Revealing Cryptic Changes of Cyanobacterial Community Structure in Two Eutrophic Lakes Using eDNA Sequencing

Harmful cyanobacterial blooms pose a risk to human health worldwide. To enhance understanding on the bloom-forming mechanism, the spatiotemporal changes in cyanobacterial diversity and composition in two eutrophic lakes (Erhai Lake and Lushui Reservoir) of China were investigated from 2010 to 2011 by high-throughput sequencing of environmental DNA. For each sample, 118 to 260 cpcBA-IGS operational taxonomic units (OTUs) were obtained. Fifty-two abundant OTUs were identified, which made up 95.2% of the total sequences and were clustered into nine cyanobacterial groups. Although the cyanobacterial communities of both lakes were mainly dominated by Microcystis, Erhai Lake had a higher cyanobacterial diversity. The abundance of mixed Nostocales species was lower than that of Microcystis, whereas Phormidium and Synechococcus were opportunistically dominant. The correlation between the occurrence frequency and relative abundance of OTUs was poorly fitted by the Sloan neutral model. Deterministic processes such as phosphorus availability were shown to have significant effects on the cyanobacterial community structure in Erhai Lake. In summary, the Microcystis-dominated cyanobacterial community was mainly affected by the deterministic process. Opportunistically dominant species have the potential to replace Microcystis and form blooms in eutrophic lakes, indicating the necessity to monitor these species for drinking water safety.


Introduction
Cyanobacterial blooms in eutrophic freshwater bodies are a globally severe environmental problem [1,2]. The massive proliferation of cyanobacteria has increased the difficulty and cost of water management of urban water supplies [3]. Off-flavor compounds are released from both growing and decaying blooms, resulting in problems with undesirable odors in ambient air and drinking water [4,5]. Moreover, the oxidation of large amounts of organic matter released from cyanobacterial cells may cause oxygen depletion, which is lethal to many aquatic animals such as fish [6]. Harmful cyanobacteria can also produce hepatotoxins, cytotoxins, neurotoxins, and dermatotoxins [7]. Some cyanotoxins are resistant to degradation and persist for a long time in natural environments [8]. As a result, cyanobacterial blooms and their metabolites pose a great risk to the health of aquatic ecosystems and drinking water safety. Therefore, it is important to monitor the dynamics of cyanobacterial communities in water bodies in order to manage cyanobacterial blooms in good time.
Microcystis blooms. Water samples from Erhai Lake were collected in 2010, as described previously [31]. The surface (0-0.5 m) water from three sites (N2, M1, and S2) in four months throughout the year (February, May, August, and December) were used in this study. Water samples at a depth of 10 m were collected from the N2 and M1 sites at the same time. Water samples from Lushui Reservoir were collected at two sites (W and E) in three months (February, May, and September) of 2011. The surface water from the two sites and water samples at a depth of 15 m from the E site were used. In total, 28 water samples were collected. All the samples were numbered according to sampling sites, depth, and month, and their information is listed in Figure 1 and Table S1. Both Erhai Lake and Lushui Reservoir are located in the subtropical monsoon climate zone, and they have been suffering from eutrophication and annual Microcystis blooms. Water samples from Erhai Lake were collected in 2010, as described previously [31]. The surface (0-0.5 m) water from three sites (N2, M1, and S2) in four months throughout the year (February, May, August, and December) were used in this study. Water samples at a depth of 10 m were collected from the N2 and M1 sites at the same time. Water samples from Lushui Reservoir were collected at two sites (W and E) in three months (February, May, and September) of 2011. The surface water from the two sites and water samples at a depth of 15 m from the E site were used. In total, 28 water samples were collected. All the samples were numbered according to sampling sites, depth, and month, and their information is listed in Figure 1 and Table S1. For each sample, cyanobacteria were collected by filtering a volume of 300 mL water through a polycarbonate membrane filter (0.45 µm pore size, Millipore, USA). The filters were preserved at −80 °C prior to DNA extraction. Surface water samples were also used to analyze the concentrations of total phosphorus (TP), total nitrogen (TN), ammonium nitrogen (NH4-N), and nitrate nitrogen (NO3-N), as described by Wu et al. (2006) [32]. The water temperature (WT) and dissolved oxygen (DO) were measured using YSI PH100 [31], and the water transparency (SD) was measured using a Secchi disk.

