Characterization of bacterial community associated with phytoplankton bloom in a eutrophic lake in South Norway using 16S rRNA gene amplicon sequence analysis

Interactions between different phytoplankton taxa and heterotrophic bacterial communities within aquatic environments can differentially support growth of various heterotrophic bacterial species. In this study, phytoplankton diversity was studied using traditional microscopic techniques and the bacterial communities associated with phytoplankton bloom were studied using High Throughput Sequencing (HTS) analysis of 16S rRNA gene amplicons from the V1-V3 and V3-V4 hypervariable regions. Samples were collected from Lake Akersvannet, a eutrophic lake in South Norway, during the growth season from June to August 2013. Microscopic examination revealed that the phytoplankton community was mostly represented by Cyanobacteria and the dinoflagellate Ceratium hirundinella. The HTS results revealed that Proteobacteria (Alpha, Beta, and Gamma), Bacteriodetes, Cyanobacteria, Actinobacteria and Verrucomicrobia dominated the bacterial community, with varying relative abundances throughout the sampling season. Species level identification of Cyanobacteria showed a mixed population of Aphanizomenon flos-aquae, Microcystis aeruginosa and Woronichinia naegeliana. A significant proportion of the microbial community was composed of unclassified taxa which might represent locally adapted freshwater bacterial groups. Comparison of cyanobacterial species composition from HTS and microscopy revealed quantitative discrepancies, indicating a need for cross validation of results. To our knowledge, this is the first study that uses HTS methods for studying the bacterial community associated with phytoplankton blooms in a Norwegian lake. The study demonstrates the value of considering results from multiple methods when studying bacterial communities.


Introduction
In recent years, reports of cyanobacterial blooms in freshwater ecosystems attributed to eutrophication and global warming have increased [1][2][3][4][5][6]. Cyanobacteria rapidly form massive water blooms under favorable conditions, causing economic, ecological and health problems [7][8][9]. Many bloom-forming species produce toxins and the types of toxin produced are dependent on the species composition within the bloom [10][11][12]. Cyanobacteria coexist and interact with heterotrophic bacterial groups commonly found in freshwater ecosystems [13,14]. They may produce various dissolved organic compounds and fix atmospheric nitrogen, supplying heterotrophic bacteria with substrates for growth [15][16][17]. Such compounds can differentially support growth of heterotrophic bacteria species [18]. Therefore, it is likely that the composition of cyanobacterial blooms will structure the associated heterotrophic bacterial populations [14].
Studies on phytoplankton cultures demonstrated a larger difference in bacterial community structure between diatom and cyanobacteria as compared to different cyanobacterial cultures. This is dependent on phytoplankton growth phase and the habitat of the associated bacterial community (free-living within bloom vs attached to the phytoplankton cell surface) [18]. In freshwater ecosystems, environmental conditions, including water temperature, residence time and mixing, light intensity and quality, nutrient concentrations and grazing pressure affect the development of cyanobacterial population [19].
Blooms in temperate regions, occur mostly during late summer and early autumn [6, 20] and are dominated by Microcystis, Aphanizomenon, Anabaena together with the invasive species e.g. Cylindrospermopsis raciborskii [21,22]. Freshwater lakes are numerous in the Nordic countries: Norway (65,000 lakes), Sweden (95,700 lakes) and Finland (187,888 lakes) [23]. Almost 40% of Norwegian water bodies are at risk from human influences such as acidification and eutrophication; climate change is likely to exacerbate these problems in the future (www. environment.no, Freshwater 2015). Reports have shown that increase in water temperature in major European lakes, coupled with decreased duration of ice cover, has caused changes in the life cycle of phytoplankton and promoted the invasion of toxic species from warmer regions [3,24].
Culture-dependent methods have shortcomings, as up to 95% of the microorganisms may be uncultivable [25][26][27]. High throughput sequencing, in contrast to culture dependent studies, can detect a large majority of microbial taxa present. This helps generate a deeper understanding when comparing bacterial communities [18,28]. In this study we have used HTS to assess the bacterial community associated with phytoplankton during the blooming season in the eutrophic Norwegian lake Akersvannet by targeting two hypervariable regions, V1-V3 and V3-V4, of the 16S rRNA gene. To the best of our knowledge, this is the first study of the bacterial community in a freshwater ecosystem in Norway using HTS 16S rRNA amplicon sequencing.

