Synergy in viral–bacterial co-infection expedites algal bloom demise

Isolation of alga, virus and BAM from a G. huxleyi bloom

G. huxleyi strain RCC6946, the virus EhV-M1 and the BAM were isolated from an induced G. huxleyi bloom during the AQUACOSM VIMS-Ehux mesocosm experiment conducted near Bergen, Norway, in 2018²⁴. For algal isolation, a water sample collected from mesocosm bag 4 during the bloom phase (day 13, 6 June 2018) was used for single-cell isolation within two weeks of sample collection at the Roscoff Culture Collection (RCC). The BAM was isolated from mesocosm bag 7, which exhibited the highest increase in viral abundance during G. huxleyi bloom demise. A water sample was collected during the demise phase (day 23, 16 June 2018), filtered through GF/F glass microfibre filters (~0.7 µm nominal pore size), and the filtrate cryopreserved with glycerol (15% final concentration). The lytic virus EhV-M1 was isolated from the same mesocosm experiment using G. huxleyi CCMP 374 as a susceptible host, and its genome was sequenced⁶².

Maintenance of G. huxleyi cultures and viral infection experiments

The G. huxleyi strain RCC6946 was maintained in batch cultures at 18 °C with a 16 h:8 h light:dark cycle. Light was provided at 100 μmol photons m⁻² s⁻¹ using cool white LEDs. The medium was composed of autoclaved filtered seawater (FSW) adjusted to 35‰ and supplemented with modified K/2 medium⁶³, modified as described in Supplementary Data 3. In parallel to the xenic culture, an axenic culture was derived by adding the antibiotics ampicillin and kanamycin at final concentrations of 100 μg ml⁻¹ and 50 μg ml⁻¹, respectively. The absence of bacteria was verified periodically by plating 100 μl on Difco Marine Agar 2216 (MA; BD Diagnostics) and enumerating bacteria using flow cytometry (Extended Data Fig. 12e). Dilutions with fresh medium were performed to maintain the algal cultures in exponential growth.

For viral infections, 20 ml of exponentially growing, axenic algal cultures at about 8 × 10⁵ cells per ml were infected with EhV-M1 (freshly propagated on RCC6946, concentrated on 100 kDa Amicon filters) at a multiplicity of particles (MOP) of 1:2 of alga:virus. For testing the effect of antibiotics on algal mortality, ampicillin and kanamycin were added to xenic RCC6946 cultures at final concentrations of 100 μg ml⁻¹ and 50 μg ml⁻¹, respectively, simultaneously with viral co-infection at an MOP of 1:5. The experiment was performed in biological quadruplicate and repeated independently twice with similar results. Statistical significance was assessed for the log₁₀-transformed data by using a two-sided linear mixed-effects model, followed by Holm-adjusted pairwise comparisons between treatments at each time point using lmerTest (v3.2-1) and emmeans (v2.0.3) in R (v4.5.1).

Viral infection experiments with additional bloom-forming marine microalgae

Xenic strains of the haptophyte P. globosa G(A), and the green algae M. commoda LAC38, and O. tauri RCC4221 were cultured and used for propagation of the respective viruses Phaeocystis globosa virus (PgV-07T), Micromonas commoda virus (McV-BIO217) and two Ostreococcus tauri viruses (specialist virus OtV5 and generalist virus OtV06-12), as described above, with some modifications to slow down the temporal dynamics of the viral infection and allow for bacterial synergy. Specifically, P. globosa was grown in f/2 medium⁶⁴ and infected at 5 × 10⁵ cells per ml with PgV-07T (filtered at 0.22 μm) at an MOP of 1:0.1. M. commoda was grown in FSW adjusted to 27‰ and supplemented with f/2 medium, and infected at 10⁶ cells per ml with McV-BIO217 (filtered at 0.45 μm) at an MOP of 1:1. O. tauri was grown in K/2 medium at 50 μmol photons m⁻² s⁻¹, and infected: (1) at 10⁷ cells per ml with OtV5 (filtered at 0.45 μm and concentrated on 100 kDa) at an MOP of 1:0.001; and (2) at 10⁶ cells per ml with OtV06-12 (filtered at 0.22 μm and concentrated on 100 kDa) at an MOP of 1:0.1. To test the effect of antibiotics, ampicillin and kanamycin were added at 100 μg ml⁻¹ and 50 μg ml⁻¹ at day −1 (P. globosa), day 0 (O. tauri with OtV5), or at days 1 and 2 (M. commoda and O. tauri with OtV06-12) depending on the microbiome response. All experiments were performed in biological triplicate and repeated independently twice with similar results. Statistical significance was assessed for the log₁₀-transformed data by using a two-sided linear mixed-effects model, followed by Holm-adjusted pairwise comparisons between treatments at each time point.

Enrichment of the BAM on G. huxleyi

To enrich the BAM on G. huxleyi, the cryopreserved BAM was grown on diverse media containing organic material derived through exudation, mechanical lysis (sonication) and viral lysis of algal cells. To collect exponential conditioned media (Exp. CM) and sonication lysates (SL), algal cultures were grown to their exponential phase at about 8 × 10⁵ cells per ml and centrifuged at 600g for 5 min. The supernatant was filtered through 0.22-μm PES filters and used as ‘Exp. CM’. The cell pellet was resuspended in 0.3 ml K/2 medium, physically disrupted through probe sonication (5 cycles of 30 s followed by 30 s rest). Once sonicated, the lysate was resuspended in 500 mL K/2 medium, filtered through 0.22-μm PES filters and used as ‘SL’. To collect viral free lysate (VFL), algal cultures were infected with EhV-M1 as described bove. After four days, the lysed cultures were filtered sequentially through 2-μm cellulose membrane filters, 0.22-μm PES filters and 100 kDa centrifugal filters, and used as ‘VFL’. All three filtrates were kept at 4 °C and used within 24 h. The BAM glycerol stock was directly inoculated into 30 ml of each filtrate in triplicate and incubated for 3 days in an orbital shaker at 18 °C. Triplicates were combined, and aliquots of 9 ml per medium were mixed and cryopreserved with glycerol (15% final concentration).

