Control of Vibrio vulnificus proliferation in the Baltic Sea through eutrophication and algal bloom management

.

The Baltic Sea is a semi-enclosed marginal sea of the Atlantic located in northern Europe, with a coastline of approximately 8000 km and covering an area of 415,266 km 2 .Saline inflows through the North Sea produce a 2000 km long lateral surface salinity gradient throughout the entire Baltic Sea, ranging from high salinities (>25) in the transition zone of the Kattegat to low salinities (<5) in the Gulf of Bothnia 1 .The Baltic Sea is characterized by an estuarine-like circulation due to the positive freshwater budget.The drainage area of the Baltic Sea encompasses a population of approximately 85 million, and is consequently heavily influenced by eutrophication 2 .In addition, annual mean sea-surface temperatures are rising 2 and the ecosystem is expected to be increasingly affected by warming in the coming decades 3 , and will be faced with extended heat wave durations 4 .
These changes favour the growth of pathogenic bacteria of the genus Vibrio and an increase in Vibrio spp.abundances, infection rates, and fatal cases along the Baltic Sea coastline has been reported [5][6][7][8] .The infections of predominantly immunodeficient humans can be associated with the consumption of raw or undercooked shellfish, but in the Baltic Sea, they frequently manifest as skin infections resulting from direct contact with coastal brackish water.Only a low number of infections are currently associated with Vibro vulnificus in the Baltic Sea 9 , but these are usually severe and often lethal 6,10,11 .
Temperature and salinity are widely accepted as the two primary regulators of V. vulnificus abundance and distribution [12][13][14][15] , but factors related to eutrophication, such as elevated dissolved organic carbon [DOC 16 ] concentration or dinoflagellate blooms 17 , have been shown to stimulate V. vulnificus growth in laboratory settings.Due to its preference for intermediate salinities and proliferation at water temperatures >15 °C18 , V. vulnificus experiences optimal growth conditions during the summer in the Baltic Sea.The expected spread of V. vulnificus does not only pose a significant threat to public health but also to the tourist, fishing, and aquaculture industries [19][20][21] .In consequence, the question arises, whether measures can be taken to regulate V. vulnificus abundances.
In a natural setting, macrophytes such as seagrass (Zostera marina) might reduce the abundance of pelagic and potentially pathogenic Vibrio spp.and other pathogens within the seagrass canopy 22,23 .The underlying mechanisms for this decrease are elusive but could include increased sedimentation rates due to hydrodynamic attenuation [24][25][26] , filter-feeding by benthic fauna 27,28 , or allopathic chemicals exuding directly from the seagrass plants 29 .Regardless of the specific causal mechanisms through which Z. marina beds potentially impact V. vulnificus abundance, they are putative nature-based solutions for reducing V. vulnificus 23 .Consequently, this study aims to elucidate important factors for mitigating V. vulnificus abundances along the Baltic coast, evaluating the potential of seagrass as a nature-based solution at the local scale and the reduction of eutrophication at the regional scale, respectively.We measured a large array of physical, biological, and hydrochemical parameters, and used three parallel methods (cultivation, amplicon sequencing, and droplet digital PCR; ddPCR) for the quantification of V. vulnificus within, adjacent and far from seagrass meadows.This enabled us to simultaneously assess the relationship between seagrass, environmental factors, such as temperature, salinity, and (in)organic nutrients, eukaryotic and prokaryotic microbial communities, and V. vulnificus over a vast salinity range, facilitating the detection of more generalizable patterns.These patterns are applicable to a wide range of environmental conditions found in estuaries and marginal seas worldwide.By combining a large number of parameters with a machine learning approach, complex and non-linear relationships are identified between the environment and V. vulnificus abundance.
V. vulnificus is likely stimulated by eutrophication-induced algal blooms; therefore, reducing eutrophication appears to be a promising strategy for mitigating the health risks associated this bacteria.Based on the results of this study, a regulatory effect of seagrass on V. vulnificus could not be determined.
Macrophytes had no significant effect on V. vulnificus abundance V. vulnificus was detected in 47% of the samples using ddPCR, in 33% through 16S rRNA gene sequencing, and both methods identified it in 20% of the cases.Across all sampled stations, local effects of Z. marina or Fucus spp. on V. vulnificus absolute or relative abundances were insignificant.V. vulnificus was equally distributed between substations A, B, and C (Wilcoxon rank sum test (Fig. 2)).This was consistent for vvhA gene copy numbers and 16S rRNA gene relative abundances, green colony forming units (CFUs), which are presumed, but not unequivocally identified, V. vulnificus colonies, and for heterotrophic bacterial cell counts (p-values > 0.05).Likewise, there was no consistent pattern between substations A, B, and C in vvhA gene copy number nor relative 16S rRNA gene abundances when observing the individual stations (Supplementary Figs. 2 & 3).V. vulnificus in sediments was also evenly distributed among substations (Supplementary Fig. 4).Based on the vvhA genes and 16 S rRNA gene sequencing, V. vulnificus was not present in Z. marina biofilms.However, one V. vulnificus strain was isolated from seagrass leaves at station BV-12.In addition, other potential pathogenic Vibrio spp.found in the Baltic Sea were also distributed equally between substations on a Baltic Sea-wide level (Supplementary Fig. 5).
Eutrophication impacted V. vulnificus abundance A positive relationship was observed between temperature and V. vulnificus abundance (approximated by vvhA gene copies mL −1 , Fig. 3).Within the temperature range of 19.5-20.5 °C, however, strongly eutrophied samples exhibited higher vvhA gene numbers compared to other samples in the same, as well as other, temperature ranges, coinciding with higher cell abundances of heterotrophic bacteria in general (Fig. 3).The average concentration of the vvhA gene in this temperature bin was higher than that at higher temperatures, and, on average, the bin exhibited higher eutrophication.The eutrophication index captured more of the variance in V. vulnificus abundance (36%) than temperature (12%) (Supplementary Fig. 6).The combination of both explained 43% while the more comprehensive and non-linear random forest (RF) including all parameters explained 59% of the variance.

