Research Article
Research Article
Monitoring harmful microalgal species and their appearance in Tokyo Bay, Japan, using metabarcoding
expand article infoSirje Sildever§, Noriko Nishi|, Nobuharu Inaba, Taiga Asakura#, Jun Kikuchi#, Yasuhito Asano¤, Takanori Kobayashi, Takashi Gojobori«, Satoshi Nagai
‡ Japan Fisheries Research and Education Agency, Fisheries Resources Institute, Yokohama, Japan
§ Tallinn University of Technology, Tallinn, Estonia
| AXIOHELIX Co. Ltd, Tokyo, Japan
¶ Civil Engineering Research Institute for Cold Region, Public Works Research Institute, Sapporo, Japan
# RIKEN Center for Sustainable Resource Science, Yokohama, Japan
¤ Toyo University, Tokyo, Japan
« King Abdullah University of Science and Technology, Thuwal, Saudi Arabia
Open Access


During the recent decade, high-throughput sequencing (HTS) techniques, in particular, DNA metabarcoding, have facilitated increased detection of biodiversity, including harmful algal bloom (HAB) species. In this study, the presence of HAB species and their appearance patterns were investigated by employing molecular and light microscopy-based monitoring in Tokyo Bay, Japan. The potential co-appearance patterns between the HAB species, as well as with other eukaryotes and prokaryotes were investigated using correlation and association rule-based time-series analysis. In total, 40 unique HAB species were detected, including 12 toxin-producing HAB species previously not reported from the area. More than half of the HAB species were present throughout the sampling season (summer to autumn) and no structuring or succession patterns associated with the environmental conditions could be detected. Statistically significant (p < 0.05, rS ranging from −0.88 to 0.90) associations were found amongst the HAB species and other eukaryotic and prokaryotic species, including genera containing growth-limiting bacteria. However, significant correlations between species differed amongst the years, indicating that variability in environmental conditions between the years may have a stronger influence on the microalgal community structure and interspecies interactions than the variability during the sampling season. The association rule-based time-series analysis allowed the detection of a previously reported negative relationship between Synechococcus sp. and Skeletonema sp. in nature. Overall, the results support the applicability of metabarcoding and HTS-based microalgae monitoring, as it facilitates more precise species identification compared to light microscopy, as well as provides input for investigating potential interactions amongst different species/groups through simultaneous detection of multiple species/genera.

Key Words

18S ribosomal RNA gene, bacteria, bloom-forming species, next-generation sequencing, phytoplankton, toxin-producing species


Microscopic algae are unicellular eukaryotic or prokaryotic organisms present in various aquatic environments. Microalgae may exist as individual cells or form chains or colonies. Some of them are important as primary producers, generating about 48% of the annual net primary production (Field et al. 1998). About 300 microalgal species are associated with harmful algal blooms (HABs; Hallegraeff 2004; Berdalet et al. 2016), a proliferation of microalgal cells, where cell abundances may reach up to millions of cells per litre (Zhang et al. 2012; Burson et al. 2014; Du et al. 2016; Viana et al. 2019). The main harmful effects are related to oxygen depletion, attenuation of light conditions, toxicity and clogging of fish gills (Hallegraeff 2004; Anderson 2009). HABs can also have a notable economic impact due to fish mortality and closure of fisheries, discolouration of seaweed (Pyropia yezoensis) cultures, negative influence on tourism, additional costs on monitoring and public health (Miyahara et al. 1996; Whyte et al. 2001; Imai et al. 2006; Lee et al. 2013; Yamaguchi et al. 2014; Clément et al. 2016; Du et al. 2016; Sanseverino et al. 2016). During recent decades, blooms, areas affected by HABs and distribution ranges of HAB species have expanded (Hallegraeff 2004; Anderson et al. 2012a; Kudela and Gobler 2012; Zhang et al. 2012; Koch et al. 2014; Shimada et al. 2016; Natsuike et al. 2019).

To prevent or minimise the adverse effects associated with HABs, information on HAB species distribution and factors supporting/hindering their growth, for example, environmental parameters and other organisms, is needed. Traditionally, microalgae, including the HAB species, have been identified and enumerated, based on their morphology using light microscopy (Nishikawa et al. 2010; Lefebvre et al. 2011; Lee et al. 2013; Yasakova 2013; Eriksen et al. 2019). However, the exact identification may be hampered by similar or variable morphology, small cell size as well as by the influence of fixatives on cells (Rhodes 1998; John et al. 2005; Zingone et al. 2006; Galluzzi et al. 2008; Gran-Stadniczeñko et al. 2017). Improved species detection can be achieved by employing high-throughput sequencing (HTS) methods, such as DNA metabarcoding. This approach allows amplifying the gene of interest from a mass collection of organisms or environmental DNA (eDNA) and thus target multiple species simultaneously (Moreno-Pino et al. 2018; Nagai 2018; Nagai et al. 2019). Metabarcoding has also been applied for investigating the diversity of microalgae, which has uncovered the presence of previously unknown HAB species/genera for the sampling area (e.g. Dzhembekova et al. 2017; Elferink et al. 2017; Nagai et al. 2017; Moreno-Pino et al. 2018; Gran-Stadniczeñko et al. 2019; Sildever et al. 2019, 2021; Sze et al. 2019; Liu et al. 2020).

Furthermore, eDNA metabarcoding has been successfully applied to reveal appearance patterns for several organism groups ranging from bacteria to mammals (Nagai et al. 2016; Hirai et al. 2017; Sigsgaard et al. 2017; Stoeckle et al. 2017; Berry et al. 2019; Zhang et al. 2020; Sildever et al. 2021; Alter et al. 2022), including the HAB species (Nagai et al. 2017, 2019; Sildever et al. 2019). As several genes or markers can be amplified from the same sample, co-appearance patterns and associations amongst species and groups can also be investigated (Lima-Mendez et al. 2015; Sawaya et al. 2019; Djurhuus et al. 2020). This may be especially useful as changes in species diversity, presence of different bacteria and parasites can be monitored continuously over a longer time scale together with the HAB species to recognise changes in those parameters/organism groups as indicators of the state of the HABs (Yang et al. 2015; Hattenrath-Lehmann and Gobler 2017; Berdjeb et al. 2018; Shin et al 2018; Hattenrath-Lehmann et al. 2019; Nagai et al. 2019; Jankowiak and Gobler 2020; Yarimizu et al. 2020).

During recent decades, several bacteria displaying growth-limiting influence on the HAB causative microalgae have been isolated from the natural environment to find a method to control the HAB occurrence and duration (Imai et al. 1993; Lovejoy et al. 1998; Nagai and Imai 1999; Skerratt et al. 2002; Kim et al. 2008; Zheng et al. 2018; Inaba et al. 2019, 2020). However, currently, no data are available to demonstrate the growth-limiting influence of bacteria on HAB causative species in nature (Mayali and Azam 2004; Meyer et al. 2017), although their influence on modifying the natural phytoplankton communities has been demonstrated in a laboratory experiment (Bigalke et al. 2019). Thus, both eukaryotes and prokaryotes were targeted in this study in Tokyo Bay from May/June to October 2017 and 2018 using metabarcoding.

Tokyo Bay is a semi-enclosed bay located in central Japan, connected with the Pacific Ocean through the Uraga Channel. The area of the Bay is 1500 km2 (Hattori 1982) with an average depth of 45 m (Kokubu et al. 2013). In total, 36 rivers flow to the bay (Nihei et al. 2009) with higher discharges between June to September (Kokubu et al. 2013). Dominant seasonal winds are from the northeast during winter and from the southwest during summer (Oda and Kanda 2009). Water exchange between the ocean and Bay is driven by the estuarine circulation and is influenced by freshwater input and solar radiation (Nakayama 2006). During summer, the water residence time is around 15 to 35 days (Takao et al. 2004). The main currents in Tokyo Bay are tidal and residual currents, whereas the offshore current from the ocean rarely influences the circulation in the inner bay (Suzuki and Matsuyama 2000). There is a notable seasonal variation in residual current due to variations in river outflow, heating and cooling of the sea surface and prevailing winds (Guo and Yanagi 1996). The water column in the inner part of Tokyo Bay is strongly stratified from June to October, whereas from November to March, the water column is mixed through (Nakane et al. 2008; Bouman et al. 2010). In summer, diurnal amplitude in sea surface temperature is > 1 °C, with the maximum measured difference of 5.5 °C (Oda and Kanda 2009).

Microalgal abundance is the highest between March and October and the community is dominated by diatoms (Nakane et al. 2008). The Bay became eutrophicated between the 1950s and 1970s, which was also reflected in the microalgal community composition as a loss of oceanic/neritic species replaced by a few dominant species (Nomura 1998). The annual average occurrence of HABs also increased with eutrophication from an average of two HABs per year until the 1940s to 19 times per year in the 1980s (Nomura 1998). The post-bloom accumulation of decaying algal cells in the seabed contributes to the development of hypoxia, resulting in losses in benthic communities (Toba et al. 2008; Kodama and Horiguchi 2011). The main HAB causative species in Tokyo Bay are a diatom Skeletonema costatum and a raphidophyte Heterosigma akashiwo (Nomura 1998) with 28 toxin-producing species also recorded from the Bay (Nomura and Yoshida 1997; Nomura 1998; Matsuoka et al. 2003; Yap-Dejeto et al. 2010; Yuki and Yoshimatsu 2012; Demura et al. 2014; Nagai et al. 2017).

This study focuses on the usefulness of metabarcoding for microalgal monitoring with emphasis on HAB species by using universal primers targeting the 18S ribosomal RNA gene. Additionally, patterns in the HAB species appearance and associations amongst HAB species, as well as between HAB species and environmental parameters, other eukaryotes and bacteria, were investigated.


Sampling and DNA extraction

Surface seawater was collected weekly between June 2017 and October 2017 (n = 21), as well as between May 2018 and October 2018 (n = 24) from Tokyo Bay, Japan (35°34.60'N, 139°65.72'E, Fig. 1) by using a bucket. Water temperature and salinity were measured during sampling by a YSI Pro30 conductivity and temperature sensor (Xylem Inc., USA). A 50 ml subsample was taken for nutrient and ion chemical analysis (NO2, NO3, PO4, Si, Na, Mg, S, Ca, K, Sr, B and Li). Inductively coupled plasma optical emission spectroscopy (ICP_OES, SPECTROBLUE, AMETEK Inc. USA), equipped with an autosampler ASX-260 (Cetac Technologies) and Smart Analyzer Vision 6.01.0945 software, was employed for ion analysis (Sekiyama et al. 2011; Ito et al. 2014; Ogawa et al. 2014; Oita et al. 2018;). The ICP-OES operating conditions were as follows: power 1.4 kW, plasma gas flow 12 l min−1, auxiliary gas flow 1.0 l min−1, nebuliser gas flow 1.00 l min−1 and peristaltic pump speed 30 rpm. From the obtained data, a matrix was built using the concentration in ppm for each element against the sampling points from the average result of the optimal dilution with the optimal wavelengths. Nutrient analysis was performed by a DR 3900 spectrometer (Hach, United States) using the following reagent kits: Total nitrogen – TNT826, ammonia – TNT830, nitrate – TNT835, nitrite – TNT839 and phosphate – TNT535. The analytical protocols followed the vendor’s manual. For Chl a analysis, 100 ml of water was filtered through GF/F filter (Whatman, GE Healthcare, Tokyo, Japan) and stored in an N,N-Dimethylformamide at −20 °C until analysed following the workflow in Holm-Hansen et al. (1965) and using a Turner Designs fluorometer (10 AU).

Figure 1.

Sampling location in Tokyo Bay, Japan. The red square symbolises Tokyo Bay on the Pacific coast of Japan and the red-filled square indicates the sampling location in Tokyo Bay.