Tripartite alga–virus–bacteria interaction experiments

Algal cultures were infected with EhV-M1, as described above, and co-inoculated with an overnight culture of BAM. BAM was inoculated from a glycerol stock into Difco Marine Broth 2216 (MB; BD Diagnostics) and incubated overnight at 28 °C with agitation (160 rpm). The overnight culture was washed three times with 0.2-µm-filtered FSW by centrifugation (13,800g, 1 min), enumerated using flow cytometry, and inoculated at 10⁴ bacterial cells per ml into the virus-infected algal cultures. The experiment was performed in n = 5 biological replicates and repeated independently at least twice with similar results. Statistical significance was assessed for the log₁₀-transformed data by using a two-sided linear mixed-effects model, followed by Holm-adjusted pairwise comparisons between treatments at each time point.

Quantification of algae, bacteria and viruses by flow cytometry

Flow cytometry analyses were performed using a CytoFLEX S flow cytometer (Beckman Coulter). Algal cells were identified as high forward scatter (FSC) and high chlorophyll population by plotting FSC (excitation: 488 nm, height threshold = 4 × 10⁴ a.u.) versus chlorophyll autofluorescence (excitation: 561 nm, emission: 665–715 nm). Only cells with a high chlorophyll autofluorescence >4 × 10⁴ a.u. were counted (Extended Data Fig. 12a,b). For bacteria and virus counts, samples were fixed with glutaraldehyde (0.5% final concentration) for at least 30 min at 4 °C, then snap-frozen in liquid nitrogen and stored at −80 °C until analysis. After thawing, samples were stained with SYBR Gold (Invitrogen; diluted 1:10⁴ in Tris-EDTA buffer) for 20 min at 80 °C and analysed by flow cytometry by plotting violet side scatter (vSSC; excitation: 405 nm, area threshold = 1 × 10³ a.u.) versus SYBR fluorescence (excitation: 488 nm, emission: 500-550 nm, area threshold = 1 × 10³ a.u.; Extended Data Fig. 12c–f). All flow cytometry data were inspected and processed in CytExpert v2.4 (Beckman Coulter), technical outliers (due to unstable timelines or shifts in signal intensity) removed, and statistical analyses were conducted in GraphPad Prism (v11.0.2) and R (v4.5.1).

Quantification of virus-infected algal cells using smFISH

Alga–virus–bacteria interaction experiments were performed as described above, with 500 ml per replicate. Cultures were sampled daily for: (1) quantification of algae, bacteria and viruses via flow cytometry as described above; and (2) smFISH analysis, as previously described³. In brief, at each time point, 25–30 ml of culture was fixed with paraformaldehyde (1% final concentration), incubated for 1 h, and centrifuged for 3 min at 3,000g. The cell pellet was incubated with a cryopreservation solution for 1 h, centrifuged for 3 min at 3,000g, and stored at −80 °C. For hybridization, chlorophyll was removed by two sequential washes with 70% and 100% ethanol, enzymes were degraded with proteinase K (10 μg ml⁻¹ final concentration) in 2× sodium chloride-sodium citrate (SSC), and the cells incubated overnight with hybridization buffer (17.5% formamide in 2× SSC) containing probes for the viral mcp transcript conjugated with the fluorophore tetramethylrhodamine (TMR; 0.1 ng ml⁻¹ final concentration). Finally, cells were washed (17.5% formamide in 2× SSC), resuspended in buffer (10 mM Tris at pH 8.0 in 2× SSC), and analysed by flow cytometry (FSC area threshold = 4 × 10³ a.u., >1,000 total events in the target gate at 30 μl min⁻¹, abort rate <7%) using the viral probe’s fluorescence (excitation: 561 nm, emission: 564–606 nm). Cells with a signal area greater than 1 × 10³ a.u. TMR were enumerated as mcp positive (mcp+), virus-infected cells (Extended Data Fig. 12g). Unstained samples were used as a reference for background fluorescence. The experiment was performed in biological quadruplicate and repeated independently twice with similar results. Statistical significance was assessed using a two-sided linear mixed-effects model, followed by Holm-adjusted pairwise comparisons per time point.

Bacterial community analysis using 16S rRNA amplicon sequencing

Alga–virus–bacteria interaction experiments were set up in five replicates at 400 ml. Aliquots of 40 ml were filtered by gravity through autoclaved GF/C filters (47 mm; Whatman), and the filtrate was collected on autoclaved 0.22–µm PVDF filters (47 mm; Durapore, Millipore). The filters were transferred directly to the lysis buffer of DNeasy PowerWater Kit (Qiagen), frozen in liquid nitrogen, and stored at −20 °C. DNA extraction was performed according to the manufacturer’s instructions. The 16S rRNA gene was amplified using the primers 515 F (GTGYCAGCMGCCGCGGTAA⁶⁵) and 806 R (GGACTACNVGGGTWTCTAAT⁶⁶) targeting the V4 region, under the following PCR conditions: initial denaturation at 94 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 45 s, annealing at 50 °C for 1 min, and elongation at 72 °C for 1 min and 30 s. A final elongation step was performed at 72 °C for 10 min to ensure complete extension of the amplified products. Sequencing was performed using a paired-end 2× 250 bp Illumina MiSeq platform at the RUSH Genomics and Microbiome Core Facility.

