Spatiotemporal changes of bacterial communities during a cyanobacterial bloom in a subtropical water source reservoir ecosystem in China

Cyanobacterial blooms, which not only threaten the health and stability of aquatic ecosystems but also influence the microbial community within, emerges as one of the most concerning problems in China. However, how cyanobacterial blooms affect the spatiotemporal variation of aquatic microbial communities remains relatively unclear. In this study, we used high-throughput sequencing to investigate how the cyanobacterial and bacterial community spatiotemporally vary along with main cyanobacterial bloom phases in upstream rivers of a eutrophicated water source reservoir. Both cyanobacterial and bacterial diversities in each river were significantly lower (P < 0.05) during the bloom outbreak phase, showing the apparent influence of cyanobacterial bloom. Dominant cyanobacterial taxa included Cyanobacteriales and Synechococcales, and dominant bacterial taxa comprised Acinetobacter, CL500-29, hgcI clade, Limnohabitans, Flavobacterium, Rhodoluna, Porphyrobacter, Rhodobacter, Pseudomonas, and Rhizobiales, whose changes of relative abundance along with the bloom indicated distinct community composition. Non-metric multidimensional scaling analysis proved that community composition had significant difference amongst bloom phases. Linear discriminant analysis (LDA) with LDA effect size analysis (LEfSe) identified unique dominant cyanobacterial and bacterial OTUs at different phases in each river, indicating spatiotemporal variations of communities. Canonical correlation analysis or redundancy analysis revealed that at different bloom phases communities of each river had distinct correlation patterns with the environmental parameters (temperature, ammonium, nitrate, and total phosphorus etc.), implying the spatial variations of microbial communities. Overall, these results expand current understanding on the spatiotemporal variations of microbial communities due to cyanobacterial blooms. Microbial interactions during the bloom may shed light on controlling cyanobacterial blooms in the similar aquatic ecosystems.

