Phytoplankton Bloom Dynamics in Incubated Natural Seawater: Predicting Bloom Magnitude and Timing

Ok, Jin Hee; Jeong, Hae Jin; You, Ji Hyun; Kang, Hee Chang; Park, Sang Ah; Lim, An Suk; Lee, Sung Yeon; Eom, Se Hee

doi:10.3389/fmars.2021.681252

ORIGINAL RESEARCH article

Front. Mar. Sci., 28 July 2021

Sec. Coastal Ocean Processes

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.681252

This article is part of the Research Topic Eutrophication, Algal Bloom, Hypoxia and Ocean Acidification in Large River Estuaries View all 11 articles

Phytoplankton Bloom Dynamics in Incubated Natural Seawater: Predicting Bloom Magnitude and Timing

$\r\nJin Hee Ok$ Jin Hee Ok¹

Hae Jin Jeong^1,2*

¹School of Earth and Environmental Sciences, College of Natural Sciences, Seoul National University, Seoul, South Korea
²Research Institute of Oceanography, Seoul National University, Seoul, South Korea
³Division of Life Science, Plant Molecular Biology and Biotechnology Research Center, Gyeongsang National University, Jinju, South Korea

Phytoplankton blooms can cause imbalances in marine ecosystems leading to great economic losses in diverse industries. Better understanding and prediction of blooms one week in advance would help to prevent massive losses, especially in areas where aquaculture cages are concentrated. This study has aimed to develop a method to predict the magnitude and timing of phytoplankton blooms using nutrient and chlorophyll-a concentrations. We explored variations in nutrient and chlorophyll-a concentrations in incubated seawater collected from the coastal waters off Yeosu, South Korea, seven times between May and August 2019. Using the data from a total of seven bottle incubations, four different linear regressions for the magnitude of bloom peaks and four linear regressions for the timing were analyzed. To predict the bloom magnitude, the chlorophyll-a peak or peak-to-initial ratio was analyzed against the initial concentrations of NO₃ or the ratio of the initial NO₃ to chlorophyll-a. To predict the timing, the chlorophyll-a peak timing or the growth rate against the natural log of NO₃ or the natural log of the ratio of the initial NO₃ to chlorophyll-a was analyzed. These regressions were all significantly correlated. From these regressions, we developed the best-fit equations to predict the magnitude and timing of the bloom peak. The results from these equations led to the predicted bloom magnitude and timing values showing significant correlations with those of natural seawater in other regions. Therefore, this method can be applied to predict bloom magnitude and timing one week in advance and give aquaculture farmers time to harvest fish in cages early or move the cages to safer regions.

Introduction

Phytoplankton are a major component of marine ecosystems and contribute approximately half of the global primary productivity (Field et al., 1998), providing organic carbon and energy for higher trophic levels in food webs (Calbet and Landry, 2004; Steinberg and Landry, 2017; Armengol et al., 2019). Many phytoplankton species cause blooms or red tides, which are water discolorations due to blooms. Harmful algal blooms often cause human illness and fish kills (Anderson et al., 2002; Hallegraeff, 2003; Glibert et al., 2018b). Understanding and predicting the outbreak of phytoplankton blooms, especially harmful algal blooms, are important to minimize the damage they may cause.

In general, the abundance of a phytoplankton species or community is typically affected by nutrient conditions (Nagasoe et al., 2010; Gobler et al., 2012; Lee et al., 2019a), as well as physical and trophic conditions such as light availability and predation (Cole and Cloern, 1984; Cloern, 1987; Buskey et al., 1994; Litchman, 1998; Turner and Granéli, 2006). The generation times of many phytoplankton species are relatively short; some phytoplankton can divide three to four times per day, and thus, their abundance rapidly increases with photosynthesis-favorable conditions such as sunny days after heavy rains (Trainer et al., 1998; Wetz and Paerl, 2008; Baek et al., 2009; Meng et al., 2017). A phytoplankton bloom outbreak cannot be easily predicted if unexpected meteorological events, such as heavy rains, typhoons, and heatwaves, occur in a short period. A few short-term forecasts (< one week) that predict the outbreaks of harmful algal blooms in the United States have been maintained by the National Oceanic and Atmospheric Administration (NOAA; Wynne et al., 2018). If the magnitude and timing of phytoplankton bloom outbreaks can be predicted quickly for other regions and announced one week prior, aquaculture farmers can harvest fish early, move the fish cages to safer areas, and use developed techniques to mitigate algal blooms. Prior to the present study, models for predicting blooms had been established based on correlations between phytoplankton abundance and various environmental factors such as nutrient concentrations and water temperature in natural aquatic environments (e.g., Pinckney et al., 1997; Onderka, 2007; Feki-Sahnoun et al., 2017; Lin et al., 2018; Kahru et al., 2020; Mahmudi et al., 2020). However, these models may not provide sufficient time for aquaculture farmers to manage their cages at sea or aqua tanks on land. In general, blooms by particular phytoplankton species kill fish in aquaculture cages when the abundance of phytoplankton exceeds a certain critical concentration (Gobler et al., 2008; Gravinese et al., 2019). Thus, to minimize losses due to harmful blooms, a model for predicting the magnitude and timing of blooms one week in advance is needed.

Gamak Bay and Yeosuhae Bay, near Yeosu city, South Korea, are ideal regions for developing and testing methods to predict the magnitude and timing of potential phytoplankton bloom outbreaks. Freshwater from the large Seomjin River enters these bays and often causes large variations in the nutrients and chlorophyll-a concentrations (Lee M. et al., 2009; Lee Y.S. et al., 2009; Noh et al., 2010). Meanwhile, commercial fish and oyster aquaculture systems are concentrated in these bays (Lee M. et al., 2009). There have been many harmful algal blooms in these bays National Institute of Fisheries Science, Korea (2020)¹ causing great financial losses in the aquaculture industry (Park et al., 2013).

