Marine particle microbiomes during a spring diatom bloom contain active sulfate-reducing bacteria

Abstract Phytoplankton blooms fuel marine food webs with labile dissolved carbon and also lead to the formation of particulate organic matter composed of living and dead algal cells. These particles contribute to carbon sequestration and are sites of intense algal-bacterial interactions, providing diverse niches for microbes to thrive. We analyzed 16S and 18S ribosomal RNA gene amplicon sequences obtained from 51 time points and metaproteomes from 3 time points during a spring phytoplankton bloom in a shallow location (6-10 m depth) in the North Sea. Particulate fractions larger than 10 µm diameter were collected at near daily intervals between early March and late May in 2018. Network analysis identified two major modules representing bacteria co-occurring with diatoms and with dinoflagellates, respectively. The diatom network module included known sulfate-reducing Desulfobacterota as well as potentially sulfur-oxidizing Ectothiorhodospiraceae. Metaproteome analyses confirmed presence of key enzymes involved in dissimilatory sulfate reduction, a process known to occur in sinking particles at greater depths and in sediments. Our results indicate the presence of sufficiently anoxic niches in the particle fraction of an active phytoplankton bloom to sustain sulfate reduction, and an important role of benthic-pelagic coupling for microbiomes in shallow environments. Our findings may have implications for the understanding of algal-bacterial interactions and carbon export during blooms in shallow-water coastal areas.


Introduction
Microalgae inject gigatons of organic carbon into coastal oceans e v ery year (Field et al. 1998 ) and phytoplankton blooms r epr esent primary productivity hotspots.It has been estimated that over 90% of algal-produced carbon is consumed by heter otr ophic bacteria in the immediate vicinity of algal cells, the phycosphere, during typical bloom situations (Seymour et al. 2017 ).Some of these bacteria ar e dir ectl y associated with living algal cells, or with sinking a ggr egates of senescent or dead algae, and ther efor e play an important role in the biological carbon pump.Despite their importance in carbon sequestration (e.g.Bligh et al. 2022 ), vertical connectivity (Mestre et al. 2018 ), and their documented complexity (e.g.Reintjes et al. 2023 ), particle-associated (PA) bacterial comm unities ar e less well understood than their fr ee-living counterparts.In fact, they are often overlooked due to the sometimes fr a gile nature of particles, in combination with the practice of prefiltration to exclude larger organisms prior to molecular analyses (e.g.Simon et al. 2002, Thiele et al. 2015, Heins et al. 2021 ).
Aggregates composed of living algae are known for a div erse micr obiome of aer obic, heter otr ophic bacteria that degr ade complex algal organic matter, such as polysaccharide-rich exudates (Enke et al. 2019, Reintjes et al. 2023 ), although anaerobic metabolism such as diazotrophy is also known to occur (Riemann et al. 2022 ).Within this micr obiome, bacteria ar e selected by factors such as host physiology and genotype (Ahern et al. 2021 ), the surrounding environment (Barreto Filho et al. 2021 ) and via stoc hastic pr ocesses (Stoc k et al. 2022 ).It is methodologicall y c hallenging to study associations between algal and bacterial taxa during natural phytoplankton blooms by direct observations, due to the transient nature of these associations (Seymour et al. 2017, Heins et al. 2021 ), as well as innate complexities and r a pid dynamics of algal and bacterial communities during bloom e v ents (Teeling et al. 2016 ).High-r esolution tempor al co-occurr ence analysis, in combination with measurement of microbial functional potential, offers an indirect way to infer algae-bacteria associations, which can facilitate generating hypotheses about specific interactions and their potential functional implications.
A defining feature of marine particles are the steep chemical and redox gradients that PA bacterial communities are exposed to, compared to free-living, planktonic bacteria (Ploug et al. 1997 ).These gradients have been studied in the context of bathypelagic sinking particles (i.e.marine snow), which harbor micronic hes enabling anaer obic metabolism, suc h as micr obial sulfate reduction (Shanks and Reeder 1993, Bryukhanov et al. 2011, Bianc hi et al. 2018 ).Dir ect micr oelectr ode measur ement of oxygen concentrations in marine particles formed in roller tanks have demonstrated the importance of particle size and sinking velocity (Ploug et al. 1997 ), surrounding oxygen concentration (Ploug and Bergkvist 2015 ), but also the species composition of diatom detritus making up the particles (Zetsche et al. 2020 ) for the formation of anaer obic nic hes.Further, ele v ated concentr ations of sulfide, inside artificial marine snow and field-collected particles compar ed to surr ounding water masses, indicates that sulfate reduction takes place within suc h nic hes (Shanks and Reeder 1993 ).Sulfate-reducing bacteria have also been detected in oxygenated surface waters in the Black Sea, complementing observations of sulfate reduction above 30 m depth in these waters (Bryukhanov et al. 2011 ).Using a modeling a ppr oac h, Bianc hi and collea gues predicted that anaerobic particle microenvironments enabling for example sulfate reduction may be more widespread in the global ocean than pr e viousl y assumed (Bianchi et al. 2018 ).In the photic zone, anaer obic micr o-nic hes hav e been less widel y inv estigated and the pr e v alence of sulfate-r educing bacteria is uncertain in the proximity of oxygen-producing living algal cells.
Here , we in vestigated microbial community dynamics of PA bacterial and eukaryotic taxa during a spring phytoplankton bloom in the southern North Sea at the shallo w-w ater long-term ecological r esearc h site Helgoland Roads (Wiltshir e et al. 2010 ) in the year 2018.We aimed to identify PA bacterial taxa co-occurring with the major eukaryotic taxa (diatoms and dinoflagellates) during the bloom.We hypothesized that bacteria co-occurring with diatoms would be compositionally and functionally distinct from those co-occurring with dinoflagellates.Using 16S and 18S rRNA gene amplicon data from a well-resolved time series (near daily sampling) collected between early March and late May in 2018, we constructed co-occurrence networks focusing exclusiv el y on bacteria-eukaryote co-occurrences in the particle fraction (larger than 10 μm).In addition, we addressed bacterial functional gene expression by analysis of metaproteomes from three selected time points during the bloom.