High-Throughput Sequencing of Cyanobacterial cpcBA-IGS
Environmental DNA (eDNA) was extracted from the filters using the method developed in an earlier study [33]. Cyanobacterial cpcBA-IGS sequences were amplified from the eDNA of all the samples using a forward primer, PCβF [34], targeting the cpcB gene, and a reverse degenerate primer, PCαR454 [26] targeting the 5′ end of the cpcA gene. Adapters and specific barcode sequences for each sample were designed and added to the 5′ side of the PCαR454 primer for 454 pyrosequencing (Table S2, in the supplementary data). As described by Jiang et al. (2017) [26], PCR reactions were performed in triplicate using Tks Gflex DNA polymerase (Takara, Japan), and the For each sample, cyanobacteria were collected by filtering a volume of 300 mL water through a polycarbonate membrane filter (0.45 µm pore size, Millipore, USA). The filters were preserved at −80 • C prior to DNA extraction. Surface water samples were also used to analyze the concentrations of total phosphorus (TP), total nitrogen (TN), ammonium nitrogen (NH 4 -N), and nitrate nitrogen (NO 3 -N), as described by Wu et al. (2006) [32]. The water temperature (WT) and dissolved oxygen (DO) were measured using YSI PH100 [31], and the water transparency (SD) was measured using a Secchi disk.

High-Throughput Sequencing of Cyanobacterial cpcBA-IGS
Environmental DNA (eDNA) was extracted from the filters using the method developed in an earlier study [33]. Cyanobacterial cpcBA-IGS sequences were amplified from the eDNA of all the samples using a forward primer, PCβF [34], targeting the cpcB gene, and a reverse degenerate primer, PCαR454 [26] targeting the 5 end of the cpcA gene. Adapters and specific barcode sequences for each sample were designed and added to the 5 side of the PCαR454 primer for 454 pyrosequencing (Table S2, in the supplementary data). As described by Jiang et al. (2017) [26], PCR reactions were performed in triplicate using Tks Gflex DNA polymerase (Takara, Japan), and the amplification products were purified and pooled to minimize random bias. High-throughput sequencing was performed using the GS-FLX Titanium platform (Roche 454 Life Sciences, Branford, CT, USA). High-quality reads were selected by discarding low-quality reads, reads with barcode and/or primer errors, sequences shorter than 200 bp, and chimeras. Afterwards, robust sequences were obtained by trimming the barcode and primer sequences from the 5 end. The sequence data were deposited under the accession number PRJNA638599 in the Sequence Read Archive of National Center for Biotechnology Information (NCBI).

Sequence Clustering
Unique sequences were picked out from the dataset of robust sequences, followed by aligning against an alignment template of reference cyanobacterial sequences [26] using Mothur v.1.33.0 [35]. The aligned cpcBA-IGS sequences were clustered into operational taxonomic units (OTUs) using the average neighbor clustering algorithm with a distance of 0.05 by Usearch v11.0.667 [36]. Prior to further analysis, the sequences of all the samples were subsampled and normalized to the smallest read number (i.e., 4560 in EHM11002). Subsequently, the OTU richness, Coverage, Chao 1 index, Shannon-Wiener index (here called Shannon index), Simpson index, Simpson's evenness index, and rarefaction values were estimated using the vegan package of R v3.6.2 (https://www.r-project.org).

Taxonomy Assignment of OTUs
Representative sequences of OTUs were aligned with cpcBA-IGS sequences from the reference cyanobacterial strains used by Jiang et al. (2017) [26]. The alignment was conducted using the ClustalW v2.0 program (Conway institute, University College Dublin, Ireland). A phylogenetic tree of cpcBA-IGS sequences was created based on the neighbor-joining algorithm using the MEGA v7.0 software [37]. Nucleotide substitution was fitted by the Kimura 2-parameter model, and the bootstrap replications were set as 1000. The OTUs were annotated by the closest related cyanobacterial species on the phylogenetic tree. The taxonomy assignment of each OTU was also verified by homologous searching against the GenBank database on the NCBI website.

