The study focused on three sub-regions of the NIRL ecosystem, i.e., the southern Mosquito Lagoon; the NIRL proper and the Banana River Lagoon (Figure 1). The three primary sampling sites were Site 1 in the southern Mosquito Lagoon, Site 2 in the Indian River Lagoon near Titusville, and Site 3 in the central Banana River lagoon (Figure 1). For broader comparisons of Top-200 biomass peaks of individual taxonomic groups in the three sub-regions, data for an additional site in each sub-region were included, i.e., Site 1b in the central Mosquito Lagoon, Site 2b in the Indian River Lagoon at Cocoa and Site 3b in the northern Banana River near Cocoa Beach (Figure 1).

The three sub-regions of the NIRL are connected by navigational canals (i.e., Haulover and Barge canals) and natural gaps in land barriers. All three sub-regions are microtidal and have long water residence times, with estimated mean water half-lives (i.e., E₅₀—50% exchange) of 75 days in the southern Mosquito Lagoon, 107 days in the NIRL, and 156 days in the Banana River Lagoon (Sheng and Davis, 2003; Steward et al., 2005). All three sub-regions included in the study are characterized by small watersheds relative to the size of the receiving basins, but vary in terms of land-use (Sigua and Tweedale, 2003; Adkins et al., 2004; Steward et al., 2005; Gao, 2009). The Banana River Lagoon watersheds have significant urban/residential areas, along with some agricultural uses and undeveloped lands. The watersheds of the southern Mosquito Lagoon and NIRL near Titusville are characterized by higher percentages of natural upland and wetland preserves (mostly associated with the Canaveral National Seashore and Merritt Island National Wildlife Refuge) than the Banana River Lagoon. The southern Mosquito Lagoon does experience some tidal and wind–driven water exchange with the northern Mosquito Lagoon, which is characterized by moderate urban and residential development.

Sampling locations and frequencies varied over the study period (i.e., September 1997 through December 2020). Site 1 was sampled twice per month from May 2006 through December 2020. Sites 2, and 3 were sampled approximately monthly from September 1997 through April 2006, then twice per month from May of 2006 through December 2020.

Salinity (practical salinity units, psu) was measured on site with YSI or Hach/Hydrolab environmental multi-probes. Water temperatures described in this paper were based on the USGS continuous monitoring site in the Haulover Canal near Titusville¹ through September 2020², after which temperature values are from a continuous monitoring YSI sonde platform located near sampling Site 2 in the NIRL near Titusville (St. Johns River Water Management District, Palatka, FL, United States). Light transmission through the water column was measured using a Secchi disk (meters).

Water was collected using a vertically integrating sample tube that captures water evenly from the surface to within 0.2 m of the bottom. Split phytoplankton samples were preserved on site, one with Lugol’s solution and the other with glutaraldehyde in 0.1 M sodium cacodylate buffer. Aliquots of water were kept on ice and processed upon return to the laboratory for analysis of picocyanobacteria densities.

Data for total phosphorus (TP) and total nitrogen (TN = total Kjeldahl nitrogen + NO_x) were obtained from the St. Johns River Water Management District, Palatka, Florida,³ which carries out contemporaneous monthly analyses of water chemistry on samples collected during this project. Data from May 1997 through April 2020 were selected.

Rainfall data for the Titusville Meteorological Station, located near the NIRL sampling Site 2, was obtained from the NOAA Climatological Data for Florida web site.⁴ Rainfall totals at Titusville over the study period are shown for 12 month groupings from fall through summer (September–August). The time interval was selected to allow for the incorporation of fall seasonal rain, which can affect nutrient loads, which in turn influence the potential for winter/spring blooms because of very long water residence times in NIRL.

