Harmful algal blooms are preceded by a predictable and quantifiable shift in the oceanic microbiome

Domain: earth_science
Published: 2025-04-29
Editorial outcome: major_revision

Abstract
Harmful algal blooms (HABs) have become a worldwide environmental and human health problem, stressing the urgent need for a reliable forecasting tool. Dynamic interactions between algae, including harmful algae, and bacteria play a large role regulating water chemistry. Free-living bacteria quickly respond to small physical and/or chemical environmental changes by adjusting their proteome. We hypothesize that this response is detectable at the peptide level and occurs before rapid phytoplankton growth characteristic of harmful bloom events. To characterize the microbiome’s physiological changes preceding bloom onset, we collected and analyzed a high-resolution metaproteomic time series of a free-living microbiome in a coastal ecosystem. We confirm that twelve candidate HAB biomarkers are detectable, quantifiable, and correlated across two pre-bloom periods. This study identifies proteomic shifts in bacterial peptides which may be used as predictive biomarkers for forecasting harmful algal bloom initiation, potentially mitigating detrimental algal bloom outcomes in the future.

Introduction
Harmful algal blooms (HABs) are sudden accumulations of photosynthetic phytoplankton that cause adverse effects to aquatic ecosystems1. These blooms pose serious threats to both human health and the economy2 with negative consequences on fishing, tourism, and real estate industries as well as mass mortalities of coastal marine mammals, seabirds, and fishes. Harmful effects caused by bloom-forming species are primarily attributed to the production of highly toxic compounds such as neurotoxins3; that said, there are many non-toxic HAB species that achieve remarkably high cell densities that can block sunlight to deeper waters, alter local water chemistry (e.g., pH, nutrients)4, generate hypoxic or anoxic zones5, and/or obstruct fish gills, all potentially leading to a swift rise in marine animal mortalities6.

A variety of human-made and environmental perturbations are capable of initiating these blooms7; however, in the majority of bloom events, the specific underlying cause for HAB initiation is unknown. In general, two fundamental goals of HAB research are (1) to determine what factors lead to the dominance of HAB-forming species and (2) to predict bloom onset, enabling warning and mitigation strategies to combat the increasing frequency and intensity of HAB events globally8.

Interactions between phytoplankton and bacteria have been shown to play a paramount role in shaping the dynamics of nutrient availability, biogeochemical cycles, and microbial diversity and function in marine ecosystems9–11. The individual and interacting factors influencing HAB initiation, including chemical, physical, and biological aspects of an ecosystem, are experienced by both HAB-forming phytoplankton as well as their bacterial neighbors. In the specific context of HABs, tight associations have been observed within the microbiome as phytoplankton advance through bloom cycles, with bacteria playing a deterministic role in the outcome of HAB progression and algal byproducts12–14. Additionally, certain consistently abundant taxa have been noted as commonly associated with different phytoplankton species15–17, lending precedent to the idea that community structure and function are driven by specific associations18. The influence these organisms have on each other’s biological potential creates a recognized opportunity19 to utilize bacterial responses as a rapid in situ predictor for phytoplankton dynamics, regardless of the causative influence that may trigger bloom initiation.

Proteomic techniques are capable of capturing this active functionality of cells through the measurement of peptides, the component pieces of proteins that are tightly regulated and rapidly modified in response to minute environmental perturbations20,21. Recent advances in mass spectrometry (MS) techniques and methods have enabled highly accurate quantification of peptides in a variety of sample types across a wide dynamic range22. Additionally, as the development of peptide-based biomarkers for diagnosing human health23 is gaining momentum in medical research, a range of potential molecular or antigen rapid tests could be deployed via citizen science24 or autonomous vehicle monitoring for environmental applications25. While proteomic analysis of eukaryotic organisms like phytoplankton remain difficult26, a predictive biomarker developed using a bacteria-centric approach would potentially eliminate the need for targeted forecasting due to taxonomic and functional overlap of marine microbiomes globally27,28.

Here we present a high-resolution time series that captures the biological and chemical characteristics of a dynamic coastal ecosystem during two pre-bloom periods of Chaetoceros socialis, a cosmopolitan and abundant HAB-forming diatom capable of deteriorating water quality and killing fish through physical irritation29. To interrogate the microbiome for biomarkers capable of predicting a HAB >24 hours in advance, we perform a metaproteomic analysis that quantifies bacterial peptides collected every four hours. A multistep framework is developed for selecting and statistically validating potential biomarker peptides using a hidden Markov model (HMM), a statistical approach that can identify non-reverting state shifts in time series data. Our investigation reveals 314 potential biomarkers for predicting bloom onset. Targeted validation of peptides from a second pre-bloom period with a different chemical profile confirms a subset of 12 peptides as quantifiable, non-reverting biomarkers. These findings provide evidence towards the use of the bacterial microbiome to generate biomarkers for predicting future HABs.

Results
Field-sampling captures two pre-bloom initiation periods with divergent nutrient profiles
A high-resolution time series of a coastal bacterial microbiome was acquired from May 27 to June 18, 2021, in East Sound, WA, USA, capturing samples every 4 hours across one full HAB cycle and the pre-bloom dynamics of a second HAB (Fig. 1a; Supplementary Fig. 1; Supplementary Data 1). On-site microscopy determined that the eukaryotic fraction of the ecosystem was dominated by Chaetoceros socialis (SoundToxins Phytoplankton Monitoring Network)30. Two bloom initiations of C. socialis were identified from the detailed chlorophyll a fluorescence trace (Fig. 1b, c). Calculated bloom start dates (BSD)31, defined as the time point with the highest positive change in chlorophyll a per change in time preceding the bloom peak, were May 31 at 23:00 and June 18 at 11:00 (Supplementary Fig. 2). Cell counts of nanoeukaryotes and picoeukaryotes confirmed the May 31st BSD30 (Fig. 1d), while Sound Toxin microscopy and satellite images confirmed the presence of a HAB on June 19 [Sound Toxins Phytoplankton Monitoring Network]30 (Fig. 2). The nutrient conditions over the two pre-bloom initiation periods revealed divergent chemistries (Fig. 1h, i). Specifically, the pre-bloom period in May experienced a 20-30 mM increase in silicate and a 15 mM increase in nitrate ~60 hours in advance of the BSD (Fig. 1h), coinciding with a 4°C fluctuation in temperature (Fig. 1j). In contrast, in the pre-bloom conditions prior to the June BSD concentrations of all nutrients measured fluctuated <5 mM, revealing a more stable chemical ecosystem (Fig. 1i).

As a result of these divergent nutrient profiles in the pre-bloom conditions, the two time series were considered as separate in situ pre-bloom experiments (discovery and validation), allowing the opportunity to characterize each ecosystem independently (Fig. 1a). June 12 to 18 was selected as the discovery period because the environmental parameters monitored were more stable (Fig. 1f). May 28 to June 2 was selected as the validation period to test candidate biomarkers (Fig. 1a) to confirm that peptides remained non-reverting at least 24 hours post-BSD.

Discovery phase: identifying peptide biomarkers through quantitative metaproteomics
A metaproteomic pipeline was utilized to 1) confidently identify peptides from this complex microbiome, 2) generate quantitative peptide data, and 3) annotate peptides using a lowest common ancestor approach. Briefly, data-dependent acquisition (DDA) MS data from the June 12 to 18 sample set (36 samples) were searched against the metagenome-generated database (DB) of 3,457,818 proteins (PRIDE identifier PXD052873). DDA analysis yielded 36,429 protein detections from the 9765 peptides confidently detected, which were used to generate the smaller DDA-refined database for quantitative data independent acquisition (DIA) analysis and interpretation (Supplementary Fig. 3; Supplementary Data 2)32,33. The resulting DDA-refined database represents 1% of the original metagenome. Next, the DIA analysis of the discovery sample set detected 4645 unique peptides. A conservative approach was taken to only include peptides that yielded a quantitative value for >50% of the 36 time points, yielding a final set of 4468 peptides for the next phase of our biomarker screen (Fig. 3; Supplementary Data 3).

For Phase 1A, a principal component analysis (PCA) was performed on the time series to identify how the quantified peptides varied by time and day of collection. A PCA of all quantified peptides from the DIA analysis of the discovery period revealed that variation across PC1 (proportion of variance = 12.67%) was driven by absence of peptide signals across the time series (Supplementary Fig. 4). Visualization of PC2 (proportion of variance = 9.17%) vs PC3 (proportion of variance = 6.67%) revealed two clusters of data separating along PC2 (Fig. 4a). Analysis of the date and time of collection represented by each data point revealed that the separation was driven by the date of collection (between June 15-16; Fig. 4a), not the time of day of collection as expected. A Mann-Whitney test of PC2 loadings between June 12 to 15 and June 16 to 18 was significant (p = 1.916e-09). A total of 447 peptides, representing the top 10% of peptides driving variation across PC2, were selected for Phase 2 analysis (Supplementary Data 4).

In Phase 1B, 2015 of the total 4468 peptides were annotated with MetaGOmics34, revealing representation of 29 total bacterial classes in the quantitative peptide dataset (Supplementary Fig. 5; Supplementary Data 5). Nearly half of those annotated peptides were from Alphaproteobacteria (956 peptides). The next six most abundant classes were Flavobacteria (431 peptides), Gammaproteobacteria (339 peptides), Verrucomicrobiae (91 peptides), Bacilli (75 peptides), Opitutae (39 peptides), and Betaproteobacteria (14 peptides) (Supplementary Table 1). Seven of the 29 taxa showed significant correlation between summed quantitative peptide abundance values and time over the pre-bloom period (Spearman’s rho cut-off = |0.5|; Fig. 4c; Supplementary Table 2). Verrucomicrobiae, Opitutae, Acidobacteriia (1 peptide), and Zetaproteobacteria (1 peptide) exhibited a strong positive trend over time (⍴≥0.7, Fig. 4c), with a distinct inflection in abundance at the same June 15-16 changepoint identified from the PCA (Fig. 4a, c). The peptide trend from the most abundant class, Alphaproteobacteria, exhibited a weak positive correlation (⍴=0.35, Fig. 4c), while peptide trends from Gammaproteobacteria, Flavobacteria, and Epsilonproteobacteria (6 peptides) exhibited a negative correlation with time (⍴≤-0.5, Fig. 4c). All peptides from classes with a ⍴≥|0.5| were selected for the next phase of biomarker discovery. Additionally, peptides from the most abundant class, Alphaproteobacteria, were also selected to proceed through biomarker development. In total, our taxa-driven discovery approach yielded 1858 peptides for consideration as biomarkers in Phase 2.

As our goal was to identify biomarkers to predict HAB initiation, a hidden Markov model (HMM) approach was applied to reveal peptides with a non-reverting shift in underlying state processes. A two-state HMM was fit to each of the 2,015 Phase 1 candidate peptides’ temporal abundances (Fig. 3, yellow). A total of 314 peptides were confirmed to undergo a non-reverting state shift >24 hours prior to the bloom start date (Fig. 4b; Supplementary Data 6). This set included 88 peptides that were identified from both Phase 1 analyses, as well as 181 that originated from only the Phase 1A PCA and 221 from only the Phase 1B correlation analysis (Fig. 3).

From the 314 peptides that passed Phase 2, 87 peptides were randomly selected to test for a quantitative response in a dilution series using targeted mass spectrometry35. Because each microbiome sample carries a complex mixture of thousands of peptides, determining the quantity of any individual peptide abundance within the sample can be difficult. In order to ensure that the selected biomarker peptides are both detectable and quantifiable when present in complex mixtures, we created serial dilutions of the sample to generate a calibration curve prior to further analysis (Fig. 2, blue). Because there are interfering ions simultaneously eluting from the chromatography column prior to MS analyses, which can influence detectability of ions and the resulting quantities reported, these serial dilutions were generated with a complex suite of peptides, rather than water. In this dilution series, the total amount of protein injected for each MS run remained the same while the ratio of microbiome sample to diluent was varied incrementally. The diluent used here consisted of a protein digest of Colwellia psychrerythraea 34H, a model bacterium with no shared endogenous peptides. By diluting our test microbiome with peptides of interest against a background of equivalent biochemical complexity with no shared endogenous peptides, we were able to determine if candidate biomarkers exhibited a linear response during targeted MS experiments and establish a limit of detection (LOD) for each peptide. This method, the matrix-matched calibration curve (MMCC), is employed for peptide biomarker development in complex systems35. Samples prepared with the MMCC method were analyzed using parallel reaction monitoring (PRM) mass spectrometry to determine the LOD for 87 peptides down to detectable threshold of 1e3 and to assess whether each peptide’s dilution series exhibited a linear response above this threshold (Supplementary Fig. 6). Due to sample limitation, a random subset of 27 peptides that passed the MMCC standards were selected to progress through to the validation phase.