The 16S rRNA amplicon sequencing data were processed using the DADA2 pipeline⁶⁷ (v1.28.0) in R (v4.3.2) to perform quality filtering, error correction, and amplicon sequence variant (ASV) inference. Raw sequencing reads were quality-checked using the plotQualityProfile function, and low-quality bases were trimmed using the filterAndTrim function with truncation parameters set based on quality score distributions. Error rates were learned from the dataset using the learnErrors function, and sequences were dereplicated using derepFastq. ASVs were inferred using the dada function, followed by chimera removal with removeBimeraDenovo. Taxonomic classification was performed using the SILVA_nr99 reference database⁶⁸ (v138.1), and ASV tables were generated for downstream analyses. The average sequencing depth per sample was 67,413 ± 26,969 (minimum 24,150 and maximum 132,743 reads). Chloroplast reads accounting for 28 ± 23% (maximum 73%) were removed. The processed data were then analysed using the Phyloseq package⁶⁹ (v1.44.0) in R to visualize microbial community composition. The nine most abundant genera are shown individually, while all remaining genera (<5%) were summed as ‘Others’.

Calculation of SI and BI

SI was calculated by comparing the AUC of virus-infected G. huxleyi growth in the presence and absence of bacteria (the entire growth curves were used). We chose to focus on AUC for four main reasons: First, it is biologically meaningful and captures the cumulative biomass change that would occur in natural bloom systems, reflecting the total ‘ecological impact’ (for example, the total amount of carbon fixed). Second, it captures complete temporal dynamics rather than relying on a single time point or summary statistic that would miss critical dynamics. Third, it is more robust against experimental variability, as it takes into account three measurements from each experiment day. Finally, normalizing AUC difference by another AUC creates a dimensionless metric that allows direct comparison across different bacterial strains and experimental conditions even if their baseline levels are different. SI was expressed as a percentage and calculated as:

where: AUC_virus represents the algal population dynamics under viral infection alone, and AUC_tripartite represents the algal population dynamics under co-infection with both viruses and bacteria.

We computed the AUC of the algal abundance for each replicate and treatment using GraphPad Prism (v11.0.2). Then, the mean AUC of each treatment (n = 3) was used to compute the SI as described above. The standard error (SE) of the SI was estimated as \({\rm{S}}{\rm{E}}\,({\rm{ \% }})=\frac{\sqrt{{\rm{V}}{\rm{a}}{\rm{r}}({\rm{S}}{\rm{I}})}}{\sqrt{{n}}}\times 100\) using analytical error propagation via the delta method, whereby the variance of the SI was determined as \(\mathrm{Var}(\mathrm{SI})=\frac{1}{{{\mathrm{AUC}}_{\mathrm{virus}}}^{2}}\,\times \) \({\rm{V}}{\rm{a}}{\rm{r}}({{\rm{A}}{\rm{U}}{\rm{C}}}_{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{p}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{i}}{\rm{t}}{\rm{e}}})+(\frac{{{\rm{A}}{\rm{U}}{\rm{C}}}_{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{p}}{\rm{a}}{\rm{r}}{\rm{t}}{\rm{i}}{\rm{t}}{\rm{e}}}}{{{{\rm{A}}{\rm{U}}{\rm{C}}}_{{\rm{v}}{\rm{i}}{\rm{r}}{\rm{u}}{\rm{s}}}}^{2}}{)}^{2}\times {\rm{V}}{\rm{a}}{\rm{r}}({{\rm{A}}{\rm{U}}{\rm{C}}}_{{\rm{v}}{\rm{i}}{\rm{r}}{\rm{u}}{\rm{s}}})\). To test for a significant treatment difference between the tripartite and virus-infected demise dynamics and to label the SI values accordingly, a repeated-measures two-way ANOVA was performed at a two-tailed significance level of α = 0.05 to evaluate the effects of time and the bacterial strain on algal abundance. To account for cell calcification and thus sedimentation, the raw data were detrended with a linear mixed-effects model (fixed slope, random intercept per day) using the robustlmm package (v3.4-3) in R (v4.5.1). Results are reported as SI ± SE.

BI quantifies bacterial growth differences over time by comparing the AUC in tripartite conditions versus bacterial growth alone. BI was expressed as a percentage and calculated as:

where AUC_tripartite represents the bacterial growth over time when co-cultured with virus-infected G. huxleyi, and AUC_bacteria represents bacterial growth over time in the presence of G. huxleyi but without viral infection.

We first computed the AUC of the bacterial abundance for each replicate and treatment between day 1 and day 6 (at day 0, bacterial inoculation density was below the detection limit of the flow cytometry acquisition settings). Then, the mean AUC of each treatment (n = 3) was used to compute the BI as described above. The SE of the BI was estimated as \(\mathrm{SE}\,( \% )=\frac{\sqrt{\mathrm{Var}(\mathrm{BI})}}{\sqrt{n}}\times 100\) using analytical error propagation via the delta method, whereby the variance of the BI was determined as \(\mathrm{Var}(\mathrm{BI})=\frac{1}{{{\mathrm{AUC}}_{\mathrm{bacteria}}}^{2}}\times \mathrm{Var}({\mathrm{AUC}}_{\mathrm{tripartite}})+{\left(\frac{{\mathrm{AUC}}_{\mathrm{tripartite}}}{{{\mathrm{AUC}}_{\mathrm{bacteria}}}^{2}}\right)}^{2}\times \) \(\mathrm{Var}({\mathrm{AUC}}_{\mathrm{bacteria}})\). To test for significant differences between tripartite and bacterial growth dynamics and to label the BI values accordingly, a repeated-measures two-way ANOVA was performed at a two-tailed significance level of α = 0.05 to evaluate the effects of time and the bacterial strain on bacterial abundance (in case of missing values, a mixed-effects model was fitted). Results are reported as BI ± SE.

