Microbial and chemical dynamics of a toxic dinoflagellate bloom

Sampling sites, filtration, and preservation

Over two days in June 2018, we sampled twelve shoreline sites near the mouth of Tampa Bay, Florida for phytoplankton counts, microbial DNA, and chemical metabolites (Fig. 1, Table S1). This time period marked the beginning of a new surge in K. brevis abundances after 8–10 weeks of low (<100 cells/mL) cell counts along Sarasota and Manatee County coastlines (Fig. S1), although the bloom was technically still in progress statewide both before and after our sampling. Six sites were north of Tampa Bay; six were south of Tampa Bay, and in or around Sarasota Bay. The sites represented two levels of a Karenia brevis bloom based on Florida Fish and Wildlife Conservation Commission (FWC) reports at the time of sampling (https://myfwc.com/research/redtide/statewide/). Surface currents and salinity distributions were obtained from the Finite Volume Community Ocean Models made available by the Ocean Circulation Group at the University of South Florida (http://ocgweb.marine.usf.edu/Models/FVCOM/month_tb_ch.php?year1=2018&mon1=jun).

Sites included docks within sheltered bays, ocean-facing beaches with high wave action, and semi-sheltered sites near bay entrances (Table S1). At each site, approximately 3 L of surface water was collected by hand from immediately offshore, either by wading from the beach into the surf zone (1–2 m depth) or from a dock. Water was stored in autoclaved polystyrene screw-top bottles at room temperature in the dark for 1–5 h before transportation to Mote Marine Laboratory for processing. Additionally, at each site two replicate samples of 20 mL unfiltered water were preserved in glass scintillation vials with 200 µL acidic Lugol’s solution (5 g I₂, 10 g KI, 100 mL DI water, 10 mL glacial acetic acid; final iodine concentration 1%) for microscopy and stored in the dark at room temperature for 2–4 days until processing at Georgia Tech.

Each 3 L sample was pre-filtered through a 30 µm pore size filter to remove any detritus or larger organisms. From this filtrate, 200 mL was processed for metabolomic analysis by sequential filtration onto a 3.0 µm nitrocellulose filter (MilliporeSigma) and a 0.2 µm polycarbonate filter (MilliporeSigma). Both filters were 47 mm in diameter and were stored in glass scintillation vials with methanol (enough to submerge filter) at room temperature until transport to Georgia Tech 1–2 days after sampling, where they were then stored at 4 °C until processing four weeks later.

The remaining ∼2800 mL of filtrate was processed for molecular analysis by sequential filtration onto 3.0 µm nitrocellulose (MilliporeSigma) and 0.22 µm Sterivex (MilliporeSigma) filters. The 3.0 um filters were 47 mm in diameter and housed in Swinnex (MilliporeSigma) filter holders. Filters were preserved with approximately three mL of DNA/RNA stabilization buffer (25 mM sodium citrate, 10 mM EDTA, 5.3 M ammonium sulfate (pH 5.2); (Patin et al., 2018) in five mL CryoTubes (ThermoFisher) and stored at −80 ° C or on dry ice until further processing at Georgia Tech.

Storage and handling time were minimized as much as possible while balancing the importance of treating all samples identically. It is possible that the samples underwent biological or chemical changes (bottle effects) during the time (1–5 h) between water collection and filtration. This is a common trade-off faced by researchers in field studies (Scofield et al., 2015). Bottle effects, if they occurred, would most likely result in a type II (false negative) error, assuming they operated relatively uniformly across samples and were larger than differences due to bloom state.

Phytoplankton cell counts

From each 20 mL sample preserved with Lugol’s solution, the settled cells were enumerated using a 10 mL settling chamber (Hydro-Bios). Cells were settled overnight before counting with an Olympus IX-50 inverted microscope (400X magnification). Phytoplankton cells were identified and quantified under the following categories: diatom single cells, diatom colonies, K. brevis, and other (non-K.brevis) dinoflagellates.

DNA extraction, PCR amplification, and sequencing

DNA was extracted from the filters according to the protocol in (Patin et al., 2018). Briefly, the preservation buffer was flushed from the Sterivex cartridges and replaced with two mL DNA lysis buffer and 40 µL of lysozyme. The membranes (3 and 30 µm pore size) were similarly rinsed and combined with five mL lysis buffer and lysozyme in 15-mL Falcon tubes. Filters were incubated and rotated at 37 °C for 45 min, after which 100 µL of proteinase K and 100 µL of SDS (20%) were added, followed by a 2-hr incubation and rotation at 55 °C. Sterivex cartridges were then flushed with approximately three mL of lysis buffer, and the combined elution volume was transferred to 15 mL centrifuge tubes. Membrane filters were removed from the Falcon tubes, and the remaining solution was used for extraction. Three mL of phenol:chloroform:isoamyl alcohol (25:24:1) were added before vortexing and centrifugation for 5 min at 3500xG. Following centrifugation, the aqueous layer was transferred to a new 15 mL centrifuge tube. Three mL of chloroform:isoamyl alcohol were added and the samples were vortexed. The samples were centrifuged at 2500xG for 5 min. For each sample, the aqueous layer was transferred to an Amicon Ultra-4 centrifuge filter (10 kDa pore size). These were centrifuged at 3500xG until the volume was under 250 µL. TE buffer ( four mL) was added as a rinse and the samples were again centrifuged at 3500xG until the final volume was below 100 µL.

Purified DNA was quantified on a Qubit fluorometer (Life Technologies) and all samples were standardized to a concentration of 4 ng/µL. Sequencing libraries were generated by PCR amplification of the V4 region of the 16S rRNA gene using Earth Microbiome Project primers 515F (Parada, Needham & Fuhrman, 2016) and 806R (Apprill et al., 2015) in a protocol adopted from Kozich et al. (2013). Each PCR contained 21 µL Platinum Taq HiFi Mix (Invitrogen), 1 µL BSA, 0.5 µL each of the forward and reverse primers (10 µM), and 2 µL template DNA. Thermal cycling consisted of initial denaturation for 1 min at 94 °C, followed by 24 cycles of 30 s at 94 °C, 30 s at 55  °C, and 60 s at 68 °C, and a final step of 60 s at 68 °C. Due to poor amplification, DNA from the 30 µm filter samples from sites 2, 3, 4, and 6 was subjected to 30 PCR cycles using a modified protocol involving a 45 s annealing step at 55 °C and a 90 s extension step at 68 °C. Successful amplification was confirmed by agarose gel electrophoresis. We failed to amplify DNA from the largest (>30 µm) size fraction from sites 2, 3, and 4. Amplicons from all other samples were combined in equimolar concentrations, and the combined library was purified using AMPure XP beads (Beckman Coulter) and sequenced on an Illumina MiSeq using a v2 500 cycle kit and a 2 x 250 bp paired-end protocol according to (Patin et al., 2019) and the manufacturer’s instructions.

Bioinformatic processing

Amplicon sequences were quality-filtered using Trimmomatic v. 0.36 (Bolger, Lohse & Usadel, 2014) as in (Patin et al., 2019) with the following parameters: removal of bases at the start of a read below quality 3 (LEADING:3), removal of bases at the end of a read below quality 3 (TRAILING:3), read scanning with a 4-base sliding window cutting when average quality per base drops below 25 (SLIDINGWINDOW:4:25), and minimum read length of 150 bp (MINLEN:150). The intermediate size fraction (3–30 µm) dataset from site 12 contained only 92 sequences after quality filtering and was excluded from further analyses (Table S1). Quality filtered sequences were analyzed with the QIIME2 pipeline (qiime2.org) (Bolyen et al., 2019) by applying Deblur algorithm (Amir et al., 2017), which uses sequencing error profiles to remove likely artifacts instead of clustering them in operational taxonomic units. The resulting individual sequence variants (SVs) represent distinct taxa. SV taxonomy was assigned using a Bayes classifier trained on the “Silva 132 99% OTUs from the 515FB/806RB region” database (QIIME2 2018.6 version), and sequences classified as chloroplasts were removed, along with sequences that could not be assigned to a domain. A phylogeny was generated with FastTree 2 and used for calculating phylogeny-based diversity metrics, including weighted and unweighted Unifrac (Lozupone et al., 2011).

To estimate eukaryotic phytoplankton composition, we further analyzed the fraction of the QIIME2-processed sequences that were classified as chloroplast 16S rRNA. For these sequences, we first extracted the corresponding full length amplicons (post-quality filtering but pre-Deblur processing). We then compared these longer sequences to the PhytoRef database (phytoref.sb-roscoff.fr) (Decelle et al., 2015) using BLAST v. 2.2.28+ (Camacho et al., 2009) with max_target_seqs set to five. Taxonomic identity was assigned based on top BLAST match (highest % sequence similarity). For queries having multiple top matches (taxa) with identical % similarity, taxonomy was assigned based on the taxon most likely to occur in Gulf coastal waters (based on literature). In three cases, a query had two top matches, both of which were equally likely to occur at the sampling site; for these, we noted the ambiguity by retaining both taxonomic names.

Statistical analyses

Alpha diversity (Shannon diversity and Pielou’s evenness) was calculated in QIIME2 based on 10,000 reads per sample. Diversity values were compared between size fractions and bloom levels using box plots generated in R using the ggplot2 package (Wickham, 2016).

SV counts were transformed using the varianceStabilizingTransformation function in the DESeq2 library (Love, Huber & Anders, 2014) in R to obtain a matrix of counts that were approximately homoscedastic. A principal components analysis (PCA) was performed with the resulting matrix based on singular value decomposition using the prcomp function in R. Points on the PCA plot were colored according to bloom level and size fraction. A two-way permutational analysis of variance (PERMANOVA) was performed to test for significant effects of both size fraction and bloom level using the vegan package in R (Oksanen et al., 2019).

To obtain an overview of taxonomic trends among size fractions and between high and low bloom sites, box plots were generated using the Seaborn package in iPython. The small size fraction dataset from site 12 was excluded from these comparisons because it featured several major taxa with outlier trends compared to all other sites (the intermediate sample from this site was already excluded due to low sequence counts). We included the small size fraction from site 12 in all statistical analyses described elsewhere.

To statistically test for SVs differentially abundant between the two bloom conditions, raw SV counts were used to generate a beta-binomial regression model with the corncob package in R (Martin, Witten & Willis, 2019). The beta-binomial model is uniquely appropriate for amplicon sequence data as it allows for zero values (absences of an SV in samples), variability in per-sample total counts, and overdispersion in relative abundance values (a common property of microbiome data). Furthermore, the model allows for detection of both differentially-abundant and differentially-variable taxa. Here, both were identified with the differentialTest function using the Wald likelihood ratio test with no parametric bootstrap algorithm and a desired Type I error rate (fdr_cutoff) of 0.05. Separate analyses were run for the small and intermediate size fractions.

To test for significant associations between taxa, correlation networks were generated in the program Maximal Information Nonparametric Exploration (MINE) (Reshef et al., 2011). MINE identifies a range of correlation function types (e.g., linear, exponential, or periodic) based on maximal information coefficient (MIC) values calculated from raw SV counts. K. brevis counts (based on microscopy) were included in this analysis along with the SV count data. Separate analyses were run for the small and intermediate size fractions. For each, all pairwise association values above the statistically significant level (as determined by sample size-based p-values) were extracted from the results. Of these, associations with an MIC value greater than 0.85 were used to generate networks in Gephi (Bastian, Heymann & Jacomy, 2009). A cutoff of 0.85 was chosen because it limited networks to visually interpretable sizes. Undirected networks were generated using a Fruchterman Reingold layout and modularity was calculated using a randomized community detection algorithm with a resolution of 1.0. Nodes colored by modularity class and sized by average weighted degree. All nodes connected to K. brevis were labeled with font size proportional to the weighted degree of the node. For the intermediate size fraction, the nodes of five modules with no connections to the main network were also labeled.

Extract processing for metabolomics

The 0.2 and 3.0 µm disk filters were stored at 4  °C upon arrival at Georgia Tech until extraction using a protocol modified from Poulson-Ellestad et al. (2014). Filters were flash frozen with liquid nitrogen and ground with a mortar and pestle, and metabolites were extracted with methanol (LC/MS grade). Filter particulates were pelleted by centrifugation and the methanol extract was transferred to a fresh vial. The filter pellet was then rinsed 3 times with methanol (LC/MS grade). Methanol extracts were filtered through 0.2 µm polytetrafluoroethylene (PTFE) disk filters to remove any remaining suspended particles and then dried in vacuo on a ThermoSavant SpeedVac vacuum concentrator. Dried methanol extracts were stored at −20 °C until next steps. The dried extracts were then re-dissolved in a chloroform/water/methanol solution and partitioned twice more with water and methanol, resulting in a 19:27:30 solvent ratio. The chloroform/methanol and methanol/water phases were each filtered through a 0.2 µm nylon disk filter and dried on a ThermoSavant SpeedVac vacuum concentrator, resulting in a lipid-soluble extract (from the chloroform/methanol phase) and an aqueous extract (from the methanol/water phase).

¹H NMR spectroscopy and metabolomics analysis

¹H NMR spectral data were acquired from extracts representing a uniform concentration of ∼3 X 10⁵ K. brevis cells/mL based on microscopy counts (Table 1). The lipid-soluble extracts were prepared for analysis in deuterated chloroform containing 1% trimethylsilane (TMS); the aqueous extracts were too dilute for NMR spectroscopy. ¹H NMR spectra of lipid-soluble extracts were acquired on an 800 MHz Bruker Ascend™ NMR spectrometer using 400 scans per sample with a spectral width of 12 ppm centered at 5.5 ppm. Spectra were aligned by setting the internal standard peak (TMS) at 0.0 ppm, manually phased, baseline corrected, and segmentally aligned. For both the small (0.2–3.0 µm) and intermediate (3.0–30 µm) size fractions, spectral regions around TMS (−0.3 to 0.3 ppm), a negative peak (5.45–5.50 ppm), the residual chloroform signal (7.235–7.285 ppm for the small size fraction, 7.20–7.325 ppm for the intermediate fraction), and the excess area at the end with no peaks (8.0–11.5 ppm) were excluded from analysis. Spectral regions around containments (7.65 to 8.00 ppm) were also excluded from the small size fraction. Before multivariate analysis, NMR spectra were binned (0.005 ppm) and preprocessed via probabilistic quotient normalization, generalized logarithmic transformation, and mean-centering. Orthogonal partial least squares discriminant analysis (PLS-DA) was used to investigate variance in the chemical profiles between lipid-soluble extracts of samples classified as “high” (>100 K. brevis cells/mL) and “low” (<100 K. brevis cells/mL) bloom (MATLAB, version 8.1.0, PLS Toolbox, version 8.0.2; Eigenvector Research).

Sparse canonical correlation analysis of microbiome and metabolome data

To identify features driving the composition of both microbiome and metabolome data, transformed SV counts and NMR peak intensities were combined in a sparse canonical correlation analysis (sCCA) using the PMA package in R (Witten & Tibshirani, 2009). This method compares sets of features across high-dimensional data tables and identifies a subset of features that capture the most covariance (i.e., signals present in both tables) (Holmes & Huber, 2019). Samples (microbiome or metabolome) with counts of zero for six or more features were filtered out of the matrices and the sCCA function was run with the following sparsity penalties: small size fraction penaltyx = 0.2 and penaltyz = 0.05, intermediate size fraction penaltyx = 0.1 and penaltyz = 0.15. Separate analyses were run for the small and intermediate size fractions.

Phytoplankton community structure

Karenia brevis population densities ranged from 0 to 1463 cells/mL among sites sampled over two days in June 2018 (Fig. 1, Table 1, Table S1). All sites south of Tampa Bay (#1–6) had levels above 100 cells/mL, corresponding to a “medium” bloom level designated by the Florida Fish and Wildlife Conservation Commission (FWC). All sites north of the entrance to Tampa Bay (#7-12) contained fewer than 100 cells/mL, with site 12 being the only site with no detectable K. brevis cells. We therefore divided the sample sites into two groups: those associated with “high” or “low” K. brevis bloom levels, corresponding to southern versus northern sites, respectively (Fig. 1). Oceanographic models showed a relatively modest, but uniform, influence from lower salinity bay waters. All sites lay between the same salinity isolines at the time of sampling (33-34 psu), with the exception of site 12, which was situated between 31 and 32 psu isolines during the sampling time (http://ocgweb.marine.usf.edu/Models/FVCOM/month_tb_ch.php?year1=2018&mon1=jun).

Other (non-K. brevis) dinoflagellates were largely absent from samples, with the exception of the site 12 sample, which contained >11,000 cells/mL belonging to the genus Gonyaulax while lacking K. brevis altogether (Table 1). Unicellular diatom abundances ranged from 73 to 1,200 cells/mL, and colonial diatom abundances ranged from 0 to 1001 cells/mL (Table 1). The abundances of phytoplankton groups as categorized here were not correlated with K. brevis levels (correlation values not shown).

Taxonomic analysis of chloroplast 16S rRNA gene sequences confirmed the presence of several groups of diatoms in all samples, including the marine families Chaetoceraceae, Cymatosiraceae, and Triceratiaceae (Table S2). Other phytoplankton identified by chloroplast sequences included chrysophytes, also called golden algae, and prymnesiophytes, which include coccolithophores. No Karenia brevis or Gonyaulux sp. chloroplast sequences were identified.

Microbiome diversity, composition, and connectivity

Sequence read counts following trimming and quality filtering ranged from 10,961 to 66,447 per sample (excluding the site 12 sample from the 3–30 µm size fraction, which was removed due to low counts) (Table S1). Post-Deblur counts ranged from 4,348 to 33,394 (Table S1).

Microbiome community composition and alpha diversity were strongly influenced by size fraction (0.2–3 µm, 3–30 µm, and >30 µm) (Fig. 2, Figs. S2, S3; PERMANOVA p = 0.001). Measures of Faith’s phylogenetic diversity were highest in the largest fraction and lowest in the free-living (0.2–3.0 µm) fraction. Similarly, in most cases, taxonomic groups in the intermediate size fraction (3–3.0 µm) were at intermediate relative abundance compared to abundance in the smallest and largest size fractions; i.e., there were only a few taxonomic groups enriched or depleted in the intermediate fraction compared to the other two fractions. These included Euryarchaeota (Marine Group II clade, hereafter referred to as MGII) and, most notably, the class Cyanobacteria (phylum Cyanobacteria, hereafter referred to as cyanobacteria), both of which were present in higher proportions in the intermediate fraction (Figs. S2, S4). In comparison to the particle-associated (3–3.0 and >30 µm) fractions, the free-living fraction contained relatively more Proteobacteria (of all classes, excluding Deltaproteobacteria, which was most abundant in the largest size fraction), Deferribacteres, and Acidobacteria, but relatively fewer cyanobacteria and members of the phylum Bacteroidetes, particularly of the class Sphingobacteria (Figs. S2, S4). The largest (>30 µm) size fraction contained higher proportions of Verrucomicrobia and Planctomycetes, notably Planctomycetes of the uncultured class OM190 (Fig. S4), but lower proportions of Actinobacteria and Archaea compared to the other fractions. The >30 µm fraction also contained a substantial portion (2–10%) of a taxon unclassified beyond the domain Bacteria (Fig. S4).

Community composition also differed significantly based on bloom level (>100 cells/mL or “high”, vs. <100 cells/mL or “low”; PERMANOVA p = 0.002). This effect was observed within each size fraction, although separation by bloom level was less robust compared to separation by size fraction, with some overlap among the high and low bloom level microbiomes (Fig. 2, Table S3). Alpha diversity, however, did not differ significantly between bloom levels (Fig. S3, Table S4).

Several broad taxonomic groups differed notably in relative abundance between high and low bloom conditions (Fig. S5), with these differences most pronounced in the particle-associated fractions. In the free-living size fraction, taxon-specific variation based on bloom level was minor and evident only for the classes Flavobacteriia (phylum Bacteroidetes, hereafter referred to as flavobacteria) and Acidobacteria (Fig. S5). In the intermediate size fraction, flavobacteria, Deferribacteres (phylum Deferribacteres), and Deltaproteobacteria were differentially abundant between high and low bloom sites (Fig. S5). Of these taxa, only flavobacteria were enriched in high bloom sites.

To identify individual SVs contributing to these broad taxonomic patterns and statistically verify apparent differences, we applied two analyses, a binomial regression-based test (using corncob) to identify SVs significantly enriched/depleted in the high bloom dataset and a network analysis (using MINE) to identify SVs significantly correlated with K. brevis cell counts. For the small size fraction, corncob identified 20 SVs enriched and 17 SVs depleted under high bloom conditions (Fig. 3). For the intermediate size fraction, corncob identified 21 SVs enriched and 21 depleted under high bloom conditions (Fig. 3). MINE analysis further identified 31 and 34 SVs significantly associated with K. brevis in the small and intermediate size fraction networks, respectively (Fig. 4, Tables S5–S7). Network modularity was 0.537 and 0.605 for the small and intermediate size fraction datasets, respectively. Modularity is a value between -1 and 1 measuring the degree to which data features (nodes, in this case representing taxa) partition into modules (communities of connected nodes). These positive values indicate a greater number of connections between nodes within modules than expected by chance. The small and large size fraction networks contained 13 and 15 modules, respectively, with modules ranging in size from two to 56 (small fraction) or two to 69 (intermediate fraction) nodes. The largest module occurred in the intermediate size fraction and included the K. brevis node.

In the smallest size fraction, 25 SVs were both significantly associated with bloom state (either high or low; as identified by corncob) and significantly associated with K. brevis counts (MINE). For 11 of these SVs, the MIC value, a measure of two-variable dependence that is approximately equal to R² but considers both linear and nonlinear associations, was greater than 0.99 (Table S8). These SVs included four flavobacteria, three Gammaproteobacteria, one member of the Alphaproteobacteria, two non-flavobacteria Bacteroidetes, and one member of the phylum Cyanobacteria.

In the intermediate size fraction, which presumably contained the bulk of K. brevis cells before filtration, 29 SVs were identified as associated with both bloom state (either high or low) and K. brevis counts; 10 of these had MIC values >0.99 (Table S8). Of these 29, 8 were enriched in high bloom conditions (Table 2). These included four flavobacteria, as well as SVs of the Puniceicoccaceae, SAR86, and Alteromonadaceae groups; SAR86 was also detected by both corncob and MINE analysis of the small size fraction. One of the 8 SVs could only be classified to the domain Bacteria. The remaining 21 SVs, which included three planctomycetes and a member of the MGII euryarchaeota, were enriched under low bloom conditions.

Metabolomes

Multivariate statistical analysis (PLS-DA) of ¹H NMR spectra representing size-fractionated plankton community metabolomes differentiated samples from high and low bloom sites (Fig. 5). High bloom metabolomes were less variable among samples compared to low bloom metabolomes for the intermediate size fraction, as shown by the tighter intra-class clustering and significantly smaller 95% confidence ellipse for samples in the high bloom class (Fig. 5B). Based on the non-overlapping confidence interval areas, the differentiation of high bloom and low bloom metabolomes was greatest in this intermediate size fraction containing K. brevis, other phytoplankton, and their attached bacteria. For the smallest size fraction containing free-living bacteria, metabolomes of the low bloom samples clustered tightly within the high bloom confidence interval (Fig. 5A). In contrast, the high bloom samples of this smallest size fraction showed greater variability. (Note that we did not generate metabolome data from extracts of the >30 µm fraction.)

¹H NMR spectral analysis of lipid-soluble bloom extracts was consistent with the presence of aromatic compounds (with peaks between 6.8 ppm and 7.4 ppm) (Figs. 6A, 7A). These features were present in some high bloom extracts of the intermediate size fraction, but not in all high bloom samples, and in fewer low bloom samples (Fig. 7A). The enhanced concentrations of aromatic compounds in the high bloom extracts was made clear by signals at 6.8, 7.1, and 7.4 ppm in the loadings plot being positive on latent variable 2 (Fig. 7B) which best separated high and low bloom extracts in the scores plot (Fig. 5B). Similar NMR signals (6.8–7.4 ppm) in the smallest size fraction were also consistent with enrichment of aromatic compounds in high bloom samples with these indicative peaks positive on latent variable 1 (Fig. 6B), which best separated high and low bloom extracts in the scores plot (Fig. 5A), and largely only present in high but not low bloom extracts (Fig. 6A). Broad singlets, at 4.3 and 4.7 ppm, and peaks with varying multiplicities at 3.2–3.3 ppm, were indicative of alcohol, amine, ester, or halogen-containing metabolites, which occurred in the intermediate size fraction in both high and low bloom samples (Figs. 6A, 7A). These signals occurred more often and with greater prominence in the high bloom extracts, as shown by the positivity on latent variable 2 (Fig. 7B), although they were not found in all samples in either bloom level (Fig. 7A). Peaks in similar spectral regions (3.1–3.2 and 4.0–4.9 ppm) were found in the smallest size fraction (Fig. 6A); as in the intermediate size fraction, their positivity on latent variable 1 indicates these peaks drove the variation between the high and low bloom extracts primarily by their prominence in the high bloom extracts (Fig. 6B).

Additional NMR signals were consistent with the presence of lipids in both size fractions. These features, which include those within the vinylic region (5.0–6.0 ppm), alkyl region (1.2–1.6 ppm), allylic region (1.8–1.9 ppm), and the region adjacent to carbonyls (2.1–2.9 ppm) (Figs. 6A, 7B), also drove the variation between high and low bloom metabolomes (Figs. 5A & 5B). The majority of statistically significant NMR signals separating the two bloom levels were attributed to the high bloom extracts, as shown by the abundance of signals in lipid-related regions that were positive on latent variables 1 (Fig. 6B) and 2 (Fig. 7B) in the loadings plots for the small and intermediate size fractions, respectively. While these features were not present in all samples, they occurred in more of the high than low bloom extracts for both size fractions. These patterns are represented in the ¹H NMR spectra, where exemplary peaks that vary between high and low bloom metabolomes are numbered both in the raw spectra (Figs. 6A, 7A) and in the loadings plots (Figs. 6B, 7B). Signals indicative of brevetoxins were not detected in any of the spectra. Overall, statistical analyses indicate high bloom metabolomes were richer in lipids and aromatic metabolites compared to low bloom metabolomes (Figs. 6B, 7B), with substantial site-to-site variation in chemical profiles (Figs. 6B, 7B).

Sparse canonical correlation analysis

We applied sCCA to identify features driving separation of combined microbiome and metabolome data. Samples of the smallest size fraction separated largely along one axis (79% of variation) into three main groups, as follows: (1) sites 2, 4, 6, 7, 8, 9, (2) sites 1, 3, 11, and (3) sites 5, 10 (Fig. 8A). Sixteen SVs and seven NMR peaks were identified as driving this separation (Table S9). Samples of the intermediate size fraction separated into two groups along one axis (94% of variation), as follows: (1) sites 1, 2, 3, 5, 10, 11, and (2) sites 4, 6, 7, 8, 9 (Fig. 8B). Six SVs and 52 NMR peaks were identified as driving this separation (Table S9).