Prevalence of V. vulnificus was correlated with indicators of blooming situations
Both RF models exhibited moderate performance, explaining 0.59 ± 0.32 and 0.59 ± 0.45 of the variance in the derived total V. vulnificus abundances Fig. 2 | V. vulnificus abundance and microbial abundances at three substations summarized across all sampled locations.a V. vulnificus abundances (vvhA gene copies mL −1 ), b Relative V. vulnificus abundance based on 16 S rRNA genes, c Presumptive V. vulnificus based on green CFUs mL −1 on TCBS agar, d Number of heterotrophic bacteria mL −1 .Station A is located within a macrophyte meadow, while station B is located 15 m, and station C is located 100 m from the meadow edge.Boxes represent the 25th and 75th percentiles.and of the variation of vvhA gene copies mL −1 , respectively.Root Mean Square Error (RMSE) was 697.16 ± 597.50 and 6.33 ± 4.66 of the target values 8,130.20 and 37.93, and the mean absolute error (MAE) was 519.76 ± 307 and 4.56 ± 2.55.

Discussion
We present a comparative study of 19 coastal Baltic Sea sampling stations in summer/early fall 2021, exploring the potential of seagrass as a local solution and the reduction of eutrophication as a more regional strategy to mitigate the proliferation of potentially pathogenic V. vulnificus.For these stations, physical-, biological-and hydrochemical parameters were measured and variables explaining V. vulnificus occurrence were identified.We found that eutrophication-related parameters such as POC/PON and high Cyanobium sp. and Synechococcus sp.abundance predicted high V. vulnificus abundances, while the occurrence and density of Z. marina showed no predictive value.This implies that reducing eutrophication on a regional level could be a promising strategy for constraining further proliferation of V. vulnificus along the Baltic Sea coast.
The distribution and abundance of V. vulnificus has been extensively studied, with temperature and salinity being consistently identified as key factors [12][13][14][15] .Our study observed temperatures exceeding 15 °C and salinities ranging from 6 to 16, which are generally considered favorable conditions for V. vulnificus 5,30 , yet abundances varied between stations.Notably, stations with a high eutrophication index tended to exhibit higher V. vulnificus abundances compared to those with similar or even higher temperatures but a lower eutrophication index (Fig. 3).Furthermore, the eutrophication index explained a greater proportion of the variance in the V. vulnificus abundance than temperature (Supplementary Fig. 6).This finding is consistent with earlier studies, which also related elevated concentrations of DOC and chl-a to high V. vulnificus abundances 16,30 .Chl-a was, besides being part of the eutrophication index, also closely associated with the predictors for V. vulnificus identified by RFE, namely POC, PON, Cyanobium sp. and Synechococcus (relative) abundance in this study (Supplementary Fig. 7).Additionally, the effectiveness of heterotrophic cell counts (LNA cells in Fig. 4) in the prediction underscored the link to eutrophication.These results align with observations from the Neuse River estuary, which faces similar problems with eutrophication as the Baltic Sea, where elevated levels of both total Vibrio spp.and V. vulnificus corresponded with high phytoplankton biomass 31 and an investigation after Hurricane Ian (October 2022), where high phytoplankton mass, approximated by chl-a concentration, and associated zooplankton abundance, was found to potentially stimulate V. vulnificus proliferation 32 .
Cyanobium sp., which thrives under eutrophic conditions and is known to form blooms and correlate positively with the decay phases of other blooms 33 , was found to be an important V. vulnificus predictor.The identified ASVs for Cyanobium sp.comprised up to 18.7% of the total prokaryotic sequence reads at some of the highest vvhA gene copy numbers.Associations between pathogenic Vibrio spp.and (harmful) algal blooms have been observed before 34 and our results substantiate this relationship.It is plausible that eutrophication, by stimulating blooms, served as an important indirect driver of the distribution and abundance of V. vulnificus via the release of DOC/DON and POC/PON.The idea that the increase in DOC/DON triggered by the bloom contributed to the proliferation of V. vulnificus prior to sampling is supported by the predictive power of the actinobacteria Candidatus Aquiluna (ASV236), ML602J-51 (ASV654), and Nocardioides (ASV1070).Actinobacteria are commonly associated with the demise of phytoplankton blooms and the generation of cyanobacterially sourced DOC 35,36 .The stimulation of V. vulnificus abundance by algal blooms was also inferred from the high V. vulnificus abundances coinciding with high filamentous Ceramium sp.abundances, peaking at 9.2% of the sequence reads during both the highest vvhA gene concentrations and derived V. vulnificus cell counts at station BV-08-A (Fig. 4).Ceramium sp. is an opportunistic fast-growing algae, with high nutrient uptake rates 37 .Interestingly enough, the eutrophication index, POC and PON are all low at this station.This could be an indication that the blooming (in terms of relative abundance) of specific plankton, for example, Ceramium sp. and Cyanobium sp. can already convey some benefit.POC/PON are still good predictors of V. vulnificus in that scenario, as they generally increase with phytoplankton blooms, including all the blooms that provide speciesspecific benefits.Regardless of the specific combination of factors that promote V. vulnificus proliferation, reducing eutrophication would limit the extent of algal blooms in the Baltic Sea, likely reducing infection risk.Among the top five instances of highest vvhA gene copy abundances and derived V. vulnificus cell numbers, either Ceramium or Cyanobium sp. was consistent, except at station BV-24-B, found at high relative abundance.
The high relative abundance of Cyanobium sp. was concurrent with other species expected to either infect or co-occur with a bloom.The parasite Aplanochytrium sp., which is known to bloom 38 and to correlate strongly with phototrophs 39 , and the fungi Rhizophydium sp.(ASV 312) 40 , reached their highest relative abundance (max 0.4% and 1%, respectively) at stations with high Cyanobium sp. and V. vulnificus (relative) abundance.Under conditions conducive to extensive bacterial growth, such as during blooms, Thaumatomastix (ASV320) might be able to temporarily compete with other organisms 41 .Thaumatomastix exceeded 0.1% of the eukaryotic community only in samples where Cyanobium sp., Synechococcus, and their potential parasites were relatively abundant, further pointing to a highly productive system.This scenario coincided with elevated V. vulnificus abundance.
We speculate that when grazing pressures on the algal bloom reached critical levels, the bloom ceased, leading to a depletion of the organic nutrients crucial for the proliferation of V. vulnificus.In instances where we observed high relative abundances of the potential predator Teleaulax (>1%), there was a concurrent decrease in Synechococcus and Cyanobium sp.(relative) abundances and a decrease in the concentrations of PON and POC (e.g.BV-18 & BV-20).This phenomenon was consistently associated with the absence or low relative abundance of V. vulnificus.This may be attributed to the potential role of Teleaulax in depleting the Synechoccocus population 42,43 and thereby reducing the food sources for V. vulnificus.
While organic nutrients were good predictors of V. vulnificus abundance, inorganic nutrients were not.Additionally V. vulnificus was generally less correlated with inorganic nutrients (e.g.PO 4 3-) than heterotrophic bacteria overall (Supplementary Fig. 8).This suggests that potential bloomforming organisms may play a pivotal role in the chain leading from inorganic nutrient input to increased V. vulnificus abundance by converting inorganic nutrients into organic forms.
In addition to the increase in organic constituents present during a phytoplankton bloom, the higher abundance of phytoplankton and the associated zooplankton community could potentially offer protection against grazing by bacterivorous protozoa, as previously demonstrated for V. cholerae 44 .V. cholerae population growth is normally balanced by protozoan predation, but during declining blooms, V. cholerae has been shown to multiply at an increased rate due to the higher availability of DOC and DON and/or change to a particle-attached lifestyle for protection from predation 45 .A similar protective effect has been shown for other species of Vibrio, which significantly benefitted from association to phytoplankton during periods of intense grazing 46 .We hypothesize that V. vulnificus might also benefit from a similar interaction.
The discussed microorganisms and eutrophication index exhibited greater explanatory power than temperature for V. vulnificus abundance.This was likely due to their sensitivity to historical and current temperatures and nutrient availability.As a result, they served as a record of the water mass that was more predictive than temperature alone.The suggested pathway through which eutrophication could stimulate V. vulnificus abundance is summarized in Fig. 5.The identity of the blooming organism and its predators might not be as critical as the fact that the bloom and its associated predators stimulate organic nutrient availability and might provide a degree of protection.
While eutrophication is often a regional problem, with regional solutions, such as the Baltic Sea action plan 47 .Intraregional differences exist between bays due to differences in nutrient sensitivity 48 and nutrient input from point sources 49 .Engineering solutions to eutrophication in semienclosed Baltic Sea bays have been tested successfully 50 , potentially providing a solution for some Baltic Sea coasts.Besides allowing for a potential (multi)-national approach to reduce the V. vulnificus hazard, the possible impact of eutrophication can be used to improve risk assessment and early warning by integrating the existing HELCOM Eutrophication Assessment Tool 48 into the currently available ECDC vibrio viewer, which relies solely on temperature and salinity, for the Baltic Sea 51 .In other estuaries or marginal seas, the integration of chl-a data from remote sensing into early warning systems may provide an improved prediction of V. vulnificus abundance.
The assessment of V. vulnificus predictors was hampered by the observed discrepancies between the quantification methods employed.They are likely due to multiple biases associated with each methodology.These are, in no particular order: CFUs by nature of the method, underestimate the real abundance, as they are unable to quantify viable but not culturable cells, resulting in the lowest abundances measured by any of the methods (Fig. 2).Differences between vvhA gene copy numbers and derived cell counts can be partly attributed to multiple copies of the 16S rRNA gene 52 , which range between 8 and 12 for V. vulnificus 53 .We did not correct for these, according to the advice of Louca et al. 54 .Additionally, not all V. vulnificus cells contain the vvhA gene, which further complicates the comparison.The discrepancy between vvhA and derived V. vulnificus abundance could be an indication that a sizeable part of the sampled V. vulnificus community lacks this pathogenicity gene.Further differences Fig. 5 | Overview of the suggested pathway through which eutrophication promotes V. vulnificus proliferation.P and N represent the total dissolved inorganic phosphorus and nitrogen.Inorganic nutrient inflow from the land induces algal blooms, providing the organic material required for V. vulnificus proliferation and potentially protection from bacterivorous protozoa.The figure also displays other organisms identified as important predictors and their potential connections to blooms.Notably, Zostera marina could not be shown to significantly affect abundance in our study.might arise from primer bias, PCR biases and biases in quantification of bacterial counts.When assessing the vvhA gene in our study, fewer instances of V. vulnificus are identified per station; however, V. vulnificus is detected at more stations.A prime example is station BV-11-C (Fig. 4), which has the second highest vvhA gene copy numbers, and where presumable V. vulnificus colonies were cultured, but no V. vulnificus was detected with the 16S rRNA amplicon sequencing, while at stations BV-11-A and BV-11-B V. vulnificus is detected by all methods.In this case, the discrepancy arose because the 16S rRNA gene sequences amplified at station BV-11-C could not unambigiously be identified as V. vulnificus.A similar degree of similarity was observed with a related Vibrio species.Additional discrepancies might arise from V. vulnificus variants not accounted for in the database used for comparison.This complicates the comparison between the methods, but overall, the association between algal blooms and V. vulnificus could still be shown across all measured parameters.
V. vulnificus was found to be absent from seagrass leaves by molecular methods and only one strain could be isolated from them, even though V. vulnificus is known to form biofilms 55 under a variety of conditions [56][57][58] .This is surprising given the release of labile DOC from seagrass leaves 59 but could be due to allelopathic chemicals exuded from the seagrass leaves hampering the growth of V. vulnificus, as observed for a general bacterial community 60 .An alternative explanation is that there is a strong association between the seagrass and their associated microbiome and that V. vulnificus was consistently outcompeted 61,62 .V. vulnificus was, however, found in the benthic and pelagic environment, in similar abundances at substations inside and outside of the meadows (Fig. 2, Supplementary Figs. 1, 2), implying that the impact of seagrass meadows on V. vulnificus abundance, relative abundance, and culturable potential within the surrounding environment was minor in the Baltic Sea.Additionally, no effect on the relative abundance of other potential Vibrio pathogens could be detected within the sampled salinity range (Supplementary Fig. 5).Similar to the findings in San Diego coastal waters, our results contrast the reduction in pathogenic bacteria abundance observed by Lamb et al. in 2017 63 .Future findings may differ for other conditions, e.g., water temperatures, salinity ranges, or less eutrophied systems, which were outside our study framework.The discrepancy between our study and earlier work on the influence of seagrass on potential pathogens may also be ascribed to a difference in methodology or study area.Reusch et al. (2021) used solely chrom-agar based blue colony counts at salinities between 14.3 and 22.1 23 , while we employed methods more specific to V. vulnificus over a different salinity gradient.The salinity of our study areas was mostly lower and may have affected the physiology of seagrass and its possible effect on V. vulnificus.Based on our results, it could not be shown that seagrass meadows are a direct effective nature-based solution for reducing V. vulnificus abundance and associated infections in the Baltic Sea.However, the current high nutrient concentrations in the Baltic Sea may have masked possible direct effects of seagrass, and probably also Fucus spp., on V. vulnificus, and indirect effects through the reduction of the nutrient load.
Our study suggests that a reduced anthropogenic nutrient input could mitigate the cascade of events that stimulate algal blooms and through that, organic nutrient availability.This in turn, would make conditions less conducive for V. vulnificus proliferation.This could reduce public health risks and be a critical management tool for V. vulnificus in the Baltic Sea, as well as in other coastal brackish water systems worldwide.

