Metagenomic Insights Into Co-Proliferation of Vibrio Spp. and Dinoagellates Prorocentrum During A Spring Algal Bloom in the Coastal East China Sea Near Xiamen

In this study, we interrogated the prokaryotic microbiomes of surface water samples collected at two neighboring segments of East China Sea that contrast greatly in terms of the intensity and frequency of Prorocentrum-dominated HAB. Mantel tests identied signicant correlations between the structural and functional composition of the microbiomes and the physicochemical state and the algal biomass density of the surface seawater, implying the possibility that prokaryotic microbiota may play key roles in coastal HAB. A conspicuous feature of the microbiomes at the sites characterized with high trophic state index and eukaryotic algal cell counts was disproportionate proliferation of Vibrio spp., and their complete domination of the functional genes attributable to the dissimilatory nitrate reduction to ammonia (DNRA) pathway signicantly enriched at these sites. The genes attributed to phosphorus uptake function were also signicantly enriched at these sites, presumably due to the P i -deciency induced by algal growth; however, the proles of the phosphorus mineralization genes lacked consistency, barring any conclusive evidence with regard to contribution of prokaryotic microbiota to phosphorous bioavailability. The results of the co-occurrence network analysis performed with the core prokaryotic microbiome also supported that the observed proliferation of Vibrio and HAB may be causally associated.


Abstract Background
Coastal harmful algal blooms (HAB), commonly termed 'Red tides', have severe undesirable consequences to the marine ecosystems and local shery and tourism industries. Increase in nitrogen and/or phosphate loading is often attributed as the major culprits of increasing frequency and intensity of the coastal HAB; however, fundamental understanding is lacking as to the causes and mechanism of bloom formation despite decades of intensive investigation.

Results
In this study, we interrogated the prokaryotic microbiomes of surface water samples collected at two neighboring segments of East China Sea that contrast greatly in terms of the intensity and frequency of Prorocentrum-dominated HAB. Mantel tests identi ed signi cant correlations between the structural and functional composition of the microbiomes and the physicochemical state and the algal biomass density of the surface seawater, implying the possibility that prokaryotic microbiota may play key roles in coastal HAB. A conspicuous feature of the microbiomes at the sites characterized with high trophic state index and eukaryotic algal cell counts was disproportionate proliferation of Vibrio spp., and their complete domination of the functional genes attributable to the dissimilatory nitrate reduction to ammonia (DNRA) pathway signi cantly enriched at these sites. The genes attributed to phosphorus uptake function were also signi cantly enriched at these sites, presumably due to the P i -de ciency induced by algal growth; however, the pro les of the phosphorus mineralization genes lacked consistency, barring any conclusive evidence with regard to contribution of prokaryotic microbiota to phosphorous bioavailability. The results of the co-occurrence network analysis performed with the core prokaryotic microbiome also supported that the observed proliferation of Vibrio and HAB may be causally associated.

Conclusions
The ndings of this study suggest a previously unidenti ed association between Vibrio proliferation and the Prorocentrum-dominated HAB in the subtropical East China Sea, and opens a discussion regarding an unlikely, but still possible, involvement of Vibrio-mediated DNRA in Vibrio-Prorocentrum symbiosis. Further substantiation of this interaction with culture-based experiments may prove crucial for understanding of the dynamics of explosive local algal growth during HAB events endemic to the region.