In the present study, we quantified the concentrations of nutrients and chlorophyll-a (Chl-a), as well as the abundance of dominant plankton species, daily or every other day for 10 days after incubating seawater collected from Gamak Bay and Yeosuhae Bay, South Korea. During the study period, the nutrient conditions of sampled seawater varied greatly; nutrient levels were low during dry periods and before the passage of a typhoon, but high after typhoons and rainfall. We used the results of enclosure experiments and tested the correlations between pairs of different variables (e.g., the initial nutrient concentrations and the initial and peak concentrations of Chl-a and ratios) to propose equations in order to determine the magnitude (i.e., Chl-a peak) and the timing of phytoplankton bloom peaks. To validate these equations empirically, we compared the predicted values of magnitude and timing with the actual values from the monitoring sites from May to October annually from 2016 to 2019. The results of this work provided a basis for understanding bloom dynamics in coastal waters and contributed to being able to predict phytoplankton blooms in advance. This can help ensure that preemptive actions are taken to safeguard aquaculture operations.

Materials and Methods

Sampling Site Description and Water Sampling

The two sampling sites, Station Gukdong and Station Gyedong, were located in the innermost areas of Gamak Bay and Yeosuhae Bay, South Korea, respectively (Figures 1A,B). These bays are located near the mouth of the Seomjin River, which is one of the largest rivers in Korea, from which freshwater runoff reaches the bays after heavy rainfall (Jeong et al., 2017). Twenty liters of seawater were collected below the water surface at St. Gukdong using a clean bucket, six times across May, July, and August 2019, and 20 L of seawater was gently collected below the surface at St. Gyedong, once in August 2019 (Table 1). The temperature, salinity, pH, and dissolved oxygen (DO) of the samples were measured using a YSI Professional Plus meter (YSI, Yellow Springs, OH, United States) immediately after collection. The seawater samples were gently filtered through a 100-μm Nitex mesh to minimize the effects of large detritus or solids in suspension on phytoplankton growth during the bottle incubation experiments (e.g., shading, flocculation, and nutrient adsorption) (Young and Barber, 1973; Sew and Todd, 2020). The collected seawater was then immediately transported to the laboratory and placed in culture chambers at the same temperature (±2°C) as the seawater temperature at collection. The amount of rainfall and the wind speed and path of Typhoon Danas during the sampling period were obtained from the Korea Meteorological Administration (2020)² and the Korea Hydrographic and Oceanographic Agency (2021).³

FIGURE 1

Figure 1. (A) Locations of the study areas in the South Sea of Korea. (B) Stations for seawater sample collection in Gamak Bay (St. Gukdong) and Yeosuhae Bay (St. Gyedong) and the station for automated seawater quality monitoring data acquisition in Gwangyang Bay (St. Gwangyang). (C) Station for the automated seawater quality monitoring data acquisition in Masan Bay (St. Masan). The dotted line in panel (A) indicates the path of the Danas typhoon on July 20, 2019.

TABLE 1

Table 1. Meteorological conditions when the seawater samples were collected from Station Gukdong and Station Gyedong in 2019 and their incubation conditions.

Experimental Procedures

Seawater collected on each sampling date was gently mixed and evenly distributed into 2 L polystyrene bottles (JetBiofiltration Co., Ltd, Guangzhou, China). All bottles incubated under each condition were set up in triplicates and placed in temperature-controlled chambers set to the ambient water temperature (±2°C; Table 1). These bottles were capped loosely and then incubated for 10 days at 100 μmol photons m^–2 s^–1 under a 14 h:10 h light:dark cycle using light-emitting diodes (LEDs; FS-075MU, 6500K; Suram Inc., Suwon, South Korea). Based on the light intensity at which the maximum growth rates of many phytoplankton occurred and the effects of photoinhibition on phytoplankton growth were minimized, a light intensity of 100 μmol photons m^–2 s^–1 was used. This ensured that light was not a limiting factor in this study. For example, the light intensities at which the maximum growth rates were observed were >10 μmol photons m^–2 s^–1 for the dinoflagellate Paragymnodinium shiwhaense (Jeong et al., 2018), >58 μmol photons m^–2 s^–1 for the dinoflagellate Alexandrium pohangense (Lim et al., 2019), >70 μmol photons m^–2 s^–1 for the diatom Leptocylindrus danicus (Verity, 1982), >80 μmol photons m^–2 s^–1 for the diatoms Skeletonema costatum, Thalassiosira minima, and Thalassiosira eccentrica, and the dinoflagellates Protodinium simplex, Scrippsiella sweeneyae, and Prorocentrum micans (Chan, 1978), and >90 μmol photons m^–2 s^–1 for the dinoflagellate Margalefidinium polykrikoides (Kim et al., 2004). Moreover, the light intensities at which the negative growth rates of phytoplankton due to photoinhibition occurred were >123 μmol photons m^–2 s^–1 for the dinoflagellate Karenia brevis (Magaña and Villareal, 2006), >135 μmol photons m^–2 s^–1 for the dinoflagellate Tripos furca (Nordli, 1957), >247 μmol photons m^–2 s^–1 for the dinoflagellate Takayama helix (Ok et al., 2019), and >317 μmol photons m^–2 s^–1 for the diatoms Ditylum brightwellii and Thalassiosira sp. (Brand and Guillard, 1981). Daylength in the study area from May to August is approximately 13–14.5 h (Sunrise-Sunset, 2021).⁴ Thus, this light condition can simulate the natural environment in the study area.

Subsampling and Component Analyses

Subsampling was conducted daily in all experimental seawaters except for the one collected on May 16, 2019 (every two days) to quantify the Chl-a and nutrient concentrations and the planktonic protist abundance.

For Chl-a analysis, a 20 mL aliquot taken from each bottle was filtered through a GF/F filter (Whatman Inc., Clifton, NJ, United States), and the filtered membrane was placed in a 15 mL conical tube and stored below −20°C. Then, 10 mL of 90% acetone was added to the conical tube containing the sample. The conical tubes were sonicated for 10 min, wrapped in aluminum foil to avoid light, stored in a 4°C chamber overnight, and centrifuged, and then 5 mL of the supernatant was taken from each tube and measured using a 10-AU Turner fluorometer (Turner Designs, Sunnyvale, CA, United States). For nutrient analysis, a 20 mL aliquot of the filtrate samples was stored in polyethylene vials below −20°C. The concentrations of nitrate plus nitrite (hereafter nitrate or NO₃), ammonium (NH₄), phosphate (PO₄), and silicate (SiO₂) in the seawater samples were measured using a 4-channel nutrient auto-analyzer (QuAAtro, SEAL Analytical GmbH, Norderstadt, Germany; Bran+Luebbe Analyzing Technologies, 2004a, b, 2005; SEAL Analytical GmbH, 2005).

