The ecological responses of bacterioplankton during a Phaeocystis globosa bloom in Beibu Gulf, China highlighted by integrated metagenomics and metatranscriptomics

The interaction between bacteria and phytoplankton during bloom events is complicated throughout the developmental processes of algal blooms. The detailed ecological roles of bacterioplankton during algal blooms still need to be investigated comprehensively. With the assistance of omics techniques, the composition and function of bacterioplankton were studied during the blooming and recession periods of Phaeocystis globosa in the Beibu Gulf, China. The transcriptionally active taxa of Vibrio, which correlated with most functional genes, were enriched in the blooming period, whereas the active microbial groups, such as Erythrobacter and Candidatus puniceispirillum, increased significantly in the recession period of the P. globosa bloom. The transcriptional analyses indicated the blooming and recession periods of P. globosa had different impacts on the function of bacterioplankton, including shifts in expression of specific metabolic pathways for variable nutritional utilization and their advantage in bacterial motility, quorum sensing, and extracellular secretion. Overall, the integrated field investigation and in-depth analysis of molecular data highlighted the difference of bacterial community during the outbreak and collapse stages of P. globosa bloom, and provided fundamental data for better understanding of the bacterial community in response to blooming and recession of the P. globosa bloom.


Introduction
Marine bacterial community (including bacteria and archaea) play fundamental roles in global biogeochemical cycles by mediating organic matter mineralization (DeLong and Karl 2005;Crespo et al. 2013). A large proportion (10-50%) of organic matter produced by phytoplankton are remineralized by marine bacteria (Cole et al. 1988), returning massive amounts of organically bound nutrients for reuptake by phytoplankton (Castberg et al. 2001;Azam and Malfatti 2007;Egan et al. 2013). It is, however, increasingly recognized that the interaction between bacteria and phytoplankton during bloom events vary throughout the developmental processes of algal blooms (Buchan et al. 2014;Chen et al. 2015). Bacteria may pose stimulative or inhibitory effects on the growth of algal cells by secreting soluble extracellular substances (Guo et al. 2016;Tan et al. 2016;van Tol et al. 2017), or even promote algal aggregation and sedimentation in seawaters (Sonnenschein et al. 2012;Powell and Hill 2013). The detailed ecological roles of Responsible Editor: Sandra Shumway.
* Caiwen Li cwli@qdio.ac.cn 1 bacterioplankton during algal blooms still need to be investigated further. Progress in next-generation sequencing is fueling a rapid increase in the number and scope of bacterial communitytargeted studies (Gilbert et al. 2008;Vila-Costa et al. 2010;Rinta-Kanto et al. 2012). Meanwhile, metagenomics and single-cell amplified genomes have provided insight into the taxonomic composition and metabolic potential of a complex bacterial community (DeLong et al. 2006;Swan et al. 2011;Rinke et al. 2013;Yooseph et al. 2013). Lately, metatranscriptomics have served to unveil the actual gene expressions and metabolic activities of the whole community at a specific time and place, and how those activities vary in response to environmental fluctuation or biotic interaction . It has been demonstrated that comparative metatranscriptomics were highly efficient in deciphering the phylogenetic composition, metabolic potential, and pathways of marine bacterial community (Gifford et al. 2013;Ottesen et al. 2013;Wu et al. 2013;Penn et al. 2014). Previous metatranscriptomic studies indicated differences in gene expressions during the peak of a marine phytoplankton bloom compared to the post-bloom period in microcosms (Gilbert et al. 2008). The metabolic activities of bacterial families, including Flavobacteria (genera Ulvibacter, Formosa, and Polaribacter), Alphaproteobacteria (SAR11 clade and Rhodobacteraceae) and Gammaproteobacteria (genus Reinekea and SAR92 clade), were higher during bloom events (Klindworth et al. 2014;Wemheuer et al. 2014Wemheuer et al. , 2015. The genes involved in cell surface adhesion and organic acid utilization (Rinta-Kanto et al. 2012), toxic microcystin (Penn et al. 2014), and transposases (Wemheuer et al. 2015) were elevated and differentially expressed during blooms. Those previous studies provide useful insight into the in situ metabolic activities and functional partition of microbes in the bloom-associated bacterial community.
The marine prymnesiophyte Phaeocystis globosa, has a complex life cycle (usually exists as a colony or a freeliving cell), was globally distributed and frequently caused large-scale algal bloom . Since 2011, blooms of P. globosa have occurred annually in the Beibu Gulf, Guangxi Province, China, and have resulted in significant negative impacts to local industry and aquaculture (Qin et al. 2016). Series of investigations have been conducted to explore the succession of marine plankton, as well as the correlations among abundance of key function groups and nutrients Qin et al. 2018); however, limited effort has been initiated to elucidate the detailed ecological responses of the bacterial community during P. globosa blooms by modern comparative metagenomics or metatranscriptomics. Thus, in the present study, coupled metagenomic and metatranscriptomic sequencing derived from environmental DNA and RNA were applied to elucidate the compositional and functional responses of bacterial community during a Phaeocystis globosa bloom event in a typical tropical coastal bay. The goals were to address the following questions: How did the bacterial community change in response to blooming and recession of the P. globosa bloom? What functional processes and specific gene transcriptional behaviors of the bacterial community were influenced by the algal bloom? And, how did the dominant metabolic processes respond to the fluctuating nutrient environment? The result will contribute to better understanding of the correlations among microbial assemblages and functional attributes, and further highlight the ecological roles of bacterioplankton across 'boom' and 'bust' period of harmful algal blooms.

