Substantial microbial community shifts in response to an exceptional harmful algal bloom in coastal Southern California

Abstract

Phytoplankton blooms create organic matter that stimulates entire marine ecosystems, including other components of the microbial community. How the ecosystem responds varies depending on the intensity, duration, and composition of the bloom. When the bloom has a direct or indirect negative impact on the ecosystem, it is termed a harmful algal bloom (HAB). HAB frequency is expected to increase in response to changing oceanic conditions and coastal nutrient supply. Characterizing the response of the bacterial and archaeal communities to HABs will improve our understanding of the ecological impacts of these phenomena. We utilized time series of chlorophyll a, phaeophytin, dissolved oxygen, flow cytometry cell counts, and microbial community structure (assessed via 16S rRNA gene sequences) maintained by several observing programs to investigate how the microbial community was affected by an exceptional bloom of Lingulodinium polyedra in coastal Southern California. These multi-year datasets allowed us to compare the microbial community response to past events, such as a smaller L. polyedra bloom the previous year. We demonstrated that the bacterial and archaeal response to the 2020 bloom was unique taxonomically, with many novel heterotrophs, and higher trophic state variance. The measured heterotrophic response to the bloom resulted in massive oxygen drawdown and may have impacted the length of the bloom and contributed to a secondary diatom bloom following the main HAB event. Taken together, these data illustrate how the massive 2020 L. polyedra bloom created unique ecological conditions for coastal Southern California.

Introduction

Phytoplankton perform many ecosystem functions with positive outcomes for coastal ecosystems. However, large-scale phytoplankton blooms can be detrimental, with negative impacts to public health, coastal economies, and the seafood industry (Anderson, 1995; Anderson et al., 2015). Large-scale, harmful algal blooms (HABs) result in the production of massive amounts of organic matter that can lead to oxygen draw-down and hypoxic or anoxic conditions through microbial respiration. In some cases, these blooms also produce phytotoxins that have a direct impact on upper trophic level consumers (Anderson, 1995; Anderson et al., 2015). Understanding the ecological consequences of such events is vital, as HABs are occurring at increased frequencies with greater impacts due to global change (Lewitus et al., 2012; Gobler, 2020).

The ecological response to large-scale harmful blooms is often widespread and includes heterotrophic members of the microbial community. These bacteria, archaea, and heterotrophic eukaryotes are particularly important for remineralizing bloom-derived organic matter (Buchan et al., 2014). Previous work has identified consistent taxonomic shifts in portions of the heterotrophic community in response to ‘typical’ sized blooms (Teeling et al., 2016; Wilson et al., 2021) and succession in the bacterial community has been observed following a HAB of Heterosigma akashiwo (Matcher et al., 2021). In particular, members of Flavobacteriales (Kirchman, 2002; Avci et al., 2020; Ferrer-Gonzalez et al., 2020; Matcher et al., 2021; Wilson et al., 2021) and Rhodobacterales (Teeling et al., 2012; Buchan et al., 2014) have both been shown to increase in response to phytoplankton-derived polysaccharides from variously sized blooms. A majority of work analyzing microbial responses to phytoplankton blooms, however, are either based on time series with ‘typical’ seasonal and small non-seasonal blooms (Wear et al., 2015; Teeling et al., 2016; Wilson et al., 2021) or based on opportunistic encounters with large-scale harmful blooms (Matcher et al., 2021), which makes comparing the magnitude of change in microbial taxonomy and function to HABs difficult. We were able to circumvent this shortcoming when an exceptional bloom of the dinoflagellate Lingulodinium polyedra occurred in coastal Southern California over an area where existing time-series data are collected by the Scripps Ecological Observatory (SEO), the Southern California Coastal Ocean Observing System (SCCOOS), and the Scripps Ocean Acidification Real-Time Monitoring Program (SOAR) (Table 1).

Table 1.

Acronym identification guide

Acronym	Full Name	Definition/Purpose
AF	Auto-fluorescent	Referencing unstained chlorophyll-containing flow cytometry groups (cells that auto-fluoresce)
AOU	Apparent oxygen utilization	An indication of the amount of oxygen biologically utilized (the difference between DO at saturation and measured DO), where negative values indicate autotrophy and positive values indicate heterotrophy
ASVs	Amplicon sequence variants	Inferred unique sequences obtained via high throughput sequencing of DNA marker genes (16S rRNA genes were used in this study)
DO	Dissolved oxygen	The concentration of oxygen present in water
HAB	Harmful algal bloom	A rapid and large-scale algal bloom that causes ecological harm/damage
MM	Module membership	Within WGCNA, how well each taxon’s abundance pattern fits with the subnetwork it was assigned (calculated by way of a Pearson correlation coefficient)
SCCOOS	Southern California Coastal Ocean Observing System	A research group that collects scientific data to understand coastal systems and HABs in Southern California, with a node at the Ellen Browning Scripps Pier
SEO	Scripps Ecological Observatory	An ecological observatory based at Scripps Institution of Oceanography that collects and leverages scientific data related to microorganisms (identity and function) from the Ellen Browning Scripps Pier
SG	SYBR Green I	Referencing non-auto-fluorescent flow cytometry groups stained with SYBR Green I (which binds to double-stranded DNA and fluoresces)
SOAR	Scripps Ocean Acidification Real-Time Monitoring Program	A monitoring program that focuses on measurements of real-time ocean acidification and DO at the Ellen Browning Scripps Pier
SOM	Self-organizing map	A machine learning technique which can be used to reduce multidimensional data to a single categorical variable for each sample event/date
TM	Taxonomic mode	The term used to represent the different taxonomic (16S rRNA gene-based) SOM community types for a given date; similar SOM map units were grouped into distinct TMs using K-means clustering
TOM	Topological overlap measure	A calculation completed during WGCNA that demonstrates network interconnectedness by assessing the connection strengths of shared 16S rRNA genes for each pair of taxa in the network
WGCNA	Weighted gene correlation network analysis	A network analysis tool used to identify groups of co-occurring microorganisms (based on 16S rRNA gene data) throughout samples and relate those groups (or subnetworks) to ecological and environmental data

Lingulodinium polyedra is a bioluminescent, red-tide-forming dinoflagellate species commonly observed along the western coast of the United States and is capable of forming HABs (Kahru and Mitchell, 1998; Kudela and Cochlan, 2000) with low oxygen conditions and the production of yessotoxin (Howard et al., 2008; Lewitus et al., 2012). The range of L. polyedra includes nearshore regions of the Southern California Bight, though blooms in this area prior to 2020 have not been classified as harmful (Allen, 1943; Holmes et al., 1967; Kahru and Mitchell, 1998; Moorthi et al., 2006). However, a massive bloom of L. polyedra resulted in the highest concentrations recorded by the Scripps Pier Station node of SCCOOS in April and May of 2020. In addition to visibly red and bioluminescent water that smelled strongly of sulfur compounds, 25% of beachgoers reported minor health symptoms often observed during HABs (including sore throat, watering eyes, and dermatitis), and there was a massive die off event for invertebrates and fishes (Anderson and Hepner-Medina, 2020).

The purpose of this study was to assess the bacterial and archaeal community response to the unprecedented 2020 Lingulodinium polyedra bloom using the existing long-running Scripps Pier time series. This included characterizing chlorophyll a, phaeophytin, dissolved oxygen (DO), apparent oxygen utilization (AOU), cell counts using flow cytometry, and microbial community structure (assessed via 16S rRNA gene sequences) before, during, and after the bloom. To place our observations in the context of past seasonal and stochastic changes, we used similar machine learning and network analysis methods as a previous time-series analysis at the site (Wilson et al., 2021). That work, which relied on the first 18 months of many of the data products presented in this study, found that a minority subnetwork correlated with all changes in phytoplankton and phytoplankton-derived products against a backdrop of seasonal changes in the microbial community. This approach allowed us to determine the novelty of the bacterial and archaeal response as we were able to compare it to a lesser L. polyedra bloom in 2019 and various other seasonal phytoplankton blooms. Our findings suggest that the 2020 bloom included a novel taxonomic bacterial and archaeal community, with implications for concurrent ecosystem function and trophic status. Following a period of succession, the microbial community returned to seasonal norms after bloom collapse.

Materials and methods

Sample collection

A list of all acronyms used in this paper can be found in Table 1. All data were collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). Scripps Ecological Observatory (SEO; Bowman et al., 2021) conducts bi-weekly (i.e., twice a week) analysis of the microbial community in conjunction with the Southern California Coastal Ocean Observing System (SCCOOS, 2022) Scripps Pier Shore Station sampling effort. As described in Wilson et al. (2021), water samples collected by SCCOOS were analyzed for bacterial and archaeal community structure, determined through sequencing of the 16S rRNA gene, and microbial abundance via flow cytometry.

Surface water was collected using an acid-washed bucket off the Scripps Pier at approximately 10:00 AM and transferred back to the laboratory in acid-washed containers for immediate analysis on a bi-weekly basis. As described in Hatch et al. (2013), bi-weekly concentrations of chlorophyll a and phaeophytin and weekly counts of Lingulodinium polyedra, dinoflagellates, diatoms, and total phytoplankton were provided by the McGowan Plankton and Chlorophyll Program and SCCOOS (2022). Continuous temperature and water depth (hydrostatic pressure) at the Scripps Pier were provided by the SCCOOS Automated Shore Station Program (SCCOOS, 2022).

