Sampling for this study took place on August 31 in Yeosu, South Sea of Korea, specifically in the northern region of Bodolbada, located between Yeosu and Goheung (Fig. 1). The survey was conducted using a vessel operated by Yeosu City for Margalefidinium red-tide control projects. At each sampling station, water temperature and salinity were recorded using a multi-parameter analyzer (Orion Star A329; Thermo Fisher Scientific, Waltham, MA, USA), and water samples for chlorophyll-a (Chl-a), particulate organic carbon (POC), and nutrient analysis were collected using bucket sampling. Raw water samples were initially pre-filtered through a 200 μm nylon sieve to remove larger plankton that could potentially interfere with subsequent analyses, before being aliquoted for biochemical analyses.

For Chl-a sample collection, 200 mL of seawater was filtered through a 0.7 μm filter (47 mm, Whatman GF/F; Cytiva, Buckinghamshire, UK) under gentle vacuum pressure (<100 mmHg). Similarly, 500 mL of water samples for POC and nutrient analyses were filtered through 0.7 μm GF/F filters, with the filter residue used for POC analysis and the 15 mL of filtrate for nutrient analysis. All samples were frozen at −20°C until further analysis.

To characterize protist communities, up to 1 L of water was collected using a bucket sampling at each sampling and filtered through a 0.45 μm cellulose nitrate membrane filter (47 mm; Whatman, Cytiva, Germany) for 18S rRNA gene sequencing. Additionally, a 0.6 L sample was preserved in Lugol’s iodine solution (final concentration of 2%) for cell counts. For metazoans, samples were collected using double 300 μm-mesh Bongo nets (60 cm mouth diameter) equipped with filtering cod ends through horizontal tows in surface waters. Samples were preserved in 5% formalin in seawater. The volume of water filtered during each tow was measured using a flow meter mounted at the center of the net mouth (model 438 115; Hydro-Bios, Kiel, Germany).

For the Chl-a analysis, 10 mL of 90% acetone was added to the tube containing 500 mL of filtered sample, followed by 10 min of sonification. The samples were then placed in a dark chamber and kept at 4°C overnight. The following day, the samples were centrifuged, and the supernatant was utilized for Chl-a measurement using a Trilogy Laboratory Fluorometer (Turner Designs, San Jose, CA, USA). Details are described by Holm-Hansen and Riemann (1978).

Filter paper samples for POC analysis were dried at 50°C for 12 h to eliminate moisture. Subsequently, acid fumigation with HCl was carried out in a desiccator for 12 h to remove inorganic carbon, and POC was measured using an element analyzer (Flash 2000; Thermo Fisher Scientific). Detailed methods are described by Choi et al. (2024).

For nutrient analysis, concentrations of nitrate + nitrite (NO₃ + NO₂), ammonium (NH₄), phosphate (PO₄), and silicic acid (Si[OH]₄) were quantified using a nutrient auto-analyzer (New QuAAtro39; SEAL Analytical, Southampton, UK). Quality control was ensured by using reference standards for each nutrient.

A plankton sample was collected from bloom waters of Yeosu, Korea, in August 2022. The sample was passed through a 154-μm Nitex mesh screen and subsequently placed in 6-well tissue culture plates. A clonal culture was established through two successive single-cell isolations. As the culture volume increased, cells were transferred sequentially to 32-, 250-, and 500-mL polycarbonate bottles. The bottles were placed on a shelf at 25°C under an illumination of 20 μE m⁻² s⁻¹ and a 14 : 10-h LD cycle. Detailed methods of isolation and culturing are described by Jang et al. (2017).

For metabarcoding analyses, DNA was extracted from the filters using the DNeasy PowerSoil Kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. The extracted DNA was quantified with a Qubit fluorometer using Quant-iT PicoGreen (Invitrogen, Carlsbad, CA, USA) and stored at −80°C until use in polymerase chain reactions (PCRs).