General phytoplankton composition was determined using the Utermöhl method (Utermöhl, 1958). Samples preserved in Lugol’s were settled in 19-mm diameter cylindrical chambers. Phytoplankton cells were identified and counted at 400× and 100× with a Leica phase contrast inverted microscope. At 400×, a minimum of 100 cells of a single taxon and 30 grids were counted. If 100 cells of a single taxon were not counted by 30 grids, up to a maximum of 100 grids were counted until 100 cells of a single taxon were reached. At 100×, a total bottom count was completed for taxa >30 μm in size. Light microscopy was aided by other techniques for proper identification, such as the squash technique and scanning electron microscopy (Badylak et al., 2004). Scanning and transmission electron microscopy and fluorescence microscopy was used to confirm the taxonomic composition of several blooms dominated by pico-nanoplanktonic species, particularly in the post-2010 period. For example, the identity of Aureoumbra lagunensis as the dominant species in the 2012 and 2013 was further confirmed using SEM/TEM (Phlips and Badylak, 2012). Fluorescence microscopy was used to enumerate picoplanktonic cyanobacteria (e.g., Synechococcus spp. and spherical picocyanobacteria spp.) at 1000x magnification (Phlips et al., 1999). Subsamples of seawater were filtered onto 0.2-μm Nucleopore filters and mounted between a microscope slide and cover slip with immersion oil.

Cell biovolumes (μm^–3 cell^–1) were estimated by assigning combinations of geometric shapes to fit the characteristics of individual taxa (Smayda, 1978; Sun and Liu, 2003). Specific phytoplankton dimensions were measured for at least 30 randomly selected cells. Species which vary substantially in size, such as many diatom species, were placed into size categories. Phytoplankton biomass as carbon values (μg carbon ml^–1) was estimated by using conversion factors for different taxonomic groups applied to biovolume estimates, i.e., 0.065 × biovolume (μm^–3 ml^–1 × 10^–6) of diatoms, 0.16 × biovolume of dinoflagellates and 0.22 × biovolume of cyanobacteria and other phytoplankton taxa (Strathmann, 1967; Ahlgren, 1983; Sicko-Goad et al., 1984; Verity et al., 1992; Work et al., 2005).

For the purpose of description and discussion, “blooms” were defined as phytoplankton carbon values which fell within the top 20% of total phytoplankton carbon values in all samples from the three main sites in NIRL examined over the study period, i.e., >2 μg carbon ml^–1.

All three sub-regions of NIRL showed similar trends in biomass peaks over the study period, with more intense (i.e., biomass peaks) and prolonged blooms in the 2011–2020 time period than from 1997 to 2010 (Figure 2). Changepoints were identified in the mean and variance of the total phytoplankton biomass time series for the average three main sites combined (Sites 1, 2, and 3) in March 2011 and October 2015 using the regression-based changepoint detection approach (Figure 3A), and in April 2011, June 2012, and September 2015, using the tests for change in mean and variance, mean only, and variance only, respectively (Figure 3B). These results highlight a major shift in bloom biomass beginning in 2011, with further intensification in subsequent years. The major trends in blooms characteristics at the three principal sampling sites include;

Surface water temperatures in NIRL typically ranged from winter lows between 10 and 15°C to summer highs between 30 and 32°C (Figure 4). However, during the winters of 2009/2010 and 2010/2011, minimum water temperatures were exceptionally low, i.e., 3.7 and 7.0°C, respectively.

Rainfall totals ranged from 100 to 180 cm (Figure 5). Totals were generally lower during periods dominated by La Niña conditions, such as 1998–2001 and 2007–2011, than during periods dominated by El Niño conditions, such as 2002–2006, 2014–2016, and 2018–2019 (Figure 6). Another notable feature of the rainfall trends is the influence of hurricanes on total rainfall (Figure 5). Certain storm events had a strong contribution to rainfall totals in the study region, including the active tropical storm periods of 2004 (Charley, August; Francis, September; Ivan, September; Jeanne, September), 2005 (Wilma, October), 2008 (Fay, August), 2012 (Beryl, May; Sandy, October), 2014 (Arthur, July), 2015 (Erika, September; Joaquin, October), 2016 (Colin, June; Julia, September; Mathew, October), 2017 (Emily, August; Irma, September), 2018 (Gordon, September), 2019 (Dorian, September), 2020 (Isaias, August)