To quantify the abundance of dominant planktonic protists, a 10 mL aliquot was taken from each bottle and fixed with 5% acidic Lugol’s solution (Sournia, 1978). Each dominant taxon was enumerated by counting >200 or all cells for each species in a Sedgewick-Rafter chamber using an inverted microscope (BX53, Olympus, Japan) at a magnification of ×100–200.

To convert the abundance (cells mL^–1) of dominant planktonic protists to biomass (ng C mL^–1), the cellular carbon content of each taxon was obtained from the literature or estimated from its biovolume according to Menden-Deuer and Lessard (2000). The biovolume of each taxon was calculated using Hillebrand et al. (1999) after measuring the cell length and width, or diameter and height of ≥5 cells of the taxon.

Data Analysis and Equation Development

The concentrations of NO₃, NH₄, PO₄, SiO₂, and Chl-a were measured every sampling day and dissolved inorganic nitrogen (NO₃+NH₄; DIN) was calculated from these data. When Chl-a increased, the magnitude of the Chl-a peak (i.e., the maximum Chl-a) and the time elapsed until the Chl-a peak occurred (Day_peak) were determined. When Chl-a decreased, the initial Chl-a was regarded as the Chl-a peak. Using the data, the ratio of the initial NO₃ to Chl-a [μM (μg L^–1)^–1] was calculated by dividing the initial NO₃ by the initial Chl-a for each incubated seawater sample. This value was also calculated as a unitless ratio, multiplied by the molecular weight of N (=14 g mol^–1) to convert μmol L^–1 (=μM) to μg L^–1. The growth rate (μ, day^–1) of the phytoplankton community was calculated as follows:

μ = \frac{Ln (\frac{Chl a peak}{initial Chl a})}{{Day}_{peak}}

When the Chl-a decreased, the growth rate was allocated to zero.

To develop equations for the magnitude of the Chl-a peak of the phytoplankton community, four different simple linear regressions were analyzed: the Chl-a peak versus the initial NO₃; the Chl-a peak versus the ratio of the initial NO₃ to Chl-a; the ratio of the peak-to-initial Chl-a versus the initial NO₃; and the ratio of the peak-to-initial Chl-a versus the ratio of the initial NO₃ to Chl-a. Furthermore, to develop an equation for Day_peak, four different simple linear regressions were also analyzed: Day_peak versus the natural log of the initial NO₃; Day_peak versus the natural log of the ratio of the initial NO₃ to Chl-a; the growth rate versus the natural log of the initial NO₃; and the growth rate versus the natural log of the ratio of the initial NO₃ to Chl-a. The growth rate was used as a variable because the growth rate equation itself includes time. Similar analyses were conducted using NH₄, DIN, PO₄ and SiO₂, as well as the ratio of the initial NO₃ to PO₄ and the ratio of the initial NO₃ to SiO₂. To apply this method at a species level, the peak abundance of Skeletonema costatum of each seawater sample during bottle incubation was determined. The goals were to develop equations for predicting bloom outbreaks of each species and to assess the contribution of each bloom species of interest to the phytoplankton community. The cellular chlorophyll-a content of S. costatum was obtained from Hitchcock (1980). The estimated Chl-a peak and Day_peak of S. costatum were determined using the analyses described above.

Linear regression equations were established by fitting the data using DeltaGraph (SPSS Inc., Chicago, IL, United States). These equations were solved to deduce the Chl-a peak and Day_peak.

Comparison of the Predicted and Field Data

To test whether the bloom prediction equations obtained from the present study (in sections “Correlations Between the Variables Related to the Chl-a Peak” and “Correlations Between Variables Related to the Timing of Chl-a Peaks”) can be used to predict bloom outbreaks in the field, in situ time-series data of the automated seawater quality monitoring system in the South Sea of Korea (St. Gwangyang and St. Masan; Figure 1) were obtained from the Marine Environment Information System (2021).⁵ Data from May to October in 2016–2019 were used because the water temperatures during this period were similar to those of our incubation experiments using the seawaters collected from St. Gukdong and St. Gyedong. The values of NO₃, Chl-a, salinity, and water temperature measured at five-minute intervals in a day were each averaged to obtain the mean value of the day (Supplementary Table 1). Based on these monitored data, the following steps were conducted (Supplementary Figure 1). First, the Chl-a peak was found when the Chl-a exceeded 5 μg L^–1 with an increasing trend, assuming the criterion of a phytoplankton bloom as 200 ng C mL^–1 and a carbon to Chl-a ratio of 40 (Peterson and Festa, 1984; Jeong et al., 2013). Second, Chl-a and NO₃ levels, prior to reaching the Chl-a peak, were selected. Third, the time elapsed until the Chl-a peak occurred (Day_peak) was counted. Fourth, using the selected NO₃ and Chl-a data and the equations developed in this study (in Sections “Correlations Between the Variables Related to the Chl-a Peak” and “Correlations Between Variables Related to the Timing of Chl-a Peaks”), the predicted values were calculated for each date and then compared with the actual Chl-a peak or Day_peak using linear regression analyses. Data were not used when the daily average salinity changed by >2 because it indicated seawater with different nutrient conditions introduced from the river or large streams.

Statistical Analysis

Pearson’s correlation was used to examine the relationships between variables (Pearson, 1895). All statistical analyses were performed using SPSS (version 25.0; IBM Corp., Armonk, NY, United States). The significance criterion was set at 0.05.

Results

Physical and Chemical Properties of the Seawater Samples

During the study period in 2019, water temperature, precipitation, and salinity varied considerably whereas pH and dissolved oxygen (DO) did not (Figure 2): the water temperature ranged from 18.2°C on May 19 to 26.6°C on August 22 (Figure 2A); the precipitation summed for 7 days was 274.1 mm on July 22 and July 24 but <100 mm on the remaining dates (Figure 2B); the salinity was 24.36 and 29.10 on July 22 and July 24, respectively, but exceeded 29.9 on the remaining dates, with a maximum of 33.78 on May 16 (Figure 2C); the pH was 7.89 on July 22 but exceeded 7.90 on the remaining dates, with a maximum of 8.10 on August 22 (Figure 2D); and the DO was 5.01 mg L^–1 on July 24 but exceeded 5.4 mg L^–1 on the remaining dates with a maximum of 6.07 mg L^–1 on May 16 (Figure 2E).