Sequencing libraries were prepared following the Illumina metagenomic sequencing library protocol for gene amplification (San Diego, CA, USA). Two nanograms of genomic DNA were PCR-amplified using the universal protist primer pair TAReuk454FWD1 (5′-CCAGCASCYGCGGTAATTCC-3′) and V4 18S Next.Rev (5′-ACTTTCGTTCTTGATYRATGA-3′), at 500 nM each, targeting the V4 region of the 18S rRNA gene (Stoeck et al. 2010, Piredda et al. 2017). Amplification was carried out with Herculase II fusion DNA polymerase (Agilent Technologies, Santa Clara, CA, USA) in a 5× reaction buffer with a 1 mM dNTP mix. The thermocycler protocol for the initial PCR included 3 min at 95°C, followed by 25 cycles of 30 s at 98°C for denaturation, 60 s at 65°C for annealing, and 90 s at 72°C for extension. The resulting products were purified using AMPure beads (Agencourt Bioscience, Beverly, MA, USA). Following purification, 2 μL of each initial PCR product underwent a second amplification for final library construction using NexteraXT Indexed Primers. The thermocycler conditions for the second PCR were 3 min at 95°C; 10 cycles of 30 s at 95°C, 30 s at 55°C, and 30 s at 72°C; followed by a final 5-min extension at 72°C. After purification using AMPure beads, the purified PCR product was quantified via quantitative PCR using KAPA Library Quantification kits for Illumina sequencing platforms (Roche, Basel, Switzerland) according to the manufacturer’s protocol and assessed using a TapeStation D1000 ScreenTape system (Agilent Technologies, Waldbronn, Germany). Paired-end (2 × 250 bp) sequencing was conducted at Macrogen (Seoul, Korea) using a single lane on the Illumina MiSeq platform (Illumina, San Diego, CA, USA).

Sequencing data were demultiplexed using QIIME 1.9.1 and Cutadapt 3.2 and then analyzed via the QIIME 1.9 pipeline (Martin 2011, Bolyen et al. 2019). The reads were trimmed, assembled, and subjected to quality control. Denoising, pairing, and chimera removal were performed on the demultiplexed paired-end reads using the DADA2 software (Callahan et al. 2016). The processed sequences were then clustered into amplicon sequence variants (ASVs). To address potential sampling depth biases, the sequencing pool was randomly subsampled 10 times using the “multiple_rarefactions_even_depth.py” script in QIIME. The resulting ASV tables were consolidated for further analysis. Taxonomic assignment of ASVs was conducted using the PR2 4.14.0 database with the default parameters (Guillou et al. 2012), and ASVs identified as metazoans were excluded from the analysis.

To further identify strains of M. polykrikoides, which dominated the HAB in Yeosu, single cells from samples fixed in Lugol’s solution and stored at 4°C were isolated in 0.2 mL PCR tubes containing 15 μL of Milli-Q water. The tubes were then frozen at −80°C for 3 min to rupture the cell membranes and subsequently thawed. Amplification was performed using 5 μL of 10× Ex Taq buffer, 4 μL of 2.5 mM dNTP mix, 0.25 μL of 5 U μL⁻¹ ExTaq polymerase (RR006A; Takara Bio, Madison, WI, USA), and 0.2 μM of primers for the 28S rRNA gene region: D1RF (Scholin et al. 1994) and LSUB (Litaker et al. 2003). The PCR thermocycler protocol consisted of an initial denaturation at 95°C for 5 min, followed by 38 cycles of 95°C for 40 s, 59°C for 2 min, and 72°C for 1 min, with a final extension step at 72°C for 5 min. PCR products showing a clear single band were then purified using the ExoSAP-IT PCR Purification Kit (Qiagen). The sequences from Yeosu’s Margalefidinium, along with multiple sequences obtained from the NCBI nucleotide database, were aligned using the native implementation of ClustalW in MEGA v.4.0 (Tamura et al. 2007) and manually refined. Phylogenetic relationships were then inferred using the maximum-likelihood algorithm in the RAxML 7.0.3 program (Stamatakis 2006) and the Bayesian method in MrBayes v.3.1 (Ronquist and Huelsenbeck 2003). Detailed methods are described by Jang et al. (2018).