Sampling and sample processing
A large volume of seawater (40 L-140 L) was sampled using a clean bucket from 1 m depth below the water surface as pr e viously described (Teeling et al. 2012, Wang et al. 2024 ) from the r esearc h v essel Aade in the morning at near dail y interv als between beginning of March and end of May in 2018 at the long term ecological r esearc h (LTER) site Helgoland Roads (50 • 11.3 N, 7 • 54.0 E; DEIMS.iD:https:// deims.org/1e96ef9b-0915-4661-849f-b3a72f5aa9b1 ).The site is located near the small island of Helgoland in the south-eastern North Sea and has a water depth of 6-10 m depending on tide (Wiltshire et al. 2010 ).For chlorophyll a (chl a ) analysis, sample filtration was carried out in a laboratory under dim light to avoid the loss of pigments during the filtr ation pr ocess.We used a combined method of Za pata et al. ( 2000 ) and Garrido et al. ( 2003 ) for chl a extraction and analysis.Pigments wer e separ ated via high-performance liquid c hr omatogr a phy (HPLC) (Waters 2695 Separation Module), and detected with a Waters 996 Photodiode Array Detector.Secchi depth was measured from the vessel on site .T he abundance of "detritus" (non-identifiable matter) was estimated micr oscopicall y on a scale from 0-6, corresponding to "none" (0), "moderate" (3) and "massive" (6) levels .Publicly a vailable wind data from the weather station on Helgoland were obtained via the website www.wetterk ontor.de .Water le v el data, collected at the harbor on the island, were obtained from http:// www.portal-tideelbe.de/ .For 16S and 18S rRNA gene amplicon sequencing and metagenome anal ysis, plankton biomass fr om a 1 L seawater subsample was filtered using 10 μm pore size polycarbonate membrane filters (47 mm diameter, Millipor e, Sc hwalbac h, German y) to separ ate PA microbes ( > 10 μm) from smaller size fractions (not analyzed in this study).At three selected time points during the bloom (J ulian da ys 107, 128 and 144, r epr esenting earl y-, mid-and late bloom phases), a separate filtration was performed using larger (142 mm diameter, Millipore) 10 μm pore size polycarbonate membrane filters for metaproteomic analysis.In order to maximize biomass harvest while avoiding clogging, filtered volumes varied between 15 and 30.5 L per filter for meta pr oteomic anal ysis .T he > 10 μm particle fraction comprises e v erything lar ger than 10 μm, and can include living phytoplankton cells, phytoplankton a ggr egates, zooplankton, fecal pellets and resuspended material of benthic origin.For the purpose of this study, we refer to microbes as particle-associated (PA) if detected in this filter fraction without making assumptions about the nature of the particle.

rRNA gene amplicon sequencing and analysis
Samples from 52 time points throughout the bloom were collected and analyzed.DN A w as extracted from the filters using the Qiagen DNeasy Po w er soil Pro kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions.Dislocation of microbial cells from the filters and mechanical lysis were ac hie v ed by bead beating in a FastPrep 24 5 G (MP Biomedicals , Irvine , C A, USA).DN A concentrations w ere measured at a Qubit 3.0 fluorometer (In vitrogen, Carlsbad, C A, USA).Extracted DN A w as amplified with primer pairs targeting the V4 region of the 16S rRNA gene [515f: 5 -GTGYCAGCMGCCGCGGTAA-3 , 806r: 5 -GGA CTA CNV GGGTWTCTAAT-3 (Walters et al. 2016 )] and the V7 region of the 18S rRNA gene [F-1183mod: 5 -AATTTGA CTCAA CRCGGG-3 , R-1443mod: 5 -GRGC ATC AC AGACCTG-3 ] (Ra y et al. 2016 ) coupled to custom adaptor-barcode constructs.PCR amplification and Illumina MiSeq (Illumina, San Diego, CA, USA) libr ary pr epar ation and sequencing (V3 chemistry) were carried out by LGC Genomics (LGC Genomics, Berlin, Germany).
Sequence r eads fr ee of ada ptor and primer sequence r emains wer e pr ocessed using the D AD A2 pac ka ge (v1.2.0) in R (Callahan et al. 2016 ).In summary, forw ar d and r e v erse Illumina MiSeq r eads wer e truncated to 200 bp, filtered (maxEE = 2, truncQ = 2, minLen = 175), dereplicated and error rates were estimated using the maximum possible error estimate from the data as initial guess.Sample sequences were inferred, paired forw ar d and rev erse r eads wer e mer ged and c himeric sequences wer e r emov ed using the r emov eBimer aDenovo function.The resulting amplicon sequence variants (ASVs) were taxonomically classified using the Silva database (nr 99 v 138.1, Pruesse et al. 2007 ) for 16S rRNA and the PR2 database (version 4.13, minboot: 50, Guillou et al. 2013 ) for 18S rRNA sequences using the build-in RDP classifier.16S rRNA gene amplicon reads classified as c hlor oplasts and mitochondria, as well as the 18S rRNA gene reads classified as Metazoa (zooplankton) were removed prior to downstream analyses (a single time point was excluded due to suspected contamination).The diatom genera Nitzschia , Navicula and Cocconeis were classified as primarily benthic based on local reference literature (Hustedt 1959 ).Corr elation anal ysis was carried out on the r elativ e abundance of these benthic diatom genera and Delsulfobacterota against water le v el, wind speed, and Secc hi depth data using linear r egr ession (function lm).
Co-occurrence networks wer e gener ated in R using Spearman r ank corr elation, as described pr e viousl y (Bengtsson et al. 2017 ).Briefly, we excluded r ar e ASVs with a total abundance < 100 reads (16S) and < 500 reads (18S) across the whole dataset.Then, pairwise correlations between all remaining ASVs were calculated using the rcorr function (Hmisc R pac ka ge), follo w ed b y Pvalue adjustment for multiple testing (function p.adjust) using the Benjamini-Hoc hber g method (Benjamini and Hoc hber g 1995 ).For the final netw ork, w e consider ed exclusiv el y corr elations between 18S and 16S ASVs, and with a correlation coefficient > 0.7 and an adjusted p < 0.01.The netw ork w as then plotted using the igr a ph R pac ka ge (Csardi andNepusz 2006 , R Cor e Team 2023 ).

Metaproteomics
Pr oteins wer e extr acted fr om biomass filters of the 17th of April, 8th of May and 24th of May 2018 (Julian days 107, 128, and 144) as described pr e viousl y (Sc hultz et al. 2020 , Sc hultz 2022 ) and analyzed with liquid chromatography-tandem mass spectrometry in triplicates.Briefly, pr oteins wer e extr acted via bead-beating follo w ed b y acetone pr ecipitation.The extr acts wer e separ ated and fractionated by 1D SDS-PAGE and in-gel trypsin digested.After desalting and concentration of the peptides using C18 Millipore ® ZipTip columns, the samples were measured with an Orbitrap Velos TM mass spectrometer (ThermoFisher Scientific, Waltham, MA, USA).After conversion into mgf file format using MS conv ert in Pr oteoWizard (P alo Alto, CA, USA), spectr a wer e matc hed against a metagenome-based database containing 14764755 entries .Mascot (Matrix Science , London, UK) and Scaffold (Proteome Software , P ortland, OR, US) were used for peptide-spectrummatc hing, pr otein identification and pr otein gr ouping.Instead of setting an FDR-threshold, identification of protein groups was based on number of peptide-matches with a minimum of two, a pr otein thr eshold of 99% and a peptide threshold of 95%.Identified pr otein gr oups wer e annotated via Pr ophane v6.2.3 (Sc hiebenhoefer et al. 2020 ), using the Uniprot-TrEMBL (as of September 2021) and NCBI nr (as of February 2022) databases for taxonomic and EggNOG v5.0.2 (Huerta-Cepas et al. 2019 ) and TIGRFAMs 15 (Haft et al. 2001 ) databases for functional annotation with Prophane default settings.Taxonomic information for Desulfobacterota proteins was manually confirmed against the most recent NCBI nr database (blastp, https://blast.ncbi.nlm.nih.gov, September 2023).TrEMBL taxonomic information for Desulfobacterales was manually curated and set to phylum level for better comparability with the Silv a database.Relativ e abundances wer e calculated as Normalized Spectral Abundance Factor (NSAF) values using the quantification method "max_nsaf" integrated in Pr ophane.Briefly, SAF v alues were calculated by dividing exclusive unique spectrum counts for each protein group by protein length of the longest sequence in that pr otein gr oup.For normalization, SAF v alues wer e then divided by the sum of all SAF values of the sample.

Da ta visualiza tion
Stacked bar plots for rRNA gene amplicon data and metaproteomics were created with R version 4.3.0 using the tidyverse pac ka ge (Wic kham et al. 2019 ) in combination with the svglite, pol yc hr ome, patc hwork, glue and ggnested pac ka ges.

Particle microbial community dynamics during the course of the bloom
The spring phytoplankton bloom in 2018 was c har acterized by an initial dominance of diatoms ( Bacillarioph yceae ), follo w ed b y an increase in dinoflagellate ( Dinophyceae ) relative abundances after the bloom peak (peak in chl a concentration, Fig. 1 A and B).The algal bloom (chl a measurements) peaked around the 26th of April (J ulian da y 116), coinciding with a decrease in water clarity, as indicated by Secchi depth rising at the onset of the bloom.16S rRNA gene amplicon sequencing r e v ealed a total of 17 615 ASVs in the particle bacterial micr obiomes, whic h consisted mainl y of Proteobacteria ( e.g.Spongibacteriaceae, Methylophagaceae, Rhodobacteriaceae) , Bacteroidetes (e.g.Polaribacter ), Verrucomicrobia (e.g.Persicirhabdus ) and Planctomycetes (e.g.Phycisphaeraceae ).Further, Actinobacteria (e.g.Illumatobacter ) , and, notably, Desulfobacterota were also r elativ el y abundant (Fig. 1 C).

Taxonomic and functional composition as detected by metaproteome analysis
Meta pr oteome sampling was performed on three time points that correspond to the early, mid and late phase of the bloom (Julian days 107, 128 and 144), as indicated in Fig. 1 A).We detected 4 584 pr otein gr oups in > 10 μm filter fr actions fr om these sampled time points.Meta pr oteomic anal yses confirmed high abundances of Proteobacteria and Bacteroidetes in the bacterial fraction (Fig. 2 A), and dominance of Bacillariopyhta in the eukaryote fraction (Fig. 2 B), with eukaryote proteins making up 90% of all proteins at the first selected time point.The proportion of eukaryotic proteins was reduced to 70% in the last selected time point (144), while bacterial proteins became comparably more abundant (25% on day 144 compared to 4% on day 107).Functional analysis revealed a predominance of eukaryotic pr oteins involv ed in metabolism, including energy production and conversion, and a shift to w ar ds expression of proteins relevant in cellular processes and signaling over the course of the bloom (Fig. 2 C).

Co-occurrence network analysis
Co-occurr ence anal ysis r esulted in distinct network modules center ed ar ound diatom-and dinoflagellate 18S rRNA gene amplicon sequence variants (ASVs).Based on significant positive correlations (Spearman Rho > 0.7, corrected p > 0.01) between eukary otic 18S rRN A gene and prokary otic 16S rRN A gene ASVs, including the 51 analyzed time points during the 2018 spring bloom, the netw ork w as dominated b y three major distinct modules (Fig. 3 A).Two of these modules were dominated by eukaryotic ASVs belonging to diatoms and dinofla gellates, r espectiv el y, while the third module mostly contained diatom ASVs and was linked to the dinoflagellate-dominated network module.Along the time line of the bloom, these modules r oughl y corr esponded to the phytoplankton taxa pr e v alent in the early stages of the bloom (module I, mainly diatoms), during the late stages of the bloom (mod-ule II, mainly dinoflagellates) and mid-bloom around peak chl a (module III, Fig. 3 A, Fig. 1 A).The composition of bacterial ASVs that co-occurred with diatoms and dinoflagellates is depicted in Fig. 3

Analysis of Desulfobacterota proteins
Out of the 4 584 protein groups detected by metaproteome analysis, 19 were classified as belonging to different Desulfobacterota orders ( Table S1 ).The r elativ e abundance of predicted Desulfobacterota pr oteins av er a ged ov er all time points in the meta pr oteome was 0.14% (expressed as normalized spectral abundance factor-NSAF).This order of magnitude corresponded well to the relative abundance of Desulfobacterota ASVs in the microbiome at the same time points (av er a ge 0.38%, Fig 1 C).Remarkably, no less than 32% of the Desulfobacterota meta pr oteome fr action (in terms of protein abundance) consisted of k e y enzymes for dissimilatory sulfate reduction (Fig. 3 C).Of these , ATP sulfurylase , both the alpha and the beta subunits of dissimilatory sulfite reductase (DSR) and aden yl yl-sulfate (adenosine-5 -phosphosulfate) reductase (APS r eductase) wer e detected.We used a conserv ativ e threshold for protein identification of at least two matching peptides to consider a protein as validly detected.

Assessment of sediment influence
In order to assess the potential influence of resuspension of anoxic sediments on our results, we analyzed local wind speed maxima, tide le v els and Secc hi depth (Fig. 1 A), as well as detritus levels, and the abundance of benthic diatom taxa (Fig. 4 ) along the course of the bloom.The r elativ e abundance of total Desulfobacterota ASVs correlated significantly with Secchi depth (R 2 = 0.33, P < 0.01) and with wind speed (R 2 = 0.11, P < 0.01), but not with water le v el (R 2 = 0.05, p > 0.05), as illustr ated in Fig. 5 .In addition, low Secchi depth coincided with high estimated le v els of "detritus", i.e. micr oscopicall y unidentifiable particular material (can have both benthic and pelag ic orig in, Fig. 4 A).The r elativ e abundance of known benthic diatom genera was low (at most < 0.05% of 18S amplicon reads, Fig. 4 B).The benthic diatom genus Cocconeis significantl y corr elated with wind speed (R 2 = 0.14, P < 0.01), indicating resuspension from benthic en vironments .Further, we searched for other known sediment-associated or ganisms, suc h as Bathyarchaeota , in the 16S amplicon dataset, which were not detected.We also did not detect any enzymes involved in denitrification (nrfA, narG, napA, nirK, nirS, nor and nosZ) in the meta pr oteomes.

Discussion
The succession we observed during the 2018 spring bloom, with an initial dominance of diatoms, follo w ed b y dinoflagellates, using 18S rRNA gene amplicon sequencing a gr ees with the typical phytoplankton community succession at Helgoland Roads (Wiltshire et al. 2008, Käse et al. 2020 ).Meta pr oteome anal ysis highlighted the dominance of eukaryotic proteins in the sampled particles, making up 70-80% of detected proteins, most of which belonged to the major diatom phytoplankton.The high pr e v alence of pr oteins involved in basic metabolism such as energy production and conversion is consistent with the sampled > 10 μm particle fraction mostly comprising of living phytoplankton cells, especially at the beginning of the bloom.Ho w e v er, bacterial pr oteins wer e pr esent The network was calculated based on Spearman correlations (r > 0.7, p < 0.01) exclusively between 18S ASVs and 16S ASVs.Desulfobacterota (y ello w cir cles) w er e onl y associated with the diatom-dominated module I.The sizes of the symbols corr espond to the number of significant correlations of the nodes (degree).(B) The bacterial taxa co-occurring with diatoms and dinoflagellates, respectively, are displayed as horizontal bars with a length proportional to the number of ASVs belonging to eac h linea ge.Desulfobacterota (y ello w) as w ell as potentially sulfur-oxidizing Ectothiorodospiraceae ( Gammaproteobacteria , blue) were positiv el y associated with diatoms but not with dinoflagellates.(C) The pie chart displays r elativ e abundances for all 19 pr otein gr oups classified as belonging to Desulfobacterota during all thr ee time points sampled for meta pr oteomics .Of these , se v en pr otein gr oups (32%, bright y ello w) r epr esented enzymes involv ed in dissimilatory sulfate r eduction.Yello w cir cles indicate which of the k e y enzymes in this metabolic pathway were detected.(A) Secchi depth and abundance of "detritus" (unidentifiable material) along the course of the bloom.Detritus le v els wer e estimated under the microscope on a scale between 0 and 6 (0:"none", 3: "moderate" 6:"massive").(B) Relative abundance of Desulfobacterota (16S rRNA gene) in relation to r elativ e abundances of kno wn benthic diatom taxa ( Cocconeis , Navicula, Nitzschia , 18S rRNA gene). in all of the three selected time points, with a notable peak at day 144, and their taxonomic composition was very similar to that observed via 16S rRNA gene amplicon sequencing.Overall, the composition of particle bacterial microbiomes agreed with other reports from similar environments (Wang et al. 2024, Crump et al. 1999, Schultz et al. 2020, Heins et al. 2021, Reintjes et al. 2023 ).
As hypothesized, bacteria co-occurring with abundant diatoms formed a distinct network module, with taxa that were different from those co-occurring with dinoflagellates in a second distinct module.A thir d netw ork module contained mostly diatom taxa, but also one dinoflagellate ASV and other phytoplankton taxa.These network patterns offer an alternative way to visualize the temporal dynamics of the bloom, and should not be inter pr eted as evidence of physical interactions between taxa (Röttjers and Faust 2018 ).Ho w e v er, one pattern that is striking in our network analysis is the exclusive co-occurrence of Desulfobacterota with diatoms.
We observed a high relative abundance of Desulfobacterota (in total 0.38% of bacterial amplicon reads) in particles ( > 10 μm), especially in the early phase of the bloom when diatoms were dominating.Despite the limited resolution of meta pr oteome data compared to DNA-based methods, desulfobacterial proteins wer e r epr esented at all selected time points (in total 0.14% of meta pr oteomes).Network anal ysis further highlighted the tempor al co-occurr ence of Desulfobacterota with se v er al diatom taxa.This raises the question of the niche filled by these anaerobic bacteria during an active phytoplankton bloom.
Diatoms, such as Thallassiosira spp.and Thalassionema spp., whic h wer e co-occurring with Desulfobacterota in this study (Fig. 3 A), are known to form aggregates (Thornton 2002 ).For example, se v er al species of Thalassiosira extrude long chitin fibrils, whic h pr e v ent sinking and bind exopol ymeric substances (EPS) also produced by the algal cells (Herth andBarthlott 1979 , Den et al. 2023 ).This creates a favor able envir onment for bacteria to attac h, whic h can in turn stimulate algal EPS production (Gärdes et al. 2011 ).EPS mak es particles adhesi ve and thus bacteria can be ca ptur ed in this sticky EPS layer.Smaller particles can a ggr egate to larger particles by collision and adhesion, in particular during phytoplankton blooms with high particle densities.Such aggregates can feature high numbers of living, photosynthesizing algal cells (Thornton 2002 ) which produce ample oxygen during the day, but at night r espir ation by the algal cells and their surrounding bacteria may deplete oxygen sufficiently for low oxygen or e v en anaer obic micr o-nic hes to form.Experiments with marine particles in laboratory roller tanks have demonstrated how oxygendepleted zones can form in diatom-derived particles (Ploug and Berggkvist 2015 ), in part due to EPS associated to the algae rendering the particles impermeable to water flow (Zetsche et al. 2020 ).Inter estingl y, sulfate r eduction was detected in suc h particles in an earlier study, and the reducing microzones where this process was pr esumabl y taking place wer e fr equentl y associated with diatom frustules (Shanks and Reeder 1993 ).
Desulfobacterota have not been identified as frequent members of diatom microbiomes in either cultures or in the field so far (Helliwell et al. 2022 ).Howe v er, a r ecent global surv ey of the diatom interactome detected positive correlations between diatoms and sulfate-reducing bacteria ( Desulfovibrio ) in samples from the Tara Oceans expedition (Vincent and Bowler 2020 ).In addition, Desulfobacterota hav e pr e viousl y been r epeatedl y detected in particleassociated communities in the photic zone (Crump et al. 1999, Liu et al. 2020, Hallstrøm et al. 2022 ).In a par allel study, we r econstructed a genome of a Desulfobacterota member from metagenomic data from material sampled during the same phytoplankton bloom (Wang et al. 2024 ).Ectothiorhodospiraceae ar e pur ple  sulfur bacteria (belonging to Gammaproteobacteria ), that oxidize reduced sulfur compounds as electron donors during anoxygenic photosynthesis and are anaerobic to microaerophilic (Imhoff et al. 2022 ).The y can o xidize H 2 S, e.g.produced via dissimilatory sulfate reduction by members of the Desulfobacterota .Our observation of Desulfobacterota as well as Ectothiorhodospiraceae co-occurring with diatoms is consistent with potential sulfur cycling during the diatom-dominated phase of this phytoplankton bloom under anoxic to very low oxygen conditions.While co-occurrence of diatom and Desulfobacterota rRNA genes does not by itself indicate that sulfate reduction is taking place in diatom-derived particles, detection of k e y functional enzymes by meta pr oteomics suggests that Desulfobacterota were actively carrying out sulfate reduction and thus gaining energy through anaerobic respiration during the bloom.The higher detection of these enzymes especially at the last time point (Julian day 144, Table S1 can likely be attributed to the lar ger pr oportion of bacterial pr oteins at this time point, when the diatom bloom was declining. Tempor al associations, suc h as detected by our co-occurrence network analyses , ha ve to be interpreted with caution as additional variable factors that were not taken into account may influence the observed correlations.Importantly, we cannot quantitativ el y assess the influence of underlying sediment microbiomes at the rather shallow sampling site (6-10 m depth depending on tide), whic h fr equentl y become r esuspended during times of heavy wind ther eby intr oducing anaer obic micr obes to the pela gic environment.Indeed, some of our results point to w ar ds a significant influence of sediment micr obiomes, suc h as the correlation between r elativ e abundance of total Desulfobacterota ASVs and secchi depth, as well as Desulfobacterota and wind speed.The detected Desulfobacterota ASVs (e.g.classifying as Desulfosarcina , Desulfocapsaceae , Desulfobulbaceae ) are related to typical benthic lineages (Rav ensc hla g et al. 1999 ), suggesting that they originated from the underlying sediments .T he diatoms co-occurring with Desulfobacteria were primarily classified as common pelagic genera such as Thalassiosira and Thalassionema .Ho w e v er, the genus Brockmanniella , which can also inhabit benthic biofilms (Hernández Fariñas et al. 2017 ), was also found within network module I.The primarily benthic diatom genera Nitzschia , Navicula and Cocconeis were not part of the network, but Cocconeis r elativ e abundances correlated with wind speed, indicating a turbulence-driven benthic-pelagic coupling.Ne v ertheless, a benthic origin of the Desulfobacterota detected in our study does not exclude physical interactions with the blooming pelagic diatoms.In shallow environments, benthic and pela gic micr obiomes ar e in close contact within the photic zone, and seeding of pelagic particles by benthic microbes is likely frequent, which may explain the observed co-occurrences.In fact, most of the enzymes involved in sulfate r eduction wer e detected in the last timepoint, when the diatom bloom was declining and wind speeds were moderate, indicating low sediment resuspension.Micr oscopic anal ysis r e v ealed high abundances of the ha ptoph yte ph ytoplankton Phaeocystis globosa were also detected at this timepoint (Wang et al. 2024 ), which may have contributed to a ggr egate formation (Schoemann et al. 2005 ) Importantly, our methodology does not allow us to determine the physical proximity of Desulfobacterota and diatom cells.We ther efor e cannot rule out that the observed co-occurrences reflect a pur el y tempor al association between pela gic diatoms and r esuspended sediment bacteria.Indeed, demonstrating a potential physical association between Desulfobacterota and diatoms would r equir e in situ microscopic investigation using for example taxonor gene-specific fluor escent pr obes.With this a ppr oac h, sulfate r educers wer e for example detected in oxygenated surface waters in the Black Sea (Bryukhanov et al. 2011 ).Likewise, it has yet to be clarified, whether any such temporal or physical associations ar e annuall y r ecurr ent, under what specific conditions they occur, and which specific diatom taxa are in volved.T hus , further studies are needed to confirm whether our results are re presentati ve for phytoplankton blooms in shallow water coastal areas.
Our results highlight the complexity of particle microbiomes and corr obor ate the need to study algae-bacteria particles as spatiall y heter ogeneous entities .T he micr obiomes of a ggr egateforming phytoplankton may indeed be similarly complex as those of animals and other multicellular organisms, featuring distinct micr o-nic hes with sub-micr obiomes (analogous to e.g.human skin vs. gut microbiomes), whose compositions depend on the c hemical envir onment on a small spatial scale.Anaer obic micr oniches in phytoplankton-derived aggregates may affect carbon cycling insofar, as sulfate reduction can chemically alter particu-late organic matter, rendering it resistant to microbial degradation (Raven et al. 2021 ).Ho w ever, this effect has been demonstrated in oxygen-deficient zones of the ocean, and it is unclear this can be expected in oxygenated surface waters.Sulfate reduction in marine particles has been predicted to be prevalent in vast hypoxic ( < 60 μM oxygen concentration) waters of the global oceans, based on modeling of particle properties and oxygen regimes (Bianchi et al. 2018 ).The model could explain observed precipitation of Cadmium sulfide (CdS) in these waters (Janssen et al. 2014 ).Surface waters in the photic zone at our sampling site are expected to feature far higher oxygen concentrations at around 300 μM (7-10 mg/l, https:// dashboard.awi.de/?dashboard=34404 ).Howe v er, anaer obic nic hes in marine phytoplankton-deriv ed particles wer e alr eady fr equentl y r eported in the context of nitrogen fixation (i.e.diazotrophy, e.g.Riemann et al. 2022 ).Recent w ork has sho wn that many Desulfobacterota are in fact capable of both diazotrophy and sulfate reduction (Liesirova et al. 2023 ) and that diazotrophic Desulfobacterota can move to w ar ds phytoplankton-deriv ed or ganic matter via chemotaxis.Similar to diazotr ophs, sulfate r educers in turn provide additional niches for bacteria which scavenge their metabolic products , i.e .sulfur oxidizers such as the observed Ectothiorhodospiraceae , leading to higher metabolic diversity and complexity within particles.We suggest that incor por ation of sulfate reducing bacteria in phytoplankton-derived aggregates may be a r ele v ant phenomenon, especiall y in shallow coastal areas where the microbiomes of underlying sediments provide a pool of abundant sulfate reducers to colonize aggregates.Considering the significance of marine particle microbiomes in carbon sequestration, it is vital to understand the consequences of their metabolic and compositional complexity, and the resulting microbial interactions in these microbiomes.

Figure 1 .
Figure 1.Pr ogr ession of the 2018 spring bloom by Helgoland, North Sea.(A) Chl a (green) peaked around the 26th of April (Julian day 116), coinciding with a drop in Secchi depth (brown).Wind speed maxima (grey) reached storm levels (German: "Sturm", > 75 km/h, dotted line) on three occasions in the first half of the bloom.Points indicate water le v els (m from global reference zero point) and tidal direction (white points: falling, black points: rising).Asterisks indicate time points (Julian days) for which metaproteomic sampling was performed.(B) 18S rRNA gene amplicon sequencing r e v ealed that diatoms ( Bacillariophyta ) were the most abundant phytoplankton lineage, although dinoflagellates ( Dinophyceae ) dominated after the chl a peak.(C) 16S rRNA gene amplicon sequencing showed the highest r elativ e abundances of Desulfobacterota (yellow) during the first half of the bloom.
B, and in more detail in supplementary Fig.1.Ele v en ASVs belonging to the Desulfobacterota were part of the main diatomdominated network module I and co-occurred exclusively with diatoms, including the genus Desulfosarcina , and the families Desulfocapsaceae and Desulfobulbaceae as well as the lineage Sva1033(Rav ensc hla g et al. 1999 ).In addition, 4 ASVs classified as Ectothiorhodospir aceae (genus Thiogr anum ) also co-occurred with diatoms (Fig.3 B).

Figure 2 .
Figure2.Taxonomic and functional annotation of particle meta pr oteomes fr om thr ee selected time points during the bloom.(A) The pr oportion of bacterial proteins increased during the course of the bloom (B) Eukaryotic proteins made up the majority of identified proteins and sho w ed an initial strong dominance of Bacillariophyta (C).Functional annotation of metaproteomes indicated a high contribution of metabolism-related eukaryotic pr oteins earl y in the bloom, while cellular pr ocesses and signaling incr eased during later time points.

Figure 3 .
Figure3.Co-occurrence of eukaryotes and bacteria, and detected protein groups of Desulfobacterota .(A) A network analysis of 18S rRNA (squares) and 16S rRN A (cir cles) gene ASV co-occurrences resulted in two major network modules containing diatom [I] and dinoflagellate [II] 18S ASVs, r espectiv el y, as well as one mixed module[III].The network was calculated based on Spearman correlations (r > 0.7, p < 0.01) exclusively between 18S ASVs and 16S ASVs.Desulfobacterota (y ello w cir cles) w er e onl y associated with the diatom-dominated module I.The sizes of the symbols corr espond to the number of significant correlations of the nodes (degree).(B) The bacterial taxa co-occurring with diatoms and dinoflagellates, respectively, are displayed as horizontal bars with a length proportional to the number of ASVs belonging to eac h linea ge.Desulfobacterota (y ello w) as w ell as potentially sulfur-oxidizing Ectothiorodospiraceae ( Gammaproteobacteria , blue) were positiv el y associated with diatoms but not with dinoflagellates.(C) The pie chart displays r elativ e abundances for all 19 pr otein gr oups classified as belonging to Desulfobacterota during all thr ee time points sampled for meta pr oteomics .Of these , se v en pr otein gr oups (32%, bright y ello w) r epr esented enzymes involv ed in dissimilatory sulfate r eduction.Yello w cir cles indicate which of the k e y enzymes in this metabolic pathway were detected.