Site description and sampling
Water samples were collected from three stations (P1, P4, and P7) during the different stages of a bloom of Phaeocystis globosa on January 15th and February 2nd 2018, in the north Beibu Gulf, Guangxi Province, China (Supporting Information (SI) Fig. S1 and Table 1). Nine water samples (triplicate for the three stations, identified as bP1, bP4, and bP7) were taken during the bloom period (bP), and the other nine (triplicate for the three stations, identified as rP1, rP4, and rP7) were taken during the recession period (rP) of the bloom event. For bacterial DNA/RNA extraction, samples (20-L per treatment) were prefiltered through a 20-μm sieve and filtered through a 0.22-μm pore size polycarbonate filter (142-mm diameter, Millipore, Germany) successively using a peristaltic pump. Samples of 2-mL were filtered through a 20-μm sieve and fixed with 10% sterile particlefree paraformaldehyde (1% final concentration, v v −1 ) for bacterial enumeration using flow cytometer (Accuri C6, Becton, Dickinson and Company, the United States). Physicochemical characteristics (temperature (T), salinity (S), chlorophyll a (Chl a), dissolved oxygen (DO), total nitrogen (TN), ammonium salt (NH 4 + ), nitrate plus nitrite (NO 2 − , NO 3 − ), soluble phosphate (PO 4 3− ), silicate (SiO 3 2− ), bacterial abundance, and colony abundance) were measured either on site or in the laboratory for each sample. Detailed protocols are provided in SI.

Processing and analysis of environmental metagenomic samples
For the water samples, total environmental DNA (triplicate for each sample, 6*3 = 18 samples in total) was extracted, processed, and analyzed by PCR amplification of the V3-V4 region of the 16S rRNA genes. The protocols are detailed in the SI. Briefly, DNA was extracted from the filter sandwich and prepared for sequencing using a series of kits following manufacture instructions. We ultimately obtained an average of 85,591 sequence reads per sample, and 80,135 high-quality sequences were obtained by quality control, with a quality control efficiency of 93.65% (Table S1). The raw reads have been deposited in the Sequence Read Archive (NCBI) under accession number PRJNA678802. Sequence analyses were performed by Uparse software (Uparse v7.0.1001; Quast et al. 2013). A total of 1262 OTUs were obtained, with sequences ≥ 97% similarity assigned to the same OTUs. A representative sequence for each OTU was screened and classified using the Silva Database (Version 132) for further annotation (Edgar et al. 2011). QIIME (Version1.7.0) was used to rarefy the OTUs table and normalize the read counts, in which the lowest sequencing depth obtained from a sample required to normalize the differences. Then, a cumulative-sum scaling (CSS) method was used to calculate the scaling normalization factors equal to the sum of counts up to a particular quantile (Wemheuer et al. 2013). Subsequent analyses of alpha diversity and beta diversity were performed with the normalized OTU abundance information.