Isolation of bacteria from BAM and induced G. huxleyi blooms

For bacterial isolation, glycerol stocks were sampled on different days of the induced G. huxleyi blooms and following various enrichments (Supplementary Data 1). In brief, seawater samples (filtered using 10 µm and GF/F filters) were collected during the bloom (day 13) and demise (day 21 and day 23) phases and cryopreserved as glycerol stocks at a final concentration of 15%, snap-frozen in liquid nitrogen, and stored at −80 °C. The BAM glycerol stock was prepared as described above. To enrich alga- and viral infection-associated taxa, the mesocosm and BAM glycerol stocks were inoculated into uninfected and virus-infected cultures of G. huxleyi before isolation. To enrich G. huxleyi-associated taxa, xenic cultures of RCC6946, which contain bacteria that were co-isolated with the algal cell from the mesocosm bloom, were infected with EhV-M1 before isolation. From overnight cultures of the glycerol stocks and enrichment cultures, 100 µl were plated at 100×, 1,000×, and 10,000× dilution on MA plates and incubated at 28 °C, allowing for single-colony isolation. Selected colonies were streaked onto fresh MA three times or until a uniform colony morphology was observed, ensuring bacterial isolate purity. Colonies were inoculated in 3 ml MB medium, grown overnight at 28 °C with agitation (160 rpm), and cryopreserved as glycerol stocks. Tripartite alga–virus–bacteria experiments were performed as described above in biological triplicate and repeated independently twice with similar results.

For sequencing the complete 16S rRNA gene, cells from pure colonies were lysed using a buffer containing TE and 0.1% Triton-X100, followed by heating at 95 °C for 5 min. Cell lysates were centrifuged at 13,300 rpm for 10 min, and 2 µl of the supernatant was used as a DNA template for PCR amplification. The 16S rRNA gene was amplified using the universal primers 16S-27F (AGAGTTTGATCMTGGCTCAG)⁷⁰ and 16S-1492R (CGGTTACCTTGTTACGACTT)⁷¹. The PCR reaction was carried out in 20 µl using 2× HY-Taq (HyLabs) or REDTaq (Sigma) Ready Mix. PCR amplification was performed using the following conditions: initial denaturation at 94 °C for 3 min, followed by 25–35 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C for 30 s, and elongation at 72 °C. A final elongation step was conducted at 72 °C for 10 min, followed by a hold at 12 °C. PCR products were visualized in 1% w/v agarose gel with ethidium bromide. Amplicons were purified with the illustra ExoProStar 1-Step enzymatic cleanup kit (GE Healthcare Life Sciences) using a 5 µl aliquot of the PCR reaction and incubated at 37 °C for 15 min, followed by enzyme inactivation at 80 °C for 15 min. For Sanger sequencing, 1–2 µl of the purified amplicons were analysed by the Life Sciences Core Facilities at the Weizmann Institute of Science. All 16S rRNA gene sequences were classified at the genus level using the SINA classifier tool⁷² (v1.2.12) within the SILVA rRNA database (Supplementary Data 1).

Quantification of algal cell death using SYTOX staining

Alga–virus–bacteria interaction experiments were performed as described above. Cultures were sampled daily for: (1) enumeration of algae, bacteria and viruses via flow cytometry as described above; and (2) quantification of dying (membrane-permeable cells) using the DNA-intercalating stain SYTOX Green (Invitrogen). Samples were stained with a final concentration of 1 μM SYTOX Green, incubated for at least 20 min in the dark at room temperature, and analysed by flow cytometry (FSC height threshold = 2.5 × 10⁴ a.u., vSSC area threshold = 1 × 10³ a.u.) using the SYTOX Green fluorescence (excitation: 488 nm; emission: 505–545 nm). For samples with ≥400 cells (parent gate), cells with a signal area >1 × 10⁵ a.u. SYTOX Green (Extended Data Fig. 12h) were counted as SYTOX⁺, membrane-permeable cells. Cells with a chlorophyll autofluorescence >2 × 10⁵ a.u. were enumerated as high-chlorophyll cells. Unstained samples were used as a reference for background fluorescence. The experiment was performed in biological quadruplicate and repeated independently twice with similar results. Statistical significance was assessed using a two-sided linear mixed-effects model, followed by Holm-adjusted pairwise comparisons per time point.

Genomic sequencing of bacterial isolates and metabolic pathway analysis

Bacterial isolates (Supplementary Data 1) were inoculated from glycerol stocks in MB medium and grown overnight at 28 °C with agitation (160 rpm). Cultures were centrifuged at 12,000 rpm at room temperature for 1 min, the cell pellet was washed 3× with FSW, and frozen in liquid nitrogen. Genomic DNA was extracted following a modified version⁷³ of a previously described protocol⁷⁴. High molecular weight of the genomic DNA was verified by gel electrophoresis or Genomic DNA ScreenTape Analysis using TapeStation 4150 (Agilent Technologies), while DNA concentration and purity were validated using a Qubit fluorometer with the dsDNA HS Assay Kit (Thermo Fisher Scientific) and NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific), respectively. If necessary, samples were concentrated using the NEBNext Ultra II FS DNA Library Prep Kit (New England Biolabs) or AMPure XP beads (Beckman Coulter).