Neutral Model for Cyanobacterial Community
The impacts of stochastic and deterministic processes on the assembly of cyanobacterial communities were evaluated by predicting the relationship between the occurrence frequency and relative abundance of OTUs using the Sloan neutral model [38,39]. Subsequently, the OTU data of Erhai Lake and Lushui Reservoir were fitted to the neutral model. Two parameters, the goodness of fit (R 2 ) and immigration rate (m), were calculated. The 95% confidence intervals around the fitting statistics were estimated by bootstrapping, and the bootstrap replicates were set as 1000. All model calculations were performed using R v3.6.2.

Statistical Analysis
The OTU abundance was the square root transformed and used for calculating the Bray-Curtis and Jaccard dissimilarity among samples in each water body. The independent effects of different environmental variables on the beta diversity of samples were estimated by hierarchical partitioning [40] and statistical significance was calculated by a randomization test. The correlations between environmental factors and cyanobacterial community structure were also evaluated by a Mantel test [41] using Pearson's coefficients. The environmental factors were ln (x + 1) transformed and Euclidean distances between samples were calculated. The effects of the ratio of TN and TP (TN/TP) on cyanobacterial community structure was also evaluated by a Mantel test. All analyses were performed using the hier.part and vegan packages of R v3.6.2.

Sequence Data of cpcBA-IGS
Twenty-eight filter samples were collected from Erhai Lake and Lushui Reservoir, respectively (Table S1). High-quality eDNA was extracted from the filters, and cpcBA-IGS was amplified and sequenced. A total of 210,806 robust cpcBA-IGS reads were obtained after sequence processing. Each sample comprised 4560 to 11,784 sequences. The length of all the cpcBA-IGS reads varied from 267 bp to 418 bp ( Figure S1 in the supplementary data), indicating a high variability of this genomic region.

OTU Composition and Diversity
Following the normalization of the sequence number for each sample, a total of 3339 OTUs were obtained with 118 to 260 OTUs per sample (Table 1). Overall, the samples from Erhai Lake contained more OTUs than those from Lushui Reservoir. For all the samples, a relatively low coverage of OTU richness was found, with values between 0.94 and 0.97. The Chao 1 index values were much higher than the OTU richness values, suggesting the existence of many rare OTUs. This result was consistent with the extremely low evenness (≤0.03) of the OTU composition. The rarefaction curves of the OTU richness and Chao 1 index did not plateau at the sequencing depth in this study ( Figure S2A,B), indicating that there were undetected OTUs, as suggested by the low coverage of OTU richness. The Shannon index (H ) and Simpson's reciprocal index (D −1 ) were calculated based on the OTU composition of each sample (Table 1). These two diversity indices varied similarly among the samples. In Erhai Lake, although the highest OTU diversity was found in one December sample, EHN2D1012 (OTU = 260, H = 3.00, D −1 = 8.90), other samples collected in December always had a lower OTU diversity despite being from different sampling sites. For both water bodies, the OTU diversity of each site was relatively high in February, with the exception of the S2 site in Erhai Lake. The rarefaction curves of the Shannon index were in close proximity to the plateau ( Figure S2C, in the supplementary data), indicating that most of the abundant OTUs had been detected, as the Shannon index is mainly affected by the percentage of abundant OTUs.

