Research Article |
Corresponding author: Satoshi Nagai ( snagai@affrc.go.jp ) Academic editor: Florian Leese
© 2022 Sirje Sildever, Noriko Nishi, Nobuharu Inaba, Taiga Asakura, Jun Kikuchi, Yasuhito Asano, Takanori Kobayashi, Takashi Gojobori, Satoshi Nagai.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Sildever S, Nishi N, Inaba N, Asakura T, Kikuchi J, Asano Y, Kobayashi T, Gojobori T, Nagai S (2022) Monitoring harmful microalgal species and their appearance in Tokyo Bay, Japan, using metabarcoding. Metabarcoding and Metagenomics 6: e79471. https://doi.org/10.3897/mbmg.6.79471
|
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.
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 (
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 (
Furthermore, eDNA metabarcoding has been successfully applied to reveal appearance patterns for several organism groups ranging from bacteria to mammals (
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 (
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 (
Microalgal abundance is the highest between March and October and the community is dominated by diatoms (
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.
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.
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
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,
The demultiplexing and trimming were performed using Trimmomatic version 0.35 (http://www.usadellab.org/cms/?page=trimmomatic). 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 (http://www.drive5.com/usearch/) with default settings (≥ 16 bp overlap, ≥ 90% similarity and mismatch ≤ 5 bp; http://www.drive5.com/usearch/manual/merge_options_html) resulting in a maximum sequence length of 584 bp. Singletons were also removed. Sequences were aligned using Clustal Omega v. 1.2.0. (http://www.clustal.org/omega/) 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 (http://www.mothur.org./wiki/MiSeq_SOP;
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 (
All data analysis was conducted in R v. 4.0.3 (
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 data analysis (
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.
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
In 2017, 673 unique eukaryotic OTUs (≥ 99% BLAST top-hit similarity) were detected, belonging to nine different supergroups and 32 phyla (Suppl. material
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
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
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.
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 (
Ordination
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
Redundancy analysis, based on HAB-associated OTUs detected by metabarcoding and HTS approach, A. 2017; B. 2018.
Correlation
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
Eukaryotes
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
Prokaryotes
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
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.
Based on the time-series analysis of selected phyla, several associations between bacteria and phytoplankton were detected (Suppl. material
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 (
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 (
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 (
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 (
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 (
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 (
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;
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 (
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 (
Another important aspect of structuring a species’ appearance is interaction with other organisms (
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 (
Bacterial communities are also influenced by the phytoplankton blooms and species composition (
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 (
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].
Figures S1–S6
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).
Tables S1–S13
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.