Protists in the collected samples were counted and identified at the species level. For microscopic enumeration, 600 mL subsamples preserved in Lugol’s solution were concentrated 5 to 10 times using a 2-d settling method. After thorough mixing, either all or a minimum of 300 cells from three to five 1 mL Sedgwick–Rafter counting chambers were counted at 100× and/or 200× magnification using a light microscope (Nikon, Tokyo, Japan).

Metazoans in the net samples were counted and identified down to the species level. Samples were subsampled using a Folsom-type splitter, with ratios ranging from 1/64–1/2. Species identification and quantification were performed under a stereomicroscope (Nikon SMZ645) within a Bogorov counting chamber, and further detailed identification was conducted using a high-magnification optical microscope (Nikon ECLIPSE 80i). Metazoan abundance was expressed as the number of individuals per cubic meter (individuals m⁻³).

The M. polykrikoides bloom was investigated on Aug 31, 2022. On the days leading up to the survey date, the average daily sea surface temperatures in the study area showed a continuous decline, reaching approximately 23°C on Aug 31 (Fig. 2A). Salinity ranged from 31.5 to 33 psu before and after the survey period, with more than 20 mm of rainfall beginning the day before the survey and continuing into the day of the survey (Fig. 2B). The wind speed remained steady at 6–9 m s⁻¹ up until the day of the survey, after which it recorded an upward trend, reaching 11–16 m s⁻¹ (Fig. 2C). Margalefidinium polykrikoides cell densities continuously increased, reaching a maximum of 3,550 cells mL⁻¹ on the day of the survey, followed by a declining trend thereafter (Fig. 2D).

At the high bloom-intensity station (H-BS), the HAB was visually distinct, with the water appearing darker and more turbid compared to the clearer waters at the adjacent, non-bloom station (N-BS) (Fig. 3A & B). The recorded maximum density of M. polykrikoides at H-BS, 2,344 cells mL⁻¹, was significantly higher than that at N-BS, only 30 cells mL⁻¹ (Table 1). These differences indicate the spatial heterogeneity of the bloom intensity across the study sites. Additionally, the total Chl-a concentrations at H-BS were over 10 times higher than those at N-BS (t-test, p < 0.01) (Fig. 3C), which is consistent with the bloom’s cell density. However, POC and NO₃ + NO₂, PO₄, and Si(OH)₄, showed relatively minor variations between the stations, while NH₄ levels were elevated at H-BS (Fig. 3D & E).

The intense bloom, predominantly composed of the dinoflagellate M. polykrikoides, was subjected to further genetic analysis by sequencing the large subunit (LSU) region of the rRNA gene from individual cells (Fig. 4). Sequence alignment revealed that the LSU rRNA of the M. polykrikoides present at Yeosu was identical to certain known strains yet showed approximately 10% divergence from others. The phylogenetic analysis indicated that our strain grouped within a major clade of the East Asian ribotype, but was distinctly separated from the strain originating from the Mediterranean, as well as the American/Malaysian and Philippine ribotypes (Fig. 4).

A total of 30 protist taxa were observed through cell counts in the surveyed region (Table 1). Of these, 15 species were detected at H-BS and 23 species at N-BS, with seven species being unique to H-BS and 15 to N-BS. In terms of total cell abundance, excluding the bloom-forming species M. polykrikoides (2,344 cells mL⁻¹ at H-BS), the overall protist cell abundance was lower at H-BS than at N-BS. However, certain species, such as the diatoms Chaetoceros decipiens and Pseudo-nitzschia pungens, the dinoflagellate Alexandrium sp., and the naked ciliate Mesodinium rubrum, exhibited remarkably higher abundances at H-BS than at N-BS (Table 1). Moreover, the protozoan dinoflagellate Noctiluca scintillans, which dominated the meso-sized plankton group, exhibited nearly double the biomass at H-BS (Table 2).