For Illumina short-read sequencing, shotgun metagenome libraries were prepared using the NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New England Biolabs) with enzymatic shearing, according to the manufacturer’s instructions. The final size-selected pool was sequenced on an Illumina NovaSeq 6000 instrument with an S2 flow cell, utilizing 2× 150 base reads. Short reads were quality trimmed using Trimmomatic (v0.39) with the following parameters: ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10, LEADING:10, TRAILING:10, SLIDINGWINDOW:4:15, MINLEN:50. Reads that passed quality trimming were used for assembly construction with SPAdes (v3.13.0), using default k-mer sizes (21 bp, 33 bp, 55 bp and 77 bp). Contigs were filtered based on mean coverage (≥3×) and length (≥300 bp) and clustered into bins using MaxBin2 (v2.2.7) with the default parameters. Assembly quality was assessed using CheckM2 (v1.0.2) with default parameters. Multiple bins of the same isolate were concatenated into single genomes if they complemented one another based on completeness and had similar taxonomic annotation. Genes were predicted from the genomes using Prodigal (v2.6.3) in metagenomic mode. Predicted genes were filtered based on length (≥300 bp), the presence of a start codon at the first trimer, and the absence of a premature stop codon. Genes that passed filtration were clustered using MMseqs2 (v13.45111) with the following parameters: coverage mode 0, minimum coverage 0.8, and minimum sequence identity 0.9. Representative genes were annotated using eggNOG-mapper (v2.1.3) and used for downstream analysis. For PacBio long-read sequencing, multiplexed microbial libraries were prepared, using SMRTBell Express Prep Kit 2.0 (Pacific Biosciences). A single PacBio sequencing library was constructed that was used for SMRTbell template preparation as described in the protocol (Pacific Biosciences). The resulting SMRTbell template was annealed to the sequencing polymerase, using a Sequel binding kit 2.0. Sequencing was conducted with a 12pM on-plate concentration using a Sequel I System in CLR mode and a single SMRT Cell (1M; PacBio) at the WIS Genomics Unit. PacBio raw reads (polymerase reads) were processed using the SMRT Link software (v8.0). Genomes were assembled using Trycycler (v0.5.5) and Flye (v2.9.6) with default parameters. Assembly quality was assessed using CheckM2, and genes were predicted, filtered, clustered, and representative genes annotated, as described above. Nanopore sequencing was performed by Plasmidsaurus following their standard bacterial genome sequencing service using the Oxford Nanopore Technology with custom analysis and annotation. All assembled genome sequences were classified at the species level using the GTDB-Tk software toolkit (v2.4.1) with GTDB database⁷⁵ v220 (Supplementary Data 1). The genome of Sulfitobacter sp. D7 was sequenced as previously described⁷⁶.

To assess the relationship between bacterial metabolic pathways and bacterial synergy, all KEGG pathways and their constituent KEGG orthologues were retrieved via KEGG API. Using the eggNOG-based genome annotations for each taxon, pathway coverage (%) was calculated as the proportion of KEGG orthologues present in a strain relative to the total KEGG orthologues in that pathway. Pathway coverage for the FAD pathway was then correlated with the SI using a Spearman’s rank correlation, and the correlation coefficient ρ with the respective P value determined.

Enumeration of algae and viruses in induced G. huxleyi blooms

The abundance of G. huxleyi and intracellular EhV throughout the succession of the induced G. huxleyi bloom in mesocosm bag 7 were enumerated as previously described²⁴. In brief, for the enumeration of G. huxleyi cells, water samples were pre-filtered using 40-µm cell strainers and immediately analysed with an Eclipse iCyt flow cytometer using ec800 software v1.3.7 (Sony Biotechnology). Calcified G. huxleyi cells were identified as a population of cells with high SSC and high chlorophyll (excitation: 488 nm, emission: 663–737 nm). For the enumeration of intracellular EhV, water samples were pre-filtered through 20-µm filters, and DNA was extracted from 2-µm filters. EhV abundance was determined by quantitative PCR using the major capsid protein (mcp) gene. Absolute abundances were calculated using a calibration curve of DNA at known concentrations.

16S rRNA amplicon sequencing of bacteria in induced G. huxleyi blooms

Planktonic (>0.2 µm and <2 µm) and particle-attached (>20 µm and <200 µm) fractions of the bacterial community were assessed throughout the succession of the induced G. huxleyi bloom in mesocosm bag 7 using 16S rRNA amplicon sequencing, as previously described²⁴. In brief, water samples were pre-filtered through 200 µm, sequentially filtered through 20-µm, 2-µm and 0.22-µm polycarbonate filters (Isopore, 47 mm; Merck Millipore), and immediately frozen in liquid nitrogen. DNA extraction of the 20-µm and 0.22-µm filters was performed with the DNeasy PowerWater kit (Qiagen), and sequenced following the EMP 16S amplicon protocol targeting the V4 region of the 16S gene using the 515F-806R primers⁷⁷ including degeneracy to reduce bias towards Thaumarchaeota and SAR11^65,66. The particle-attached fraction was sequenced using a paired-end 2× 153 bp Illumina MiniSeq platform at the RUSH Genomics and Microbiome Core Facility. Data analysis for planktonic bacteria was performed using the DADA2 pipeline and ASVs annotated with the RDP classifier. The average sequencing depth per sample and chloroplast read proportion for the planktonic community was previously reported²⁴. Data analysis for particle-attached bacteria was performed using the DADA2 pipeline (setting ‘justConcatenate = TRUE’ to merge the forward and reverse reads introducing 10x ‘N’ in the middle) and ASVs annotated with the SILVA_nr99 reference database. The average sequencing depth for the particle-attached size fraction in bag 7 was 20,734 ± 6,517 (minimum 8,057 and maximum 29,629) reads per sample. Chloroplast reads accounting for 7 ± 8% (minimum 1% and maximum 36%) in the particle-attached fraction were removed. Kruskal–Wallis tests were performed to compare planktonic and particle-attached bacterial fractions for each family and bloom phase (pre-bloom, days 1–6; bloom, days 7–14; and viral infection phase, days 15–23). Pairwise Welch’s t-tests with the Holm–Bonferroni correction for multiple comparisons were performed comparing the log₂-transformed fold change of particle-attached planktonic bacterial fractions between the bloom phases. For missing samples (day 2 of the planktonic and day 1 of the particle-attached size fraction), the relative abundances were interpolated. Absence of Pseudoalteromonadaceae in some planktonic samples prevented log₂ fold change computation.

