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

Sampling area and environmental factors determination

This study was carried out along the four upstream rivers, namely Huangshan River (HS), Jinfuzhai River (JFZ), Sihe River (SH), and Xihua River (XH), of Yankou Reservoir basin (29°17'25"~29°18'49" N, 119°54'11"~119°55'13" E) in Yiwu City, Zhejiang province, middle east of China (Fig. 1). Considering the dynamics of an algal bloom, the outbreak course was primarily divided into four phases, including before the algal bloom (Before AB), followed by early phase (Pre AB), outbreak phase (During AB), and late phase (Post AB). Additionally, each river has at least one tributary pool (HS has three pools) connecting to mainstream of the river; therefore, water samples would also be collected from the pools. A total number of 104 surface water samples (0.5 m depth) were collected, respectively in May, July, August, and September 2020, representing the four algal bloom phases. The surface water samples were collected within a 3-day period in all the phases. After immediate collection, the samples were kept in an ice box and transported back to the laboratory for subsequent treatment. Water samples (~3 000 mL) were first passed through a 200-µm pore-size sieve to remove any debris for subsequent determination of water quality parameters and microbial community analysis. For bacterial and cyanobacterial community analyses, water samples (~1 500 mL) were further filtered through 0.22-µm pore-size polycarbonate membranes (50 mm, Jinteng®, Tianjin, China), and those membranes were stored at -80℃ until DNA extraction. For water quality determination, the water samples passed through 200-µm sieve were further filtered through 0.45-µm membranes (50 mm, Jinteng®, Tianjin, China), and the filtration was used subsequently.

The physiochemical parameters of the water samples were determined according to standard methods described previously. Briefly, a Hydrolab HQ30D multiparameter water quality meter (HACH Company Co., Loveland, the USA) was used in situ to monitor the water temperature (TEMP), oxidation reduction potential (ORP), pH value, and chlorophyll-a (Chl-a). Other water quality parameters, including chemical oxygen demand (COD), total nitrogen (TN), dissolved total nitrogen (DTN), ammonium (NH₄⁺), nitrite (NO₂⁻), nitrate (NO₃⁻), total phosphorus (TP), dissolved total phosphorus (DTP) and orthophosphate (MPO₄⁻) were determined using the standard methods published by China’s State Environmental Protection Administration. In total, 13 environmental factors were determined in this study.

DNA extraction and metagenomic sequencing

The total DNA of bacterial and cyanobacterial communities was extracted directly from the membrane filters using a FastDNA^TM Spin Kit (MP Biomedicals, Santa Ana, CA, the USA) according to the manufacturer’s instructions. The purified total DNA was processed and sequenced by Shanghai Majorbio Bio-Pharm Technology Co.,Ltd, and the downstream alignment and analyses were performed using the integrated and free online platform of Majorbio Cloud Platform (www.majorbio.com).

The hypervariable V3-V4 regions of prokaryotic 16S rRNA gene were first amplified through the universal primer set 338F and 806R (Chao et al. 2021). The 20-µL PCR mixture contained 4 µL of 5×FastPfu Buffer, 0.25 mM dNTPs, 10 U of FastPfu polymerase (TransStart®, TransGen Biotech, Beijing, China), 0.2 µM of forward and reverse primers, and 10 ng of the sample DNA. PCR reactions were performed as follows: initial denaturation at 95℃ for 4 min, followed by 30 cycles of 30 s at 94℃, 30 s at 53℃, and 45 s at 72℃, the amplicons of which were subjected to a final extension at 72℃ for 10 min.

Triplicated PCR products for each of 104 samples were examined through gel electrophoresis, and purified using AxyPrepDNA Gel Extraction Kit (Boyao Biotechnology, Shanghai, China). Sequencing libraries were constructed using the TruSeq^TM DNA Sample Prep Kit for Illumina (Illumina Inc., San Diego, CA, the USA) according to manufacturer’s instructions, and index codes were correspondingly supplemented. The quality of constructed library was assessed using the Aglient Bioanalyzer 2100 system and Qubit 2.0 Fluorometer. Finally, the bar-coded amplicons of each sample were mixed and sequenced using an Illumina MiSeq platform (Illumina Inc., San Diego, CA, the USA) according to the manufacturer’s protocols.

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 reads of low-quality. After quality control process, sequences were clustered into operation taxonomic units (OTUs) using UPARSE (v.7.0.1090) with 97% sequence similarity threshold (Edgar 2013). Representative sequence from each OTU were run against the databases of SILVA (Release 138), RDP (11.5), and GreenGenes (Release 135) for reference alignment using the RDP classifier (v.2.11). To avoid mistaken taxonomic alignments, the taxonomic classifications were double-checked against the reference prokaryotes. Those unclassified OTUs were discarded before subsequent analyses, and the sequence data were finally normalized to minimize sequencing biases and allow for appropriate comparison of community variations amongst those samples.

For comparisons of microbial taxa variations across the algal 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 (Chen et al. 2019).

Statistical analyses

After sequencing, standard quality control and data optimization, the sequencing data was normalized for better comparison of variances of microbial communities amongst samples. In this study, the algal bloom was mainly attributed to the cyanobacterial overproduction, therefore, cyanobacterial community would be further extracted from the normalized data for each sample, and alpha- and beta-diversities were analyzed 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 the same river as well as variance of communities at the same phase amongst different rivers was performed.

The alpha-diversity indices including OTU numbers (Sobs), abundance-based coverage estimator (ACE), Shannon, and Chao 1, as well as Venn diagrams compared different samples were calculated using 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 significant level at 0.05. For beta-diversity analysis, the non-metric 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, significance test was first performed using Kruskal-Wallis test coupled with Tukey-Kramer post-hoc 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 biomarkers (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 significant level of 0.05.

Finally, a variance inflation factor analysis (VIF) was first performed to determine which physiochemical parameters would be selected for subsequent correlation analyses between microbial communities and these parameters. Canonical correlation analysis (CCA) or redundancy analysis (RDA) was conditionally chosen based upon results of detrended correspondence analysis (DCA) to reveal the correlations.

Data Availability

The original sequencing data was updated to the NCBI database under the accession number from SAMN23614736 to SAMN23614843.

Environmental parameters of each river at different phases

Water samples were collected from corresponding water bodies, important environmental parameters including water temperature, pH, Chl-a, total nitrogen, and total phosphorus were summarized in Table 1, and complete environmental parameters were summarized in Table S1. In general, the Chl-a in the pools was higher than that in the corresponding rivers, and the maximum concentration reached 535.68µg/L in JFZ’s pool. Further statistical analyses on the complete parameters of each river at different algal bloom phases indicated that for each river, there was no significant difference amongst the parameters at different phases (indicating as N.S.), except ORP and Chl-a (P<0.05, Table S1). For these two parameters, detailed significances between relevant phases were included.

General statistics of sequencing data and alpha-diversity comparison

Comparisons of cyanobacterial and bacterial community composition

In general, most diverse cyanobacterial OTUs were assigned to phylogenetic orders of Chloroplast, 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 three taxa took turns as the most dominant community, for example, in the tributary pool of HS, the most dominant community was Chloroplast, Synechococcales, and Cyanobacteriales across each of the phases; whilst JFZ had Chloroplast and Cyanobacteriales as the most dominant communities. For the cyanobacteria in rivers, Chloroplast was the most dominant taxa in HS river, but Synechococcales was the most dominant in the rest of three rivers along the algal bloom (Fig. 3A).

For bacterial communities, abundant phylogenetic genus to which most diverse OTUs were assigned were sorted out in each river. Based upon 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, community of hgcI clade was always abundant in HS, JFZ, and SH (Fig. 3B-D), whilst was conditionally varied in XH (Fig. 3E), and its relative abundance was greatly reduced during the algal bloom, especially in XH. Communities of Pseudarcicella and Limnohabitans in HS were always abundant and their relative abundance increased during the algal bloom (Fig. 3B), however, their relative abundances were conditionally varied and noticeably reduced in the rest three rivers during the same phase (Fig. 3C-E). Similarly, Comamonadaceae was always abundant in each river, but its relative abundance clearly reduced during the algal bloom, and Allobaculum was always rare in HS, SH, and XH (Fig. 3B, C, E), but conditionally varied in JFZ (Fig. 3D) that its relative abundance greatly increased during the algal bloom.

Variations of cyanobacterial and bacterial communities were further illustrated through NMDS with PERMANOVA test. Cyanobacterial communities at the “Before AB” phase exhibited significant variation amongst the studied rivers (PERMANOVA, df = 3, F model = 3.220, R² = 0.305, P = 0.001), however not at the rest three phases (Fig. 4A1-A4). Across the algal bloom phases, cyanobacterial communities in SH and XH showed significant variations (PERMANOVA, df = 3, F model = 2.252 and 2.298, R² = 0.252 and 0.256, P = 0.011 and 0.02, respectively), whilst those in HS and JFZ did not (Fig. 4B1-B4). For bacterial communities, they only showed significant difference at the “Post AB” phase amongst the studied rivers (df = 3, F model = 1.729, R² = 0.191, P = 0.039) (Fig. 4C1-C4). However, except for HS, bacterial communities in JFZ, SH, and XH collectively revealed significant discrepancy across the algal bloom phases (df = 3, F model = 1.884, 3.731, and 4.525, R² = 0.220, 0.359, and 0.404, P = 0.006, 0.001, and 0.001, respectively) (Fig. 4D1-D4).

Significant difference in cyanobacterial and bacterial biomarkers (OTUs)

Above results indicated that the discrepancies of both cyanobacterial and bacterial communities could be more attributed to the difference of algal bloom phases instead of the spatial heterogeneity of rivers. Therefore, the significantly different OTUs of cyanobacterial and bacterial communities were analyzed amongst samples at different phases in each individual 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 were mainly identified at the rest three phases except the “During AB” phase (Fig. 5). Cyanobacterial OTUs at different phases were mainly classified into Chloroplast (OTU2648, OTU7367, OTU8492, OTU6617, OTU5009, and OTU867) and Synechococcales (OTU7163, OTU6821, OTU8339, OTU8374, and OTU6991) in JFZ (Fig. 5B), whilst the counterparts were mainly categorized into Chloroplast (OTU8417, OTU7367, OTU2648, and OTU983, etc.), Cyanobacteriales (OTU1843), and Leptolyngbyales (OTU3191) in SH (Fig. 5C), and were pooled into Chloroplast (OTU8492, OTU 6675, OTU7166, OTU2082, and OTU983, etc.) and Cyanobacteriales (OTU6199) in XH (Fig. 5D), respectively. Additionally, several OTUs, including OTU2648, OTU7367, OTU8417, and OTU867, were found to be widely distributed amongst these three rivers, which indicated that they might share a core cyanobacterial community.

Similarly for bacterial community, more OTUs were identified in JFZ, SH, and XH, however were universally identified across the algal bloom phases (Fig. 6). Members of Rhodobacteraceae (OTU 853 and OTU 5851), Limnohabitans (OTU4177, OTU216, and OTU3910), Porphyrobacter (OTU6265), and Sporichthyaceae (OTU2594) were mainly screened out in JFZ (Fig. 6B), whilst Pseudarcicella (OTU3296), Comamonadaceae (OTU3244), Sporichthyaceae (OTU2594), Limnohabitans (OTU3910), Acinetobacter (OTU8318, OTU2460, OTU2685, OTU5230, OTU2703, OTU5164, OTU483, and OTU2184) were broadly identified in common in SH and XH (Fig. 6C and D). Notably, members of Rhodobacteraceae (OTU853 and OTU5851), Porphyrobacter (OTU6265), and Acinetobacter (OTU2460, OTU2685, OTU2184, and OTU5230) in JFZ, SH, and XH exhibited significantly higher (P<0.05) relative abundances at “During AB” phase than the rest three phases (Fig. 6B-D), which may indicate their relationships with the cyanobacterial communities causing algal bloom. Correlations between the dominant cyanobacterial taxa and the above identified biomarkers were subsequently analyzed (Fig. 7). In general, the results indicated that same bacterial taxa in different aquatic ecosystem showed different correlations with the same cyanobacterial taxa. For example, OTU853 in JFZ was negatively correlated with Chloroplast and Cyanobacteriales along the whole bloom phases (Fig. 7A), however, it showed contrasting correlation patterns with Chloroplast and Cyanobacteriales in XH (Fig. 7C). OTU6265 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 on correlation analyses indicated that the bacteria-cyanobacteria interaction may additionally be affected by environmental factors.

Correlations of cyanobacterial and bacterial communities with the physiochemical parameters

The VIF analysis (Table S4) on the physiochemical parameters against the cyanobacterial and bacterial communities 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₄⁻, NO₃⁻ and Chl-a, and communities of each individual river had specific correlation with additional parameters (Fig. 8A1-A4). For example, cyanobacteria in JFZ also had close correlation with TP, TN, and NH₄⁺ (Fig. 8A2), and the ones in SH had close correlation with TN, NH₄⁺, and NO₂⁻ (Fig. 8A3). Additionally, cyanobacterial communities of these rivers also exhibited variance on correlations with the main parameters, which indicated the discrepancy of the cyanobacterial communities amongst these rivers. For instance, most of the samples of HS across algal bloom were positively correlated with ORP, TN, NH₄⁺, NO₃⁻, MPO₄⁻, TP, TN, and were either not or negatively correlated with the rest main parameters, including pH, TEMP, NO₂⁻, and Chl-a (Fig. 8A1), whilst most samples in JFZ but the ones of “Before AB” showed positive correlation with TEMP, ORP, NH₄⁺, and MPO₄⁻, and negative correlation with TN, NO₃⁻, TP, pH, and Chl-a (Fig. 8A2). Cyanobacterial communities in SH were universally negatively correlated with the main parameters (Fig. 8A3), whilst the ones in XH largely illustrated positive correlation with TEMP, pH, NO₃⁻, and MPO₄⁻ (Fig. 8A4). The correlation between XH communities and the main parameters were reflected through RDA according to the DCA analysis.

For bacterial communities in all rivers at different phases, they were mainly correlated with COD, TN, DTN, NO₃⁻, ORP, TEMP, and pH, and communities of each individual river also showed specific correlation with additional parameters (Fig. 8B1-B4). Bacterial communities in HS additionally exhibited close correlation with TP, NO₂⁻ and Chl-a (Fig. 8B1), the counterparts in JFZ were with NH₄⁺ and Chl-a (Fig. 8B2), the ones in SH and XH were both with TP, MPO₄⁻ and NH₄⁺, whilst XH showed additional correlation with DTP, and NO₂⁻ (Fig. 8B3&B4). 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 positive correlation with pH, TEMP, and TP, and negatively with TN, DTN, NO₃⁻, and ORP, except HS, the rest three rivers also collectively showed negative correlations with Chl-a (Fig. 8B1-B4). Furthermore, communities of “During AB” in these rivers consistently exhibited contrasting correlation patterns compared with their counterparts either in “Post AB” or “Before AB” (Fig. 8B1-B4).

Algal bloom is a major factor influencing diversities of microbial community

Although the diversity indices indicated that both the cyanobacterial and bacterial communities varied not only across the algal bloom phases in each individual river of HS, JFZ, SH, and XH, but also amongst rivers at each of the same phases (Fig. S2), the significance tests clearly demonstrated that comparing with other factors (Table S1) the algal bloom was a major factor influencing the diversities of both cyanobacterial and bacterial communities (Fig. 2). Only several cases of significant differences of diversities were determined between rivers at the same phases, such as HS and JFZ at “Before AB” phase (Fig. 2A1) as well as JFZ and SH at “Post AB” phase (Fig. 2A4) for cyanobacteria, and JFZ and XH at “Post AB” phase for bacteria (Fig. 2C4), more universal significant difference cases were shown in rivers of JFZ, SH, and XH across the algal bloom, which all indicated that both cyanobacterial and bacterial communities were significantly lower during algal bloom (Fig. 2B2-B4 and D2-D4). It has been universally revealed that the main harmful algal blooms could significantly reduce the diversities and functions of other microbial communities, including planktonic communities of coast and freshwater (Chai et al. 2020, Amorim and Moura 2021), microbial communities associated with algae (Zhou et al. 2020, Zhu et al. 2021), as well as those grew with macrophytes (Jiang et al. 2019), and sedimentary dwellers (Yang et al. 2021). During the algal bloom phase, cyanobacterial community diversities in JFZ, SH, and XH also significantly reduced, which inferred that the algal bloom was due to the overproduction of a narrowed range of cyanobacterial taxa. Most studies concerning algal blooms focused on Microcystis and related species, which were reported ubiquitous in life-depending rivers and lakes of China, including Yangtze River, Taihu Lake, and Chaohu Lake (Huang et al. 2018, Shi et al. 2018, Tang et al. 2018, Huo et al. 2021), it is of importance to identify the dominant cyanobacterial taxa, their dynamic changes, and how the bacterial community compositions were correspondingly affected along the course of algal bloom in rivers of this study.

Distinct spatial and temporal variance of cyanobacterial community composition along the algal bloom course

Instead of Microcystis, the dominant cyanobacterial taxa were Chloroplast, Cyanobacteriales, and Synechococcales both in the concerned rivers and their tributary pools (Fig. 3A). Even though these water bodies are close in locality and have similar environmental conditions, they revealed distinct cyanobacterial community compositions and patterns of dynamic changes in relative abundances along the algal bloom phases. Additional beta-diversity analyses on OTU level compared both different rivers at the same phase and each individual river along the algal bloom course, the results further illustrated that cyanobacterial community significances were shown in rivers at different algal phases (Fig. 4A1-B4), especially for SH, and XH. Beside the algal bloom influence, the differences of community composition amongst rivers at each corresponding phase were also shown (Fig. 4A1). The spatial heterogeneity might be due to human activities and riverine characteristics as suggested in other studies (Su et al. 2014, Mao et al. 2015). 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 to some extent.

The cyanobacterial biomarkers (OTUs) were further screened out through significance tests and LEfSe analyses on each river at different algal bloom phases (Fig. 5). The identified biomarkers belonged to cyanobacterial taxa of Mychonastes, Trachydiscus, Oscillatoriales, Chamaesiphon, Nannochloropsis, Monoraphidium, and Cyanobium, whose ecological functions included feeds for livestock (Rezanka et al. 2015, Saadaoui et al. 2020), constitutive members of biofilms (Nowicka-Krawczyk and Zelazna-Wieczorek 2017, Leiser et al. 2021), biodiesel derivation (Holbrook et al. 2014, Shi et al. 2021), and cyanobacterial early-stage bloom (Li et al. 2020, Fernandes et al. 2021). No biomarkers were identified from the “During AB” phase of the rivers, which may imply that relevant cyanobacterial OTUs in each river were relatively similar in abundance and composition during the algal bloom.

These results together may shed lights 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 algal 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 and implications concerning their ecological functions

The dominant bacterial genera along the algal bloom course in each river were identified and categorized into groups of AAT, CVT, AMT, and ART based upon dynamic changes of their relative abundance along the course (Fig. 3B-E). Comparisons of bacterial genera composition in each river at different phases and amongst rivers at the same phases collectively revealed the spatial and temporal heterogeneity of bacterial communities. Similar with cyanobacteria, the beta-diversity analyses illustrated that bacterial community significances were more apparent in rivers along the algal bloom course (Fig. 4C1-D4), which additionally proved the major influence imposed by algal bloom. Furthermore, the identified bacterial significant biomarkers in each river at each phase (Fig. 6) further demonstrated the responses of exact members of the dominant bacterial taxa along the algal bloom course, indicating their close interactions with the algal bloom.

The dominant bacterial taxa identified along the algal bloom course included Acinetobacter, CL500-29, hgcI clade, Limnohabitans, Flavobacterium, Rhodoluna, Porphyrobacter, Rhodobacter, Pseudomonas, and Rhizobiales (Fig. 3B-E), which could also be traced in other relevant studies. In a study concerning Tianmuhu Lake, CL500-29 and hgcI clade were found as dominant genera, which also closely associated with cyanobacteria, whilst in the lake’s river system, Flavobacterium, Limnohabitans, and Rhodoluna were identified as the dominant genera (Xie et al. 2021). Porphyrobacter, Rhodobacter, Pseudomonas, and Rhizobiales were also found dominant in other studies concerning either cyanobacterial bloom or Microcystis bloom in different aquatic ecosystems (Shi et al. 2017, Li et al. 2019, Kim et al. 2020, Pineda-Mendoza et al. 2020, Liu et al. 2021). Additionally, it was reported that species of Acinetobacter exhibited algicidal activity against Microcystic aeruginosa (Su et al. 2016), Flavobacterium and Limnohabitans also showed cyanobactericidal and growth-inhibiting activities (Li et al. 2019, Liu et al. 2021), and Rhizobiales could degrade the cyanobacteria-derived particulate organic matters (Shi et al. 2017), all of which could be used as biological measures to alleviate the ecological influences of algal bloom. The relative abundances of Limnohabitans were significantly reduced at “During AB” phase in JFZ, SH, and XH laterally proved its cyanobacteria-inhibiting activity, and noteworthily, Pseudarcicella, Sporichthyaceae, and other taxa identified in this study may also serve the similar function.

Discrepancies of correlations between microbial communities and environmental parameters

In general, the correlation patterns of cyanobacterial communities with the environmental parameters additionally indicate the discrepancies of community composition amongst samples at different algal bloom phases in each river (Fig. 8A1-A4). The environmental parameters, including temperature, nitrogen, and phosphorus, were universally considered to be the crucial factors affecting algal growth and bloom (Paerl et al. 2011, Xu et al. 2015, Li et al. 2021, Liu et al. 2021). This study additionally showed that environmental parameters even with minor difference (Table S1) could shape distinct cyanobacterial community composition, which could be categorized as nitrogen-limited and phosphorus-limited growth types (Xu et al. 2010, Paerl et al. 2011, Paerl et al. 2014, Xu et al. 2015).

The algal bloom clearly influenced the correlation patterns of bacterial communities of each river with the environmental parameters (Fig. 8B1-B4), that bacterial communities during the algal bloom phase generally had contrasting correlation pattern compared with the rest phases. The difference may infer that functional bacterial community on certain environmental parameters were additionally impaired due to the algal bloom, and which metabolic pathways were impaired will be included in future studies.

The interaction between cyanobacteria and bacteria

This study reveals the dynamic changes of cyanobacteria and bacteria along the algal bloom, and it is believed that the algal bloom forming cyanobacteria have profound interaction with the bacterial communities, which has long been the research of interest (Kouzuma and Watanabe 2015, Ramanan et al. 2016, Cirri and Pohnert 2019, Jiang et al. 2021, Sial et al. 2021). Cyanobacterial species can act as bactericidal or growth-inhibiting factors, influencing the bacterial communities (Su et al. 2017, Kim et al. 2020, Zhou et al. 2020, Carpine and Sieber 2021), similar as various microbial agents act in other natural environments and artificial constructions (Feng et al. 2016, Yang et al. 2021). On the other hand, various bacterial species, including those included in this study, have also been reported to be cyanobactericidal or growth-inhibiting factors (Su et al. 2016, Zheng et al. 2018, Wang and Coyne 2020, He et al. 2021), which are further studied to develop microbial reagents to alleviate or control the effects of algal bloom. Additionally, many cyanobacteria species are integral members of biofilms in various environments (Nowicka-Krawczyk and Zelazna-Wieczorek 2017, Xia et al. 2020, Leiser et al. 2021). These facts jointly reveal that cyanobacteria as a whole is the indispensable part of microbial community. Based upon this cognition, to alleviate or control the algal 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 significant biomarkers relevant to each algal bloom phase is taking one key step towards the comprehension of cyanobacteria-bacteria interactions and the behind mechanisms.