A series of MiniDOT optical dissolved oxygen (DO) and temperature sensors (Precision Measurement Engineering Inc., Vista, CA, USA), maintained by the Scripps Ocean Acidification Real-Time Monitoring Program (SOAR, 2022), were used to measure DO and temperature at a depth of 5 m. The time intervals of measurements ranged from 5 to 15 minutes from August 9, 2018, through October 31, 2020. A total of four different dataloggers were utilized to minimize gaps in data collection as sensors were cleaned and data downloaded throughout the time series (with the maximum deployment lasting 84 days from February 14, 2020, through May 8, 2020, and the minimum deployment lasting 5 days from August 9, 2018, through August 14, 2018). Sensors were checked once or twice a month for biofouling, with minor biofouling observed in early September 2018 and, to a lesser degree, November 2018 (Figure S1). Sensors were swapped at the first signs of growth during regular maintenance checks and the introduction of copper mesh coverings helped to prevent further biofouling starting October 10, 2018. SOAR applies a constant salinity of 33.5 to correct DO and DO equilibrium saturation (DO_sat), which was calculated with the MiniDOT datalogger software using sea level, temperature, and salinity. Apparent oxygen utilization (AOU; DO_sat−DO) was calculated for each time point, after which a daily average AOU was calculated. The daily average was selected to integrate across the day-night cycle and provide a net AOU value for a 24-hour period. Raw (unaveraged) DO and AOU time-series data along with sensor deployments are depicted in Figure S1.

Flow cytometry

Flow cytometry data were analyzed bi-weekly by SEO from October 28, 2019, through December 28, 2020. Flow cytometry samples were prefiltered (60 µm), aliquoted into 2-mL sample tubes, and fixed with 10 µl of 0.25% glutaraldehyde before running on a Guava easyCyte 11HT (Luminex, Austin, TX, USA). The Guava 11HT interrogates cells with both a blue (488 nm) and violet (405 nm) laser. The fluorescence signals of unstained chlorophyll-containing (autofluorescent; AF) samples were measured for forward scatter, side scatter, red emission (405 nm excitation/695 nm emission), and yellow emission (488 nm excitation/583 nm emission). Duplicate samples were stained with SYBR Green I (SG-stained; Molecular Probes, Inc., Eugene, OR, USA) at the manufacturers recommended working concentration to visualize all particles that possessed double-stranded DNA. These SG-stained samples were measured for forward scatter, side scatter, and green emission (488 nm excitation/525 nm emission). Absolute cell counts were determined by spiking a standard volume of 1:10 diluted 123-count beads (ThermoFisher, Waltham, MA, USA) to each sample and blank.

Subgroups were identified using a self-organizing map (SOM) from forward scatter, side scatter, red emission, and yellow emission for AF communities, and from forward scatter, side scatter, and green emission for SG-stained communities following Wilson et al. (2021). This analysis was conducted using the ‘kohonen’ package in R (Wehrens and Kruisselbrink, 2018). In brief, a training set was constructed with data from randomized sample days (December 28, 2020, and July 20, 2020, were selected). These data were trained using a toroidal map with a grid size of 41 x 41. Populations were identified using k-mean clustering, and k was chosen through the visual evaluation of a within-cluster sum of squares scree plot and a priori knowledge of populations (AF k = 7, SG-stained k = 5). The two AF and SG-stained k-mean cluster models were then used to predict the clustering of events in each AF and SG-stained flow cytometry sample, respectively.

DNA collection

Bi-weekly characterization of the microbial community via 16S rRNA gene sequencing took place from January 4, 2018, through July 2, 2020, and was a continuation of the time series presented in Wilson et al. (2021), which lasted from January 4, 2018, until June 20, 2019. Seawater was filtered through a sterile 0.2 µm Supor membrane disc filter (Pall Corporation, Port Washington, NY, USA) and stored at −80ºC until extraction. Through January 8, 2019, filters were extracted manually using the MoBio DNEASY PowerWater Kit (Qiagen, Venlo, Netherlands), after which filters were extracted using the KingFisher^TM Flex Purification System and MagMax Microbiome Ultra Nucleic Acid Extraction kit (ThermoFisher Scientific, Waltham, Massachusetts, USA). Extracted DNA was sent to Argonne National Laboratory for amplicon library preparation and sequencing using the Illumina MiSeq platform, universal primers 515F and 806 R (Walters et al., 2016), and a 2 x 151 bp library architecture. Reads were filtered, denoised, and merged with dada2 (Callahan et al., 2016) using the Bowman Lab protocols (2022), and analyzed with paprica v0.7.0 (Bowman and Ducklow, 2015). Paprica utilizes phylogenetic placement with EPA-ng (Barbera et al., 2019), Infernal (Nawrocki and Eddy, 2013), and Gappa (Czech et al., 2020) to place query reads on a reference tree constructed from the full-length 16S rRNA genes from all completed genomes in GenBank (Haft et al., 2018). All unique reads were assigned to internal branches or terminal branches on the reference tree and placements were used to estimate 16S rRNA gene copy number, genome size, and GC content. Sequences were submitted to NCBI SRA at BioProject PRJNA662174. Diversity was calculated for each sample using the inverse Simpson Index.

Determination of taxonomic modes and subnetworks

The entire 30-month 16S rRNA gene time series was assessed via microbial community segmentation (Bowman et al., 2017) using a pre-existing SOM that was created utilizing the original 18-month dataset (from January 4, 2018, through June 20, 2019) (Wilson et al., 2021). SOMs are an effective way to reduce multidimensional taxonomic datasets to a single categorical variable for each sample (‘kohonen’ package in R version 3.0.8; Wehrens and Kruisselbrink, 2018). Existing and new data were (re)assigned to the existing map units (which were arranged in a 2-D toroidal lattice) by identifying the shortest Euclidean distance between the Hellinger-transformed relative abundances of each sample and the codebook vector associated with each map unit. The codebook vector is the “representative community” that was assigned to each map unit during the training phase. Put another way, we identified which representative Hellinger-transformed relative abundance community (map unit) was most similar to each sample and grouped samples accordingly. Because all samples were reassigned to the existing SOM, the Euclidean distances indicate how similar each sample is to the codebook vector for the map unit to which it was assigned. The final clustering of map units into taxonomic modes via k-means clustering was as described in Wilson et al. (2021), and taxonomic mode assignments were propagated to samples via sample association with map units.

Only samples that had >5,000 reads were used for the SOM-based analysis. Analysis was carried out in two ways: 1) utilizing all amplicon sequence variants (ASVs); and 2) eliminating putative chloroplasts (ASVs identified as Cyanobacteria that were not from the genera Synechococcus or Prochlorococcus). Potential chloroplast sequences were included in the SOM training data of Wilson et al. (2021). Masking those ASVs allowed us to examine the microbial community response associated exclusively with heterotrophs and Cyanobacteria.

The entire 30-month time series was also re-evaluated using weighted gene correlation network analysis (WGCNA) as implemented in the R package ‘WGCNA’ (Langfelder and Horvath, 2007; 2008). This network analysis tool was used in Wilson et al. (2021) to characterize variation in subsets (or subnetworks) of the microbial community over the original 18-month dataset. WGCNA shows how co-occurring groups within a sample vary over time and in relation to environmental data. Reanalysis allowed us to assess if the same subnetworks continued to exist during the 2020 bloom. Following Wilson et al. (2021), we again used a Hellinger-transformed relative abundance matrix (that factored in all 6,757 unique taxa across all 240 samples) but limited the matrix to the 1,451 most abundant taxa so that the topological overlap measure (TOM) would fit. We utilized a signed adjacency matrix and selected a soft thresholding power of 7 which gave an R² of 0.97 for the TOM model fit and a mean connectivity between taxa of 18.4. We set our minimum module (subnetwork) size to 50, and the co-occurrence abundance profiles for each subnetwork were then related to other variables using the first principal component of each subnetwork by way of the Pearson correlation coefficient. Similarly, each taxon was related to its subnetwork using the Pearson correlation coefficient (a term called module membership or MM). Wilson et al. (2021) utilized a Spearman rank correlation for these two measures, but Pearson showed stronger correlations with certain environmental variables during the bloom and yields a more interpretable correlation. The correlation method (Pearson) for grouping taxa into subnetworks was identical between this work and Wilson et al. (2021).

We also used WGCNA to identify subnetworks from the 30-month time series excluding the 2020 Lingulodinium polyedra bloom (excluding dates after April 1, 2020), making it a 27-month time series. The same Hellinger-transformed relative abundance matrix of the 1,451 most abundant taxa went into creating the subnetworks (sans data after April 2020), though 15 ASVs were dropped initially because they were only present during and after the 2020 bioluminescent bloom and 11 additional ASVs (that were originally in the ‘bloom’ and ‘post-bloom’ subnetworks) could not be sorted due to a lack of association with subnetworks above the minimum module size threshold of 50. All model parameters were identical except for a soft thresholding power of 10, which gave an R² of 0.94 for the TOM model fit and a mean connectivity between taxa of 3.99. The higher soft thresholding power was necessary to get the R² of the TOM model fit above 0.8.

Results

Phytoplankton variation

Phytoplankton standing stocks exhibited temporal variation throughout the time series, with chlorophyll a and phaeophytin concentrations typically peaking in the spring and summer (Figure 1A and B). Chlorophyll a possessed particularly high peaks in June 2019 and April/May 2020, both periods that possessed visibly bioluminescent red tides caused by Lingulodinium polyedra (Figure 1C). During these periods, chlorophyll a and L. polyedra both peaked on June 3, 2019, and on April 30, 2020, respectively. Phaeophytin had a peak in 2020 (also on April 30, 2020) but did not have a peak to match the chlorophyll a or L. polyedra peaks in 2019 (Figure 1B).

Figure 1.

