Genomic analysis of red-tide water bloomed with Heterosigma akashiwo in Geoje

Phytoplankton-specific universal primer set optimized for the MiSeq platform

A universal primer set was designed to increase specificity as well as taxon coverage of phytoplankton (Table 1). A total of 1,473 23S rDNA sequences (997 from proteobacteria and 476 from phytoplankton and cyanobacteria) obtained from the public databases including GenBank (https://www.ncbi.nlm.nih.gov/genbank/) & BOLD (http://www.barcodinglife.org) were compared using the Clustal omega program (http://www.ebi.ac.uk/Tools/msa/clustalo/). A new 23S universal forward primer (P23MISQF1) was designed by several modifications of the previously designed ones (Table S1). Briefly, we added the degenerate sequence at the seventh position from the 5′ end of P23MISQF1 to increase the taxon coverage for the Heterokonts (A/T), which was previously adenine (A) in the A23SrVF1 (Yoon et al., 2016; hereafter referred to as Yoon’s 23S universal primer) and p23SrV-f1 primers (Sherwood & Presting, 2007; hereafter referred to as Sherwood’s 23S universal primer). To increase specificity for phytoplankton, guanine (G) was added to the 3′ end of P23MISQF1 (Table S1). We also introduced two changes in the reverse primer, P23MISQR1 (Table 1). First, the nitrogenous base in the fifth nucleotide from its 5′ end was replaced by pyrimidine (Y) bases to increase the coverage of algal sequences (Table S2). The second modification was the addition of two nucleotides at the 3′ end of P23MISQR1 (Table S2), which increased both its melting temperature (T_m) and the sequence specificity for phytoplankton sequences. The expected sizes of amplified products by the newly modified 23S universal primer set (hereafter referred to as Kang’s 23S universal primer set) ranged from 407 to 414 bps, which is optimized for the MiSeq platform. In order to evaluate the designed plastid 23S universal primer set, In-silico PCR was performed on LSU-132 database (RefNR sequence collection) using the SILVA TestPrime tool with zero mismatches (https://www.arb-silva.de/search/testprime/).

Sample collection and DNA extraction

We tested the reliability of Kang’s 23S universal primer set using two seawater samples collected from the East/Japan Sea in 2014 as part of the “Long-term change of structure and function in marine ecosystems of Korea” project funded by the Ministry of Oceans and Fisheries, Korea. To analyze the bloom, water samples were collected on Aug 20, 2015 from three sample sites in Geoje, Korea (Fig. 1). A water sample collected from a site distantly located from the bloom (N34°82′825″, E128°56′533″) was used as the control. Two water samples were collected from the center of the bloom (N34°80′110″, E128°52′572″) and at its edge site, which was close to the bloom, but did not exhibit a water-color change (E site; N34°80′245″, E128°52′175″). From each sample site, 1L of surface water was collected and stored in an ice bucket before filtration through a 0.45 µm GH polypro membrane filter (Pall Corporation, New York, NY, USA). The membrane filters were then cut into small pieces using autoclaved dissecting scissors and completely grinded in a mortar and pestle with liquid nitrogen. Genomic DNA was extracted using a DNeasy^® plant mini kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. Isolated genomic DNA was quantified using an ND-1000 nanodrop spectrophotometer (Thermo Scientific, Waltham, MA, USA) and stored at −70 °C until used for library construction.

Library preparation and sequencing

Isolated genomic DNA was used as a template for the library construction of MiSeq sequencing. Libraries derived from the 16S (Bakt_341F and Bakt_805R) and Kang’s 23S universal primer sets (P23MISQF1 and P23MISQR1) were used for microbial and phytoplankton communities, respectively (Table 1). Additionally, two 23S universal primer sets (Sherwood’s and Yoon’s 23S universal primer sets) were used to test the reliability of Kang’s 23S universal primer set for phytoplankton community analyses (Table 1). The library was constructed using the Nextera XT index kit (Illumina, San Diego, CA, USA) according to the manufacturer’s manual. First, PCR amplification was done using the universal primer sets (NexBakt_341F and NexBakt_805R, NexP23MISQF1 and NexP23MISQR1), which overhang the adapter sequence on forward and reverse primers, respectively (Table 1). The PCR reaction (total volume 20 µL) contained 10 ng of template, 1 µL of each primer (10 pmol), 2 µL of dNTPs (10 mM), 0.2 µL Phusion High Fidelity DNA polymerase (New England Biolabs, Hitchin, UK), and 4 µL 5X buffer. The first PCR condition was an initial denaturation at 94 °C for 3 min, followed by 15 cycles at 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 3 min. The PCR products from the first amplification were purified using the AccuPrep^® PCR purification Kit (Bioneer, Daejeon, Republic of Korea) and eluted with 20 µL of elution buffer. The same conditions including PCR cycles and volume of components were employed for the second PCR amplification, except that 4 µL of purified first PCR product was used as a template and the indexing primers for the MiSeq platform. The second PCR amplicon was separated by 1.5% agarose gel electrophoresis and stained with loading star dye (Dynebio, Seoul, Republic of Korea). PCR products with the expected sizes (approx. 580 bp for analysis of 16s rRNA sequences and approximately 540 bp for analysis of 23s rRNA sequences) were cut from the gel and purified using an AccuPrep^® gel purification Kit (Bioneer, Daejeon, Republic of Korea). The quality and quantity of the libraries were measured using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Finally, constructed libraries were loaded with a MiSeq 600-cycle Reagent Kit v3 (Illumina, San Diego, CA, USA) to perform 300-bp paired-end sequencing on a MiSeq instrument.

Bioinformatics analysis of NGS data

The raw reads with a low quality (QV < 20) and shorter than 100 nucleotides were eliminated from further analysis using the CLC Genomic Workbench v.8.0 (CLC Bio, Cambridge, MA, USA). The reads were merged with longer than 6 bp overlapping sequences without any mismatches. The merged read with the expected size ranges (400∼500 for 16S and 350∼450 for 23S) were selected and their primer sequences were trimmed using Mothur software v.1.35.0 (Schloss et al., 2009). The obtained merged reads were clustered into operational taxonomic units (OTUs) at 99% similarities and chimeras were removed using UCHIME software v.8.1 (Edgar et al., 2011). Operational taxonomic units (OTUs) with less than 10 merged reads or below 0.1% of the total merged reads were eliminated from further analysis. The species name for each OTU was assigned by the similarity search using a blastn search of BLAST +2.2.30 (Camacho et al., 2009) on the NCBI non-redundant nucleotide database (ftp://ftp.ncbi.nlm.nih.gov/blast/db/; accession date: 04/04/2017). Top-scored species name was assigned for each OTU with higher than 98% sequence identity to the database. The OTUs with 90–98% identities in the database were described as “Genus name with highest score” followed by “sp.” OTUs with less than 90% identity were classified as “Unknown”. A phylogenetic tree was constructed by the Minimum Evolution algorithm using Molecular Evolutionary Genetics Analysis (MEGA ver 6.0) (Tamura et al., 2013).

Quantitative PCR for microbial communities

To quantify total microorganisms and phytoplankton communities, qPCR with two different universal primer sets (Bakt_341F and Bakt_805R and P23MISQF1 and P23MISQR1, respectively) were employed (Table 1). It was performed using a DNA Engine Chromo 4 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) under the following conditions: initial template denaturation (94 °C for 3 min); 40 amplification cycles (94 °C for 30 s; 55 °C for 30 s; 72 °C for 30 s) and a final extension step 72 °C for 3 min. A 20 µL volume of the qPCR mixture contained 10 µl of 2 X SYBR Green premix Ex Taq II (Takara Bio Inc., Kuratsu, Japan), 4 µl of template, 1 µl of forward and reverse primers (10 pmol), and 4 µl of purified PCR grade water. Standard curves were constructed to confirm the efficiency of each primer set and quantify copy numbers.

Comparative analysis of phytoplankton community structures generated by three 23S universal primer sets

As the result of SILVA TestPrime tool, all matched sequences (7,749) were photosynthetic groups without any heterotrophic bacterial sequences indicating high specificity of Kang’s 23S universal primer set to photosynthetic phytoplankton. To determine the reliability of the modified Kang’s 23S universal primer set for the analysis of the phytoplankton community, the NGS results of the same seawater sample with three different 23S universal primer sets (Sherwood’s, Yoon’s, and Kang’s) were compared (Table 2). After trimming and clustering the raw reads, 103,359, 175,854, and 54,129 of the merged reads were obtained by Sherwood’s, Yoon’s, and Kang’s 23S primer sets, respectively (Table 2). The highest OTU numbers were identified in the results of Kang’s primer set (98 OTUs) followed by Yoon’s (81 OTUs) and Sherwood’s primer set (28 OTUs). Only 16 eukaryotic algal OTUs were obtained using Sherwood’s 23S primers, while 60 and 67 phytoplankton OTUs (one cyanobacteria and 59 eukaryotic algae by Yoon’s and one cyanobacteria and 66 eukaryotic algae by Kang’s) were identified, respectively (Table 2). In contrast, the highest heterotrophic bacterial OTUs were identified by Sherwood’s primer set (eight), followed by Yoon’s primer set (six OTUs). Only three heterotrophic bacterial OTUs were identified by Kang’s primer set (Table 2). Proportions of bacterial reads were also highest in Sherwood’s primer set (70.2%) followed by Yoon’s primer (59.34%). Only 0.86% of the bacterial sequences was identified by Kang’s 23S primer set. These results showed that Kang’s 23S universal primer set is a reliable tool for analyzing a phytoplankton community because of its high taxon-specificity, excluding bacterial sequences.

The taxon coverage of three 23S universal primer sets were also compared (Table 2). Five eukaryotic phytoplankton phyla, Bacillariophyta, Chlorophyta, Haptophyta, Miozoa, and Rhodophyta, were identified by Sherwood’s primer set, while two and three additional eukaryotic phytoplankton phyla were identified by Kang’s (Ochrophyta and Streptophyta) and Yoon’s primer sets (Ochrophyta, Streptophyta, and Cryptophyta), respectively (Table 2). Cyanobacterial sequences were identified only by Yoon’s and Kang’s primer sets. Among the 20 most abundant OTUs, 70.21% (7 OTUs) and 52.36% (3 OTUs) were occupied by the bacterial OTUs of Sherwood’s and Yoon’s primer sets, respectively (Table S3), which was not suitable for phytoplankton community analysis presenting dominant bacterial OTUs. In contrast, only one bacterial OTU ranked at 20th with a negligible proportion (0.76%) by Kang’s 23S universal primer set supporting that Kang’s 23S universal primer set is a reliable tool for analyzing the phytoplankton community structure from the entire microbial community with a high taxon specificity and coverage.

Changes in total microbial communities during the bloom

To determine the microbial community changes caused by algal bloom, an NGS analysis was conducted using the 16S universal primer set (Table 1). After trimming and clustering the raw reads, 6,588 reads from the control station, 21,190 from the edge, and 32,461 from the bloom were generated using the 16S universal primer set (Table 3). Clustered with a 99% sequence identity, 161 microbial OTUs were identified from three water samples in the coastal waters of Geoje in 2015. The highest OTU numbers were identified at the edge (103), followed the bloom (89), and control (61) sites (Table 4). All OTUs showed more than a 98% sequence identity indicating the high quality of the 16S database (Table 3). Of 161 OTUs, 23 OTUs were identified as ‘Uncultured Bacteria’ their phyla were determined by phylogenetic analysis. (Fig. 2A).

Thirteen ‘uncultured’ OTUs were bacteria, including eight Bacteroidetes, three Proteobacteria, one Actinobacteria, and one Verrucomicrobia, while seven belonged to Archaea and the final one was the eukaryotic algal species, Bacillariophyta. Finally, it was difficult to determine the taxonomic rank of two OTUs (OTU109 and OTU137) by phylogenetic analysis and were therefore named ‘unclassified OTUs’ (Table 4). Besides these two ‘unclassified’ OTUs, the 159 obtained OTUs amplified by the 16S universal primer set were further classified into 11 phyla, including five prokaryotic heterotrophs (Actinobacteria, Bacteroidetes, Proteobacteria, Verrucomicrobia, and Archaea), 1 eukaryotic protist (Foraminifera), one photosynthetic prokaryote (Cyanobacteria), and 4 photosynthetic eukaryote (Mioza, Ochrophyta, Bacillariophyta, and Cryptophyta) (Table 4 and Fig. 2A).

Community structures of three sample sites generated by the 16S universal primer set were compared (Table 4). In all three sites, OTUs in Proteobacteria were predominant in all three sites and occupied more than 50% of the total OTU numbers and proportions (Fig. 3A), whereas only 24 photosynthetic phytoplankton OTUs of all the microorganisms (cyanobacteria and eukaryotic algae) were identified and their proportions were very low (4.36% and 11.86% at the control and edge sites, respectively) (Table 4). At the control site, heterotrophic bacterial OTUs occupied 80.07%, followed by the eukaryotic protists, Foraminifera (15.57%), and photosynthetic phytoplankton (4.36%). Unlike the microbial community structure at the control site, those at both the bloom and edge sites were highly similar (Table 4). In both bloom and edge sites, proportions of heterotrophic bacteria were 87.53% and 90.5%, respectively, which was higher than those in the control site (Table 4). In heterotrophic bacteria, proportions of Bacteroidetes and Verrucomicrobia were among the most significantly increased phyla during the bloom occurred; their proportions increased by 7.35 and 5.98 folds at the edge site and 8.94 and 5.42-folds at the bloom site, respectively (Table 4). In contrast, Actinobacteria abundance decreased during the bloom. Among the phytoplankton phyla, Cryptophyta, Miozoa, and Ochrophyta, were only identified at the bloom sites, while the proportions of Bacillariophyta were low at both edge and bloom sites (Table 4). Proportions of phytoplankton at bloom and edge sites were 11.86% and 9.23%, respectively, which was also higher than at the control site (4.36%). Collectively, proportions of both photosynthetic algae and heterotrophic bacteria increased as the bloom occurred. However, proportions of Foraminifera, the amoeboid protist phylum, considerably decreased from 15.57% at the control site to 0.27% and 0.13% at the edge and bloom sites, respectively (Fig. 3 and Table 4).

Operational taxonomic units in Cryptophyta, Ochrophyta and Miozoa, were only identified at the bloom and edge sites (Table 4). Interestingly, all the OTUs exclusively identified at bloom sites were those responsible for algal blooming. One Miozoa (OTU10) at the bloom and edge sites was Karenia mikimotoi (GenBank Number: AB027236). Dinoflagellate K. mikimotoi is one of the common species responsible for harmful algal bloom (HAB) causing massive fish mortality and human health risks (Anderson et al., 2010; Chen et al., 2011; Gentien et al., 2007). Two Ochrophyta (OTU12, OTU107) were Heterosigma akashiwo, which is also a well-known species responsible for HAB as well as K. mikimotoi (Nagasaki & Yamaguchi, 1997). Operational taxonomic unit 26 (Teleaulax amphioxeia) and OTU88, (Plagioselmis sp.) in Cryptophyta were also known as the Cryptophyta bloom (Seoane et al., 2012; Šupraha et al., 2014).

To determine the changes in the community structure due to algal bloom, we analyzed the commonly found and bloom-specific OTUs (Table 5, Fig. 4A). Among these 161 OTUs obtained by the 16S universal primer set, 18 OTUs were commonly identified, but their proportions in each site were 64.59%, 64.50%, and 62.93%, respectively (Table 5, Fig. 4A). Except for one cyanobacterial OTU (Synechococcus sp., OTU 15), all other commonly identified OTUs were heterotrophic bacterial sequences. Candidatus pelagibacter (GenBank: LN850161) was identified as the most abundant OTU in all three sample sites, making up more than 30% of the populations in all three sample sites (Supplementary 5). The control site had the highest site-specific OTU numbers (41), followed by the edge (30), and bloom sites (16). Only one OTU was shared between the control, edge, and bloom sites, whereas there were 54 OTUs supporting community structures of the edge and bloom sites, which were highly similar each other (Table 4 and Supplementary 5). Among the 41 control-specific OTUs, Virgulinella fragilis and Rhodobacteraceae sp. occupied about 50% of their proportions (Supplementary 4). The other two site-specific OTUs in both bloom and edge sites occupied only small proportions (3.61% and 6.38%, respectively), which suggests that the bloom did not originate from an outbreak of new species, but proportions of preexisting OTUs changed considerably thereby changing the community structure (Table 5).

Operational taxonomic units with more than two-fold changes were analyzed (Table 6). Proportions of 10 OTUs (seven heterotrophic bacteria and three photosynthetic phytoplankton) increased by more than two fold at the edge and bloom sites, respectively, unlike those from the control sites. Moreover, as shown in the community structure, changes in OTUs at the edge and bloom sites were highly similar (Table 6). Increased heterotrophic bacterial OTUs belonged to phyla Bacteroidetes, Verrucomicrobia, and Proteobacteria (Table 6). Although two phytoplankton OTUs, K. mikimotoi and H. akashiwo were identified in the bloom sites generated by the 16S universal primer set, their proportions were so low due to the outnumbered bacterial sequences (Table 6). Eight OTUs were identified as species that were highly decreased by the bloom (Table 7). The decreased OTUs in both bloom and edge sites were also similar as shown in the increased OTUs. Foraminifera, V. fragilis OTUs decreased the most, followed by two bacterial OTUs, Proteobacteria and Actinobacteria (Table 7). Interestingly, one Rhodobacteraceae sp. (GenBank Number: KU382430) increased while the proportions of the other OTU, which showed a 99% identity to KU382430, decreased at bloom sites (Tables 6 and 7).

Changes in phytoplankton community during the bloom

After trimming and clustering, 25,085, 31,236, and 23,341 reads were finally generated by the Kang’s 23S universal primer set at the control, bloom, and edge sites, respectively (Supplementary 4). At 99% sequence identity, 140 OTUs were obtained and no heterotrophic bacterial OTUs were identified suggesting that Kang’s 23S universal primer set is specific for phytoplankton species (Table 4). The quality of the 23S region database was not as good as that of the 16S region in which species names could not be assigned for about 50% of the OTUs generated by Kang’s 23S primer set at 98% sequence identity (Table 3). Therefore, phylum names were assigned for those exhibiting sequence identities between 90% and 98% as in the previous study (Yoon et al., 2016). Of the total OTUs, 2.15% showed a less than 90% sequence identity to the database and were assigned as “Unknown”. The total 140 OTUs were classified into nine phytoplankton phyla including Haptophyta (27.86%), Chlorophyta (18.57%), Bacillariophyta (14.29%), Cyanobacteria (12.14%), Miozoa (11.43%), Ochrophyta (8.57%), Cryptophyta (3.57%), Cercozoa (0.71%), and Rhodophyta (0.71). Operational taxonomic units in four phyla including Haptophyta, Chlorophyta, Cercozoa, and Rhodophyta were identified exclusively in the results generated by Kang’s 23S universal primer set, which reinforces its importance for microbial community study (Table 4). The most abundant phytoplankton OTU was Ostreococcus sp. (32.78%, GenBank number: KF285533), which belongs to Chlorophyta in the control site. Heterosigma akashiwo (GenBank number: EU168191) was the most abundant OTU in both the bloom and edge sites and occupied 38.31% and 35.67% of the total reads, respectively (Supplementary 5).

Community structures of the three sample sites generated by Kang’s 23S universal primer set were compared (Table 4, Fig. 3B). In the control site, OTUs in Chlorophyta occupied 65.87%, followed by those in Bacillariophyta (14.33%) and Haptophyta (8.44%). Considering the low proportions of phytoplankton (4.36%), which included Bacillariophyta (3.70%) and Cyanobacteria (0.66%), the 16S universal primer set was not as efficient as the 23S universal primer set in presenting phytoplankton community structures (Table 4). In both the bloom and edge sites, proportions of Ochrophyta were highest (43.18% for edge and 40.05% for bloom) followed by Chlorophyta (14.74% for edge and 19.34% for bloom). The phytoplankton community structures by Kang’s 23S universal primer set were largely similar between the bloom and edge sites as shown in the microbial community structure created by 16S (Table 4 and Fig. 3).

Among phytoplankton OTUs, the highest difference in Chlorophyta between the control and bloom sites was identified. Chlorophyta occupied 65.87% in control site, whereas 14.74% and 19.34% was shown at the bloom and edge sites, respectively. By contrast, OTUs in Ochrophyta occupied only 0.18% at the control site, whereas its proportion was 43.18% at edge site and 40.05% at bloom site, respectively. The proportions of Cryptophyta and cyanobacteria at the control site were also 4.48-folds and 6.46-folds and 6.27-folds and 7.69-folds higher than those at the edge and bloom sites, respectively (Table 4). Although Miozoa decreased at the edge and bloom sites by 2.83 and 3.34 folds, the proportion of K. mikimotoi (OTU17) increased. Besides these changes, another algal bloom species, Alexandrium affine (OTU65) was identified only in the bloom site. Three ‘unclassified’ OTUs (OTU27, OTU124, OTU132) were also detected at both the bloom and edge sites, but not the control site (Table 4). Collectively, analysis by Kang’s 23S universal primer set was more sensitive to recent changes in phytoplankton communities during the bloom than those by 16S universal primer set (Table 4).

Commonly identified and site-specific OTUs in all three sites were analyzed (Table 5 & Fig. 4B). Among the 140 OTUs obtained by the 23S universal primer set, 23 OTUs were commonly identified in all three sites at 71.19%, 72.81%, and 76.43%, respectively, which was similar to the results obtained by the 16S universal primer set (Table 5). The commonly identified OTUs were eight in Haptophyta, five in Chlorophyta, four in Bacillariophyta, three in Cryptophyta, and one in Cyanobacteria, Cercozoa, and Ochrophyta. Among the 23 commonly identified OTUs, a small chlorophyte, Ostreococcus sp. (GenBank numbers: KF285533) was the most abundant species, but was not identified by the 16S universal primer set (Supplementary 6). From the 16S universal primers, the control site exhibited the highest site-specific phytoplankton OTU numbers (37), followed by the edge (21) and bloom sites (15). As shown in the results of the 16S universal primer set, more than 95% of the reads shared between the edge and bloom sites by the 23S universal primer set (Table 5 and Fig. 4B). Control-specific OTUs occupied about 14% of their proportions and the other two site-specific OTUs occupied only marginal proportions (Table 5).

To determine the changes in the phytoplankton community structure caused by the algal bloom, OTUs with more than two-fold changes were analyzed (Table 6). At the edge and bloom sites, nine OTUs were significantly higher than in the control site (Table 6). With the exception of four OTUs, all highly increased phytoplankton OTUs generated by Kang’s 23S primer set were responsible for the algal bloom in the ocean (Table 6). Heterosigma akashiwo was identified as the most highly increased phytoplankton species in both the bloom and edge sites, whose proportions were more than 200-fold higher than those in the control site. This result indicated that the major bloom species in the Geoje in 2015 was H. akashiwo. Additionally, we identified changes in other species responsible for the bloom including three Teleaulax spp. and K. mikimotoi. Interestingly, two OTUs in the genus Micromonas exhibited a different pattern in which one OTU increased while the other one decreased (Tables 6 and 7). Among the eight OTUs that decreased, chlorophytes including Ostreococcus sp. were the most highly decreased OTUs in edge site, followed by Chaetoceros sp. (Table 7). Interestingly, Dinophysis acuta was the most highly decreased OTU at the bloom site unlike at the control site, and is also known as the species partly responsible for HAB (Table 7).