Background
Today, many coastal regions across the globe suffer from harmful algal blooms (HAB), more widely known as the layman's term 'red tides' [1]. The East China Sea off the southeast coast of China has been one of the regions most severely affected by HAB since the 1980s [2]. Since the turn of millennium, the major culprits of HAB in the East China Sea have shifted to the dino agellates Prorocentrum spp., a group of fast-growing phototrophic algae widespread in coastal waters across the globe [3]. The species prevalent in HAB in East China Sea, which had been named P. donghaiense, is neither morphologically nor genotypically clearly distinguishable from other Prorocentrum species with worldwide distribution, e.g., P. dentatum and P. obtusidens [4,5]. Although toxin production by this particular species has not been reported, the possibility cannot be totally ruled out, as several Prorocentrum species have been observed to produce okadaic acid, dinophysistoxin-1, and prorocentrolide [6,7]. Even apart from toxin production, intense Prorocentrum blooms cause substantial detrimental impacts to the local marine and estuarine ecosystems by means of oxygen depletion, alteration of pH, and/or light attenuation, which are typical for high-biomass algal blooms [1,8,9]. Such ecological damage, as well as aesthetic deterioration due to coloration and putrefaction accompanied with the blooms, have caused immense damage to the local tourism, shing, and aquaculture industries [10].
The annually repetitive nature of the Prorocentrum algal bloom in the East China Sea, as well as explosive growths of P. donghaiense observed during progressions of the algal blooms endemic to the region, suggests that in situ growth is more important than vertical and/or horizontal migration in understanding the mechanisms of bloom dynamics in the region [11,12]. As P. donghaiense is primarily a photoautotrophic organism with carbon needs mostly ful lled from uptake of inorganic carbon, the environmental factors regarded to be crucial to its overgrowth during the bloom events are presumably availability of nitrogen and phosphorous, as suggested by previous laboratory observations [1,13,14].
Although several published data collected from the East China Sea suggest potential correlations between nitrogen and phosphorous availability and P. donghaiense blooms, conclusive experimental evidences are still lacking, and it is yet unclear whether or which speci c forms of nitrogen and phosphorous support rapid growth of P. donghaiense [15,16]. As such, the decades-long question as to whether or how relief of nitrogen and/or phosphorous limitation initiates P. donghaiense bloom remains unanswered.
In this study, we approached this same old question from a different perspective, exploring the possibility that the sea surface microbiome may play a key role in formation and progression of Prorocentrum blooms. Diverse bacterial and archaeal constituents of the sea surface microbiota are involved with cycling of nitrogen and phosphorous [17]. In the euphotic zone, generally characterized by supple sunlight and oxygen availability, photosynthetic cyanobacteria x atmospheric nitrogen to ammonium and organic nitrogen, and diverse bacteria and archaea catalyze aerobic mineralization of organic N and oxidation of NH 3 and NO 2 − [18 -20]. Typically, NO 3 − is often the most abundant form of nitrogen at the sea surface due to its oxidizing condition [21,22]. The microbial reactions perceived as predominantly anaerobic, e.g., denitri cation, anaerobic ammonia oxidation (anammox), and dissimilatory nitrate/nitrite reduction to ammonium (DNRA), may also occur during periodic hypoxia, e.g., nocturnal anoxic spells during high-biomass algal blooms and periods of massive biomass decay following demise of algal blooms [23,24]. Such microbial mediated nitrogen transformation reactions may affect nitrogen supply to phytoplankton, as they may alter the overall nitrogen availability, as well as nitrogen speciation. As nitrogen is often regarded as the limiting nutrient for algal growth during rapid blooms and previous laboratory experiments have shown that certain nitrogen species are favored by Prorocentrum spp. over others, investigations into these microbially mediated nitrogen-transformations may help improve understanding of the bloom dynamics in the East China Sea [3,13,25].
Less predictable is the in uence of microbiota on HAB via alteration of the marine phosphorus cycle, given the patchy ecophysiological understanding of marine microbial phosphorous metabolism. Previous metagenomic analyses have revealed abundance and omnipresence of diverse phoA-, phoD-, and phoXtype alkaline phosphatases (AP) and C-P lyases and hydrolases in seawater microbiomes, which have been known to be involved with disintegration of inorganic phosphorus (P i ) from phosphoesters and phosphonates [26][27][28]. Expressions of these organic phosphorus metabolizing enzymes are generally perceived to be under Pho regulon control, such that they are expressed only under P i -de ciency; however, P i -insensitive expression and activity of alkaline phosphatases and C-P hydrolases have been experimentally demonstrated in both in situ and laboratory experiments [29,30]. Thus, bacterial mineralization of organic phosphorus may enrich the inorganic phosphorus pool in situ and increase phosphorous bioavailability to phytoplanktons, contributing to formation of HAB [27,31,32]. Marine eukaryotic phytoplankton, including Prorocentrum spp., also synthesize active AP [32][33][34]; however, a recent eld study conducted in the East China Sea during a late-spring P. donghaiense bloom showed that the bacterial fraction accounted for the majority of AP activity and that the contribution of eukaryotic algae was minimal [35].
Here, a comparative survey of the sea surface microbiomes was performed at two neighboring coastal seas off the coastline of Fujian Province in South China. The sea surrounding Pingtan Island, hereafter referred to as Pingtan Sea, has been regularly and frequently affected by P. donghaiense blooms. In contrast, Xiamen Bay ~ 200 km to the southwest has stayed mostly free of Prorocentrum spp., although occasional occurrences of moderate diatom and dino agellate blooms have been reported [36]. During a typical spring dino agellate bloom event in the Pingtan Sea in 2019, surface seawater samples were collected from ve locations in Pingtan Sea and three locations in Xiamen Bay, varying in the intensity of HAB at the time of sampling, and the physicochemical properties of the seawater and population densities and compositions of the eukaryotic algae were analyzed. The microbiomes of these seawater samples were then analyzed with shotgun metagenome sequencing. The compositional and functional features of the microbiomes were analyzed, with particular emphasis on those relating to nitrogen and phosphorous metabolism and statistically associated with the severity of HAB in the region. The ndings from the comparative metagenomics analysis warrant further investigation into interactions between key microbial groups, e.g., Vibrio spp., and eukaryotic algae, for a more complete understanding of HAB dynamics in coastal seas.