Three-year time series of pigments and phytoplankton including the red tide in 2020. Concentrations/abundances of (A) chlorophyll a, (B) phaeophytin, (C) Lingulodinium polyedra, (D) dinoflagellates, and (E) diatoms before, during, and after the 2020 bioluminescent red-tide collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). Points represent sample days, with the bloom period denoted in red (from the beginning of April through May 18, 2020—based on when the abundance of L. polyedra was elevated and water was visibly red), and the post-bloom period denoted in orange (from May 19 through June 4, 2020—based on when water no longer had red pigmentation but remained difficult to filter and possessed a unique microbial community). Vertical lines indicate April 1 and June 3; insets expand the scales for this period.

[Figure omitted. See PDF]

While Lingulodinium polyedra blooms occurred in both 2019 and 2020, the 2020 bloom possessed cell concentrations almost 50x higher than the 2019 bloom (Figure 1C). According to light microscopy counts performed during regular SCCOOS sampling (Hatch et al., 2013; Wilson et al., 2020), L. polyedra concentrations were on the order of 10² and 10³ cells L^–1 in the spring of 2020, before they increased to 10⁴ cells L^–1 on March 30, 2020. Seawater had visible red bands by April 7, 2020, and filters had visible red pigmentation beginning April 16 through May 14, 2020. Water became universally red (no more banding) by April 20, 2020, when L. polyedra cells reached concentrations greater than 10⁶ cells L^–1. L. polyedra peaked on the order of 10⁷ cells L^–1 (making up almost all enumerated cells; Figure 1D) and chlorophyll a reached concentrations of more than 1,300 mg m^-3. The concentration of L. polyedra remained high until May 11, though the water remained red and other slime-producing dinoflagellate species were present afterward (Figure 1C and E). Diatoms possessed a very strong peak following the 2020 HAB, reaching a maximum on May 26, 2020, as the concentration of L. polyedra was still decreasing (Figure 1E). Diatoms did not possess a corresponding peak following the 2019 L. polyedra bloom.

Temperature, DO, and AOU variation

Temperature, DO, and AOU varied throughout the time series prior to the 2020 Lingulodinium polyedra bloom (Figure 2). From August 9, 2018 (the start of the temperature/DO time series), through March 31, 2020, daily average temperature followed a seasonal cycle with maxima in the summer and minima in the winter. Daily average DO concentrations had a mean value of 241 µmol L^–1, a median value of 242 µmol L^–1, and standard deviation of 14.6 µmol L^–1. During this period, DO reached a maximum of 332 µmol L^–1 on June 2, 2019, and a minimum of 179 µmol L^–1 on September 17, 2018. Similarly, daily AOU over this period had a mean value of 1.82 µmol L^–1, a median value of 3.79 µmol L^–1, and a standard deviation of 16.0 µmol L^–1. AOU reached a minimum of −89.1 µmol L^–1 on June 2, 2019, and a maximum of 50.2 µmol L^–1 on September 17, 2018. These results show a seasonal pattern of high O₂ production/low O₂ consumption accompanied by a peak in DO in the late spring/early summer, followed by high O₂ consumption and a dip in DO at the end of the summer.

Figure 2.

Time series of water temperature, dissolved oxygen, and apparent oxygen utilization. Time series of (A) daily average water temperature, (B) daily average dissolved oxygen concentration, and (C) daily average apparent oxygen utilization from August 9, 2018, through October 31, 2020. All data were collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). Points represent sample days, with the bloom period denoted in red (from the beginning of April through May 18, 2020), and the post-bloom period denoted in orange (from May 19, 2020, through June 4, 2020). Vertical lines indicate April 1 and June 30.

[Figure omitted. See PDF]

The period from April 1, 2020, through June 30, 2020 (which encompasses the Lingulodinium polyedra bloom and post-bloom period), possessed DO concentrations and AOU values that were well outside and more variable than the preceding 20 months (Figure 2). There was an increase in DO and a decrease in AOU at the beginning of the bloom, which reached a maximum DO concentration of 331 µmol L^–1 and a minimum AOU value of −83.9 µmol L^–1 on April 17, 2020. The 2020 peak mimicked early June the previous year but DO decreased and AOU increased in late April midway through the 2020 bloom. DO concentrations reached an all-time low of 22.1 mmol L^–1 while AOU reached an all-time high of 216 µmol L^–1 on May 5, 2020. This period also had increased variance (standard deviation for DO and AOU were 56.2 and 54.9 mmol L^–1, respectively). DO and AOU subsequently returned to pre-bloom levels in June and then decreased (for DO) and increased (for AOU) sharply again in July and in August.

Flow cytometry variation

Six of the seven AF populations were elevated at some point during the 2020 Lingulodinium polyedra bloom, with five populations spanning across all size classes reaching all-time maxima during the bloom (Figure 3). The largest AF size group with the highest chlorophyll a fluorescence (AF group 1 and likely diatoms) reached a maximum on May 21, 2020, shortly after the sudden drop in L. polyedra abundance. This group had an earlier local peak on April 30, 2020, in addition to others later in the summer. All other AF groups except medium-sized AF cells with high chlorophyll a fluorescence (AF group 3 which were likely small eukaryotic phytoplankton) and small-sized AF cells with low chlorophyll a fluorescence (AF group 5 which was a population of Cyanobacteria) reached maxima on either April 30 or May 4, 2020. However, it should be noted that AF group 3 was still elevated during the bloom, with a local maximum on April 30, 2020. AF group 5 (a different population of Cyanobacteria) was the only AF group to reach an all-time minimum during the bloom (on May 11, 2020) and was generally lowest during the 2020 L. polyedra bloom.

Figure 3.

Time series and chlorophyll a florescence of autofluorescent groups. Time series of (A–G) different autofluorescent (AF) populations (AF groups 1–7) from October 28, 2019, through December 28, 2020; and (H) a forward scatter vs chlorophyll a fluorescence plot of the different AF populations for combined data from July 20 and December 28, 2020. All data were collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). For the time series, points represent sample days, with the bloom period denoted in red (from the beginning of April through May 18, 2020), and the post-bloom period denoted in orange (from May 19, 2020, through June 4, 2020). Vertical lines indicate April 1 and June 30.

[Figure omitted. See PDF]

SG-stained populations showed similar trends to various AF populations, with three of the five SG-stained populations reaching maxima during the 2020 Lingulodinium polyedra bloom and one of the remaining two reaching an all-time minimum during the bloom (Figure 4). Specifically, the medium- and larger-sized SG-stained cells (SG-stained groups 1–3) reached maxima on either May 4 or May 18, 2020. Meanwhile, small SG-stained cells with medium/high nucleic acid content (SG-stained group 4) peaked on June 25, 2020, and small SG-stained cells with low nucleic acid content (SG-stained group 5) reached an all-time minimum on May 4, 2020. Finally, both medium- and small-sized SG-stained cells with medium/high nucleic acid content (SG-stained groups 3 and 4) reached minima on March 30, 2020, as L. polyedra abundance was increasing.

Figure 4.

Time series and SYBR Green I fluorescence of SG-stained groups. Time series of (A–E) different SYBR Green I (SG)-stained populations from October 28, 2019, through December 28, 2020; and (F) a forward scatter vs SG fluorescence plot of the different SG-stained populations for combined data from July 20 and December 28, 2020. All data were collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). For the time series, points represent sample days, with the bloom period denoted in red (from the beginning of April through May 18, 2020), and the post-bloom period denoted in orange (from May 19 through June 4, 2020). Vertical lines indicate April 1 and June 30.

[Figure omitted. See PDF]

Microbial community variation

We next assessed the microbial community via analysis of 16S rRNA gene sequences. Sequencing produced 3,721,874 copy number corrected reads (with a mean library size of 15,508 ± 11,838) across 240 sampling days. The mean and standard deviation for 16S rRNA gene copy number were each elevated during the 2020 Lingulodinium polyedra bloom (with a maximum on May 4, 2020; Figure 5A). Meanwhile, genome size, GC content, and diversity were each more typical in terms of the magnitude of peaks/troughs for the spring/summer, with the timing varying slightly year-to-year (Figure 5B–D). Note that the spring/summer peak for genome size and trough for diversity aligned with the 2020 L. polyedra bloom. The highest genome size in 2020 was observed on May 4 (during the height of the bloom) and the lowest diversity in 2020 was observed on May 21 (following the sudden drop in L. polyedra on May 18).

Figure 5.

Time series of 16S rRNA gene copy number, genome size, GC content, and diversity. Time series of (A) mean 16S rRNA gene copy number, (B) mean genome size, (C) mean GC content, and (D) inverse Simpson index from January 4, 2018, through July 2, 2020, collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). Points represent sample days, with the bloom period denoted in red (from the beginning of April through May 18, 2020), and the post-bloom period denoted in orange (from May 19, 2020 through June 4, 2020). Vertical lines indicate April 1 and June 30.

[Figure omitted. See PDF]

Microbial community segmentation with an existing SOM (Wilson et al., 2021) demonstrated increased novelty during and after the 2020 Lingulodinium polyedra bloom. Because the original SOM includes chloroplast ASVs for a holistic view of microbial community structure, we evaluated the distances between all samples and codebook vectors both with and without chloroplast ASVs. Including chloroplasts, there were a similar number of map units observed during and immediately after the 2020 bloom when compared with similar months in previous years (Figure 6A). Because similar samples are assigned to the same map units, a greater number of map units for a given time period means that samples are more varied in composition. April 1 through June 30 each possessed 9 map units in 2018 and 2019 vs 7 in 2020. However, the 2020 bloom samples occupied several map units that were not occupied in previous years from April through June (Figure 6A). Additionally, Euclidean distances between each sample and the codebook vector for the closest map unit were higher during the 2020 bloom (Figure 6B). Distances generally increased during the bloom, reaching higher values than previously observed (Figure 6B). The highest value occurred on May 21, 2020, on the first sampling date of the post-bloom period, before they gradually decreased post-bloom (Figure 6B). Novelty was also observed at the level of taxonomic mode (TM). The end of the bloom and the beginning of the post-bloom period occupied a TM only observed later in the summer in previous years (Figure 6C).

