Synchronous dynamics and correlations between bacteria and phytoplankton in a subtropical drinking water reservoir

Both phytoplankton and bacteria are key and abundant components of aquatic ecosystems and play pivotal roles in maintaining ecosystem structure and function. However, the extent to which phytoplankton community succession inﬂu-ences changes in bacterial community composition (BCC) is largely unknown. In this study, we evaluated the correlations between bacteria and phytoplankton communities and determined the relative contribution of phytoplankton community succession to temporal variation of BCC in a subtropical drinking water reservoir (Tingxi Reservoir, southeast China). Bacterial communities were investigated by quantitative PCR and 454 pyrosequencing of 16S rRNA genes, while phytoplankton communities were analyzed by light microscopy. A remarkable seasonal succession from Cyanophyta to Bacillariophyta was observed during the study period, and this succession can accurately predict the distribution and abundance of the bacterial OTUs based on the discriminant function analysis. Association networks revealed that 38 of the 46 abundant bacterial OTUs exhibited signiﬁcant correlations with phytoplankton. More interestingly, the positive correlations dominated the associated network, which may suggest that facilitative correlations between phytoplankton and bacteria are more important than inhibitory correlations in the Tingxi Reservoir. In addition, some bacterial OTUs were closely correlated with the dynamics of Microcystis , and they were afﬁliated with the divisions Acidobacteria , Actinobacteria , and Proteobacteria . Structural equation model showed that succession of phytoplankton community explained the largest part of temporal variation in BCC. Therefore, our data suggest that the distinct succession of phytoplankton community may mediate the temporal dynamics of bacterial community in the Tingxi Reservoir.


Introduction
Bacteria and phytoplankton constitute the base of aquatic microbial food webs and play essential roles in the structure and function of aquatic ecosystems (Reynolds, 1984;Newton et al., 2011;Dickerson & Williams, 2013). In aquatic microbial food webs, bacteria and phytoplankton can interact with each other. These interactions can include mutualism, competition and parasitism, influencing the microbial abundance, diversity, community composition, and ecological function (Cole, 1982;Amin et al., 2012). Previous work has proposed the hypothesis that phytoplankton succession may influence the dynamics of the bacterial community, as different phytoplankton communities lead to difference in quantity and quality of organic matter which the bacterial community can utilize . Various studies attempted to test this hypothesis, for example, Kent et al. (2007) using automated ribosomal intergenic spacer analysis (ARISA), demonstrated that the bacterial community composition (BCC) had tight correlations with the succession of the phytoplankton community in six humic lakes. Niu et al. (2011) observed that several specific bacterioplankton groups tightly coupled to different phytoplankton blooms or dynamics and suggested that the succession of the phytoplankton community was a major factor that drove the dynamics of bacterioplankton community in Lake Taihu. Recently, Paver et al. (2013) investigated the correlations between bacteria and phytoplankton in three humic lakes, and their mesocosm experiment showed that the tight interactions between specific phytoplankton and bacteria affected the temporal dynamics of bacterial community. However, there are potential limitations in these studies. First, PCR-based fingerprinting which was used in these studies finds it difficult to amplify microorganisms with abundances in the ecosystem of < 1% of the total community (Pedr os-Ali o, 2012), and some DGGE bands or ARISA fragment lengths can represent multiple, diverse taxa (Diez et al., 2001;Newton et al., 2006). Second, Paver et al. (2013) investigated the correlations of bacteria with phytoplankton using bacterial relative abundance data, which may lead to spurious correlations. An apparent decrease in a population may in fact be the opposite if the total communities increase in size at the same time.
Most importantly, these studies did not directly test the relative contribution of phytoplankton community succession to temporal variation of BCC. To better understand the relationships of bacteria and phytoplankton, in this study, we used quantitative real-time PCR and 454 pyrosequencing to analyze the bacterial community and used classical microscopy to investigate the phytoplankton community. The relative contribution of phytoplankton community succession was tested using structural equation model (SEM). In addition, the specific relationships between these two communities were performed using association network at species level. This approach has been proposed for exploring the complex interactions that are likely to take place among microorganisms (Ruan et al., 2006;Fuhrman & Steele, 2008;Fuhrman, 2009;Steele et al., 2011). However, the studies published to date have focused on the bacteria, archaea, and protists. There is little information on freshwater bacteria and phytoplankton. The aims of this study were (1) to investigate the temporal dynamics of bacterial and phytoplankton communities in a subtropical drinking water reservoir; (2) to examine the possible relationships between bacteria and phytoplankton at species level; and (3) to assess the relative importance of environmental and phytoplankton community factors affecting the temporal variation of BCC.