FIGURE 2

Figure 2. Background environmental conditions of water samples from May 16 to August 22, 2019. (A) Temperature (T, °C), (B) amount of rainfall for 7 days [Rainfall_7d, mm (7 d)^{− 1}], (C) salinity, (D) pH, and (E) dissolved oxygen (DO, mg L^{− 1}).

There was 65.8 mm of rainfall between May 17 and May 18, 274.1 mm between July 17 and July 21, and 31.9 mm between August 21 and 22. Typhoon Danas passed through Korea and affected the study area on July 20 (Figure 1A). The maximum wind speed during the passage of the typhoon in the study area was 21.3 m s^–1.

Dominant Protist Species at the Beginning of Incubation

At the beginning of each bottle incubation experiment, the most dominant planktonic protist species based on carbon biomass in chronological order was the diatom Cylindrotheca closterium on May 16, the photosynthetic dinoflagellate T. furca on May 19, the tintinnid ciliate Tintinnopsis sp. on July 17 and July 22, the diatoms Pseudo-nitzschia pungens on July 24 and Chaetoceros curvisetus on August 20, and the photosynthetic ciliate Mesodinium rubrum on August 22 (Figures 3A–G).

FIGURE 3

Figure 3. (A–G) Biomass of the five dominant planktonic protists at the beginning of the incubation of the seven different seawater samples collected from St. Gyedong (August 22) or St. Gukdong (the remaining dates) in 2019. To convert the unit from cell abundance to biomass, the cellular carbon content of each dominant species was obtained from Goldman (1993), Olenina et al. (2006), Yoo et al. (2010); Lee et al. (2014), Harrison et al. (2015), and Lim et al. (2019). The cellular carbon content of Alexandrium catenella was calculated from the equivalent spherical diameter in Kang et al. (2018). The cellular carbon contents of Chrysochromulina sp. (0.004 ng C cell^{− 1}), Chaetoceros sp. (0.001 ng C cell^{− 1}), Thalassiosira cf. minima (0.02 ng C cell^{− 1}), Thalassiosira sp. (0.10 ng C cell^{− 1}), Protoperidinium sp. (1.81 ng C cell^{− 1}), Strobilidium sp. (1.00 ng C cell^{− 1}), and Tintinnopsis sp. (3.78 ng C cell^{− 1}) were calculated from the biovolume in this study.

Initial Nutrient and Chl-a Concentrations

The initial concentrations of NO₃, NH₄, and DIN of the experimental seawater ranged from 0.6 μM (August 20) to 39.7 μM (July 22), from 0.1 μM (July 17) to 10.2 μM (July 22), and from 1.5 μM (July 17) to 49.9 μM (July 22), respectively (Figures 4A–C). The initial PO₄ concentrations ranged from 0.1 μM (July 24) to 1.1 μM (July 22; Figure 4D), while the initial SiO₂ concentrations varied from 4.4 μM (August 20) to 93.4 μM (July 22; Figure 4E). The initial Chl-a concentrations ranged from 0.8 μg L^–1 (May 19) to 24.5 μg L^–1 (July 24; Figure 4F).

FIGURE 4

Figure 4. Initial concentrations of nutrients (A–E; NO₃, NH₄, DIN, PO₄, and SiO₂; μM) and chlorophyll-a (F; Chl-a; μg L^{− 1}) of each seawater collected from St. Gyedong (August 22) or St. Gukdong (the rest dates) in 2019. Vertical symbols represent treatment means ± 1 SE.

Variations in Nutrient and Chl-a During Incubation

Chl-a showed two different patterns as incubation time progressed (Figures 5A–G); first, it decreased in the seawater samples from May 16, July 17, and August 20 (Figures 5A–C); second, it increased to a peak, but then decreased later in the rest of the samples as NO₃, NH₄, DIN, PO₄, or SiO₂ became depleted (Figures 5D–G).

FIGURE 5

Figure 5. Variations in the concentrations of chlorophyll-a (Chl-a) and nutrients. Chl-a (A–G), nitrate (NO₃) (H–N), ammonium (NH₄) (O–U), dissolved inorganic nitrogen (DIN) (V–AB), phosphate (PO₄) (AC–AI), and silicate (SiO₂) (AJ–AP) during incubation of seawaters. A blue arrow indicates the highest Chl-a in the seawater where any increase in Chl-a did not occur during incubation. A red arrow indicates the Chl-a peak in the seawater where increases in Chl-a occurred during incubation. The date in each box indicates the sampling date of the seawater collected from St. Gyedong (August 22) or St. Gukdong (the remaining dates) in 2019. Vertical symbols represent treatment means ± 1 SE.

Chlorophyll-a decreased during incubation when the initial NO₃ concentrations were ≤2.4 μM (Figures 5A–C,H–J and Supplementary Figure 2A). In contrast, Chl-a increased when the initial NO₃ concentrations were ≥3.5 μM (Figures 5D–G,K–N and Supplementary Figure 2B). Chl-a decreased when the initial NH₄ concentrations were ≤1.0 μM, while Chl-a increased when the concentrations were ≥2.1 μM (Figures 5O–Q,R–U). Furthermore, Chl-a decreased when the initial DIN concentrations were ≤3 μM but increased when the concentrations were ≥6.4 μM (Figure 5V–X,Y–AB). Chl-a decreased during incubation when the initial PO₄ concentrations were ≤0.3 μM (Figures 5A–C,AC–AE) with one exception: Chl-a increased when the initial PO₄ concentration was 0.1 μM in the seawater on July 24 (Figures 5F,AH). Chl-a also decreased when the initial SiO₂ concentrations were ≤24.5 μM (Figures 5A–C,AJ–AL) with one exception: Chl-a increased when the initial SiO₂ concentration was 18.7 μM in the seawater of May 19 (Figures 5D,AM).

Correlations Between the Variables Related to the Chl-a Peak

The Chl-a peak of the phytoplankton community was in between 4.6 and 114 μg L^–1 (Figures 6A,B). The ratio of the peak-to-initial Chl-a of the phytoplankton community varied from 1.0 to 61.4 (Figures 6C,D). The ratio of the initial NO₃ to Chl-a of the phytoplankton community ranged from 0.1 to 24.2 μM (μg L^–1)^–1 (equivalent to 0.8–339 in a unitless ratio, respectively; Figures 6B,D).