Figure 6.

Time series of each sample's map unit, Euclidian distance to its assigned unit, and taxonomic mode. Time series showing sample assignment of the microbial community (assessed via the 16S rRNA gene) to the original self-organizing map (SOM) that was created using data from January 4, 2018, through June 20, 2019 (Wilson et al., 2021). All samples (from January 4, 2018, through July 2, 2020) were collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). The time series shows (A) the map unit for each sample, (B) the Euclidean distance for each sample to its assigned map unit, and (C) the taxonomic mode for each sample. Each map unit possesses a “representative community” that was derived during the machine learning phase of analysis (completed using the original 18-month dataset in Wilson et al., 2021). Because all of the samples were reassigned to the existing SOM, the Euclidean distances show how well each sample fits with the map unit community to which it was assigned. We were able to further cluster similar map units into like groups that we termed taxonomic modes using k-means clustering. Points represent sample days, with the bloom period denoted in red (from the beginning of April through May 18, 2020), and the post-bloom period denoted in orange (from May 19 through June 4, 2020). Vertical lines indicate April 1 and June 30.

[Figure omitted. See PDF]

Rerunning the SOM analysis excluding chloroplast ASVs allowed us to focus on changes in the heterotrophic community, with similar overall results. From April 1 through June 30, 2018, the heterotrophic and cyanobacterial community had 9 map units versus 8 in 2019 and 2020 (Figure S2A). The 2020 bloom samples still occupied several map units that were not occupied in previous years (Figure S2A), and the Euclidean distances between each sample and the codebook vector for the closest map unit were again highest during the 2020 bloom (Figure S2B). However, while the distances were still elevated, removing chloroplasts reduced the distances during the post-bloom phase (Figure S2B). Similar novelty with regard to TMs was observed by excluding chloroplasts (Figure S2).

WGCNA analysis of microbial subnetworks reinforced the typical seasonal patterns observed in Wilson et al. (2021) and also demonstrated the novelty of the microbial community during and after the 2020 Lingulodinium polyedra bloom. Seasonally dominant subnetworks corresponded to ‘fall/winter’, ‘winter/spring’, ‘early-summer’, and ‘late-summer’ (similar to Wilson et al., 2021; Figures 7A and 8). There were clear differences between the 18-month and 30-month time series for the less-dominant subnetworks. Specifically, no ‘spring’ subnetwork (originally identified for the 18-month time series) was observed for the 30-month time series, and a low abundance winter subnetwork appeared during the 2018/2019 winter that was not present in the 18-month time series. Additionally, a ‘chlorophyll a’ subnetwork observed in the 18-month time series—which remained non-dominant but strongly correlated with chlorophyll a, dinoflagellate abundance, and multiple autofluorescent size classes—was not present in the 30-month time series. Finally, there were two aseasonal subnetworks that were prevalent during and after the 2020 L. polyedra bloom (the ‘bloom’ and ‘post-bloom’ subnetworks Figure 7A and B). The ‘bloom’ and ‘post-bloom’ subnetworks were each made up of more than 40% new ASVs that were not present abundantly enough to be included in the initial 18-month analysis (Figure 8; Table S1).

Figure 7.

Time series of microbial subnetworks. Time series showing the sum of the Hellinger-transformed relative abundances of 16S rRNA gene sequences sorted into different weighted gene correlation network analysis subnetworks (with a module membership p-value <0.05) from January 4, 2018, through July 2, 2020. All data were collected from the Ellen Browning Scripps Pier in San Diego, CA (32.8663° N, 117.2546° W). Both (A) the entire time series and (B) the first six months of 2020 are shown. Vertical lines indicate April 1 and June 30.

[Figure omitted. See PDF]

Figure 8.

Heatmap time series of most abundant ASVs in each of the microbial subnetworks. Heatmap time series showing the Hellinger-transformed relative abundances of the top 7 amplicon sequence variants (ASVs) for each of the 8 subnetworks that were identified with weighted gene correlation network analysis, with the lowest possible taxonomic identification listed for each ASV. The 8 subnetworks are color-coded on the left axis (top to bottom): early-summer (light green), late-summer (dark green), fall/winter (blue), fall/winter 2019 (brown), winter/spring (light blue), other (black), bloom (red), and post-bloom (orange). The blue scale bar indicates the Hellinger-transformed relative abundance value. The 16S rRNA sequences and the other subnetwork members can be found in Table S1. Module membership values for displayed ASVs ranged from 0.66 to 0.96, with higher values indicating ASVs had similar abundance profiles to their subnetwork.

[Figure omitted. See PDF]

The ‘bloom’ and ‘post-bloom’ subnetworks correlated with specific environmental drivers (Figure 9) and were taxonomically distinct from one another (Figure 8; Table S1). The ‘bloom’ subnetwork consisted of 66 bacterial ASVs, with top members belonging to Rhodobacterales, Flavobacteriales, and other orders often associated with phytoplankton blooms. Meanwhile, the ‘post-bloom’ subnetwork consisted of 114 ASVs, with a high proportion of Cyanobacteria/chloroplasts (47%, with only 3% identified as Synechococcus or Prochlorococcus; Table S1). The ‘bloom’ subnetwork correlated with chlorophyll a, phaeophytin, Lingulodinium polyedra, and dinoflagellates (Pearson’s r = 0.32−0.44), while the ‘post-bloom’ subnetwork correlated with diatoms (r = 0.6). The ‘bloom’ subnetwork also correlated with low DO, high AOU, and high average 16S gene copy number and genome size (r = −0.58, 0.57, 0.79, and 0.41, respectively; Figure 9). Of note, the ‘bloom’ subnetwork was elevated during the less intense 2019 red tide, though it remained non-dominant (Figure 7A).

Figure 9.

Pearson correlation coefficients for each microbial subnetwork and the measured environmental variables. Pearson correlation coefficients for each microbial subnetwork, determined by weighted gene correlation network analysis (shown on the vertical axis), and measured environmental variables (shown on the horizontal axis). The 8 subnetworks are color-coded on the left axis (top to bottom): early-summer (light green), late-summer (dark green), fall/winter (blue), fall/winter 2019 (brown), winter/spring (light blue), other (black), bloom (red), and post-bloom (orange). The first principal component of each subnetwork’s abundance profile (based on Hellinger-transformed relative abundances) was related to environmental and ecological variables when the time series overlapped. Positive Pearson relationships are shown in red and negative in blue, with a color bar ranging from –1 to 1 on the right for reference. In each cell, the top number is Pearson’s r and the number in the parentheses is the p-value.

[Figure omitted. See PDF]

The novelty of the ‘bloom’ and ‘post-bloom’ subnetworks was further demonstrated by reassessing the longer time series with WGCNA excluding bloom dates (eliminating data from April 1, 2020, onward), making it a 27-month time series. We observed that 30% of the ‘bloom’ subnetwork’s ASVs and 53% of the ‘post-bloom’ subnetwork’s ASVs were either absent or had an MM p-value below 0.05 in the 27-month time series. Additionally, the resulting 27-month time series again possessed a ‘chlorophyll a’ subnetwork, with over 67% of ASVs from the original 18-month time series in this group. The 27-month ‘chlorophyll a’ subnetwork tended to be most abundant from April to June each year and positively correlated with chlorophyll a, Lingulodinium polyedra, dinoflagellates, DO, 16S rRNA gene copy number, genome size, and GC content (r = 0.25–0.69), and negatively correlated with AOU (r = −0.36; Figure S3). When we extended the 27-month time-series subnetworks out through the end of sampling, the 27-month ‘chlorophyll a’ subnetwork also peaked during and after the 2020 bloom but did not reach the same magnitude as the ‘bloom’ subnetwork (Figure S4).

Discussion

During the spring of 2020, several populations of primary producers bloomed in the coastal Southern California Bight, with Lingulodinium polyedra reaching concentrations orders of magnitude higher than previously observed. The coastal water quality conditions brought on by the 2020 L. polyedra bloom were unprecedented, as the spike in L. polyedra cells, chlorophyll a, phaeophytin, and strong sulfur smell all indicated that an unprecedented amount of organic matter was introduced and then broken down. We demonstrated some of the ecological effects of this bloom by combining multiple long-standing biological and oceanographic time-series datasets. Specifically, we characterized bacterial and archaeal taxonomic identity, DO, and oxygen utilization. Further, we were able to compare these data to an earlier study based on the first 18 months of the 30 months presented here. These observations allowed us to determine the novelty of the 16S rRNA gene community and other conditions and to place them in the context of typical seasonal patterns. We were also able to compare the ecological response to the 2020 L. polyedra bloom to that of a smaller and more typical red tide event in the summer of 2019. Doing so demonstrated that the 2020 bloom and post-bloom communities were taxonomically unique as many taxa were only observed at these times. The 2020 bloom and post-bloom communities were likely functionally unique as well.

