ABSTRACT
Harmful algal blooms (HABs) are serious ecological disasters in coastal areas, significantly influencing biogeochemical cycles driven by bacteria. The shifts in microbial communities during HABs have been widely investigated, but the assembly mechanisms of microbial communities during HABs are poorly understood. Here, using 16S rRNA gene amplicon sequencing, we analyzed the microbial communities during an early-spring diatom bloom, in order to investigate the dynamics of microbial assembly processes. Rhodobacteraceae, Flavobacteriaceae, and Microbacteriaceae were the main bacterial families during the bloom. The 30 most abundant operational taxonomic units (OTUs) segregated into 4 clusters according to specific bloom stages, exhibiting clear successional patterns during the bloom process. The succession of microbial communities correlated with changes in the dynamics of algal species. Based on the β-nearest taxon distance, we constructed a simulation model, which demonstrated that the assembly of microbial communities shifted from strong heterogenous selection in the early stage of the bloom to stochasticity in the middle stage and then to strong homogeneous selection in the late and after-bloom stages. These successions were driven mainly by chlorophyll a contents, which were affected mainly by Skeletonema costatum. Moreover, functional prediction of microbial communities showed that microbial metabolic functions were significantly related to nitrogen metabolism. In summary, our results clearly suggested a dominant role of determinacy in microbial community assembly in HABs and will facilitate deeper understanding of the ecological processes shaping microbial communities during the algal bloom process.
IMPORTANCE Harmful algal blooms (HABs) significantly influence biogeochemical cycles driven by bacteria. The shifts in microbial communities during HABs have been studied intensively, but the assembly mechanisms of microbial communities during HABs are poorly understood, with limited investigation of the balance of deterministic and stochastic processes in shaping microbial communities in HABs. In this study, the dynamics and assembly of microbial communities in an early-spring diatom bloom process were investigated. Our data both confirm previously observed general microbial successional patterns and show new detailed mechanisms for microbial assembly in HABs. These results will facilitate deeper understanding of the ecological processes shaping microbial communities in HABs. In addition, predictions of metabolic potential in this study will facilitate understanding of the influence of HABs on nitrogen metabolism in marine environments.
INTRODUCTION
Phytoplankton blooms (e.g., diatoms, haptophytes, and autotrophic dinoflagellates) account for a significant fraction of marine primary productivity. Diatoms alone are estimated to contribute 20% to 40% of global carbon fixation (1). Harmful algal blooms (HABs) induced by diatoms in the open ocean or in coastal areas can be a serious problem for aquaculture and fisheries (2). Extensive efforts have been made to reveal the factors that trigger HABs, as it is crucial to understand their dynamics and thus control them. The occurrence and development of HABs have been explained by a complex combination of environmental conditions and microbial factors (1). Because of their roles in biogeochemical cycling during phytoplankton blooms, prokaryotes are now considered significant factors influencing the development of HABs (3).
Bacterial communities are considered a fundamental component of marine ecosystems, and they are a major driver of nutrient cycles in coastal waters. Moreover, bacteria have extensive interactions with phytoplankton (4, 5), and these interactions may influence microbial community compositions (MCCs). Bacterial clades, responding sensitively to phytoplankton blooms, belong to Gammaproteobacteria, Flavobacteriia (phylum Bacteroidetes), and the Roseobacter clade (6). This response is not so specific and includes distinct clades of taxa that bloom one after another. Study of the MCCs in mosaic phytoplankton blooms showed that Alphaproteobacteria, Gammaproteobacteria, and Flavobacteria represented the majority of the bacterial communities (7). SAR11 and the Roseobacter clade, within the class Alphaproteobacteria, are of particular importance for the turnover of dissolved organic matter (DOM) in the bloom (7). This trophic connection leads to synchronized blooms of bacteria during phytoplankton blooms (8, 9). Studies of the variation of MCCs during the HAB process have been widely reported, but few studies have revealed MCC assembly mechanisms during the bloom process.
Elucidating the processes regulating MCCs would make a significant contribution toward a comprehensive understanding of bacterial dynamics. It is well recognized that determinacy and stochasticity are two major types of assembly processes influencing MCCs (10, 11). The deterministic process involves selection from the abiotic environment and species interactions (12). In contrast, the stochastic process is induced by unpredictable disturbances, probabilistic dispersal of microbes, and random birth and death events (13). Although deterministic and stochastic processes have been found in animal and plant communities, studies of these two processes in microbial communities have been reported only in recent years (12, 14). Based on current findings, the deterministic process may either sort species adapted to a given environment or exclude species that may otherwise enhance ecosystem function (15). The balance of the stochastic process has a significant role in defining the assembly of microbial communities (16). Deterministic and stochastic processes have been widely discussed in relation to various ecological processes, but little is known about their roles in the assembly of microbial communities during natural HAB processes.
Xiangshan Bay, located in Zhejiang Province of China, is one of the most important aquaculture bases with poor water exchange ability. Thousands of fish cages, combined with the input of residual feed and fish excretion, lead to severe eutrophication in this semienclosed harbor (17). Consequently, HABs occur frequently, resulting in MCC disturbances (9). An algal bloom is a strong environmental filtering process that affects MCCs, but the assembly mechanisms of microbial communities during the bloom process need to be elucidated. In the present study, we investigated microbial communities during a diatom bloom process in Xiangshan Bay. We sought to elucidate (i) the composition and succession patterns of microbial communities during the HAB process and (ii) changes in the relative influences of stochastic and deterministic processes in the assembly of microbial communities during distinct bloom stages. We also evaluated the predicted metabolic functions of microbial communities during the HAB process, which few studies have reported previously. We think that the results of this study will facilitate understanding of the importance of bacterial roles in the fate of HABs and the mechanisms for the microbial assembly process during HABs.
RESULTS
Bloom characteristics.The sudden early-spring bloom in Xiangshan Bay was induced by Skeletonema costatum and Guinardia delicatula. Chlorophyll a (Chl a) contents reached 40 μg/liter in the early stage (ES) (Fig. 1B). Cell densities of S. costatum and G. delicatula ranged from 5.7 × 106 to 8.1 × 106 cells/liter in the ES; the dominant algal species changed from S. costatum to G. delicatula during this stage. After the ES, the abundance of S. costatum decreased dramatically, but the G. delicatula cell density reached its maximum in the middle stage (MS), contributing to the large contents of Chl a. In the late stage (LS), the abundance of S. costatum was 2.1 × 105 cells/liter, which was 30 times lower than that in the ES. The contents of Chl a were mainly accompanied by the occurrence of S. costatum. In the after-bloom (AB) stage, the concentration of Chl a decreased to 0.98 μg/liter and we could not precisely quantify S. costatum and G. delicatula; therefore, the abundance data are missing.
Map of Xiangshan Bay, indicating the sampling sites (A), and the abundance of algal species (B). ES1 and ES2 indicate samples from 10 February 2017 and 15 February, respectively, and MS1 and MS2 indicate samples from 19 February and 23 February. LS indicates samples from 27 February. AB indicates samples from 8 March. ES1 and ES2 were integrated into the ES of the bloom according to the Chl a contents, and MS1 and MS2 were integrated into the MS of the bloom. The map was created using Surfer 13 (Golden Software).
Microbial α-diversity.We observed clear differences in richness and phylogenetic diversity at distinct bloom stages, gradually decreasing from the ES to the AB stage (see Fig. S1 in the supplemental material). There were no remarkable changes in the Shannon index and evenness during the whole bloom process. Richness and phylogenetic diversity showed the strongest positive correlations with Chl a contents (r = 0.426 [P < 0.05] and r = 0.494 [P < 0.01], respectively), followed by dissolved oxygen (DO) (r = 0.386 [P < 0.05] and r = 0.430 [P < 0.05], respectively). Richness and phylogenetic diversity showed the strongest negative correlations with salinity (r = −0.491 [P < 0.01] and r = −0.531 [P < 0.01], respectively), ammonium (r = −0.467 [P < 0.05] and r = −0.475 [P < 0.01], respectively), and water temperature (r = −0.423 [P < 0.05] and r = −0.465 [P < 0.05], respectively) (Table S1). The environmental variables were not correlated with the Shannon index and evenness. These results indicated that richness and phylogenetic diversity were susceptible indices during algal bloom processes.
MCCs during the bloom process.The dominant bacterial classes were Alphaproteobacteria (relative abundance, 32.9%), Flavobacteriia (28.8%), Gammaproteobacteria (18.6%), Actinobacteria (8.2%), Betaproteobacteria (3.7%), and Verrucomicrobiae (1.9%), which constituted 94.1% of the total number of sequences (Fig. 2A). Furthermore, Rhodobacteraceae (23.8%), Flavobacteriaceae (20.4%), Microbacteriaceae (8.1%), Cryomorphaceae (8.1%), Surface1 (SAR11) (6.6%), Oceanospirillaceae (5.7%), SAR86 (3.9%), Verrucomicrobiaceae (1.8%), and Puniceicoccaceae (1.2%) were the main bacterial families, accounting for 79.6% of the total number of sequences (Fig. 2B).
Relative abundances of microbial lineages at the class (relative abundance of >0.5%) (A) and family (relative abundance of >5%) (B) levels during the diatom bloom process.
MCCs changed greatly during the bloom process (Fig. 2). Rhodobacteraceae was the most abundant taxon (average relative abundance, 37.1%) in the ES but decreased significantly in the AB stage (15.2% abundance). Compared with the bloom stages (ES, MS, and LS), the proportions of SAR11 (from 2.9% to 17.9%), Comamonadaceae (from 0.11% to 6.4%), SAR86 (from 1.7% to 7.3%), and Vibrionaceae (from 0.03% to 1.6%) increased significantly at the AB stage. Oceanospirillaceae and Verrucomicrobiaceae reached maximum abundance in the MS. Flavobacteriaceae had relatively stable abundance throughout the bloom process. In addition, the 30 most abundant operational taxonomic units (OTUs) were affiliated mainly with Rhodobacteraceae, Flavobacteriaceae, Oceanospirillaceae, Verrucomicrobiaceae, and SAR11 (Fig. 3). With only these OTUs involved in the random forest analysis, the overall out-of-bag (OOB) error was only 3.45% (Table S2), with 100% classification accuracy in the MS, the LS, and the AB stage, suggesting that the 30 most abundant OTUs could indicate the distinct algal bloom stages precisely. The cluster of columns in the heatmap (Fig. 3) revealed approximately 4 clusters according to the bloom stages, exhibiting clear successional patterns during the bloom process.
Abundance of the 30 most abundant OTUs in 4 bloom stages. Samples from the same stage clustered together. Microbial abundance was scaled with log transformation in the heatmap.
We also investigated the correlation of the 30 most abundant OTUs with the abundance of S. costatum. Fourteen of the 30 most abundant OTUs were strongly correlated with S. costatum (Fig. S2). Five OTUs showed significant (P < 0.001) positive correlations with S. costatum, i.e., OTU12736 (R2 = 0.827), OTU7254 (R2 = 0.696), OTU5946 (R2 = 0.678), OTU10900 (R2 = 0.646), and OTU12316 (R2 = 0.450). Negative linear correlations with S. costatum were observed for OTU2635 (R2 = 0.598 [P < 0.001]), OTU4140 (R2 = 0.418 [P < 0.001]), OTU3566 (R2 = 0.373 [P < 0.01]), OTU10614 (R2 = 0.306 [P < 0.01]), OTU8289 (R2 = 0.305 [P < 0.01]), and OTU6302 (R2 = 0.240 [P < 0.05]). Three OTUs (OTU1301, OTU5613, and OTU10568) showed negative nonlinear correlations with S. costatum. In summary, 9 of the 30 most abundant OTUs were closely associated with the decline of S. costatum, which indicated that those OTUs might inhibit the growth of S. costatum or promote its death during the algal bloom process.
Successional patterns of microbial communities during the bloom process.The nonmetric multidimensional scaling (NMDS) plot showed that samples clustered together according to distinct bloom stages (Fig. 4A), suggesting that there were variations in community structure over the bloom stages. Permutational multivariate analysis of variance (PERMANOVA) revealed that the structures of the microbial communities varied significantly (R2 = 0.527 [P < 0.0001]) among distinct bloom stages (Table 1), which indicated that the succession of microbial communities was induced mainly by the diatom bloom. Moreover, there were significant correlations between microbial communities and all environmental variables (Table 1; also see Fig. S3 in the supplemental material), as evidenced by Mantel tests and PERMANOVA (Table 1). In addition, the succession of microbial communities occurred mainly along the NMDS1 axis, which was significantly correlated with the Chl a content (R2 = 0.478 [P < 0.01]) (Fig. 4B). Therefore, even if this sudden algal bloom occurred during dramatic environmental changes, the dynamics of algal species might be a key factor driving the succession of microbial communities.
NMDS ordinations of microbial community dissimilarity in the bloom process (A), and correlation between Chl a contents and NMDS1 (B). Colors correspond to different sampling stages, indicating the succession of microbial communities together with NMDS1.
PERMANOVA (Adonis) analysis based on Bray-Curtis dissimilarity with 999 permutations and simple Mantel tests for correlations between environmental variables (Euclidean distance) and β-diversity of microbial communities (Bray-Curtis dissimilarity) with 999 permutations
Community assembly processes in successive stages.The β-nearest taxon index (βNTI) was used to provide deep insight into the potential roles of deterministic or stochastic factors in phylogenetic community succession of microbial communities. Based on the analysis of phylogenetic turnover, βNTI values were greater than 2 in the ES (Fig. 5A). After that, βNTI values decreased significantly throughout the bloom process. In the MS, most βNTI values ranged from −2 to 2, suggesting a role of stochasticity. In the LS and the AB stage, almost all of the βNTI values were lower than −2, indicating a deterministic process. These results revealed that βNTI distributions progressively shifted along the bloom stages, and the mechanisms driving the succession of microbial communities might have been different.
Patterns of distribution of βNTI values across bloom stages (A), and effect of algal biomass (indicated by Chl a contents) on the distribution of βNTI values (B). In panel A, horizontal dashed lines indicate lower and upper significance thresholds at −2 and 2, respectively; the asterisks denote the βNTI average in each group.
Based on the pattern of βNTI values and on environmental heterogeneity in different bloom stages, we assembled microbial communities at 4 bloom stages in a simulation model (Fig. 6). Environmental heterogeneity was the highest in the ES, indicating that selection by environmental variables ruled the deterministic process. In the MS, there was a transition from selection of heterogenous variables to stochasticity, as environmental heterogeneity decreased significantly. In the LS and the AB stage, environmental heterogeneity continuously decreased, and the deterministic assembly process of the microbial communities was induced mainly by homogeneous selection. The distribution of βNTI values changed in the simulation, corresponding closely to the shifts in Chl a contents (Mantel ρ = 0.264 [P = 0.004]) (Table 2 and Fig. 5B) and silicate levels (Mantel ρ = 0.202 [P = 0.003]). The simulation showed that assembly of microbial communities shifted from strong heterogenous selection (βNTI of >2, in the ES) to stochasticity (∣βNTI∣ of <2, in the MS) and then to strong homogeneous selection (βNTI of less than −2, in the LS and the AB stage), accompanied mainly by the dynamics of algal species.
Simulation model of the successional pattern of microbial communities based on βNTI distributions. Ecological selection is weak from −2 to 2 (green area), while selection is strong toward positive and negative values on the y axis, resulting from strong variable selection (VS) and homogeneous selection (HS), respectively. The boxplot shows the environmental heterogeneity in each bloom stage. The bar plot shows the abundance of 2 harmful algal species. The prebloom stage, without bloom disturbance, may have ecological selection similar to that of the AB stage.
Simple Mantel tests for correlations between environmental variables (Euclidean distance) and βNTI distance matrix for microbial communities with 999 permutations
Predicted metabolic potential during the bloom process.The nearest sequenced taxon index (NSTI) values ranged from 0.056 to 0.113; the mean value was 0.090 ± 0.013 for all samples (see Data Set S2 in the supplemental material). The Bray-Curtis distance-based redundancy analysis (db-RDA) showed a clear shift in the composition of functional gene families throughout the 4 bloom stages (Fig. S4). The most frequent predicted functions were amino acid metabolism (12.2%), membrane transport (11.4%), carbohydrate metabolism (10.4%), replication and repair (7.1%), and energy metabolism (5.8%) (Fig. S5).
Nutrient levels, including nitrate, nitrite, ammonium, phosphate, and silicate levels, were positively correlated with predicted metabolic potential in the MS and LS, indicating greater relative potential of nitrogen metabolic dynamics than in the other stages (Fig. S4 and Table S3). Therefore, gene families of nitrogen metabolism were investigated in detail (Fig. 7). Seventeen genes segregated into 4 clusters according to their predicted abundance (Fig. 7). Cluster I included mainly NO/N2O-reducing-related genes (norB, norC, and nosZ). Cluster II included mainly nitrogen-fixation-related genes (nifK, nifH, and nifD), which were particularly frequent in most samples from the ES. Cluster III genes, including nitrite-reducing-related genes (nirB and nirD), had greater relative abundance than genes in other clusters. Genes in cluster IV (napB, narG_narZ_nxrA, narl_narV, and narH_narY_nxrB) were related to nitrate reduction.
Relative abundances of nitrogen metabolism-related genes during the bloom process, as predicted by PICRUSt. The predicted genes could be segregated into 4 clusters in a row according to their functions.
Regarding the contributions of bacterial families to the selected genes (Fig. S6), genes in cluster I were predicted mainly to occur in Flavobacteriaceae (>43.9% of contributions), which contributed only to this cluster. The high prevalence of cluster II genes in the ES was attributed mainly to Puniceicoccaceae (Verrucomicrobia; >49.1% of contributions). The predicted source of cluster III genes was primarily the Rhodobacteraceae family (>54.4% of contributions). The predicted sources of cluster IV genes were Comamonadaceae (>49.8% of contributions) and Halomonadaceae (47.8% of contributions in napB). The taxonomic metagenomic contributions identified here revealed that each of the gene clusters was contributed mainly by a single bacterial family.
DISCUSSION
Microbial α-diversity variations during the bloom process.Microbial α-diversity decreased significantly during the diatom bloom process (see Fig. S1 in the supplemental material), as observed previously in Xiangshan Bay (9). Factors correlated with bacterial richness and growth would be expected to determine the microbial succession. Chl a is one of those factors, which can be used to estimate phytoplankton biomass directly (18). The spatial and temporal distributions of Chl a are closely related to marine bacterial abundance (19, 20). Moreover, a microcosm study found that the presence of phytoplankton could affect bacterial richness and production significantly (21). In line with these studies, we found that microbial α-diversity was strongly influenced by diatom biomass (Table S1). Here we also showed that the dominant algal species shifted during the bloom process (Fig. 1B), which may also affect microbial α-diversity, because distinct algal species could release DOM, facilitating differential bacterial growth (22).
We should note that environmental variables, including salinity and water temperature, were related to microbial α-diversity. During this early-spring bloom, the water temperature increased significantly from the ES to the AB stage, which would evidently affect bacterial abundance (19). Salinity is generally associated with dramatic changes in bacterial abundance and community composition (23, 24). Consequently, microbial α-diversity was influenced mainly by bloom disturbances, together with water temperature and salinity.
MCCs associated with the diatom bloom.Members of the Alphaproteobacteria/Gammaproteobacteria and Flavobacteria constituted the major parts of MCCs in this early-spring bloom (Fig. 2A). Flavobacteriaceae and Rhodobacteraceae were the most dominant bacteria during the bloom process (Fig. 2B). The relative abundances of these taxa were similar to those in previously reported HABs in the Xiangshan Bay (9) and in other locations (21, 25). Flavobacteriaceae, which can utilize a broad range of biopolymers (particularly polysaccharides and proteins) (7), are known for their ability to thrive under bloom conditions in temperate waters, as well as in polar systems (7, 25). The relative abundance of Flavobacteria in the present study was stable, possibly as a result of eutrophication and input of organic matter from aquaculture in Xiangshan. Marine Rhodobacteraceae taxa are major vitamin suppliers for B12-auxotrophic prokaryotes and for eukaryotic primary producers (e.g., diatoms, dinoflagellates, and brown algae) (26, 27). Interestingly, about 90% of Rhodobacterales species possess the complete vitamin B12 synthesis pathways (28), indicating that Rhodobacterales might directly support the diatom bloom by supplying vitamins. Hence, a decrease in Rhodobacteraceae abundance during the bloom process might promote the decline in diatoms.
The bacterial taxa SAR11 and SAR86 are oligotrophic bacterial clades whose abundance may increase with diatom biomass and with Chl a contents (29, 30). Increases in Vibrio abundance have been associated with algal blooms and, in some cases, with specific species of algae (31). High Vibrio abundance may occur in diatom-dominated phytoplankton assemblages (32). Moreover, Vibrio abundance was related to water temperature (31). Therefore, diatom blooms and water temperature may contribute to the progression of Vibrio. Oceanospirillaceae was tightly linked to DOM (33), so that DOM derived from algae would significantly influence their growth. Verrucomicrobiaceae were also enriched in this bloom, because members of this taxon could use diverse glucoside hydrolases to degrade polysaccharides during the bloom process (34).
According to an analysis of the 30 most abundant OTUs, the abundance of Rhodobacteraceae was closely linked to that of S. costatum, or the decrease in Rhodobacteraceae was unfavorable for the growth of S. costatum. Nine OTUs with negative correlations with S. costatum belong to Oceanospirillaceae and Flavobacteriaceae. Members of Oceanospirillaceae and Flavobacteriaceae have long been reported to have algicidal activity (35, 36), which may inhibit the growth of S. costatum. Therefore, our results indicate that some dominant bacterial members may change the course of diatom blooms, according to the interactions between bacteria and diatoms.
Linking of microbial community variations with diatom blooms.Algal blooms could drive the variations in MCCs significantly, because DOM, allelochemical matter, and other growth factors derived from algal species could strongly influence bacteria growing in alga-related circumstances. Our results revealed that samples could cluster into the corresponding bloom stages (Fig. 3 and 4A; also see Fig. S3), and the diatom bloom could explain the largest MCC variations (Table 1), indicating that MCC variations were closely associated with bloom dynamics. Among the environmental variables analyzed, DO and pH were major abiotic factors governing the dynamics of microbial community structures (Table 1; also see Fig. S3). DO and pH are closely linked to algal growth, due to photosynthesis and the respiration of algal cells. Previous studies proved that DO and pH are major factors affecting microbial community structures in various ecosystems (37, 38). Taken together, our results demonstrated that diatom bloom dynamics would significantly determine the MCCs in Xiangshan Bay.
MCCs varied dramatically over the whole process, which supported the idea of strong temporal shifts in microbial assemblages, and clustered according to bloom stages, which suggested a potential repeatable bloom succession. The cluster of the 30 most abundant OTUs and random forest analysis revealed that MCCs shifted precisely with bloom dynamics (Fig. 3; also see Table S2). Moreover, random forest classification achieved a high level of prediction accuracy with the 30 most abundant OTUs (Table S2), suggesting the great bloom-indicating ability of these OTUs. This result has evident ecological significance for improving the accuracy of the forecast model for HABs. There are many heuristic models capable of simulating various aspects of HABs (39), but most of them have a relatively low level of forecast accuracy. The input parameters in these models are mainly physiochemical variables and not biological factors (e.g., bacteria or algae) that are inherent factors in HABs (2). Therefore, bacteria with bloom-indicating abilities should be integrated into forecast models to improve their accuracy in the future.
Relative roles of deterministic versus stochastic processes.Variations in the importance of deterministic and stochastic assembly processes are prevalent in natural ecosystems (12, 14, 40). To the best of our knowledge, the present study is somewhat novel in estimating the relative strengths of deterministic and stochastic processes in defining microbial communities across whole HABs. Our results revealed that diatom biomass was the most important factor driving microbial assembly throughout the HAB (Table 2 and Fig. 5B), but there might have been other factors driving microbial assembly in specific bloom stages.
According to βNTI values (Fig. 5A) and the simulation model (Fig. 6), we observed a greater influence of determinacy in the initial successional stages. This result is consistent with previous studies that reported microbial primary successions in various ecosystems (14, 16, 40). The high level of determinacy in the ES may be the result of eutrophication. The nutrient loads from catchments and excessive aquaculture activities exacerbate the degree of eutrophication in Xiangshan Bay (41). Consequently, strong selective forces imposed by high nutrient inputs may drive succession of the microbial communities in the initial or prebloom stages. In the MS, stochasticity mainly drove bacterial assembly. We inferred that physiochemical conditions are not as extreme as in the ES, because environmental heterogeneity decreased significantly. In contrast to conditions in other bloom stages, there was a strong wind on the MS sampling day, which would promote the immigration of bacteria by strong waves. Thus, random dispersal could have been an important ecological factor in the MS. In the LS and the AB stage, the deterministic process again governed bacterial assembly, but it was induced mainly by environmental homogeneity (Fig. 6). In the LS, many of the S. costatum organisms were lysed, and DOM derived from this algal species was released into the seawater, leading to strong homogeneous selection (5). In the AB stage, we observed an evident increase in temperature, which might also be an important factor driving MCCs, as reported previously (42). We inferred that the influence of the rapidly rising temperatures might override the effects of other environmental variables, inducing homogeneous selection at the AB stage.
Although our analyses explained the MCC successional patterns in the whole HAB, the repeatability in Xiangshan Bay and other algal bloom processes still needs to be tested. Nonetheless, the present data helped to elucidate the relative magnitude of forces driving MCCs, which are crucial components of both natural and agricultural ecosystems. In the future, more samples will be integrated into this model to test its repeatability and to improve its accuracy.
Predicted metabolic potential in the diatom bloom.Evidence for variation in MCCs in HABs has been widely reported, but few studies focused on the metabolic potential of the microbial communities. To gain more insights into the metabolic potential of bacteria in the HAB process, PICRUSt was used to predict the metagenome of each bloom stage. The NSTI values (see Data Set S2 in the supplemental material) were either similar to or lower than the mean NSTI value reported for bacterial communities in samples collected from various environments (43–45). Therefore, we inferred that our predicted metabolic potential is supported by the microbial genomic reference database.
The present study revealed that the metabolic functions of microbial communities were significantly related to nitrogen (Table S3 and Fig. S4), and distinct bacterial families contributed to the different metabolic potentials (Fig. 7; also see Fig. S6), which may be induced by the nitrogen nutrient levels. Cluster I consisted mainly of genes encoding nitric oxide and nitrous oxide reductases. This metabolic potential might be due to the family Flavobacteriaceae (46). The metabolic potential in cluster II was related to nitrogen fixation, with Puniceicoccaceae (phylum Verrucomicrobia) being the main contributor. A novel bacterium of the family Puniceicoccaceae, with nitrogen-fixing ability, was found last year (47). In addition, nitrogen fixation by heterotrophic bacterial communities could be significantly stimulated by high nutrient loads (48); therefore, there may be a relatively greater nitrogen fixation potential in the ES, in comparison to other stages (Fig. 7). Rhodobacteraceae, the main clade in cluster III, is an important denitrifying family (49) related to nitrite reduction. Members of Comamonadaceae and Halomonadaceae in cluster IV have been observed to catalyze the reduction of nitrate (50, 51). These results indicated that nitrogen metabolism is mediated by a diverse group of bacteria, even though some of them are not the most abundant taxa. Moreover, anthropogenic and agricultural activities resulted in large nitrate amounts entering Xiangshan Bay, and HABs may significantly influence the removal of nitrate from the ecosystem.
In the present study, we studied the microbial metabolic function potential in the HAB process, based on 16S rRNA gene sequencing and PICRUSt, and we realized that there might be some bias in this method. Gene ontology in the PICRUSt script is based on the GreenGenes v13.5 database, which is not as accurate as the newest SILVA database. Moreover, the annotation is implemented based on well-characterized human-pathogenic species, which may not be suitable for environmental samples. A more accurate input metagenome would require a newer database and more annotated species from environmental samples. In addition, quantitative data, such as absolute bacterial abundance, and functional genes should be integrated. Therefore, more reliable analytical techniques, including quantitative PCR, metagenomic sequencing, and metatranscriptomic sequencing, will be applied to study metabolic function potential in future work.
Overall, microbial communities in this early-spring bloom exhibited changes in diversity and composition during the bloom process. In particular, the 30 most abundant OTUs responded significantly to the algal bloom, exhibiting clear successional patterns in distinct stages. Moreover, our exploration of the diatom bloom in the Xiangshan Bay also provided novel insights about the relative roles of the deterministic and stochastic processes in driving microbial community structures. Both ecological processes played important roles in the assembly of microbial communities, but the deterministic process was of greater importance during the bloom process. In addition, metabolic potential prediction suggested that the outbreak of diatoms might have a significant impact on the denitrification process. In the future, detailed information about the interaction between microbial communities and phytoplankton, as well as metabolic functional characteristics, is essential to better understand the mechanisms influencing the formation of HABs.
MATERIALS AND METHODS
Sample collection.Samples were collected in February and March 2017 from 5 stations in Xiangshan Bay during the diatom bloom process, in the ES, MS, LS, and AB stage (Fig. 1). For the ES, samples were collected on 10 February and 15 February. For the MS, samples were collected on 19 February and 23 February. For the LS and the AB stage, samples were collected on 27 February and 8 March, respectively. Twenty-nine water samples were collected from the surface layer (0.5 m deep), using a 5-liter water sampler. For phytoplankton quantification, 1 liter of each water sample was fixed with 1% Lugol's solution. After filtration with a 100-μm mesh, approximately 600 ml of seawater was filtered using 0.2-μm filters (47-mm-diameter polycarbonate; Millipore, USA); the filters were then stored at −80°C prior to DNA extraction. To prevent contamination between samples, the water sampler and filtration systems were washed carefully with sterile water before water collection and filtration.
After homogenization by gentle shaking, phytoplankton samples were characterized and enumerated using a counting chamber with a light microscope (Axioplan; Carl Zeiss, Jena, Germany). Water temperature, salinity, and pH were measured on board. DO, Chl a, chemical oxygen demand (COD), phosphate, silicate, nitrate, nitrite, and ammonium levels were measured using standard methods (52).
DNA extraction and PCR amplification.DNA extraction was performed using the Power Soil DNA isolation kit (Mo Bio, USA), according to the manufacturer's instructions. DNA concentrations were quantified using a Nanodrop 2000 instrument (Thermo Fisher Scientific, Wilmington, DE, USA). The V4 hypervariable region of the 16S rRNA gene was amplified in triplicate PCRs (20-μl reaction volume, with 5× FastPfu buffer, 2.5 mM deoxynucleoside triphosphates [dNTPs], 5 μM combined primers, and FastPfu polymerase), using 10 ng DNA template and the primer pair 515f and 806r (53). PCR cycling conditions were as follows: denaturation at 95°C for 3 min, 27 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 45 s, and a final extension at 72°C for 10 min. Triplicate PCR amplicons were mixed after purification using a TaKaRa purification kit (TaKaRa, Japan). Libraries were then generated using the TruSeq DNA sample preparation kit for Illumina (Illumina, San Diego, CA, USA), following the manufacturer's instructions. The libraries were sequenced using a MiSeq platform (Illumina) in a paired-end 250-bp sequence read run at MajorBio Co. Ltd. (Shanghai, China).
Data processing.QIIME v1.9 (54) and USEARCH v6.1 (55) pipelines were used to perform quality control and to cluster sequences into OTUs. Overlapping reads were merged using FLASH, with default parameters (56). Raw sequences were demultiplexed and filtered to remove low-quality sequences using QIIME. Chimeras were detected and removed with USEARCH using the SILVA 128 database (57). OTUs were clustered at 97% sequence similarity using UCLUST (58). Taxonomic assignments were achieved by UCLUST with the SILVA 128 reference database. OTUs affiliated with chloroplast, mitochondrion, unassigned, and unclassified sequences and singletons were removed from the data set before downstream analysis. The numbers of final clean reads in each sample ranged from 16,618 to 82,467 reads, yielding a total of 3,529 OTUs (see Data Set S1 in the supplemental material).
Statistical analyses.For estimation of α-diversity, sequences in each sample were normalized at the lowest sequence depth and rarefied at 16,610 reads per sample. The different α-diversity estimations (richness, Shannon index, phylogenetic diversity, and evenness) were calculated using QIIME and visualized using the ggplot2 package (59). The relationships between α-diversity indices and environmental variables were calculated by Spearman's rank correlation using the psych package (60). The MCCs were analyzed at the class (relative abundance of >0.5%) and family (relative abundance of >5%) levels. The 30 most abundant OTUs were selected for study of the shifts in MCCs at distinct bloom stages, and the abundance of each OTU was scaled by color in a heat map using the pheatmap package (61). Random forest analysis was performed to test whether the 30 most abundant OTUs could predict the distinct bloom stages using the randomForest package (62). To track the successional trajectories of microbial communities in the bloom process, β-diversity was calculated using the Bray-Curtis distance and visualized by a two-dimensional NMDS plot. The db-RDA (63) was calculated using the capscale function in vegan (64), to explore the effects of environmental variables on MCCs. Mantel tests (65) and PERMANOVA (66) were used to verify the explanation of environmental variables in the succession of MCCs. To determine environmental heterogeneity, the homogeneity of multivariate dispersions, applied to the environmental variables (67), was calculated using the betadisper function in vegan. All of the environmental variables were standardized before the calculation of environmental heterogeneity. All statistical analyses were performed with R statistical software (https://www.r-project.org).
Analysis of assembly processes of microbial communities.The phylogenetic turnover (phylogenetic β-diversity) was quantified using the βNTI (12). Based on the random shuffling of OTU labels across the tips of the phylogeny, the βNTI was computed as the number of standard deviations that observed β-mean nearest taxon distance (βMNTD) values deviated from the mean of the null distribution (999 null replicates) (12, 68). The absolute magnitude of the βNTI reflects the influence of the deterministic process (∣βNTI∣ of >2); the greater the magnitude is, the greater the influence of the deterministic process is. Values of ∣βNTI∣ of <2 indicate that the stochastic process may play an important role in structuring microbial communities during the bloom process.
Metagenomic prediction of microbial communities.Functional metagenomes during the bloom stages were predicted from clean reads from 16S rRNA sequencing using PICRUSt (43). Briefly, OTUs were picked based on the closed reference OTU-picking method, at 97% sequence similarity, against the GreenGenes v13.5 database in QIIME. The OTU table was submitted to local PICRUSt v1.1.2 for metagenome prediction, which then created the final metagenomic functional predictions based on KEGG pathways, at classification levels 2 and 3. To ensure the accuracy of metagenomic predictions, the NSTI value, which indicates the distance of the 16S rRNA-sequenced microbes from the genome-sequenced microbes, was calculated. Lower NSTI values reflect a shorter distance between our sequenced microbes and the genome-sequenced microbes in the reference database (43).
A db-RDA was used to estimate the relationships between metabolic gene ontologies (Bray-Curtis distances from predicted gene abundance) and environmental variables, including nitrate, nitrite, ammonium, phosphate, and silicate levels. The genes predicted to be associated with nitrogen metabolism were clustered (genes with total abundance of <7,000 copies in all samples were discarded) and visualized using the pheatmap package in R (61). The contribution of OTUs to nitrogen metabolism-related genes was computed using metagenome_contributions.py in PICRUSt. Bacterial families with contributions to nitrogen metabolism of >1% were selected and graphed.
Accession number(s).Raw sequences were submitted to the NCBI Sequence Read Archive with the accession number SRP130859.
ACKNOWLEDGMENTS
This work was funded by the National Natural Science Foundation of China (grant 41706132), the Zhejiang Provincial Natural Science Foundation of China (grant LQ17D060001), the Natural Science Foundation of Ningbo (grant 2016A610094), the Open Fund of Key Laboratory of Integrated Marine Monitoring and Applied Technologies for Harmful Algal Blooms, State Oceanic Administration (grant MATHAB201708), and the K. C. Wong Magna Fund of Ningbo University.
FOOTNOTES
- Received 27 April 2018.
- Accepted 9 July 2018.
- Accepted manuscript posted online 13 July 2018.
- Address correspondence to Demin Zhang, zhangdemin{at}nbu.edu.cn, or Xiangyu Zhu, xiangshg{at}126.com.
Citation Zhang H, Wang K, Shen L, Chen H, Hou F, Zhou X, Zhang D, Zhu X. 2018. Microbial community dynamics and assembly follow trajectories of an early-spring diatom bloom in a semienclosed bay. Appl Environ Microbiol 84:e01000-18. https://doi.org/10.1128/AEM.01000-18.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.01000-18.
REFERENCES
- Copyright © 2018 American Society for Microbiology.