The 23-year study period was characterized by two broad trends in salinity related to multi-year shifts in rainfall, i.e., (1) increasing salinities during generally below average rainfall periods, including 1998 through the summer of 2001, 2006 through the summer of 2008, 2010 through the summer of 2014, and 2016 and (2) declining salinities during generally above average rainfall periods, including fall 2001 through 2005, 2008–2009, fall of 2014 through 2015, 2017, and 2020 (Figure 7). The trends were subject to episodic departures related to shorter-term rainfall events, e.g., the strong Tropical Storm Fay in the late summer of 2008, which resulted in a dip in salinity during a period of generally increasing values. Salinity ranges and amplitudes of variability also varied between sites. Salinities were highest at Site 1, with peaks up to 46 psu, followed by Site 2, and lowest at Site 3 (Figure 7). Differences in salinity patterns between the three sub-regions of NIRL reflect differences in hydrologic characteristics (e.g., mean depths, water residence times, and watershed characteristics) and fine spatial scale variability in rainfall patterns.

Secchi disk depths varied between 0.5 and 2.5 m in all three sub-regions of NIRL over the first 13 years of the study (Figure 8). After 2010, extended periods of Secchi depths below 0.5 m were observed in all three sub-regions in association with major phytoplankton bloom events, with values as low as 0.1 m.

Overall patterns of TP concentrations were similar among sub-regions, with a strong rise in 2010–2013, followed by multiple major peaks in subsequent years (Figure 9). From 1997 to 2010, TP concentrations ranged from 0.004 to 0.076 mg L^–1 at Site 1 in southern Mosquito Lagoon, 0.006 to 0.105 mg L^–1 at Site 2 in the northern IRL near Titusville, and 0.008 to 0.125 mg L^–1 at Site 3 in the central Banana River Lagoon. Mean values for 1997 to 2010 were 0.031, 0.036, and 0.051 mg L^–1, for Sites 1, 2 and 3, respectively (Table 3). From 2011 to 2020, TP concentrations ranged from 0.010 to 0.205 mg L^–1 at Site 1 in southern Mosquito Lagoon, 0.007 to 0.185 mg L^–1 at Site 2 in the northern IRL near Titusville, and 0.011 to 0.227 mg L^–1 at Site 3 in the central Banana River (Figure 9). Mean values for 2011 to 2020 were 0.073, 0.081, and 0.084 mg L^–1, for Sites 1, 2, and 3, respectively (Table 3). Mean values for the 2011–2020 period were higher than the 1997–2010 period by 2.35-times at Site 1, 2.25-times at Site 2, and 1.65-times at Site 3.

Total nitrogen (TN) concentrations ranged from below 0.5 to near or above 3.0 mg L^–1 at all three sites in the study (Figure 9). Peak concentrations gradually declined over the first 13 years of the study at all three sites, then rose again after 2010. From 1997 to 2010, TN concentrations ranged from 0.4 to 2.4 mg L^–1 at Site 1 in southern Mosquito Lagoon, 0.2 to 2.3 mg L^–1 at Site 2 in the northern IRL near Titusville, and 0.2 to 2.5 mg L^–1 at Site 3 in the central Banana River (Figure 7). Mean values for 1997 to 2010 were 1.32, 1.31, and 1.50 mg L^–1, for Sites 1, 2, and 3, respectively (Table 3). From 2010 to 2020, TN concentrations ranged from 0.6 to 2.9 mg L^–1 at Site 1 in southern Mosquito Lagoon, 0.2 to 3.3 mg L^–1 at Site 2 in the northern IRL near Titusville, and 0.2 to 3.2 mg L^–1 at Site 3 in the central Banana River (Figure 9). Mean values for 2011 to 2020 were 1.51, 1.46, and 1.40 mg L^–1, for Sites 1, 2, and 3, respectively (Table 3). Mean values for the 2011–2020 period were higher than the 1997–2010 period by 1.14-times at Site 1, and 1.11-times at Site 2, but declined by 0.93-times at Site 3.