Validation of biomarkers as predictors of HABs
In order to determine whether our final subset of 27 biomarkers could predict a HAB, targeted MS analyses were performed on the May 28-June 2 pre-bloom sample set (i.e., the validation period, Fig. 1c). The goal of the validation experiment was to evaluate whether the 27 biomarkers identified in the discovery bloom dataset also exhibited a non-reverting state-shift prior to the June 18 BSD and whether the resulting temporal patterns were consistent across multiple blooms. All peptides identified as potential pre-bloom biomarkers were above the desired limit of detection (LOD) in the validation time series. A two-state HMM was fit to each peptide, revealing a total of 12 peptides that exhibited a non-reverting state shift prior to the BSD (Fig. 5a, b; Supplementary Figs. 7-8). The state shift for all peptides occurred at least 24 hours prior to the observed BSD for both blooms (Figs. 1c, 5b; Supplementary Figs. 7-8; Supplementary Table 3). Additionally, a Spearman correlation analysis of the discovery and validation blooms showed a strong positive correlation for eight of the quantified peptides when peptide abundance was aligned to the BSD (correlation value≥0.6 (e.g., Fig. 5a, b; Supplementary Fig. 9; Supplementary Data 7). A detailed annotation of the peptides with Basic Local Alignment Search Tool (BLAST) protein analysis revealed that, of the final 12 validated biomarkers, 8 peptides were from Gammaproteobacteria, 2 peptides were from Flavobacteriia, and 2 peptides were from Verrucomicrobiae. Functionally, 6 of the biomarker peptides were TonB-dependent receptors, 2 peptides were ribosomal, and others are indicative of an active response to the environment (Table 1; Supplementary Data 8).

Methods
Sample collection and water column analysis
Sample collection methods are detailed in Nunn et. al., 202430. Briefly, field sampling was conducted in East Sound, WA, USA from May 27 to June 18, 2021, by pumping water from 100 m offshore at a constant depth of 2 m and flow rate of 5 L min−1 to obtain 20 L. From the 20 L, multiple peristaltic pumps were set up to simultaneously filter between 500–2000 mL total volume per sample, averaging 1 L per microbiome filter. Whole water was filtered through a 200 µm nylon mesh followed by three subsequent filters at 100 µm, 10 µm, 1 µm before final collection onto a 0.22 µm polyethersulfone (PES) 47 mm filter. Final metaproteomics samples were stored long term in -80 °C on 0.22 µm polyethersulfone (PES) filters. Additional samples collected included nutrients, flow cytometry, prokaryotic DNA/proteomics, eukaryotic DNA/RNA/proteomics, and dissolved metabolites. Water chemistry was also monitored at the collection site with a YSI Exo1 Sonde probe capturing chlorophyll a concentration, salinity, temperature, dissolved oxygen, and pH every 10 minutes.

Metagenome assembly and database generation
DNA extraction was performed according to methods detailed in Nunn et al.30 on samples from all time points collected. DNA from all 132 samples were combined to generate a pooled sample with equal amounts of DNA (50 ng) from each sample for whole genome sequencing (WGS) using 150 bp paired-end Illumina sequencing (NextSeq 500 Mid Output v2 Kit-300 cycle) through the Northwest Genomics Center (NWGC) at the University of Washington. DNA library preparation, quality controls, and WGS were performed by the NWGC. The resulting FASTQs were assembled into a metagenome using MEGAHIT76. Samtools77 1.14 was used with Python78 3.7.7 to filter only mapped reads, and Prodigal79 2.6.3 was used to convert contigs into a FASTA file consisting of 3,457,818 proteins (PRIDE identifier PXD052873).

Proteomic sample processing
Samples collected for microbiome proteomic analysis (collected on 0.22 µm PES) were manually agitated within a 3×4 cm Ziplock baggie with 100 µl S-Trap solubilization buffer (5% SDS, 50 mM TEAB, 2 mM MgCl2, 1X HALT protease and phosphatase inhibitors in nanopure water) to ensure equal spreading of the buffer across the filter surface. The cell lysate was collected via pipette into a 2 mL microfuge tube. This was followed by an additional rinse of the filter with 100 µl nanopure water, followed by manual agitation and immediate transfer of the solution to the pooled fractions in the 2 mL tube. Two additional rinses of the filter with 100 µl nanopure water without manual agitation were performed and the solution added to the pool. The resulting 400 µl whole cell solution was then sonicated (Branson 250 Sonifier; 20 kHz, 3×10 s on ice). Samples were then evaporated in a SpeedVac (Thermo Fisher Scientific SRF110) at room temperature to reduce the volume to 100 µL (5% SDS). Protein concentrations for each sample were quantified in technical triplicates using a Bicinchoninic Acid protein assay kit (Pierce Thermo Scientific) using a microplate reader. At total of 20 µg protein per sample was isolated for enzymatic digestion and MS-based identification. Enolase protein standard (0.16 µL of 100 ng µl−1 enolase per 1 µg protein) was added to each sample to track digestion efficiency. Each sample (20 µg protein) was treated with benzonase to degrade DNA (0.5 µL 250 unit µl−1 10 minutes at 95 °C), and proteins were reduced with 20 mM dithiothreitol for 10 minutes at 60 °C followed by a four minute cool down to room temperature, alkylated with 40 mM iodoacetamide for 30 minutes in the dark, acidified to pH <2 (using 1.2% aqueous phosphoric acid), and then processed on an S-trap micro column according to manufacturer’s recommendations (Protifi S-Trap micro protocol). Proteins were digested with Promega modified trypsin (2 µg for 1:10 ratio, 1 hour 37 °C; Thermo Fisher Scientific). Resulting peptides were evaporated to dryness and resuspended in 2% acetonitrile (ACN), 0.1% formic acid with final concentration of 0.5 µg protein µl−1. Prior to analysis on the mass spectrometer, each sample received 50 fmol µg−1 protein of Pierce Retention Time Calibration (PRTC) as an external standard to monitor MS performance across experiments.

Proteomic mass spectrometry overview
All mass spectrometry analyses were performed on an Orbitrap Lumos Fusion Tribrid instrument with an inline EASY-nLC 1200 (Thermo Fisher Scientific). Reverse phase chromatography was achieved using a manufactured PicoTip fused silica capillary column (75 µm i.d., 40 cm long) packed in-house with C18 beads (Dr. Maisch ReproSil-Pur; C18-Aq, 120 Å, 3 µm) and maintained at 50 °C. The analytical column was fitted with a 4 cm long, 150 µm i.d. in-house packed kasil-frit pre-column (Dr. Maisch ReproSil-Pur; C18-Aq, 120 Å, 3 µm). Peptides were eluted using an acidified (formic acid, 0.1% v/v) water-acetonitrile gradient (2-45% acetonitrile). Each sample was analyzed with data-dependent acquisition (DDA) and data-independent acquisition (DIA) in random order. Skyline was used to evaluate the consistency of the resulting peak areas and retention times of QC standards through all analyses80. The DDA MS data are deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD052873. The DIA and PRM MS data are deposited to Panorama Public81,82: (https://panoramaweb.org/Orcas-Island-Time-Series.url).

DDA analyses
Samples from June 12, 2021 (17:00) to June 18, 2021 (13:00) were analyzed on the mass spectrometer operating in DDA mode. The top 20 most intense ions were selected for MS2 acquisition from precursor ion scans of 400–1200 m/z. Centroid full MS resolution data was collected at 70,000 with AGC target of 1×106 and centroid MS2 data was collected at resolution of 35,000 with an AGC target of 5 × 104. Dynamic exclusion was set to 10 seconds and +2, +3, +4 ions were selected for MS2. Comet 2022.0183,84 was used to search the DDA files; parameters included: reverse concatenated sequence database search, 10 ppm precursor mass tolerance, semi-tryptic specificity with 3 allowed missed cleavages, cysteine modification of 57 Da (resulting from the iodoacetamide) and modifications on methionine of 15.999 Da (oxidation). Percolator (version 3.02)85,86 was used to select peptides with 1% FDR. An open-source Nextflow workflow (https://github.com/mriffle/nf-filter-fasta) was used to reduce the metagenomic FASTA from 3,457,818 proteins into a smaller DDA-refined database representing all possible proteins from the peptides detected in the samples for proteomic interpretation. These peptides were mapped back to the full metagenomic FASTA to compile a database of all proteins with at least 1 peptide identified by the workflow. The resulting DDA-refined database contained 36,429 protein sequences (Supplementary Data 2).

DIA analyses
Samples from June 12, 2021 (17:00) to June 18, 2021 (13:00) were analyzed on the mass spectrometer using DIA mode. For each DIA MS experiment, 4 µg digested peptides with 7 fmol PRTC were analyzed. Centroid full MS resolution data was collected at 120,000 and MS2 data was collected at a resolution of 15,000 with an AGC target of 4×105 across isolation windows of 8 m/z. Prosit87 was used to process the DDA-refined database using an HCD intensity and ion retention time prediction model according to methods from Searle et. al., 202088. Six gas phase fractionations (GPFs) spanning 100 m/z for each analysis were completed with a narrow isolation window (4 m/z) to generate a chromatogram library using Encyclopedia89 v. 1.4.01 following Pino et al., 202090. The 36 individual time series MS experiments were searched against the resulting chromatogram library using Encyclopedia v. 1.12.31. Encyclopedia-exported results were visualized using Skyline-daily version 22.2.0.527, including quality assurance that every peptide contained at least three quantifiable transitions and external standards (PRTC peptides) had a coefficient of variation below 20%. Quantitative values for each peptide were generated using Skyline, which calculated the area under the curve (AUC) for each peptide after total ion current (TIC) normalization using Encyclopedia-inferred peak boundaries.

Taxonomic and Functional Annotations
MetaGOmics34 software (version 0.2.0; https://www.yeastrc.org/metagomics) identifies the lowest common ancestor on each peptide’s taxonomic and functional Gene Ontology (GO) assignment using a provided site-specific metagenomic-generated protein sequence database. Here, we used the DDA-refined database as our background proteome and provided a complete list of all peptides identified at each time point with their respective peptide quantities (DIA-MS). The MetaGOmics reports were used to annotate peptides to their respective GO terms and taxonomic classes.

Framework for selecting peptides
Phase 1A: Data dimensionality reduction - PCA. PC1 was driven by missingness of data across the 36 time points and was not included in further exploratory analysis. PC2 and PC3 were separated by day of collection, identifying a critical time point of change. The top 10% of peptides driving variation across PC2 were selected as candidate peptides (Supplementary Data 4).

Phase 1B: Taxonomic trends in summed peptide abundances. Peptide abundances resulting from each class as annotated by MetaGOmics, were summed at each time point and plotted through time (Supplementary Data 5). A Spearman’s correlation coefficient rho (⍴) was then calculated (Supplementary Table 2; Supplementary Fig. 5) comparing the summed quantitative peptide abundance values for each class and ranked time. All peptides from the bacterial classes with a correlation coefficient ⍴>=|0.5| were selected for further development and testing as biomarkers (i.e., Flavobacteriia, Gammaproteobacteria, Verrucomicrobiae, and Opitutae). Alphaproteobacteria was additionally added to the analysis as it represented the most abundant taxonomic class.

Phase 2: HMM State-shift time series analysis. A two-state hidden Markov model (HMM) was then fit to the quantitative values associated with each candidate peptide to identify those that undergo a non-reverting state shift >24 hours prior to the BSD (R package depmixS491; Supplementary Fig. 7-8). The resulting probability matrix was used to determine the latent state of the observations, and the determined state outcomes were plotted through time. Of the 314 peptides determined to undergo a non-reverting state shift, 87 were selected for the next phase of biomarker development.

Phase 3: Matrix-Matched Calibration Curve (MMCC). To determine if the candidate peptides are quantitative by mass spectrometry, a matrix-matched calibration curve (MMCC) was created35. MMCC analysis provides the user with the ability to determine a peptide’s detectability within a complex matrix and whether the peptide yields a linear response when analyzed with targeted parallel reaction monitoring (PRM) methods in a dilution series. A pooled sample consisting of 2 µg aliquots of digested peptides from 36 samples from our time series served as a representative water microbiome. This microbiome sample was serially diluted with a matched matrix of peptides resulting from the marine bacterium Colwellia psychrerythraea 34H92. The resulting dilution series consisted of six microbiome-to-matrix ratios that spanned the expected dynamic range of targeted peptides: 0%, 0.3%, 3%, 50%, 70%, 100% microbiome.

PRM analyses on MMCC
The targeted peptide list was generated using the Skyline analysis of the DIA data. Peptides were trapped and separated on a near-identical chromatography setup detailed earlier with the exception of the analytical column length (30 cm). Mass spectrometry experiments were completed with 2 µg digested peptides with 150 fmol PRTC. Peptides were eluted using the same gradient as DDA/DIA experiments followed by a 30 minute wash with 80% ACN and 10 minute re-equilibration with 3% ACN. Targeted ions were selected for MS2 acquisition from precursor ion scans of 400–1000 m/z with an isolation window of 1.6 m/z, flow rate of 300 nL/min, and resolution of 30,000. Centroid full MS and MS2 data was collected with an AGC target of 4×105. Data was analyzed with Skyline80 software to identify peptide peaks and chromatographic time boundaries. Skyline-calculated peptide quantification from AUC were exported and analyzed with R following Pino et al.35. A total of 27 peptides were identified to yield linear responses via the MMCC with resulting limit of detection (LOD) >1e6.

Biomarker Peptide Validation
To validate the final 27 candidate peptides as predictive HAB biomarkers, PRM MS analyses were performed on the set of samples collected from May 28, 2021 (01:00) to June 3, 2021 (01:00) using the same target peptide list. MS experiments were completed on 2 µg peptides with 150 fmol PRTC external standard. Skyline (version 23.0.9.187) was used to visualize the peptide peak boundaries, which were visualized and manually adjusted based on the parameters established by the MMCC methods. The resulting AUC report was exported from Skyline and evaluated in R to analyze the temporal trends for each peptide. A two-state HMM was fit to each peptide to select for those that undergo a non-reverting state shift, resulting in the final panel of 12 successful biomarker peptides. A Spearman correlation analysis was performed on the 12 peptides to compare peptide abundances across the discovery and validation datasets and assess for consistency of pre-bloom temporal trends. Taxonomic and functional annotations were assigned to each of the 12 peptides by first mapping each peptide back to its respective parent proteins within the DDA-refined database. Subsequently, a standard protein search was performed against the publicly available NR database using the Basic Local Alignment Search Tool (BLAST; 2/6/2024) (Table 1).

Statistics and Reproducibility
Sample analyses on the MS were randomized to reduce batch effects and quality control (QC) peptide mixtures were analyzed every sixth injection to monitor chromatography and MS sensitivity. Quantitative peptide AUC values were first median normalized and log transformed prior to performing the PCA. A Mann-Whitney test was used to test for significance of the PCA results between the first and second half of the time series. For testing correlation strength throughout the study, a Spearman’s rank-order correlation test was used to calculate Spearman’s rho (ρ) and significance was assessed by calculating a two-tailed p-value. The two-state hidden Markov model (HMM) used throughout the study was fit using the depmixS4 package. Time series data was log transformed prior to model fitting. The model was specified to have a Gaussian distribution and transition probabilities were assumed to be stationary. Model parameters were estimated using an Expectation-Maximization (EM) algorithm, with a random seed set to ensure reproducibility. Data are presented as either raw values or the mean ± SD of at least duplicate measures. No data were excluded from the analyses. All statistical analyses were performed using R version 4.3.2.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an impressive body of work demonstrating alterations in the metaproteome of the microbiome in the >24h preceding a harmful algal bloom in a coastal system of the US – a discovery that has potential for development as a bloom forecasting tool. It was apparent that a lot of thought went into how best to present the various data streams and complicated process for data analysis, which was greatly appreciated. The work is novel, interesting, and an enjoyable read, and I expect will be influential within the HAB community and general phytoplankton ecology.\nIntroduction:\nLine 64: please specify “animal” mortalities for clarification since the paragraph discusses both human and animal impacts.\nLine 92-94: Since here is the first indication of the model HAB organism, add references that describe it as a HAB for readers. Also, considering adding the negative consequences typically associated with this particular HAB.\nResults:\nFigure 1, Supplemental Figure 3, and associated text (Lines 111-116): The earlier bloom is obvious in your figure and no further data are needed to confirm you captured the initiation of the 1st bloom; however, more details in the text or figure(s) would strengthen your argument that there was a 2nd bloom after your data collection ceased on June 18. Is it possible to add a timeseries panel for SoundToxins microscopy counts over the whole time period and into the 2nd bloom, so that you can extend Figure 1 beyond June 18? Or is there complimentary chl data they can provide for the entire sampling period and into the second bloom that will 1) ground truth the satellite imagery you provide under Suppl. Figure 3 while also (2) allowing you to show elevated chl after June 18? Or alternatively, can you provide a quantitative measure (that can be similar between both blooms) that can be used to confirm a 2nd bloom initiated on June 18, e.g., a level of chl or cell concentrations reached/exceeded during both events to classify them as blooms, or rate of growth/division reached? Further explanation of which method was chosen to determine BSD in reference 26 (Brody et al., 2013) will also add clarity. The confirmation of a second bloom, and comparison of its duration and/or intensity, relative to the first bloom, will provide additional information regarding how robust this prediction tool might be, e.g., does there need to be a bloom of a certain magnitude (chl, cell count) to induce a detectable response in the preceding microbiome shift. Bloom and BSD criteria will also explain why you didn’t call the elevated chl levels earlier in May “blooms” which would have added more BSDs and validation periods to your analysis, and will allow for others to repeat your approach with other HAB genera.\nLine 123 – 126: Can you please explain why the validation period, May 28 – June 2, extended beyond the BSD. What is the advantage of this approach? This is inconsistent with the discovery period which does not include any dates past the BSD.\nLine 284-296: Can you hypothesize if these peptides (or the tool in general) would be more relevant to (1) diatom-based blooms (e.g., Pseudo-nitzschia spp), as opposed to dinoflagellate or raphidophyte blooms, (2) blooms of autotrophs or mixotrophs, (3) blooms only of a certain density (and therefore might exclude blooms of a low concentration but high potency, e.g., Dinophysis spp.), or (4) toxin-producing vs non-toxin producing blooms?\nJuliette L. Smith (VIMS)", "reviewer_2": "Paper: Mudge et. al:\nHarmful algal blooms preceded by a predictable and quantifiable shift in the oceanic microbiome\nMudge et. al provide evidence that a harmful algal bloom can be detected by non-reverting peptide changes in marine bacteria. By identifying two blooms within their time series, they were able to not only identify peptides of interest, but reconfirm them within different environmental conditions leading up to a bloom, showing the pertinence of the study. Their metaproteomic pipelines were well designed to provide minimal bias and ensure accuracy in their findings. The specific peptides of interest are also connected to maximizing metabolic functions that would be beneficial to bacteria experiencing a bloom, supporting the hypothesis that proteomic shifts can predict an oncoming bloom.\nThe authors perform a solid and detailed analysis that has potential to advance the field of HAB detection. The proteomics approach is well formulated and we especially appreciate that the PRM data and analysis was made publicly available (and is high quality). Our main questions concern the formulation of the study, specifically that the focus on bacterial proteins is not rounded out by consideration of peptides from the bloom forming diatom nor other environmental drivers such as nutrient concentrations or even the concentration of total bacterial cells. While these peptides seem to be good enough for predicting the bloom, where is the evidence that they are the best option? We describe our major comments in more detail below.\nThis review was produced collaboratively by the assigned reviewer and an early career researcher.\nMajor Comments:\nThe authors could do more to explain the motivation for the study. There are many approaches for monitoring HABs - why is it worth trying metaproteomics over other methods such as in situ imaging systems? Specifically, we have questions about the 24-48 hour time frame that the authors identify. How does this time frame compare to other HAB monitoring/reporting methods? Is 24-48 hours enough time to combat economic/health effects of a HAB in this region? Providing this context would support the method of targeting peptides within this time frame and not earlier before the bloom. We also suggest discussing cost-effectiveness of the proteomics approach versus other common approaches.\nThere is consensus in the literature about what makes a good biomarker. These include simplicity (e.g. a relatively direct causal mechanism/relationship between the process being studied and the biomarker), sensitivity, and measurability. The authors do a very good job of describing the sensitivity and measurability of peptide biomarkers, but we have questions about the mechanistic interpretations. The authors decided to look at the bacterial microbiome but did not discuss peptides/proteins from the blooming phytoplankton Chaetoceros socialis. What was the motivation for this? We wonder whether changes in the bacterial community precede changes in the bloom forming algae (suggesting that bacterial peptides respond first) or rather if the bacterial community is changing due to a change in phytoplankton activity or abundance (suggesting that algal peptides might allow for earlier detection). Of course, this will be related to the mechanisms of bloom initiation in the first place. Is this known for this region and if not, can the proteomics data generated in this study be used to clarify this? That would be great! Being able to describe this causality is an important part of biomarker validation, and we therefore highly recommend that the authors include Chaetoceros socialis peptides and proteins in their analysis, which may require analyzing the larger size fractions.\nRelatedly, it seemed to us that the biomarker peptides were diagnosing a shift in the bacterial community composition, seemingly more than they diagnosed a functional shift within that community. This suggests that metaproteomics is able to resolve microbial community composition and shifts in it (that is useful to point out in itself), but also perhaps that other methods for quantifying bacterial composition such as qPCR, which might be faster and/or less expensive, would work just as well. Can the authors comment on this?\nLastly, we are not sure that the manuscript as such is best served by framing it as a ‘new framework/pipeline’ for biomarker discovery, because the methodology is not sufficiently different from standard biomarker or metaproteomic approaches. We recommend focusing more on the environmental or biological findings (including addressing the questions detailed above, that we think the DDA and DIA proteomics data would be helpful for), and removing some of this phrasing in order to make the manuscript stronger.\nMinor Comments:\nMain Text\nLine 39. We think the abstract is well done and all encompassing.\nLine 80-81. The second half of this sentence is a bit convoluted, could it be rephrased?\nLine 96/Line 353. We are a bit confused on what the actual sample volume was. Line 96 mentions collection was every four hours and in the methods, the flow rate is 5L per minute, does this mean each collection was 1,200L?\nLine 198: We need more introduction to the matched matrix for readers who are not proteomics experts.\nLine 236. Could you please cite an example of cell counts being unreliable to help establish why you dismissed it for chlorophyll concentrations?\nLine 245-247. What is a double-dipping method? We are not familiar and it was not mentioned in the reference. Likewise, the reference specifically mentioned proteomics pipelines, which can be quite distinct from metaproteomics. How widespread are these methods?\nLine 265: We don’t have data in the figures to show us the second bloom. How was the second bloom confirmed?\nLine 298: I would avoid statements of “first-ness”\nLine 296-301. These findings are fascinating!\nLine 762. The github page is a bit difficult to navigate. We recommend having a main README that describes the differences between Metagenomics_Sci_Data and Metaproteomics_Nat_Commun. Then separate README files in each folder detailing the specifics of the contained files. This would make approaching the code and data more manageable.\nLine 784. Do you mean blue and pink arrows? The figure does not appear green.\nLine 794. Figure 3b and 3c should have letters switched so the figure description follows the normal right to left reading path and minimizes confusion.\nSupplementary Figures\nLine 38. We recommend moving Supplementary Figure 3, or other visual representation of the location of the experiment, into the main text for geographic context.\nLine 59. Similarly to Figure 3 in the main text, shaded regions representing loess smoothing should be added to Supplementary Figure 6 to aid in visualizing the proteomic response. This should also aid in supporting the bacterial classes chosen for Figure 3.\nLine 66. Similarly to the previous comment, is there a statistical analysis that can be used for graph visualization that will aid in distinguishing the more significant peptides?\nSupplementary Data\nThe various excel sheets offered as supplementary data are well done, but organizationally a bit confusing. Could you potentially combine the excels and have different sheets that are descriptively named? Currently, some are descriptive and others are “Supplementary_data_[number].” It may also be helpful to have a readme file of sorts to explain the different datasets. A description of the columns, especially including units, would also be useful.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We would like to thank the reviewer for their recognition of the work here and its potential to greatly impact the field. Our responses to individual comments are detailed below in blue.\nIntroduction:\nLine 64: please specify “animal” mortalities for clarification since the paragraph discusses both human and animal impacts.\nRESPONSE: We appreciate the reviewer’s suggestion and have amended the text as requested.\nOriginal L61-64: “that said, there are many non-toxic HAB species that achieve remarkably high cell densities that can block sunlight to deeper waters, alter local water chemistry (e.g., pH, nutrients)4, generate hypoxic or anoxic zones5, and/or obstruct fish gills, all potentially leading to a swift rise in mortalities6.”\nRevised L61-64: “that said, there are many non-toxic HAB species that achieve remarkably high cell densities that can block sunlight to deeper waters, alter local water chemistry (e.g., pH, nutrients)4, generate hypoxic or anoxic zones5, and/or obstruct fish gills, all potentially leading to a swift rise in marine animal mortalities6.”\nLine 92-94: Since here is the first indication of the model HAB organism, add references that describe it as a HAB for readers. Also, considering adding the negative consequences typically associated with this particular HAB.\nRESPONSE: We thank the reviewer for their suggestion and have added the consequences of a Chaetoceros socialis HAB along with a reference to support this assertion.\nOriginal L92-94: “The high-resolution time series presented here captured the biological and chemical characteristics of a dynamic coastal ecosystem during two pre-bloom periods of the HAB-forming diatom Chaetoceros socialis.”\nRevised L98-101: “The high-resolution time series presented here captured the biological and chemical characteristics of a dynamic coastal ecosystem during two pre-bloom periods of Chaetoceros socialis, a cosmopolitan and abundant HAB-forming diatom capable of deteriorating water quality and killing fish through physical irritation29.”\nReference 29: Li, Xiao-Yu, Ren-Cheng Yu, Anthony J. Richardson, Chaojiao Sun, Ruth Eriksen, Fan-Zhou Kong, Zheng-Xi Zhou, Hui-Xia Geng, Qing-Chun Zhang, and Ming-Jiang Zhou. 2023. “Marked Shifts of Harmful Algal Blooms in the Bohai Sea Linked with Combined Impacts of Environmental Changes.” Harmful Algae 121 (January):102370. https://doi.org/10.1016/j.hal.2022.102370.\nResults:\nFigure 1, Supplemental Figure 3, and associated text (Lines 111-116): The earlier bloom is obvious in your figure and no further data are needed to confirm you captured the initiation of the 1st bloom; however, more details in the text or figure(s) would strengthen your argument that there was a 2nd bloom after your data collection ceased on June 18. Is it possible to add a timeseries panel for SoundToxins microscopy counts over the whole time period and into the 2nd bloom, so that you can extend Figure 1 beyond June 18? Or is there complimentary chl data they can provide for the entire sampling period and into the second bloom that will 1) ground truth the satellite imagery you provide under Suppl. Figure 3 while also (2) allowing you to show elevated chl after June 18? Or alternatively, can you provide a quantitative measure (that can be similar between both blooms) that can be used to confirm a 2nd bloom initiated on June 18, e.g., a level of chl or cell concentrations reached/exceeded during both events to classify them as blooms, or rate of growth/division reached?\nRESPONSE: We appreciate the feedback from the reviewer concerning bloom capture. We hope that the following response and additional information clarifies to the readership and reviewers the presence of the 2 blooms mentioned. We have chosen to amend Figure 1 by adding a line beneath the full timeseries chlorophyll a trace that shows when Chaetoceros was detected by SoundToxins, Puget Sounds’ HAB monitoring Agency, via microscopy during the field season. We have also taken the advice of another reviewer to move Supplementary Figure 3 to the main text. This figure shows satellite imagery of the region on June 21, three days after the bloom start date, to confirm that the area was experiencing a HAB.\nBecause the type of high-resolution monitoring we performed on the area is not performed outside of this study, it is not possible to provide a measure of elevated chlorophyll outside of the satellite imagery. That said, we have provided additional satellite images (all of which are also available via the provided MODIS reference), microscopy images from our team, and all of the available SoundToxins observations from the field site for select timepoints in our manuscript that is currently in review with Scientific Reports (now available in BioRxiv: https://doi.org/10.1101/2024.09.27.614549). We ask that the reviewer refer to Supplemental Figure 1 of this pre-print (note that we have included below). Sound Toxins confirmed the second bloom on June 22, 2021 via a phytoplankton analysis by volunteer scientists performing microscopy and noting the relative abundance of the phytoplankton types in the water. This data is available with the Scientific Data publication (and in BioRxiv) but is also available here through a private link for reviewers on figshare (https://figshare.com/s/5ef2b91998217cabaab5, file: Nunn_OrcasIsland_Data_SoundToxins). For ease of reviewing this manuscript, we have included the Supplementary Figure 1 below.\nWe believe that these three lines of evidence, in conjunction with our chlorophyll measurements directly at the sample collection site, support the narrative that multiple Chaetoceros blooms occurred in the area.\nReference: Nunn, Brook L., Emma Timmins-Schiffman, Miranda C. Mudge, Deanna Plubell, Gabriella Chebli, Julia Kubanek, Michael Riffle, et al. 2024. “Microbial Metagenomes Across a Complete Phytoplankton Bloom Cycle: High-Resolution Sampling Every 4 Hours Over 22 Days.” bioRxiv. https://doi.org/10.1101/2024.09.27.614549.\nSupplemental Figure 1 (from BioRxiv manuscript accepted to Sci Data Reports). MODIS satellite imagery and on-site microscopy confirm the presence of two distinct phytoplankton blooms. Representative satellite and microscopic images confirm pre-bloom and bloom start for two phytoplankton blooms, as well as mid-bloom and post-bloom of the first phytoplankton bloom. MODIS satellite images of near-surface chlorophyll a concentration (mg m-3) in waters surrounding East Sound, WA, USA. The MODIS Chlorophyll a product is available from both the Terra and Aqua satellites. The sensor and imagery resolution is 1 km, and the temporal resolution is daily. References: MODISA_L2_OC doi:10.5067/AQUA/MODIS/L2/OC/2018. Microscopy images of whole water collected on-site confirm presence of Chaetoceros spp. throughout field season.\nFurther explanation of which method was chosen to determine BSD in reference 26 (Brody et al., 2013) will also add clarity.\nRESPONSE: We are happy to clarify. The topic of bloom initiation and determining a specific BSD remains difficult due to the unpredictable nature of blooms, thus the importance of the research presented here. Many studies of bloom progression occur at a low resolution, with sampling occurring on the scale of days to months. Because of this, it can be difficult to accurately determine the best way to measure a BSD on the scale of measurements occurring every 4 hours as presented here.\nWe decided to utilize the rate of change method utilized in Brody et al., 2013 due to the unique high-resolution nature of our data and the timing of our field season across the blooms. The BSD was determined as the time point with the highest change in chlorophyll a per change in time preceding the bloom peak. As addressed in the previous comment, due to field site time-constraints we were unable to capture the full bloom cycle of the second bloom using the direct measure of chlorophyll a via fluorescent probe. We concede that the exact BSD may differ slightly from what we state in the paper, however, from all the data and evidence that we do have, we believe that the calculated BSD coincides with the initiation of the bloom. The Threshold and Cumulative Sum methods in Brody et al., 2013 are not applicable to our data since we don’t know the bloom peak for the second bloom. We believe the combination of three lines of evidence (chlorophyll measure, microscopy, and satellite imagery) confirm that multiple blooms occurred in the area.\nOriginal L111: “Calculated bloom start dates (BSD)27 were May 31 at 23:00 and June 18 at 11:00 (Supplementary Fig. 2).”\nRevised L118-120: “Calculated bloom start dates (BSD)27, defined as the time point with the highest positive change in chlorophyll a per change in time preceding the bloom peak, were May 31 at 23:00 and June 18 at 11:00 (Supplementary Fig. 2).”\nThe confirmation of a second bloom, and comparison of its duration and/or intensity, relative to the first bloom, will provide additional information regarding how robust this prediction tool might be, e.g., does there need to be a bloom of a certain magnitude (chl, cell count) to induce a detectable response in the preceding microbiome shift. Bloom and BSD criteria will also explain why you didn’t call the elevated chl levels earlier in May “blooms” which would have added more BSDs and validation periods to your analysis, and will allow for others to repeat your approach with other HAB genera.\nRESPONSE: We acknowledge the reviewer’s suggestion and our primary goal is to make this approach clear and concise to simplify the work for other researchers to replicate and test. In reference to the duration and intensity of the bloom, we hope that the detailed SoundToxins data provided on the figshare portal will be sufficient to contextualize the bloom for readers (https://figshare.com/s/5ef2b91998217cabaab5). The file Nunn_OrcasIsland_Data_SoundToxins provides evidence that Chaetoceros socialis was blooming on June 22, 2021 and was common in the water on June 29, 2021. Additionally, satellite imagery on June 19, 2021 in Supplementary Figure 1 shows the scale of the HAB in the East Sound region of Orcas Island. At this time, we cannot comment on the magnitude of the bloom required to induce a detectable response in the microbiome. This is a recognized limitation of the study that we hope to address in future research at different sites for longer durations. We have added a sentence to the discussion to mention this caveat.\nOriginal L343-345: “Additionally, further study is needed in both oceanic and coastal marine ecosystems to validate if these findings are translatable to blooms of other planktonic types, including cyanobacteria and toxin-producing blooms.”\nRevised L367-: “Further study is needed in both oceanic and coastal marine ecosystems to validate if these findings are translatable to blooms of other planktonic types, including cyanobacteria and toxin-producing blooms. While both blooms reported here were documented to affect the greater East Sound area, characterization of the magnitude of the bloom required to elicit a detectable response in the preceding microbiome shift would be beneficial for determining the robustness of these biomarkers across bloom systems of varying duration and intensity.”\nLine 123 – 126: Can you please explain why the validation period, May 28 – June 2, extended beyond the BSD. What is the advantage of this approach? This is inconsistent with the discovery period which does not include any dates past the BSD.\nRESPONSE: We are happy to clarify our reasoning. During the Discovery phase, our data limitations prevented us from extending the analysis beyond the BSD. However, we recognized the practical value of determining whether the biomarker peptides remained non-reverting after the bloom. By incorporating post-BSD time points in the Validation phase, we confirmed that the biomarkers maintained their non-reverting state. This information is crucial for further application, as the biomarkers operate in a binary manner: even if a scientist conducts their analysis after the BSD but before any visual or environmental signs of the bloom, the biomarker remains valid. Without these additional time points, there would be a risk that the biomarkers could revert, narrowing the window for detecting bloom initiation. This would be impractical for researchers unable to sample as frequently or at the high temporal resolution employed in our study. We have amended the text to include our rationale.\nOriginal L123-126: “June 12 - 18 was selected as the discovery period because the environmental parameters monitored were more stable (Fig.1f). May 28 - June 2 was selected as the validation period to test candidate biomarkers (Fig. 1a).”\nRevised L131-136: “June 12 - 18 was selected as the discovery period because the environmental parameters monitored were more stable (Fig.1f). May 28 - June 2 was selected as the validation period to test candidate biomarkers (Fig. 1a) to confirm that peptides remained non-reverting at least 24 hours post-BSD.”\nLine 284-296: Can you hypothesize if these peptides (or the tool in general) would be more relevant to (1) diatom-based blooms (e.g., Pseudo-nitzschia spp), as opposed to dinoflagellate or raphidophyte blooms, (2) blooms of autotrophs or mixotrophs, (3) blooms only of a certain density (and therefore might exclude blooms of a low concentration but high potency, e.g., Dinophysis spp.), or (4) toxin-producing vs non-toxin producing blooms?\nRESPONSE: We appreciate the reviewer’s interest in the relevance of this tool for predicting blooms outside the context of this study. We believe part of the excitement of the findings presented here is their potential to be applicable in a variety of contexts as marine microbiomes reported in metagenomic literature have a high degree of similarity in taxonomic content. Additionally, there is often functional redundancy across the microbiome with many different bacteria expressing similar proteins that execute the same function. With further research, the biomarkers presented here may have the potential to reveal unique metabolic pathways specific to a pre-bloom response. We do not wish to generalize the potential limitations of the biomarkers to the confines of this study until further research has been conducted to validate the biomarkers on blooms of different types as mentioned here. We have restructured the discussion and added the following text below to address this comment.\nRevised L372-377: “Because marine microbiomes globally share similarities in taxonomic structure15,16,17, and bacteria commonly associated with HAB-forming phytoplankton exhibit a functional redundancy that creates overlap in the metabolic potentials of unique taxa76,28, we hypothesize that the bacterial microbiome could generate biomarkers applicable for blooms of various types (diatom, dinoflagellate, raphidophyte, toxin vs non-toxin producing, etc.).”\n-------------------------------------------------------------------------------------------------------------------------------\nReviewer #2 (Remarks to the Author):\nPaper: Mudge et. al:\nHarmful algal blooms preceded by a predictable and quantifiable shift in the oceanic microbiome\nMudge et. al provide evidence that a harmful algal bloom can be detected by non-reverting peptide changes in marine bacteria. By identifying two blooms within their time series, they were able to not only identify peptides of interest, but reconfirm them within different environmental conditions leading up to a bloom, showing the pertinence of the study. Their metaproteomic pipelines were well designed to provide minimal bias and ensure accuracy in their findings. The specific peptides of interest are also connected to maximizing metabolic functions that would be beneficial to bacteria experiencing a bloom, supporting the hypothesis that proteomic shifts can predict an oncoming bloom.\nThe authors perform a solid and detailed analysis that has potential to advance the field of HAB detection. The proteomics approach is well formulated and we especially appreciate that the PRM data and analysis was made publicly available (and is high quality). Our main questions concern the formulation of the study, specifically that the focus on bacterial proteins is not rounded out by consideration of peptides from the bloom forming diatom nor other environmental drivers such as nutrient concentrations or even the concentration of total bacterial cells. While these peptides seem to be good enough for predicting the bloom, where is the evidence that they are the best option? We describe our major comments in more detail below.\nThis review was produced collaboratively by the assigned reviewer and an early career researcher.\nRESPONSE: We appreciate the reviewer feedback and have provided our response below by breaking down the major comments to address each point. We wish to express that we don’t intend to make the claim that the methodology presented here is the best possible approach to bloom prediction. Rather, we present evidence that this novel way of thinking about bloom prediction has demonstrated usefulness in identifying robust, detectable, and quantifiable biomarkers in the form of peptides that exhibited a statistical shift in underlying state in two distinct blooms contexts. Whether there exists another, potentially better, biomarker within the complex molecular signals of these ecosystems is a question for a separate manuscript. As we revised the manuscript, however, we kept this over-arching point raised by Reviewer 2 in our minds with the goal of addressing it throughout.\nMajor Comments:\nThe authors could do more to explain the motivation for the study. There are many approaches for monitoring HABs - why is it worth trying metaproteomics over other methods such as in situ imaging systems? Specifically, we have questions about the 24-48 hour time frame that the authors identify. How does this time frame compare to other HAB monitoring/reporting methods? Is 24-48 hours enough time to combat economic/health effects of a HAB in this region? Providing this context would support the method of targeting peptides within this time frame and not earlier before the bloom.\nRESPONSE: We appreciate the reviewer’s suggestion to clarify the study motivation. Current monitoring processes for HABs often utilize a chain of communication that proceeds from site-specific observations to reporting of environmental events to the appropriate authorities who can disseminate information. Current monitoring approaches detect the toxin in the water column, then alert authorities (pers. comm with WA State SoundToxins and WA State Board of Health). A change in biomarker that occurs 24-48 hours in advance of bloom development would be enough time for researchers to report their findings to the relevant agency who would in turn initiate the process of warning the public. We have direct reporting from SoundToxins, one of the on the ground phytoplankton monitoring groups in the state of Washington, that the development of a tool that could forecast blooms within this time frame would be of great value. We have added additional text to contextualize the benefit of our biomarkers in comparison to in situ imaging or other monitoring methods.\nRevised L357-362: “Current monitoring techniques (i.e., in situ imaging, toxin analysis) require a bloom to reach a critical mass or produce dangerous toxins before HAB detection is confirmed; in contrast, the peptide biomarkers detailed here undergo a detectable and significant shift prior to bloom development thereby allowing sufficient time to alert critical reporting agencies of the oncoming dangers of bloom development prior to public exposure and safety issues.”\nWe also suggest discussing cost-effectiveness of the proteomics approach versus other common approaches.\nRESPONSE: We appreciate the reviewer addressing the importance of cost in conducting field studies. While we agree that proteomics studies can be expensive, we view this research as a vital first step in identifying effective biomarkers that could predict blooms. We have added sentences to address the rapid developments taking place in the field of proteomics research that are continuing to reduce the cost of MS runs and also emphasize that future directions for this study would ideally include the development of a cost-effective tool utilizing these peptides in a way that may not require the use of mass spectrometry. Future research can narrow down a suite of appropriate non-reverting biomarkers (we propose many rather than a few) and those peptides can then be made into a rapid assay-based system similar to a COVID test (a familiar example for us all).\nRevised L265-267: “Presently, mass spectrometry analysis can pose a challenge for researchers unfamiliar with the field; however, as metaproteomic studies gain traction, such methods are becoming more accessible and affordable41.”\nReference 41: Salvato, Fernanda, Robert L. Hettich, and Manuel Kleiner. 2021. “Five Key Aspects of Metaproteomics as a Tool to Understand Functional Interactions in Host-Associated Microbiomes.” PLoS Pathogens 17 (2): e1009245. https://doi.org/10.1371/journal.ppat.1009245.\nRevised L362-365: “the unique state shift occurring in the biomarker peptide abundances lends itself to the development of an assay in which the biomarkers act as an on/off switch for determining bloom initiation, eliminating the need for high-resolution monitoring and cutting development costs for stakeholders invested in using the tool.”\nThere is consensus in the literature about what makes a good biomarker. These include simplicity (e.g. a relatively direct causal mechanism/relationship between the process being studied and the biomarker), sensitivity, and measurability. The authors do a very good job of describing the sensitivity and measurability of peptide biomarkers, but we have questions about the mechanistic interpretations. The authors decided to look at the bacterial microbiome but did not discuss peptides/proteins from the blooming phytoplankton Chaetoceros socialis. What was the motivation for this? We wonder whether changes in the bacterial community precede changes in the bloom forming algae (suggesting that bacterial peptides respond first) or rather if the bacterial community is changing due to a change in phytoplankton activity or abundance (suggesting that algal peptides might allow for earlier detection).\nRESPONSE: We appreciate the reviewer’s perspective on the biomarker discovery and validation process and are happy to clarify our motivation for the study design as presented.\nImportantly, before starting the work we recognized that there are a range of toxic phytoplankton species, yet, at the class-level, the microbiome in the water column has very similar players present with many highly conserved functional domains. Since our expertise lies in metaproteomics, specifically the analysis of bacterial proteomes and how changes in the metabolic pathways observed reflect how these organisms respond to and interact with their environment, it made the most sense to us to first approach the problem by looking at the surrounding microbiome’s response. In our introduction, we provide reasoning and references to support that within the microbiome, bacteria and phytoplankton interact closely and these interactions play a deterministic role in the outcome of bloom progression. We agree with the reviewer that the bacterial community could be changing for a variety of reasons and that more research is necessary to elucidate the cause. We have current collaborators (who are authors on this manuscript) who will be exploring dissolved metabolites at the site, and another who will focus on the phytoplankton fraction. These studies are ~2 years out still. Here, we posit that because this change is robust and precedes the HAB, the main value of these peptide biomarkers is in their predictive potential. We have revised the following text to clarify.\nRevised L71-83: “Interactions between phytoplankton and bacteria have been shown to play a paramount role in shaping the dynamics of nutrient availability, biogeochemical cycles, and microbial diversity and function in marine ecosystems9,10,11. The individual and interacting factors influencing HAB initiation, including chemical, physical, and biological aspects of an ecosystem, are experienced by both HAB-forming phytoplankton as well as their bacterial neighbors. In the specific context of HABs, tight associations have been observed within the microbiome as phytoplankton advance through bloom cycles, with bacteria playing a deterministic role in the outcome of HAB progression and algal byproducts12,13,14. Additionally, certain consistently abundant taxa have been noted as commonly associated with different phytoplankton species15,16,17, lending precedent to the idea that community structure and function are driven by specific associations18. The influence these organisms have on each other's biological potential creates a recognized opportunity19 to utilize bacterial responses as a rapid in situ predictor for phytoplankton dynamics, regardless of the causative influence that may trigger bloom initiation.”\nFollowing this logic, we hypothesize that as bacteria respond to the environmental or biological triggers that precede bloom development that this change can be detected through measurement of their proteomes. A robust, detectable, and significant change in the proteome, more specifically within individual bacterial peptides, could serve as a predictive biomarker for forecasting future bloom events. The goal of the study is not to characterize bloom development or identify the mechanistic cause of the bloom, but rather to isolate signals of interest (from the bacterial microbiome) that could be generalizable for predicting blooms independent of HAB type, toxicity level, location, environment, etc.\nRevised L372-377: “Because marine microbiomes globally share similarities in structure15,16,17, and bacteria commonly associated with HAB-forming phytoplankton exhibit a functional redundancy that creates overlap in the metabolic potentials of unique taxa73,26, we hypothesize that the bacterial microbiome could generate biomarkers applicable for blooms of various types (diatom, dinoflagellate, raphidophyte, toxin vs non-toxin producing, etc.).”\nOf course, this will be related to the mechanisms of bloom initiation in the first place. Is this known for this region and if not, can the proteomics data generated in this study be used to clarify this?\nResponse: Thank you for the kind remarks, we agree that this information would be of great value! While a deep understanding of the mechanisms behind bloom initiation would be an incredible resource to researchers, this information remains unknown. This is true for the vast majority of HAB-forming phytoplankton, in part due to the myriad of potentially causative factors that vary across HABs.\nWhile we have not isolated the direct mechanisms of bloom development, we did identify multiple biomarkers related to TonB-dependent receptors (Table 1). These proteins are involved in bacterial uptake of scare vital resources, yielding clues to changes in ecosystem that could influence bloom initiation. We also wish to emphasize that this research requires follow-up, ideally in the form of metabolomics studies to integrate how changes in the biology of the microbiome correlate with changes in the environment.\nOne major problem with focusing on isolating the causative mechanisms of one bloom type is that we don’t know what type of bloom will occur in any given area or when. Even within the dataset presented here, we acknowledge that two separate blooms of the same type develop within the same spatial area but completely distinct nutrient and temperature profiles. We show that certain known forcing variables like nutrient concentrations don’t undergo a significant change in abundance suitable for predicting bloom events (Fig. 1). Per the correspondence of SoundToxins (data provided above), other bloom types develop in this area regularly that include Pseudo-nitzschia spp., Alexandrium and Heterosigma akashiwo, all of which are toxin-producing. Because the environmental conditions that researchers generally attempt to correlate with bloom formation vary greatly between the two HABs presented here, it is likely that the phytoplankton dynamics will also vary. However, the bacterial peptides show a consistent response despite these differences.\nWe provide the following reference to address the difficulties in developing a mechanistic model for HAB development:\n“In general, the many biological processes that contribute to HAB development remain poorly defined and present major sources of uncertainty in process-based models.” (Ralston and Moore, 2020)\n“Process-based (or mechanistic) models use mathematical equations to explicitly simulate key physical and biological processes that govern HABs and HAB outcomes. Their development requires detailed knowledge of critical life history characteristics and the factors that modulate them as well as transport pathways. … Because process-based models are more comprehensive than statistical models, they take more time and effort to develop and are more computationally expensive to run.” (Ralston and Moore, 2020)\nOriginal L345-348: “While capturing a bloom system prior to onset remains difficult, isolating this critical time frame in future studies could advance our understanding of what triggers a bloom and how community structure is impacted by environmental factors prompting bloom initiation.”\nRevised L380-383: “While capturing and characterizing a bloom system prior to onset remains difficult76, isolating this critical time frame in future studies could advance our understanding of what triggers a bloom and how community structure is impacted by environmental factors prompting bloom initiation.”\nReference: Ralston, David K., and Stephanie K. Moore. 2020. “Modeling Harmful Algal Blooms in a Changing Climate.” Harmful Algae 91 (January):101729. https://doi.org/10.1016/j.hal.2019.101729.\nBeing able to describe this causality is an important part of biomarker validation, and we therefore highly recommend that the authors include Chaetoceros socialis peptides and proteins in their analysis, which may require analyzing the larger size fractions.\nResponse: In contrast to the bacterial metaproteome, phytoplankton peptides from mixed natural assemblages are more challenging to extract and accurately analyze and quantify. To expand on this point, a benefit of identifying predictive biomarkers using a bacterial microbiome approach is that proteomic/metaproteomic pipelines for eukaryotic organisms are exceedingly more difficult to the point that this would prove to be a barrier of entry for other researchers in the field. While our motivation stems from our own expertise in processing and analyzing thousands of bacterial samples, there are many studies and ongoing efforts that substantiate the difficulty in analyzing mixed communities of phytoplankton. For example, proteomics is reliant on accurate metagenomes to be completed. In the case for phytoplankton, we would need accurate transcriptomics and downstream transcriptomes would need to be time matched to each proteomic samples. Further, phytoplankton produce chlorophyll which makes determining protein concentration challenging due to widespread use of colorimetric methods for quantifying protein. Methods for improving protein quantification would need to be developed. We will be working on this. We hope the revision below adequately addresses this concern.\nRevised L92-96: “While proteomic analysis of eukaryotic organisms like phytoplankton remain difficult26, a predictive biomarker developed using a bacteria-centric approach would potentially eliminate the need for targeted forecasting due to taxonomic and functional overlap of marine microbiomes globally27,28.”\nReference 26: Keeling, Patrick J., Fabien Burki, Heather M. Wilcox, Bassem Allam, Eric E. Allen, Linda A. Amaral-Zettler, E. Virginia Armbrust, et al. 2014. “The Marine Microbial Eukaryote Transcriptome Sequencing Project (MMETSP): Illuminating the Functional Diversity of Eukaryotic Life in the Oceans through Transcriptome Sequencing.” PLOS Biology 12 (6): e1001889. https://doi.org/10.1371/journal.pbio.1001889.\nIn addition to the physical efforts, time, and expertise required to perform the reviewer’s request, a significant monetary investment would be required to both sequence the eukaryotic metagenome (to follow our pipeline for database construction) and analyze the eukaryotic metaproteome via mass spectrometry (which the reviewer recognized in a later comment is a costly analysis). Finally, all of these efforts would result in a major influx of information that is outside the scope of this study. While Chaetoceros socialis peptides may undergo a detectable change prior to HAB formation, there is a chance they would only be applicable to C. socialis blooms. In contrast, the bacterial microbiome surrounding the C. socialis HAB contains many bacterial classes that are commonly found surrounding HABs of other types, increasing the likelihood that biomarkers would be applicable across different HAB types. Shifting the emphasis from the microbiome to the phytoplankton would limit the applicability of the research to forecasting efforts.\nWhile we agree that more validation is necessary before this tool can be confirmed for future use, we believe this manuscript presents strong evidence that the biomarker peptides identified have predictive potential and are worth disseminating to the Nature Communications readership, along with a detailed pipeline for other researchers to replicate in their own studies.\nRelatedly, it seemed to us that the biomarker peptides were diagnosing a shift in the bacterial community composition, seemingly more than they diagnosed a functional shift within that community. This suggests that metaproteomics is able to resolve microbial community composition and shifts in it (that is useful to point out in itself), but also perhaps that other methods for quantifying bacterial composition such as qPCR, which might be faster and/or less expensive, would work just as well. Can the authors comment on this?\nRESPONSE: We thank the reviewer for their comment on alternative means for characterizing the microbiome. While we agree that other methods may be equally effective at isolating shifts in community composition, the goal in focusing on peptides was to identify biomarker potential in both taxonomic and functional contexts. When the study was initiated, we anticipated that we would discover details on how the microbiome was functioning (i.e., what enzymes were increased or decreased so we could potentially determine if there were key metabolites present in the water column). In the discovery bloom, the trends experienced by classes like Verrucomicrobiae hint at potential taxonomic drivers of blooms in that the many of peptides annotated to the class in Phase 1 fit a non-reverting two-state hidden Markov model in Phase 2 (Supplementary Table 2). In contrast, Gammaproteobacteria as a class exhibited a gradual decline in peptide abundance, and the biomarker peptides isolated from the study had both positive and negative abundance trends specific to their protein of origin (Table 1). This is why we argue that the framework established here of confirming biomarkers with the HMM is important, as trends in taxa alone don’t convey the pre-bloom dynamics in the same detail as the individual peptides.\nWe also want to emphasize that we entered this research endeavor with the idea that we could glean a greater mechanistic understanding of the factors influencing bloom initiation by utilizing metaproteomics. Peptides capture both the taxonomic composition and the active functionality of cells as they alter their proteome in response to changes in their environment. In this case, if we are interested in eventually understating key chemical metabolites or signals that trigger or cause blooms, quantitative metaproteomics is superior to metagenomics, metatranscriptomics, etc. for isolating and identifying the time-specific response during complex pre-bloom dynamics (reference provided in text). While the trends revealed in the discovery bloom reveal potential bacterial classes of interest, a suite of peptides make ideal biomarkers for the reasons the reviewer emphasized in a previous comment (measurability, sensitivity, ability to be made into rapid assay, etc.).\nThe proteomic study conducted here identified peptides that exhibit a significant shift in state that is more useful as a specific and precise biomarker for HAB initiation than metagenomic content based on a study conducted in parallel with the work presented here. A metagenomic analysis of relevant gene count for the classes of bacteria represented across the time series (one metagenome per timepoint) was conducted as part of a separate manuscript in prep (note that proteomics were started two years prior to the metagenomes as a PhD thesis). The metagenomic data shows that while certain classes of bacteria do undergo observable increases/decreases in relative gene count, the trends are not consistent prior to both blooms. Additionally, the timing of the increase for classes that undergo a significant change in abundance prior to the discovery bloom (i.e., Verrucomicrobiae) do so gradually and within a timeframe that does not make clear the incoming BSD. The gene counts gradually increase across a six day period, while prior to the validation bloom relative gene counts are consistent prior to the BSD. In this way, peptides provide a precise measure of identifying bloom initiation even outside of obvious changes in metagenomic content. We chose not to include this information here as the details surrounding the preparation of the metagenomic data are substantial and independent of the metagenomic sequencing we used in the manuscript to create the DDA-refined metagenome. We agree that follow-up studies to identify the best way to utilize the exciting results detailed in this data are necessary. Our goal is to release the data to the public as quickly as we can (e.g., BioRxiv has the metagenomes and all JGI data is open source, and this manuscript provides all raw and searched proteomic data for users to interrogate and model). More is on the way!\nLastly, we are not sure that the manuscript as such is best served by framing it as a ‘new framework/pipeline’ for biomarker discovery, because the methodology is not sufficiently different from standard biomarker or metaproteomic approaches. We recommend focusing more on the environmental or biological findings (including addressing the questions detailed above, that we think the DDA and DIA proteomics data would be helpful for) and removing some of this phrasing in order to make the manuscript stronger.\nRESPONSE: We appreciate the reviewer’s interest in the biological findings presented in the manuscript. While we agree that a new mass spectrometry technique was not invented here, we purposefully formatted the manuscript as a framework for biomarker discovery due to the novel nature of the pipeline for evaluating metaproteomic data from a bacterial microbiome. Further, we would argue that no one has used time series statistics to explore how a microbiome’s peptides change through time in order to understand a dynamic ecosystem (in particular in the ocean, but we could also find no references on Google Scholar, nor after speaking to multiple time series statistics experts). Further, as addressed in a response below, we utilized a new approach to curating a protein database for DIA-MS analysis using DDA results that is statistically robust and unbiased. We also utilized a hidden Markov model to identify which peptides undergo a state shift prior to the BSD and the specific timing of the shift. Again, to our knowledge, this is the first research utilizing an HMM in the context of biomarker development.\nAdditionally, there is the possibility that the microbiome associated with different types of HABs would undergo similar state shifts but in peptides with a different taxonomic or functional origin. By creating and detailing the pipeline here, researchers could identify subsets of peptides that have relevance to the blooms occurring in their area. This would make the manuscript relevant to the field even if the biomarker peptides identified here are not tractable to other bloom systems. Last, we did not want to overinterpret these results as we were focused on using accurate quantitative data to discover trends in an unbiased manner and follow up with an unbiased test on the 2nd bloom.\nMinor Comments:\nMain Text\nLine 39. We think the abstract is well done and all encompassing.\nRESPONSE: We thank the reviewer for including this comment.\nLine 80-81. The second half of this sentence is a bit convoluted, could it be rephrased?\nRESPONSE: We are happy to clarify and have made the following change:\nOriginal L80-81: “The influence these organisms have on each other's biological potential creates a recognized opportunity19 to utilize rapid bacterial active responses to individual and interacting chemical, physical, and biological factors as an in situ predictor for phytoplankton dynamics.”\nRevised L81-83: “The influence these organisms have on each other's biological potential creates a recognized opportunity19 to utilize bacterial peptides as a rapid in situ predictor for phytoplankton dynamics, regardless of the chemical, physical, and biological factors that may trigger bloom initiation.”\nLine 96/Line 353. We are a bit confused on what the actual sample volume was. Line 96 mentions collection was every four hours and in the methods, the flow rate is 5L per minute, does this mean each collection was 1,200L?\nRESPONSE: The sample volume was approximately 1.2 L per sample, but always between 500 mL - 2000 mL. The text has been edited to reflect this in the Methods section.\nOriginal L351-353: “Briefly, field sampling was conducted in East Sound, WA, USA from May 27 - June 18, 2021, by pumping water from 100 m offshore at a constant depth of 2 m and flow rate of 5 L min-1…”\nRevised L386-389: “Briefly, field sampling was conducted in East Sound, WA, USA from May 27 - June 18, 2021, by pumping water from 100 m offshore at a constant depth of 2 m and flow rate of 5 L min-1 to obtain 20 L. From the 20L, multiple peristaltic pumps were set up to simultaneously filter 500 – 2000 mL total volume per sample, averaging 1L per microbiome filter.”\nLine 198: We need more introduction to the matched matrix for readers who are not proteomics experts.\nRESPONSE: Thank you for this suggestion. We agree that the MMCC experiment can be confusing if not well explained and have added extra text to explain the concept of a matrix matched calibration curve.\nOriginal L184-196: “Because metaproteomic measurements using MS techniques are prone to confounding matrix effects, we constructed a matrix-matched calibration curve to ensure that candidate peptides could be accurately and reliably quantified prior to further analysis (Fig. 2, blue). The MMCC consisted of a dilution series of this study’s microbiome peptides against a matched-matrix of Colwellia psychrerythraea 34H peptides to determine if candidate biomarkers exhibited a linear response during targeted MS experiments and establish a limit of detection (LOD) for each peptide. The MMCC method followed by parallel reaction monitoring (PRM) allowed evaluation of the LOD for the 87 peptides down to a detectable threshold of 1e3, and the dilution series demonstrated a linear trend above the LOD (Supplementary Fig. 7).”\nRevised L195-208: “Because each microbiome sample carries a complex mixture of thousands of peptides, determining the quantity of any individual peptide abundance within the sample can be difficult. In order to ensure that the selected biomarker peptides are both detectable and quantifiable when present in complex mixtures, we created serial dilutions of the sample to generate a calibration curve prior to further analysis (Fig. 2, blue). Because there are interfering ions simultaneously eluting from the chromatography column prior to MS analyses which can influence detectability of ions and the resulting quantities reported, these serial dilutions were generated with a complex suite of peptides, rather than water. In this dilution series, the total amount of protein injected for each MS run remained the same while the ratio of microbiome sample to diluent was varied incrementally. The diluent used here consisted of a protein digest of Colwellia psychrerythraea 34H, a model bacterium with no shared endogenous peptides. By diluting our test microbiome with peptides of interest against a background of equivalent biochemical complexity with no shared endogenous peptides, we were able to determine if candidate biomarkers exhibited a linear response during targeted MS experiments and establish a limit of detection (LOD) for each peptide. This method, the matrix-matched calibration curve (MMCC), is employed for peptide biomarker development in complex systems36. The MMCC method followed by parallel reaction monitoring (PRM) allowed evaluation of the LOD for the 87 peptides down to a detectable threshold of 1e3, and the dilution series demonstrated a linear trend above the LOD (Supplementary Fig. 6).”\nLine 236. Could you please cite an example of cell counts being unreliable to help establish why you dismissed it for chlorophyll concentrations?\nRESPONSE: We are happy to clarify. We chose to utilize chlorophyll concentrations because cell counts rely on human input, making more prone to human error and more labor intensive. A probe can monitor chlorophyll concentration constantly (every 10 minutes) and results are immediate, while cell counts are limited due to the need for human input and time constraints. Cell counts also can vary depending on the phytoplankton present in the water column, the size of those phytoplankton, and the discretion of the person performing the cell count. Additionally, the use of a fluorescent probe in monitoring chlorophyll concentration as a measure for phytoplankton abundance is widely accepted as a valid proxy. We have included a reference validating the use of a fluorescent probe in monitoring phytoplankton quantity. Last, we were looking for methods that are rapid and easy to employ for a range of sites for the future.\nOriginal L236: “Additionally, probe-based hourly chlorophyll a abundance was used to calculate the BSD27, streamlining the assessment of HAB onset independent of laborious, and sometimes unreliable, methods used to quantify cell counts.”\nRevised L248-250: “Additionally, probe-based hourly chlorophyll a abundance was used to calculate the BSD31, streamlining the assessment of HAB onset independent of laborious, and sometimes unreliable, methods used to quantify cell counts39.”\nReference 39: Gregor, J, and B Maršálek. 2004. “Freshwater Phytoplankton Quantification by Chlorophyll a: A Comparative Study of in Vitro, in Vivo and in Situ Methods.” Water Research 38 (3): 517–22. https://doi.org/10.1016/j.watres.2003.10.033.\nLine 245-247. What is a double-dipping method? We are not familiar and it was not mentioned in the reference. Likewise, the reference specifically mentioned proteomics pipelines, which can be quite distinct from metaproteomics. How widespread are these methods?\nRESPONSE: We are happy to clarify for readers. In this section, we use the term double-dipping in reference to the statistical problem created when one data set is used to both create a hypothesis and also test that hypothesis. We have added more detail to the text as well as a reference that explains double-dipping. We agree that metaproteomics methods are distinct and more complicated than traditional proteomics. Many labs have developed their own method for analyzing data, particularly for environmental samples where a reference database must be independently constructed. We assert in the text that our approach is unique in that it is unbiased and robust. This is also part of the reason we drafted the paper with an emphasis on the framework for biomarker development, as the choice to curate our DIA data using DDA data first is a new approach in the field. An analogy to double dipping: Imagine you're searching for a needle in a haystack, and you decide to take a handful of hay, sift through it, and pull out what looks like a needle. Then, without discarding that same handful, you sift through it again, hoping to find another needle.\nIn statistics, this \"double dipping\" happens when you reuse the same data to test hypotheses after already selecting them. Just like repeatedly searching the same handful of hay increases the chances of finding something that just looks like a needle, reanalyzing the same data increases the likelihood of finding patterns that appear significant but are really just random noise. This can lead to misleading conclusions because you're effectively giving yourself multiple chances to stumble upon something by chance (or stumbling on the thing you were looking for), rather than through meaningful evidence.\nOriginal L245-247: “By doing the database reduction using DDA data that is distinct from the DIA data used in the discovery Phases, this approach avoids potential biases incurred by double-dipping methods widespread among other metaproteomic pipelines29.”\nRevised L259-268: “By doing the database reduction using DDA data that is distinct from the DIA data used in the discovery Phases, this approach avoids potential biases incurred by double-dipping methods that assess the value of signals of interest from the same data used to identify those signals, a widespread occurrence among other metaproteomic pipelines33,40,41. This novel approach allowed for confident detection and quantification of DIA peptide signals in an unbiased, statistically robust manner, while accurately controlling the false discovery rate (FDR).”\nReference 33: Everett, L. J., Bierl, C. & Master, S. R. Unbiased Statistical Analysis for Multi-Stage Proteomic Search Strategies. J. Proteome Res. 9, 700–707 (2010).\nReference 40: Ball, Tali M., Lindsay M. Squeglia, Susan F. Tapert, and Martin P. Paulus. 2020. “Double Dipping in Machine Learning: Problems and Solutions.” Biological Psychiatry. Cognitive Neuroscience and Neuroimaging 5 (3): 261–63. https://doi.org/10.1016/j.bpsc.2019.09.003.\nReference 41: Salvato, Fernanda, Robert L. Hettich, and Manuel Kleiner. 2021. “Five Key Aspects of Metaproteomics as a Tool to Understand Functional Interactions in Host-Associated Microbiomes.” PLoS Pathogens 17 (2): e1009245. https://doi.org/10.1371/journal.ppat.1009245.\nLine 265: We don’t have data in the figures to show us the second bloom. How was the second bloom confirmed?\nRESPONSE: We agree that this is important to clarify. The second bloom was confirmed using two separate lines of evidence. The first confirmation was microscopy from SoundToxins Phytoplankton Monitoring Network through the Washington Sea Grant program: weekly measurements examining the phytoplankton content of the water column at the field site confirmed a bloom of Chaetoceros socialis on June 22, 2021. This data is included in a separate manuscript that has been accepted for publication at Scientific Data as Supplemental Figure 1. Additionally, the SoundToxins data is available at the following private repository via figshare (https://figshare.com/s/5ef2b91998217cabaab5, file: Nunn_OrcasIsland_Data_SoundToxins). The data will be made publically available prior to publication.\nThe second confirmation was satellite imagery on June 19, 2021 showing an increase in chlorophyll a levels at the field site. For ease of review, we have chosen to publish this manuscript to Bioarxiv (https://doi.org/10.1101/2024.09.27.614549). We ask that the reviewer refer to Supplemental Figure 1 which we have also included below. We hope that this data will adequately satisfy the reviewer’s request.\nSupplemental Figure 1. MODIS satellite imagery and on-site microscopy confirm the presence of two distinct phytoplankton blooms. Representative satellite and microscopic images confirm pre-bloom and bloom start for two phytoplankton blooms, as well as mid-bloom and post-bloom of the first phytoplankton bloom. MODIS satellite images of near-surface chlorophyll a concentration (mg m-3) in waters surrounding East Sound, WA, USA. The MODIS Chlorophyll a product is available from both the Terra and Aqua satellites. The sensor and imagery resolution is 1 km, and the temporal resolution is daily. References: MODISA_L2_OC doi:10.5067/AQUA/MODIS/L2/OC/2018. Microscopy images of whole water collected on-site confirm presence of Chaetoceros spp. throughout field season.\nLine 298: I would avoid statements of “first-ness”\nRESPONSE: We believe the reviewer is referring to line 278 and have revised the sentence to remove the word “first” as requested.\nOriginal L278: “This is the first study to statistically identify and isolate this critical time frame in pre-bloom dynamics for future research and funding mechanisms to focus on.”\nRevised L295-296: “Here we statistically identify and isolate the critical time frame in pre-bloom dynamics for future research and funding mechanisms to focus on.”\nLine 296-301. These findings are fascinating!\nRESPONSE: We agree with the reviewer and are excited to share this research with the Nature Communications readership.\nLine 762. The github page is a bit difficult to navigate. We recommend having a main README that describes the differences between Metagenomics_Sci_Data and Metaproteomics_Nat_Commun. Then separate README files in each folder detailing the specifics of the contained files. This would make approaching the code and data more manageable.\nRESPONSE: The GitHub page has been amended to include a README page that clearly navigates the viewer to the files for the paper and distinguishes between folders for the different manuscripts contained within.\nLine 784. Do you mean blue and pink arrows? The figure does not appear green.\nRESPONSE: We apologize for the confusion and thank the reviewer for bringing this to our attention. A previous iteration of the figure contained green points. The colors noted in the text have been corrected to reflect the figure as presented.\nLine 794. Figure 3b and 3c should have letters switched so the figure description follows the normal right to left reading path and minimizes confusion.\nRESPONSE: We agree that this was organizationally confusing and have switched the lettering on Fig. 3B and 3C. We have also updated the text so all references to the figure are accurate.\nSupplementary Figures\nLine 38. We recommend moving Supplementary Figure 3, or other visual representation of the location of the experiment, into the main text for geographic context.\nRESPONSE: We thank the reviewer for this suggestion and agree that Supplementary Figure 3 will add valuable context for readers to understand the location and impact of the blooms studied on the surrounding area. We have moved Fig. S3 to be the new Figure 2 and corrected all figure labeling in the text.\nLine 59. Similarly to Figure 3 in the main text, shaded regions representing loess smoothing should be added to Supplementary Figure 6 to aid in visualizing the proteomic response. This should also aid in supporting the bacterial classes chosen for Figure 3.\nRESPONSE: We agree with the reviewer that this change will provide a valuable visual aid to readers. Loess smoothing has been added to Supplementary Figure 6 for all classes.\nLine 66. Similarly to the previous comment, is there a statistical analysis that can be used for graph visualization that will aid in distinguishing the more significant peptides?\nRESPONSE: While we appreciate the reviewers request, we don’t want to suggest to readers that one peptide might be better for finding and predicting blooms than another. We specifically think that using a suite of biomarkers like a biosignature is more appropriate than putting all your faith on a single biomarker. Sorry, but we hope the reviewers understand.\nSupplementary Data\nThe various excel sheets offered as supplementary data are well done, but organizationally a bit confusing. Could you potentially combine the excels and have different sheets that are descriptively named? Currently, some are descriptive and others are “Supplementary_data_[number].” It may also be helpful to have a readme file of sorts to explain the different datasets. A description of the columns, especially including units, would also be useful.\nRESPONSE: We appreciate the reviewer’s suggestion and have added a key to the Data Availability section with a brief description of each numbered dataset. We have also included details on each file, including column descriptions, within the Supplementary Information file. All datasets have been uploaded to figshare along with a key to help readers navigate each file. Datasets containing spreadsheets have been condensed so that each file is now one tab of a single excel file. Please note that because the files are now uploaded to figshare we made the decision to transition all names from “Supplementary Dataset” to “Dataset”."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors addressed all of my questions and concerns, by either adding rationale or identifying future research directions based on their findings.", "reviewer_2": "We thank the authors for addressing our comments. However we still have some lingering questions/comments. Overall, this is a nice body of work and an impressive dataset. The below comments are made because we see opportunities to increase the long-term relevance of this paper, which has potential to become an important work in HAB research and marine metaproteomics, if it is carefully framed.\nIn our prior review, we asked for clarification of the data/code (e.g. a README). This is especially poignant given that the authors frame their work as a new pipeline. We see that a README was created, but it is not sufficient. A complete README would include a delineation of the files in the folder and column names, with brief descriptions of each, especially since abbreviations are used.\nWe are still unsure about the reasoning for considering only bacterial proteins/peptides and avoiding analysis of the eukaryotic phytoplankton peptides. While we acknowledge that it is not easy, a number of studies have successfully analyzed the metaproteomes of eukaryotic phytoplankton (McCain et al., 2021, Chen et al, 2021, Zhao et al., 2024, Cohen et al., 2021). The authors argue that focusing on the bacteria might allow phytoplankton blooms to be predicted without knowing a-priori which phytoplankton are/will be blooming. However, bacterial populations can also shift and indeed the evidence from this study points to this (e.g. shifts in bacterial community composition were resolved in the metaproteomic data). So, how do the authors know that these bacterial peptides will work for blooms of other eukaryotic phytoplankton in this region? Note: we do not mean to totally dismiss the approach, but rather point out that the phytoplankton themselves could and probably be considered in a complete analysis. This would enhance the validity of the biomarkers by helping to establish their mechanistic basis. Even though this may be difficult, we think that the benefit would warrant this effort.\nWe disagree with the statement that cell counts require human input; there are a variety of means (automated flow cytometry with appropriate controls, imaging based methods) that can provide more objective results. The point being, it is unclear to us whether the metaproteomics based methods are actually able to provide faster information than other HAB monitoring methods (we say this as fans/advocates of metaproteomics in general, who are concerned about/want to avoid overselling of metaproteomics-based results which might later hinder future efforts to develop these methods by the authors and others).\nLastly we still think that this work would be better served by softening some of the language about this being a new pipeline or approach. While this might be the first time this pipeline is used for marine HABs, we can think of more than a few examples where DDA data is used to inform DIA searches (this is the basis of spectral library generation), or where HMM models are used to identify key features of a proteome (e.g. Herod et al., 2020 Nat Comms). Again, this is not to dismiss the work, but rather to encourage further strengthening and focus on the scientific discoveries. If the authors really want to establish a new pipeline, the best way to demonstrate this would be to test the pipeline for a variety of biomarker discovery systems.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ">Thank you for your thoughtful feedback and for recognizing our efforts to address your questions and concerns. We greatly appreciate your insightful comments, which have helped strengthen our manuscript.\n> We apologize that we provided README that is not sufficient. We have now included a complete README file on the GitHub site; it includes delineation of the files in each folder and column names with descriptions, as requested. It is additionally attached for simpler review.\n>Thank you for recognizing the potential significance of this paper and our research. While we would love to expand our analysis to include eukaryotic phytoplankton peptides, this lies beyond the scope of the currently funded work for several reasons. Additionally, we aim to address the reviewer’s request for clarification regarding the differences between predictive and diagnostic biomarkers (see response below). Briefly, predictive biomarkers are used to identify the likelihood of a system responding to a particular treatment or environmental condition, whereas diagnostic biomarkers typically require a mechanistic understanding derived from the findings. In our case, this distinction suggests that developing a detailed biological mechanism is not essential at this time but remains a promising avenue for future exploration.\nWhile we acknowledge the contributions of the authors mentioned above who have conducted eukaryotic-focused metaproteomic analyses, this was not the primary objective of the funded pilot project. We agree that such analyses are possible if an accurate, time- and site-specific protein database is generated. However, the aim of this project was to determine whether the bacterial microbiome exhibited detectable and reproducible changes in its metaproteome prior to a phytoplankton bloom. Our goal was to identify microbiome-specific peptides that could serve as potential predictive biomarkers for bloom initiation using a (biologically and functionally) agnostic approach that relies on 1) accurate mass assignments, 2) accurate peptide sequence identifications, 3) reproducibility and 4) validation on another bloom occurrence.\nMoreover, the Nunn lab also adheres to a rigorous standard for metaproteomics, requiring a time- and site-specific metatranscriptome or metagenome to accurately assign peptides to spectra generated by the mass spectrometer. As such, performing metatranscriptomics would be a necessary prerequisite for further phytoplankton proteomic analyses -without the transctipome we would not be able to assign peptide sequences to the MS data with high confidence. These objectives are planned for future work, but currently, we lack the necessary funding and personnel to undertake this additional effort. Notably, the lead graduate student on this project, Miranda Mudge, has defended her PhD and since moved on. The lab does anticipate publishing additional papers resulting from this work, as we believe it has significant potential. While we are not opposed to completing the eukaryotic fraction, it remains outside the current scope of the project and capacity of our lab.\nThe additional work requested by the reviewer would take years to complete. The need for tools to understand and assess HAB risk is immediate. We have produced results that will move the field of HAB research forward and other researchers and environmental managers should have access to these results on a reasonable timescale. We feel that 1) demonstrating the potential of using the microbiome as a predictor for blooms and 2) that our ability to detect chemical/molecular changes days before bloom initiation should be provided to the science community sooner than later. As these HAB events are occurring more frequently and more intensely, it emphasizes that we, as a scientific community, are running out of time. We hope new investigators are intrigued and tackle this long-standing problem with new perspectives, technology, and ideas. We are excited to be part of that future work.\n>This is an excellent point! We don’t know if these bacterial peptides will work for blooms of other eukaryotic phytoplankton in this region or others. In future work, we expect to explore other blooms to explore the generality of these biomarkers in predicting bloom onset, vs. other changes in the plankton community.\n>To be clear, we are suggesting that we have found predictive biomarkers, not necessarily diagnostic biomarkers, which do require mechanistic understanding. In environmental microbiome research, predictive biomarkers are often identified without a full understanding of the underlying mechanisms. For example, certain microbial taxa have been associated with the health of coral reefs, serving as indicators of reef condition (Baquiran et al., 2024). However, the specific roles these microbes play in reef health are not entirely clear, making them predictive biomarkers without a defined mechanism. Similarly, in soil ecosystems, the presence of particular bacterial groups has been linked to soil fertility and plant productivity, yet the exact processes through which these microbes influence soil health remain unidentified. These instances illustrate how microbial indicators can predict environmental outcomes even when their functional roles are not fully understood.\nAdditionally, a biomarker paper by Califf 2018 (> 1000 citations) provided clarity on the definition: “The basic definition of a biomarker is deceptively simple: 'A defined characteristic that is measured as an indicator of normal biological processes, pathogenic processes or responses to an exposure or intervention (Cagney et al 2018).'” Additionally, within the Calif manuscript, they have a figure that demonstrates the steps in the evaluation framework for biomarkers that we based many of our developments on.\nBaquiran, Jake Ivan P., John Bennedick Quijano, Madeleine JH van Oppen, Patrick C. Cabaitan, Peter L. Harrison, and Cecilia Conaco. \"Microbiome stability is linked to coral thermotolerance.\" bioRxiv (2024): 2024-10.\nCagney, D.N., Sul, J., Huang, R.Y., Ligon, K.L., Wen, P.Y. and Alexander, B.M., 2018. The FDA NIH Biomarkers, EndpointS, and other Tools (BEST) resource in neuro-oncology. Neuro-oncology, 20(9), pp.1162-1172.\nCaliff, R.M., 2018. Biomarker definitions and their applications. Experimental biology and medicine, 243(3), pp.213-221.\nWe have revised the original text to avoid overselling the applicability of these biomarkers and hope that this satisfies Reviewer #2. This includes removing any instances of the word “diagnostic” from the text.\nOriginal L51: “This study proposes and verifies that the changing molecular signals within the ocean microbiome can be used as diagnostic biomarkers for predicting HAB initiation, potentially mitigating detrimental algal bloom outcomes in the future.”\nRevised L51: “This study proposes and verifies that the changing molecular signals within the ocean microbiome can be used as predictive biomarkers for forecasting HAB initiation, potentially mitigating detrimental algal bloom outcomes in the future.”\nOriginal L338-342: “While a diverse screen is advantageous for capturing the heterogeneity of potential HAB systems, having multiple biomarkers that respond differently to the same underlying process adds valuable redundancy to a detection system. If one biomarker fails or provides ambiguous results, others may provide diagnostic information, enhancing the overall reliability of the detection system59.”\nRevised Sentence: “While a diverse screen is advantageous for capturing the heterogeneity of potential HAB systems, having multiple biomarkers that respond differently adds valuable redundancy, thus enhancing the overall reliability of the detection system59.”\nWe have also added the following sentence to emphasize the need for further studies to identify the molecular mechanisms linking our biomarkers to an environmental or biological cause.\nNew Addition L290: “We acknowledge that despite finding predictive biomarkers, we have not revealed a mechanistic understanding or causal factors that lead to these detections. Further research is needed to identify mechanisms that lead to these detections.”\n>We thank you for your concern, but we have taken care to ensure that we are not overstating or overselling what was discovered.\nCurrent monitoring methods are not able to identify biological changes occurring in the microbiome. They do not target the bacterial fraction and sample on a coarse timescale that likely misses key ecosystem inflection points. What was revealed through this investigation does demonstrate the power of proteomic technology detections (i.e, we used 1-2 L of water per sample) and flexibility in platforms for downstream validation (i.e., from discovery to quantitative to targeted MS approaches). We have attempted to soften the language in a few of our statements throughout the manuscript (see above revisions).\nOne of the important, yet overlooked, components of the work was defining the time frame prior to a bloom on which future work should be focused for collecting and analyzing when looking for new forecasting methods. Cell counting, satellite imaging, or detecting the toxin work well, but require the bloom to already be on the rise. We note that even in the point made by the reviewer, HAB “monitoring” methods are the standard to beat, while this research has the potential to move the standard to “predicting” methods, which is a tremendous leap forward in the field. Prior to this work, no research has revealed that there are detectable and quantifiable changes in anything in the water days before the bloom is initiated. This finding fundamentally changes the playing field. We see that there are detectable changes in the community and likely in the water (regardless of exploring the phytoplankton physiology). Further, to our knowledge, this is the only study that has discovered predictive (not diagnostic) biomarkers in the environment and tested them on another bloom to see if they are valid and predictive. We had no idea that would work.\nMillions of dollars from multiple organizations around the world have been devoted to moving the field beyond the current monitoring techniques. We believe this is a solid first step towards creating an assay that could be developed commercially, but because we don’t make claims that this is the end all be all of HAB forecasting, we feel that the conclusions made in this manuscript are valid. Additionally, since we have accomplished this for a single bloom system and made all our methods and results open with this publication, others will be able to relatively easily follow the same pipeline for easy comparison of results.\n>The pipeline - as a specific series of steps conducted in order - was developed in-house. The methods detailed in the manuscript build upon previous developments in our lab by combining them with other approaches. If we were to remove the word “developed” and replace it with “utilized” we would be required to include a citation, for which there is none. Additionally, each of the metaproteomic papers cited above by the reviewer contains its own detailed methodology for collecting, processing, and analyzing the metaproteomic data. We were not trying to oversell the biology as more field sites need to be collected, yet the framework we used for making downstream decisions on this complex time series was generated by our author group as there are no other high resolution metaproteomic time series papers that reveal state shifts across days on a 4 hour timescale. If we were to follow a published pipeline, we would be eliminating steps in our process, which would change the final list of identified biomarkers significantly.\nWe were unable to find the Herod paper.\nPer your suggestion, we have softened the language in the manuscript:\nOriginal Sentence L 137: A new metaproteomic pipeline was developed to 1) confidently identify peptides from this complex microbiome, 2) generate quantitative peptide data, and 3) annotate peptides using a lowest common ancestor approach.\nRevised Sentence L 137: A metaproteomic pipeline was developed to 1) confidently identify peptides from this complex microbiome, 2) generate quantitative peptide data, and 3) annotate peptides using a lowest common ancestor approach."}, {"round": 3, "reviewer_comments": {"reviewer_2": "We thank the authors for addressing our additional comments. We reviewed the new README and find it to be much more informative - thank you for making that change. One more minor point is that we wondered why in Fig 4A/Supp Fig 4 the authors do not note the percent variance explained by the PCs (it could be informative to note this). Thank you also for addressing the difference between predictive and diagnostic biomarkers - we think this helps to clarify the overall goal of the work.\nOn a broader note, we acknowledge that our comments about scope (e.g. including phototrophs) and framing (e.g. as a specific new method) were frustrating. We are still concerned about framing this as a new pipeline as we do not think this is necessary, true, nor strengthens the scientific message. We think the scientific results are interesting enough themselves. Avoiding this language would help the authors to avoid the ongoing conversation we have raised throughout this review process and rather let this interesting science stand for itself. We therefore still highly recommend making this change and (and note it is not necessary to say we 'utilized' a framework and make a citation - just tell us what you did).", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We’re pleased that our edits have made the manuscript more informative. We also thank the reviewer for their comments and suggestions. Per the suggestion on percent variation, we agree these values are informative and have ensured they are included in the text as detailed below.\nL170-174: “A PCA of all quantified peptides from the DIA analysis of the discovery period revealed that variation across PC1 (proportion of variance = 12.67%) was driven by absence of peptide signals across the time series (Supplementary Fig. 4). Visualization of PC2 (proportion of variance = 9.17%) vs PC3 (proportion of variance = 6.67%) revealed two clusters of data separating along PC2 (Fig. 4a).”\nWe appreciate the reviewer’s help in homing in on the main scientific results of the manuscript. One large change that should help to remove emphasis from the framework aspect of the manuscript is an amendment to the abstract – we have deleted text referring to the framework from this section. We also deleted a sentence emphasizing the pipeline in the Discussion section. With this, reference to our development of a pipeline is reduced to one occurrence within the Results section, Figure 3 (a sequential flow of the study design), and the Methods section. While extensive details remain within the methods section, we believe reporting this information is important for encouraging reproducibility. We hope that these changes are enough to soften the focus on the pipeline.\nDeleted L59: “We developed an agnostic, comprehensive, and statistically robust framework to query the microbial metaproteome…”\nDeleted L279-281: “To further optimize data analyses, we developed a statistically robust metaproteomic pipeline to increase the accuracy of peptide identifications through an agnostic scaling of the MS search space.”\nWe once again thank the reviewer for the time and guidance that went into reviewing our manuscript."}]

Decision letter