Processing and analysis of environmental metatranscriptomic samples
Environmental total RNA (triplicate samples from each station were pooled together to acquire higher amount of environmental RNA, thus six in total) was extracted (equivalent to 60-L seawater processed), purified, and prepared for Illumina high-throughput sequencing as described in detail in SI. Sequencing libraries used were NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, USA) generated according to the manufacturer recommendations, and index codes were added to attribute sequences to each sample. With RNA-seq, after quality filtering by removing reads containing adapter, reads containing ploy-N, and low-quality reads from raw data, an average of 35,568,957 paired reads per sample were obtained (Table S2). The raw data have been deposited in the NCBI under accession number PRJNA686139. Ribosomal RNA reads were filtered out using SortMeRNA (version 2.1, 2/2/2016). The remaining non-ribosomal sequences (an average of 30,443,470) were assembled de novo with the stitching software Trinity (version: r20140413p1). The sequences were integrated, and redundancy eliminated using CD-HIT-EST with the sequence identity threshold of 0.95, and obtained 1,585,720 unigenes for further analyses (Table S2). The cDNA sequences were identified by BLAST against the NCBI nonredundant protein database (Nr, Version: 2016-11-05) using DIAMOND software (Conesa et al. 2016) with pair (blastp, evalue ≤ 1e −5 ). Bacterial and archaeal protein-coding gene sequences were extracted from Nr (the LCA algorithm), and the system classification of the MEGAN software (Huson  (Table S2).

Statistical methodology
All statistical analyses were conducted in the R environment (v.3.6.1; http:// www.r-proje ct. org/) unless otherwise indicated. Alpha diversity and beta diversity based on Bray-Curtis distance were calculated with multiple indices (observed species, Chao1, and Shannon-Wiener index) and the Bray-Curtis distance between samples. Nonmetric multidimensional scaling (NMDS) based on Bray-Curtis distance was performed using the "vegan" package (v.2.5-6) to investigate the general distributions of bacterial community composition with feature significance confirmed by permutational multivariate analysis of variance (PERMANOVA; Oksanen et al. 2016). Linear least-squares regression analysis was used to model relationships between the significant environmental conditions (colony abundance of P. globosa) and taxonomic features (alpha diversity variables) with Spearman's correlation test.
To compare the differences in taxonomic and functional composition of bacterial taxa, taxonomic ranks (at genus level) combined with their functions annotated with KEGG databases were comprehensively analyzed (Carrion et al. 2019). To determine the taxonomic composition of functions-associated taxa, custom R commands were used to aggregate all the differential expressed genes and the correlated taxonomic information in the Nr database, and a metatranscriptomic matrix based on the read abundance per functional taxon and per sample was generated. The relative abundance of bacterial taxa (at DNA, RNA, and the functional genus level) was used to construct chord diagram using "circlize" package (v.0.4.6; Gu et al. 2014) and heatmap with hierarchical clustering performed by "pheatmap" package (v.1.0.12). Two-tailed t-tests was performed by the STAMP (v.2.1.3; Parks et al. 2014) to assess whether the abundance of specific microbial taxa was different between groups. Obtained P values were adjusted using the Benjamini-Hochberg correction method. A Pearson's correlation relationship between various environmental variables was performed by "ggcor" package (v. 0.9.6; Huang 2019). To evaluate their associations with the structures of bacterial community (OTU level), constrained analysis of principal coordinates (CAP) based on Bray-Curtis distance was performed through 999 times forward selection. The statistical significance of relationship between each bacterial taxon (based on DNA dataset, RNA dataset, and the functional taxa) at genera level and environmental factors was evaluated with Mantel test, with transformed (ln (x + 1)) and normalized environmental factors (Yang et al. 2017).
Hierarchical clustering analysis of KEGG databases was performed with the correlation distance method based on a correlation matrix (1 − Pearson correlation coefficient) using the "pheatmap" package (v.1.0.12). Differential expression analysis of two groups was performed using the "DESeq2" package (v.1.16.1), genes with an adjusted P value < 0.05 were assigned differentially (Love et al. 2014). A cluster Profiler R package was used to test the statistical enrichment of differential expressed genes in KEGG pathways based on hypergeometric test (Yu and Qing 2016), genes with corrected P value < 0.05 (corrected with Benjamini and Hochberg method) were considered significant differentially expressed. The degree of enrichment in KEGG pathway was analyzed with "ggplot2" package (v.3.1.1), and measured by Rich factor (the ratio of the number of differentially expressed genes in specific pathway to the total number of annotated genes in the pathway), P value (corrected by multiple hypothesis testing), and the number of genes enriched in this pathway.
The differentially expressed KEGG orthologues (KOs) of KEGG pathway were identified through the two-tailed T test using the STAMP (v.2.1.3) with adjusted P values. The relative abundance of differentially expressed KOs associated with key functional activities were used to construct bubble diagram using "ggplot2" package. Simplified schematic of key metabolite pathways in a bacterioplankton cell was illustrated with Pathway builder tool (v.2.0).

Correlations between environment factors and bacterioplankton community
During the blooming period (bP), the colony abundance of colloidal vesicles of P. globosa was 4.1-7.3 m −3 , with average particle size ranging from 1.5 mm to 7.0 mm in diameter. In the recession period (rP), the colony of P. globosa gradually disappeared, whereas bacterial abundance increased from 8.91 to 12.20 × 10 5 cells mL −1 . Environmental parameters were listed in Table 1. Pearson rank correlation analysis showed that the abundance of bacteria had a significant negative correlation with the colony abundance of P. globosa, temperature, dissolved oxygen, and Chl a. (P < 0.05; Fig. 1).
The influences of environment factors on the structures of bacterial communities (OTU level) during the blooming and recession period were further evaluated with CAP analysis (Fig. S2). Results indicated that 36.47% of overall variability of bacterial composition was explained by first two principal components (CAP1 and CAP2). DO, S, and colony abundance were significantly related to changes between bP and rP communities based on 999 times forward model selection (P < 0.01), and explained 12.24% of the variation in the bacterial communities (PERMANOVA test). Additionally, Mantel test was used to examine the potential correlation between environmental variables and the dissimilarity of taxonomic composition, transcriptionally composition, and gene function associatedactive taxon composition (Fig. 1). Bacterial abundance, colony abundance, Chl a, DO, T, and PO 4 3− were strongly correlated with the structures of bacterial community at RNA level (R > 0.5, P < 0.05); DO and colony abundance exerted significant effect on the taxonomic composition of bacterial community at DNA and gene function level, respectively (P < 0.05; Fig. 1 and Fig. S2).

Bacterioplankton community composition and transcriptionally active taxa
Total and active bacterioplankton community structures were assessed by high throughput sequencing-based analyses of the 16S rRNA amplicon from environmental DNA and metatranscriptomic from environmental RNA, respectively. A total of ~ 513 Gbp and ~ 32 Gbp sequences were obtained for the metagenomic and metatranscriptomic data sets, respectively (Table S1 and S2). Taxonomic assignments of the DNA and cDNA protein-coding gene sequences revealed more than 100 phyla in each of the six communities, with majority of the bacterial community was recovered by the surveying effort. The Proteobacteria were the most abundant bacterial phyla across all samples (~ 51.70% of sequences; Fig. S3a), most of which were attributed to Alphaproteobacteria (~ 36.96%) and Gammaproteobacteria (~ 14.45%). The other major groups included Bacteroidetes (~ 23.33%), Cyanobacteria (~ 16.90%), and Actinobacteria (~ 6.91%; Fig. S3a). The relationships of bacterial community were visualized with nonmetric multidimensional scaling (NMDS) analysis (Fig. 2a). The water samples from the blooming period and the recession period formed distinct clusters as confirmed by PERMANOVA (R 2 = 0.2742, P < 0.01). And the bacterial richness and diversity were negatively correlated with the colony abundance of P. globosa (spearman's correlation test, P < 0.01; Fig. 2b, c). At the genus level, the bP community was dominated by Glaciecola and Balneola members. The rP community possessed a high abundance of Candidatus aquiluna, Lentibacter and Candidatus actinomarina (Fig. S4  and Fig. 3). Significant shifts in specific taxa were observed between the bP and rP community. The operational taxonomic units (OTUs) corresponding to Porticoccus, Altererythrobacter, Candidatus puniceispirllum, Aurantivirga, and Ilumatobacter increased significantly in the rP community (two-tailed T test, P < 0.05; Fig. S4).
Taxonomic assignments of cDNA protein-coding gene sequences revealed the Proteobacteria phyla were with the highest transcripts across all samples (Fig. S3b). Among them, Gammaproteobacteria phyla vibrio were the most transcriptionally active populations in the bP community, as more active gene transcription per cell in comparison with other taxa according to the DNA abundance ( Fig. 3 and Fig.  S5). The bP community also had a high transcriptionally active of Pseudoalteromonas and Pseudomonas. Moreover, Aphanizonmenon were dominated the rP community ( Fig. 3 and Fig. S5). The transcriptional activities of Erythrobacter and Candidatus puniceispirillum were significant higher (two-tailed T test, P < 0.05) during the recession period of P. globosa bloom compared to its blooming period (Fig.  S5). Interestingly, Lentibacter and Candidatus actinomarina showed extremely low transcriptional activity in all samples, even the two genera were considerable abundant and classified as dominant taxa across the bacterial communities (Fig. 3).
To determine which members of the bacterial community were metabolism-associated during different bloom periods, significant differentially expressed genes were mapped to contigs with phylogenetic assignments, and the abundances of functions-associated taxa derived from metatranscriptomic sequencing were determined for each period (Fig. 3 and Fig. S6). Overall, 179 genus taxa could be assigned to these genes. The transcriptionally active taxa of Vibrio, which correlated with most functional genes were significant enriched in the bP communities. While the diversity Fig. 2 General distributions of the bacterial community in water samples from the blooming period (bP) and the recession period (rP) of P. globosa bloom. a CSS (cumulative-sum scaling) transformed reads were used to calculate Bray-Curtis distances, and a Nonmetric multidimensional scaling (NMDS) analysis (with 95% confidence ellipses) showed the structure of bacterial community among the samples of different bloom period. Significance of the unconstrained analysis was assessed through PERMANOVA (P < 0.01). The size of each point is proportional to the colony abundance of P. globosa. Drivers of the distributions of bacterial community were estimated via linear least-squares regression analysis between the diversity (b)/richness index (c) and the colony abundance of P. globosa (with 95% confidence intervals). *P < 0.05; **P < 0.01; ***P < 0.001 of contributors to differential genes increased during recession of P. globosa bloom, such as Pseudoalteromonas, Pseudomonas, and Alcanivorax were significantly increased in the rP communities (two-tailed T test, P < 0.05; Fig. S6).

Metatranscriptomic comparison of bacterioplankton community
The total numbers of sequence reads (normalized by library size) assigned to KEGG databases revealed the metabolic potentials and functional activities of bacterial community. The majority of annotated transcripts (~ 72%) were assigned to genes related to metabolism, particularly to three KEGG categories: energy metabolism (~ 14.62%), carbohydrate metabolism (~ 10.94%), and amino acid metabolism (~ 8.40%; Fig. S7a and S7b). The translation pathways (~ 14.51%) were common in the community transcriptome, specifically for the ribosomal pathway. Hierarchical clustering based on protein-coding gene sequences showed the bP and rP community exhibited discrepant metabolic potentials ( Fig. S7a and S7b), there was highly elevated bacterial translation activities in the blooming period of P. globosa while highly elevated metabolism pathways during the recession of P. globosa bloom, mirroring the relatively discrepant environmental conditions (Fig. 1) and community composition ( Fig. 2 and Fig. S2).
At gene transcriptional level, 4,121 of 1,048,575 KEGG genes indicated significant difference across the six communities (P < 0.05; Fig. 4a). More genes had higher levels of expression during the blooming period (2879 genes) than that of the recession period (1242 genes). The number of differentially expressed genes mapped to central metabolic pathways varied across communities (Fig. 4b).
Among the 277 significant differentially expressed KEGG metabolic pathways, 29 were better represented in the bP community (e.g., carbon metabolism, biosynthesis of amino acids, and nitrogen metabolism), and 21 were better represented in the rP community (e.g., ribosome, 2-component system, and quorum sensing; Fisher's exact test, P < 0.05).

Fig. 3
Total taxonomic composition (based on DNA dataset), and transcriptionally active taxa (based on RNA dataset), and gene functions-associated taxa (based on significantly different genes) of the bacterial community. In which taxonomic assignments were binned at the genus level (only the top 15 of each dataset were shown) between bP (blooming period) and rP (recession period). The chord diagram arranged the nodes radially (the upper nodes and lower nodes represents taxonomic composition and bacterioplankton, respectively), drawing thick curves between nodes. Each bacterial taxon had a link for each taxonomic composition in the bP or the rP community, and the thicker the links, the more relative abundant distribution of bacterial genes or transcripts

Diversity and function of cell motility, membrane transport and two-component systems
Expression of cell motility (including flagellar assembly and bacterial chemotaxis), membrane transport (including ABC transporters, phosphotransferase system and bacterial secretion system) and two-component systems within the major bloom-associated bacterial taxa (Figs. 4b, 5) indicated its sensing of environmental resources and adaption to particular ecological niches (Penn et al. 2014). KEGG orthologues (KOs) with significant differential expression between different bloom periods, based on their relative abundance in key KEGG metabolic pathways, were identified through the twotailed T test with adjusted P values. Within the significant differentially expressed KOs, five cell motility-associated transcripts, 36 membrane transport-associated transcripts, and 30 two-component systems-associated transcripts were the most comprehensive and representative pathways, that inducing the sensing and uptake of nutrients, cell motility with flagellar assembly, extracellular secretion of proteins, and other functional activities (Figs. 4, 5, Tables S3-S7, and Fig. S8, S9).

Nutrient assimilation
In the bP community, the genes of general L-amino acid transport system aapM and aapP were highly expressed ABC transporter (Fig. 5a, b), which were correlated with the transcriptionally active taxa of Vibrio ( Fig. S8; and Table S4). In contrast, multiple strategies for utilizing both inorganic and organic sources were observed in the rP community. The utilization of phosphate was significantly increased in the recession period (two-tailed T test, P < 0.05), supported by sensing and responding to phosphorus assimilation among two-component system transcripts (senX3). The utilization of inorganic and organic phosphorous was demonstrated by the up-regulated expression of phosphate transporter (pstC, pstS) and phosphonate transporter (phnC, phnD), respectively (Fig. 5a, b). Further, the uptake and transport of metallic cations (afuA, cusA, cusB, znuB), amino acids (proV, glnH, lysY, aspQ), nitrate/nitrites (frdA, ntrY), monosaccharides (gtsC, xylF, frcA, rhaS), multiple sugar (gguB), and short chain fatty acids (atoB, atoE) were also highly expressed in the diversiform rP microbes (two-tailed T Fig. 4 The genes with significant differential expression levels in the bP (blooming period) and rP (recession period) communities. a Hierarchical clustering based on their relative abundance of differential expressed genes in each community, was performed with a correla-tion matrix (1-Pearson correlation coefficient). b The functional categories of KEGG pathways in the bP community with significantly higher or lower expression levels than those in the rP community. Only the P < 0.05 in each community was displayed test, P < 0.05; Fig. 5a, b) including Pseudomonas, Vibrio, Candidatus puniceispirillum, and Psychrobacter ( Fig. S8 and Table S5).

Cell motility and quorum sensing
Higher activities of cell mobility were indicated in the rP community compared to the bP community (Fig. 4),   Fig. 5 Functional assignment of bacterial metatranscriptome sequences to key metabolic pathways. a Bubble diagram showing differentially expressed KEGG orthologues (KOs) in cell motility, twocomponent system, and membrane transport pathways. Orthologues were represented by circles whose size is proportional to read number. The foldchange value was the ratio of the average relative abundance of bP (bloom period) to rP (recession period) community. The value of fold change > 0, indicated the expressed KOs in bP community higher than those in the rP community. Orthologues with significantly difference (P < 0.05) between bP and rP were shown above the dotted horizontal line (P = 0.05). b The table contains a description of the function followed the gene names for the system based on KEGG annotation. The typical transported substrates. The color of the gene was corresponded to the three metabolic pathways in the bubble diagram. c Simplified schematic representing key metabolite pathways in a bacterioplankton cell. The subcellular localization of gene products highlighted in color was determined by reference to annotation data and the literature (see text for details). Only the genes with significant differential expression between the communities were shown (P < 0.05). Abbreviation: DMS, dimethylsulfide; DMSP, dimethylsulfoniopropionate; MeSH, methanethiol; PO 4 3− , orthophosphate; NO 2 − , nitrite; NO 3 − , nitrate; DOM, dissolved organic matter; POM, particulate organic matter especially in the active heterotrophic Pseudomonas, Labrenzia, Vibrio, and Deferrisoma (Fig. S8). This result was manifested by the up-regulated flagella regulon (fliA) in the twocomponent system, as well as the flagellar hook-basal body complex protein (fliE), flagellar protein (fliS), and flagellar synthesis protein (flgM, flgN) in flagellar assembly pathway (two-tailed T test, P < 0.05; Fig. 5a-c). Further, autoinducers-2 (AI-2) is a type of interspecific signaling molecule that is generally released by exogenous microorganisms and responsible for intercellular communication in quorum sensing (QS). The AI-2 binding periplasmic protein (luxP), which regulates bioluminescence, exhibited higher expression in the bP community. Whereas, ribose transport system substrate-binding protein (rbsB) that transports AI-2 inside the cell through the ABC transport system, was significantly increased in the rP community (two-tailed T test, P < 0.05; Fig. 5a, b). Quorum sensing were better represented in the rP community than bP community (Fisher's exact test, P < 0.05; Fig. 4).

Extracellular secretion
The transcripts of Type II secretory system, which relied on the participation of Sec and GSP protein to secrete virulence factors from the cytoplasm to the extracellular domain, were highly abundant in the cDNA data sets (Table S3). The Sec proteins (secY) were dominant in the bP community (Fig. 5a, b) associated with Vibrio, Alcanivorax, Acinetobacter, and Pseudomonas (Fig. S8), suggesting the extracellular toxins may cross the intima into the periplasm first, then assemble into a mature protein after folding and further modification (e.g., formation of sulfur bonds). The protein will be secreted from the periplasm to the extracellular passage through the outer membrane with the assistance of a complex of GSP proteins (ftsY, gspE, gspF). The expression of GSP proteins were significantly up-regulated in the rP community (Fig. 5a, b), particularly associated with Pseudohongiella, Shewanella, Psychrobacter and ten less active taxa (Fig. S8).

Other functional activities
The bP community harbored one significantly higher orthologue for dimethylsulfoniopropionate (DMSP) demethylase (dmdA; two-tailed T test, P < 0.05; Fig. 5c) associated with Rhodobacteraceae (Table S7), indicating highly active behavior of DMSP catabolism to form methanethiol (MeSH). In addition, arylsulfatase B (ARSB), one of the glycosaminoglycan degradation enzymes, was detected in the rP community (two-tailed T test, P < 0.05; Fig. 5c) associated with the less active taxa of Niabella (Table S7), suggesting the presence of a functional response toward the accumulation of mucopolysaccharide during P. globosa bloom. The exopolysaccharide (EPS) biosynthetic polymerization and transportation genes, such as UDP-glucose 6-dehydrogenase (ugd, key precursor in the synthesis of polysaccharides (Li et al. 2010)), and UDP-glucose 4-epimerase (galE, polysaccharide repeat unit polymerase (Li et al. 2019)), were highly enriched in the six communities (Fig. S9). Among them, UDP-glucosyltransferase (rfaG, a key enzyme in the synthesis of microbial EPS) was significantly increased in rP samples ( Fig. S6 and Fig. 5c).

The difference between total bacterial communities and transcriptionally active taxa
Marine bacterial communities in the outbreak and decline of P. globosa bloom were isolated and amplified, respectively, to perform their taxonomic composition and functional response analyses using an integrated metagenomic and metatranscriptomic approach. The bacterial abundance and diversity increased along with the decline of the P. globosa bloom when massive organic matter was released from Phaeocystis cells and colonies at the peak of the bloom and in the following recession period (Bratbak and Thingstad 1985;Danger et al. 2007). Similar to previous studies (Bagatini et al. 2014;Wemheuer et al. 2014), there were significant differences among bacterioplankton community structures (both of the total and active) during the bloom and recession periods. Active bacterial communities were dominated by only a few marine groups (e.g., Vibrio, Pseudoalteromonas) in the bloom period, whereas there was a higher richness and diversity of bacterioplankton in the recession period of the P. globosa bloom, demonstrating that the bacterioplankton experienced distinct metabolic conditions across the bloom and recession stages.
Heterotrophic bacteria differed in their ability to of consume and remineralize various substrates (Stepanauskas et al. 2003;Alonso-Saez and Gasol 2007), and key gene function-associated bacterial composition showed Pseudoalteromonas, Pseudomonas, and Alcanivorax with a significant increase in the communities during the recession period. The link between community structure and function in these communities suggested that these communities were functionally distinct during different periods of the algal bloom. In addition, an overt difference between community composition (16S rRNA) and the active subset (metatranscriptome) was observed in all communities, which also found in the microbial communities of the mesoand bathypelagic realm in North Pacific Ocean (Wu et al. 2013), revealed differential relative transcriptional activities per cell. Particularly, the taxa Lentibacter and Candidatus actinomarina had high abundances but possessed low transcriptional activity, indicating that the dominant taxa may not necessarily be among the great functional contributors of the bacterial community . Compared to conventional methods, omics techniques provided less bias for microbial gene expression in situ, and provided novel insights into relative transcriptional activity in the microbial community. Moreover, the discrepancies between the metagenomic and metatranscriptomic libraries had been also attributed to the standards of how metabolically active bacterial species were identified, as well as database bias for sequenced genomes in metatranscriptome annotation (Gifford et al. 2011;Vikram et al. 2016).

Bacterial utilization of inorganic and organic nutrients
It is speculated that phosphate limitation was one of the major causes for the decay of the P. globosa bloom in the Beibu Gulf . In this study, phosphate content in water exerted strongly effect on the abundance of P. globosa colonies and the taxonomic composition of transcriptionally active bacteria. Compared the bacterial gene transcription patterns of ABC transporters and two-component systems between the bloom and recession communities, the genes involved in the uptake and utilization of inorganic phosphorus organic phosphate have different expression strategies during the bloom process. With subsiding of P. globosa bloom, the transcripts of phnC and phnD genes which were involved in the transport of organic phosphate, were significantly upregulated, indicating that bacterial uptake and utilization of organic phosphate were increased significantly, likely to satisfy utilization of dissolved organic matter for secondary production (Kirchman and Wheeler 1988). The transcriptional activity of organic phosphate utilization by bacteria was up-regulated in the recession period, which should result in increase of inorganic phosphorus content in water. At the same time, the pstC and pstS genes involved in the uptake and utilization of inorganic phosphorus system, and the senX3 gene which senses and responds to inorganic phosphorus assimilation, were significantly up-regulated during the bloom period, indicating that bacterioplankton compete with phytoplankton for limited dissolved inorganic phosphorus during this stage. And an increase of phosphate content was also detected during the recession period of the P. globosa bloom in the present study, which may due to the bacteria absorbed only 50% of the inorganic phosphorus produced by its decomposition during phosphate-limited environment (Bjorkmen and Karl 1994). The decrease of inorganic phosphorus ingested by Phaeocystis or exogenous supplementation from the estuary was also a possible explanation for the results described here.
Uptake of organic carbon reveals the potential for niche overlap and specialization. In a nutrient-limiting environment, the release of labile, low molecular weight, dissolved organic matter (DOM) (e.g., amino acids, organic acids) from living P. globosa was often induced or increased (Bratbak and Thingstad 1985;Danger et al. 2007). These substances could be utilized by bacteria, thereby relieving the competitive relationship between P. globosa and bacterionplankton during the algal bloom (Buchan et al. 2014). In addition, the released substances may also stimulate the beneficial bacteria to produce phytoplankton growth-promoting compounds, such as vitamins and trace elements. And as evidence, the transcripts of cusA and cusB, involved in copper/silver efflux of bacteria (e.g., Vibrio, Pseudoalteromonas, and Micromonospora) in the two-component system, was indeed up-regulated significantly during the blooming period. During the recession period, the compactness of Phaeocystis colonies decreased along with its growth in size and age (Albright et al. 1986), resulted in release of particulate organic matter (POM) from the colony matrix. At this stage, the transcripts of bacterial genes associated with transport of organic phosphate (phnC and phnD) and uptake of glutamine (glnH), monosaccharides (gtsC, xylF, frcA and rhaS) and polysaccharides (gguB) were significantly upregulated, indicating the release of the mucus polymers and the simple sugar components from the colony matrix (Thingstad and Billen 1994), which in turn enhanced conversion of organic matters. Among the major bacterioplankton taxa, there were high and variable proportions of transcripts associated with sensing and uptake of inorganic nutrients, DOM, and POM suggesting the presence of resource partitioning and active recycling of nutrients within the bloom.

Bacterial motility and extracellular secretion during the P. globosa bloom
Large Phaeocystis colonies with decreased compactness could stimulate bacterial activities in clumping, attachment, and invasion during recession of an algal bloom (Albright et al., 1986). During recession of the P. globosa bloom, the genes associated with flagellum-mediated motility (flgM, flgN, fliA, fliE, and fliS) were highly transcribed, showed that the motility of related bacterial groups (such as Pseudomonas, Alcanivorax, and Candidatus puniceispirillum) is enhanced during period, which may promote the chemotaxis of bacteria and directly affected the colonization of bacteria in the algal surface (Manson 1992). The ribose transport system substrate-binding protein (rbsB), an AI-2 binding protein, might regulate cell motility by acting on target genes related to cell chemotaxis and bacterial flagellum synthesis (Karunakaran and Esther 2012). The rbsB mutant strain could inhibit the formation of bacterial biofilm and further reduce its colonization ability (Armbruster et al. 2011). Here, the transcripts of AI-2 binding protein rbsB were up-regulated expressed in the bacterial communities during bloom recession period. The related microbial groups (Vibrio, Planktomarina, Labrenzia, and Rhizobium) may mediate the surface translocation of those bacteria by quorum sensing (QS) in a competitive environment. In addition, the QS system regulated the synthesis of virulence factors (e.g., gelatinase, lipase, lecithinase, and hemolysin; Bai and Zhang 2010) and motility regulons were also able to control the expression of virulence determinants (Mattick 2002). In the present study, seven genes of Type II secretory system were expressed significantly high across the bloom and recession stages, through two processes of extracellular toxins entering the periplasmic space (secY) or secreted proteins from the periplasm to the outer membrane (ftsY, gspE, gspF, etc.), respectively. And these genes could induce secretion of extracellular virulence factors during both periods, which was closely related to the subsequent degradation of algal cell walls, growth inhibition and its subsequent mortality (Pokrzywinski et al. 2012). Although the extracellular enzyme content was not continuously monitored during bloom, the elevated transcripts of several related genes (such as, rfbB, galE, and rfaG) corresponding to bacterial EPS synthesis in both periods could be evidence for their presence in environmental water.
In the laboratory and simulation in natural seawater, the taxa of Pseudomonas (Dakhama et al. 1993), Vibrio (Zhao et al. 2009), Pseudoalteromonas (Seymour et al. 2009), and Pseudohongiella (Cai et al. 2009) were confirmed to secrete toxic substances into the environment inhibiting the growth of certain algae such as cyanobacteria, dinoflagellates, and diatoms. In particular, bacterial groups of Bacillus and Streptomyces were had strong algicidal activity against P. globosa, through overproduction of extracellular active substances (Guan et al. 2015;Hou et al. 2016;Tan et al. 2016) and resulting in growth inhibition of P. globosa. In this study, the transcriptionally active genes in the recession period were mainly associated with the bacterial taxa of Bacillus and Streptomycesxia, indicated that these bacteria were closely related to inhibition of P. globosa growth and recession of its bloom.
In summary, bacterioplankton were able to decompose the colony and cells of P. globosa and release the intracellular substances of Phaeocystis (such as monosaccharides, polysaccharides, amino acids, etc. in the colony matrix), thus provided nutrient basis and conditions for the growth of bacterioplankton in bloom environment (Zheng et al. 2013;Zhang et al. 2015). In addition, the genes putatively involved in flagellum-associated chemotaxis, AI-2 quorum sensing signal, and participated in regulation of polysaccharide secretion, which may further induce bacterial adhesion, survival, movement and synthesis of virulence factors (Manson 1992;Stoodley et al. 2002). Therefore, the high expression of genes involved in nutrient uptake and key cellular processes (e.g., bacterial motility, quorum sensing, and extracellular secretion) facilitated their advantage in finding and competing for carbon sources and energy sources during the recession period (Harshey 2003), which also promoted recession of the algal bloom.

Conclusions
As the decline of P. globosa bloom, the pivotal roles of bloom-associated bacterioplankton in transforming phytoplankton-derived organic matter and gene transcriptional responses were explored in the present study. Taxonomic analysis uncovered the bacterial communities experienced distinctly different metabolic conditions across the bloom onset and collapse stages. Transcriptome analysis highlighted a potential mechanism for variable nutritional utilization, which was mediated through cell motility, membrane transport, and two-component systems. Bacterioplankton was deeply involved in the uptake of dissolved and particulate organic matter during recession of P. globosa bloom. Additonaly, the bacterial community transcribed more copies of genes associated with bacterial motility, quorum sensing and extracellular secretion during the recession period, indicating their advantage in location, competition and degradation of carbon sources. The results demonstrate that the gene transcriptional profiles of bloom-associated bacterioplankton are closely related to the site physiochemical characteristics, providing insights into the bacterial response and adaptation mechanisms that dictate their behaviors across 'boom' and 'bust' period of harmful algal blooms.