Chemotaxis assay and 16S rRNA amplicon sequencing of chemotactic bacteria in induced G. huxleyi blooms

To prepare cell extracts as chemoattractants, G. huxleyi CCMP 2090 was grown in 19 l modified K/2 medium in 2 carboys with stirring and aeration. One exponentially growing culture (1.2 × 10⁶ cells per ml) was infected with EhV 201 at an MOP of 1:1 while the other culture served as an uninfected control. The viruses were propagated on CCMP 2090, washed and concentrated on 100 kDa Amicon-ultra filters (EMD Millipore), and stored at 4 °C until use within two months. Cells were counted daily using a Multisizer 4 Coulter Counter (Beckman Coulter). After 48 h, the uninfected (4.4 × 10⁶ cells per ml) and infected (0.9 × 10⁶ cells per ml) cultures were filtered onto pre-combusted GF/C filters (47 mm; Whatman) using gentle vacuum (600 mbar under-pressure). Filters were immediately transferred to 500 ml pre-cooled (−20 °C) methanol (HPLC grade, J.T. Baker) and extracted overnight at 4 °C while shaking. Extracts were sonicated for 30 min and filtered through pre-combusted GF/C filters before being evaporated to 20 ml using a Rotovap (30 °C, 70–100 rpm). Extracts were centrifuged at 20,800g for 10 min at 4 °C, and the supernatants were evaporated to complete dryness using a SpeedVac (28 °C, 1,980 rpm). Dried extracts were re-dissolved in methanol to reduce the salt load, aliquoted into eight centrifuge tubes, dried using a SpeedVac, and stored at −80 °C until further use. The dry weight was 57.5 ± 0.3 mg and 38.3 ± 0.5 mg for the uninfected and infected extracts, respectively.

In situ chemotaxis assays (ISCAs) containing the G. huxleyi cell extracts were deployed five times throughout the mesocosm bloom succession, namely at days 3, 9, 12, 18 and 20. For each deployment, 12.5 ml of bloom water were collected per mesocosm bag (1–4; 50 ml in total), mixed and subjected to a quadruple-filtration process. Specifically, the mixed bloom water was first filtered through a 0.2-μm Millex FG (Millipore) to remove large particles, twice through a 0.2-μm Sterivex filter (Millipore), and finally through a 0.02-μm Anatop filter (Whatman). This quadruple filtration aimed to remove all microorganisms from the seawater. We confirmed the efficiency of the filtration procedure using flow cytometry and DNA sequencing⁷⁸. The filtered bloom water was used as a control in the ISCA and was also used to resuspend the dried G. huxleyi extracts at a final concentration of 1 mg ml⁻¹. The resuspended extracts were filtered through a 0.2-μm Millex FG (Millipore) to remove undissolved particles.

The working concentration of G. huxleyi extracts (1 mg ml⁻¹ of dry weight) is representative of environmental hotspots. Indeed, chemical concentrations decay exponentially with distance away from the ISCA opening (port)⁴¹. This means that concentrations sensed by bacteria in the surrounding bloom water will be substantially lower than those in the ISCA wells and within the same order of magnitude as the concentrations arising from the lysis of a phytoplankton cell⁷⁹.

Treatments (filtered bloom water and the two G. huxleyi cell extracts) were randomly allocated to an ISCA row (consisting of 5 wells with 110 μl each). All wells in a row contained the same treatment, and each treatment was replicated on three different ISCA devices, which were deployed in parallel to act as biological replicates (n = 3). Each ISCA was placed into a flow-damping acrylic enclosure, following a pre-established protocol⁸⁰. Each enclosure was filled with water coming from the first four mesocosm bags (bags 1–4), previously mixed at equal volume. The flow-damping enclosures were filled slowly with water from the inlet in the lower surface and were sealed with a plug once full. The enclosures were then incubated in situ at 1 m depth for 2 h (ref. ⁷⁸). Following deployments, enclosures were retrieved and slowly drained. The volume of each ISCA well was retrieved using a sterile 1 ml syringe and 27G needle.

ISCA samples for DNA extraction were snap-frozen immediately. DNA extraction was performed using a microvolume extraction approach (physical lysis method⁸¹). Libraries for 16S rRNA gene sequencing were generated for all samples using the Nextera XT DNA Sample Preparation Kit (Illumina) following a previously described protocol for low-input DNA libraries⁸². Primers targeting the V4–V5 regions of the 16S rRNA gene were used, namely 515F and 926R (CCGYCAATTYMTTTRAGTTT⁶⁵). Samples were then sequenced on an Illumina MiSeq platform with a paired-end 2× 300 bp output. Libraries were pooled on an indexed sequencing run and yielded ~300 bp per sample. The composition of the chemotactic bacterial communities was analysed as described above using the DADA2 pipeline, and ASVs were annotated with the SILVA_nr99 reference database. For day 18, a Welch’s ANOVA followed by two-sided pairwise t-tests with Benjamini–Hochberg correction were performed to compare the fraction of Alteromonadaceae across chemoattractants using the R ‘stats’ package (v4.5.1; Shapiro–Wilk test and Levene test were performed to assure normality and homogeneity of variances).