Materials and Methods
Sampling stations and SCUBA sampling strategy Samples from fifteen Z. marina and, for comparison, four Fucus spp.stations were collected from July 25th to September 2nd, 2021, covering the German, Estonian, Finnish, Polish, Swedish, and Danish coasts (Fig. 1).Water and sediment samples were collected inside the macrophyte meadows (substation A), as well as the leaves/fronds of the macrophyte.Analogously, water and sediment samples were collected at control stations without macrophytes, located 15 m (substation B) and 100 m (substation C) from the edge of the meadow (Fig. 1).Sampling was carried out by SCUBA or snorkel divers at water depths ranging from 0.6 to 4.7 m.Sampling was only performed up to a moderate breeze of 4 on the Beaufort scale with wave heights <0.5 m to prevent heavy mixing between substations.The prevailing current direction was visually determined at the sampling depth by releasing food colourant (Dr.Oetker 4 Back & Speisefarbe).Divers approached the sampling site against the current and maintained buoyancy to prevent sediment resuspension.
Water samples were collected following two methods.First, for DNA extraction, flow cytometry, and plating, 100 mL of water was collected in 9 replicates using rinsed syringes.The samples were collected ca. 5 cm from a macrophyte and ca.20 cm above the sediment at substation A, and ca.20 cm above the sediment at substations B and C. Second, sterile 5 L plastic bags, used for chl-a, inorganic nutrients, POC and PON, DOC and DON determination, were filled by swimming against the current while keeping the bag open 20 cm above the sediment.This water was split between 1 L plastic and glass bottles after triplicate rinsing.Sediment was collected in nine sterile 50 mL Falcon tubes by scraping the top 1 cm of sediment.At substation A, Z. marina/Fucus spp.samples were collected in nine replicates by pinching the connection between the root and stalk.Subsequently, the leaves/fronds were separated by hand for further processing.All samples were transferred to a 4 °C cooler immediately and stored (maximally 8 h) until processing.