The large-scale 2020 Lingulodinium polyedra bloom disrupted typical patterns of local seasonal ecology and resulted in taxonomically unique microbes, most likely through the introduction of organic matter and potential toxins. The large Euclidean distances for SOM map units starting in mid-April 2020 were a function of this disruption, as they were caused by new ASVs that were not present in the original 18 months of training data. WGCNA further revealed the degree to which natural microbial community ecology was interrupted. The seasonal ‘winter/spring’ subnetwork was increasing in early April of 2020, but was effectively interrupted from being the dominant subnetwork by the rise of the novel ‘bloom’ subnetwork. Throughout April, the ‘winter/spring’ subnetwork decreased in abundance and the non-seasonal ‘bloom’ subnetwork increased, becoming the dominant subnetwork on April 30, 2020. The ‘bloom’ subnetwork remained dominant until May 21, 2020, when the ‘post-bloom’ subnetwork became dominant. The ‘post-bloom’ subnetwork then dominated until June 8, 2020, at which point ‘typical’ seasonal succession resumed and the ‘early-summer’ subnetwork increased. However, the ‘post-bloom’ subnetwork remained above background levels through the end of sampling. Despite the fact that the ‘bloom’ subnetwork was putatively responding to phytoplankton-derived organic matter, many top ASVs were not in the original 18-month time series (Wilson et al., 2021) or the longer time series excluding 2020 bloom dates. The fact that the ‘chlorophyll a’ subnetwork (which responded to more typical-sized blooms) was not present during the time series that included the 2020 bloom is because only a portion of its members responded to the L. polyedra bloom, likely due to the unique combination of extremely high organic matter, dinoflagellate-produced toxins (Paz et al., 2004), and extremely low DO. Additionally, when we extended the ‘chlorophyll a’ subnetwork through the 2020 bloom, it did not have a strong peak because novel ASVs were so important to the ‘bloom’ subnetwork.

Locally unique and specialized microbes likely took advantage of the niches created by the massive amount of organic matter produced during the 2020 Lingulodinium polyedra bloom. The ‘bloom’ subnetwork largely consisted of non-photosynthetic r-type growth strategists, as indicated by the increase in medium to large non-autofluorescent microbial cells (SG-stained groups 1–3), the increase in 16S rRNA gene copy number, and the taxonomic make-up of the ‘bloom’ subnetwork. Slow-growing bacteria from low-nutrient aquatic environments have been found to possess a low number of 16S rRNA genes (Button et al., 1998; Fegatella et al., 1998) and multiple studies have shown that bacteria that respond rapidly to substrates tend to have more 16S rRNA genes (Kappenbach et al., 2000; Chen et al., 2020). Additionally, many top members of the ‘bloom’ subnetwork were identified as Flavobacteriaceae and Rhodobacteraceae. The genomic properties of Flavobacteriales and Rhodobacterales, such as large genomes and higher gene content, enable a fast response to transient nutrient pulses (Buchan et al., 2014). As such, both Flavobacteriales and Rhodobacterales have often been associated with blooms, with Flavobacteriales shown to specialize in the degradation of complex organic matter (Pinhassi et al., 2004; Edwards et al., 2010; Thomas et al., 2011; Gómez-Pereira et al., 2012; Teeling et al., 2012) and Rhodobacterales thought to consume low molecular weight organic matter (Landa et al., 2017; Ferrer-Gonzalez et al., 2020).

The life histories of ASVs and assigned taxa within the ‘bloom’ subnetwork suggest that heterotrophs controlled the ultimate fate of organic matter and may have contributed to the length of the bloom itself. While many of the taxa observed in the ‘bloom’ subnetwork were not observed during previous blooms in the area (Fandino et al., 2001; Mayali et al., 2011; Wilson et al., 2021), they have been observed in other high organic matter environments. For example, three top module members (MM > 0.70) of the ‘bloom’ subnetwork (two that were not in the original 18-month time series) had 100% similarity to homologs found during the Tara Oceans cruise from either the surface or deep chlorophyll maximum in locations with particularly high primary production, chlorophyll a, and particulate organic carbon (OM-RGC.v1.039482917, OM-RGC.v1.001612655, OM-RGC.v1.001602887; Villar et al., 2018). The Ribosomal Database Project (RDP) (Wang et al., 2007) also refined the identity of a Rhodobacterales ASV (MM = 0.83) to either Lentibacter spp. or Litoreibacter spp. Lentibacter spp. have been found to specialize in dissolved organic matter transformations and are an effective competitor for many substrates found during blooms (Han et al., 2021). Polaribacter spp. were also identified as part of the ‘bloom’ subnetwork (MM = 0.45) and were detected during both a 2005 (Mayali et al., 2011) and 2019 (Wilson et al., 2021) Lingulodinium polyedra bloom off the Scripps Pier. Polaribacter spp. are known to compete with Lentibacter during blooms (Han et al., 2021), where they have been observed to degrade polysaccharides (Teeling et al., 2012). The importance of bacteria in controlling the fate of L. polyedra was further indicated by Marivivens spp. (MM = 0.81), which are found in the nutrient-rich phycosphere of dinoflagellates and diatoms (Fu et al., 2020) and can supply L. polyedra with vitamins or other ecological benefits (Cruz-Lopez and Maske, 2016; Fu et al., 2020). Finally, the ASV with the highest membership score (MM = 0.96) for the ‘bloom’ subnetwork had an RDP seqmatch score of 1.0 to the Rhodobacterales, Jannaschia cystaugens, also known as Thalasobacter stenotrophicus (Pujalte et al., 2005). This bacterium is known to promote cyst formation in other HAB species (e.g., Alexandrium) and thus plays a significant role in the disintegration of some blooms. Because of this role, it tends to be a dominant taxon during bloom peak and decay (Adachi et al., 1999; Adachi et al., 2004); and it was only present during the 2020 L. polyedra bloom in our time series. However, it should be noted that the environment also likely contributed to the break-up of the L. polyedra bloom. For example, temperatures decreased towards the end of the bloom (Figure 2A) and dinoflagellates often favor warmer stratified waters following relaxation of upwelling (Mantyla et al., 2008; Smayda and Trainer, 2010; McGowan et al., 2017).

The ‘bloom’ subnetwork also likely impacted the environment. For example, ‘bloom’ subnetwork ASVs were statistically linked to DO concentrations reaching an all-time low. Previous work has shown that dissolved organic carbon released from red-tide-forming algae can trigger the growth and respiration of heterotrophic bacteria (Wada et al., 2018) and that subsets of the microbial community correlate with community respiration (Wilson et al., 2018), DO (Beman et al., 2021), and carbon export (Guidi et al., 2016). In this study, the ‘bloom’ subnetwork’s inferred heterotrophic activities likely created circumstances that stimulated another bloom following the decline in Lingulodinium polyedra via a combination of high nutrients and open niches that other photosynthetic organisms were able to occupy. Many phytoplankton blooms are characterized by succession of several different phytoplankton species (Chang et al., 2003; Teeling et al., 2012; Shao et al., 2020), and in our study the secondary bloom was made up primarily of diatoms that began increasing after L. polyedra declined dramatically. Multiple analyses demonstrated this secondary diatom bloom, including the fact that: 1) chloroplasts contributed significantly to the novelty of post-bloom samples in the SOM analysis, 2) the ‘post-bloom’ subnetwork possessed many chloroplasts/Cyanobacteria with high MM, and 3) both diatoms and AF group 1 (which were likely diatoms) were high following the L. polyedra bloom. This specific taxonomic succession has been seen in other coastal areas (Shao et al., 2020). However, the opposite situation of dinoflagellates following diatoms has been observed most commonly (Smayda and Trainer, 2010), with dinoflagellates sometimes attributed with breaking up a diatom bloom (Tiselius and Kuylenstierna, 1996). In our study, however, the high concentration of diatoms (versus some other photosynthesizer seizing the available niches) in the post-HAB bloom was likely linked to the cooler water.

In addition to increased taxonomic novelty, there was increased community/subnetwork variance during and after the 2020 Lingulodinium polyedra bloom. Rather than the typical singular dominant (and two subdominant) subnetwork(s) from April through June of 2020, there were five different dominant subnetworks. This increased taxonomic variance followed clear succession as the system transitioned from a typical seasonal pattern (though the switch between the ‘fall/winter’ and ‘winter/spring’ subnetworks occurred later than usual) to ‘bloom’ conditions, ‘post-bloom’ conditions, and back to a typical seasonal pattern. The ‘bloom’/’post-bloom’ succession we observed is consistent with the literature that has shown succession linked to different bloom stages (Teeling et al., 2012; Teeling et al., 2016; Bunse and Pinhassi, 2017; Shao et al., 2020). The degree to which all other subnetworks decreased when the ‘bloom’ subnetwork initially increased shows the selective power of the 2020 L. polyedra bloom to initiate the two novel successional stages.

We hypothesize that our results also indicate increased variance in microbial functioning during and after the Lingulodinium polyedra bloom, which we base on correlations between subnetworks and measured oxygen cycling as well as gene predictions based on taxonomic identifications (i.e., measured genetic potential) rather than measured functional capacity. Nonetheless, despite low DO conditions throughout the ‘bloom’ and ‘post-bloom’ period, functional oxygen cycling was not constant and flickered between trophic states even as the ‘bloom’ and ‘post-bloom’ communities remained. This trophic state variability occurred to the point that the ‘post-bloom’ subnetwork did not correlate with DO or AOU. Many ecological studies have found that increased variance in ecological function (here, AOU) indicates a critical transition and overall regime shift (Carpenter and Brock, 2006; Scheffer et al., 2009; Scheffer et al., 2012; Wang et al., 2012). Though this concept has not been widely applied to natural microbial systems due to a lack of high frequency data that can capture variance on microbial timescales, one study using high frequency data from a marine lake saw ecosystem metabolism flicker as the epilimnion changed depth and a new normal was established (Wilson et al., 2019). If such a situation applies to our system, the L. polyedra bloom and the ‘bloom’ and ‘post-bloom’ subnetworks can be thought of as regime shifts on timescales relevant to microorganisms. Variables that could have contributed to this potential flickering include top-down forces (grazing and viral lysis) and bottom-up forces (nutrient or substrate limitation) which could have selected for different subsets and altered function (Wilhelm and Suttle, 1999; Chang et al., 2003; Buchan et al., 2014). However, purely physical forces (e.g., resuspension of organic matter or advection) or self-limitation by aerobic organisms may also have contributed or caused this flickering without increased variation in functional capacity, only in measured function. The flickering in oxygen cycling that was observed during and after the 2020 L. polyedra bloom persisted into the summer (with similar bloom values in July and August), well after the red tide event ended. Although microbial data after July 2, 2020, is not yet available, we hypothesize that portions of the ‘bloom’ and ‘post-bloom’ community may have increased again during the DO and AOU late-summer flickers or that the typical seasonal community was reestablishing and thus displayed functional variation.