Differences in mean TN and TP concentrations from the 1997–2010 and 2011–2020 time periods were further reflected in the TN/TP ratios (Table 3). Mean TN/TP ratios were 53.5, 42.5, and 33.1 at Sites 1-3, respectively, in 1997–2010, then declined to 23.5, 19.9, and 18.7 at the three sites in 2011–2020.

One of the prominent cyclical patterns observed over the 23-year time series was higher mean total phytoplankton biomass values during El Niño periods, (i.e., 2002–2006, 2009–2010, 2014–2016 and 2018–2019) than La Niña periods (1999–2001, 2007–2008, 2011–2014, and 2017–2018) (Table 4). Mean biomass values for the pre- and post- initial changepoint (2011) periods exhibited differences in range, but the same relationships between La Niña and El Niño periods. El Niño periods are generally associated with declining salinities, which reflect overall elevated levels of rainfall and external watershed inputs. La Niña periods exhibit the opposite trend.

The dramatic change observed in bloom characteristics in 2011 fit the general criteria for a regime shift (aka, alternative stable state change), i.e., a “stark change” in the structure and function of an ecosystem inferred from long-term empirical data (MacNally et al., 2014). The results indicate two important differences in the magnitude of phytoplankton blooms and their relationship to nutrient concentrations for the pre-2011 (1997–2010) and post-2010 (2011–2020) time periods, i.e., (1) higher range of TP concentrations in the post-2010 period than pre-2011 period in all three sub-regions of NIRL, and a similar pattern for TN at Site 1, but less pronounced at Sites 2 and 3 and (2) higher phytoplankton biomass maxima at any given nutrient level in the post-2010 period for both TP and TN (Figure 10). These relationships are further exemplified by the regression relationships between annual average TP and TN concentrations and annual average phytoplankton biomass in the 1997–2010 and 2011–2020 time periods, which show steeper slopes for the post-2010 period (Figure 11). Slopes for both the pre- and post- change point periods regressions were significantly different than zero (p < 0.01) for all of the regressions except that for TN in the pre-period (Figure 11C). R² values for the regressions with TP were higher than for TN.

From a seasonal perspective, the highest peaks in phytoplankton biomass in NIRL typically occur during the warmer months of the year, which coincides with the wet season (May–October), when external nutrient loads are elevated (Phlips et al., 2010). However, due to the location of NIRL in the subtropical region of east Florida, relatively modest seasonal variability in water temperature and incident irradiance enhance the potential for winter bloom events. During this study, winter bloom events were often dominated by cosmopolitan species with tolerance to lower temperatures, such as the diatom C. pelagica, and the dinoflagellate A. sanguinea, which have broader temperature preferences than the sub-tropical/tropical species found in NIRL (Badylak and Phlips, 2004; Phlips et al., 2010, 2015; Badylak et al., 2014). The potential for cold season blooms in NIRL is further illustrated by the intense A. lagunensis brown tides of 2016 and 2018, both of which began in the winter. By contrast, blooms in the warm season were often dominated by the tropical dinoflagellate species P. bahamense (Phlips et al., 2006).

In the future, progressive warming associated with climate change may result in further expansion of the peak bloom season (Jahan and Choi, 2014). The warming trend may also increase the frequency and intensity of tropical storms and hurricanes (Phlips et al., 2020). High rainfall and winds associated with major storms can increase external and internal (e.g., sediment resuspension and benthic biomass disruption) nutrient loads that support bloom development. This is illustrated by the occurrence of two major hurricanes (Emily and Irma) in Florida in 2017, which was followed by a major bloom event in NIRL starting in the winter of 2017/18 and extending into 2018.