Macrophyte characteristics
Z. marina plant densities were determined by counting 20 × 20 cm squares in triplicate.The length of 30 leaves was measured per meadow.Fucus spp.densities were not determined due to the large heterogeneity within the individual beds.

Environmental parameters
Salinity, temperature, and depth were measured using a CTD48M (Sea & Sun Technology, Trappenkamp, Germany) attached to a diver.DO and pH were measured using a handheld multimeter (HQ40D Portables 2-Channel Multimeter) (HACH, Iowa, USA) with Hach Intellical LDO101 and Intellical PHC101 probes, respectively.
Water from the plastic bottles was filtered through 25Ø GF/F filters (Whatman plc, Maidstone, UK), with the filters stored for chl-a determination at −80 °C (250-500 mL) and the filtrate for nutrient analysis stored at −20 °C.Water from the glass bottles was filtered through pre-combusted (450 °C, 4 h) 25Ø GF/F filters for POC and PON determination and the filtrate for DOC and DON analysis, all stored at −20 °C.PO , and SiO 2 concentrations were measured using a Seal Analytical QuAAtro automated constant flow analyzer (SEAL Analytical Ltd, Nordestedt, Germany), with detection limits of 0.1, 0.2, 0.05, 0.5 and 2 μM, respectively 64 .DOC and DON were analyzed with a TOCL-CPH/ TOC-VCPH TOC-Analyzer (Shimadzu, Germany) 65 , and POC/PON measurements were conducted with a varioMICRO cube element analyzer (Elementar Analysensysteme, Germany).Chl-a was measured fluorometrically using a 10-AU-005-CE fluorometer (Turner, San Jose, USA) according to the HELCOM 66 guidelines and corrected for phaeopigment.