FIGURE 6

Figure 6. (A,B) Chl-a peak (μg L^{− 1}) or (C,D) peak-to-initial Chl-a ratio of the phytoplankton community as a function of the initial NO₃ (μM) or the ratio of initial NO₃ to Chl-a [μM (μg L^{− 1})^{− 1}]. The ratio of initial NO₃ to Chl-a [μM (μg L^{− 1})^{− 1}] in panels (B,D) can be described as a unitless ratio by multiplying it with 14. Data were obtained from the seven experiments incubating seawater samples collected from St. Gukdong and St. Gyedong in 2019.

To develop equations for the magnitude of the Chl-a peak using the initial NO₃ or the ratio of the initial NO₃ to Chl-a, four different linear regressions between two variables were analyzed (see section “Data Analysis and Equation Development”). The Chl-a peak of the phytoplankton community was significantly correlated with the initial NO₃ concentration and the ratio of the initial NO₃ to Chl-a (Figures 6A,B). The ratio of peak-to-initial Chl-a of the phytoplankton community was also significantly correlated with the initial NO₃ and with the ratio of initial NO₃ to Chl-a (Figures 6C,D). Equations (Eqs 1–4 in Table 2) for calculating the Chl-a peak of the phytoplankton community were obtained from linear regressions in Figure 6.

TABLE 2

Table 2. Equations for calculating the peak chlorophyll-a concentration (Chl-a peak, μg L^–1) of the phytoplankton community and the day elapsed until the Chl-a peak occurred (Day_peak; day).

These types of correlations were significant when NH₄, DIN, PO₄, and SiO₂ were used (Supplementary Figures 3, 4). However, the Chl-a peak or the ratio of the peak-to-initial Chl-a of the phytoplankton community was not significantly correlated with the ratio of the initial NO₃ to PO₄ and the ratio of the initial NO₃ to SiO₂ (Supplementary Figure 5).

Correlations Between Variables Related to the Timing of Chl-a Peaks

In these incubation experiments, the elapsed day until the Chl-a peak of the phytoplankton community occurred (Day_peak) ranged from days 0 to 4, (Figures 7A,B). Furthermore, the growth rates of the phytoplankton community ranged from 0.0 to 1.4 day^–1 (Figures 7C,D).

FIGURE 7

Figure 7. (A,B) The elapsed day until peak chlorophyll-a occurred (Day_peak, day) or (C,D) growth rate of the phytoplankton community (day^{− 1}) as a function of the natural log of the initial NO₃ (Ln μM) or the natural log of the ratio of the initial NO₃ to Chl-a [Ln μM (μg L^{− 1})^{− 1}]. Data were obtained from the seven experiments incubating seawater collected from St. Gukdong and St. Gyedong in 2019.

To develop equations for determining Day_peak using the natural log of the initial NO₃ or the natural log of the ratio of initial NO₃ to Chl-a, four different linear regressions between two variables were analyzed (see section “Data Analysis and Equation Development”). Day_peak of the phytoplankton community was significantly correlated with the natural log of the initial NO₃ and with the natural log of the ratio of the initial NO₃ to Chl-a (Figures 7A,B). The growth rate of the phytoplankton community was also significantly correlated with the natural log of the initial NO₃ and with the natural log of the ratio of the initial NO₃ to Chl-a (Figures 7C,D). Equations (Eqs 5–8 in Table 2) for calculating Day_peak were obtained from the equations of linear regressions in Figure 7.

These types of correlations were significant when NH₄, DIN, and SiO₂ were used, however, this was not always the case when PO₄ was used (Supplementary Figures 6, 7). Moreover, Day_peak or growth rate of the phytoplankton community was significantly correlated with the natural log of the ratio of initial NO₃ to PO₄ but not with the ratio of the initial NO₃ to SiO₂ (Supplementary Figure 8).

Comparison of the Predicted and Field Values of Chl-a Peak and Day_peak

The ranges of the daily average values of the water temperature, salinity, NO₃, and Chl-a in St. Gwangyang and St. Masan at the selected dates in 2016–2019 were different from those of St. Gukdong and St. Gyedong (Table 3 and Supplementary Table 1). The daily average water temperature ranges at St. Gwangyang (17.9–31.2°C) and St. Masan (16.5–30.6°C) were greater than those at St. Gukdong and St. Gyedong (18.2–26.6°C). The maximum values of the daily average salinity at each station were slightly different, but the minimum values at St. Gwangyang (5.1) and St. Masan (19.1) were much lower than those at St. Gukdong and St. Gyedong (24.4). The ranges of the daily average NO₃ concentration at St. Gukdong and St. Gyedong (0.6–39.7 μM) were much greater than those at St. Gwangyang (0.6–16.1 μM) and St. Masan (0.1–7.7 μM). However, the ranges of the daily average Chl-a concentration at St. Gwangyang (1.1–49.5 μg L^–1) and St. Masan (0.3–36.8 μg L^–1) were much wider than those at St. Gukdong and St. Gyedong (0.8–24.5 μg L^–1).

TABLE 3

Table 3. Comparisons of water temperature (T, °C), salinity, NO₃ (μM), and initial chlorophyll-a (Chl-a, μg L^–1) in the seawater samples used for the incubation experiments (St. Gukdong and St. Gyedong) and the data obtained from the monitoring sites (St. Gwangyang and St. Masan) in the selected periods between 2016 and 2019 in the South Sea of Korea.

Applying Eqs 1–3 from Table 2, the actual peak Chl-a value measured at St. Gwangyang or St. Masan and the combined actual Chl-a peak values from both stations were not always significantly correlated with the predicted peak Chl-a value (Figures 8A–I). By contrast, using Eq. 4, the predicted values of Chl-a peak showed significant correlations with the actual Chl-a peak (Figures 8J–L). The root mean square errors (RMSE) were 10.6, 9.2, and 11.6 μg L^–1 for St. Gwangyang, St. Masan, and the combined data of the two stations, respectively. Thus, Eq. 4 was chosen as a best-fit equation to predict the Chl-a peak in the field.

FIGURE 8