The environmental DNA revealed more about the species composition and relative proportions of the two protist communities (Fig. 5, Supplementary Table S1). A total of nine species were uniquely identified at H-BS, while a significantly higher number, 152, was found at N-BS (Fig. 5A). At H-BS, M. polykrikoides dominated the community in terms of relative abundance, whereas at N-BS, the Syndiniales group occupied a notable portion as the second most abundant group, followed by diatoms and cryptophytes (Fig. 5B). Within the Dinoflagellata group, Alexandrium affine showed significant relative abundance at H-BS, while Noctiluca scintillans, Gyrodinium spp., and A. affine were all prominent at N-BS (Fig. 5C). The relative abundances of taxonomic groups within Syndiniales were similar between the two stations. Additionally, within the diatom group, Chaetoceros spp., primarily comprised of C. affinis, showed the highest relative abundance in both communities, followed by Coscinodiscus spp. Of the ciliate community, the non-loricate group Strombidiida was the most dominant, followed by Choreotrichida and the loricate group Tintinnida, with no significant differences in composition between the stations.

A total of 14 metazoan taxa, along with the meso-sized dinoflagellate N. scintillans, were observed through microscopic examination in the surveyed region (Table 2). All species were detected at H-BS, while only five species were found at N-BS. Notably, the relative abundances of the Cladoceran species Evadne tergestina, the copepods Acartia ohtsukai and Tortanus forcipatus, and the Chaetognathan Sagitta crassa were significantly higher at H-BS compared to N-BS. Overall, the total abundance of metazoans was higher at H-BS, where the M. polykrikoides bloom occurred.

The environmental parameters during the M. polykrikoides bloom highlight several key factors that may have played a role in promoting the bloom’s intensity. First, the observed sea surface temperature of 23°C and salinity of around 32 psu aligns with previous studies that identified the maximum growth ranges of M. polykrikoides (Supplementary Table S2). While the temperature and salinity ranges for bloom occurrences in the southern coastal waters of Korea tend to be broader, they generally focus around 25°C and the low 30s psu range, respectively (Supplementary Table S3).

While this study did not measure nutrient concentrations before the red tide event, presenting limitations in interpreting nutrient availability, rainfall occurring before and during the survey likely facilitated nutrient influx from terrestrial runoff. Indeed, the concentrations of nitrate and phosphate measured on the survey day exceed the half-saturation constant (Ks) values for M. polykrikoides (Ks = 2.10 and 0.57 μM, respectively) reported by Kim et al. (2001). This suggests, albeit indirectly, that the nutrient availability may have partially supported the bloom development process (Gobler et al. 2012, Al-Azri et al. 2014). Furthermore, assuming the two surveyed stations are within the same water mass, the elevated ammonium concentrations at the H-BS station are likely driven by secondary processes, including the decomposition of organic matter resulting from the bloom (Lim et al. 2021).

Another significant observation was the area’s wind speeds. Before the bloom, the wind speed remained moderate at 6–9 m s⁻¹, providing stable conditions conducive to the accumulation of surface waters and nutrients. However, as wind speeds increased to 11–16 m s⁻¹ following the peak of the bloom, it is possible that enhanced mixing and turbulent water conditions may have contributed to the subsequent decline in bloom density. Previous studies, such as Lim et al. (2015a), suggest that maximum wind speeds exceeding 14 m s⁻¹ could disrupt HABs by promoting vertical mixing and the dispersal of bloom-forming species.

The spatial heterogeneity observed between the high bloom-intensity H-BS and non-bloom N-BS stations also underscores a physical effect, namely wind-driven transport. Despite both stations being within the same general water mass, the pronounced difference in bloom intensity suggests that local hydrodynamic conditions may have contributed to the aggregation and retention of M. polykrikoides at H-BS, while they dispersed cells at N-BS.

Furthermore, the Chl-a concentrations, which were over 10 times higher at H-BS than at N-BS, confirm the presence of a dense algal community at the bloom site. Interestingly, despite the sharp contrast in Chl-a, the POC levels between the two stations showed minimal differences. This suggests that the contribution of the M. polykrikoides bloom to the overall POC within the community was low at the time of observation. This could be attributed to factors such as the contribution of bacterial biomass, which was not examined in this study, or the relatively low bloom-derived contribution due to the high background organic carbon levels characteristic of coastal areas.