On a longer time scale, another element of cyclical variability of phytoplankton blooms is the influence of El Niño and La Niña periods. A number of recent studies of coastal ecosystems have linked El Niño/La Niña cycles to changes in phytoplankton biomass and composition (Zhu et al., 2017; Zhao et al., 2018; Jiménez-Quiroz et al., 2019; Phlips et al., 2020). In the Southeastern United States, El Niño and La Niña periods are associated with wetter and dryer than average conditions, respectively (Piechota and Dracup, 1996; Schmidt and Luther, 2002). The pattern results in above and below average watershed inflows to coastal aquatic ecosystems (Kahya and Dracup, 1993; Dracup and Kahya, 1994). Nutrients enter NIRL from a variety of external sources, including surface water runoff, tributary inflows, groundwater inputs, rainfall, septic system seepage, and periodic permitted and accidental releases from sewage treatment systems (Swarzenski et al., 2001; Adkins et al., 2004; Steward and Green, 2007; Gao, 2009; Lapointe et al., 2015; Trefry and Fox, 2021). All of these inputs can be exacerbated by high rainfall levels. El Niño periods are noted for elevated rainfall levels in the dry season (i.e., late fall through early spring). Given the very long water residence times in NIRL (E₅₀ > 3 months), elevated nutrient loads induced by rainfall in the dry season can enhance the potential for blooms from winter through early summer.

In this study, El Niño periods were characterized by higher peaks in bloom biomass (e.g., 2002–2006, 2009–2010, 2014–2016, 2018–2019), than in La Niña periods (e.g., 1999–2001, 2007–2008, 2011–2014, 2017–2018). One of the most dramatic examples of this relationship is the intense and protracted bloom in NIRL in 2015–2016, which began in the winter of 2015/16. This relationship is further illustrated by the positive correlation between the ENSO index and peak biomass of one of the dominant bloom-forming species in NIRL, the dinoflagellate P. bahamense (Phlips et al., 2006, 2015, 2020). A similar positive relationship between periods of elevated rainfall and blooms of P. bahamense has been reported in the Indo-Pacific region (Usup et al., 2012). Similar trends have been reported for rates of primary production in other coastal ecosystems (Belgrano et al., 1999; Paerl et al., 1999). By contrast, high rainfall El Niño periods in Baffin Bay in Texas appear to be associated with less intense brown tide blooms due to elevated flushing rates and reduced salinities (Cira et al., 2021). These contrasting impacts of El Niño/La Niña periods highlights the importance of recognizing regional differences in key ecosystem characteristics when determing the effects on HAB dynamics.

In addition to the aforementioned cyclical patterns, a sudden, stark and persistent shift in the biomass, composition and timing of phytoplankton blooms in NIRL was observed in 2011. The highest peaks in biomass in the post-2010 period were up to five or more times higher than prior to 2011, and the blooms expanded into the winter and spring seasons. From the perspective of phytoplankton community structure, an important feature of the post-shift period was the prominence of pico/nano-planktonic algae, including the brown tide species A. lagunensis, pico/nano-planktonic cyanobacteria and nano-planktonic eukaryotes. By contrast, in the pre-shift period microplanktonic dinoflagellates and diatoms were the dominant taxa in blooms. The results of this study provide insights into the drivers of the shift and how it fits into the concept of a regime shift.

Regime shifts in marine ecosystems have been attributed to a range of drivers of varying scale and character (Conversi et al., 2015; McMahon et al., 2015; Möllmann et al., 2015). Reports of regime shifts in planktonic communities have stressed the influences of climatic changes or events, such as well-documented long-term studies of temperature changes in the NE Atlantic Ocean (Beaugrand, 2004; Molinero et al., 2013; Alheit et al., 2019; Bode et al., 2020). Many reports of shifts also highlight the role that regional drivers, such as anthropogenic eutrophication and hydrologic alterations, can play in driving major and persistent changes in the structure and function of planktonic communities (Jahan and Choi, 2014; Phlips et al., 2015, 2020; Peperzak and Witte, 2019; Bode et al., 2020). There are frequently multiple drivers of change that operate synchronously (Beaugrand et al., 2015; Möllmann et al., 2015), and can lead to a cascade of trophic level impacts (Pershing et al., 2015). This is illustrated by the multi-year intense blooms of A. lagunensis in Laguna Madre lagoon in Texas that negatively impacted seagrasses and certain faunal communities (Buskey et al., 1998, 2001).