Microalgae were identified and enumerated from 1 ml of raw seawater and also from 1 ml of Lugol (2%) fixed sample concentrated from 200 ml of water to 8–10 ml by using an inverted microscope (Eclipse Ti-U, Nikon, Tokyo, Japan) and a Sedgewick Rafter counting chamber. The cells in the entire chamber were identified and enumerated. Morphologically similar species, for example, belonging to Chaetoceros, Pseudo-nitzschia and Skeletonema, were identified to genus level. Microscopic count data on species/genus level is only available from 2018, whereas total microalgal cell count data is available for 2017. For metabarcoding, 200–1000 ml of seawater was filtered through 8 and 1 µm (in 2017) or 1 µm (in 2018, 8 µm filtration was discontinued to reduce the number of steps in the DNA extraction) and 0.22 µm pore-size filters (Nuclepore membrane, GE Healthcare, Tokyo, Japan). DNA was extracted after filtration or stored at −20 °C until the extraction by using a 5% Chelex buffer (Chelex 100, Molecular Biology Grade Resin, Bio-Rad, Hercules, CA, USA), following the protocol described by Nagai et al. (2012). First, 250 µl of the 5% Chelex buffer was added to the 1.5 ml tubes containing the filters. A pellet pestle motor (Kontes Glass, Vineland, NJ, USA) was used for 30 s to break the cells on the filters. As the final step in the extraction, the samples were heated for 20 min at 97 °C.

Paired-end library preparation, sequencing and bioinformatics

For 2017 samples, DNA extracted from 1 and 8 µm filters, was mixed in equal volumes and used as a template for paired-end library preparation. For 2018, DNA extracted from 1 µm filters, was used as a template. Eukaryotic species were targeted by universal primers for the 18S rRNA gene (V7-V9 region, Dzhembekova et al. 2018). For prokaryotic species, universal 16S rRNA gene primers (V3-V4 region, Herlemann, et al. 2011) were used with DNA extracted from both filter fractions (0.22 µm – 1 µm and ≥ 1 µm) as a template. Two-step PCR for the construction of paired-end libraries and HTS on the Illumina MiSeq platform (MiSeq Reagent Kit v.3, 600-cycle; Illumina, San Diego, CA, USA) followed the protocol in Dzhembekova et al. (2017). Treatment of obtained sequences, selection of operational taxonomic units (OTUs) and taxonomic identification of OTUs were done according to the workflow described by Dzhembekova et al. (2017) with the exception that sequences longer than 300 bp were truncated to 300 bp by trimming the 3′ tails.