Overall, the environmental factors driving the M. polykrikoides bloom appear to be a combination of optimal temperature, nutrient enrichment from rainfall, and localized wind conditions that initially facilitated bloom formation but later contributed to its dissipation by increasing water column mixing. These findings emphasize the intricate balance of physical and chemical processes that influence HAB dynamics in coastal regions.

The genetic analysis of the M. polykrikoides strain present in the Yeosu area, particularly the high divergence of its LSU rRNA sequence from other strains, reaching 10% dissimilarity, suggests potential adaptation to local environmental conditions. This bloom-forming strain clustered within the East Asian ribotype, genetically distinct from those in regions such as the Philippines and America. These findings indicate that local environmental conditions, such as regional nutrient availability and temperature and salinity fluctuations, could drive the evolutionary divergence of local variants (Moore et al. 1998, Iwataki et al. 2008). This highlights the importance of considering regional genetic diversity when predicting potential future blooms and their environmental impacts.

Generally, the presence of cysts, which are known to contribute to the persistence and sudden outbreak of dinoflagellate blooms, can be a key factor in bloom dynamics (Tang and Gobler 2012, Shin et al. 2017). It is plausible that our strain of M. polykrikoides may have remained in a dormant cyst form in the sediments until favorable conditions occurred, as reported in other studies in the region (Li et al. 2015, Park et al. 2016).

Moreover, the observed genetic divergence raises questions about the strain’s competitive ability and allelopathic potential, particularly in terms of how it interacts with other species within the local ecosystem. Local adaptation may influence the strain’s capacity to outcompete other protists, thereby shaping the bloom’s dynamics and its broader ecological impacts. These insights highlight the critical need for the ongoing monitoring of genetic variability within HAB species, especially in regions prone to frequent bloom events (Lebret et al. 2012).

The significant differences in protist species richness and community composition observed between the two stations highlight key ecological dynamics that drive bloom formation and its effects on community structure. Blooms often lead to reduced diversity as bloom-forming species monopolize available resources, such as nutrients and light, creating less favorable conditions for other taxa (Koch et al. 2014, Jeong et al. 2015). The notable reduction in the overall protist abundance at H-BS further supports the hypothesis of resource depletion or potential allelopathic effects induced by M. polykrikoides (Tang and Gobler 2010, Band-Schmidt et al. 2020). This competitive strategy may allow M. polykrikoides to dominate the community, thereby reducing the abundance and diversity of other protists in its proximity.

Interestingly, certain species, such as C. affinis, C. decipiens, and P. pungens, exhibited higher abundances at H-BS, suggesting they may possess unique characteristics that enable them to persist or even thrive under bloom conditions. These diatoms are known for their resilience to nutrient stress and chain-length plasticity, which may enhance their resistance to grazing (Eppley et al. 1969, Bergkvist et al. 2012). Also, some of these diatoms are specifically recognized for inhibiting the growth of M. polykrikoides through allelopathy (Lim et al. 2014). Their success in the bloom-dominated environment indicates that some diatoms can coexist with bloom-forming species, at least under specific conditions, possibly due to their ability to utilize different nutrient sources or occupy distinct ecological niches.

The presence of the meso-sized dinoflagellate N. scintillans, with H-BS having nearly double the biomass of N-BS, is particularly noteworthy. N. scintillans is a known grazer of smaller protists and organic particles, and its elevated abundance suggests that it may benefit from the high concentrations of particulate organic matter and smaller prey species associated with blooms (Kiørboe and Titelman 1998, Kiørboe 2003). While direct evidence of N. scintillans feeding on M. polykrikoides is lacking, the frequent co-occurrence of N. scintillans in other reported M. polykrikoides bloom communities implies a potential for a predator-prey interaction between the two species (Lim et al. 2017, Xiaodong et al. 2023). Similarly, the significant abundances of dinoflagellate Gyrodinium species within the community suggest their potential role as predators (Kang et al. 2020).