Conclusions

Taken together, these data illustrate how the massive 2020 Lingulodinium polyedra bloom created unique ecological conditions for coastal Southern California. A very particular combination of bacterial and archaeal ASVs quickly responded before following a clear successional trajectory. Many of the ASVs that responded had not been detected previously in the area. However, these unique local ASVs have been recorded in other extremely high organic matter environments and following other HABs, demonstrating a specific set of functional life histories. The inferred function (based on taxonomic identities) of the responding microbes with regard to oxygen cycling likely depended on many different factors, including available substrates and other unassessed variables such as physical forcing, predators, or nutrients. Due to some combination of these factors, this function (and potentially functional capacity) flickered throughout the ‘bloom’ subnetwork’s dominance, as the system attempted to acclimate to the new normal. A similar situation occurred during the ‘post-bloom’ time period, where abundance and activity of the post-HAB diatom bloom and associated microorganisms were ultimately limited, either by the community’s own members or other external factors. Even as the ‘bloom’ subnetwork persisted, the system was in a state of previously unobserved flux such that measured function continued to flicker or shift in a novel way for the existing time series. The shift/flicker in measured function likely included changes in carbon and nutrient cycling as the bloom/post-bloom microbial succession and functional flickering indicated a substantial change in substrates (which many of the life histories of microbes in the relevant subnetworks also suggest). Through these data, we were able to show that the 2020 L. polyedra bloom elicited a novel microbial response taxonomically and (potentially) functionally. This work also gives insight into how microbes in coastal Southern California may respond to a similar event. Further studies regarding the effects of HABs can use these results as a springboard to target the microbial groups (e.g., bloom degraders like Thalasobacter stenotrophicus) and conditions that most influence ecosystem responses during such events.

Data accessibility statement

All data that were generated (sequences and flow cytometry data) are publically available. 16S rRNA gene sequences were submitted to NCBI SRA at BioProject PRJNA662174. Flow cytometry data are available through the UC San Diego Library Digital Collections: Bowman, JS, Wilson, J, Connors, B. 2021. Scripps Ecological Observatory. UC San Diego Library Digital Collections. https://doi.org/10.6075/J0348KHJ.

Phytoplankton counts, chlorophyll a concentration, and phaeophytin concentration came from Southern California Coastal Ocean Observing System (SCCOOS): https://sccoos.org/; DOI: 10.18436/S6159F.

All water temperature and oxygen data came from Scripps Ocean Acidification Real-Time (SOAR) Monitoring Program. 2022. https://coralreefecology.ucsd.edu/research/scrippsoceanacidificationreal-timesoarmonitoringprogram/.

Supplemental files

The supplemental files for this article can be found as follows:

Figures S1–S4. PDF

Table S1. xlsx

Acknowledgments

The authors would like to thank Avishek Dutta, Benjamin Klempay, Melissa Hopkins, Gabriella Berman, and Mackenzie Davey for their assistance with microbial and flow cytometry sample collection and processing. We would also like to thank Kristi Seech for her assistance with SCCOOS phytoplankton cell identification and enumeration and John A. McGowan for support and leadership over the decades with the McGowan Plankton and Chlorophyll Program. Finally, the authors would like to thank Sarah Abboud for her assistance with the manuscript.

Funding

The Southern California Coastal Ocean Observing Harmful Algal Bloom Monitoring Program was supported by the National Oceanic and Atmospheric Administration (NOAA NA16NOS0120022, NA11NOS120029, and NA17RJ1231). Funding for the McGowan Plankton and Chlorophyll Program was provided by private donors and the MacArthur Foundation. This work was also supported by a Simons Foundation Early Career Marine Microbial Ecology and Evolution award to JSB.

Competing interests

The authors report no competing interests. JSB is an associate editor at Elementa. He was not involved in the review process of this manuscript.

Author contributions

Completed field and laboratory work: JMW, NE, EC, EJC, SMC, MLC, JES, JSB.

Completed statistical/data analysis: JMW.

Wrote the paper: JMW, NE, EC, EJC, JSB.

Reviewed, revised, and approved the manuscript: All authors.

References

Adachi, M, Kanno T, Matsubara, T, Nishijima T, Itakura, S, Yamaguchi, M. 1999. Promotion of cyst formation in the toxic dinoflagellate Alexandrium (Dinophyceae) by natural bacterial assemblages from Hiroshima Bay, Japan. Marine Ecology Progress Series 191: 175–185. DOI: https://dx.doi.org/10.3354/meps191175.

Adachi, M, Kanno, T, Okamoto, R, Shinozaki, A, Fujikawa-Adachi, K, Nishijima, T. 2004. Jannaschia cystaugens sp. nov., an Alexandrium (Dinophyceae) cyst formation-promoting bacterium from Hiroshima Bay, Japan. International Journal of Systematic Evolutionary Microbiology 54(5): 1687–1692. DOI: https://dx.doi.org/10.1099/ijs.0.03029-0.

Allen, WE. 1943. “Red Water” in La Jolla Bay in 1942. Transactions of the American Microscopical Society 62(3): 262–264. DOI: https://dx.doi.org/10.2307/3223029.

Anderson, C, Hepner-Medina, M. 2020. SCCOOS red tide bulletin: Spring 2020. Available at https://sccoos.org/california-hab-bulletin/red-tide/. Accessed June 6, 2021.

Anderson, CR, Moore, SK, Tomlinson, MC, Silke, J, Cusack, CK. 2015. Living with harmful algal blooms in a changing world: Strategies for modeling and mitigating their effects in coastal marine ecosystems, in Ellis, JT, Sherman, DJ, Shroder, JF eds., Coastal and marine hazards, risks, and disasters. Waltham, MA: Elsevier: 495–561.

Anderson, DM. 1995. Toxic red tides and harmful algal blooms: A practical challenge in coastal oceanography. Reviews of Geophysics 33(S2): 1189–1200. DOI: https://dx.doi.org/10.1029/95RG00440.

Avci, B, Kruger, K, Fuchs, BM, Teeling, H, Amann, RI. 2020. Polysaccharide niche partitioning of distinct Polaribacter clades during North Sea spring algal blooms. The ISME Journal 14: 1369–1383. DOI: https://dx.doi.org/10.1038/s41396-020-0601-y.

Barbera, P, Kozlov, AM, Czech, L, Morel, B, Darriba, D, Flouri, T, Stamatakis, A. 2019. EPA-ng: Massively parallel evolutionary placement of genetic sequences. Systematic Biology 68(2): 365–369. DOI: https://dx.doi.org/10.1093/sysbio/syy054.

Beman, JM, Vargas, SM, Vazquez, S, Wilson, JM, Yu, A, Cairo, A, Perez-Coronel, E. 2021. Biogeochemistry and hydrography shape microbial community assembly and activity in the eastern tropical North Pacific Ocean oxygen minimum zone. Environmental Microbiology 23(6): 2765–2781. DOI: https://dx.doi.org/10.1111/1462-2920.15215.

Bowman, JS. 2022. Paprica lab protocol. Available at https://github.com/bowmanlab/seq_data_scripts/blob/master/dada2.r. Accessed February 12, 2021.

Bowman, JS, Amaral-Zettler, LA, Rich, JJ, Luria, CM, Ducklow, HW. 2017. Bacterial community segmentation facilitates the prediction of ecosystem function along the coast of the western Antarctic Peninsula. The ISME Journal 11(6): 1460–1471. DOI: https://dx.doi.org/10.1038/ismej.2016.204.

Bowman, JS, Ducklow, HW. 2015. Microbial communities can be described by metabolic structure: A general framework and application to a seasonally variable, depth-stratified microbial community from the coastal West Antarctic Peninsula. PLoS One 10(8): e0135868. DOI: https://dx.doi.org/10.1371/journal.pone.0135868.

Bowman, JS, Wilson, J, Connors, B. 2021. Scripps Ecological Observatory. UC San Diego Library Digital Collections. DOI: https://doi.org/10.6075/J0348KHJ.

Buchan, A, LeCleir, GR, Gulvik, CA, Gonzalez, JM. 2014. Master recyclers: Features and functions of bacteria associated with phytoplankton blooms. Nature Reviews Microbiology 12: 686–698. DOI: https://dx.doi.org/10.1038/nrmicro3326.

Bunse, C, Pinhassi, J. 2017. Marine bacterioplankton seasonal succession dynamics. Trends in Microbiology 25(6): 494–505. DOI: https://dx.doi.org/10.1016/j.tim.2016.12.013.