Study area and sampling
Tingxi Reservoir is located in the north of Xiamen city, southeast China. This reservoir was completed in 1956 with a drainage area of 100.8 km 2 and a total storage capacity of 48.45 million m 3 . It is one of the most impor-tant sources for drinking water supply of Xiamen city. The area has a subtropical humid monsoon climate with an annual mean precipitation of 1600 mm and an annual mean temperature of 21°C. Twelve surface water samples (upper 50 cm) were collected monthly in the middle of the reservoir (24°48 0 N, 118°8 0 E) from May 2010 to April 2011. Bacterial communities (400 mL water) were collected on 0.22-lm-pore-size polycarbonate filters (47 mm diameter, Millipore). The filters were stored at À80°C until DNA extraction. For the phytoplankton analysis, a total of 2.5 L surface water samples were fixed in Lugol's solution and passed through a series of settling and siphoning steps to produce a 25 mL concentrate (Lv et al., 2014).

Physico-chemical analysis
Water temperature, pH, dissolved oxygen, electrical conductivity (EC), and chlorophyll a were measured in situ with a Hydrolab DS5 multi-parameter water quality analyzer (Hach). Water transparency was determined with a 30-cm Secchi disk. Total nitrogen (TN), ammonium nitrogen (NH 4 -N), nitrite and nitrate nitrogen (NO x -N), total phosphorus (TP), and phosphate phosphorus (PO 4 -P) were measured following methods used in our previous study (Liu et al., 2013).

Phytoplankton analysis
Phytoplankton were identified and counted under an inverted microscope (Motic, China) according to Shen & Zhang (1990), Zhang & Huang (1991), and Hu & Wei (2006). Three subsamples were investigated for each sample, and more than 500 individuals were counted for each sample. Biomass was estimated from cell numbers and cell size measurements, assuming that 1 mm 3 of algal volume equal 1 mg of fresh weight biomass (Lv et al., 2014).

DNA extraction and quantitative real-time PCR
Total DNA was extracted directly from the filter using a Fast DNA spin kit (Bio101) according to the manufacturer's instructions. The extracted DNA was dissolved in 50 lL TE buffer, quantified by spectrophotometer, and stored at À20°C until further use.
The partial 16S rRNA genes were amplified to prepare the quantitative PCR standard curve with the primers 341F (5 0 -CCTACGGGNGGCWGCAG-3 0 ) (Herlemann et al., 2011) and 515R (5 0 -ATTCCGCGGCTGGCA-3 0 ) (L opez- Guti errez et al., 2004). The 50-lL PCR mixture contained 1 lL of the primer set (25 pmol each), 0.5 lL (1.25 U) of ExTaq DNA polymerase (Takara, Japan), 5 lL of Ex Taq buffer, 4 lL of deoxyribonucleotide triphosphate (dNTP) mixture (2.5 mM each), 100 ng of DNA template, and 33.5 lL of ddH 2 O. The PCR program began with a 5 min denaturation at 94°C; this was followed by 25 cycles of 94°C for 30 s, 51°C for 30 s, and 72°C for 60 s. The final cycle was extended at 72°C for 10 min. The purified PCR products (Promega) were cloned into the pMD18vector (Takara) and transformed into Escherichia coli DH5a competent cells (Takara). The successfully inserted plasmids were sequenced using an automated sequencer (ABI3730). The plasmids DNA then were extracted using the MiniPrep kit (Qiagen, Hilden, Germany), and the plasmid concentrations were determined by spectrophotometry using a BioPhotometer (Eppendorf, Germany). Standard curves were prepared from linearized plasmid serial dilutions containing between 10 4 and 10 10 16S rRNA genes copies calculated directly from the concentration of extracted plasmid. A standard curve was generated by plotting the threshold cycle values vs. log 10 of the gene copy numbers. The amplification efficiency (E) was estimated using the slope of the standard curve through the following formula: E = (10 À1/slope ) À1.
The quantitative PCR was carried out in a volume of 20 lL comprising 10 lL SYBR Premix ExTaq TM (Takara), 0.25 lM of each primer (341F and 515R), 2 lL of total DNA, and RNase-free water. The thermocycling programs of quantitative PCR included an initial 30 s at 95°C, followed by 40 cycles of 5 s at 95°C and 34 s at 60°C. All the measurements were taken in triplicate. In this study, a standard curve (r 2 = 0.991) was used as the reference to calculate the concentrations of environmental 16S rDNA. The efficiency of PCR amplification of the 16S rRNA genes was 105.8%. The amplified products were visualized on a BioRad Gel Documentation System after running in 1% agarose gels.