The dominance of Alexandrium species at H-BS hints at possible competitive dynamics with M. polykrikoides. Although M. polykrikoides was the primary bloom-former, the high relative abundance of Alexandrium suggests that it may also play a significant role in the bloom ecology, as a competitor resilient to bloom-related stress. For instance, Lim et al. (2019) suggested that Alexandrium species exhibit a relative competitive advantage over M. polykrikoides under higher temperatures and lower salinities. Furthermore, certain Alexandrium species has been reported to gain a competitive edge by directly feeding on M. polykrikoides (Lim et al. 2015b). Moreover, frequent reports of Alexandrium blooms either occurring in proximity to or co-occurring with M. polykrikoides blooms in time series data (Jeong et al. 2017, Lim et al. 2019, Wolny et al. 2020) highlight the need for further research into their interactions, which could provide crucial insights into the mechanisms driving bloom dynamics.

There is limited direct research specifically focusing on the relationship between Syndiniales and M. polykrikoides, but some studies have explored potential parasitic interactions involving Syndiniales in HABs (Hattenrath-Lehmann et al. 2019, Cui et al. 2020). Syndiniales, particularly species from Dino-Group II, are widely recognized for parasitizing a variety of protists, with dinoflagellates being the primary hosts (Sehein et al. 2022). Given that Dino-Group II species were collectively the highest relative abundance group within the Syndiniales during this bloom event, further research into the parasitic interactions of Syndiniales could provide critical insights into the mechanisms governing bloom regulation.

The higher species richness and abundance of metazoans at H-BS suggest that these organisms may possess greater resistance to the adverse effects typically associated with HABs, such as the production of toxins or the depletion of oxygen. Several studies have indicated that certain metazoans, such as copepods and cladocerans, display a higher tolerance than protists to the toxic conditions induced by dinoflagellate blooms (Marcoval et al. 2013, Turner 2014). The remarkably higher relative abundances of species such as E. tergestina (Cladocera), A. ohtsukai (Copepoda), and S. crassa (Chaetognatha) at H-BS support the hypothesis that these species may have adaptations enabling them to cope with or even thrive under bloom conditions. For instance, Jiang et al. (2011) proposed that copepod adaptation to M. polykrikoides could result from the selection of resistant genotypes during a persistent bloom event. Additionally, prior research suggesting that metazoans can tolerate HABs due to physiological mechanisms that mitigate the impacts of toxins or hypoxic conditions (Turner 2014, Yang et al. 2024) also partially explains their abundance within bloom conditions.

Furthermore, the increased abundance of metazoans in the bloom-dominated waters at H-BS suggests potential trophic interactions between these species and M. polykrikoides. It is plausible that certain metazoans, such as Acartia spp. and E. tergestina, exhibit positive associations with the bloom due to the availability of prey that thrive in the bloom environment (Jeong et al. 2008, Kim et al. 2013, Lee et al. 2017b, Song et al. 2024). While the precise relationship between M. polykrikoides and metazoans remains largely unclear, the potential for these species to feed on other organisms within the bloom suggests a complex web of interactions. Further investigation into the specific feeding behaviors of these metazoans during bloom events could provide valuable insights into the ecological roles they play in controlling bloom dynamics or mitigating bloom intensity. This study provides critical insights into the ecological dynamics of M. polykrikoides blooms in the southern coastal waters of Korea. Our comparison between high-intensity-bloom and adjacent non-bloom waters highlights key differences in environmental factors, community structure, and species interactions. The bloom site (H-BS) showed a higher abundance and diversity of metazoans, likely due to their resilience against the stressors typically associated with HABs. In contrast, the non-bloom site (N-BS) exhibited greater protist species richness, suggesting that M. polykrikoides may outcompete or inhibit other protists during blooms, potentially through resource competition or allelopathy. Thus, the distinct community structures observed at H-BS and N-BS reflect the complex interplay of competition, predation, and parasitism that defines bloom dynamics. This study shows that gaining a deeper understanding of the interactions between M. polykrikoides and co-occurring planktonic species is essential for effectively predicting and mitigating future HABs in coastal ecosystems. Ongoing monitoring and investigation into species-specific adaptations, competition, and trophic dynamics will be critical for minimizing the ecological and economic consequences of blooms.