The regime shift in the intensity of phytoplankton blooms in NIRL involves both climatic and regional drivers. The increases in bloom biomass in 2011 coincided with a large increase in nutrient concentrations. Annual average TP concentrations after 2010 were two to three times higher than before 2011 and annual average TN concentrations were as high as any observed before 2010. The large increases in nutrients throughout NIRL in 2011 do not appear to be related to sudden increases in external load since it coincided with a largely below average rainfall period. Normally, low rainfall periods in NIRL are characterized by reduced external nutrient loads and comparatively low nutrient concentrations and phytoplankton biomass (Steward et al., 2005; Steward and Green, 2007; Phlips et al., 2010, 2015). This conundrum warrants a closer examination of changes in sources and sinks of nutrients underlying the regime shift in bloom characteristics in 2010. One major change in NIRL that coincided with the increases in nutrient levels and bloom intensities was widespread losses in seagrass communities, drift macroalgae and some pelagic fauna, including zooplankton and fish, in 2009–2011 (Phlips et al., 2015; Stevens et al., 2016). Major biomass losses in these communities appear to have precipitated a cascade of impacts, beginning with a re-distribution of bioavailable nutrients toward the phytoplankton community, and culminating in the regime shift in bloom intensities (Phlips et al., 2015). Estimates of the decline in seagrass, seagrass epiphytes and drift macroalgae biomass from 2009 to 2011 approximate the observed increase in peak phytoplankton biomass levels in 2011, providing supporting evidence for the latter hypothesis (Phlips et al., 2015).

Shifting of nutrient pools from benthic primary producer biomass into phytoplankton biomass has been observed in other eutrophic aquatic ecosystems subject to persistent phytoplankton blooms (Burkholder et al., 2007; Waycott et al., 2009; Schmidt et al., 2012). Losses of benthic primary producer biomass can enhance nutrient availability for phytoplankton blooms in several other ways. Bioavailable forms of nutrients entering the water column are competed for by planktonic and benthic primary producers, most prominently drift macroalgae, microphytobenthos (e.g., microalgae on the sediment surface) and epiphytes associated with seagrass communities (Phlips et al., 2015), therefore, losses in the latter communities can enhance nutrient availability for phytoplankton. Seagrass and microphytobenthos communities provide stabilization of surface sediments, therefore their demise decreases the benthic nutrient filter effect, resulting in increased fluxes of nutrients from sediments into the water column (Macintyre et al., 2004; Sundbäck et al., 2004; Murrell et al., 2009; Eyre et al., 2011; Anderson et al., 2014). Sediment destabilization also can enhance the potential for re-suspension of nutrient-rich muck sediments (Trefry et al., 1990; Aoki et al., 2020). In addition, periodic reductions in oxygen concentrations near the sediment surface associated with major phytoplankton blooms can enhance mobilization of nutrients from sediments into the water column (Reay et al., 1995; Engelsen et al., 2008; Cook et al., 2009). The senescence of blooms also contributes to the accumulation of nutrients in sediments, thereby enhancing legacy loads in the ecosystem (Lemley et al., 2021). All of these factors likely contributed to the additional nutrients needed for the exceptionally intense phytoplankton blooms observed in NIRL after 2010.