pyrosequencing
The V3-V5 region of bacterial 16S rRNA genes were amplified and sequenced by 454 pyrosequencing. Total DNA was sent to the Personal Biotechnology Co., Ltd in Shanghai, China, for high-throughput sequencing on a Roche-FLX+ pyrosequencer. The V3-V5 region of 16S rRNA genes were amplified with the primers 357F (5 0 -CCT ACG GGA GGC AGC AG -3 0 ) (Muyzer et al., 1998) and 926R (5 0 -CCG TCA ATT CMT TTR AGT -3 0 ) (Lane, 1991). Each DNA sample was individually PCR amplified in triplicated 25 lL reactions included an initial denaturation at 94°C for 4 min, followed by 25 cycles of 30 s at 94°C, 45 s at 50°C, and 60 s at 72°C. At the end of the amplification, the amplicons were subjected to final 8 min extension at 72°C. Each reaction contained 19 PCR buffer, 2.5 mM dNTPs, 0.625 U of Taq DNA polymerase, 10 lM of each primer, and 20 ng of target DNA.

Sequence analysis
Raw sequence data were processed using QIIME software package (Caporaso et al., 2010). Sequences were quality controlled using the Split_Libraries.py with following settings: 200 < sequence length < 1000, mean quality > 25, ambiguous bases < 1, and homopolymer length < 6. The remaining sequences were then clustered using the pick_otus.py script with the Uclust method (97% sequence similarity) (Fortunato et al., 2012). The sequencing data were denoised by running the MOTHUR implementation of the AmpliconNoise algorithm (shhh.seqs) (Schloss et al., 2009;Quince et al., 2011). Potentially chimeric sequences were identified and removed by the Chimera Slayer using a 50 000-column wide SILVA reference alignment (Schloss et al., 2009;Haas et al., 2011), and a total of 7044 sequences composing 584 OTUs were eliminated from the data set. Further, the OTUs which contained only 1 read were not used to avoid possible biases. We randomly selected identical sequences from each sample to standardize sequencing effort across samples (Zinger et al., 2012). Then, we constructed a bacterial absolute abundance matrix, based on the bacterial relative abundance 9 16S rRNA gene copy number. Finally, all reads identified as Archaea or chloroplasts were removed, and to avoid artificially increasing the correlations of Cyanophyta between microscopy and 454 sequencing, all Cyanophyta OTUs of 454 sequencing were eliminated.
Six bacterial OTUs were significantly correlated with the dynamics of Microcystis flosaquae based on the extended local similarity analysis (eLSA). The represent sequences of these OTUs were aligned with the CLUSTAL W program and compared with sequences available in Gen-Bank databases using BLAST (Thompson et al., 1994).