Isolation and chemotactic specificity of bloom-associated bacteria

To isolate chemotactic bacteria throughout the succession of the induced G. huxleyi bloom, we prepared agar plates made with 0.22-μm-filtered environmental bloom water and amended with 1 mg ml⁻¹ of the corresponding G. huxleyi cell extracts. The agar plates were inoculated with the chemotactic bacterial community from the ISCA wells containing uninfected and virus-infected G. huxleyi extracts and incubated at fjord water temperature (10 °C) until colonies appeared. Individual colonies were then grown in 10% MB in autoclaved seawater overnight at room temperature and at 10 °C. Each isolate went through two additional rounds of purification via successive growth on plates and in liquid cultures. Isolates were identified by 16S rRNA gene Sanger sequencing⁷⁸. Selected strains (Supplementary Data 1) were further used for genome sequencing and characterization of their chemotaxis towards different G. huxleyi cell extracts.

For genome sequencing, overnight cultures were inoculated in MB medium at 0.1 OD₆₀₀. Exponentially growing cultures at 1.0 OD₆₀₀ were centrifuged, and cell pellets of 15–30 mg were re-dissolved in 0.5 ml Zymo DNA/RNA Shield. Genomic DNA extraction and Nanopore sequencing were performed by Plasmidsaurus following their standard bacterial genome sequencing service. Assembled genome sequences were classified at the species level using the GTDB-Tk software toolkit.

For laboratory-based chemotaxis assays, the bacteria were plated from cryostocks on MA. Multiple colonies from each strain were inoculated in 1% MB in artificial seawater (33 g l⁻¹ Instant Ocean, Aquarium Systems) and incubated at 180 rpm and 27 °C for 4 h to promote bacterial motility. An aliquot was fixed with glutaraldehyde (1% final concentration; Sigma Aldrich), stained for 15 min with SYBR Green I (ThermoFisher), and bacterial abundance enumerated using a CytoFLEX S flow cytometer (Beckman Coulter) as described above. Chemotaxis was characterized using ISCA devices that were each loaded with four different 0.2-μm-filtered treatments and incubated at room temperature for 1 h in artificial seawater inoculated with 10⁶ bacteria per ml. As chemoattractants, we used: (1) extract from uninfected G. huxleyi CCMP 2090 cultures (1 mg ml⁻¹ in artificial seawater); (2) extract from virus-infected G. huxleyi CCMP 2090 cultures (1 mg ml⁻¹ in artificial seawater); (3) 10% MB (acting as a positive control), and artificial seawater (acting as a negative control to account for cells entering wells by random motility). Chemoattractants were inoculated in ISCA wells using a sterile 1 ml syringe (Codan Medical ApS) and needle (27 G, BD Microlance). In each ISCA, each row of five wells served as technical replicates for an individual treatment. Each chemoattractant was deployed simultaneously in four ISCA replicates. After incubation, contents of the ISCA wells were recovered using a sterile syringe and needle, technical replicates were pooled⁸⁰, and 200 µl per sample were fixed, stained and enumerated using flow cytometry. Chemotactic indices (I_c) were determined as ratios of the cell concentrations in treatment wells relative to the average artificial seawater control. For statistical analysis, I_c values were log1p-transformed (ln(1 + x)) and the Levene test of homogeneity of variances and Shapiro–Wilk test of normality were performed using R. All bacteria showed significant differences between the artificial seawater control and the treatments, including the MB control (Dunnett’s test with P < 0.001). Differences between the two cell extracts were assessed using a two-sided Welch’s t-test with the Holm–Bonferroni correction for multiple comparisons. For fold change comparison, the transformed I_c values were further log₂-transformed, the normality per strain and treatment validated using a Shapiro–Wilk test, and treatment differences per strain assessed using a two-sided Welch’s t-test with Benjamini–Hochberg correction for multiple comparisons. The log₂-transformed fold changes (LFC) between the infected and control extracts were computed as LFC = mean(log(Inf)) − mean(log(Ctl)) with s.e.m. = √ (Var_Inf/n_Inf + Var_Ctl/n_Ctl).

16S rRNA amplicon sequencing of coccolithophore blooms in the North Atlantic

Coccolithophore blooms were sampled in the North Atlantic during the NA-VICE cruise in 2012 with large spatial coverage¹⁸ and the Tara Breizh cruise in 2019 with high temporal resolution⁴³.

During the NA-VICE cruise, coccolithophore blooms across a 2,000-mile transect between the Azores and Iceland were identified via satellite remote sensing and sampled at 30 stations onboard the RV Knorr between 13 June and 16 July 2012, as previously described¹⁸. Samples for DNA extraction and 16S rRNA amplicon sequencing were collected between the surface and 150 m depth with Niskin bottles, and filtered through a 200-μm prefilter onto 0.8-µm polycarbonate filters. Genomic DNA was extracted and amplified in a dual PCR strategy using primers targeting the V1–V3 region of the 16S rRNA gene⁸³ linked to adaptors, namely CS1_27F_YM (ACACTGACGACATGGTTCTACA-AGAGTTTGATYMTGGCTCAG) and CS2_534R_Y (TACGGTAGCAGAGACTTGGTCT-TYACCGCGGCTGCTGG), and subsequently linked to unique sample barcodes for multiplexing. For PCR amplification, the AccuPrime SuperMix II (Life Technologies) was used using the following cycling conditions: denaturation for 5 min at 95 °C, 28 cycles of 30 s at 95 °C, 30 s at 56 °C and 5 min at 68 °C, and final 7 min elongation at 68 °C. The PCR amplicon quality was validated by gel electrophoresis. In a second PCR step, sample barcodes were amplified with the following cycling conditions: 5 min at 95 °C, 8 cycles of 30 s at 95 °C, 30 s at 60 °C and 30 s at 68 °C, and final 7 min elongation at 68 °C. Samples were pooled and purified using AMPure XP beads, and concentrations were estimated using Qubit. Sequencing was performed using a paired-end 2× 300 bp Illumina MiSeq technology at the DNA Services Facility, University of Illinois at Chicago.