Sediment grain size
The grain size was analysed after 5 d of lyophilization with a Delta 1-24 LSCplus (Martin Christ Gefriertrocknungsanlagen, Osterode am Harz, Germany).Sediment samples were sieved through a 3.5 mm sieve and visible biological elements were removed.The remaining sediment chunks were broken by grinding for a minimum of 90 s with mortar and pestle.Samples were measured using the dry cell of a Mastersizer 3000 with a range from 0.01 to 3.500 μm (Malvern Panalytical, Malvern, UK).

Microbial analyses
Vibrio spp.CFUs and V. vulnificus isolates were obtained from water, sediment, and macrophytes from six independent replicate samples.For this, water aliquots of 50, 100 or 200 μL were (a) plated directly onto Vibrio selective thiosulfate citrate bile sucrose (TCBS) agar (Merck, Darmstadt, Germany) and (b) aliquots of 2, 5, 10 or 25 mL were filtered onto 0.2 μm PC-filters (Merck-Millipore, Burlington, USA) and placed onto TCBS agar.The sediment samples were homogenized after the removal of the overlying water, a subsample of approximately 10 g (dry-weight determined accurately after lyophilization) was transferred from six sediment samples to sterile 50 mL falcon tubes, where 40 mL of double 0.2 μm-filtered station water was added.To detach bacteria, five ultrasonic pulses of 10 s at 25% capacity at 5 s intervals using the Bandelin SONOPULS HD 2200.2 (Bandelin, Berlin, Germany) were applied.After subsequent vortexing and settling of the sediment, water aliquots of 50, 100, or 200 μL were plated on TCBS agar in six biological replicates.The same method was applied to the macrophytes.The only modification was that 5 − 20 mL of the supernatant was filtered over a PC-filter to obtain colonies.After 24 h incubation at 37 °C, CFUs of green colonies were determined for all plates.
To isolate V. vulnificus, green colonies from stations BV-05, BV-06, BV-09, BV−10, BV−12, BV−15, BV-20, and BV-21 were further cultivated on CHROMagar_vibrio™ (Chromagar Ltd.Paris, France) for 24 h at 37 °C and blue-colored colonies were restreaked on TCBS agar, CHROMagar, and Columbia sheep blood agar (Oxoid, Basingstoke, UK).DNA of isolates which were green on TCBS, blue on CHROMagar, and confirmed to be pure cultures on blood agar, was extracted (DNeasy Blood and Tissue Kit, Qiagen, Hilden, Germany), and purity and concentration were determined using a NanoDrop Spectrophotometer (Thermo Fisher, Waltham, USA).To confirm V. vulnificus, the species-specific vvhA gene sequence was targeted using multiplex real-time PCR (5´nuclease assay).In the same assay, an internal amplification control (KOMA) was detected 67 .Primers for the detection of the vvhA gene 45 are provided (Supplementary Table 3).