Button, DK, Robertson, BR, Lepp, PW, Schmidt, TM. 1998. A small, dilute-cytoplasm, high-affinity, novel bacterium isolated by extinction culture and having kinetic constants compatible with growth at ambient concentrations of dissolved nutrients in seawater. Applied Environmental Microbiology 64(11): 4467–4476. DOI: https://dx.doi.org/10.1128/aem.64.11.4467-4476.1998.

Callahan, BJ, McMurdie, PJ, Rosen, MJ, Han, AW, Johnson, AJA, Holmes, SP. 2016. DADA2: High-resolution sample inference from Illumina amplicon data. Nature Methods 13: 581–583. DOI: https://dx.doi.org/10.1038/nmeth.3869.

Carpenter, SR, Brock, WA. 2006. Rising variance: A leading indicator of ecological transition. Ecology Letters 9(3). DOI: https://dx.doi.org/10.1111/j.1461-0248.2005.00877.x.

Chang, FH, Zeldis, J, Gall, M, Hall, J. 2003. Seasonal and spatial variation of phytoplankton assemblages, biomass and cell size from spring to summer across the north-eastern New Zealand continental shelf. Journal of Plankton Research 25(7): 737–758. DOI: https://dx.doi.org/10.1093/plankt/25.7.737.

Chen, Y-J, Leung, PM, Bay, S, Hugenholtz, P, Kessler, A, Shelley, G, Waite, D, Cook, P, Greening, C. 2020. Metabolic flexibility allows generalist bacteria to become dominant in a frequently disturbed ecosystem. bioRxiv. DOI: https://dx.doi.org/10.1101/2020.02.12.945220.

Cruz-López, R, Maske, H. 2016. The vitamin B1 and B12 required by the marine dinoflagellate Lingulodinium polyedrum can be provided by its associated bacterial community in culture. Frontiers in Microbiology 7(May): 560. DOI: https://dx.doi.org/10.3389/fmicb.2016.00560.

Czech, L, Barbera, P, Stamatakis, A. 2020. Genesis and Gappa: Processing, analyzing and visualizing phylogenetic (placement) data. Bioinformatics 36(10): 3263–3265. DOI: https://dx.doi.org/10.1093/bioinformatics/btaa070.

Edwards, JL, Smith, DL, Connolly, J, McDonald, JE, Cox, MJ, Joint, I, Edwards, C, McCarthy, AJ. 2010. Identification of carbohydrate metabolism genes in the metagenome of a marine biofilm community shown to be dominated by Gammaproteobacteria and Bacteroidetes. Genes 1(3): 371–384. DOI: https://dx.doi.org/10.3390/genes1030371.

Fegatella, F, Lim, J, Kjelleberg, S, Cavicchioli, R. 1998. Implications of rRNA operon copy number and ribosome content in the marine oligotrophic ultramicrobacterium Sphingomonas sp. strain RB2256. Applied and Environmental Microbiology 64(11): 4433–4438. DOI: https://dx.doi.org/10.1128/AEM.64.11.4433-4438.1998.

Fandino, LB, Riemann, L, Steward, GF, Long, RA, Azam, F. 2001. Variations in bacterial community structure during a dinoflagellate bloom analyzed by DGGE and 16S rDNA sequencing. Aquatic Microbial Ecology 23(2): 119–130. DOI: https://dx.doi.org/10.3354/ame023119.

Ferrer-Gonzalez, FX, Widner, B, Holderman, NR, Glushka, J, Edison, AS, Kujawinski, EB, Moran, MA. 2020. Resource partitioning of phytoplankton metabolites that support bacterial heterotrophy. The ISME Journal 15(3): 762–773. DOI: https://dx.doi.org/10.1038/s41396-020-00811-y.

Fu, H, Uchimiya, M, Gore, J, Moran, MA. 2020. Ecological drivers of bacterial community assembly in synthetic phycospheres. Proceedings of National Academy of Sciences USA 117(7): 3656–3662. DOI: https://dx.doi.org/10.1073/pnas.1917265117.

Gobler, CJ. 2020. Climate change and harmful algal blooms: Insights and perspective. Harmful Algae 91: 101731. DOI: https://dx.doi.org/10.1016/j.hal.2019.101731.

Gómez-Pereira, PR, Schüler, M, Fuchs, BM, Bennke, C, Teeling, H, Waldmann, J, Richter, M, Barbe, V, Bataille, E, Glöckner, FO, Amann, R. 2012. Genomic content of uncultured Bacteroidetes from contrasting oceanic provinces in the North Atlantic Ocean. Environmental Microbiology 14(1): 52–66. DOI: https://dx.doi.org/10.1111/j.1462-2920.2011.02555.x.

Guidi, L, Chaffron, S, Bittner, L, Eveillard, D, Larhlimi, A, Roux, S, Darzi, Y, Audic, S, Berline, L, Brum, JR, Coelho, LP, Tara Oceans Consortium Coordinators, Stemmann, L, Not, F, Hingamp, P, Speich, S, Follows, M, Karp-Boss, L, Boss, E, Ogata, H, Pesant, S, Weissenbach, J, Wincker, P, Acinas, SG, Bork, P, de Vargas, C, Iudicone, D, Sullivan, MB, Raes, J, Karsenti, E, Bowler, C, Gorsky, G. 2016. Plankton networks driving carbon export in the oligotrophic ocean. Nature 532(7600): 465–470. DOI: https://dx.doi.org/10.1038/nature16942.

Haft, DH, DiCuccio, M, Badretdin, A, Brover, V, Chetvernin, V, O’Neill, K, Li, W, Chitsaz, F, Derbyshire, MK, Gonzales, NR, Gwadz, M, Lu, F, Marchler, GH, Song, JS, Thanki, N, Yamashita, RA, Zheng, C, Thibaud-Nissen, F, Geer, LY, Marchler-Bauer, A, Pruitt, KD. 2018. RefSeq: An update on prokaryotic genome annotation and curation. Nucleic Acids Research 46(D1): D851–D860. DOI: https://dx.doi.org/10.1093/nar/gkx1068.

Han, Y, Jiao, N, Zhang, Y, Zhang, F, He, C, Liang, X, Cai, R, Shi, Q, Tang, K. 2021. Opportunistic bacteria with reduced genomes are effective competitors for organic nitrogen compounds in coastal dinoflagellate blooms. Microbiome 9(1). DOI: https://dx.doi.org/10.1186/s40168-021-01022-z.

Hatch, MBA, Schellenberg, SA, Carter, ML. 2013. Ba/Ca variations in the modern intertidal bean clam Donax gouldii: An upwelling proxy? Palaeogeography, Palaeoclimatology, Palaeoecology 373: 98–107. DOI: https://dx.doi.org/10.1016/j.palaeo.2012.03.006.

Holmes, RW, Williams, PM, Eppley, RW. 1967. Red water in La Jolla Bay, 1964–1966. Limnology and Oceanography 12(3): 503–512. DOI: https://dx.doi.org/10.4319/lo.1967.12.3.0503.

Howard, MDA, Silver, M, Kudela, RM. 2008. Yessotoxin detected in mussel (Mytilus californicus) and phytoplankton samples from the U.S. west coast. Harmful Algae 7(5): 646–652. DOI: https://dx.doi.org/10.1016/j.hal.2008.01.003.

Kahru, M, Mitchell, BG. 1998. Spectral reflectance and absorption of a massive red tide off southern California. Journal of Geophysical Research Ocean 103(C10): 21601–21609. DOI: https://dx.doi.org/10.1029/98JC01945.

Kirchman, DL. 2002. The ecology of Cytophaga-Flavobacteria in aquatic environments. FEMS Microbiology Ecology 39(2): 91–100. DOI: https://dx.doi.org/10.1111/j.1574-6941.2002.tb00910.x.

Klappenbach, JA, Dunbar, JM, Schmidt, TM. 2000. rRNA operon copy number reflects ecological strategies of bacteria. Applied Environmental Microbiology 66(4): 1328–1333. DOI: https://dx.doi.org/10.1128/AEM.66.4.1328-1333.2000.

Kudela, RM, Cochlan, WP. 2000. Nitrogen and carbon uptake kinetics and the influence of irradiance for a red tide bloom off southern California. Aquatic Microbial Ecology 21(1): 31–47. DOI: https://dx.doi.org/10.3354/ame021031.

Landa, M, Burns, AS, Roth, SJ, Moran, MA. 2017. Bacterial transcriptome remodeling during sequential co-culture with a marine dinoflagellate and diatom. The ISME Journal 11: 2677–2690. DOI: https://dx.doi.org/10.1038/ismej.2017.117.

Langfelder, P, Horvath, S. 2007. Eigengene networks for studying the relationships between co-expression modules. BMC Systems Biology 1: 54. DOI: https://dx.doi.org/10.1186/1752-0509-1-54.

Langfelder, P, Horvath, S. 2008. WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics 9: 59. DOI: https://dx.doi.org/10.1186/1471-2105-9-559.

Lewitus, AJ, Horner, RA, Caron, DA, Garcia-Mendoza, E, Hickey, BM, Hunter, M, Huppert, DD, Kudela, RM, Langlois, GW, Largier, JL, Lessard, EJ, RaLongde, R, Rensel, JE, Strutton, PG, Trainer, VL, Tweddle, JF. 2012. Harmful algal blooms along the North American west coast region: History, trends, causes, and impacts. Harmful Algae 19: 133–159. DOI: https://dx.doi.org/10.1016/j.hal.2012.06.009.

Mantyla, AW, Bograd, SJ, Venrick, EL. 2008. Patterns and controls of chlorophyll-a and primary productivity cycles in the Southern California Bight. Journal of Marine Systems 73(1–2): 48–60. DOI: https://dx.doi.org/10.1016/j.jmarsys.2007.08.001.