Figure 8. Comparison between the actual and predicted values of the Chl-a peak (μg L^{− 1}) at St. Gwangyang (closed circle), St. Masan (open square), and both stations. The predicted Chl-a peak values calculated using concentrations of NO₃ and Chl-a at each station and Eq. 1 (A–C), Eq. 2 (D–F), Eq. 3 (G–I), and Eq. 4 (J–L) in Table 2.

Based on Eqs 5, 7, and 8 from Table 2, the actual Day_peak data obtained from St. Gwangyang or St. Masan and the combined data of the two stations were not always significantly correlated with the predicted Day_peak (Figures 9A–C,G–L). On the other hand, with Eq. 6, the actual Day_peak data were significantly correlated with the predicted Day_peak (Figures 9D–F). The RMSE was all 1 day for St. Gwangyang, St. Masan, and the combined data of the two stations. Thus, Eq. 6 was selected as a best-fit equation to predict Day_peak in the field.

FIGURE 9

Figure 9. Comparison between the actual and predicted values of the day elapsed until the Chl-a peak occurred (Day_peak) at St. Gwangyang (closed circle), St. Masan (open square), and both stations. The predicted Day_peak values calculated using concentrations of NO₃ and Chl-a at each station and Eq. 5 (A–C), Eq. 6 (D–F), Eq. 7 (G–I), and Eq. 8 (J–L) in Table 2. Predicted Day_peak in panels (G–L) were calculated using the predicted Chl-a peak from Eq. 4.

These comparisons were also made on Eqs 9–16 for in situ NH₄ data and Eqs 17–24 for in situ DIN data (Supplementary Tables 1–3). The actual peak Chl-a values measured at St. Gwangyang or St. Masan and the combined actual peak Chl-a values from both stations were significantly correlated with the predicted values calculated using Eq. 11, 19, and 20 (Supplementary Figures 9, 10). The predicted Day_peak values were significantly correlated with the actual values when Eqs. 13 or 14 and Eqs. 21 or 22 were used (Supplementary Figures 11, 12). Therefore, any type of nitrogen data (NO₃, NH₄, or DIN concentrations) can be implemented to predict Day_peak values at St. Gwangyang or St. Masan.

Correlations Between Variables Related to the Magnitude and Timing of Skeletonema costatum

The Chl-a peak estimated for Skeletonema costatum ranged from 0.02 to 12.2 μg L^–1 (Supplementary Figure 13A,B). The ratio of the peak-to-initial Chl-a estimated for S. costatum showed a spectrum of 1.0–90.4 (Supplementary Figures 13C,D). The ratio of the initial NO₃ to Chl-a estimated for S. costatum had values in the range 2.7–1225.6 μM (μg L^–1)^–1 (equivalent to 38 and 17,158, in a unitless ratio, respectively; Supplementary Figures 13B,D).

Among the four types of linear regression (see section “Data Analysis and Equation Development”), there were significant correlations between the Chl-a peak estimated for S. costatum against the initial NO₃ concentration and the ratio of the latter to the Chl-a concentration estimated for S. costatum against the ratio of the peak-to-initial Chl-a concentration estimated for S. costatum (Supplementary Figures 13A–D). These two correlations were always significant when NH₄, DIN, PO₄, or SiO₂ were used (Supplementary Table 4). However, the Chl-a peak or the ratio of the peak-to-initial Chl-a concentration estimated for S. costatum was not significantly correlated with the ratio of initial NO₃ to PO₄ concentrations nor the ratio of initial NO₃ to SiO₂ concentrations (Supplementary Table 4).

In these incubation experiments, it took 0–4 days to reach the Chl-a peak (Supplementary Figures 13E,F). The growth rates of S. costatum ranged from 0.0 to 2.3 day^–1 (Supplementary Figures 13G,H). Four different linear regressions between two variables were all significant when NO₃, NH₄, and DIN were applied (Supplementary Table 5 and Figures 13E–H). However, when PO₄ was used, there were no significant correlations (Supplementary Table 5). Moreover, Day_peak or the growth rate of S. costatum was not always significantly correlated with the natural log of the ratio of initial NO₃ to PO₄ or initial NO₃ to SiO₂ (Supplementary Table 5).

Discussion

In this study, we successfully obtained equations to predict the magnitude and timing of peak phytoplankton blooms by following several critical steps: measuring daily fluctuations in the concentrations of nutrients and Chl-a in seawater enclosures; performing linear regression analyses between pairs of diverse variables described in section “Data Analysis and Equation Development.” The regression analyses led to the development of the eight equations for predicting the magnitude and timing of bloom peaks that are presented in this paper. Among these, Eq. 4 for magnitude and Eq. 6 for timing were significantly correlated with the observed actual data from St. Gwangyang and St. Masan. Thus, the magnitude and timing of bloom peaks can be easily identified a few days in advance. One week may be a critical time for aquafarmers to reduce economic losses due to red tides or harmful algal blooms. Moreover, NO₃ and Chl-a in natural seawater samples can be evaluated in a few hours using a nutrient measuring instrument and a Chl-a sensor. Therefore, the method developed in this study is a quick way to predict the outbreaks of phytoplankton blooms.

The results of the present study showed that Chl-a decreased when the initial NO₃ and NH₄ concentrations were ≤2.4 and 1 μM (on May 16, July 17, and August 20), while it increased when the initial NO₃ and NH₄ concentrations were ≥3.5 and 2.1 μM (on the other dates), respectively. When NO₃ concentration was ≤2.4 μM, the diatoms C. closterium (May 16) and C. curvisetus (August 20) most dominated the phytoplankton assemblages. The half-saturation constants of C. closterium and C. curvisetus for NO₃ uptake are 0.4 and 0.6 μM, respectively (Anderson and Roels, 1981; Li et al., 2020). The half-saturation constant of C. closterium for NH₄ uptake ranged from 0.1 to 0.3 μM, which was lower than the ambient NH₄ concentration of 0.6 μM in the seawater where this diatom was dominant (May 16; Williams, 1964; Sunlu et al., 2006; Kingston, 2009). C. closterium and C. curvisetus were able to grow under these circumstances; however, the SiO₂ concentrations on May 16 (5.1 μM) and August 22 (4.4 μM) were the lowest among the seawater samples collected during the study period. This meant that deficiencies in SiO₂ concentration might be one of the reasons why these diatoms did not grow during incubation on those dates.