. Diagram of Yankou Reservoir Basin indicating the sampling points on each river. The sampling sites were located using the global positioning system (GPS) and were nominated after the names of corresponding rivers. Map was created with ArcGIS software (V10.6, https:// deskt op. arcgis. com/ zh-cn/ deskt op/). Field photos were taken at each sampling point (pins only indicated some of the locations for concise presentation), and representative ones were selected to show the bloom course from May to September. Bioinformatic treatment. The sequenced paired-end reads were merged with FLASH (v. 1.2.11), and the raw reads were processed and analyzed using QIIEM (v.1.9.1) to remove low quality reads. After quality control, sequences were clustered into operation taxonomic units (OTUs) using UPARSE (v.7.0.1090) with a 97% sequence similarity threshold 16 . A representative sequence from each OTU was run against the databases of SILVA (Release 138). To avoid mistaken taxonomic alignments, the taxonomic classifications were doublechecked against the reference prokaryotes. Those unclassified OTUs were discarded, followed by normalization of the filtered sequence data to minimize sequencing biases and allow for appropriate comparison of community variations.
For comparisons of microbial taxa variations across the cyanobacterial bloom phases, we defined those taxa with relative abundance less than 0.01% as rare taxa, more than 1% as abundant taxa, and between 0.01% and 1% as moderate taxa, respectively. Based upon the definition, those relevant taxa in subsequent analyses were categorized into groups as follows: always abundant taxa (AAT), always rare taxa (ART), always moderate taxa (AMT), and conditionally varied taxa (CVT), and the taxa with more than 0.01% relative abundance at certain phases in each sample were also referred as "dominant taxa" according to previous studies 17 . Statistical analyses. In this study, the bloom was mainly attributed to cyanobacterial overproduction, therefore, the cyanobacterial community would be further extracted from the normalized data for each sample, and alpha-and beta diversities were estimated both on cyanobacterial and the rest of bacterial communities (will be referred as "bacterial communities" below). Additionally, a comprehensive analysis on the changes of bacterial communities at different phases in each river as well as the variance of communities at the same phase amongst different rivers was performed.
To better compare the difference of communities on relative abundances and community compositions, both alpha-and beta-diversity indices were calculated. The alpha-diversity indices including OTU numbers (Sobs), abundance-based coverage estimator (ACE), Shannon, and Chao 1, as well as Venn diagrams compared to different samples were calculated using the vegan package in R language (v.4.0.0). Rarefaction curves and Good's coverage were performed with MOTHUR (v. 1.30.2). Significance amongst alpha-diversity indices were calculated using one-way ANOVA and Student's t-test with the significant level at 0.05. For beta-diversity analysis, the nonmetric multidimensional scaling analysis (NMDS) with weighted unifrac similarity coefficient and permutational multivariate analysis of variance test (PERMANOVA) was used to evaluate both the discrepancy of bacterial communities at different phases and the variance of communities at the same phase in different rivers.
To identify the cyanobacterial and bacterial OTUs contributing to the difference amongst samples described above, a significance test was first performed using the Kruskal-Wallis test coupled with Tukey-Kramer posthoc examination, based upon OTUs' relative abundances. Subsequently, the linear discriminant analysis (LDA) coupled with the LDA effect size (LEfSe) technique, with a discriminant analysis score of 2.0, was performed to identify the significantly different dominant OTUs amongst the samples. Additionally, cyanobacterial and bacterial OTUs with LDA scores over 3.0 and 4.0, respectively, were screened out to double-check with the significance tests. The significance tests and LEfSe analyses were all examined at a significant level of 0.05.
Finally, a variance inflation factor analysis (VIF) was first performed on all of the environmental parameters, only when VIF scores of the parameters were under 10, would the parameters be selected for subsequent correlation analyses. In this study, the VIF scores of DTN, DTP, and COD were over 10, which would be eliminated. Canonical correlation analysis (CCA) or redundancy analysis (RDA) was conditionally chosen based upon results of detrended correspondence analysis (DCA) to reveal the correlations.

Alpha-diversity comparisons and statistics.
For cyanobacterial communities in pools of HS, JFZ, SH, and XH across the bloom phases, the Shannon index maximized at the "During AB" phase for JFZ and XH. Although it maximized at the "Post AB" and "Before AB" phase for HS and SH, respectively, the diversities at the "During AB" phase were also relatively high for both pools (Fig. S2A). Whilst Chao 1 index of these cyanobacterial communities showed a concurrent pattern that community diversities minimized at "During AB" and maximized at "Post AB" phase. For bacterial communities, both Shannon and Chao 1 indicated that community diversity decreased at "During AB" and increased afterward. On the contrary to pools, cyanobacterial communities in all rivers minimized at "During AB" phases according to Shannon and Chao 1 index (Fig. S2B), except for Shannon of HS and JFZ. Similar to pools, bacterial communities in all rivers were minimized at "During AB" phases and increased afterward based upon both indices estimates.
Differences on alpha-diversity were further determined amongst different rivers at the same phase as well as each river along the bloom phases (Fig. 2). For the cyanobacterial communities at each phase, there was, in general, no significant difference amongst the studied rivers (P > 0.05), except for the comparisons between HS and JFZ at "Before AB", as well as JFZ and SH at "Post AB" phases (P < 0.05) ( Fig. 2A). Whereas, community diversities in individual rivers exhibited apparent variances across the consecutive phases (Fig. 2B). Cyanobacterial communities in HS showed no significant difference across the phases, whilst communities in JFZ, SH, and XH demonstrated significant differences amongst these phases, and especially, communities at the "During AB" phase were significantly lower (P < 0.05) than the rest phases according to Chao 1 index for JFZ, SH, and XH.
Similarly, for bacterial communities at the same phase, there was, in general, no significant difference amongst the rivers (P > 0.05), except for the comparison between JFZ and XH according to Shannon index (Fig. 2C). Bacterial communities in individual rivers also showed similar variance with the cyanobacteria across the phases (Fig. 2D). HS communities had no significant difference amongst the phases (P > 0.05), whilst JFZ, SH, and XH communities showed significant differences, and communities at the "During AB" phase were also significantly lower (P < 0.05) than some of the rest phases in JFZ, SH, and XH, respectively.
OTU comparisons of both cyanobacteria and bacteria in tributary pools and their rivers at corresponding phases were performed (Fig. S3). In general, for both cyanobacterial and bacterial OTUs, the comparisons either amongst pools and rivers at certain phases ( Fig. S3) or in each pool and river at different phases (Fig. S3), only a small portion of OTUs was shared in common amongst the four compared counterparts. It should be noted that only a small number of cyanobacterial OTUs was classified in pools of JFZ, SH, and XH during the bloom (Fig. S3), which supported the alpha-diversity analyses that the cyanobacteria diversities during the bloom were significantly lower than the rest phases for JFZ, SH, and XH (Fig. 2).
Comparisons of cyanobacterial and bacterial community composition. In general, the most diverse cyanobacterial OTUs were assigned to phylogenetic orders of Cyanobacteriales and Synechococcales in both rivers and their tributary pools. However, the cyanobacterial community composition exhibited spatial variance amongst these ecosystems. In each pool, the 2 taxa took turns as the most dominant community, for example, in the tributary pool of HS, the most dominant community was Synechococcales and Cyanobacteriales across each of the phases; whilst JFZ had Cyanobacteriales as the most dominant community. And Synechococcales was the most dominant in the rest of the three rivers along with the bloom (Fig. 3A).
For bacterial communities, the abundant genus to which most diverse OTUs were assigned were sorted out in each river. Based upon the relative abundance of each taxon, the phylogenetic genera were further categorized into groups of AAT, ART, AMT, and CVT, respectively ( Fig. 3B-E). Although these rivers had some genera in common in each group (AAT, ART, AMT, and CVT), such as Acinetobacter and Flavobacterium the bacterial community composition also exhibited spatial variance amongst the studied rivers. For example, a community of hgcI clade was always abundant in HS, JFZ, and SH, whilst was conditionally varied in XH, and its relative abundance was greatly reduced during the cyanobacterial bloom, especially in XH. Communities of Pseudarcicella and Limnohabitans in HS were always abundant, and their relative abundance increased during the cyanobaterial www.nature.com/scientificreports/ bloom, however, they were conditionally varied and noticeably reduced in the rest three rivers during the same phase. Similarly, Comamonadaceae was always abundant in each river, but its relative abundance clearly reduced during the bloom, and Allobaculum was always rare in HS, SH, and XH, but conditionally varied in JFZ that its relative abundance greatly increased during the bloom.
Variations of cyanobacterial and bacterial communities were further illustrated through NMDS with PER-MANOVA test. Cyanobacterial communities at the "Before AB" phase exhibited significant variation amongst the studied rivers (PERMANOVA, df = 3, F model = 3.220, R 2 = 0.305, P = 0.001), however not at the rest phases (Fig. 4A). Across the bloom phases, cyanobacterial communities in SH and XH showed significant variations (PERMANOVA, df = 3, F model = 2.252 and 2.298, R 2 = 0.252 and 0.256, P = 0.011 and 0.02, respectively), whilst those in HS and JFZ did not (Fig. 4B). For bacterial communities, they only showed significant difference at the "Post AB" phase amongst the studied rivers (df = 3, F model = 1.729, R 2 = 0.191, P = 0.039) (Fig. 4C). However, except for HS, bacterial communities in JFZ, SH, and XH collectively revealed significant discrepancy across the cyanobacterial bloom phases (df = 3, F model = 1.884, 3.731, and 4.525, R 2 = 0.220, 0.359, and 0.404, P = 0.006, 0.001, and 0.001 respectively) ( Fig. 4D1-D4). www.nature.com/scientificreports/ The significant difference in cyanobacterial and bacterial dominant species (OTUs). The above results indicated that the discrepancies of both cyanobacterial and bacterial communities could be more attributed to the difference of cyanobacterial bloom phases instead of spatial heterogeneity of rivers. Therefore, the significantly different OTUs of cyanobacterial and bacterial communities were analyzed amongst samples at different phases in each river (containing corresponding tributary pools) ( Table S2). The significance test was subsequently used for LEfSe analysis (Table S3) to identify those cyanobacterial and bacterial OTUs at different phases contributing to the significance. For cyanobacteria, more OTUs were identified in JFZ, SH, and XH, and these OTUs have mainly identified at the rest three phases except the "During AB" phase ( Fig. 5). Cyanobacterial OTUs at different phases were mainly classified into Cyanobacteriales (OTU1843), and Leptolyngbyales (OTU3191) in SH (Fig. 5C), and were pooled into Cyanobacteriales (OTU6199) in XH (Fig. 5D), respectively. Additionally, several OTUs, including OTU2648, OTU7367, OTU8417, and OTU867, were widely distributed amongst these three rivers, indicating that a core cyanobacterial community was shared in common in these rivers.
Correlations between the dominant cyanobacterial taxa and the above-identified dominant species were subsequently analyzed (Fig. 7). In general, the results indicated that some bacterial taxa in different aquatic ecosystems showed different correlations with the same cyanobacterial taxa. For example, OTU853 (Rhodobacter) 25  www.nature.com/scientificreports/ in JFZ was negatively correlated with Cyanobacteriales along with the whole bloom phases (Fig. 7A), however, it showed contrasting correlation patterns with Cyanobacteriales in XH (Fig. 7C). OTU6265 (Porphyrobacter) in JFZ was positively correlated with Cyanobacteriales and Synechococcales for the first three phases, but negatively correlated in the last phase (Fig. 7A), however, it was mostly negatively correlated with these cyanobacterial taxa in SH (Fig. 7B). These discrepancies indicated that the bacteria-cyanobacteria interaction may additionally be affected by environmental parameters.

Correlations of cyanobacterial and bacterial communities with the environmental parameters.
The VIF analysis (Table S4) concluded that all parameters expect DTN, DTP, and COD were selected for cyanobacterial CCA/RDA, whilst all the parameters were used for bacterial CCA/RDA, and analyses were separately performed on microbial communities of each river at different phases (Fig. 8).
In general, cyanobacterial communities in these rivers at different phases were mainly correlated with ORP, pH, TEMP, MPO 4 − , NO 3 − and Chl-a, and communities of each river had a specific correlation with additional parameters (Fig. 8A). For example, cyanobacteria in JFZ also had a close correlation with TP, TN, and NH 4 + (Fig. 8A2), and the ones in SH had a close correlation with TN, NH 4 + , and NO 2 − (Fig. 8A3). Additionally, the cyanobacterial communities of these rivers also exhibited variance in correlations with the main parameters, which indicated the discrepancy of the cyanobacterial communities amongst these rivers.  For bacterial communities in all rivers at different phases, they were mainly correlated with COD, TN, DTN, NO 3 − , ORP, TEMP, and pH, and communities of each river also showed specific correlation with additional parameters (Fig. 8B). Bacterial communities in HS additionally exhibited a close correlation with TP, NO 2 − and Chl-a (Fig. 8B1), the counterparts in JFZ were with NH 4 + and Chl-a (Fig. 8B2), the ones in SH and XH were both with TP, MPO 4 − and NH 4 + (Fig. 8B3), whilst XH showed additional correlation with DTP, and NO 2 − (Fig. 8B4). Interestingly, different from cyanobacteria, the bacterial communities of these rivers demonstrated relatively consistent correlation patterns with the main parameters. For example, communities of "During AB" in these rivers universally showed a positive correlation with pH, TEMP, and TP, and negatively with TN, DTN, NO 3 − , and ORP, except HS, the rest three rivers also collectively showed negative correlations with Chl-a. Furthermore, communities of "During AB" in these rivers consistently exhibited contrasting correlation patterns compared with their counterparts either in "Post AB" or "Before AB".

Discussion
Cyanobacterial bloom is a major factor influencing diversities of microbial community. The significance tests amongst the environmental parameters demonstrated that most of the parameters were not significantly different across the bloom phases (Table S1, P > 0.05), however, both the cyanobacterial and bacterial communities varied not only across the bloom phases in each individual river, but also amongst rivers at each of the same phases (Fig. S2). Additionally, more universal significant differences of cyanobacterial and bacterial communities were determined in rivers of JFZ, SH, and XH across the bloom (Fig. 2B,D), and the community diversity was significantly lower at the "During AB" phase (P < 0.05). The results as a whole indicate that the cyanobacterial bloom is a major factor affecting the diversities of both cyanobacterial and bacterial community.
It has been universally revealed that the main harmful algal blooms could significantly reduce the diversities and functions of other microbial communities, including coast and freshwater planktonic communities 5,6 , microbial communities associated with algae 18,19 , as well as those grew with macrophytes 20 , and sedimentary dwellers 7 . At the "During AB" phase, diversities of cyanobacterial community in JFZ, SH, and XH also significantly reduced, which inferred that the bloom was due to the overproduction of a narrower range of cyanobacterial taxa. Most studies concerning algal blooms focused on Microcystis and related species, which were ubiquitous   www.nature.com/scientificreports/ in life-depending rivers and lakes in China, including Yangtze River, Taihu Lake, and Chaohu Lake [21][22][23][24] . It is therefore of importance to identify the dominant cyanobacterial taxa, their dynamic changes, and how the bacterial community were correspondingly affected along the course of cyanobacterial bloom in rivers of this study.
The distinct spatial and temporal variance of cyanobacterial community composition along the bloom course. Instead of Microcystis, the dominant cyanobacterial taxa were Cyanobacteriales and Synechococcales both in the rivers and their tributary pools (Fig. 3A). Even though these water bodies are close in the locality and have similar environmental conditions, they revealed distinct cyanobacterial community compositions and patterns of dynamic changes in abundances along the bloom phases. Additional beta-diversity analyses on OTU level compared both different rivers at the same phase and each river along the bloom course, the results further illustrated that cyanobacterial community significances were shown in rivers at different bloom phases (Fig. 4A,B), especially for SH, and XH. Beside the bloom influence, the differences of community composition amongst rivers at each phase were also shown, indicating the spatial heterogeneity, which might be due to human activities and riverine characteristics as suggested in other studies 25,26 . Because of the relatively slow flowrate and large surface area, the tributary pools may therefore possess distinct cyanobacterial communities relative to the mainstream rivers, rendering the pools relatively independent ecosystems from rivers to some extent. The dominant cyanobacterial species (OTUs) were further screened out through significance tests and LEfSe analyses on each river at different bloom phases (Fig. 5). The identified dominant species belonged to cyanobacterial genus of Mychonastes, Trachydiscus, Oscillatoriales, Chamaesiphon, Nannochloropsis, Monoraphidium, and Cyanobium, whose ecological functions included feeds for livestock 27,28 , constitutive members of biofilms 29,30 , biodiesel derivation 31,32 , and cyanobacterial early-stage bloom 12,33 . No dominant species were identified from rivers at the "During AB" phase, which may imply that relevant cyanobacterial OTUs in each river were relatively similar in abundance and composition during the bloom.
These results together may shed light on the ongoing treatment campaigns on eutrophication and algal blooms in China, that corresponding strategies should be taken into consideration when dealing with rivers and pools (ponds, reservoirs, and lakes). In general, the cyanobacterial bloom is the major factor that leads to the dynamic variation of community composition, which not only reduces the microbial diversities, assimilating cyanobacterial communities but also would affect the functions of relevant bacterial communities.

Dynamic changes of dominant bacterial community, implications of their ecological functions and cyanobacteria-bacteria interaction.
Collective comparisons of bacterial genera composition in each river at different phases and amongst rivers at the same phases revealed the spatial and temporal heterogeneity of bacterial communities (Fig. 3B-E). The dominant bacterial genus including Acinetobacter, CL500-29, hgcI clade, Limnohabitans, Flavobacterium, Rhodoluna, Porphyrobacter, Rhodobacter, Pseudomonas, and Rhizobiales could also be traced in other relevant studies. For example, in a study concerning Tianmuhu Lake, CL500-29 and hgcI clade were dominant genera and closely associated with cyanobacteria, whilst in the lake's rivers, Flavobacterium, Limnohabitans, and Rhodoluna were the dominant genera 34 . Porphyrobacter, Rhodobacter, Pseudomonas, and Rhizobiales were also dominant concerning either cyanobacterial or Microcystis bloom in different aquatic ecosystems [35][36][37][38][39] .
The LEfSe analysis identified dominant bacterial species (OTUs) in each river at each phase, which demonstrated the dynamic responses of exact bacterial species to the bloom, implying their close interactions with the bloom-forming cyanobacteria (Fig. 6). For example, abundance of Porphyrobacter (OTU6265), Rhodobacter (OTU853 and OTU5851), and Acinetobacter (OTU2460, OTU2184, and OTU2685) were significantly higher at "During AB" phase. It was reported that species of Acinetobacter exhibited algicidal activity against Microcystic aeruginosa 40 , therefore, the outbreak of cyanobacterial bloom may also be likely to give rise to its adverse communities to keep the dynamic balance.
Furthermore, it is believed that the bloom-forming cyanobacteria have profound interaction with the bacterial communities, which has long been the research of interest [9][10][11]41,42 . Cyanobacterial species can act as bactericidal or growth-inhibiting factors, influencing the bacterial communities 18,35,43 , similar as various microbial agents act in other natural environments and artificial constructions 7,44 . On the other hand, various bacterial species, including those identified in this study, have also been reported to be cyanobacterial or growth-inhibiting factors 40,[45][46][47] , which can be further studied to develop microbial reagents to alleviate or control the effects of cyanobacterial bloom. Additionally, many cyanobacteria species are integral members of biofilms in various environments 29,30,48 . These facts jointly reveal that cyanobacteria as a whole are indispensable part of the microbial community. Based upon this cognition, to alleviate or control the cyanobacterial bloom is to restore microbial equilibrium in ecosystems rather than to eliminate cyanobacteria, and understanding about the dynamic changes of cyanobacterial and bacterial communities, as well as those microbial taxa and dominant species relevant to each cyanobacterial bloom phase is taking one key step towards the comprehension of cyanobacteria-bacteria interactions and the behind mechanisms.
Discrepancies of correlations between microbial communities and environmental parameters. The correlation patterns of cyanobacterial communities with the environmental parameters additionally indicated the discrepancies of community composition amongst samples at different bloom phases in each river (Fig. 7). Temperature, nitrogen, and phosphorus were universally considered to be the crucial factors affecting cyanobacterial growth and bloom 2,3,49,50 , and this study showed that environmental parameters, including temperature, nitrogen (ammonium, nitrate, etc.), and phosphorus (MPO 4 − and TP), were also important to the cyanobacterial community compositions [2][3][4]51 , and parameters even with not significant difference (Table S1) could shape distinct cyanobacterial community. www.nature.com/scientificreports/ The cyanobacterial bloom influenced the correlation patterns of bacterial communities of each river with the environmental parameters (Fig. 7B), that bacterial communities at the "During AB" phase generally had contrasting correlation patterns compared with the rest phases. The difference may infer that functional bacterial community on certain environmental parameters were additionally impaired due to the cyanobacterial bloom, and which metabolic pathways were impaired will be included in future studies.

Conclusions
This study investigated the dynamic changes of cyanobacterial and bacterial communities in four upstream rivers, namely Huangshan River (HS), Jinfuzhai River (JFZ), Sihe river (SH), and Xihua River (XH), of a eutrophicated water source reservoir in the subtropical area of China. The major conclusions were drawn as follows.
• The bloom has a major influence on cyanobacterial and bacterial communities, rendering the alpha-diversities of bacterial communities significantly lower at "during cyanobacterial bloom" phase compared with other relevant phases; • The studied rivers commonly shared dominant cyanobacterial taxa, including Cyanobacteriales and Synechococcales, and commonly shared several dominant bacterial genera, including Acinetobacter, Flavobacterium, Rhodobacteraceae, Comamonadaceae, etc. However, the corresponding compositions varied along with the bloom, making cyanobacterial communities significant difference amongst phases in SH and XH, and bacterial communities significant difference amongst phases in JFZ, SH, and XH; • The studied rivers had distinct cyanobacterial and bacterial dominant species (OTUs) along with the bloom, and the bacterial dominant species with significantly higher relative abundances at "during cyanobacterial bloom" phase than other phases, including Rhodobacter, Rhodobacteraceae, Porphyrobacter, Acinetobacter, and Pelomonas, may indicate their close correlations with cyanobacterial bloom in each river; • Cyanobacterial communities at different bloom phases in each river were mainly correlated with ORP, pH, TEMP, MPO 4 − , NO 3 − , and Chl-a, and the bacterial counterparts were mainly with ORP, pH, TEMP, COD, TN, DTN, and NO 3 − . However, communities of different rivers at different phases showed specific correlations with the above parameters.
The dynamic shifts of cyanobacterial and bacterial communities may imply the interactions between cyanobacteria and bacteria, to fully understand and verify the interaction, it is suggested that future studies should expand to network analysis and field experiments investigating the effects of bacterial species identified in this study on the cyanobacterial blooms.

Data availability
The original sequencing data was updated to the NCBI database under the accession number from SAMN23614736 to SAMN23614843. The direct link to the uploaded data is as follows: https:// www. ncbi. nlm. nih. gov/ sra? linkn ame= biopr oject_ sra_ all& from_ uid= 786052.