Matcher, G, Lemley, D, Adams, J. 2021. Bacterial community dynamics during a harmful algal bloom of Heterosigma akashiwo. Aquatic Microbial Ecology 86: 153–167. DOI: https://dx.doi.org/10.3354/ame01963.

Mayali, X, Franks, PJS, Burton, RS. 2011. Temporal attachment dynamics by distinct bacterial taxa during a dinoflagellate bloom. Aquatic Microbial Ecology 63(2): 111–122. DOI: https://dx.doi.org/10.3354/ame01483.

McGowan, JA, Deyle, ER, Ye, H, Carter, ML, Perretti, CT, Seger, KD, de Verneil, A, Sugihara, G. 2017. Predicting coastal algal blooms in southern California. Ecology 98(5): 1419–1433. DOI: https://dx.doi.org/10.1002/ecy.1804.

Moorthi, SD, Countway, PD, Stauffer, BA, Caron, DA. 2006. Use of quantitative real-time PCR to investigate the dynamics of the red tide dinoflagellate Lingulodinium polyedrum. Microbial Ecology 52(1): 136–150. DOI: https://dx.doi.org/10.1007/s00248-006-9030-3.

Nawrocki, EP, Eddy, SR. 2013. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics 29(22): 2933–2935. DOI: https://dx.doi.org/10.1093/bioinformatics/btt509.

Paz, B, Riobó, P, Fernández, ML, Fraga S, Franco, JM. 2004. Production and release of yessotoxins by the dinoflagellates Protoceratium reticulatum and Lingulodinium polyedrum in culture. Toxicon 44(3): 251–258. DOI: https://dx.doi.org/10.1016/j.toxicon.2004.05.021.

Pinhassi, J, Sala, MM, Havskum, H, Peters, F, Guadayol, Ò, Malits, A, Marrasé, C. 2004. Changes in bacterioplankton composition under different phytoplankton regimens. Applied Environmental Microbiology 70(11): 6753–6766. DOI: https://dx.doi.org/10.1128/AEM.70.11.6753-6766.2004.

Pujalte, MJ, Macián, MC, Arahal, DR, Garay, E. 2005. Thalassobacter stenotrophicus Macián et al. 2005 is a later synonym of Jannaschia cystaugens Adachi et al. 2004, with emended description of the genus Thalassobacter. International Journal of Systematic and Evolutionary Microbiology 55(5): 1959–1963. DOI: https://dx.doi.org/10.1099/ijs.0.63617-0.

Scheffer, M, Bascompte, J, Brock, WA, Brovkin, V, Carpenter, SR, Dakos, V, Held, H, van Nes, EH, Rietkerk, M, Sugihara, G. 2009. Early-warning signals for critical transitions. Nature 461(7260): 53–59. DOI: https://dx.doi.org/10.1038/nature08227.

Scheffer, M, Carpenter, SR, Lenton, TM, Bascompte, J, Brock, W, Dakos, V, van de Koppel, J, van de Leemput, IA, Levin, SA, van Nes, EH, Pascual, M, Vandermeer, J. 2012. Anticipating critical transitions. Science 338(6105): 344–348. DOI: https://dx.doi.org/10.1126/science.1225244.

Scripps Ocean Acidification Real-Time (SOAR) Monitoring Program. 2022. Available at https://coralreefecology.ucsd.edu/research/scrippsoceanacidificationreal-timesoarmonitoringprogram/. Accessed April 15, 2021.

Shao, Q, Lin, Z, Zhou, C, Zhu, P, Yan, X. 2020. Succession of bacterioplankton communities over complete Gymnodinium-diatom bloom cycles. Science of the Total Environment 709: 135951. DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.135951.

Smayda, TJ, Trainer, VL. 2010. Dinoflagellate blooms in upwelling systems: Seeding, variability, and contrasts with diatom bloom behaviour. Progress in Oceanography 85(1–2): 92–107. DOI: https://dx.doi.org/10.1016/j.pocean.2010.02.006.

Southern California Coastal Ocean Observing System (SCCOOS). 2002. https://sccoos.org/; DOI: https://dx.doi.org/10.18436/S6159F. Accessed June 12, 2021.

Teeling, H, Fuchs, BM, Becher, D, Klockow, C, Gardebrecht, A, Bennke, CM, Kassabgy, M, Huang, S, Mann, AJ, Waldmann, J, Weber, M, Klindworth, A, Otto, A, Lange, J, Bernhardt, J, Reinsch, C, Hecker, M, Peplies, J, Bockelmann, FD, Callies, U, Gerdts, G, Wichels, A, Wiltshire, KH, Glöckner, FO, Schweder, T, Amann, R. 2012. Substrate-controlled succession of marine bacterioplankton populations induced by a phytoplankton bloom. Science 336(6081): 608–611. DOI: https://dx.doi.org/10.1126/science.1218344.

Teeling, H, Fuchs, BM, Bennke, CM, Krüger, K, Chafee, M, Kappelmann, L, Reintjes, G, Waldmann, J, Quast, C, Glöckner, FO, Lucas, J. 2016. Recurring patterns in bacterioplankton dynamics during coastal spring algae blooms. eLife 5: e11888. DOI: https://dx.doi.org/10.7554/eLife.11888.

Thomas, F, Hehemann, J-H, Rebuffet, E, Czjzek, M, Michel, G. 2011. Environmental and gut Bacteroidetes: The food connection. Frontiers in Microbiology 2(May): 93. DOI: https://dx.doi.org/10.3389/fmicb.2011.00093.

Tiselius, P, Kuylenstierna, M. 1996. Growth and decline of a diatom spring bloom: Phytoplankton species composition, formation of marine snow and the role of heterotrophic dinoflagellates. Journal of Plankton Research 18(2): 133–155. DOI: https://dx.doi.org/10.1093/plankt/18.2.133.

Villar, E, Vannier, T, Vernette, C, Lescot, M, Cuenca, M, Alexandre, A, Bachelerie, P, Rosnet, T, Pelletier, E, Sunagawa, S, Hingamp, P. 2018. The Ocean Gene Atlas: Exploring the biogeography of plankton genes online. Nucleic Acids Research 46(W1): W289–W295. DOI: https://dx.doi.org/10.1093/nar/gky376.

Wada, M, Takano, Y, Nagae, S, Ohtake, Y, Umezawa, Y, Nakamura, S, Yoshida, M, Matsuyama, Y, Iwataki, M, Takeshita, S, Oda, T. 2018. Temporal dynamics of dissolved organic carbon (DOC) produced in a microcosm with red tide forming algae Chattonella marina and its associated bacteria. Journal of Oceanography 74(6): 587–593. DOI: https://dx.doi.org/10.1007/s10872-017-0455-8.

Walters, W, Hyde, ER, Berg-Lyons, D, Ackermann, G, Humphrey, G, Parada, A, Gilbert, JA, Jansson, JK, Caporaso, JG, Fuhrman, JA, Apprill, A. 2016. Improved bacterial 16S rRNA gene (V4 and V4-5) and fungal internal transcribed spacer marker gene primers for microbial community surveys. mSystems 1(1): e00009–e000015. DOI: https://dx.doi.org/10.1128/msystems.00009-15.

Wang, R, Dearing, JA, Langdon, PG, Zhang, E, Yang, X, Dakos, V, Scheffer, M. 2012. Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature 492(7429): 419–422. DOI: https://dx.doi.org/10.1038/nature11655.

Wang, QG, Garrity, M, Tiedje, JM, Cole, JR. 2007. Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Applied Environmental Microbiology 73(16): 5261–5267. DOI: https://dx.doi.org/10.1128/AEM.00062-07.

Wear, EK, Carlson, CA, James, AK, Brzezinski, MA, Windecker, LA, Nelson, CE. 2015. Synchronous shifts in dissolved organic carbon bioavailability and bacterial community responses over the course of an upwelling-driven phytoplankton bloom. Limnology and Oceanography 60(2): 657–677. DOI: https://dx.doi.org/10.1002/lno.10042.

Wehrens, R, Kruisselbrink, J. 2018. Flexible self-organizing maps in kohonen 3.0. Journal of Statistical Software 87(7): 1–18. DOI: https://dx.doi.org/10.18637/jss.v087.i07.

Wilhelm, SW, Suttle, CA. 1999. Viruses and nutrient cycles in the sea: Viruses play critical roles in the structure and function of aquatic food webs. Bioscience 49(10): 781–788. DOI: https://dx.doi.org/10.2307/1313569.

Wilson, J, Ucharm, G, Beman, JM. 2019. Climatic, physical, and biogeochemical changes drive rapid oxygen loss and recovery in a marine ecosystem. Scientific Reports 9(1): 16114. DOI: https://dx.doi.org/10.1038/s41598-019-52430-z.

Wilson, JM, Carter, ML, Mühle, J, Bowman, JS. 2020. Using empirical dynamic modeling to assess relationships between atmospheric trace gases and eukaryotic phytoplankton populations in coastal Southern California. Marine Chemistry 227: 103896. DOI: https://dx.doi.org/10.1016/j.marchem.2020.103896.

Wilson, JM, Chamberlain, EJ, Erazo, N, Carter, ML, Bowman, JS. 2021. Recurrent microbial community types driven by nearshore and seasonal processes in coastal Southern California. Environmental Microbiology 23(6): 3225–3239. DOI: https://dx.doi.org/10.1111/1462-2920.15548.

Wilson, JM, Litvin, SY, Beman, JM. 2018. Microbial community networks associated with variations in community respiration rates during upwelling in nearshore Monterey Bay, California. Environmental Microbiology Reports 10(3): 272–282. DOI: https://dx.doi.org/10.1111/1758-2229.12635.