The dominant phototrophic protists in the seawaters used in this study, S. costatum, C. curvisetus, C. closterium, P. pungens, and T. furca, and the photosynthetic ciliate, M. rubrum, are commonly found in the seawater of many other regions (e.g., Tilstone et al., 1994; Marshall and Nesius, 1996; Zhang et al., 2015). Furthermore, they are known to often cause blooms in the study region and many other regions (e.g., Lips and Lips, 2017; Eom et al., 2021; Jeong et al., 2021). The naked ciliate Strobilidium spp. and the tintinnid ciliate Tintinnopsis spp. have also been reported to be common in the seawaters of many other regions (e.g., Dolan, 1991; Mironova et al., 2009). Thus, the results of the present study can be applied to other regions in which the dominant phototrophic and heterotrophic protists are similar to those of this study area. The method using the concentrations of nutrients and Chl-a, developed in this study, was able to predict the potential phytoplankton blooms even when microzooplankton such as Protoperidinium sp. and Tintinnopsis sp. were present in experimental seawaters. Moreover, although mesozooplankton might be removed through filtration of experimental seawaters, their grazing pressure on phytoplankton blooms has been reported to be low (Fiedler, 1982; Huntley, 1982; Stoecker and Sanders, 1985; Kim et al., 2013; Lee et al., 2017). Thus, grazing by mesozooplankton was not considered in this study.

The present study used a conversion factor of 14 to convert the unit from moles of nitrogen to gram because this factor has been most commonly used in previous studies (e.g., Baldock et al., 2004; Cosme et al., 2015). Some studies have empirically proved that 1 μM of nitrogen is equal to 1–2 μg L^–1 of chlorophyll in Scottish coastal waters and San Francisco Bay Delta (Gowen et al., 1992; Glibert et al., 2014). Thus, 0.5–1 can also be used as a conversion factor to convert the unit depending on the region or the dominant phytoplankton species.

The reported growth rates of phytoplankton assemblages in temperate coastal waters showed a spectrum of 0.0–3.4 day^–1 (Furnas, 1990 and references therein); the phytoplankton community growth rate in Loch Creran in Scotland (0.2–0.5 day^–1), Celtic Sea (0.3–1.2 day^–1), Southern California in US (0.3–1.3 day^–1), Dabob Bay in US (0.0–0.9 day^–1), and the east coast of the Izu Peninsula in Japan (0.9 day^–1) were within the range of that of the present study (0.0–1.4 day^–1). Moreover, the in situ growth rate of S. costatum in many other regions has been reported to have values between 1 and 4.1 day^–1 (Furnas, 1990 and references therein); the in situ growth rates of this species in the Great Barrier Reef in Australia (2.2 day^–1), Tampa Bay in US (1.0 day^–1), and Trondhiem Fjord in Norway (1.0 day^–1) were within the range of that in the present study (0.0–2.3 day^–1). Thus, the method for predicting the timing of a bloom peak based on the growth rate can be used in regions where the growth rate of the phytoplankton community or S. costatum are similar.

In the present study, S. costatum was chosen as a representative species to determine Chl-a peaks and Day_peak at a species level. This species was present in all seawater samples collected during the study period. Furthermore, S. costatum is the most common red-tide diatom in the coastal waters of the countries that experience red tides (Jeong et al., 2021). Finally, S. costatum can inhibit the growth of the harmful dinoflagellate Margalefidinium polykrikoides, which has caused an annual economic loss of up to USD 60 million in Korea (Lim et al., 2014, 2015). Thus, if S. costatum blooms can be predicted accurately, the population dynamics of M. polykrikoides can also be determined. Taking this point into consideration, the initial nutrient-to-Chl-a ratio and the peak-to-initial Chl-a ratio of each species should be determined at the species level to develop methods for predicting blooms for each species. The range of the ratio of the initial NO₃ to Chl-a estimated for S. costatum [2.7–1,226 μM (μg L^–1)^–1] was much wider than that of the total phytoplankton community [0.1–24.2 μM (μg L^–1)^–1]. This is because S. costatum was part of the phytoplankton community and diatoms such as S. costatum had lower concentrations of Chl-a than other types of phytoplankton (Hitchcock, 1982; Moal et al., 1987). Moreover, the range of peak-to-initial Chl-a ratio of S. costatum (1.0–90.4) was much wider than that of the phytoplankton community (1.0–61.4). This species grows faster than the other phytoplankton species that make up the phytoplankton community (Finenko and Krupatkina-Akinina, 1974).

The seawater samples for experiments collected from St. Gukdong and St. Gyedong had slightly different water temperature ranges than those at the validation sites (St. Gwangyang and St. Masan). However, the ranges of other environmental factors—such as salinity and Chl-a—were largely different between the sites for the experiments (St. Gukdong or St. Gyedong) and validation (St. Gwangyang or St. Masan). Moreover, interannual variations in diverse factors (e.g., salinity, temperature, mixing, or turbidity) were not considered in this study. Nevertheless, predictions of the magnitude and timing of the peaks of the blooms at the validation sites corresponded well with the actual measurements from May to October, for each year from 2016 to 2019. Thus, the method developed in this study can be used in other regions and periods having different ranges of environmental factors from those of the study sites.

Other types of nutrients (e.g., NH₄, DIN, PO₄, or SiO₂) can also be used to predict peak bloom periods using the method described in this study, depending on the characteristics of the site of interest and the dominant bloom taxa. Some urbanized coastal and estuarine regions (e.g., Delaware Estuary, Neuse Estuary, San Fransico Estuary, Manila Bay, Scheldt Bay, Rhine Bay) or periods in a region, in particular, show high NH₄ concentrations (Middelburg and Nieuwenhuize, 2000; Burkholder et al., 2006; Yoshiyama and Sharp, 2006; Chang et al., 2009; Parker et al., 2012a; Seeyave et al., 2013; Glibert et al., 2016). NH₄ also contributes more to the growth of phytoplankton communities than NO₃ in some regions (e.g., Dokai Bay, Santa Monica Bay, and Kaneohe Bay; Harvey and Caperon, 1976; Eppley et al., 1979; Tada et al., 2009), and many phytoplankton taxa that cause noxious blooms are associated with the reduced forms of nitrogen, including NH₄ (Bates et al., 1993; Leong et al., 2004; Trainer et al., 2007; Seeyave et al., 2009; Killberg-Thoreson et al., 2014). In these cases, NH₄ should be used as the main nutrient component to predict bloom peaks. The present study mainly used equations that were based on NO₃ because DIN in seawater samples from the study area was dominated by NO₃. Moreover, NO₃ is a good indicator for predicting the magnitude and timing of algal bloom peaks for regions where diatoms are one of the most common species in regions like the sampling sites of this study because these species are preferentially NO₃-users (Lomas and Glibert, 1999a, b, 2000; Kudela and Dugdale, 2000; Glibert et al., 2014, 2016).