Enumeration of heterotrophic bacterial cells and autofluorescent phytoplankton
Four mL of water from three replicate 100 mL syringes, which were also used for the DNA extraction, was pipetted in duplicate into 5 mL sterile cryovials (VWR, Radnor, USA), and 200 μL formaldehyde (37%) was added.After homogenization and 1 h incubation at 8 °C, the samples were shock frozen and stored at −80 °C until flow cytometry analysis.Autofluorescent and heterotrophic cells were counted separately using a Cyto-FLEX S flow cytometer (Beckman coulter, Brea, USA) and unstained and SYBR Green I (Invitrogen, Waltham, USA) stained samples, respectively 68 .Autofluorescent cells were analysed with CytExpert software (Beckman coulter, Brea, USA) and grouped into Synechococcus, picoeukaryotes, and nanoeukaryotes.Stained samples were gated and clustered using FlowClust 69 after quality control by FlowAI 70 based on green and red fluorescence and side scatter 68 .The minimum Bayesian Information Criterion was reached for three clusters, showing beads, HNA, and LNA cells.The flow rate was monitored by adding beads to each sample (NFPPS-52-4K, Spherotech, Lake Forest, USA).

Molecular prokaryotic and eukaryotic community composition
Water (92 mL) from three syringes was filtered through 0.22 μm polyvinylidene fluoride membrane filters (Merck, Darmstadt, Germany), which were shock-frozen in liquid nitrogen and stored at −80 °C for downstream molecular analysis.Three falcon tubes with sediment from each substation and three falcon tubes with macrophytes from substation A were frozen at −80 °C.DNA extractions were performed using the DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany).The sediment and macrophyte samples were thawed and homogenized.A subsample of 500 mg of sediment or 4-6 pieces of a macrophyte of 2-3 cm lengths were transferred into bead-beating tubes.These were placed on ice and sonicated twice for 7 min and beadbeaten for 30 s at 4 m/s.Subsequently, the manufacturer's instructions were followed and DNA yield was quantified using Picogreen (Thermo Fisher, Waltham, USA).
16S rRNA genes were amplified using primers, covering the prokaryotic V3-V4 hypervariable region 71 , and 18S rRNA genes by targetting the eukaryotic V4 region (Supplementary Table 3) 72 , with a PCR protocol derived from Latz et al. 73 .(Supplementary Table 4).Illumina sequencing adapters were included in the 5' ends of the 16S and 18S primers (Supplementary Table 3).Phased primers 74 were used to increase the complexity of the sequencing libraries.For 16S rRNA genes, forward and reverse primers were phased (CTAGAGT, TAGAGT, etc. for the forward and ACTACTG, CTACTG, etc. for the reverse primer), for 18S rRNA only the forward primer (ATG, TG, G, or no base) was phased.PCR thermal conditions and master mix details are provided (Supplementary Tables 4 &  5).Leftover adapters were removed using the MagSi-NGS PREP Plus Kit (MDKT00010075, magtivio BV., Nuth, the Netherlands).The purified product was indexed through a second PCR (Supplementary Table 4) following the Adapterama indexing scheme 75 , pooled in equimolar ratios, and sequenced on MiSeq for 16S and 18S rRNA gene metabarcoding by Sci-LifeLab/NGI (Solna, Sweden).In addition to the MiSeq sequencing, the pelagic 16S rRNA gene libraries were deep-sequenced on NovaSeq 6000 (Illumina Inc, San Diego, CA, US).
Phased primer sequences were removed from the reads using a snakemake pipeline 76,77 .The pipeline (https://github.com/biodiversitydata-se/amplicon-multi-cutadapt) encompassed: the removal of read-pairs that contain Illumina adapters, exclusion of read-pairs lacking the expected primer sequences located at the 5' ends of the reads, removal of the primer sequences from the remaining reads and elimination of read-pairs that have misplaced primer sequences.DADA2 78 was used for denoising, concatenating paired-end reads, and chimera removal.The resulting amplicon sequence variants (ASVs) were taxonomically assigned with DADA2 using PR2 79 (V.4.14.0) as a training set for 18S rRNA gene amplicons and Silva 138.1 80 for 16 S rRNA gene amplicons.
Species-level classification of Vibrio ASVs was achieved by sequence comparison (BLAST, V. 2.13.0 81 ) to a custom database.This database includes 16S rRNA gene sequences of complete Vibrio genomes from RefSeq 82 (51 species, 317 strains-including 22 V. vulnificus strains), 41 draft V. vulnificus genomes from clinical isolates from the Baltic Sea region 11,83 , and 84 draft V. vulnificus genomes from environmental Baltic Sea isolates.In order for an ASV to obtain a species-level assignment to a Vibrio spp., we required perfect (100% identity) alignment of the full ASV to 16 S rRNA genes of a single species in the custom database.None of the ASVs that perfectly matched to V. vulnificus 16S rRNA genes also matched perfectly to 16S genes of other species.
In order to obtain derived V. vulnificus cells mL −1 , the integrated relative abundace of V. vulnificus ASVs was multiplied by the flowcytometry measured heterotrophic cells per mL −1 .
Quantification of V. vulnificus V. vulnificus was quantified using ddPCR (QX200™ Droplet Digital PCR System, Bio-Rad, München, Germany) targeting the vvhA gene 84 (Supplementary Table 3).PCR thermal conditions and master mix details are provided (Supplementary Tables 4 & 5).Droplets were generated using the QX100 droplet generator (Bio-Rad, Hercules, USA).Emulsified samples were transferred to a 96-well plate and sealed by a pierceable foil hot seal (BioRad, 181-4040) using PX1 PCR Plate Sealer™ (Bio-Rad, Hercules, USA) (5 s at 180 °C).The PCR was performed using a Bio-Rad C1000 Touch™ thermal cycler (Supplementary Table 4).Subsequently, the plate was analyzed with the QX100 droplet reader using the Quantasoft 1.74.09.17 software.Positive and negative controls were 50 ng of DNA of a V. vulnificus isolate and of a V. harveyi isolate, respectively.In the contamination control, template DNA was substituted with DEPC water.

Statistical analyses and machine learning
To assess the Baltic Sea-wide differences in V. vulnificus abundance between substations A (N = 42), B (N = 45), and C (N = 42) of seagrass stations, a two-sided Wilcoxon rank sum test, assuming independence and equal variance, using the rstatix (V.0.7.2) package in R (V. 4.3.0) 85was performed and holm correction for multiple testing was applied.Additionally, heterotrophic bacterial cell abundances and green CFUs, inside and outside of the seagrass meadows, were compared.Sample sizes for comparison between substations are consistent for the water column and sediment.No repeat measurements were performed, every data point is a distinct sample.
Three separate linear models for V. vulnificus were used to compare the explanatory power (R 2 ) of the traditional predictor temperature with the more integrated predictor eutrophication index and their combination, using averaged vvhA gene copies mL −1 (log10 transformed) per substation as a response.
To identify the environmental conditions associated with high V. vulnificus abundance, two random forests (RF) in combination with a recursive feature elimination (RFE) algorithm were employed in caret (V.6.0.94) 86predicting both the vvhA gene copies mL −1 and the derived V. vulnificus cells mL −1 .These models used sequencing and environmental data as predictors.Combining RF and RFE addressed correlated variables in this high-dimensional dataset and was able to detect complex and nonlinear relationships between measured variables and V. vulnificus 87,88 .
Prior to model training, the 16S and 18S rRNA relative gene proportions were pre-processed: Rare (absent in at least 25 samples) for both datasets were removed.All ASVs taxonomically classified as Vibrio were removed from the predictors.After pre-processing, the data set included 964 predictors, namely 314 eukaryotic, 627 prokaryotic, 14 physicochemical, 7 biological ones, substation and macrophyte type.Samples collected from the same station were grouped, to avoid data leakage between test and train data.
Both RF models consisted of 2000 trees with 30 randomly sampled variables as candidates at each split.The model was trained and evaluated using 10-fold cross-validation on a dataset of 52 observations, each representing the average of three biological replicates.Performance metrics used for evaluation included the mean absolute error (MAE), average absolute difference between predicted and actual values, root mean squared error (RMSE), a measure of the differences between the values predicted by the model and the actual values, and coefficient of determination (R 2 ) of the 10fold cross-validation.The top 10 predictors of both models are discussed.

Eutrophication index
A eutrophication index, defined as the organic matter availability in an ecosystem 89 , was derived by performing a principal component analysis (PCA), using the "FactoMineR (V 2.8)" R package 90  .All organic nutrients strongly aligned with principal component one, explaining 63%, which was chosen as the eutrophication index (Supplementary Fig. 9).

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 1 |
Fig. 1 | Map of the Baltic Sea showing sampling stations and salinity levels.A schematic zoomed-in view of one station (BV-12) is provided to illustrate the three substations sampled at every station: A, B, and C. Station A is located within a macrophyte meadow, while station B is located 15 m, and station C is located 100 m from the meadow edge.
41.76 μM, DOC between 180.10 and 681.30μM, dissolved organic nitrogen (DON) between 8.30 and 30.60 μM and chlorophyll-a (chl-a) between 0.63 and 21.16 μg L −1 , with the highest concentrations found along the German coast.The eutrophication index was highest on the German coast (BV-11) and lowest on the Danish coast (BV-26).

Fig. 3 |
Fig. 3 | Relationship between temperature, eutrophication index, and pelagic V. vulnificus (vvhA copies mL −1 ).The x-axis consists of 1-degree temperature bins for the bars and is continuous for the points, the y-axis is the log10 transformed vvhA gene copy numbers.Colours represent the eutrophication status per sample and are averaged per bin for the bars.The point shape indicates heterotrophic bacterial cell abundance, grouped into five bins, ranging from 3.4 × 10 6 to 2.7 × 10 7 cells mL −1 .Error bars represent standard deviation.

Fig. 4 |
Fig. 4 | Heatmap of the most important predictors for the vvhA gene copy number, derived V. vulnificus abundance, and the eutrophication index.Predictors are ordered from high (left) to low (right) importance for the respective models.Observations are ordered according to decreasing vvhA gene copy numbers.Values of each predictor are scaled from 0 to 1 for color representation in the heatmap, while the original values are printed in each cell.