Data analysis
Two Bray-Curtis similarity matrices were constructed with the bacterial absolute abundance data and phytoplankton absolute biomass data generated from each site, respectively. The nonmetric multidimensional scaling (MDS) ordination was used to investigate differences in microbial communities among samples (Clarke & Gorley, 2001). A measure of goodness of fit of the ordination is given by a stress value (Kruskal's stress formula) that should be < 0.20 to minimize misinterpretation (Kruskal, 1964). To evaluate the significant differences between groups, we used the randomization/permutation procedure analysis of similarities (ANOSIM). The ANOSIM statistic R is calculated by the difference of the between-group and within-group mean rank similarities, thus it displays the degree of separation between groups. Complete separation is indicated by R = 1, whereas R = 0 suggests no separation (Clarke & Gorley, 2001). The RELATE works by relating the bacterial community similarity matrix to phytoplankton community similarity matrix, by calculating rank correlation between these two matrices, thus provides a significance test with the matching coefficient qm, which is equivalent to the Mantel's test (Clarke & Gorley, 2001). We referred to positive and negative correlation levels between 0.5 and 1 as strong relations at P < 0.01.
The Cyanophyta dominant period (May and June 2010) and non-Cyanophyta dominant period (July 2010 to April 2011) were arbitrarily defined based on the absolute and relative biomass of Cyanophyta. Similarly, the Bacillariophyta dominant period (October 2013-January 2014) and non-Bacillariophyta dominant period (May-September 2010 and February-April 2011) were arbitrarily defined based on the absolute and relative biomass of Bacillariophyta. To display the differences between the different periods (i.e. Cyanophyta dominant period vs. non-Cyanophyta dominant period), discriminant function analysis (DFA) was performed on the abundant bacterial OTUs or the seven phytoplankton phyla (Manly, 2005). We use discriminant scores as our measure of bacterial diversity that best distinguished different periods . Stepwise multiple regression analyses (MRA) were performed to examine the relationship between discriminant scores (bacterial diversity) and each groups of biological and physicochemical variables (Scheiner & Gurevitch, 1993). We used the software package SPSS version 19 for these analyses.
We analyzed the associations between abundant bacteria and phytoplankton using eLSA based on bacterial absolute abundance and phytoplankton absolute biomass (http://meta.usc.edu/softs/lsa/) (Ruan et al., 2006;Xia et al., 2011). Abundant bacteria were defined as the OTUs with a representation ≥ 1% within a sample, and abundant phytoplankton were defined as the species with a representation ≥ 1% within a sample (Pedr os-Ali o, 2012). Finally, we included 74 variables in the analysis: 28 phytoplankton species and 46 bacterial OTUs. Parameters for the eLSA were set to 1000 permutations. We used CYTO-SCAPE 2.6.3 (Shannon et al., 2003) to create networks showing the 62 variables with 144 significant local similarity correlations between bacteria and phytoplankton (P < 0.01) with a false discovery rate (that is, Q-value) < 0.004.
We used a forward selection procedure ('ordistep' function from package 'VEGAN') to select physico-chemical variables that related significantly (P < 0.05) with the variation of BCC and phytoplankton community composition (Blanchet et al., 2008). The forward selection procedure analysis was performed using the R software. Directed dependencies between the response variables (BCC and phytoplankton community composition) and all groups of significant relevant contextual physicochemical parameters were assessed in a SEM using path analysis (Legendre & Legendre, 1998;Shipley, 2002;Bienhold et al., 2012). The Spearman coefficient was used to derive a correlation matrix between groups of variables. We started with an initial model that included all plausible pathways between BCC, significant relevant contextual physico-chemical parameters (water temperature, TP, and EC), and phytoplankton community composition. Subsequently, the significance of each path coefficient was tested by its critical ratio (P < 0.05), and nonsignificant paths were removed in a stepwise fashion until all remaining paths were significant (Kline, 2005). The overall fit of the final model was evaluated with the goodnessof-fit index (GFI), the adjusted goodness-of-fit index (AGFI), and the Bentler comparative fit index (CFI) (Schermelleh-Engel et al., 2003). The structural equation analysis was performed using the software package AMOS version 17.

Nucleotide sequence accession number
All sequence data from this study have been deposited in the public NCBI database (http://www.ncbi.nlm.nih.gov/) under the accession number SRX400186. The represent sequences of the six bacterial OTUs that were significantly correlated with the dynamics of M. flosaquae based on the eLSA have been deposited to GenBank under the Accession Number KF984315-KF984318, KF984320 and KJ872512.

Environmental and biological variables
All of the 13 environmental and biological parameters exhibited wide ranges and showed diverse seasonal patterns (Table 1, Fig. 1).
The highest concentrations of TP, EC and chlorophyll a and the lowest transparency occurred in June, while the maximum copy number of 16S rRNA genes and TN, and minimum dissolved oxygen were detected in July. In this period, the water temperature was high.
The minimum TP and EC was detected in January, while the maximum DO was detected in February. The lowest concentrations of chlorophyll a and the maximum transparency occurred in December. During this period, the water temperature was relatively low.
The biomass of Cyanophyta was high in May and June 2010 and then decreased dramatically. In contrast, the biomass of Bacillariophyta was low from May to August and then increased remarkably. The biomass of Bacillariophyta achieved a high value from October 2010 to January 2011 and then decreased substantially (Fig. 1).

Synchronous succession of bacterial and phytoplankton communities
In general, there was no significant variation in bacterial communities at phylum level. Actinobacteria was the most dominant group in all twelve samples (the relative abundance ranged from 50.2% to 74.1%) during our sampling period, followed by Proteobacteria, unclassified bacteria, and Bacteroidetes (Fig. 2). In contrast, phytoplankton communities showed a distinct seasonal shift at phylum level. Cyanophyta and Bacillariophyta were the two major dominant phyla in the Tingxi Reservoir. Cyanophyta dominated the phytoplankton communities from May to June 2010, while Bacillariophyta dominated the communities from September 2010 to April 2011 based on the relative biomass (Fig. 2). 2.02 9 10 9 1.09 9 10 8 1.04 9 10 10 At the community level, however, the dynamics of bacterial and phytoplankton communities were strongly concordant according to the nonmetric MDS (Fig. 3). Interestingly, all samples were grouped into two clusters for both bacterial and phytoplankton communities. One cluster included five samples from May to September 2010 and associated with the warmer season and higher nutrients. The other cluster, from the colder season with lower nutrients, consisted of seven samples from October 2010 to April 2011. The ANOSIM analyses revealed a significant separation of these two clusters (for bacterial community R = 0.524, P = 0.003; for phytoplankton community R = 0.873, P = 0.001). The RELATE analysis showed the matching coefficient qm of 0.685 between bacterial and phytoplankton communities at P = 0.001, indicating their community compositions were significantly correlated.
The discriminate functions can accurately predict the months of phytoplankton blooms, as Cyanophyta and Bacillariophyta dominated the calculation of the discriminant function. Also, the bloom of Cyanophyta and Bacillariophyta can accurately predict the distribution and abundance of the bacterial OTUs. MRA showed that the bacterioplankton diversity, as indexed by DFA scores, was consistently predictable from the phytoplankton biomass and environmental variables (Table 2).

Association network
The number of positive correlations between bacteria phylotypes and phytoplankton species was 130, while the number of negative correlations was 14 (Supporting Information, Fig. S1). In addition, six bacterial OTUs (Acidobacteria OTU1610, Actinobacteria OTU764, Actinobacteria OTU2937, Betaproteobacteria OTU787, Gammaproteobacteria OTU1015 and OTU2930) were significantly correlated with the dynamics of a common bloom- forming species, M. flosaquae (Fig. 4). According to the BLAST analyses, these six OTUs were closely related to Microcystis species (Table 3).

Drivers of temporal dynamics of bacterial community
Water temperature, TP, and EC were significantly related to the variation of BCC by forward model selection (P < 0.05). Both water temperature and TP were significantly correlated with the succession of phytoplankton community (P < 0.05).
To gain further understanding of the linkages between BCC, significant relevant contextual environmental parameters and phytoplankton community succession, we developed a structural equation model (SEM). Our initial model included all plausible path ways between BCC, significant environmental parameters and phytoplankton community composition (Fig. S2). We trimmed the initial model by removing all nonsignificant pathways, until we arrived at a final model in which all path ways were significant. The AGFI of the trimmed model was 0.93, indicating a good fit of the model to the original data.
TP concentration and succession of phytoplankton community had direct effects on BCC. However, water temperature did not display a significant direct effect on BCC, but had indirect effect on BCC due to their effect on succession of phytoplankton community (Fig. 5).

Discussion
The MDS ordination revealed a distinct succession in the phytoplankton community from May 2010 to April 2011. Interestingly, the bacterial community yielded a similar seasonal succession pattern when compared to the phytoplankton community. This indicated that there were potential correlations between the bacteria and phytoplankton taxa. Although previous studies have noted these links between bacterial and phytoplankton communities using aggregate measures or DNA fingerprinting techniques (Kent et al., 2004, we demonstrated the complex correlations between successions of bacterial and phytoplankton community assessed by 454 high-throughput pyrosequencing, which is particularly helpful to describe the microbial diversity at a fine scale based on large number of sequences (Fortunato et al., 2012).
The release of dissolved organic matter (DOM) by phytoplankton has long been recognized as an important source of high quality carbon to bacteria (Cole, 1982). Meanwhile, the DOM evolves in quality and quantity during phytoplankton community succession . Thus, changes in composition of phytoplankton community can drive changes in BCC through DOM production. It has been reported that phytoplankton plays a key role in regulating BCC in both mesocosms and natural ecosystems (H€ ofle et al., 1999;Pinhassi et al., 2004;Rooney-Varga et al., 2005;Tian et al., 2009). In the present study, the biomass of Cyanophyta increased rapidly from May to June 2010 and then was replaced by Bacillariophyta. The biomass of Bacillariophyta achieved dominant position from October 2010 to January 2011. Therefore, we can infer that the diverse DOM originated from Cyanophyta and Bacillariophyta might have markedly influenced the composition of bacterial community. Furthermore, we used DFA to discriminate the phytoplankton dominant period and nonphytoplankton dominant period. Our DFA results demonstrated that the functions can accurately predict the months of phytoplankton blooms, as Cyanophyta and Bacillariophyta dominated the calculation of the discriminant functions. More interestingly, the bloom of Cyanophyta and Bacillariophyta can accurately predict the distribution and abundance of the bacterial OTUs. Conversely, some bacterial OTUs were detected as indicator species of phytoplankton blooms in the Tingxi Reservoir. For example, 8 OTUs dominated the calculation of the discriminant functionsthey were affiliated with the divisions Actinobacteria, Proteobacteria, and unclassified bacteria. Additionally, MRA showed the DFA scores of bacterial community were predictable primarily from the environmental variables of pH, EC, chlorophyll a, biomass of Cyanophyta and Bacillariophyta, TP, and TN/TP. The biomass of Cyanophyta and Bacillariophyta and TP  . 1). Cyan -Cyanophyta, Baci -Bacillariophyta, Chl achlorophyll a, ECelectrical conductivity, TPtotal phosphorus, TN/TPtotal nitrogen/total phosphorus (mass ratio). The red rectangles represent phytoplankton species, while the blue ellipses represent bacterial OTUs (For full network see Fig. S1).
explained the highest proportion of the variance in diversity. Therefore, our results suggest that the bacterial communities are likely to be regulated by the succession of phytoplankton in the Tingxi Reservoir. It has been reported that correlations between succession of bacterial and phytoplankton communities are extremely complex in aquatic ecosystems, which include completion, symbioses or parasitism, etc. (Cole, 1982;Amin et al., 2012). In some mesocosms and natural experiments, specific bacterioplankton taxa that were in association with different phytoplankton taxonomic groups were also found (Pinhassi et al., 2004;Li et al., 2011). However, few studies attempt to test the statistical significance of correlations between bacteria and phytoplankton of whole communities at species level. In this study, we used eLSA based on the absolute quantity data to examine the complex correlations between bacterial and phytoplankton communities at species level. More importantly, we found that positive correlations dominated our association network (130 positive correlations vs. 14 negative correlations). Positive correlations may indicate microorganisms that thrive under similar envi-  ronmental conditions. On the other hand, positive correlations suggest the possibility of functional interdependencies between bacteria and phytoplankton. Further, phytoplankton are able to release pheromones to construct their own phycosphere , and these pheromones can only be perceived by specific bacteria, which form consistent associations and engage in specific interactions with phytoplankton. Phytoplankton can benefit from these interactions, as they require the bacterial production of vitamins. On the other hand, bacteria are able to acquire dissolved organic carbon from the autotrophic phytoplankton (Maurin et al., 1997;Kent et al., 2006;Sadro et al., 2011). Therefore, specific phytoplankton can interact with certain bacteria and thus affect lake bacterial community succession (Paver et al., 2013).
Recently, various studies focused on the bloom-forming Microcystis and found some specific bacterioplankton appear to strongly positively associate with Microcystis blooms (Niu et al., 2011;Cai et al., 2013;Parveen et al., 2013). For example, in our study, M. flosaquae was the dominated species from May to June. We also found six bacterial OTUs which were affiliated with the divisions Acidobacteria, Actinobacteria, and Proteobacteria were significantly correlated with the dynamics of M. flosaquae (Fig. 4). More interestingly, according to the BLAST analyses, these six OTUs might be closely related to Microcystis species. Of the twelve clones listed in  (Li et al., 2007;Wang et al., 2012). In contrast, our association network also showed some negative correlations between bacteria and phytoplankton. Bacteria can utilize inorganic nutrients in freshwater ecosystems. Through this progress, bacteria could have a profound influence on phytoplankton which growth depends on the supply of phosphorus. Currie & Kalff (1984) demonstrated that bacterioplankton are a superior competitor for phosphorus relative to the phytoplankton (algae). They found a freshwater bacterium (Pseudomonas paucimobilis), which utilized the excretion of phosphoruslimited alga with a carbon source, can outcompete that alga for phosphorus under widely varied phosphorus supply rates (Currie & Kalff, 1984). Besides phosphorus, heterotrophic bacterioplankton are more efficient at acquiring nitrogen from organic compounds than phytoplankton such as leucine (Løvdal et al., 2008). Moreover, specific bacteria are able to release sufficient molecules to kill phytoplankton (Mu et al., 2007;Amin et al., 2012). These algicidal bacteria have potential application as biocontrol agents of harmful algal blooms, because they increase in abundance concurrently with the decline of algal blooms (Mayali & Azam, 2004;Kodama et al., 2006).
We used a forward selection procedure to select physico-chemical variables that related significantly (P < 0.05) with the variation of BCC and phytoplankton community composition. In addition, we validated the causal relationship between BCC, significant physico-chemical parameters (water temperature, TP, and EC), and phytoplankton community succession using SEM. Bienhold et al. (2012) indicated that forward selection procedure and structural equation are complementary ecological modeling methodswith the former, the response and explanatory variables are defined a priori, while with the latter, different ecological correlations among variables are evaluated within a causal modeling framework. Therefore, we selected the significant variables based on the forward selection procedure and described the plausible interactions between the variables in the initial model. In this study, the SEM showed an indirect effect of water temperature on the bacterial community, acting through phytoplankton community succession. As one of the main seasonal factors, water temperature may influence abundance and composition of the phytoplankton community (Lv et al., 2014). For example, Cyanophyta may benefit from high temperature, as they have high optimum growth temperatures (Kosten et al., 2012). Hence, the differences of water temperature in Tingxi Reservoir could select different phytoplankton and thereby affecting the variation of BCC.
In conclusion, our results revealed a remarkable seasonality of bacterial community which was closely related to phytoplankton community succession. Furthermore, the association networks showed complex correlations between bacteria and phytoplankton communities at species level. Although, there is still a long way to go to show that those correlations are real interactions, our study provided an incomplete list of potentially interacting bacteria and phytoplankton that can be tested experimentally in future. The dominated positive correlations may suggest facilitative correlations between phytoplankton and bacteria are more important in the Tingxi Reservoir, while bacterial communities also have the potential negative correlations with phytoplankton community dynamics; for example, bacterioplankton compete with phytoplankton for inorganic nutrients, or bacteria kill phytoplankton through attachment and/or release of enzymes. Furthermore, we found phytoplankton succession explained the largest part of the temporal dynamics of the bacterial community. Our results provided new information on the microbial interactions as well as the possible mechanism that driving the variation of microbial plankton community composition in an aquatic ecosystem. Finally, additional long-term and extensive investigations from more freshwater ecosystems are needed to better understand the correlations of bacteria and phytoplankton communities.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Fig. S1. Association network showing all significant correlations (P < 0.01, Q < 0.004) between phytoplankton species and bacterial OTUs detected in the Tingxi Reservoir. Fig. S2. Initial structural equation model illustrating all plausible interaction pathways between bacterial community composition, significant relevant contextual physico-chemical parameters (water temperature, total phosphorus and electrical conductivity), and phytoplankton community. Fig. S3. Rarefaction curves of similarity-based operational taxonomic unit (OTUs) at cluster distance value of 0.03. Left -After identical subsampling, right -After removing Cyanobacteria, Archaea and chloroplasts (see Materials and methods). Table S1. Similarity-based OTUs (cluster distance values = 0.03) and species richness estimates (see Materials and methods).