Enclosure experiments have enabled the investigation of phytoplankton responses to environmental factors and the calculation of growth or nutrient uptake rates of phytoplankton, preventing physical mixing with surrounding waters (Pitcher et al., 1993; Kudela and Dugdale, 2000; Parker et al., 2012b; Lee et al., 2019b; Ferreira et al., 2020). However, results from the enclosure experiments are limited by differences in chemical (e.g., nutrients, dissolved gases), physical (e.g., wind-induced turbulence, light intensity), and biological variables (e.g., predators and prey) from nature, the so-called “bottle effects” (Veldhuis and Timmermans, 2007; Nogueira et al., 2014). The predicted Day_peak value may be lower than the actual Day_peak value because equations developed in this study may be affected by the absence of mixture or dilution in nature. Comparisons between the predicted and actual Day_peak values showed one-day differences at St. Gwangyang and St. Masan. Thus, a one-day time lag should be considered when using this equation (Eq. 6). Residence time is another important parameter for determining the timing of bloom outbreaks in nature, and it varies among regions (Delesalle and Sournia, 1992; De Sève, 1993; Phlips et al., 2012; Chen et al., 2013; Wan et al., 2013; Ralston et al., 2015; Glibert et al., 2018a). The residence time in Gamak Bay is approximately 11 days (Kim et al., 2016); this is greater than the highest Day_peak of 4 days obtained in the present study. Thus, residence time might not largely affect the equations derived from the results of the enclosure experiments in this study.

Models for predicting phytoplankton blooms have been developed using technological advances and knowledge accumulated over decades (McGillicuddy, 2010; Anderson et al., 2015; Franks, 2018; Zohdi and Abbaspour, 2019). However, there are still many challenges with this due to the complexity of the models involved and the difficulty of validating them; therefore, works combining laboratory and field studies are needed (Flynn and McGillicuddy, 2018). The present study combined laboratory and field work to improve the prediction of the magnitude and timing of phytoplankton blooms. Moreover, the main method used in this study was a simple linear regression analysis that has been successfully used to predict the magnitude of phytoplankton blooms in aquatic environments for decades (e.g., Blum et al., 2006; Phillips et al., 2008; Kahru et al., 2020). These simple models may help understand and interpret ecosystem dynamics, thereby demonstrating the relationship between two ecological variables (Glibert et al., 2010; Anderson et al., 2015). Thus, the equations developed in the present study only need two variables: the concentrations of Chl-a and the main nutrient source in a target region. This means that it can be used easily by those in the aquaculture industry to predict daily bloom peaks.

In summary, this study suggests the following method for predicting phytoplankton bloom dynamics. First, seawater from a target region should be incubated for ten days in a closed system; the concentrations of the main nutrients and Chl-a should be measured daily. Second, the linear regressions described in Figures 6, 7 should be used to analyze the data collected. Third, in situ Chl-a and nutrient concentrations have to be introduced into the linear regression equations to allow calculation of the predicted values for Chl-a peak and Day_peak. By demonstrating how this can be done, the present study could provide a basis for understanding phytoplankton bloom dynamics via the use of modeling to predict harmful algal blooms in various regions worldwide.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

JO and HJ designed the study conception and drafted the manuscript. JO, HJ, JY, HK, SP, AL, SL, and SE obtained the data, conducted the experiments, and approved the final submitted manuscript. JO, HJ, and JY performed data analyses. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Useful Dinoflagellate program of Korea Institute of Marine Science and Technology Promotion (KIMST) funded by the Ministry of Oceans and Fisheries (MOF) and the National Research Foundation (NRF) funded by the Ministry of Science and ICT (NRF-2017R1E1A1A01074419; NRF-2020M3F6A1110582) award to HJ.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a past co-authorship with one of the authors HJ.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to Eun Chong Park for technical support. We also thank the editor and reviewers for their valuable comments.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.681252/full#supplementary-material

Abbreviations

Day_peak, the time elapsed until the Chl-a peak occurred.

Footnotes

References

Anderson, C. R., Moore, S. K., Tomlinson, M. C., Silke, J., and Cusack, C. K. (2015). “Living with harmful algal blooms in a changing world: strategies for modeling and mitigating their effects in coastal marine ecosystems,” in Coastal and Marine Hazards, Risks, and Disasters, eds J. F. Shroder, J. T. Ellis, and D. J. Sherman (Amsterdam: Elsevier), 495–561. doi: 10.1016/B978-0-12-396483-0.00017-0

ORIGINAL RESEARCH article

Phytoplankton Bloom Dynamics in Incubated Natural Seawater: Predicting Bloom Magnitude and Timing

Introduction

Materials and Methods

Sampling Site Description and Water Sampling

Experimental Procedures

Subsampling and Component Analyses

Data Analysis and Equation Development

Comparison of the Predicted and Field Data

Statistical Analysis

Results

Physical and Chemical Properties of the Seawater Samples

Dominant Protist Species at the Beginning of Incubation

Initial Nutrient and Chl-a Concentrations

Variations in Nutrient and Chl-a During Incubation

Correlations Between the Variables Related to the Chl-a Peak

Correlations Between Variables Related to the Timing of Chl-a Peaks

Comparison of the Predicted and Field Values of Chl-a Peak and Daypeak

Correlations Between Variables Related to the Magnitude and Timing of Skeletonema costatum

Discussion

Data Availability Statement

Author Contributions

Funding

Conflict of Interest

Publisher’s Note

Acknowledgments

Supplementary Material

Abbreviations

Footnotes

References

Comparison of the Predicted and Field Values of Chl-a Peak and Day_peak