Study site and sampling
Lake Akersvannet is situated in Stokke Municipality, Vestfold County (South Norway). The lake has a mean depth of 6 m, maximum depth of 13 m and a total area of 2.4 km 2 (Fig 1). It is an important wetland habitat for birds during the spring and autumn migration. During summer and fall cyanobacterial blooms occur frequently [29]. The lake is ice-covered during the winter. The lake is surrounded by agricultural land and is used for recreational activities including water sports and fishing. The first documented bloom at this lake was identified as Aphanizomenon flos-aquae (NIVA 1959). Lowering of the water level due to canalization during the year 1968 led to a shift in species composition. By the year 1994, the phytoplankton community was dominated by species of Anabaena and Aphanizomenon and the dinoflagellate Ceratium hirundinella [30].
This study was conducted in close cooperation with Vestfold County administration and the group of land owners surrounding the lake. In accordance with "Section IV" of the "Preservation regulations for Lake Akersvannet nature reserve" (www.lovdata.no) and "Section 3.3" of the "Management plan for Lake Akersvannet" (www.fylkesmannen.no), specific permission for taking water samples for research and educational purposes is not needed. Water samples were collected monthly from June to August 2013. Samples were taken at the sampling point indicated in Fig 1 (59.25˚N 10.33˚E). Separate one liter water samples were collected at depths of 0 m, 2 m and 4 m respectively. Net samples (25 μm pore size) for microscopic analysis were collected at the same time and location. Samples were immediately transferred to a cooled chamber and stored at 4˚C.

Toxin analyses
Enzyme-Linked Immunosorbent Assay technique (ELISA) with Microcystin ADDA ELISA Kit (product No. 520011) from Abraxis (Abraxis LLC, Warminster PA) was used for the microcystin analyses. All samples were freeze-thawed twice to lyse the cells and release the microcystins into the water. An AccuReader (Metertech, Taiwan) model 965 was used to read the absorbance at 450 nm.

Qualitative and quantitative analyses of phytoplankton
Net samples were fixed with Lugol solution and examined under 100 and 400 x magnification for qualitative determination of the phytoplankton community (S1 Table). Phytoplankton enumeration was performed by Utermöhl inverted microscope technique using a 10ml sedimentation chamber. Species identification was performed according to Tikkanen and Willén [31] (S2 Table).

DNA extraction
For each sampling date, three 1L samples, collected at depths of 0 m, 2 m and 4 m were available. Each sample was filtered through four to five 1.2 μm pore size, 47 mm diameter GF/C filters (Whatman, Norway). Each set of filters was cut into fragments in mm-sized pieces with a disposable sterile scalpel blade and incubated in XS buffer for 120 min and subjected to DNA extraction using the Xanthogenate-SDS (XS) method [32]. DNA pellets were washed twice with 70% ethanol, air-dried and suspended in 50μl TE buffer (10 mM Tris-HCL, 1 mM EDTA pH 8). The DNA extract from each set of filters were pooled to prepare one DNA pool per sampling depth. DNA pools from each sampling depth (0, 2 and 4 m) were then combined, giving an integrated sample. This integrated sample was then split into two equal volumes (replicates), and each replicate was used as a template for amplification and sequencing of V1-V3 and V3-V4 target regions of the 16S rRNA gene.