Results And Discussion
Characterization of the sampling sites Notable differences were observed between the physicochemical characteristics of the surface water samples collected from the Pingtan Sea (PIN1-PIN5) and the Xiamen Bay (XIA1-XIA3) (Fig. S1, Table S1).
The surface water temperature was higher at the Xiamen Bay sites by approximately 3 °C at the time of sampling, and the pH and salinity were slightly but signi cantly higher (p<0.05) at the Pingtan Sea sites. The dissolved oxygen (DO) concentration of the surface water ranged between 4.7 and 6.9 mg/L, but no signi cant difference was found between the two regions (p>0.05). The measured oxidation reduction potential (ORP) values were 43.9±2.5 mV and 22.1±6.5 mV at the Pingtan Sea and the Xiamen Bay sites, respectively (Fig. S1). Overall, Xiamen Bay was more nutrient-rich than the Pingtan Sea. The dissolved inorganic nitrogen (DIN) concentration was signi cantly higher (p<0.05) in the Xiamen Bay samples  (Table S2). The only other site with a substantial Prorocentrum spp. population was PIN5, where 3.6 x 10 5 cells/L, out of the total algal population of 5.2 x 10 5 cells/L, were identi ed as Prorocentrum spp. The total 16S rRNA gene copy numbers, as an indicator of prokaryotic population, ranged between 9.9±1.3 x10 7 and 4.2±0.2 x10 8 copies/L in Pingtan Sea and between 9.1±0.9 x10 7 and 2.9±0.9 x10 8 copies/L in Xiamen Bay (Fig 1b and Table S2).