The demultiplexing and trimming were performed using Trimmomatic version 0.35 ( Nucleotide sequences were demultiplexed according to the 5´-multiplex identifier (MID) tag and primer sequences were demultiplexed according to the default format in MiSeq. Sections containing: (1) palindrome clips longer than 30 bp and (2) monopolymers longer than 9 bp were trimmed from the sequences at both ends. 3´ tails with an average quality score lower than 30 at the end of the final 25-bp window were also trimmed from each sequence. 5´ and 3´ tails with an average quality score lower than 20 at the end of the final window were also trimmed. The remaining sequences were merged into paired reads using Usearch version 8.0.1517 ( with default settings (≥ 16 bp overlap, ≥ 90% similarity and mismatch ≤ 5 bp; resulting in a maximum sequence length of 584 bp. Singletons were also removed. Sequences were aligned using Clustal Omega v. 1.2.0. ( and only sequences that were contained in more than 75% of the read positions were extracted. Filtering and a part of the multiple alignment process were performed using the screen.seqs and filter. seqs commands in Mothur, as described in the MiSeq SOP (; Schloss et al. 2011). Erroneous and chimeric sequences were detected and removed using the pre.cluster (diffs = 4) and chimera.uchime (minh = 0.1;; Edgar et al. 2011) commands in Mothur, respectively. Using the unique.seqs command of Mothur, the same sequences were collected into OTUs. The contig sequences were counted as OTUs by count.seqs and used for the subsequent taxonomic identification analysis. Demultiplexed, filtered, but untrimmed sequence are available in the DDBJ Sequence Read Archive (Accession number: DRA013265). Sequences were clustered to OTUs at ≥ 99% similarity level. The sequence database, used for assigning taxonomy to OTUs, was downloaded from GenBank on 22.03.2019. OTUs associated with multiple records from the same genus were merged with the OTU associated with a single species from the same genus if the multiple records consisted of a single species and other records identified as sp. of the same genus. For eukaryotes, only OTUs that had ≥ 99% similarity with the best match from the database were used for further analysis. For prokaryotes, only OTUs with > 50 total sequence reads for each fraction and sampling year were included in the further analysis.

HAB species selection

The presence of HAB species amongst the detected eukaryotic OTUs was investigated by comparing the list of OTUs associated with phytoplankton species against the IOC-UNESCO Taxonomic Reference List of Harmful Micro Algae (Moestrup et al. 2022), as well as with the known non-toxic HAB species from the Tokyo Bay (Takano 1956; Nomura 1998). To avoid the inclusion of OTUs with ambiguous identities, representative sequences of all OTUs associated with the HAB species were also manually BLAST-searched from the GenBank online database ( The validity of the taxonomic names was checked against the AlgaeBase ( and the World Register of Marine Species ( Only OTUs that had the ≥ 99% BLAST top hit similarity with the associated single species were included in the further analysis. For OTUs associated with Alexandrium catenella and Alexandrium tamarense, the nomenclature described by Litaker et al. (2018) was followed. To facilitate the comparison with the LM data, OTUs are referred to as species or genera. For example, species belonging to Pseudo-nitzchia were detected to the genus level by LM; however, they are also shown as single species detected by the metabarcoding and HTS to provide more detailed information. Species belonging to genera/higher taxonomic levels reported by Nomura (1998; Chaetoceros, Euglenophyceae, Gymnodinium, Naviculaceae, Prorocentrum, Pyramimonas, Rhizosolenia and Thalassiosira) are shown on a genus/higher taxonomic level. An exception to this is Cryptomonadaceae spp., which is represented by Goniomonas aff. amphinema as it was the only species from this group detected by metabarcoding and the HTS-approach, whereas no species belonging to this group were detected by LM in this study.

Statistical analysis

All data analysis was conducted in R v. 4.0.3 (R Core Team 2021) using the “vegan” (v.2.5-7, Oksanen et al. 2021) or “psych” packages (correlation analysis; Revelle 2020). The number of shared OTUs amongst different sampling years and methods was analysed and visualised as a Venn diagram. As LM did not allow species-level identification for some genera, the number of shared OTUs is compared both at species and genus levels. Species richness was calculated as the number of OTUs. Ordination analyses were applied to reveal the relationship between the presence/absence of OTUs associated with HAB species and several environmental parameters. First, a Detrended Correspondence Analysis was performed to choose between linear and unimodal ordination methods. The length of the first DCA axis was < 3 standard deviations for 2017 and 2018 data (Table 1) and, thus, linear methods were applied (PCA: Principal Component Analysis and RDA: Redundancy Analysis; Ter Braak and Prentice 1988). PCA was used to investigate general trends in the HAB species appearances with environmental parameters added as vectors using the envfit function. RDA was used to analyse the influence of environmental parameters on the HAB species’ appearances. Monte Carlo Permutation tests with 999 unrestricted permutations were applied to test the significance of the variation in community composition explained by the environmental variables under a global model. When the result was statistically significant, further Monte Carlo permutation tests with 999 unrestricted permutations were performed to test the variation explained by the individual axes and environmental variables. Spearman correlations between environmental parameters and presence/absence of OTUs associated with each HAB species, between OTUs associated with HAB species, other eukaryotic species detected and with bacteria were calculated. P-value was adjusted for multiple comparisons by the Benjamini-Hochberg correction (Benjamini and Hochberg 1995). The correlation analysis results were visualised by heatmap.2 (“gplots” v.3.1.1; Warnes et al. 2020) between eukaryotes and by Gephi (v. 0.9.2; Bastian et al. 2009) between pro- and eukaryotes.

Table 1.

Correlations of environmental parameters with PCA axes and Monte Carlo permutation test results for RDA, based on eukaryotic OTUs associated with HAB species/genera and environmental parameters in 2017 and 2018. Sal: salinity, WT: water temperature, Mg: magnesium.

Dataset DCA 1st axis length Sign. correlation of env. param. with the PCA axes 1 and 2 (R2/p-value) Global test (F/p-value) Sign. axis (F/p-value) Sign. env. param. (terms) (F/p-value) Sign. env. param. (margin) (F/p-value)
2017 1.62 WT (0.46/ 0.01) 1.55/0.002 RDA1 (5.70/0.003) Sal (1.85/0.036) WT (1.84/0.046)
WT (4.20/0.001)
Mg (2.80/0.001)
2018 1.8 0.405

Association rule-based time-series analysis

Association rule-based time-series data analysis (Asano et al. 2019) was applied for investigating the changes in phytoplankton abundance following the changes in bacterial abundance. The method allows us to find association rules considering the potential lag effect (Agrawal et al. 1993), for example, changes in cell abundance or number of live cells may display a lagged response to a change in a parameter (Collos 1986; Nagai and Imai 1999; Śliwińska-Wilczewska et al. 2019). The number of sequence reads for Gymnodinium spp., Pseudo-nitzschia spp. and Skeletonema spp. were analysed, together with the bacterial sequence reads from 2017 and 2018 (non-rarefied data). Those groups were chosen due to their presence throughout most of the sampling season in both years. To obtain a minimum of four association rules per group analysed, the characteristic change in each group was set as follows: p = 0.1 for Gymnodinium spp. and Pseudo-nitzschia spp. and p = 0.3 for Skeletonema spp. The analysis was implemented in R v. 4.0.3 (R Core Team 2021).


Environmental parameters

During both years, salinity and water temperature displayed a similar pattern ranging from 24.5 to 32 and 18.2 °C to 27.9 °C in 2017 and from 23.9 to 32.6 and 17.4 °C to 28.9 °C in 2018, respectively (Fig. 2). The chlorophyll a pattern was more variable between the years with the values from 0.7 to 271.35 µg l-1 in 2017 and from 1.6 to 77.59 µg l-1 (Fig. 2). Nitrite concentrations ranged from 0.015 to 0.042 mg l-1 in 2017, but were mainly below the detection limit for 2018. Nitrate (0.23–0.54 mg l-1) and phosphate (0.02–0.58 mg l-1) data were only available for 2018. Silicate concentrations were almost constant throughout the sampling seasons, probably reflecting an analytical error and were thus excluded (data not shown). Ion concentrations fluctuated throughout the sampling seasons (Fig. 2). Na, Mg and Ca concentrations followed the fluctuations in salinity. Li concentrations had two peaks at the beginning of 2017 with maximum values reaching 103 parts per billion.

Figure 2.

Overview of the environmental parameters, nutrient and ion concentrations (ppb = parts per billion) during sampling periods, blue line: 2017, red line 2018.

Overview of detected eukaryotes and prokaryotes

For eukaryotes, an average of 26,781 (SD: 3,801) sequences were detected per sample in 2017, whereas in 2018, it was 14,508 (SD: 3,593) sequences (after bioinformatics treatment, Suppl. material 2: Tables S2, S3). For prokaryotes, the average sequence reads were 30,974 (SD: 5,729) and 6,695 (SD: 2,276) for samples analysed from 1 µm filters in 2017 and 2018 (Suppl. material 2: Tables S4, S6). For samples from the 0.22 µm filters, the average number of sequences per sample was 43,532 (SD: 4,036) and 13,676 (SD: 2,116) in 2017 and 2018 (Suppl. material 2: Tables S5, S7). The OTU accumulation curves (Suppl. material 1: Figs S1–S2) indicate that the number of OTUs is nearly saturated for samples from 2017 and 2018, based on the 18S marker, but in the case of the 16S marker, the low number of OTUs and sequences for 2018 samples are visible.

In 2017, 673 unique eukaryotic OTUs (≥ 99% BLAST top-hit similarity) were detected, belonging to nine different supergroups and 32 phyla (Suppl. material 2: Table S2). In 2018, 616 eukaryote OTUs from nine supergroups and 34 phyla were present (Suppl. material 2: Table S3). In both years, the majority of the OTUs belonged to Metazoa, Dinophyceae or Bacillariophyta (Suppl. material 2: Tables S2, S3). The minimum number of unique OTUs were 120 and 105 and the maximum number of OTUs were 265 and 234 in 2017 and 2018, respectively (Fig. 2). The number of cells detected by microscopy ranged from 3 to 10,954 cells ml-1 in 2017 and 2018 (Fig. 2). Unfortunately, species-level information, based on light microscopy, is not available for 2017, but in 2018, the number of species detected per sample ranged from 17 to 50 with an average of 34 species (Fig. 2).

In prokaryotes, 715 unique OTUs (≥ 50 sequences) were present in the 0.22–1 µm bacteria fraction () in 2017, whereas 1092 OTUs were detected from the ≥ 1 µm bacteria fraction (Suppl. material 2: Tables S4, S5). In 2018, the number of OTUs detected was lower for both fractions, 424 and 413, respectively (Suppl. material 2: Tables S6, S7). In 2017, the prokaryote OTUs belonged to 19 and 33 phyla, respectively (the 0.22–1 µm / ≥ 1 µm fractions; Suppl. material 2: Tables S4, S5), whereas in 2018, the respective values were 17 and 22 (Suppl. material 2: Tables S6, S7). In 2017, the ≥ 1 µm bacteria were mainly dominated by Bacteroides/Chlorobi group and Gammaproteobacteria, whereas in the 0.22–1 µm group, bacteria associated with Alphaproteobacteria and Actinobacteria were most abundant. In 2018, the main groups represented by the bacterial OTUs were the same for both fractions: Bacteroides/Chlorobi group, followed by Gamma- and Alphaproteobacteria.

Number of HAB species/genera detected and appearance patterns

In total, 29 unique genera associated with HAB species and 40 unique HAB species were detected from Tokyo Bay during two years of monitoring, based on metabarcoding and the HTS approach and morphology-based identification under LM (Suppl. material 2: Table S1). Of the 40 species detected, 30 are known as toxin-producers (Moestrup et al. 2022) and 12 of these have not been previously reported from Tokyo Bay (Suppl. material 2: Table S1). Metabarcoding and the HTS-approach allowed the detection of 26 genera, whereas 17 genera were detected, based on morphology using LM. However, three of those (Dactylisolen, Eucampia and Leptocylindrus) were not detected by metabarcoding and HTS (Fig. 3, Suppl. material 2: Table S1). The genera detected in 2017 were also present in 2018, whereas four genera were only registered in 2018 (Goniomonas, Phaeocystis, Pseudochattonella and Rhizosolenia; Suppl. material 2: Table S1). Based on metabarcoding and HTS data from both years, 36 HAB species were detected compared to 14 HAB species identified under LM in 2018 (Fig. 3, Suppl. material 2: Table S1). In 2018, five additional species (Alexandrium catenella, A. pacificum, Goniomonas aff. amphinema, Phaeocystis globosa and Pseudochattonella verruculosa) were detected by HTS compared to 2017. The total number of species/genera detected by metabarcoding and the HTS approach was 43, whereas LM allowed the detection of 20 species/genera.

Figure 3.

A. Number of unique genera; B. Number of unique species detected, based on light microscopy (LM) and by metabarcoding and the HTS approach in 2017 and 2018. As LM did not allow species-level identification for some genera, the number of shared OTUs is compared both at species and genus levels.

Based on metabarcoding and the HTS approach, 65% of the detected HAB species/genera (n = 37) were present in more than half of the sampling occasions in 2017 (June – October, n = 21) with Karlodinium veneficum, Chaetoceros spp., Prorocentrum spp., Skeletonema spp. and Thalassiosira spp. detected from all the samples (Fig. 4). At the same time, A. minutum and Azadinium poporum were detected only once or twice, respectively (Fig. 4). In 2018, 42% of the detected species/genera (n = 43) were present in more than half of the samples (May – October, n = 24). Similar to 2017, K. veneficum, Chaetoceros spp. Skeletonema spp. and Thalassiosira spp. were detected from all the samples, whereas A. catenella, Goniomonas spp., Mesodinium rubrum, P. verruculosa and Rhizosolenia spp. were present in one sample and A. pacificum in two samples (Fig. 4). In morphology-based identification by LM, 50% of the detected species/genera were present in more than half of the sampling occasions in 2018 (Suppl. material 1: Fig. S3). Chaetoceros species were detected from all the samples and species belonging to Pseudo-nitzschia were present in more than 90% of the samples (Suppl. material 1: Fig. S3). Karenia mikimotoi, K. papilionaceae, Polykrikos hartmanii, Phaeocystis globosa and Prorocentrum triestinum were detected once, whereas, based on metabarcoding and the HTS approach, the species were present in 19, 4, 12, 3 and 21 samples, respectively (Fig. 4, Suppl. material 1: Fig. S3). On the contrary, Mesodinium rubrum was identified from 15 samples using LM, whereas, based on metabarcoding and HTS, the species was detected from only one sample in 2018.

Figure 4.

Presence-absence of HAB species/genera detected by metabarcoding and the HTS approach from 2017 and 2018. The black line indicates the division between the two years, the red colour indicates presence and the white absence. The left colour panel indicates months. Toxin-producing species (Moestrup et al. 2022) are indicated in bold. C: Cryptophyceae; CI: Ciliophora; RA: Raphidophyceae.

Statistical analysis


In the PCA, the first two axes explained 35.23% (2017) and 30.62% (2018) of the variability in the OTU relative abundances (Suppl. material 1: Fig. S4). In 2017, the water temperature was significantly (p < 0.05) associated with the variability represented by the first two axes, whereas no significant associations were detected, based on the 2018 dataset (Table 1). In the RDA, the environmental parameters explained 69.98% of the variability in the 2017 dataset and 55.79% in the 2018 dataset. The Monte Carlo permutation test indicated that, in 2017, significant amount of variation (p = 0.002) in the species composition can be explained by the environmental parameters; however, for the 2018 dataset, no significant influence was detected (Table 1). The permutation test for individual axes in the 2017 dataset indicated the importance of the first axis (p = 0.003). Salinity, water temperature and magnesium (Mg) concentration were identified as parameters significantly influencing the community composition (p < 0.05; Table 1). In both years, the majority of OTUs clustered to the centre of the plot and some of those were placed together with a change in the environmental parameter in the RDA. A high number of OTUs clustered together with an increasing water temperature in both years, whereas for other parameters, the pattern was not as clear (Fig. 5). In the RDA from 2017, Prorocentrum cordatum was placed close to increasing Mg concentrations and Noctilluca scintillans with increasing water temperature (Fig. 5A). In 2018, Phaeocystis globosa appeared together with increasing phosphate concentrations, whereas Phalacroma rotundatum and Pseudo-nitzschia pungens were placed together with increasing salinity and nitrate concentrations (Fig. 5B).

Figure 5.

Redundancy analysis, based on HAB-associated OTUs detected by metabarcoding and HTS approach, A. 2017; B. 2018.


Environmental parameters

The majority of the OTUs associated with HAB species/genera had no statistically significant (p > 0.05) correlations with the environmental parameters during both years (data not shown). An exception was Cerataulina pelagica displaying a significant positive correlation with the water temperature in 2017 (p < 0.05, rS = 0.71). Additionally, Pseudo-nitzschia delicatissima had significant negative correlations with silicate (p = 0.02, rS = −0.75) and lithium (p = 0.02, rS = −0.74). In 2018, a significant positive correlation (p = 0.004, rS = 0.81) was detected between Phaeocystis globosa and silicate concentrations.

Other HAB species

In 2017, 14 OTUs associated with the HAB species/genera out of 37 were significantly (p < 0.05) correlated with other HAB-associated OTUs and in 2018, 19 OTUs from 43 were also significantly correlated with other HAB-associated OTUs (Suppl. material 2: Tables S8–S9, Suppl. material 1: Fig. S5). The highest number of significant correlations per OTU was similar amongst the years (4 and 3 per OTU, respectively); however, OTUs showing the highest number of correlations differed amongst the years (Suppl. material 2: Tables S8–S9). In both 2017 and 2018, significant positive correlations were detected between P. australis and P. multiseries (2017: p < 0.0001, rS = 0.84, 2018: p = 0.024, rS = 0.63) as well as between P. cuspidata and P. turgidula (2017: p = 0.02 rS = 0.67, 2018: p = 0.001, rS = 0.74).


In both years, all of the OTUs associated with the HAB species/genera had significant (p < 0.05) correlations with other OTUs representing various eukaryotes (Suppl. material 2: Tables S8–S9, Suppl. material 1: Fig. S5). From eukaryotic OTUs, around 17% and 21% had statistically significant correlations with HAB-associated OTUs in 2017 and 2018, respectively (Suppl. material 2: Tables S8–S9). A maximum of 30 statistically significant correlations were detected per HAB-associated OTU (Goniomonas sp.; Suppl. material 2: Tables S8–S9). In both 2017 and 2018, P. australis was significantly correlated with Eucampia spp. (2017: p = 0.009, rS = 0.77, 2018: p = 0.01, rS = 0.72) and also with Thalassiosira rotula (2017: p = 0.009, rS = 0.77, 2018: p = 0.01, rS = 0.72).


In 2017, more than half of the OTUs associated with HAB species had significant (p < 0.05) correlations with bacteria belonging to 0.22–1 µm and ≥ 1 µm fractions (Suppl. material 2: Tables S10–S11, Fig. 6). In 2018, more than half of the HAB-associated OTUs were significantly correlated with bacteria from the ≥ 1 µm fraction, whereas 43% displayed a significant correlation with the 0.22–1 µm bacteria fraction. The correlations between OTUs associated with HAB and bacteria were between different OTUs in different years. Based on both size fractions and years, 35 OTUs associated with genera containing growth-limiting bacteria were detected that were significantly (p < 0.05) correlated with HAB-associated OTUs (Suppl. material 2: Table S12). In the case of the 0.22–1 µm bacteria fraction, most correlations with HAB species were positive in both years, whereas, in the case of the ≥ 1 µm bacteria fraction, the majority of the significant correlations were negative (Fig. 6). For most of the bacterial OTUs, the previously reported target HAB species did not match with the HAB-associated OTUs that had a significant correlation.

Figure 6.

Statistically significant (p < 0.05) correlations between prokaryotes and eukaryotes. Red colour lines indicate positive and blue colour negative correlation. GI: growth-inhibiting bacterial genera. A,B. Data from 2017 and 2018 with bacterial data from the 0.22–1 µm fraction; C, D. data from 2017 and 2018 with bacterial data from the ≥ 1 µm fraction.

Association rule-based analysis of time-series data

Based on the time-series analysis of selected phyla, several associations between bacteria and phytoplankton were detected (Suppl. material 2: Table S13). From the investigated genera, the highest number of associations was found between bacteria and Skeletonema spp.: 11 associations with two of those negative. For Gymnodinium spp. and Pseudo-nitzschia spp., nine and four associations were detected respectively and for both genera, one of the associations was negative.


Two years of microalgae monitoring by metabarcoding and HTS in combination with light-microscopy allowed the detection of 29 unique genera with 40 unique species associated with harmful algal blooms (Nomura 1998; Moestrup et al. 2022) with 12 toxin-producing species recorded for the first time from Tokyo Bay. The majority of the harmful algal species/genera detected did not display significant correlations with the environmental parameters investigated, whereas significant correlations were detected with other HAB species/genera, other eukaryotes and prokaryotes, including genera containing growth-limiting bacteria.

HAB species/genera detected by metabarcoding and HTS and by light microscopy

Overall, the metabarcoding and HTS-based approach registered 28 HAB-associated species/genera not identified by LM, whereas LM allowed the detection of four species and three genera not found by metabarcoding and HTS. Differences in species/genera detected by the metabarcoding and HTS-approach vs. LM have been previously explained by a poor match between the sequences of targeted taxa and the primers used, a low number of rRNA gene copies per cell in some genera, pre-filtration of samples and availability of sequences in the reference databases (Majaneva et al 2012; Eiler et al. 2013; Xiao et al. 2014; Abad et al. 2016; Banerji et al. 2018; Gran-Stadniczeñko et al. 2019). In the present study, the majority of species not detected by LM may have been unnoticed due to their small size (≤ 20 µm, for example, Azadinium spp., K. veneficum, Pyramimonas spp., Cryptomonadaceae represented by Goniomonas aff. amphinema; Martin-Cereceda et al. 2010; Moro et al. 2011; Wang et al. 2011; Luo et al. 2017) or lack of morphological features identifiable under LM (e.g. Alexandrium spp., Pseudo-nitzschia spp., T. mala; Rhodes 1998; Anderson et al. 2012b; Lim et al. 2013; John et al. 2014; Razali et al. 2016; Bates et al. 2018; Prasad et al. 2018). Furthermore, species present in low abundance may have been missed by LM (Rodriguez-Ramos et al. 2014; Engesmo et al. 2018), although both raw and concentrated phytoplankton samples were analysed. As the concentrated samples were fixed by Lugol’s iodine solution, identification may have also been hampered due to the changes in cell morphology (e.g. Fibrocapsa japonica, Vrieling et al. 1995) or due to discolouration of cells, which could mask the morphological features (Steidinger and Tangen 1996).

From the four species identified only by LM, Dinophysis caudata and Eucampia zodiacus were detected by metabarcoding and HTS, but the sequences also matched with other species from the database and were, thus, not included in further analysis. At the same time, Dactyliosolen fragilissimus and Leptocylindrus minimus were only detected by LM. As both species have been previously reported from Tokyo Bay as HAB causative species (Nomura 1998), the lack of detection by metabarcoding may be explained by the lack of sequences present in the public databases as no 18S rRNA gene sequences were available for D. fragilissimus and only one sequence covering the target region was available for L. minimus. The lack of sequences covering the target region in the database could also explain why Pseudo-nitzschia galaxiae and P. subfraudulenta that have been previously reported from the Bay (Nagai et al. 2017) were not detected. In a recent study from the Gulf of Naples, the Mediterranean Sea, P. galaxiae was successfully detected from eDNA samples by targeting the 18S rRNA gene V4 region (Ruggiero et al. 2022). Thus, until the database is improved, the V4 region of the 18S rRNA gene and other genes (e.g. 28S rRNA) could be used for improved detection of species belonging to Pseudo-nitzschia.

Two toxin-producing dinoflagellate species, Margalefidinium polykrikoides and Lingulodinium polyedra were detected by the region targeted, but their similarity scores with the best database match were lower than the 99% criteria used in this study. Thus, the availability of more sequences covering the target region might also be beneficial for improved identification of those species. Interestingly, another toxin-producing dinoflagellate species, Protoceratium reticulatum, was not detected in the samples, although the cysts of this species have been reported from the sediments throughout Tokyo Bay (Matsuoka et al. 2003). In a study from the Baltic Sea, the species could be successfully identified from eDNA samples using the same marker (Sildever et al. 2021). Thus, the lack of detection in this study could result from collecting data from only a single station, which does not capture the entire diversity of the Bay. In future studies, more stations should be included to detect as much of the diversity present as possible.

Some of the species identified by both metabarcoding and HTS, as well as by LM, were detected in notably different numbers of sampling occasions: Karenia mikimotoi, K. papilionaceae, Noctiluca scintillans, Prorocentrum cordatum and P. triestinum were registered in more samples by metabarcoding and HTS, whereas Mesodinium rubrum and Rhizosolenia spp. were recorded in a higher number of samples by LM. In the case of Karenia and Prorocentrum species, improved detection by metabarcoding and HTS may result from more precise taxonomic identification compared to LM, whereas for N. scintillans, the difference may be due to its low abundance in some sampling occasions as the species is large and easily recognisable. Further, the discrepancy in species identification between the two methods has also been explained by the differences in the sample volume analysed (Rodriguez-Ramos et al. 2014; Gran-Stadniczeñko et al. 2019): 200 to 1,000 ml by metabarcoding and HTS compared to 1 ml of raw and 1 ml of concentrated sample (total volume of 8–10 m from 200 ml) analysed by LM in this study.

In the case of M. rubrum and Rhizosolenia spp., the potential influence of the gene and the region targeted on the detection efficiency by metabarcoding and HTS is not supported by the continuous detection of M. rubrum by the metabarcoding and HTS in 2017 from Tokyo Bay, as well as by the several Rhizosolenia species recorded by the same primers from another locality in Japan (S. Nagai, unpubl. data). Interestingly, in 2017, Rhizosolenia spp. was identified from several samples by metabarcoding and HTS, but had a lower similarity (94%) with the best taxonomic match from the database. Therefore, targeting several genes and usage of genus-specific primers may improve species detection (Gran-Stadniczeñko et al. 2017; Nagai et al. 2017; Smith et al. 2017; Sildever et al. 2019, 2021). The low identification of those taxa could also be related to the overall lower number of sequences obtained per sample in 2018 (average 14,508 ± 3,592 sequences after bioinformatics treatment) compared to 2017 (26,781 ± 3,709 sequences after bioinformatics treatment). This is further supported by the OTU accumulation curve displaying lower accumulation levels for the 2018 data. Thus, re-sequencing the samples from 2018 might improve the detection of those species.

HAB species/genera appearance patterns and correlation with environmental parameters

The majority of detected HAB species were present throughout the sampling season, indicating a wide tolerance range for salinity and temperature. This is further exemplified by Pseudo-nitzschia species appearing from March to October in variable temperature and salinity conditions (this study; Yap-Dejeto et al. 2010; Nagai et al. 2017). An exception to this was P. pseudodelicatissima, present mainly in July-October (this study; Yap-Dejeto et al. 2010), potentially due to the preference for higher surface water temperatures as its growth rate increases with rising temperature (Lundholm et al. 1997). In Tokyo Bay, an increase in water temperature and input of nutrients are mainly responsible for short-term fluctuations in the microalgal community composition and abundance as phosphate is limiting growth during the stratification period from summer to autumn (Nakane et al. 2008; Kubo et al. 2019). This was also visible, based on the ordination analyses from both years, where around half of the HAB-associated OTUs grouped with the increasing water temperature. In addition, based on the permutation test, changes in water temperature could explain a significant amount of variability in the HAB community composition in 2017. However, as the statistically significant (p < 0.05) correlations with the water temperature and chemical elements were not continuous between the years, a general pattern in the HAB species appearances, based on environmental parameters, could not be detected.

OTUs associated with A. catenella, Goniomonas aff. amphinema and P. verruculosa were detected only once during the study period. For example, A. catenella was registered in the second week of July 2018, when the water temperature was 26.6 °C and, two weeks later, the species was also detected from another sampling location in Tokyo Bay (data not shown), where it also remained the single detection in both years. In different coastal areas of Japan, the vegetative cells of A. catenella have been reported from a maximum of 20 °C (Sekiguchi et al. 1986; Itakura et al. 2002; Yamamoto et al. 2017). In the coastal waters of north-eastern Japan, the cells are present from January to June (Kaga et al. 2006) with the highest cell densities at 8–9 °C (Ichimi et al. 2001). In Tokyo Bay, the species formed a bloom in June of 1984 (Han and Terazaki 1993); however, the presence of vegetative cells has not been registered in later studies (Nomura and Yoshida 1997). Furthermore, the resting cysts detected from the Bay and belonging to either A. catenella or A. pacificum, lacked cell content, potentially indicating the absence of a local population (Matsuoka et al. 2003; Kotani et al. 2006). Thus, the detection of A. catenella only once throughout the sampling seasons and at high surface water temperature is potentially explained by the transport by currents from the neighbouring coastal areas as the species has caused paralytic shellfish poisoning events on the Pacific coast close to the entrance to Tokyo Bay (Kotani et al. 2006).

From other species registered only once during the study period, P. verruculosa was detected in the second week of May in 2018, when the water temperature was 18.7 °C and salinity 31.9. Interestingly, the temperature and salinity were in the previously-reported ranges suitable for the growth of Japanese P. verruculosa throughout the majority of sampling occasions in both years (Yamaguchi et al. 1997; Skjelbred et al. 2013). In southern Japan and the Seto Inland Sea, the species has been reported from January and February (Yamaguchi et al. 1997; Orita 2016), whereas in the Sea of Okhotsk, the species appears throughout the year with the highest number of occurrences in February-March and in August (Sildever et al. 2019). As there are no previous published results on the presence of P. verruculosa in Tokyo Bay, further study is needed to investigate if the species could have been transported to the area as suggested for A. catenella or if the species occurs earlier in spring than the start of sampling periods in this study. In the case of Cryptomonadaceae, presence from June to August with the highest abundances in August in Tokyo Bay has been reported (Nomura 1998). In this study, metabarcoding and HTS could detect only one representative of this family, Goniomonas aff. amphinema, in October 2018, whereas other species from the same class (Cryptophyta) were registered continuously throughout the sampling season. Since it is not certain which species were considered as Cryptomonadaceae by Nomura (1998), for example, if other species from Cryptophyta were also included, it is not possible to compare the detected occurrence of Cryptomonadaceae only based on G. aff. amphinema with the previous information on the appearances of Cryptomonadaceae in Tokyo Bay.

HAB species/genera appearance patterns and correlation with other species/genera

Another important aspect of structuring a species’ appearance is interaction with other organisms (Lima-Mendez et al. 2015; Sourisseau et al. 2017; Zhou et al. 2018; Sassenhagen et al. 2020). In both sampling seasons, HAB species displayed significant positive and negative associations with other HAB, eukaryotic and prokaryotic species; however, the majority of the detected correlations differed amongst the years, potentially explained by the variability in environmental conditions influencing the interspecific interactions (Mikhailov et al. 2019a). Interestingly, significant positive correlations between the same species in both years were detected for P. australis and P. multiseries, P. cuspidata and P. turgidula, as well as between P. australis and Eucampia spp. and T. rotula. For all of those, presence in a wide range of temperature and/or salinity conditions has been reported (Krawiec 1982; Yap-Dejeto et al. 2010; Nishikawa et al. 2011; Nagai et al. 2017; IOC-UNESCO 2021). In this study, the species appeared throughout the sampling seasons and no significant correlations with the environmental parameters were detected. An exception was P. cuspidata, showing a significant positive correlation with the water temperature in 2017, whereas in 2018, no significant correlation was detected and the species appeared already at lower temperatures (21.2 °C) than in 2017 (26.9 °C). Thus, currently, no clear explanation for the co-appearances in both years can be provided.

Detection of 35 OTUs linked with bacteria belonging to genera associated with growth-limiting effects on HAB species and detection of significant correlations between bacteria and HAB species further confirms the usefulness of metabarcoding and the HTS approach for investigating interactions amongst bacteria, phytoplankton and environmental parameters (Yang et al. 2015; Bunse et al. 2016, 2019; Wurzbacher et al. 2017; Zhou et al. 2018; Cui et al. 2019). Interestingly, the majority of significant correlations found between HAB and bacteria from the genera containing growth-limiting bacteria were positive. However, as the detected correlations are not between HAB species and the bacterial strains with known algicidal activity (Suppl. material 2: Table S12), but with other bacterial strains/species from the same genus, their growth-limiting activity is not certain as it can vary between the strains, species within the same genus and also in different environmental conditions (Skerratt et al. 2002; Meyer and Pohnert 2019). Thus, the significant positive correlations detected in this study may instead indicate mutualism/commensalism (Paver et al. 2013; Mikhailov et al. 2019b) or co-appearance due to indirect reasons, for example, environmental preferences (Weiss et al. 2016).

Bacterial communities are also influenced by the phytoplankton blooms and species composition (Riemann et al. 2000; Teeling et al. 2012; Klindworth et al. 2014; Sison-Mangus et al. 2016; Camarena-Gómez et al. 2018). Bacteria belonging to Alpha- and Gammaproteobacteria, as well as to Flavobacteria are often dominant during the blooms (Buchan et al. 2014). Interestingly, in this study, an OTU associated with Marinobacter spp. was detected only once during the entire sampling season in 2017 (n = 21). In a laboratory study, a strain belonging to this genus impeded the growth of a toxin-producing dinoflagellate species Karenia mikimotoi (Zheng et al. 2018). At the same time, K. mikimotoi was not detected from the sample in which Marinobacter spp. was present although it was present in previous samples. In addition, by the same sampling date, the number of sequences associated with Gammaproteobacteria (class including Marinobacter spp.) had increased ≥ 10x compared to the previous week (435 to 7356 sequences/1% to 16% of the total sequence reads). Gammaproteobacteria have been reported to react later to the algal decay, with some species achieving their highest relative abundances two to four weeks after the peak of algal bloom (Teeling et al. 2012). Thus, the appearance of Marinobacter spp. might be related to the decay of an algal bloom possibly containing a high number of K. mikimotoi cells (hypothesised, based on sequence abundances, no cell count data available). This is further supported by a decrease in Chl a values (12.03 to 1.96 µg/l-1) coinciding with the decrease in K. mikimotoi sequence numbers. To confirm the growth-limiting activity, metabarcoding and HTS data should be combined with laboratory experiments, for example, most-probable number technique, microcosm studies, isolation and identification of the bacteria displaying growth-limiting activities towards the HAB species (Imai et al. 1993; Nagai and Imai 1999; Kim et al. 2008; Inaba et al. 2017; Bigalke et al. 2019).

Another useful method to detect potential growth-limiting bacteria is association rule-based time-series data analysis that also considers the time lag between the change in species abundances (Asano et al. 2019). In the present study, the method was tested with three phytoplankton genera and negative associations with bacteria were detected for all genera. Although further investigation and experiments are needed to confirm the negative influence of those bacteria are needed (e.g. Imai et al. 1995; Kim et al. 2008; Shi et al. 2018; Inaba et al. 2019), it provides a good indication of which bacteria to target. This is further exemplified by the decline in Skeletonema sp. one week after an increase in Synechococcus sp., based on sequence abundances. This coincides with the previously reported negative influence of Synechoccus sp. on Skeletonema sp. (S. marinoi) resulting in reduced growth rates, physical cell damage and a decline in photosynthetic capacity (Śliwińska-Wilczewska and Latała 2018; Śliwińska-Wilczewska et al. 2019). Thus, monitoring the increase of Synechococcus sp. by molecular tools could be used for forecasting the decline of Skeletonema sp. abundances. Identification of similar associations between bacteria and other HAB species may enable a near-real-time bloom forecasting system by combining molecular detection with hydrodynamical and ecological models (McGillicuddy 2010; Brown et al. 2013; Tian and Huang 2019).


A combination of light-microscopy-based monitoring with metabarcoding and the HTS approach allowed the detection of 40 HAB species from Tokyo Bay, including 12 toxin-producing species previously not reported from the area. Metabarcoding and the HTS approach allowed detection of twice as many HAB-associated species than light-microscopy; however, four species were detected only based on morphology under LM. This indicates the importance of using several markers to account for the differences in their identification power. Numerous statistically significant positive and negative associations could be detected amongst the HAB species, as well as amongst HAB, eukaryotic and prokaryotic species, including genera containing growth-limiting bacteria. As the majority of significant correlations were between different OTUs in different years, the interactions between different species were probably more influenced by the variability in environmental conditions between the years than within the sampling season. To understand the interactions between growth-limiting genera and microalgae, including HAB species, a further study combining time-series environmental monitoring by metabarcoding and the HTS approach and laboratory experiments are needed. The results of this study further confirm the applicability of metabarcoding and HTS-based microalgal monitoring, exemplified by the detection of several morphologically similar, small or fragile species previously not reported from Tokyo Bay.


A. Kondo, R. Kubota, K. Kamiya and S. Ori are thanked for their help with the molecular work. This study was supported by grants from the “Technological developments for characterization of harmful plankton in the seawater”, Ministry of Agriculture, Forestry and Fisheries, Japan (16808839) [SN]; Japan Society for the Promotion of Science Short-term Postdoctoral Fellowship (PE18028) [SN, SS]. The writing of this manuscript was supported by the grant from the European Regional Development Fund and the programme Mobilitas Pluss (MOBTP160) and by the Estonian Research Council grant (PSG735) [SS].


  • Abad D, Albaina A, Aguirre M, Laza-Martínez A, Uriarte I, Iriarte A, Villate F, Estonba A (2016) Is metabarcoding suitable for estuarine plankton monitoring? A comparative study with microscopy. Marine Biology 163(7): 1–13.
  • Alter SE, King CD, Chou E, Chin SC, Rekdahl M, Rosenbaum HC (2022) Using Environmental DNA to Detect Whales and Dolphins in the New York Bight. Frontiers in Conservation Science 3: 820377.
  • Anderson DM, Cembella AD, Hallegraeff GM (2012a) Progress in understanding harmful algal blooms (HABs): Paradigm shifts and new technologies for research, monitoring and management. Annual Review of Marine Science 4(1): 143–176.
  • Anderson DM, Alpermann TJ, Cembella AD, Collos Y, Masseret E, Montresor M (2012b) The globally distributed genus Alexandrium: Multifaceted roles in marine ecosystems and impacts on human health. Harmful Algae 14: 10–35.
  • Asano Y, Oikawa H, Yasuike M, Nakamura Y, Fujiwara A, Yamamoto K, Nagai S, Kobayashi T, Gojobori T (2019) Mining of knowledge relating to factors involved in the aberrant growth of plankton. In: Gojobori T, Wada T, Kobayashi T, Mineta K (Eds) Marine Metagenomics. Technological Aspects and Applications. Springer, Singapore, 249–271.
  • Banerji A, Bagley M, Elk M, Pilgrim E, Martinson J, Santo Domingo J (2018) Spatial and temporal dynamics of a freshwater eukaryotic plankton community revealed via 18S rRNA gene metabarcoding. Hydrobiologia 818(1): 71–86.
  • Benjamini Y, Hochberg Y (1995) Controlling the False Discovery Rate - A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B. Methodological 57(1): 289–300.
  • Berdalet E, Fleming LE, Gowen R, Davidson K, Hess P, Backer LC, Moore SK, Hoagland P, Enevoldsen H (2016) Marine harmful algal blooms, human health and wellbeing: Challenges and opportunities in the 21st century. Journal of the Marine Biological Association of the United Kingdom 96(1): 61–91.
  • Berdjeb L, Parada A, Needham DM, Fuhrman JA (2018) Short-term dynamics and interactions of marine protist communities during the spring–summer transition. The ISME Journal 12: 1907–1917.
  • Berry TE, Saunders BJ, Coghlan ML, Stat M, Jarman S, Richardson AJ, Davies CH, Berry O, Harvey ES, Bunce M (2019) Marine environmental DNA biomonitoring reveals seasonal patterns in biodiversity and identifies ecosystem responses to anomalous climatic events. PLoS Genet 15(2): e1007943.
  • Bigalke A, Meyer N, Papanikolopoulou LA, Wiltshire KH, Pohnert G (2019) The algicidal bacterium Kordia algicida shapes a natural plankton community. Applied and Environmental Microbiology 85(7): 1–12.
  • Bouman HA, Nakane T, Oka K, Nakata K, Kurita K, Sathyendranath S, Platt T (2010) Environmental controls on phytoplankton production in coastal ecosystems: A case study from Tokyo Bay. Estuarine, Coastal and Shelf Science 87(1): 63–72.
  • Brown CW, Hood RR, Long W, Jacobs J, Ramers DL, Wazniak C, Wiggert JD, Wood R, Xu J (2013) Ecological forecasting in Chesapeake Bay: Using a mechanistic-empirical modeling approach. Journal of Marine Systems 125: 113–125.
  • Buchan A, LeCleir GR, Gulvik CA, Gonzalez JM (2014) Master recyclers: Features and functions of bacteria associated withphytoplankton blooms. Nature Reviews. Microbiology 12(10): 686–698.
  • Bunse C, Bertos-Fortis M, Sassenhagen I, Sildever S, Sjöqvist C, Godhe A, Gross S, Kremp A, Lips I, Lundholm N, Rengefors K, Sefbom J, Pinhassi J, Legrand C (2016) Spatio-temporal interdependence of bacteria and phytoplankton during a Baltic Sea spring bloom. Frontiers in Microbiology 7: 517.
  • Bunse C, Israelsson S, Baltar F, Bertos-Fortis M, Fridolfsson E, Legrand C, Lindehoff E, Lindh MV, Martínez-García S, Pinhassi J (2019) High frequency multi-year variability in Baltic Sea microbial plankton stocks and activities. Frontiers in Microbiology 10: 3296.
  • Burson A, Matthijs HCP, de Bruijne W, Talens R, Hoogenboom R, Gerssen A, Visser PM, Stomp M, Steur K, van Scheppingen Y, Huisman J (2014) Termination of a toxic Alexandrium bloom with hydrogen peroxide. Harmful Algae 31: 125–135.
  • Camarena-Gómez MT, Lipsewers T, Piiparinen J, Eronen-Rasimus E, Perez-Quemaliños D, Hoikkala L, Sobrino C, Spilling K (2018) Shifts in phytoplankton community structure modify bacterial production, abundance and community composition. Aquatic Microbial Ecology 81(2): 149–170.
  • Clément A, Lincoquero L, Saldivia M, Brito CG, Muñoz F, Fernández C, Perez F, Maluje CP, Correa N, Moncada V, Contreras G (2016) Expectional Summer Conditions and HABs of Pseudochattonella in Southern Chile Create Record Impacts on Salmon Farms. Harnful Algae News 53: 1–3. [accessed: 21.11.2021]
  • Collos Y (1986) Time-lag algal growth dynamics: Biological constraints on primary production in aquatic environments. Marine Ecology Progress Series 33: 193–206.
  • Cui Y, Chun SJ, Baek SH, Lee M, Kim Y, Lee HG, Ko SR, Hwang S, Ahn CY, Oh HM (2019) The water depth-dependent co-occurrence patterns of marine bacteria in shallow and dynamic Southern Coast, Korea. Scientific Reports 9(1): 9176.
  • Demura M, Nakayama T, Kasai F, Kawachi M (2014) Genetic structure of Japanese Chattonella marina (Raphidophyceae) populations revealed using microsatellite markers. Phycological Research 62(2): 102–108.
  • Djurhuus A, Closek CJ, Kelly RP, Pitz KJ, Michisaki RP, Starks HA, Walz KR, Andruszkiewicz EA, Olesin E, Hubbard K, Montes E, Otis D, Muller-Karger FE, Chavez FP, Boehm AB, Breitbart M (2020) Environmental DNA reveals seasonal shifts and potential interactions in a marine community. Nature Communications 11(1): 1–9.
  • Du X, Peterson W, Fisher J, Hunter M, Peterson J (2016) Initiation and development of a toxic and persistent Pseudo-nitzschia bloom off the Oregon coast in spring/summer 2015. PLoS ONE 11(10): e0163977.
  • Dzhembekova N, Urusizaki S, Moncheva S, Ivanova P, Nagai S (2017) Applicability of massively parallel sequencing on monitoring harmful algae at Varna Bay in the Black Sea. Harmful Algae 68: 40–51.
  • Dzhembekova N, Moncheva S, Ivanova P, Slabakova N, Nagai S, Dzhembekova N (2018) Biodiversity of phytoplankton cyst assemblages in surface sediments of the Black Sea based on metabarcoding. Biotechnology, Biotechnological Equipment 32(6): 1–7.
  • Eiler A, Drakare S, Bertilsson S, Pernthaler J, Peura S, Rofner C, Simek K, Yang Y, Znachor P, Lindström ES (2013) Unveiling Distribution Patterns of Freshwater Phytoplankton by a Next Generation Sequencing Based Approach. PLoS ONE 8(1): e53516.
  • Elferink S, Neuhaus S, Wohlrab S, Toebe K, Voß D, Gottschling M, Lundholm N, Krock B, Koch BP, Zielinski O, Cembella A, John U (2017) Molecular diversity patterns among various phytoplankton size-fractions in West Greenland in late summer. Deep-Sea Research Part I 121: 54–69.
  • Engesmo A, Strand D, Gran-Stadniczeñko S, Edvardsen B, Medlin LK, Eikrem W (2018) Development of a qPCR assay to detect and quantify ichthyotoxic flagellates along the Norwegian coast, and the first Norwegian record of Fibrocapsa japonica (Raphidophyceae). Harmful Algae 75: 105–117.
  • Eriksen RS, Davies CH, Bonham P, Coman FE, Edgar S, McEnnulty FR, McLeod D, Miller MJ, Rochester W, Slotwinski A, Tonks ML, Uribe-Palomino J, Richardson AJ (2019) Australia’s long-term plankton observations: The integrated marine observing system national reference station network. Frontiers in Marine Science 6: 1–17.
  • Galluzzi L, Bertozzini E, Penna A, Perini F, Pigalarga A, Graneli E, Magnani M (2008) Detection and quantification of Prymnesium parvum (Haptophyceae) by real-time PCR. Letters in Applied Microbiology 46(2): 261–266.
  • Gran-Stadniczeñko S, Supraha L, Egge ED, Edvardsen B (2017) Haptophyte Diversity and Vertical Distribution Explored by 18S and 28S Ribosomal RNA Gene Metabarcoding and Scanning Electron Microscopy. The Journal of Eukaryotic Microbiology 64(4): 514–532.
  • Gran-Stadniczeñko S, Dunthorn Egge E, Hostyeva V, Logares R, Eikrem W, Edvardsen B (2019) Protist Diversity and Seasonal Dynamics in Skagerrak Plankton Communities as Revealed by Metabarcoding and Microscopy. The Journal of Eukaryotic Microbiology 66(3): 494–513.
  • Guo X, Yanagi T (1996) Seasonal variation of residual current in Tokyo Bay, Japan - Diagnostic numerical experiments. Journal of Oceanography 52(5): 597–616.
  • Hallegraeff GM (2004) Harmful algal blooms: a global overview. In: Hallegraeff GM, Andersson DM, Cembella AD (Eds) Manual on Harmful Marine Microalgae. UNSECO, Paris, 25–50.
  • Hattenrath-Lehmann TK, Gobler CJ (2017) Identification of unique microbiomes associated with harmful algal blooms caused by Alexandrium fundyense and Dinophysis acuminata. Harmful Algae 68: 17–30.
  • Hattenrath-Lehmann TK, Jankowiak J, Koch F, Gobler CJ (2019) Prokaryotic and eukaryotic microbiomes associated with blooms of the ichthyotoxic dinoflagellate Cochlodinium (Margalefidinium) polykrikoides in New York, USA, estuaries. PLoS ONE 14(11): e0223067.
  • Herlemann D, Labrenz M, Jürgens K, Bertilsson S, Waniek JJ, Andersson AF (2011) Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. The ISME Journal 5(10): 1571–1579.
  • Hirai J, Katakura S, Kasai H, Nagai S (2017) Cryptic Zooplankton Diversity Revealed by a Metagenetic Approach to Monitoring Metazoan Communities in the Coastal Waters of the Okhotsk Sea, Northeastern Hokkaido. Frontiers in Marine Science 4: 1–13.
  • Holm-Hansen O, Lorenzen CJ, Holmes RW, Strickland JDH (1965) Fluorometric Determination of Chlorophyll. ICES Journal of Marine Science 30(1): 3–15.
  • Ichimi K, Yamasaki M, Okumura Y, Suzuki T (2001) The growth and cyst formation of a toxic dinoflagellate, Alexandrium tamarense, at low water temperatures in northeastern Japan. Journal of Experimental Marine Biology and Ecology 261(1): 17–29.
  • Imai I, Ishida Y, Hata Y (1993) Killing of marine phytoplankton by a gliding bacterium Cytophaga sp., isolated from the coastal sea of Japan. Marine Biology 116(4): 527–532.
  • Imai I, Ishida Y, Sakaguchi K, Hata Y (1995) Algicidal Marine Bacteria Isolated from Northern Hiroshima Bay, Japan. Fisheries Science 61(4): 628–636.
  • Imai I, Yamaguchi M, Hori Y (2006) Eutrophication and occurrences of harmful algal blooms in the Seto Inland Sea, Japan. Plankton & Benthos Research 1(2): 71–84.
  • Inaba N, Trainer VL, Onishi Y, Ishii K, Wyllie ES, Imai I (2017) Algicidal and growth-inhibiting bacteria associated with seagrass and macroalgae beds in Puget Sound, WA, USA. Harmful Algae 62: 136–147.
  • Inaba N, Trainer VL, Nagai S, Kojima S, Sakami T, Takagi S, Imai I (2019) Dynamics of seagrass bed microbial communities in artificial Chattonella blooms: A laboratory microcosm study. Harmful Algae 84: 139–150.
  • Inaba N, Kodama I, Nagai S, Shiraishi T, Matsuno K, Yamaguchi A, Imai I (2020) Distribution of harmful algal growth-limiting bacteria on artificially introduced Ulva and natural macroalgal beds. Applied Sciencea 10(16): 5658.
  • IOC-UNESCO (2021) Oceanic Biogeographic System [WWW Document]. [accessed: 17.11.2021]
  • Jankowiak JG, Gobler CJ (2020) The Composition and Function of Microbiomes Within Microcystis Colonies Are Significantly Different Than Native Bacterial Assemblages in Two North American Lakes. Frontiers in Microbiology 11: 1016.
  • John U, Medlin LK, Groben R (2005) Development of specific rRNA probes to distinguish between geographic clades of the Alexandrium tamarense species complex. Journal of Plankton Research 27(2): 199–204.
  • John U, Litaker RW, Montresor M, Murray S, Brosnahan ML, Anderson DM (2014) Formal revision of the Alexandrium tamarense species complex (Dinophyceae) taxonomy: The introduction of five species with emphasis on molecular-based (rDNA) classification. Protist 165(6): 779–804.
  • Kaga S, Sekiguchi K, Yoshida M, Ogata T (2006) Occurence and toxin production of Alexandrium spp. (Dinophyceae) in coastal waters of Iwate Prefecture, Japan. Nippon Suisan Gakkaishi 72(6): 1068–1076.
  • Kim MJ, Jeong SY, Lee SJ (2008) Isolation, identification, and algicidal activity of marine bacteria against Cochlodinium polykrikoides. Journal of Applied Phycology 20(6): 1069–1078.
  • Klindworth A, Mann AJ, Huang S, Wichels A, Quast C, Waldmann J, Teeling H, Glöckner FO (2014) Diversity and activity of marine bacterioplankton during a diatom bloom in the North Sea assessed by total RNA and pyrotag sequencing. Marine Genomics 18, Part B, 185–192.
  • Koch F, Kang Y, Villareal TA, Anderson DM, Gobler CJ (2014) A novel immunofluorescence flow cytometry technique detects the expansion of brown tides caused by Aureoumbra lagunensis to the Caribbean sea. Applied and Environmental Microbiology 80(16): 4947–4957.
  • Kokubu Y, Yamazaki H, Nagai T, Gross ES (2013) Mixing observations at a constricted channel of a semi-closed estuary: Tokyo Bay. Continental Shelf Research 69: 1–16.
  • Kotani Y, Matsuyama Y, Hayashi M, Matsuoka K (2006) Distribution and Abundance of resting cysts of Alexandrium tamarense and/or A. catenella (Dinophyceae) in Tokyo Bay, Japan. Plankton & Benthos Research 1(3): 147–154.
  • Krawiec RW (1982) Autecology and clonal variability of the marine centric diatom Thalassiosira rotula (Bacillariophyceae) in response to light, temperature and salinity. Marine Biology 69(1): 79–89.
  • Kubo A, Hashihama F, Kanda J, Horimoto-Miyazaki N, Ishimaru T (2019) Long-term variability of nutrient and dissolved organic matter concentrations in Tokyo Bay between 1989 and 2015. Limnology and Oceanography 64(S1): S209–S222.
  • Kudela RM, Gobler CJ (2012) Harmful dinoflagellate blooms caused by Cochlodinium sp.: Global expansion and ecological strategies facilitating bloom formation. Harmful Algae 14: 71–86.
  • Lee CK, Park TG, Park YT, Lim WA (2013) Monitoring and trends in harmful algal blooms and red tides in Korean coastal waters, with emphasis on Cochlodinium polykrikoides. Harmful Algae 30: S3–S14.
  • Lefebvre A, Guiselin N, Barbet F, Artigas FL (2011) Long-term hydrological and phytoplankton monitoring (1992–2007) of three potentially eutrophic systems in the eastern English Channel and the Southern Bight of the North Sea. ICES Journal of Marine Science 68(10): 2029–2043.
  • Lim HC, Teng ST, Leaw CP, Lim PT (2013) Three novel species in the Pseudo-nitzschia pseudodelicatissima complex: P. batesiana sp. nov., P. lundholmiae sp. nov., and P. fukuyoi sp. nov. (Bacillariophyceae) from the Strait of Malacca, Malaysia. Journal of Phycology 49(5): 902–916.
  • Lima-Mendez G, Faust K, Henry N, Decelle J, Colin S, Carcillo F, Chaffron S, Ignacio-Espinosa JC, Roux S, Vincent F, Bittner L, Darzi Y, Wang J, Audic S, Berline L, Bontempi G, Cabello BM, Coppola L, Cornejo-Castillo FM, d’Ovidio F, De Meester L, Ferrera I, Garet-Delmas M-J, Guidi L, Lara E, Pesant S, Royo-Llonch M, Salazar G, Sánchez P, Sebastian M, Souffreau C, Dimier C, Picheral M, Searson S, Kandels-Lewis S, Tara Oceans coordinators , Gorsky G, Not F, Ogata H, Speich S, Stemann L, Weissenbach J, Wincker P, Acinas SG, Sunagawa S, Bork P, Sullivan MB, Karsenti E, Bowler C, de Vargas C, Raes J (2015) Determinants of community structure in the global plankton interactome. Science 348(6237): 1262073.
  • Litaker RW, Fraga S, Montresor M, Brosnahan ML, Hoppenrath M, Murray S, Wolny J, John U, Sampedro N, Larsen J, Calado AJ (2018) A practical guide to new nomenclature for species within the “Alexandrium tamarense species complex”. Harmful Algae News 61: 13–15.
  • Lovejoy C, Bowman JP, Hallegraeff GM (1998) Algicidal Effects of a Novel Marine Pseudoalteromonas Isolate (Class Proteobacteria, Gamma Subdivision) on Harmful Algal Bloom Species of the Genera Chattonella, Gymnodinium, and Heterosigma. Applied and Environmental Microbiology 64(8): 2806–2813.
  • Lundholm N, Skov J, Pocklington R, Moestrup Ø (1997) Studies on the marine planktonic diatom Pseudo-nitzschia. 2. Autecology of P. pseudodelicatissima based on isolates from Danish coastal waters. Phycologia 36(5): 381–388.
  • Luo Z, Krock B, Neil K, Nézan E, Chomérat N, Bilien G, Tillmann U, Gu H (2017) Adding new pieces to the Azadinium (Dinophyceae) diversity and biogeography puzzle: Non-toxigenic Azadinium zhuanum sp. nov. from China, toxigenic A. poporum from the Mediterranean, and a non-toxigenic A. dalianense from the French Atlantic. Harmful Algae 66: 65–78.
  • Majaneva M, Rintala JM, Piisilä M, Fewer DP, Blomster J (2012) Comparison of wintertime eukaryotic community from sea ice and open water in the Baltic Sea, based on sequencing of the 18S rRNA gene. Polar Biology 35(6): 875–889.
  • Martin-Cereceda M, Roberts EC, Wootton EC, Bonaccorso E, Dyal P, Guinea A, Rogers D, Wright CJ, Novarino G (2010) Morphology, ultrastructure, and small subunit rDNA phylogeny of the marine heterotrophic flagellate Goniomonas aff. amphinema. The Journal of Eukaryotic Microbiology 57(2): 159–170.
  • Matsuoka K, Joyce LB, Kotani Y, Matsuyama Y (2003) Modern dinoflagellate cysts in hypertrophic coastal waters of Tokyo Bay, Japan. Journal of Plankton Research 25(12): 1461–1470.
  • Meyer N, Pohnert G (2019) Isolate-specific resistance to the algicidal bacterium Kordia algicida in the diatom Chaetoceros genus. Botanica Marina 62(6): 527–535.
  • Meyer N, Bigalke A, Kaulfuß A, Pohnert G (2017) Strategies and ecological roles of algicidal bacteria. FEMS Microbiology Reviews 41(6): 880–899.
  • Mikhailov IS, Bukin YS, Zakharova YR, Usoltseva MV, Galachyants YP, Sakirko MV, Blinov VV, Likhoshway YV (2019a) Co-occurrence patterns between phytoplankton and bacterioplankton across the pelagic zone of Lake Baikal during spring. Journal of Microbiology (Seoul, Korea) 57(4): 252–262.
  • Mikhailov IS, Zakharova YR, Bukin YS, Galachyants YP, Petrova DP, Sakirko MV, Likhoshway YV (2019b) Correction to: Co-occurrence Networks Among Bacteria and Microbial Eukaryotes of Lake Baikal During a Spring Phytoplankton Bloom. Microbial Ecology 77(1): 96–109.
  • Miyahara K, Nagai S, Itakura S, Yamamoto K, Fujisawa K, Iwamoto T, Yoshimatsu S, Matsuoka S, Yuasa A, Makino K, Hori Y, Nagata S, Nagasaki K, Yamaguchi M, Honjo T (1996) First Record of a Bloom of Thalassiosira diporocyclus in the Eastern Seto Inland Sea. Fisheries Science 62(6): 878–882.
  • Moestrup Ø, Akselmann R, Fraga S, Hoppenrath M, Iwataki M, Komárek J, Larsen J, Lundholm N, Zingone A (2022) IOC-UNESCO Taxonomic Reference List of Harmful Micro Algae. [accessed 21.11.21]
  • Moreno-Pino M, Krock B, De la Iglesia R, Echenique-Subiabre I, Pizarro G, Vásquez M, Trefault N (2018) Next Generation Sequencing and mass spectrometry reveal high taxonomic diversity and complex phytoplankton-phycotoxins patterns in Southeastern Pacific fjords. Toxicon 151: 5–14.
  • Moro I, La Rocca N, Valle LD, Moschin E, Negrisolo E, Andreoli C (2011) Pyramimonas australis sp. nov. (Prasinophyceae, Chlorophyta) from Antarctica: Fine structure and molecular phylogeny. European Journal of Phycology 37(1): 103–114.
  • Nagai S (2018) Marine eukaryote and HAB monitoring in Japan with next generation technology. Harmful Algae News 60: 1–3.
  • Nagai S, Imai I (1999) Possibility for bio-control of harmful diatom blooms in Coscinodiscus wailesii by marine bacteria. Microbes and Environments 14(4): 253–262.
  • Nagai S, Yamamoto K, Hata N, Itakura S (2012) Study of DNA extraction methods for use in loop-mediated isothermal amplification detection of single resting cysts in the toxic dinoflagellates Alexandrium tamarense and A. catenella. Marine Genomics 7: 51–56.
  • Nagai S, Hida K, Urusizaki S, Takano Y, Hongo Y, Kameda T, Abe K (2016) Massively parallel sequencing-based survey of eukaryotic community structures in Hiroshima Bay and Ishigaki Island. Gene 576(2): 681–689.
  • Nagai S, Urusizaki S, Hongo Y, Chen H, Dzhembekova N (2017) An attempt to semi-quantify potentially toxic diatoms of the genus Pseudo-nitzschia in Tokyo Bay, Japan by using massively parallel sequencing technology. Plankton & Benthos Research 12(4): 248–258.
  • Nagai S, Chen H, Kawakami Y, Yamamoto K, Sildever S, Kanno N, Oikawa H, Yasuike M, Nakamura Y, Hongo Y, Fujiwara A, Kobayashi T, Gojobori T (2019) Monitoring of the toxic dinoflagellate Alexandrium catenella in Osaka Bay, Japan using a massively parallel sequencing (MPS)-based technique. Harmful Algae 89: 101660.
  • Nakane T, Nakaka K, Bouman H, Platt T (2008) Environmental control of short-term variation in the plankton community of inner Tokyo Bay, Japan. Estuarine, Coastal and Shelf Science 78(4): 796–810.
  • Natsuike M, Kanamori M, Shimada H (2019) Red tide and seasonal occurence of the harmful raphidophyte Heterosigma akashiwo in Hakodate Bay, Hokkaido. Scientific Reports of Hokkaido Fisheries Research Insitutes 95: 11–17.
  • Nihei Y, Shigeta K, Ito M, Hoshino A, Fukuda M, Kato Y (2009) Sediment Transports and Qualities of Sediment in Rivers Flowing into Tokyo Bay. Journal of Japan Society of Civil Engineers, Ser. B2 (Coastal Engineering) 65(1): 1171–1175. [In Japanese]
  • Nishikawa T, Hori Y, Nagai S, Miyahara K, Nakamura Y, Harada K, Tanda M, Manabe T, Tada K (2010) Nutrient and Phytoplankton Dynamics in Harima-Nada, Eastern Seto Inland Sea, Japan During a 35-Year Period from 1973 to 2007. Estuaries and Coasts 33(2): 417–427.
  • Nishikawa T, Hori Y, Nagai S, Miyahara K, Nakamura Y, Harada K, Tada K, Imai I (2011) Long time-series observations in population dynamics of the harmful diatom Eucampia zodiacus and environmental factors in Harima-Nada, eastern Seto inland Sea, Japan during 1974–2008. Plankton & Benthos Research 6(1): 26–34.
  • Nomura H (1998) Changes in Red Tide Events and Phytoplankton Community Composition in Tokyo Bay from the Historical Plankton Records in a Period between 1907 and 1997. Oceanography in Japan 7(3): 159–178.
  • Nomura H, Yoshida M (1997) Recent occurence of phytoplankton in hyper-eutrophicated inlet, Tokyo Bay, central Japan. La Mer (Paris) 35: 107–121.
  • Oda R, Kanda M (2009) Observed sea surface temperature of Tokyo Bay and its impact on urban air temperature. Journal of Applied Meteorology and Climatology 48(10): 2054–2068.
  • Ogawa DM, Moriya S, Tsuboi Y, Date Y, Rieto-da-Silva ARB, Radis-Baptista G, Yamane T, Kikuchi J (2014) Biogeochemical Typing of Paddy Field by a Data-Driven Approach Revealing Sub-Systems within a Complex Environment - A Pipeline to Filtrate, Organize and Frame Massive Dataset from Multi-Omics Analyses. PLoS ONE 9(10): e110723.
  • Oita A, Tsuboi Y, Date Y, Oshima T, Sakata K, Yokoyama A, Moriya S, Kikuchi J (2018) Profiling physicochemical and planktonic features from discretely/continuously sampled surface water. The Science of the Total Environment 636: 12–19.
  • Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, Stevens MHH, Szoecs E, Wagner H (2021) vegan: Community Ecology Package. R package version 2.5-6. [accessed 21.11.2020]
  • Orita K (2016) Fish-kill by Pseudochattonella verruculosa (Dictyochophyceae). In: Imai I, Yamaguchi M, Matsuoka K (Eds) Advances in Harmful Algal Bloom Research. Kouseisha Co. Lmt., Tokyo, 226–231.
  • Paver SF, Hayek KR, Gano KA, Fagen JR, Brown CT, Davis-Richardson AG, Crabb DB, Rosario-Passapera R, Giongo A, Triplett EW, Kent AD (2013) Interactions between specific phytoplankton and bacteria affect lake bacterial community succession. Environmental Microbiology 15(9): 2489–2504.
  • Prasad AKSK, Nienow JA, Lochner E (2018) Thalassiosira mala (Bacillariophyta), a potentially harmful, marine diatom from Chilka Lake and other coastal localities of Odisha, India: Nomenclature, frustule morphology and global biogeography. Journal of Biosciences 43(1): 59–74.
  • R Core Team (2021) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. [accessed 21.11.2020]
  • Razali RM, Leaw C-P, Lim H, Ayub MNA, Rahim MA, Lim P-T (2016) First Report of a Marine Dinoflagellate, Alexandrium andersonii (Dinophyaceae) in Malaysian Waters. Malaysian Journal of Science 35: 304–315.
  • Revelle W (2020) psych: Procedures for Personality and Psychological Research. Northwestern University, Evanston, Illinois. R package version 2.0.9. [accessed 21.11.2020]
  • Rhodes LL (1998) Identification of potentially toxic Pseudo-nitzschia (Bacillariophyceae) in New Zealand coastal waters, using lectins. New Zealand Journal of Marine and Freshwater Research 32(4): 537–544.
  • Riemann L, Steward GF, Azam F (2000) Dynamics of bacterial community composition and activity during a mesocosm diatom bloom. Applied and Environmental Microbiology 66(2): 578–587.
  • Rodriguez-Ramos T, Dornelas M, Maranon E, Cermeno P (2014) Conventional sampling methods severely underestimate phytoplankton species richness. Journal of Plankton Research 36(2): 334–343.
  • Ruggiero MV, Kooistra WHCF, Piredda R, Sarno D, Zampicinini G, Zingone A, Montresor M (2022) Temporal changes of genetic structure and diversity in a marine diatom genus discovered via metabarcoding. Environmental DNA 00(4): 1–13.
  • Sanseverino I, Conduto D, Pozzoli L, Dobricic S, Lettieri T (2016) Algal bloom and its economic impact. Joint Research Centre Eur 27905 En.
  • Sassenhagen I, Irion S, Jardillier L, Moreira D, Christaki U (2020) Protist Interactions and Community Structure During Early Autumn in the Kerguelen Region (Southern Ocean). Protist 171(1): 125709.
  • Sawaya NA, Djurhuus A, Closek CJ, Hepner M, Olesin E, Visser L, Kelble C, Hubbard K, Breitbart M (2019) Assessing eukaryotic biodiversity in the Florida Keys National Marine Sanctuary through environmental DNA metabarcoding. Ecology and Evolution 9(3): 1029–1040.
  • Sekiguchi K, Watanabe S, Shimizu M, Saito S (1986) Occurence of Protogonyaulax tamarensis on the coast of Iwate Prefecture, in relation to toxification of the scallop Patinopecten yessoensis. Bulletin of Tohoku Regional Fisheries Research Laboratory 48: 115–123.
  • Sekiyama Y, Chikayama E, Kikuchi J (2011) Evaluation of a Semipolar Solvent System as a Step toward Heteronuclear Multidimensional NMR-Based Metabolomics for 13 C-Labeled Bacteria, Plants, and Animals. Analytical Chemistry 83(3): 719–726.
  • Shi X, Liu L, Li Y, Xiao Y, Ding G, Lin S, Chen J (2018) Isolation of an algicidal bacterium and its effects against the harmful-algal- bloom dinoflagellate Prorocentrum donghaiense (Dinophyceae). Harmful Algae 80: 72–79.
  • Shimada H, Kanamori M, Yoshida H, Imai I (2016) First record of red tide due to the harmful dinoflagellate Karenia mikimotoi in Hakodate Bay, southern Hokkaido, in autumn 2015. Nippon Suisan Gakkaishi 82(6): 934–938.
  • Shin H, Lee E, Shin J, Ko S-R, Oh H-S, Ahn C-Y, Oh H-M, Cho B-K, Cho S (2018) Elucidation of the bacterial communities associated with the harmful microalgae Alexandrium tamarense and Cochlodinium polykrikoides using nanopore sequencing. Scientific Reports 8(1): 5323.
  • Sigsgaard EE, Nielsen IB, Carl H, Krag MA, Knudsen SW, Xing Y, Holm-Hansen TH, Møller PR, Thomsen PF (2017) Seawater environmental DNA reflects seasonality of a coastal fish community. Marine Biology 164(6): 1–15.
  • Sildever S, Kawakami Y, Kanno N, Kasai H, Shiomoto A, Katakura S, Nagai S (2019) Toxic HAB species from the Sea of Okhotsk detected by a metagenetic approach, seasonality and environmental drivers. Harmful Algae 87: 101631.
  • Sildever S, Laas P, Kolesova N, Lips I, Lips U, Nagai S (2021) Plankton biodiversity and species co-occurrence based on environmental DNA – a multiple marker study. Metabarcoding and Metagenomics 5: 175–197.
  • Sison-Mangus MP, Jiang S, Kudela RM, Mechic S (2016) Phytoplankton-Associated Bacterial Community Composition and Succession during Toxic Diatom Bloom and Non-Bloom Events. Frontiers in Microbiology 7: 1433.
  • Skerratt JH, Bowman JP, Hallegraeff G, James S, Nichols PD (2002) Algicidal bacteria associated with blooms of a toxic dinoflagellate in a temperate Australian estuary. Marine Ecology Progress Series 244: 1–15.
  • Skjelbred B, Edvardsen B, Andersen T (2013) Environmental optima for seven strains of Pseudochattonella (Dictyochophyceae, Heterokonta). Journal of Phycology 49(1): 54–60.
  • Śliwińska-Wilczewska S, Latała A (2018) Allelopathic activity of the bloom-forming picocyanobacterium Synechococcus sp. on the coexisting microalgae: The role of eutrophication. International Review of Hydrobiology 103(3–4): 37–47.
  • Śliwińska-Wilczewska S, Felpeto AB, Mozdzeń K, Vasconcelos V, Latała A (2019) Physiological effects on coexisting microalgae of the allelochemicals produced by the bloom-forming cyanobacteria Synechococcus sp. and Nodularia spumigena. Toxins 11(12): 712.
  • Smith KF, Kohli GS, Murray SA, Rhodes LL (2017) Assessment of the metabarcoding approach for community analysis of benthic-epiphytic dinoflagellates using mock communities. New Zealand Journal of Marine and Freshwater Research 51(4): 555–576.
  • Sourisseau M, Le Guennec V, Le Gland GL, Plus M, Chapelle A (2017) Resource competition affects plankton community structure: Evidence from trait-based modeling. Frontiers in Marine Science 4: 1–14.
  • Stoeckle MY, Soboleva L, Charlop-Powers Z (2017) Aquatic environmental DNA detects seasonal fish abundance and habitat preference in an urban estuary. PLoS ONE 12(4): 1–15.
  • Suzuki T, Matsuyama M (2000) Numerical experiments on stratified wind-induced circulation in Tokyo Bay, Japan. Estuarine, Coastal and Shelf Science 50(1): 17–25.
  • Sze Y, Miranda LN, Sin TM, Huang D (2019) Characterising planktonic dinoflagellate diversity in Singapore using DNA metabarcoding. Metabarcoding and Metagenomics 2: 1–14.
  • Takano H (1956) Harmful Blooming of Minute Cells in of Thalassiosira decipiens in Coastal Water in Tokyo Bay. Journal of the Oceanographical Society of Japan 12(2): 64–67.
  • Takao T, Okada T, Nakayama K, Furukawa K (2004) Residence time of sea water in Tokyo Bay at 2002 and reproducing residence time with multi-box vertical one-dimensional model. Proceedings of Hydraulic Engineering 48: 1243–1248.
  • Teeling H, Fuchs BM, Becher D, Klockow C, Gardebrecht A, Bennke CM, Kassabgy M, Huang S, Mann AJ, Waldmann J, Weber M, Klindworth A, Otto A, Lange J, Bernhardt J, Reinsch C, Hecker M, Peplies J, Bockelmann FD, Callies U, Gerdts G, Wichels A, Wiltshire KH, Glöckner FO, Schweder T, Amann R (2012) Substrate-Controlled Succession of Marine Bacterioplankton Populations Induced by a Phytoplankton Bloom. Science 336(6081): 608–612.
  • Tian Y, Huang M (2019) An integrated web-based system for the monitoring and forecasting of coastal harmful algae blooms: Application to Shenzhen city, China. Journal of Marine Science and Engineering 7(9): 1–17.
  • Toba M, Kosemura T, Yamakawa H, Sugiura Y, Kobayashi Y (2008) Field and laboratory observations on the hypoxic impact on survival and distribution of short-necked clam Ruditapes philippinarum larvae in Tokyo Bay, central Japan. Plankton & Benthos Research 3(3): 165–173.
  • Viana TV, Fistarol GO, Amario M, Menezes RB, Carneiro BLR, Chaves DM, Hargreaves PI, Silva-Lima AW, Valentin JL, Tenenbaum DR, Arruda EF, Paranhos R, Salomon PS (2019) Massive blooms of Chattonella subsalsa Biecheler (Raphidophyceae) in a hypereutrophic, tropical Estuary-Guanabara Bay, Brazil. Frontiers in Marine Science 6: 1–14.
  • Vrieling EG, Koeman RPT, Nagasaki K, Ishida Y, Peperzak L, Gieskes WWG, Veenhuis M (1995) Chattonella and Fibrocapsa (Raphidophyceae): First Observation of, Potentially Harmful, Red Tide Organisms in Dutch Coastal Waters. Netherlands Journal of Sea Research 33(2): 183–191.
  • Wang H, Lu D, Huang H, Göbel J, Dai X, Xia P (2011) First observation of Karlodinium veneficum from the East China Sea and the coastal waters of Germany. Acta Oceanologica Sinica 30(6): 112–121.
  • Warnes GR, Bolker B, Bonebakker L, Gentleman R, Huber W, Liaw A, Lumley T, Maechler M, Magnusson A, Moeller S, Schwartz M, Venables B (2020) gplots: Various R Programming Tools for Plotting Data. R package version 3.0.4. [accessed: 23.11.2021]
  • Weiss S, Van Treurenv W, Lozupone C, Faust K, Friedman J, Deng Y, Xia LC, Xu ZZ, Ursell L, Alm EJ, Birmingham A, Cram JA, Fuhrman JA, Raes J, Sun F, Zhou J, Knight R (2016) Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. The ISME Journal 10(7): 1669–1681.
  • Whyte JNC, Haigh N, Ginther NG, Keddy LJ (2001) First record of blooms of Cochlodinium sp. (Gymnodiniales, Dinophyceae) causing mortality to aquacultured salmon on the west coast of Canada. Phycologia 40(3): 298–304.
  • Wurzbacher C, Attermeyer K, Kettner MT, Flintrop C, Warthmann N, Hilt S, Grossart HP, Monaghan MT (2017) DNA metabarcoding of unfractionated water samples relates phyto-, zoo- and bacterioplankton dynamics and reveals a single-taxon bacterial bloom. Environmental Microbiology Reports 9(4): 383–388.
  • Xiao X, Sogge H, Lagesen K, Tooming-Klunderud A, Jakobsen KS, Rohrlack T (2014) Use of high throughput sequencing and light microscopy show contrasting results in a study of phytoplankton occurrence in a freshwater environment. PLoS ONE 9(8): e106510.
  • Yamaguchi M, Itakura S, Nagasaki K, Matsuyama Y (1997) Effects of temperature and salinity on the growth of the red tide flagellates Heterocapsa circularisquama (Dinophyceae) and Chattonella verruculosa (Raphidophyceae). Journal of Plankton Research 19(8): 1167–1174.
  • Yamaguchi H, Minamida M, Matsubara T, Okamura K (2014) Novel blooms of the diatom Asteroplanus karianus deplete nutrients from Ariake Sea coastal waters. Marine Ecology Progress Series 517: 51–60.
  • Yamamoto K, Nakajima M, Imai I (2017) Expansion of blooming in the toxic dinoflagellate Alexandrium tamarense and environmental fluctuation analyzed from long-term monitoring data in Osaka Bay, eastern Seto Inland Sea, Japan. Bulletin of the Plankton Society of Japan 64: 11–21.
  • Yang C, Li Y, Zhou B, Zhou Y, Zheng W, Tian Y, Van Nostrand JD, Wu L, He Z, Zhou J, Zheng T (2015) Illumina sequencing-based analysis of free-living bacterial community dynamics during an Akashiwo sanguine bloom in Xiamen sea, China. Scientific Reports 5(1): 1–11.
  • Yap-Dejeto LG, Omura T, Nagahama Y, Fukuyo Y (2010) Observations of eleven Pseudo-nitzschia species in Tokyo Bay, Japan. La Mer (Paris) 48: 1–16.
  • Yarimizu K, Fujiyoshi S, Kawai M, Norambuena-Subiabre L, Cascales EK, Rilling JI, Vilugrón J, Cameron H, Vergara K, Morón-López J, Acuña JJ, Gajardo G, Espinoza-González O, Guzmán L, Jorquera MA, Nagai S, Pizarro G, Riquelme C, Ueki S, Maruyama F (2020) Protocols for monitoring harmful algal blooms for sustainable aquaculture and coastal fisheries in Chile. International Journal of Environmental Research and Public Health 17(20): 1–24.
  • Yasakova ON (2013) The seasonal dynamics of potentially toxic and harmful phytoplankton species in Novorossiysk Bay (Black Sea). Russian Journal of Marine Biology 39(2): 107–115.
  • Yuki K, Yoshimatsu S (2012) Occurence of Alexandrium minutum Halim and Alexandrium ostefeldii (Paulsen) Balech et Tangen (Dinophyceae) in Yashima Bay, eastern Seto Inland Sea, Japan. Bulletin of Akashiwo Research Institute of Kagawa Prefecture 8: 1–6.
  • Zhang QC, Qiu LM, Yu RC, Kong FZ, Wang YF, Yan T, Gobler CJ, Zhou MJ (2012) Emergence of brown tides caused by Aureococcus anophagefferens Hargraves et Sieburth in China. Harmful Algae 19: 117–124.
  • Zhang Y, Pavlovska M, Stoica E, Prekrasna I, Yang J, Slobodnik J, Zhang X, Dykyi E (2020) Holistic pelagic biodiversity monitoring of the Black Sea via eDNA metabarcoding approach: From bacteria to marine mammals. Environment International 135: 105307.
  • Zheng N, Ding N, Gao P, Han M, Liu X, Wang J, Sun L, Fu B, Wang R, Zhou J (2018) Diverse algicidal bacteria associated with harmful bloom-forming Karenia mikimotoi in estuarine soil and seawater. The Science of the Total Environment 631–632: 1415–1420.
  • Zhou J, Richlen ML, Sehein TR, Kulis DM, Anderson DM, Cai Z (2018) Microbial Community Structure and Associations During a Marine Dinoflagellate Bloom. Frontiers in Microbiology 9: 1201.
  • Zingone A, Siano R, Alelio DD, Sarno D (2006) Potentially toxic and harmful microalgae from coastal waters of the Campania region (Tyrrhenian Sea, Mediterranean Sea). Harmful Algae 5(3): 321–337.

Supplementary materials

Supplementary material 1 

Figures S1–S6

Sirje Sildever, Noriko Nishi, Nobuharu Inaba, Taiga Asakura, Jun Kikuchi, Yasuhito Asano, Takanori Kobayashi, Takashi Gojobori, Satoshi Nagai

Data type: Species data, statistical data visualization

Explanation note: Overview of HAB species occurrences based on light-microscopy detection; ordination analysis (PCA), correlation analyses (heatmaps).

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (9.85 MB)
Supplementary material 2 

Tables S1–S13

Sirje Sildever, Noriko Nishi, Nobuharu Inaba, Taiga Asakura, Jun Kikuchi, Yasuhito Asano, Takanori Kobayashi, Takashi Gojobori, Satoshi Nagai

Data type: OTU data, statistical data, species data

Explanation note: OTU information used in this study, overview of harmful algal bloom species in Tokyo Bay,overview of growth-inhibiting bacterial genera, correlation analysis results, association-based time-series analysis results.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (129.96 kb)
login to comment