Data analysis
The sequence reads were analyzed using UPARSE v8.0.1623 [36], which allows accurate identification of Operational Taxonomic Units (OTU's) (http://drive5.com/usearch/manual/ uparse_pipeline.html). The raw paired-end sequence reads were merged using the command usearch -fastq_mergepairs with maximum expected error threshold set at 1.0. The merged reads were filtered on the basis of sequence quality using fastq_filter program, setting maximum expected error threshold to 0.5 and minimum sequence length to 250 bp. After quality filtering, the reads were aligned against the 16S rRNA sequences of the RDP database using the Aligner tool [37] in the RDP pipeline available at Ribosomal Database Project (RDP Release 11, Update 4), and endpoints were identified in order to allow their sorting into V1-V3 and V3-V4 reads. For both classes of reads, dereplication was performed to collapse identical reads into one single sequence. Following dereplication, abundance sorting was done and singletons were discarded. OTU's were clustered at 3% divergence threshold and chimeras were eliminated using UCHIME [38]. Representative sequences for each OTU were selected and RDP classifier (16S rRNA training set 14) (S1 and S2 Texts) [39] was used to assign taxonomy. OTU's assigned to Archaea and Chloroplast were filtered out and OTU tables were generated. Random subsampling of the OTU tables for both target regions (S3 and S4 Texts) was performed using script single_rarefaction.py in QIIME (Quantitative Insights Into Microbial Ecology package 1.9.1)[40]. Shared OTU's between replicates were computed using script shared_phylotypes.py (S5 and S6 Texts). Diversity indices (Chao1, ACE and Shannon diversity) were calculated for the subsampled OTU tables using script alpha_diversity.py. The OTU's were further classified using BLAST analysis [41] using nucleotide database (nt/nr) with uncultured/environmental sample sequences excluded. OTU's were only assigned to a species, if the sequence similarity was ! 98% and no other species showed the same level of similarity. Where this was not achieved, the RDP taxonomic assignment was retained.

Accession numbers
The raw 16S rRNA dataset from Lake Akersvannet sample are available at Sequence Read Archive (SRA) under study SRP095055.

Results
The environmental conditions and phytoplankton biomass in Lake Akersvannet during the sampling season 2013 are shown in Table 1. TN ranged from 1.5 to 1.8 mg/l and TP varied between 90 and 148 μg/l. The highest water temperature (20.7˚C) was recorded in August. The highest concentration of microcystin (0.5 μg/l) was detected in August. In June and July, the concentration of microcystin was below the detection limit of 0.15 μg/l.

Phytoplankton community composition: Microscopic identification
Total phytoplankton biomass varied between 0.4 and 16 mg/l ( Table 1) with mean biomass of 5.9 mg/l. Microscopic examination revealed variation in species composition of the phytoplankton community throughout the sampling season. Cyanobacteria represented 77% of total phytoplankton biomass in July, whereas the dinoflagellate Ceratium hirundinella dominated the community in August (96%). Species level morphological identification of Cyanobacteria showed the presence of Aphanizomenon flos-aquae (June-August), Microcystis aeruginosa (June-August), Chroococcous turgidus (June) and Woronichinia naegeliana (August). Aphanizomenon dominated the cyanobacterial community throughout the season. The proportion of Microcystis was increased in August (S2 Table).
Further analysis using BLAST allowed identification of 14 V1-V3 OTU's at species level and three at genus level; similarly, 66 and nine V3-V4 OTU's could be identified at species and genus level respectively (S3 Table). Six V1-V3 and 30 V3-V4 OTU's, were highly similar (! 99%) to unclassified bacterial strain reported in other studies (S4 Table). Taxonomic classification using RDP classifier and BLAST were largely in agreement. There were two discrepancies at phylum level.
Bacterial richness and diversity indices (Alpha diversity) are shown in Table 2. Diversity was highest in July. Diversity was higher for the V3-V4 region than for V1-V3. Replicates gave comparable diversity estimates for each month, except August; August A had fewer OTU's and lower diversity estimates. Comparison of phylotypes revealed that 76% of V1-V3 and 82% of V3-V4 OTU's detected in August A were shared with August B (S5 and S6 Texts).
After completion of taxonomic analysis, 255 V1-V3 OTU's were assigned to 11 phyla and one group of "Unclassified Bacteria", and 502 V3-V4 OTU's were assigned to 17 phyla and one group of "Unclassified Bacteria". A comparison of taxonomic distribution of reads between the two 16S rRNA target region libraries is shown in Tables 3 and 4. This comparison showed that 10 phyla were shared between both target regions, whereas respectively 1 and 7 phyla were exclusively identified in V1-V3 and V3-V4 sequence library. Phylum level distribution of reads for both the target regions is shown in Fig 2. The proportion of reads belonging to "Unclassified Bacteria" and Proteobacteria-alpha is higher in V1-V3 than V3-V4, whereas the proportion of Verrucomicrobia and Bacteriodetes is higher in V3-V4 sequence library.
The phylum and genus level bacterial community distribution is shown in Figs 3 and 4 respectively. The bacterial community was dominated by 5 phyla; namely Actinobacteria, Proteobacteria, Cyanobacteria, Verrucomicrobia and Bacteriodetes along with "Unclassified Bacteria" in varying proportions. Cyanobacteria were detected throughout the sampling period, but dominated both sequence libraries in August. Genera identified were Aphanizomenon (GpI), Microcystis (GpXI), Synechococcus (GpIIa), Pseudanabaena (GpVI), Planktothrix and Woronichinia. Species level distribution of total cyanobacterial reads is shown in Fig 5. The same dominant cyanobacterial genera and species were identified from both the 16S rRNA target regions (data in S3 Table).

Comparison of microscopic and sequence-based cyanobacterial species distribution
The cyanobacterial species distribution derived from microscopic and sequence-based analysis is compared in Table 5. The same major species were found with both methods, but relative species abundance estimates differed.

Discussion
The mean phytoplankton biomass over the study period was 5.9 mg/l, making this lake hypereutrophic according to Brettum and Andersen [42]. Many factors regulate the composition of phytoplankton communities, but TN, TP and the ratio between TN and TP are considered particularly important [43,44]. These parameters fluctuated moderately around means of TN 1.6 mg/l (SD ± 0.1 mg/l) and TP 116.3 μg/l (SD ± 23 μg/l) in this study. Although the low TN: TP ratios observed (10 to 18), were favorable for cyanobacterial blooms, the blooming organism was instead C. hirundinella, which also favors low TN: TP (i.e. < 7) as has been reported earlier in the literature [45][46][47][48]. It plausible that once a Ceratium bloom was established, it suppressed cyanobacterial bloom by competition for light and nutrients or by grazing [48,49]. The increased proportions of Cyanobacteria observed in the 16S sequence library in August, could be due to the increase in water temperature (20.7˚C) ( Table 1). This temperature is optimal for growth of Aph. flos-aquae but suboptimal for M. aeruginosa. The lowest temperature at which Aph. flos-aquae can grow is reported to be 8˚C [50]; optimum being 15 to 28˚C [51], whereas M. aeruginosa can grow slowly at 15˚C [52]; optimum being 27.5 to 32˚C [51,53]. Aph. flos-aquae has higher photosynthetic ability at lower temperatures as compared to M. aeruginosa [52]. In August, ammonium levels rose markedly and this was accompanied by an increase of M. aeruginosa (which favors ammonium as nitrogen source) [54,55]. This would in turn explain the increase of microcystin in August.

Comparison between the two target regions
Although less discriminatory than full length 16S rRNA sequencing, high throughput sequencing of shorter hypervariable regions (V1 to V9) using universal primers has become common practice in microbial community studies [35,56]. Different hypervariable regions have been targeted in different studies with no one region being preferred over others [57,58]. In our study we have utilized two of the most commonly used primer pairs, targeting the V1-V3 and V3-V4 region of the 16S rRNA gene respectively [59]. Of the ten phyla detected by both primer sets, five (Cyanobacteria, Proteobacteria, Actinobacteria, Acidobacteria and Parcubacteria) showed similar abundance profiles in both target sequence libraries, while the other five (Bacteriodetes, Verrucomicrobia, Armatimonadetes, Firmicutes and Planctomycetes) were markedly more abundant in the V3-V4 library. Eight phyla were detected by only one of the primer sets. This suggests that the primers used in this study selectively detect certain phyla and neither is truly universal [34]. For cyanobacteria, however, the results were highly concordant, suggesting that both primers sets have similar coverage for this phylum. Although the overall proportion of phylum Proteobacteria in both libraries was comparable, Alphaproteobacteria were higher in proportion in the V1-V3 library. This variation might be explained by higher coverage of V1-V3 primers for class Alphaproteobacteria as seen in studies conducted in other environment, but will require further experimentation to confirm these findings for freshwater environment [60]. The proportion of unclassified Bacteria and overall unclassified genera was greater in the V1-V3 library. Thus, use of different sets of universal primers can lead to different apparent community structures. During PCR amplification, target sequences that do not precisely match the universal primers will be either less efficiently amplified or lost completely [61] and use of particular sets of universal primers can lead to under-or over-estimation of taxa [35]. Comparison of community richness and diversity estimates showed greater diversity for V3-V4 than for V1-V3 libraries. Such discrepancies have been reported in other studies [35,61]. Replicates gave comparable diversity estimates, except for samples August A and B. Even after random subsampling was performed, August A had fewer OTU's than August B. Such discrepancies have been reported in other studies and attributed to errors in PCR amplification, library preparation or sequencing steps [61]. In our case as both target regions were equally depleted of OTUs in August A, the discrepancy is most probably due to differences in sequencing efficiency as V1-V3 and V3-V4 were sequenced together but amplified separately.
The bacterial community was dominated by five phyla; Proteobacteria, Bacteriodetes, Cyanobacteria, Actinobacteria and Verrucomicrobia in varying proportions throughout the study period. Genus level taxonomic identifications showed the presence of a great diversity of low abundance taxa. Previous studies refer to these low abundance taxa as "rare biosphere". It has been suggested that members of the rare biosphere may play important roles in the bacterial community, despite their scarcity [18,62]. In addition to the rare biosphere, many unclassified OTU's were detected, suggesting that as-yet undescribed bacterial groups with unknown metabolic capabilities are an important part of the lake microbiota, a matter which warrants further investigation. The large proportion of unclassified bacterial taxa also suggests that taxa peculiar to the Norwegian freshwater environment may be absent from the cultured and sequenced species represented by commonly used 16S rRNA databases.

Bacterial community distribution
Bacteriodetes, predominantly Flavobacteriia and Sphingobacteriia dominated the bacterial community throughout the sampling season. These are mostly chemoorganotrophs, that are well adapted to bloom conditions due to their ability to degrade complex biomolecules [18, 63,64] such as algal-derived material [65] including cyanobacterial toxins [13]. Flavobacteriia and Sphingobacteriia have been found in high numbers during episodes of phytoplankton blooms [66], particularly cyanobacterial blooms in freshwater ecosystems [14,67]. A study conducted in Sweden in 2007 reported multiple populations of Flavobacteriia in a bloom of cyanobacteria, without any seasonal succession and suggested that resource availability (organic matter) might be the driving force behind Flavobacteriia community structuring [64]. The genus Flavobacterium, is known to contain species that can lyse Microcystis cells [68] and degrade . Member of this genus are also able to degrade complex organic molecules [64]. The dominant Flavobacterium taxon in June (OTU_11, V3-V4 library) showed high sequence similarity to a microcystin-degrading strain (Flavobacterium sp. AKB-2008 TE28) isolated in a previous study [13]. Fluviicola, another dominant genus of class Flavobacteriia detected in June and August is reported to contain species responsible for denitrification [69].
Actinobacteria have been found in various aquatic habitats and are considered to be ubiquitous members of the freshwater ecosystem, especially in the lake epilimnia [66,70,71]. In this study, the phylum was represented by Actinomycetales and Acidimicrobiales, two of the most abundant freshwater Actinobacterial orders [66,72]. The ratio of allochthonous to autochthonous carbon sources, which may vary greatly during period of phytoplankton bloom, can be selective for different Actinobacterial species [73]. Actinobacteria have been previously associated with blooms of diatoms and cyanobacteria [14,18,71,74,75]. In this study, their abundance declined markedly from June to July while the cyanobacterial abundance increased. This might be due to water temperature increasing above the optimum for Actinobacteria [75] and/or competition from other, faster growing, bacterial phyla [76]. Ilumatobacter, one of the dominant Actinobacteria genera in June (V3-V4 library), are known to thrive in environments with elevated phosphorous levels [77], as found in L. Akersvannet.
In this study, Verrucomicrobia comprised 8.6% of the bacteria detected in the V3-V4 library. Although most studies report low abundance of Verrucomicrobia, much higher abundance was reported in a recent culture independent study [66,78]. Verrucomicrobia abundance may have previously been underestimated due to poor cultivability and methodological biases [66,[78][79][80]. They are known degraders of algal polysaccharides and organic matter [81][82][83]. Hence their presence during a phytoplankton bloom is not unexpected [84]. A representative of Subdivision3 (OTU_8, V3-V4 library), which dominated in July, showed high sequence similarity to a strain known to degrade laminarin, one of the most abundant algal polysaccharides in nature [83].
Most of the OTU's belonging to Actinobacteria and Verrucomicrobia, detected in this study were highly similar to isolates previously identified in North American lakes [83,85], which Phytoplankton bloom associated microbial communities suggests that these strains are geographically widespread. Comparable reports, suggesting that unique capabilities of these groups allow them to compete successfully in lakes with diverse trophic states have been published [70,86]. Proteobacteria represented a significant portion of the bacterial community, which is in agreement with other studies [74,87,88]. Many of the most prevalent and well-studied freshwater bacterial groups belong to Betaproteobacteria [66,89]. Paucibacter (Betaproteobacteria) is represented by a single dominant OTU (OTU_6) which showed high sequence similarity to Paucibacter toxinivorans, a novel species isolated from lake sediment. P. toxinivorans can degrade the cyclic cyanobacterial hepatotoxins and mineralize organic phosphorous compounds [90]. However, the abundance of P. toxinivorans was higher in July (with toxin level below detection threshold) than in August (high toxin concentration detected), suggesting that the toxin presence is not critical for the growth of this species. Limnohabitans, one the betaproteobacterial genus identified in July (V1-V3 library), are opportunistic freshwater bacteria that prefer non-acidic environments [91][92][93]. Their ability to use low molecular weight algal exudates as sole organic carbon source (DOC) and rapid growth rate gives them a competitive advantage during phytoplankton bloom [76,91,94].
Alphaproteobacteria were the dominant proteobacterial group in this study. In contrast to Betaproteobacteria, Alpha and Gammaproteobacteria rarely dominate freshwater bacterial communities, except under conditions of elevated organic and inorganic inputs during phytoplankton blooms [13, 66,76].
One of the Gammaproteobacteria OTU's found in this study was a Rheinheimera species. Rheinheimera are abundant in marine, freshwater and estuarine environments and are able to grow on and degrade organic matter [95][96][97]. Other Gammaproteobacteria were Pseudomonas and Xanthomonadaceae. These groups include opportunistic pathogens of humans, fish and plants [98,99]. Pseudomonas and Rheinheimera are thought to enhance the growth of Microcystis by modulating phosphate exchange in the cyanobacterial mucilage capsule [13,100].

Comparison of microscopic and sequence based results
The microscopic findings and 16S rRNA sequence analysis were in agreement in regard to the major species present. The quantitative results however were not in concordance. Similar findings have been reported previously [112,113]. This discrepancy needs to be addressed, as it is commonly assumed that read counts are semi-quantitative and will provide meaningful relative abundance estimates [114]. A recent study has shown that Cyanobacteria were overrepresented in HTS data as compared to microscopic biovolume data. Thus, for the moment both methods are needed to obtain a picture of the cyanobacterial community structure [115].
The discrepancies observed indicate that there is significant bias in one or both of the methods used. Microscopy is a subjective method, dependent on the skill and experience of the microscopist and overestimation or underestimation of phytoplankton diversity may occur [112,116,117]. In addition, differential losses of cells may occur during fixation, depending on the species and cell size [113,118,119]. Smaller cells that do not sink fast enough within the counting chamber can be underrepresented or missed completely [115]. 16S amplicon read counts are subject to biases caused by cell lysis efficiency, DNA extraction efficiency, PCR amplification efficiency, varying 16S rRNA gene copy number, quality filtering criteria and the sequencing method used [14,113,114,[120][121][122][123]. Our results illustrate the importance of using multiple methods while characterizing freshwater microbiota. Additional experimentation is required to resolve the observed discrepancies.

Conclusion
To our knowledge, this study is the first use of 16S rRNA targeted amplicon sequencing to characterize the bacterial communities associated with phytoplankton blooms in a eutrophic lake in Norway. Our study revealed that choice of 16S rRNA region can affect the apparent community structure and this needs to be taken into account in the planning and interpretation of such studies. Our results showed a highly complex and fluctuating microbial flora, mostly comprised of taxa previously found in association with phytoplankton blooms. A significant proportion of the community was represented by low abundance bacterial groups and unclassified taxa that might be performing hitherto unidentified ecological functions. Our study also showed that the resolving power and the quantitative results of microscopic and 16S rRNA targeted sequence analysis may differ, which needs further investigation.
The 16S rRNA results, provided a snapshot of the phytoplankton-associated bacterial community, at three individual time points. As phytoplankton communities are known to be highly dynamic, more frequent sampling will be needed to identify trends.
Supporting information S1