Cyanobacterial Community Structure
The normalized abundance of each OTU was calculated. Subsequently, the OTUs with an abundance accounting for more than 1% of the 4560 in the whole dataset were selected and defined as abundant OTUs. Fifty-two abundant OTUs accounted for 95.2% of the total sequences (Additional file 1). Therefore, these OTUs represented the majority of the cyanobacterial community. Through a phylogenetic analysis, the 52 OTUs were classified into nine groups of cyanobacteria, including Synechococcus group I (11 OTUs), Synechococcus group II (9 OTUs), Synechococcus-related species (6 OTUs), Synechocystis (2 OTUs), Microcystis (12 OTUs), Raphidiopsis (1 OTU), mixed Nostocales species (7 OTUs), Planktothrix (2 OTUs), and Phormidium (2 OTUs) ( Figure S3, in the supplementary data). The mixed Nostocales species were composed of Dolichospermum, Anabaena, Aphanizomenon, and Cuspidothrix, which could not be discriminated in the phylogenetic tree of this study.
As displayed in Figure 2, Microcystis was the most abundant cyanobacterial genus in both water bodies. The relative abundance of Microcystis varied between 29 and 99% in Erhai Lake and between 52 and 99% in Lushui Reservoir. In Erhai Lake, mixed Nostocales species were also dominant and their relative abundance even exceeded Microcystis in two samples, EHN2D1002 and EHS21008, with 36 and 47% mixed Nostocales species, respectively. Phormidium was an opportunistically dominant species in Erhai Lake, with an abundance of 50% in the EHS21005 sample. The predominance of Synechococcus group II only occurred in the EHN2D1012 sample with a relative abundance of 27%. However, in Lushui Reservoir, a high relative abundance of Synechococcus group II existed in two samples, 22 and 33% in LSW1102 and LSED1102, respectively. Synechococcus group I was dominant in LSW1102, with a relative abundance of 22%.
The Jaccard dissimilarities among samples were calculated and used to construct cluster dendrograms on the basis of shared abundant OTUs. According to the presence of OTUs, the samples from different water bodies tended to form independent clusters, with the exception of the EHN2D1005 and EHN2D1008 samples, which were closely related to the samples from Lushui Reservoir ( Figure 3A). However, several samples from Erhai Lake were clustered together with the samples from Lushui Reservoir based on the abundance of OTUs ( Figure 3B). The divergence of cyanobacterial communities at different depths was site-and month-dependent. The surface and deep water samples were obviously different for the N2 site in Erhai Lake in May, August, and December. In each water body, seasonal differences were not obvious for both the membership and structure of the cyanobacterial community.

Application of the Neutral Model for Cyanobacterial Community
The Sloan neutral model describes the stochastic assembly of the microbial community. For the relationship between the OTU occurrence frequency and relative abundance in Erhai Lake and Lushui Reservoir, the values of R 2 and m of the neutral model were low (Figure 4). These results indicated that the stochastic process was not the main factor driving the assembly of cyanobacterial communities in these two water bodies. Moreover, the abundant OTUs, with greater abundances than that predicted by the neutral model, belonged to Microcystis in both water bodies.  The Jaccard dissimilarities among samples were calculated and used to construct cluster dendrograms on the basis of shared abundant OTUs. According to the presence of OTUs, the samples from different water bodies tended to form independent clusters, with the exception of the EHN2D1005 and EHN2D1008 samples, which were closely related to the samples from Lushui Reservoir ( Figure 3A). However, several samples from Erhai Lake were clustered together with the samples from Lushui Reservoir based on the abundance of OTUs ( Figure 3B). The divergence of cyanobacterial communities at different depths was site-and month-dependent. The surface and deep water samples were obviously different for the N2 site in Erhai Lake in May, August, and December. In each water body, seasonal differences were not obvious for both the membership and structure of the cyanobacterial community.

Application of the Neutral Model for Cyanobacterial Community
The Sloan neutral model describes the stochastic assembly of the microbial community. For the relationship between the OTU occurrence frequency and relative abundance in Erhai Lake and The Sloan neutral model describes the stochastic assembly of the microbial community. For the relationship between the OTU occurrence frequency and relative abundance in Erhai Lake and Lushui Reservoir, the values of R 2 and m of the neutral model were low (Figure 4). These results indicated that the stochastic process was not the main factor driving the assembly of cyanobacterial communities in these two water bodies. Moreover, the abundant OTUs, with greater abundances than that predicted by the neutral model, belonged to Microcystis in both water bodies.

Correlations Between Environmental Factors and Cyanobacterial Community
As shown in Figure 5, the TP had significant (p < 0.05, randomization test) independent effects for both Bray-Curtis and Jaccard dissimilarities of the cyanobacterial community in Erhai Lake. In addition, the effects of TP were larger than those of TN, NH 4 -N, NO 3 -N, WT, DO, and SD. For the cyanobacterial community in Lushui Reservoir, NH 4 -N had significant (p < 0.05, randomization test) independent effects for their dissimilarities, and the effects were larger than those of TP, TN, NO 3 -N, WT, and DO.

