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

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 NO3 or the ratio of the initial NO3 to chlorophyll-a. To predict the timing, the chlorophyll-a peak timing or the growth rate against the natural log of NO3 or the natural log of the ratio of the initial NO3 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; Abbreviations: Day peak , the time elapsed until the Chl-a peak occurred. 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) 1 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.

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 . 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) 2 and the Korea Hydrographic and Oceanographic Agency (2021) 3 .

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 temperaturecontrolled 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 , >58 µmol photons m −2 s −1 for the dinoflagellate Alexandrium pohangense , >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 , 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) 4 . 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 3 ), ammonium (NH 4 ), phosphate (PO 4 ), and silicate (SiO 2 ) in the seawater samples were measured using a 4-channel nutrient auto-analyzer (QuAAtro, SEAL Analytical GmbH, Norderstadt, Germany; Bran+Luebbe Analyzing Technologies, 2004aTechnologies, ,b, 2005SEAL 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 3 , NH 4 , PO 4 , SiO 2 , and Chl-a were measured every sampling day and dissolved inorganic nitrogen (NO 3 +NH 4 ; 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 3 to Chl-a [µM (µg L −1 ) −1 ] was calculated by dividing the initial NO 3 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: µ = Ln 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 3 ; the Chl-a peak versus the ratio of the initial NO 3 to Chl-a; the ratio of the peak-to-initial Chl-a versus the initial NO 3 ; and the ratio of the peak-to-initial Chl-a versus the ratio of the initial NO 3 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 3 ; Day peak versus the natural log of the ratio of the initial NO 3 to Chl-a; the growth rate versus the natural log of the initial NO 3 ; and the growth rate versus the natural log of the ratio of the initial NO 3 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 4 , DIN, PO 4 and SiO 2 , as well as the ratio of the initial NO 3 to PO 4 and the ratio of the initial NO 3 to SiO 2 . 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) 5 . 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 3 , Chl-a, salinity, and water temperature measured at fiveminute 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 3 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 3 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.

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).
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).

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 3 , NH 4 , DIN, PO 4 , or SiO 2 became depleted (Figures 5D-G).
Chlorophyll-a decreased during incubation when the initial NO 3 concentrations were ≤2.4 µM (Figures 5A-C,H-J and Supplementary Figure 2A). In contrast, Chl-a increased when the initial NO 3 concentrations were ≥3.5 µM (Figures 5D-G,K-N and Supplementary Figure 2B). Chl-a decreased when the initial NH 4 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 4 concentrations were ≤0.3 µM (Figures 5A-C,AC-AE) with one exception: Chla increased when the initial PO 4 concentration was 0.1 µM in the seawater on July 24 (Figures 5F,AH). Chl-a also decreased when the initial SiO 2 concentrations were ≤24.5 µM (Figures 5A-C,AJ-AL) with one exception: Chl-a increased when the initial SiO 2 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-toinitial Chl-a of the phytoplankton community varied from 1.0 to 61.4 (Figures 6C,D). The ratio of the initial NO 3 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).
To develop equations for the magnitude of the Chl-a peak using the initial NO 3 or the ratio of the initial NO 3 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 3 concentration and the ratio of the initial NO 3 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 3 and with the ratio of initial NO 3 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.
These types of correlations were significant when NH 4 , DIN, PO 4 , and SiO 2 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 3 to PO 4 and the ratio of the initial NO 3 to SiO 2 (Supplementary Figure 5).

Correlations Between Variables Related to the Timing of Chl-a Peaks
In these incubation experiments, the elapsed day until the Chla 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).
To develop equations for determining Day peak using the natural log of the initial NO 3 or the natural log of the ratio of initial NO 3 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 3 and with the natural log of the ratio of the initial NO 3 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 3 and with the natural log of the ratio of the initial NO 3 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 4 , DIN, and SiO 2 were used, however, this was not always the case when PO 4 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 3 to PO 4 but not with the ratio of the initial NO 3 to SiO 2 (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 3 , 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 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.
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.
These comparisons were also made on Eqs 9-16 for in situ NH 4 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 3 , NH 4 , 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 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). 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 3 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).

Eq. Equation
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 3 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 4 , DIN, PO 4 , or SiO 2 were used (Supplementary Table 4). However, the Chla 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 3 to PO 4 concentrations nor the ratio of initial NO 3 to SiO 2 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 3 , NH 4 , and DIN were applied (Supplementary Table 5 and Figures 13E-H). However, when PO 4 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 3 to PO 4 or initial NO 3 to SiO 2 (Supplementary Table 5).  Table 2.

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 3 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 3 and NH 4 concentrations were ≤2.4 and 1 µM (on May 16, July 17, and August 20), while it increased when the initial NO 3 and NH 4 concentrations were ≥3.5 and 2.1 µM (on the other dates), respectively. When NO 3 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 3 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 4 uptake ranged from 0.1 to 0.3 µM, which was lower than the ambient NH 4 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; Masan (open square), and both stations. The predicted Day peak values calculated using concentrations of NO 3 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. however, the SiO 2 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 2 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 Chla, 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 . 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(Lim et al., , 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 3 to Chla 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 4 , DIN, PO 4 , or SiO 2 ) 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 4 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 4 also contributes more to the growth of phytoplankton communities than NO 3 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 4 (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 4 should be used as the main nutrient component to predict bloom peaks. The present study mainly used equations that were based on NO 3 because DIN in seawater samples from the study area was dominated by NO 3 . Moreover, NO 3 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 3 -users Glibert, 1999a,b, 2000;Kudela and Dugdale, 2000;Glibert et al., 2014Glibert et al., , 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., windinduced 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 oneday 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.