A number of factors likely contributed to the 2009–2011 declines in benthic primary producer communities in NIRL, including exceptionally low winter temperatures in 2009/2010 and 2010/2011 (i.e., temperature minima near 5°C) and above average salinities (i.e., salinities > 35 psu) due to below average rainfall levels. These observations also highlight a more general issue involving the response of subtropical ecosystems to low temperature excursions, because of the presence of large numbers of flora and fauna species of tropical origin sensitive to low temperature extremes. Analogous observations have been made in relation to uncharacteristically high temperature conditions in temperate ecosystems (Rhodes et al., 1993; Nehring, 1998; Cloern et al., 2005; Edwards et al., 2006). Low light availability for benthic primary producer communities caused by the protracted phytoplankton bloom in 2011 exacerbated losses of seagrass and drift macroalgal biomass (Hall et al., 2021). Continued intense phytoplankton blooms after 2011 have prolonged the problem of light limitation. The intense phytoplankton blooms in NIRL during the post-2010 period were associated with very low Secchi disk depths, i.e., 0.1–0.5 m. Based on the average depths of the sub-regions of NIRL (i.e., 2 m), such low Secchi depths create a condition where incident irradiances near the sediment surface are below that needed to support growth and survival of benthic primary producer communities, i.e., near 10–15% of surface irradiance (Dennison, 1987; Duarte, 1991; Steward et al., 2005; Duarte et al., 2007; Lee et al., 2007; Breininger et al., 2016). It is therefore not surprising that little to no recovery of seagrasses has been observed in annual surveys of distribution and biomass from 2011 to 2020 (Phlips et al., 2015). Similar issues have been raised for the 8-year brown tide event in Laguna Madre, Texas (Onuf, 1996).

The other major attribute of the regime shift in 2011 was the dominant role played by pico-nano planktonic phytoplankton in blooms. A similar dominance by nanoplanktonic species was observed during the eight-year brown tide event in Laguna Madre, Texas (Buskey et al., 2001). As observed in NIRL, the bloom in Texas was associated with fish kills and widespread losses of seagrass biomass. It has been hypothesized that elevated levels of ammonium and organic forms of nitrogen and phosphorus associated with recycled biomass (e.g., decaying seagrass, benthic algae, and senescing phytoplankton bloom biomass) provide a competitive advantage for brown tide and other pico-nanoplanktonic species in competition for nutrients (DeYoe and Suttle, 1994; Liu et al., 2001; Berg et al., 2002; Buskey et al., 2003; DeYoe et al., 2007; Glibert et al., 2007; Muhlstein and Villareal, 2007; Sun et al., 2012; Kang et al., 2015; Wetz et al., 2017). A similar argument may apply to pico/nano-planktonic bloom-forming species in NIRL post-2010 (Phlips et al., 2015; Papacek, 2018). Many species in these groups are flexible in their methods of obtaining nutrients, such as the use of organic forms, nitrogen fixation, and mixotrophic strategies (e.g., heterotrophy, osmotrophy, and phagotrophy; Phlips et al., 1989; Stolte and Riegman, 1996; Chung et al., 2003; Moore et al., 2005; Dyhrrnan and Ruttenberg, 2006; Burkholder et al., 2008; Sunda and Hardison, 2010; Kathuria and Martiny, 2011; Lindehoff et al., 2011; Yuan et al., 2012; Stoecker et al., 2017; Papacek, 2018). Many pico/nano-phytoplankton species are also known to be efficient in the use of phosphorus (Chung et al., 2003; Moore et al., 2005; Dyhrrnan and Ruttenberg, 2006; Sunda and Hardison, 2010; Kathuria and Martiny, 2011), including the brown tide species A. lagunensis that dominated many of the blooms in NIRL after 2011 (Kang et al., 2015). This observation helps to explain why some of the peaks in bloom biomass from 2012 to 2020 were high relative to TP concentrations, compared to blooms of micro-phytoplankton blooms prior to 2011.

Another factor that may have contributed to the exceptionally high phytoplankton biomass levels in the post-2010 period is reduction in top-down control. An ability to deter grazing has been reported for a range of bloom species (Turner and Tester, 1997; Buskey, 2008; Smayda, 2008; Sunda and Hardison, 2010), including the brown tide species A. lagunensis that has been a dominant feature of blooms in NIRL after 2010 (Buskey et al., 1997; Kang et al., 2015; Galimany et al., 2020). The effect may also extend to other taxa, such as certain cyanobacteria species (Urrutia-Cordero et al., 2016). The impact of variability in top-down versus bottom-up pressures on phytoplankton biomass has been observed in a number of marine ecosystems (Pershing et al., 2015).