Correlations Between Environmental Factors and Cyanobacterial Community
As shown in Figure 5, the TP had significant (p < 0.05, randomization test) independent effects for both Bray-Curtis and Jaccard dissimilarities of the cyanobacterial community in Erhai Lake. In addition, the effects of TP were larger than those of TN, NH4-N, NO3-N, WT, DO, and SD. For the cyanobacterial community in Lushui Reservoir, NH4-N had significant (p < 0.05, randomization test) independent effects for their dissimilarities, and the effects were larger than those of TP, TN, NO3-N, WT, and DO. The results of the Mantel test revealed that TN/TP had significantly (p < 0.01, Mantel test) effects on the cyanobacterial community composition in Erhai Lake (Table 2), whereas the effects of TP was relatively weak (0.05 < p < 0.10, Mantel test). The cyanobacterial community composition in Lushui Reservoir was correlated (p ≤ 0.05, Mantel test) with the WT and had a weak correlation with the DO (0.05 < p < 0.10, Mantel test).

Erhai Lake
Lushui Reservoir The results of the Mantel test revealed that TN/TP had significantly (p < 0.01, Mantel test) effects on the cyanobacterial community composition in Erhai Lake (Table 2), whereas the effects of TP was relatively weak (0.05 < p < 0.10, Mantel test). The cyanobacterial community composition in Lushui Reservoir was correlated (p ≤ 0.05, Mantel test) with the WT and had a weak correlation with the DO (0.05 < p < 0.10, Mantel test).

Discussion
In this study, the diversity and composition of cyanobacterial communities in two eutrophic freshwater lakes were described using OTUs obtained based on the results of next-generation sequencing. Correlations between OTU occurrence and abundance were tested using the Sloan neutral model. The effects of environmental factors on cyanobacterial community structure were also evaluated.
The results of this study revealed that the two eutrophic lakes have different cyanobacterial community structures, although they both displayed Microcystis blooms. Compared to Lushui Reservoir, a larger cyanobacterial diversity was observed in Erhai Lake. In a previous investigation, an extremely large number of genotypes were detected for Microcystis in Erhai Lake [9]. A negative correlation was determined between the Microcystis diversity and eutrophication levels [9]. These findings were verified in the present study because the eutrophication level was lower in Erhai Lake than in Lushui Reservoir, which had higher nutrient concentrations (Table S3 in the supplementary data). Therefore, it seems that the genetic diversity of both cyanobacterial communities and populations of individual species were affected by the eutrophication levels.
In most cases, Microcystis had a relatively higher abundance than mixed Nostocales species, although the latter also occurred throughout the year. Earlier studies found that predation by zooplankton is a key factor controlling the relative abundance of different cyanobacterial species in aquatic ecosystems [42,43]. Compared to the filaments of Nostocales species, Microcystis colonies are more resistant to grazers [42,43]. In addition, there are interspecific interactions, such as allelopathic effects between Microcystis and Nostocales species. Under nutrient-replete conditions, allelochemicals released from Microcystis can suppress the nitrogen fixation of Dolichospermum and lead to cytotoxicity [21]. In co-culture experiments, the growth of Aphanizomenon was inhibited by some Microcystis strains, and this effect was caused by inducible unknown chemicals secreted by Microcystis [44]. A previous study also found that biological interactions probably contribute to a large proportion of the inter-annual variability of cyanobacterial bloom [23]. Therefore, the cyanobacterial community structure is influenced by both biological interactions and environmental factors [22].
The difference of cyanobacterial community structure between surface and deep water was site-and month-dependent. Microcystis, mixed Nostocales species, Synechococcus, and Phormidium constituted the major portion of the cyanobacterial community from surface to deep water in Erhai Lake, whereas Phormidium was scarce in some samples. However, Microcystis and Synechococcus dominated the cyanobacterial community throughout the water column in Lushui Reservoir. These results are in contrast with previous findings in another deep lake, Dongzhen Reservoir, where distinct cyanobacterial community structures were found between the surface and deep water [26]. Synechococcus was completely dominant in the surface water of Dongzhen Reservoir [26]. The proliferation of Synechococcus was also recorded in other freshwater lakes [27] and lagoons [45]. Due to the tiny size of this cyanobacterium, little information is available about the ecological risk of its blooms. However, Synechococcus blooms may be lethal to sponges and seagrass, and may have a long-term impact on aquatic ecosystems by increasing baseline chlorophyll a concentrations [45].
As an opportunistically dominant species, Phormidium proliferated occasionally in Erhai Lake. However, reports on bloom-forming Phormidium are scarce. In the 1980s, Phormidium blooms, which damaged tap water by releasing off-flavor compounds, occurred annually in a deep lake (Biwa Lake in Japan) [46]. Phormidium was also one of the dominant cyanobacterial genera causing cyanobacterial blooms in the Greater Sudbury Area lakes of Canada [47]. However, the global proliferation of Phormidium is rising in rivers, where it forms benthic mats [48]. In addition, several species of Phormidium produce neurotoxic anatoxins [49]. Therefore, it is also necessary to monitor the abundance of Phormidium and its toxin content in freshwater lakes where this species is known to occur frequently.
Results from the Sloan neutral model suggested that the assembly of cyanobacterial communities in the two eutrophic lakes was affected by both stochastic and deterministic processes, and the latter may play a more important role. According to the results of hierarchical partitioning and the Mantel test, the cyanobacterial community composition of Erhai Lake was correlated with the TP. Previous studies have also suggested the key role of phosphorus in the succession and dominance of different cyanobacterial populations [11,21,22,28]. Thus, phosphorus concentration is likely to be a deterministic factor of the cyanobacterial community structure. For Lushui Reservoir, the results of hierarchical partitioning and the Mantel test were not consistent and significant environmental factors affecting the cyanobacterial community were not determined. Compared to TP (0.02 ± 0.01 mg L −1 ) in Erhai Lake, the TP (0.03 ± 0.01 mg L −1 ) in Lushui Reservoir was higher, indicating that phosphorus may not be a limiting factor for cyanobacteria in this water body. Moreover, in Lushui Reservoir, the OTU diversity was relatively low and the OTU01 of Microcystis had an extremely high relative abundance of 88-95% in most of the samples, indicating a robust cyanobacterial community structure dominated by Microcystis.
According to the TN/TP criteria [50], both Erhai Lake and Lushui Reservoir had high TN/TP values (Table S3) and seemed to be phosphorus-limited (TN/TP > 17). However, a significant correlation between the cyanobacterial community composition and TN/TP was observed in Erhai Lake but not Lushui Reservoir. This result may be ascribed to the higher nutrient concentrations in Lushui Reservoir. A previous study also suggested that TN/TP was an effective environmental factor for cyanobacterial bloom only within a certain nutrient level [31].
In addition, our findings revealed that only abundant OTUs belonging to Microcystis were found to have occurrence frequencies above the neutral model prediction. The dominance of Microcystis across the year was consistent with the strong acclimation of this species to different water temperatures [51,52]. In terms of the life-cycle of Microcystis, early studies proposed that this cyanobacterium overwinters in the sediment [53][54][55][56]. Recently, in a shallow eutrophic lake, a Microcystis bloom was found to persist in winter although at a relatively lower abundance than in summer [52]. The phenotypic plasticity of Microcystis in response to variations in CO 2 concentration could also be beneficial for the adaptation of this species to different environments [57]. Thus, the occurrence and strength of Microcystis blooms will probably be intensified under the conditions of climate warming and elevated CO 2 concentrations [22,57,58].