During the Tara Breizh cruise, a G. huxleyi bloom patch in the Celtic Sea was identified by satellite remote sensing and sampled onboard the SV Tara between 27 May and 2 June 2019, as previously described^21,43. Samples for DNA extraction and sequencing were collected at the maximum turbidity depth (bloom depth) by pumping seawater through a 20-μm prefilter and sequential filtration onto 3-μm and 0.2-μm PC filters. Following DNA extraction, amplicon sequencing of the V4–V5 region of the 16S rRNA gene was performed twice using Illumina MiSeq paired-end 2× 250 bp sequencing technology. Fastq files were obtained from the European Nucleotide Archive (ENA: PRJEB50692). Intracellular virus abundances were quantified from DNA extracts of 3-µm filters using droplet digital PCR and specific primers for the mcp gene of EhV. Samples for flow cytometry were collected at the bloom depth with a Niskin bottle and immediately analysed using a BD Accuri C6 Plus Personal Flow Cytometer (BD Biosciences). Cell abundances of calcified G. huxleyi were determined by gating high chlorophyll (excitation: 488 nm; emission: 663–737 nm) and high SSC (488 nm) particles.

All paired-end 16S rRNA gene amplicon sequences were processed using QIIME 2⁸⁴. ASVs were inferred using the DADA2 paired-end denoising workflow, with trimming and truncation parameters adjusted according to amplicon type. For the NA-VICE samples (V1–V3 amplicons), 20 bp were trimmed from the 5′ ends of both forward and reverse reads, leaving 280 bp (forward) and 230 bp (reverse). Only samples that retained at least 10,000 non-chimeric reads after denoising were considered for downstream analysis. Taxonomic assignment was performed on representative sequences from DADA2, using the QIIME 2 classify-sklearn classifier with the SILVA v138 reference database. ASV tables were aggregated to the genus level, and relative abundances were calculated from the full dataset. The Tara Breizh data (V4–V5 amplicons) were analysed using QIIME 2, as previously described⁴³ with the following modifications. After removal of empty reads, the demultiplexed read pairs were repaired using SeqKit2⁸⁵ (v2.3.1). Primer and orientation handling was performed with CutAdapt (v4.2) in two branches per sample, where each branch represents a primer-orientation-specific read set corresponding to one of the two detected amplicon orientations: one in the original R1/R2 orientation and one in the swapped R2/R1 orientation. Reads lacking the primer sequences, that were <75 bp, and had insertions or deletions were removed. For each sequencing run and orientation branch, denoising was done using QIIME 2 DADA2 with pseudo-pooling. Reads were truncated to 215 bp (forward) and 190 bp (reverse) for the first branch, and to 190 bp (forward) and 215 bp (reverse) for the second branch. Before merging, sequences of the second branch were reverse-complemented. ASV identifiers were normalized by sequence content and feature tables merged by summing identical ASV sequences across the sequencing runs and branches. Chimeras were removed using the de novo UCHIME algorithm of vsearch⁸⁶ (v2.22.1) and non-chimeric ASV counts of the sequencing duplicates were summed. Counts were summed at the ASV level before conversion to relative abundance. Taxonomy was assigned in QIIME 2 using the SILVA v138 99% Naive Bayes classifier. ASV tables were aggregated to the genus level and relative abundance computed. For Fig. 5, a pseudocount of 1 × 10⁻⁵ was assigned to samples with relative abundance of exactly ‘0’. A map of the sample locations was created using the cartopy (v0.25.0) package in python (v3.10) with the Natural Earth map data.

Comparisons of bacterial relative abundance distributions between oceanic datasets

To compare relative bacterial abundances from the North Atlantic coccolithophore blooms with diverse oceanic environments, taxonomic profiles from Tara Oceans⁴⁴ were obtained (https://www.ocean-microbiome.org/). Per sample, relative abundance at the genus and family levels was calculated by dividing read counts for each taxon by the total bacterial and archaeal reads in the corresponding sample. Only samples representing the planktonic bacterial fraction (0.2 to 3 µm and 0.2 to 1.6 µm filter size fraction) were selected. NA-VICE samples were filtered by collection method to keep only cast samples and omit trap samples, while for Tara Breizh all samples labelled as outside the bloom were excluded. Pairwise comparisons of relative abundance distributions were performed at matched taxonomic level and depth range. For the Tara Oceans and NA-VICE comparison, we compared the Alteromonas genus from surface to 150 m depth. Due to limited sequencing depth, which affected our ability to reliably detect Alteromonas, the Tara Oceans and Tara Breizh comparison was performed for the Alteromonadaceae family from surface to 25 m depth. Differences between groups were evaluated using a two-sided Mann–Whitney U test. For Extended Data Fig. 11, a pseudocount of 1 × 10⁻⁵ was assigned to samples with relative abundance of exactly ‘0’. A map of the sample locations was created using the cartopy (v0.25.0) package in Python (v3.10) with the Natural Earth map data.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.