Analysis of prokaryotic community and core microbiome designation
Analysis of the prokaryotic community structures of PIN1-5 and XIA1-3 identi ed Proteobacteria as the dominant phylum, with its relative abundance ranging from 76.4% to 88.2% (Fig. 2a). Alphaproteobacteria and Gammaproteobacteria constituted 18.0 -34.6% and 41.5 -70.0% of the prokaryotic communities, respectively. Bacteroidia were also abundant, constituting 6.3-17.0% of the total prokaryotic community. Although archaeal OTUs a liated to Thermoplasmata and Nitrososphaeria were found in the metagenomes, the cumulative relative abundance of the archaeal OTUs was below 1.5% in any of the samples analyzed. The beta diversity analysis of the sea surface microbiomes clearly set the Pingtan sites apart from the Xiamen sites ( Fig. 2c-e). The non-metric multidimension scaling (NMDS) plots constructed via taxonomic grouping of 16S rRNA gene sequences and functional assignments of translated enzyme sequences to KEGG Orthology (KO) groups and CAZymes groups invariably separated PIN1-5 from XIA1-3 ( Fig. 2c-e), suggesting disparities between the microbiomes of these two regions, in terms of both taxonomic and functional make-ups. The OTUs identi ed in the samples were assigned to 379 bacterial and 5 archaeal taxa (according to genus-level classi cation), and 47 of these taxa, shared by all 8 metagenomes, were further considered as the members of the core microbiome (Table S3). The core microbiome consisted of 46 bacterial taxa and one archaeal taxon, most of which are uncultured and understudied. The members of the core microbiome constituted vast majority of the total prokaryotic community in each of the analyzed samples, such that their cumulative relative abundance ranged between 70.4 and 76.9% at the Pingtan Sea sites and between 62.6 and 73.6% at the Xiamen Bay sites ( however, its low abundance in PIN4 (1.3%) suggests relevance of this taxon to eutrophication or algal bloom unlikely.

Metagenomic analysis of nitrogen and phosphorus metabolism
The metagenomic pro les of nitrogen metabolism-related genes in the microbiomes suggest possible involvement of the prokaryotic microbiome in altering the availability of labile nitrogen to Prorocentrumdominated algal populations ( Fig. 3 and Table S4). Apart from the genes involved with urea uptake and mineralization (ureABC and urtABCDE), all of the most abundant nitrogen functional genes were relevant to dissimilatory and assimilatory NO 3 --to-NO 2 reduction (napAB and nasA) and dissimilatory nitrite reduction to ammonium (DNRA; nrfA and nirB). The genes encoding for the periplasmic nitrate reductase, napAB, had substantially higher relative abundance at the Pingtan Sea sites (10.2-39.8 RPKM for napA and 6.0-27.3 RPKM for napB) than in the Xiamen Bay sites (6.9-17.1 RPKM for napA and 8.0-14.2 RPKM for napB), although the differences were not statistically signi cant (p>0.05). The nirB and nrfA genes, encoding the catalytic subunits of the two distinct forms of nitrite reductases central to the DNRA pathway, showed high relative abundance at the Pingtan Sea and Xiamen Bay sites [39]. The relative abundance of nirB ranged between 26.6 and 87.5 RPKM in the Pingtan Sea metagenomes, and 17.5 and 50.4 RPKM in the Xiamen Bay metagenomes. The nrfA genes, recovered in lower abundance than nirB in all metagenomes, were signi cantly more abundant in the Pingtan Sea metagenomes (12.7±8.9 RPKM) than in Xiamen Bay metagenomes (3.1±3.5 RPKM). Disproportionately large shares of napA, nirB, and nrfA genes were phylogenetically a liated with Vibrio spp., considering that the relative abundance of Vibrio spp. did not exceed 33% of the total prokaryotic population according to the 16S rRNA-based community composition analyses. At the Pingtan Sea sites where this tendency was more pronounced, >93.2%, >95.8% and >56.1% of the reads mapped onto napA, napB, and nirB, respectively, were phylogenetically a liated with the genus Vibrio. The gene-coding sequences identi ed as nrfA were invariably phylogenetically a liated with Vibrio spp. in all metagenomes, and the relative abundance of Vibrio spp. exhibited signi cant correlation with that of nrfA (p<0.05 according to the Spearman's rank correlation).
Compared to the functional genes involved with NO 3 --to-NO 2 reduction and DNRA, the functional genes encoding key enzymes for nitri cation and denitri cation were, in general, recovered in relatively low abundance (Fig. 3c) 1 and 532.3, respectively). This selective enrichment of nrfA in the Pingtan Sea metagenomes appear too heavily biased to be solely attributable to the high relative abundance of Vibrio spp.
Metagenomic interrogation of phosphorus uptake and metabolism also identi ed several notable differences between the Pingtan Sea and Xiamen Bay microbiomes ( Fig. 3a and Fig. 3d). The ugpABCE genes, encoding ATP-binding cassette (ABC) transporters for uptake of glycerol-3-phosphate and glycerophocholine, were found in signi cantly higher relative abundance in the Pingtan Sea metagenomes than in the Xiamen Bay metagenomes. The phosphonate transporter gene phnD and the high-a nity phosphate binding protein pstS, were both also signi cantly more abundant in the Pingtan Sea microbiomes (51.7±15.0 RPKM and 75.5±21.6 RPKM, respectively) than in the Xiamen Bay microbiomes (27.4±5.7 RPKM and 50.9±11.8 RPKM, respectively). As these phosphorus uptake genes are known to be expressed under P i -de cient conditions to cope for phosphorus starvation, their abundance in Pingtan Sea metagenomes was consistent with the low concentrations of soluble reactive phosphorous (7.9±4.3x10 -3 mg/L) at the Pingtan Sea sites [40][41][42]. Substantial proportions of phnD (35.5±11.5%), pstS (43.3±15.0%), and ugpB (34.8±15.8%) sequences recovered from the Pingtan Sea metagenomes were a liated with Vibrio spp., suggesting that possession of these diverse means of alternative phosphorus uptake may have contributed to the disproportionate proliferation of Vibrio spp. observed in Pingtan Sea.
The genes encoding three distinct forms of alkaline phosphatases, phoA, phoX, and phoD, were found in all metagenomes, and phoA and phoD, both with RPKM values above 50, were among the most abundant nitrogen-or phosphorus-related functional genes in the metagenomes (Fig. 3d). The phoA genes were signi cantly more abundant in the Pingtan metagenomes (68.7±8.4 RPKM) than in the Xiamen Bay metagenomes (49.5±2.6 RPKM; p<0.05); however, the abundances of phoX and phoD were not signi cantly different between the two regions (p>0.05). In the Pingtan metagenomes, Vibrio-a liated sequences accounted for the majority of phoX (62.3±13.4%), while only 13.6±7.0% and 0.8±0.5% of phoA and phoD, both at least 3.4-fold more abundant than phoX, were a liated with Vibrio spp. Therefore, the abundances of alkaline phosphatase genes appear to have little relevance to Vibrio proliferation or Prorocentrum bloom.
Vibrio-dino agellate relationship as inferred from the microbial network analysis of the core microbiome Microbial network analysis of the core microbiome further substantiated the alleged association between abundance of Vibrio spp. and dino agellate bloom (Fig. 4). The relative abundances of two core microbiome taxa, including Pseudoalteromonas and an unidenti ed genus of the family Cryomorphaceae, showed signi cant positive correlation with the algal biomass, and the relative abundances of four other core microbiome taxa exhibited signi cant negative correlation (Spearman's |r| >0.7, p<0.05; Fig. 4). The core microbiome taxon a liated to the genus Vibrio was positively but not signi cantly correlated with algal biomass (r=0.26, p=0.53); however, its relative abundance had strong signi cant positive correlation (r=0.81, p=0.015) with the relative abundance of Pseudoalteromonas, and signi cant negative correlations with those of three core microbiome taxa negatively correlated with algal biomass, including Marine Group II of the phylum Euryarchaeota, an unidenti ed genus of the family Parvibaculaceae, and the IS-44 subgroup of the family Nitrosomonadaceae (r=-0.71, p=0.047; r=-0.71, p=0.047; r=-0.79, p=0.021, respectively). Of note, co-abundance of Pseudoalteromonas spp. with eukaryotic algae has been previously observed and mechanistically explained with basis on the capability of the microorganisms belonging to this taxon to prey on eukaryotic algae [43]. Alongside the observation that Vibrio spp. was the single most abundant taxon in the PIN4 microbiome most severely affected by HAB, these co-occurrence patterns surrounding Vibrio in the core microbiome network corroborate the possibility of a causal association of Vibrio spp. with HAB, despite the lack of statistically signi cant correlation between the relative abundance of the Vibrio and the algal biomass.
Ecological association of Vibrio spp. with dino agellate algal bloom As these metagenome analyses suggest, association of Vibrio spp. to the Prorocentrum-dominated algal blooms was apparent in the examined segment of coastal East China Sea. Strong signi cant correlations between the abundance of Vibrio spp. and the environmental metrics for eutrophication (TSI; r=0.833, p=0.011) as well as the chlorophyll a concentration (r=0.762, p=0.028) also support that the Vibrio and algal proliferation may be mechanistically correlated (Table S5). These observations are in line with a number of previous studies that have reported occurrences of Vibrio blooms in coastal seas simultaneous with blooms of phytoplanktons, e.g., cyanobacteria, diatoms, and dino agellates, or immediately following their demise [44][45][46][47][48]. In a study that monitored microbial population dynamics over a six-year period off the coast of Plymouth, UK, the time point at which a single Vibrio sp. dominated the prokaryotic population, representing 54% of the entire pool of 16S rRNA sequences, coincided with a peak in the population of a diatom species Chaetoceros compressus [44]. Another study performed during raphidophyte blooms in Delaware's inland bays identi ed signi cant correlation between the abundances of particle-associated Vibrio spp. and the raphidophyte population [45].
These previous studies have attributed the Vibrio-algae association to protection from zooplankton predation and/or provision of algal exudates or detritus as organic substrates to opportunistic fastgrowing Vibrio spp. [45][46][47]. Although Vibrio spp. may form bio lms attached to microalgal cells and detritus, physically and chemically protected from grazing, recent culture-independent studies have repeatedly identi ed particle-associated fractions as minor fractions of Vibrio population in seawater [49]. Thus, supply of labile organics, photosynthetically produced and released by algal counterparts, is deemed more likely as the primary mechanism via which algal blooms accommodate disproportionate proliferation of Vibrio spp. [50]. Further, several Vibrio species have been physiologically characterized as alginolytic, and thus may be able to utilize dead algal biomass or even be algicidal. In this study, core microbiome analysis identi ed a strong signi cant positive correlation between Vibrio and Pseudoalteromonas, a taxon widely-known to include algicidal bacteria [51]. This correlation may be interpreted as two groups of specialists sharing a habitat abound with their common substrates, i.e., live or dead algal biomass and/or algae-derived organic materials. It may also be possible that smaller organic compounds produced from Pseudoalteromanas-mediated algal biomass digestion may serve as more labile organic substrates for fast-growing Vibrio spp. [50]. Either hypothesized mechanism would be consistent with the co-occurrence of Vibrio proliferation and Prorocentrum-dominant algal bloom observed in this study.
Whether the presumed Vibrio-Prorocentrum association may be a commensal relationship, bene tting only Vibrio spp., or a mutually-bene cial symbiotic relationship, has not been previously discussed. The comparative metagenomic analysis of the functional genes related to nitrogen metabolism suggest possible scenarios how Vibrio proliferation may provide positive feedback to growth of Prorocentrum spp. and other eukaryotic microalgae. The proportion of Vibrio-associated functional genes were, in general, higher at the Pingtan Sea sites than at the Xiamen Bay sites, probably due merely to the higher Vibrio abundance; however, Vibrio's exclusive ownership of DNRA-related functional genes, i.e., napAB and nrfA, and the near complete absence of denitri cation-related functional genes (nirK and nirS) that are often found in genomes of marine heterotrophic bacteria (e.g., algicidal Pseudoalteromonas spp.) are not explicable merely with the high Vibrio abundance. The DO concentrations of the bulk seawater at the Pingtan Sea sites, including the PIN4 site with the highest Prorocentrum population and algal biomass, indicated that the bulk seawater was oxic at the time of sampling, which, due to logistical reasons, took place during the daytime. Under oxic conditions, DNRA is not likely to occur, as it is known as a strictly anaerobic process; however considering the high populations of the fast-metabolizing Vibrio and Prorocentrum that may be potent O 2 sinks at dark hours, it is plausible that periodic oxic-anoxic shifts may occur in the region, especially at highly eutrophicated locales such as the PIN4 and PIN5 sites [52,53]. Spatial microaerobic or anoxic niches may also exist, such as in bio lms attached to particle surfaces, microbial aggregates, and intracellular spaces within zooplankton or phytoplankton where Vibrio spp. may be able to colonize [52,[54][55][56]. Rapid decay of algal carcass, and/or upwelling of hypoxic bottom water may also cause temporal anoxia that may not be captured by snapshots of DO measurement [57,58]. Further, the redox potentials measured in situ at the Pingtan Sea (43.9±2.5 mV) sites were in fact, substantially lower than the standard redox potentials (at pH 7.0) of NO 3 --to-NO 2 -(420 mV) and NO 2 --to-NH 4 + reduction (440 mV), despite the DO concentrations ranging between 4.73 -6.89 mg L -1 . Thus, it is highly likely that at least some of the napABand nrfApossessing Vibrio spp. are DNRA-active, producing NH 4 + from NO 3 and possibly supplying Prorocentrum spp. with NH 4 + , generally known as a more e cient nitrogen source for fast-growing dino agellates than NO 3 - [59]. Although any decisive conclusion cannot be reached solely based on the metagenome analysis, the ndings from this study are certainly compelling indications that NH 4 +production via DNRA may be key to the hypothesized symbiotic association between Vibrio spp. and Prorocentrum spp. and warrant further attention to the role that DNRA may play in Prorocentrumdominated HAB.

Conclusions
Occurrences of severe marine algal blooms are often attributed to nitrogen and/or phosphorus enrichment relieving nutrient limitations to algal growth; however conclusive evidence of the presumed correlation between the nutrient conditions and algal population has remained elusive. In this study, we attempted a different avenue of approach, using comparative metagenomics to investigate any association between surface seawater microbiomes and Prorocentrum-dominated algal blooms in East China Sea. The most prominent feature that distinguished the microbiomes of Pingtan Sea with frequent severe Prorocentrum blooms from that of Xiamen Bay much less frequently affected was enrichment of organisms a liated with Vibrio spp. and disproportionately biased a liation of napAB and nrfA to this taxon. The alleged association between DNRA-capable Vibrio spp. and algal bloom in this locale was further substantiated via core microbiome network analysis. The ndings suggest a previously unidenti ed causal association between Vibrio proliferation and dino agellate bloom in subtropical marine environment and broach the possibility that DNRA may play a key role in mutually bene cial symbiosis between the two groups of organisms.

Sample collection and in situ treatment and characterization
Sampling was performed at ve locations along the coast of Pingtan Island (119°49′-119°52′E, 25°28′-25°36′N) and three coastal locations in the Xiamen Bay (118°0′-118°10′E, 24°25′-24°36′N) in the mornings of April 27 and 29, 2018, respectively (Fig. 1a). Bulk seawater samples were collected from the surface layer (0-1 m depth) using hydrophores. The physicochemical properties of the seawater samples, i.e., water temperature, dissolved oxygen concentration, conductivity, salinity, oxidation-reduction potential and pH, were measured on site immediately after sampling, using Professional Plus handheld multiparameter meter (YSI Inc., Yellow Springs, OH). Secchi disk transparency (SD) was determined on site according to the standard protocol [60]. Ten liters of each seawater sample was ltered through a phytoplankton net with 20 μm mesh size (Xuanmingyu, Wuhan, China) and one liter of the ltrate was

Microbial community analysis and metagenomic analysis of nitrogen and phosphorus metabolism
Low-quality reads were removed from the raw sequence data using Trimmomatic v0.36 with the parameters set as follows: LEADING: 3, TRAILING: 3, SLIDINGWINDOW: 4:15, and MINLEN:70 [71]. The quality-trimmed reads were screened for putative 16S rRNA gene fragments, using the hmmsearch-based Meta-RNA 3 software, with the parameters molecule type and e-value set to "ssu" and "1E-10", respectively [72]. These reads were reconstructed to full-length 16S rRNA gene sequences using EMERGE through 100 iterations with Silva SSU database release 132 as reference [73]. A QIIME-compatible OTU table was constructed with these full-length 16S rRNA gene sequences and the relative abundances of the OTUs were computed using an in-house python script, which converted the normalized coverage of each EMIRGE-synthesized 16S rRNA gene to an OTU count. The OTUs were assigned taxonomic classi cation using RDP classi er v2.10.2 with Silva SSU database 132 as the reference dataset and the parameter min_con dence set to 0.8. The OTUs assigned to chloroplasts were removed for the downstream analysis [74].
For metagenomic analysis of the functional genes relevant to nitrogen-and phosphorous-metabolism, the quality-trimmed sequence reads were assembled de novo using metaSPAdes v3.14.0 with the parameters set to default values, and gene-coding sequences in the contigs were predicted using Prodigal v2.6.3 [75,76]. The predicted coding sequences were assigned functional and taxonomic annotations using DIAMOND v0.9.31.132 with NCBI's non-redundant protein database (accessed on 1/2/2020) and the Uniref90 database (accessed on 4/9/2020) as reference and the parameters set to default values [77]. For each gene-coding sequence, only the hit that returned the highest score was taken for downstream analyses. GhostKOALA (https://www.kegg.jp/ghostkoala/) assigned KEGG Orthology (KO) numbers to the gene-coding sequences [78]. For annotation of carbohydrate-active enzymes (CAZymes), the predicted gene-coding sequences were searched against dbCAN HMMdb release 9.0 using the hmmscan command in HMMER v3.1b2 [79]. Assembled contigs from the individual samples were pooled together, and for each sample, quality-trimmed raw sequence reads were mapped onto these contigs using the BWA mem v0.7.17 software (default parameters). Read counts were performed using HTSeq v0. 12.4 [80, 81]. The raw counts of the coding sequences were normalized with the nucleotide length (in kilobases) and the total read count of the sample analyzed (in million reads), yielding gene coverage values in RPKM (reads per kilobase per million mapped reads).

Statistical analysis
Statistical analyses in this study were performed using nonparametric methods, due to the small sample size (n≤8). All statistical analyses were performed using R v4.0.2 [82]. Pairwise comparisons of the environmental parameters between sampling sites were statistically evaluated with Student's t-test using wilcox.test function implemented in the "stats" package. Spearman's rank correlation tests were performed to evaluate statistical signi cance of pairwise correlations between the environmental parameters. Mantel tests evaluated pairwise correlations among physicochemical parameter matrices, the algal biomass data, the microbial community pro les, the KEGG Orthology (KO) number-based functional gene abundance pro les, and the CAZymes pro les. The dissimilarity matrices for these datasets were constructed using either Bray-Curtis dissimilatory metrics (vegdist function in "vegan" package) or Euclidean distance metrics (dist function in "stats" package), and these dissimilatory matrices were used as the inputs for the Mantel tests performed using the mantel function implemented in "vegan" package, with the 'correlation method' and 'permutation' options set to spearman and 9999, respectively.
The non-metric multidimensional scaling (NMDS) analyses of microbial communities and functional gene compositions were performed based on the Bray-Curtis dissimilarity matrices constructed using metaMDS function embedded in "vegan" package. The co-occurrence networks were constructed based on the pairwise Spearman's correlation analyses performed with the relative abundances of the selected Availability of data and materials The raw sequence data generated in the current study have been deposited in NCBI SRA database under accession number of PRJNA739281.

Competing interests
The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.