Conclusions
A similar cyanobacterial community structure may distribute across the water column of eutrophic lakes. The diversity and composition of cyanobacterial communities are affected by the geochemical characteristics of individual water bodies. Deterministic processes, such as nutrient competition, have more important effects on the succession of cyanobacterial communities than stochastic processes. Although Microcystis is predominant in some freshwater ecosystems, opportunistically dominant species, such as Phormidium, may also proliferate and should be considered in the daily monitoring of cyanobacterial blooms.
Supplementary Materials: The following are available online at http://www.mdpi.com/1660-4601/17/17/6356/s1, Figure S1: Length distribution of high-quality sequences, Figure S2: Rarefaction curves of OTU richness (A), Chao 1 index (B), and Shannon index (C) against the number of sequences, Figure S3: Unrooted neighbor-joining tree of representative cpcBA-IGS sequences from 52 OTUs and reference cyanobacterial strains. Genbank accession number of each reference sequence was displayed before the species name. Bootstrap values above 50% are indicated at the nodes of the tree, Table S1: Detailed information of samples and number of robust reads, Table S2: Reverse primers used for amplification and pyrosequencing of cpcBA-IGS, Table S3: Water quality parameters in the two eutrophic lakes, Additional file 1: The normalized number of sequences of dominant OTUs in the samples.