Bacterioplankton dynamics driven by interannual and spatial variation in diatom and dinoflagellate spring bloom communities in the Baltic Sea

In parts of the Baltic Sea, the phytoplankton spring bloom communities, commonly dominated by diatoms, are shifting toward the co‐occurrence of diatoms and dinoflagellates. Although phytoplankton are known to shape the composition and function of associated bacterioplankton communities, the potential bacterial responses to such a decrease of diatoms are unknown. Here we explored the changes in bacterial communities and heterotrophic production during the spring bloom in four consecutive spring blooms across several sub‐basins of the Baltic Sea and related them to changes in environmental variables and in phytoplankton community structure. The taxonomic structure of bacterioplankton assemblages was partially explained by salinity and temperature but also linked to the phytoplankton community. Higher carbon biomass of the diatoms Achnanthes taeniata, Skeletonema marinoi, Thalassiosira levanderi, and Chaetoceros spp. was associated with more diverse bacterial communities dominated by copiotrophic bacteria (Flavobacteriia, Gammaproteobacteria, and Betaproteobacteria) and higher bacterial production. During dinoflagellate dominance, bacterial production was low and bacterial communities were dominated by Alphaproteobacteria, mainly SAR11. Our results suggest that increases in dinoflagellate abundance during the spring bloom will largely affect the structuring and functioning of the associated bacterial communities. This could decrease pelagic remineralization of organic matter and possibly affect the bacterial grazers communities.

During the last decades, shifting phytoplankton spring bloom communities from diatom-dominated blooms toward higher abundances of dinoflagellates have been reported in some sub-basins of the Baltic Sea (Klais et al. 2011). This tendency seems to be due to certain climate change-related effects, such as increases in water temperature and in wind speed and decreases in the thickness of the ice-cover, which appear to favor the excystment of the coldwater dinoflagellates resting cysts that may dominate until the early summer (BACC II Author Team 2015; Legrand et al. 2015). Diatoms and dinoflagellates have different ecological traits, in terms of growth rate (Spilling and Markager 2008), quality and quantity of dissolved organic matter (DOM) release (Myklestad 1995), and sedimentation patterns (Heiskanen 1998). Heterotrophic bacteria are the main consumers of the DOM produced by phytoplankton in the Baltic Sea Lindh et al. 2015), and both diatoms and dinoflagellates have been associated to different bacterial communities (Pinhassi et al. 2004;Camarena-Gómez et al. 2018). Thus, understanding how the predicted changes in the nature of the Baltic Sea spring bloom may impact bacterial community composition and activity is essential to estimate any potential biogeochemical consequences of this global warming driven process.
The nature and the timing of the spring bloom in the Baltic Sea vary between the different sub-basins (Kahru and Nômmann 1990). The bloom typically starts with light promoting net primary production (PP) in the southernmost Baltic Sea in February/March, reaching the Gulf of Finland (GoF) in April, and the Gulf of Bothnia in May. The bloom reaches the peak phase when inorganic nutrients have been depleted; N-limitation prevails in most of the Baltic Sea except for the P-limited Bay of Bothnia (Andersson et al. 1996;Tamminen and Andersen 2007). The subsequent decline phase is characterized by the rapid sinking of the phytoplankton cells (Heiskanen 1998). Therefore, understanding the bacterial responses to changes in phytoplankton communities requires considering the large environmental and spatiotemporal heterogeneity between sub-basins, that is, the timing of the bloom and the pronounced salinity (2-20) and temperature (0-20 C) gradients (Hagström et al. 2001).
Several studies in the Baltic Sea have shown that the taxonomic composition of bacterial communities varies along the salinity gradient with Alphaproteobacteria and Gammaproteobacteria dominating in more saline waters, whereas Actinobacteria, Betaproteobacteria, and Planctomycetes are favored by lower salinity (Rieck et al. 2015;Bunse et al. 2016). The strong seasonality in the Baltic Sea also shapes bacterial community structure and dynamics: Alphaproteobacteria are present throughout the year Rieck et al. 2015), Bacteroidetes (mainly class Flavobacteriia), and some taxa from Beta-and Gammaproteobacteria classes are associated with phytoplankton spring blooms (Hugerth et al. 2015;Bunse et al. 2016), Verrucomicrobia, Actinobacteria, and Planctomycetes appear in summer/autumn . Although both the physicochemical environment and biotic interactions affect bacterial communities, most of these studies have been restricted to specific sub-basins or did not capture the full development of the spring bloom, thus preventing a comprehensive understanding of bacterial responses at the whole Baltic Sea scale. A recent experimental study in the Baltic Sea showed that diatomdominated treatments promoted the growth of copiotrophic bacterial groups such as Flavobacteriia or Gammaproteobacteria as well as increases in bacterial heterotrophic production compared to dinoflagellate-dominated treatments (Camarena-Gómez et al. 2018). However, so far, no study has addressed this issue in natural conditions or across the naturally occurring environmental gradients in the Baltic Sea.
The aim of this study was to investigate how differences in the taxonomic composition of phytoplankton spring bloom communities (diatom-or dinoflagellate-dominated) affect the community structure and dynamics of the associated bacterioplankton across different areas of the Baltic Sea. For that, we collected samples in four consecutive years (2013)(2014)(2015)(2016) during different phases of the spring bloom (Growth, Peak, Decline, and Postbloom), from the northern Gulf of Bothnia to the southern Baltic Proper (BP). We explored changes in bacterioplankton communities and their heterotrophic production in relation to environmental variables (abiotic and phytoplankton communities). Based on our previous experimental studies (Camarena-Gómez et al. 2018), we hypothesized that dinoflagellate dominated communities will lead to less productive bacterial communities with reduced number of copiotrophic bacteria compared with diatom-dominated communities.

Study area, sampling design, and environmental variables
The water samples were collected during four cruises (April-May; 2013-2016) onboard R/V Aranda. In total, 127 stations were sampled, and the area spanned from the southern BP (55 22 0 N) to the northern Bay of Bothnia (BoB, 65 53 0 N), also covering GoF, Åland Sea (ÅS), Archipelago Sea (ArS), and Bothnian Sea (BS) sub-basins (Table 1; Fig. 1A). Four samples were collected in the Kvarken region (bordering BoB and BS) and these were added to the BoB stations. The number of samples differed by sub-basin (Table 1). At all stations, the water was collected from 3 m depth either using an oceanographic rosette with Niskin bottles of 5 L (n = 122) or from the flowthrough system on board (n = 5) in 2013.
Salinity, temperature, and depth were recorded in situ. Inorganic nutrients (NO 2 + NO 3 -N, NH 4 -N, PO 4 -P, and dissolved silicate (DSi)) were analyzed by standard colorimetric methods (Grasshoff et al. 1983). Samples for chlorophyll a (Chl a) determination (100-200 mL) were filtered onto glass fiber filters (GF/F, pore size: 0.7 μm, Whatman) in duplicates. Chl a was extracted in 10 mL of ethanol (96% v/v) (Jespersen and Christoffersen 1987) and stored at −20 C until further determination with a fluorescence spectrophotometer (Agilent Cary Eclipse). Dissolved organic carbon (DOC) and dissolved organic nitrogen (DON) were analyzed from 0.2 μm filtered samples using a Shimadzu TOC-V CPH analyzer and following the temperature catalytic oxidation method (Benner et al. 1993).

Plankton community composition, biomass, and productivity
For the determination of nano-and microplankton, including phytoplankton, water samples (200 mL) were fixed with acid Lugol's solution on board and stored in darkness at 4 C until enumeration. Depending on the Chl a concentration, subsamples of 10-50 mL were used to enumerate the plankton species abundance and to measure cells biovolume with an inverted light microscope (Leitz DM IRB, Leica). The abundance was converted to carbon biomass as described by Lipsewers and Spilling (2018). PP, defined as gross production during 2 h incubations, was determined by measuring the incorporation of 14 C-labeled sodium bicarbonate with a scintillation counter (Wallac 1414, Perkin Elmer) according to Steemann-Nielsen (1952). Total dissolved inorganic carbon (DIC) was determined using a high-temperature combustion IR carbon analyzer. The PP was calculated from the incorporated 14 C after incubation knowing the total amount of added isotope and the DIC concentration according to Gargas (1975). A more detailed description of the method can be found in Spilling et al. (2019).
The DOC released by phytoplankton cells was measured by filtering the total community onto 0.2 μm pore size polycarbonate filters after 24 h incubation with NaH 14 CO 3 in a climate control room at surface water temperature (2-6 C). The percent extracellular release of 14 C (PER) was calculated based on the dissolved organic fraction (DO 14 C < 0.2 μm) of the total 14 C fixation after incubation. Since substantial heterotrophic consumption of DOC (on average, 30-50%) might take  Camarena-Gómez et al.
Diatoms, dinoflagellates, and bacterial dynamics place during the incubation period (Morán and Estrada 2002), the results are regarded as net DOC production rates. Different phases of the spring bloom, named as Growth (Gr), Peak (Pe), Decline (De), and Postbloom (PB) were defined for each station and sampling year, as described in detail in Spilling et al. (2019). In brief, the phases were defined based on the concentration of inorganic nutrients (PO 4 and NO 3 + NO 2 ) and Chl a (Supporting Information Table S1). The NO 3 + NO 2 was used for the N-limited sub-basins GoF, BP, ArS, ÅS, and BS; and PO 4 3− for the P-limited BoB sub-basin.
High Chl a was defined as within 20% of the average peak concentration during spring, which is different in the different sub-basins (Supporting Information Table S1).

Bacterial heterotrophic production
Bacterial production estimated as thymidine (BPT) and leucine (BPL) uptake were measured in triplicates (1 mL sample) by using the dual labeling with [metilhyl-3 H]-thymidine and [ 14 C (U)]-leucine at final saturating concentrations of 20 nmol L −1 and 166 nmol L −1 , respectively. One sample was fixed with formaldehyde (final concentration 1.85%) and was used as a blank. The samples were incubated at in situ temperature during 2 h in darkness. The incubation was stopped by the addition of formaldehyde. The samples were processed using the cold trichloroacetic acid extraction method (Fuhrman and Azam 1982) and the centrifugation method (Smith and Azam 1992). The thymidine incorporation was converted to carbon production (μg C L −1 h −1 ), using a cell conversion factor of 1.4 × 10 9 cells nmol −1 (HELCOM 2008) and a carbon conversion factor using the formula of 0.12 pg C × (μm 3 cell −1 ) 0.7 (Norland 1993). We assumed a bacterial cell volume of 0.06 μm 3 , which corresponds to the mean cell volume of bacteria previously measured during spring blooms in the Baltic Sea (Kuparinen and Kuosa 1993). The leucine incorporation was converted to carbon production using a factor of 1.5 kg C mol −1 (Simon and Azam 1989).

Bacterial abundance
Bacterial abundances (BA) were determined in duplicates (1.2 mL sample) by flow cytometry according to Gasol and del Giorgio (2000). The samples were fixed with paraformaldehyde (final concentration 1%) and incubated 10 min in darkness and stored at −80 C until further analysis. The counts were done by staining the cells with SYBER Green I (Sigma-Aldrich) and counting them with a Partec-CUBE flow cytometer as described in Camarena-Gómez et al. (2018).

Bacterial community composition
Samples for bacterial community analysis (500 mL) were gravity-filtered onto sterile 0.2 μm pore-size Whatman cellulose ester filters, which were stored at −80 C until further analysis. DNA was extracted using the Power Soil DNA isolation kit (Mo Bio Laboratories). We amplified the 16S ribosomal RNA (rRNA) gene, V1-V3 hypervariable region. We used a two-step polymerase chain reaction with the universal primers F27 (AGAGTTTGATC[ACTG]TGGCTCAG) and R519 (GTATTACCGCGGCTGCTG) (Lane 1991). Amplicons were subsequently paired-end sequenced on an Illumina MiSeq platform, using multiplexing at the Institute of Biotechnology, University of Helsinki, Finland. The processing of the sequences was done according to the pipeline described by Logares (2017). In brief, primers were removed with Cutadapt (Martin 2011). The paired-end reads were merged with PEAR (Zhang et al. 2014). Quality filtering (> 400 bp, maximum expected errors = 1), chimera checking, and operational taxonomic unit (OTU) clustering (99% similarity) were done with the UPARSE pipeline (Edgar 2013). In total, 4.22 million sequences passed the quality filtering. After the removal of singletons as well as chloroplasts and mitochondria, based on the classification with SILVA v123, 2247 OTUs and 2.25 million reads were retained. Samples with less than 7000 reads (four samples) were removed. The remaining samples (122) were normalized by rarefaction using rrarefy in R (R Development Core Team 2008) to the lowest number of reads (n = 7024 reads). Downstream analyses were carried out using the rarefied data, including 122 samples and 2128 OTUs.
The DNA sequencing data has been submitted to ENA (European Nucleotide Archive) under the accession number PRJEB37659.

Statistical analyses
Bacterial community turnover was analyzed comparing the similarity between samples. A Bray-Curtis dissimilarity matrix was generated and used for nonmetric multidimensional scaling (NMDS) analyses using the metaMDS function of the R package vegan. Taxonomic richness was estimated with the Chao 1 index and taxonomic diversity with Shannon index. Significant differences in bacterial heterotrophic production, bacterial abundance, and Shannon index between the different sampling years and phytoplankton bloom phases were assessed applying nonparametric tests: Kruskal-Wallis for comparing the samples and Wilcoxon rank-sum test as a post hoc test.
Heterogeneity in the measured local environmental variables (biotic and abiotic) was compared against community turnover (Bray-Curtis dissimilarity) by means of Mantel linear correlations (mantel.rtest function, R package ade). The distance matrices of the environmental variables, z-score transformed, were constructed by computing the Euclidean distances of either all abiotic variables (salinity, temperature, dissolved oxygen, inorganic nutrients, dissolved inorganic and organic carbon), or single abiotic variables (salinity and temperature), or by computing the Bray-Curtis distances of the biotic variables: whole phytoplankton assemblages, or only the diatom communities, or only the dinoflagellate communities.
The effect of the environmental variables (explanatory variables) on the structure of the studied bacterioplankton communities was further assessed by applying a redundancy analysis (RDA, maximum 200 permutations) with forward selection (R package vegan) to obtain the abiotic or biotic variables that significantly explain the bacterial community variance (response variable). The variables selected for the RDA analysis were all the environmental variables included in the Mantel tests. Prior to the analysis, the variables were normalized by calculating the z.scores. The RDA forward selection analysis included an ANOVA test in each step (permutation test for RDA under reduced model, permutations: free, number of permutations: 999) to identify the significant environmental variables (p < 0.05). Collinearity between the explanatory variables was tested before and after the analysis to exclude the highly correlated variables by using the variance inflation factors (VIFs); only variables with VIF < 10 were retained for further analysis. The significant explanatory variables were used to construct the final model and the results were presented in an RDA plot. A variation partitioning analysis was performed to explore which portion of the variation in bacterial community structure could be explained by either the biotic (plankton communities) or the abiotic variables. This was done considering all the data together (N = 122 samples) or each year individually (2013-2016) and constructing two explanatory matrices representing the abiotic component and the planktonic (biotic) component, which included the variables that had been previously selected with the forward selection analysis. Spearman correlations were used to compare the relative abundances of specific bacterial groups with bacterial production, bacterial abundance, and variations in environmental variables (abiotic or biotic) using the R package ggpubr. The results of the correlation analysis were presented in a Heatmap by using the R package ComplexHeatmap. All figures and analyses were performed in R (R Development Core Team 2008).

Environmental variation and plankton bloom dynamics
The environmental conditions varied largely within and between the sampling years due to the different sub-basins and phytoplankton bloom phases covered (Table 1; Fig. 1A, Spilling et al. 2019). Salinity ranged from 8 in the southernmost station of the BP to 2 in the northernmost part of the BoB sub-basin (Table 1). Temperature was generally low (< 2 C) in 2013, which had high inorganic nutrient concentrations and all the sub-basins in this year were mainly in Growth and Peak bloom phases (Table 1). The stations where the temperature was > 4 C had low inorganic nutrient concentrations and were in Decline and Postbloom phase (e.g., in 2014 and 2015). Exceptionally, BoB sub-basin had high nitrite + nitrate (NO 2 − + NO 3 − ) concentrations (4-7.5 μmol L −1 ), similar to the concentrations found during Camarena-Gómez et al. Diatoms, dinoflagellates, and bacterial dynamics the Growth phase of the bloom, and the phosphate (PO 4 ) was almost depleted after the bloom. The lowest DSi concentration (2.27 μmol L −1 ) was measured in the GoF in 2014, while the maximum DSi concentration ( 55 μmol L −1 ) was detected in the Kemi river plume in the northernmost BoB (Table 1) (Table 2) and coinciding with the stations that presented high phytoplankton carbon biomass (Fig. 1B). The highest PP was also measured in the GoF and largely in 2013 (61 μg L −1 h −1 , Table 2). The percentage extracellular release (PER) was opposite to the PP, increased in the Decline and Postbloom phases and was higher in the BP compared to the other sub-basins (Table 2).
Overall, the planktonic communities, including mostly phytoplankton and nano-microzooplankton, were largely dominated by diatoms, dinoflagellates, the mixotrophic Mesodinium rubrum, and heterotrophic ciliates, contributing > 80% to the total carbon biomass (Fig. 1B, Supporting Information Fig. S1). Diatom biomass decreased from the Growth bloom phase (38% total carbon) to the Postbloom phase (< 5% total carbon), whereas dinoflagellate biomass was similar in all bloom phases ( 33% total carbon) and heterotrophic organisms increased toward the postbloom phase (Supporting Information Fig. S1A). When comparing the different sampling years, a much higher total biomass (> 300 μg C L −1 ) was observed in 2013 than in the other 3 yr. In addition, the biomass in 70% of the communities in 2013 was largely dominated by diatoms (> 140 μg C L −1 ) and dinoflagellates (110 μg C L −1 ; Fig. 1B), coinciding with the early spring bloom (mainly Growth and Peak phases). In the other years, dinoflagellates had higher biomass than diatoms (Fig. 1B). For instance, the dinoflagellate biomass in 2014 and 2016 were > 90 μg C L −1 in 40% of the stations, which also displayed more heterotrophic organisms such as ciliates and the heterotrophic dinoflagellate group Gymnodiniales, named also as Dinoflagellates HT, coinciding with a later bloom phase (  Fig. S1). The silicoflagellate Ebria tripartita (fam. Ebriidae) represented a small fraction of the biomass but contributed more to the total biomass in 2013 than in the other years (Fig. 1B) and mainly in the sub-basins with a high abundance of diatoms (GoF and ArS; Supporting Information Fig. S1B). The contribution of other autotrophic plankton organisms, such as Cryptophyceae, Chrysophyceae, Euglenonophyceae, Prasinophyceae, Chlorophyceae, and Nostocophyceae, to the total biomass was generally low (< 1%; Supporting Information Fig. S1).
Within diatoms, Achnanthes taeniata, Skeletonema marinoi, Thalassiosira levanderi, and Chaetoceros spp. were the dominant species (Fig. 1B). In the BS sub-basin, only the Peak phase was sampled (Table 1) and the single diatom Thalassiosira baltica dominated the community (Fig. 1B). Autotrophic dinoflagellates were dominated by the dinoflagellate complex (a group Table 2. Biological variables measured along the four sampling transects (2013-2016): Chl a, PP, percentage of extracellular release (PER), bacterial heterotrophic production, measured as leucine (BPL) and thymidine (BPT) incorporations rates, bacterial abundance, and Shannon index of the bacterial community. The table shows the range between the minimum and maximum values (indicated in bold), recorded for each variable and measured in each sub-basin. The sub-basin abbreviations are defined in Table 1.

Bacterioplankton community structure and its environmental drivers
The taxonomic composition of the studied bacterial communities varied largely across the stations sampled ( Fig. 2A). Overall, Alphaproteobacteria, mainly SAR11, dominated in all the sampling years, being present in all the sub-basins and bloom phases ( Fig. 2A; Supporting Information Fig. S2). The classes belonging to the Phylum Bacteroidetes (Flavobacteriia, Sphingobacteriia and Cytophagia) had higher relative abundance during 2013 and 2014 than in the other years. The class Flavobacteriia dominated within Bacteroidetes and comprised on average 25% of the total sequences ( Fig. 2A). Within Flavobacteriia, the genus Flavobacterium contributed 8-12% to the total reads in 2013. Other Flavobacteriia genera such as Fluviicola, Polaribacter, and Owenweeksia also contributed to the community. Gammaproteobacteria dominated by the genus Crenothrix contributed 10-15% of the communities sampled in 2013 ( Fig. 2A) and, in general, in the GoF and ÅS sub-basins (Supporting Information Fig. S2B) but decreased, even more acutely than Flavobacteriia, as the bloom progressed (Supporting Information Fig. S2A). Betaproteobacteria, largely dominated by the BAL58 marine group, showed higher relative abundance (5-7%) in 2013 and 2014 than in 2015 and 2016 ( Fig. 2A). The classes Actinobacteria (hgcl clade) and Acidimicrobiia (CL500-29 marine group) had higher relative abundance in 2014 compared to the other years, reaching values of 12.73% and 9%, respectively ( Fig. 2A). Actinobacteria were found mainly to be associated with the late bloom phases (Supporting Information Fig. S2A). The phylum Phycisphaerae (CL500-3) showed generally low abundances in all the bloom phases but reached its highest values (11.03%) in the BoB ( Fig. 2A; Supporting Information Fig. S2B), sampled only in 2015.
Based on their overall taxonomic structure, the bacterial communities clustered into two groups differing largely in taxonomic richness (Chao 1; Fig. 2B) and diversity (Shannon index; Table 2). The group with significantly higher taxonomic diversity comprised all communities sampled during 2013 and 2014 (Wilcoxon; p < 0.0001; Table 2), which were characterized by higher relative abundances of Flavobacteriia, Sphingobacteriia, Cytophagia, Betaand Gammaproteobacteria, and Actinobacteria. The low taxonomic richness group, sampled in 2015-2016, was clearly  Table 1 and in Fig. 1. The Y-axis indicated the relative abundance of the bacterial taxa. Only the bacterial groups that contributed more than 0.25% of the total sequences are indicated in the legend and the rest are grouped as "other bacteria." The classification was performed at the class or genus level in most cases, but the class Alphaproteobacteria was split into the main orders Rhodobacterales and SAR11. The lines separate the different sampling years. (B) NMDS of bacterial community structure based on Bray-Curtis dissimilarity (N = 122, stress = 0.134). The color indicates the sampling years and the dot size is proportional to the OTU richness (Chao 1) in each community. The vectors indicate the main bacterial groups found in the 4 yr, corresponding to those shown in panel (A). dominated by Alphaproteobacteria (SAR11 and Rhodobacteriaceae) and Cyanobacteria (Fig. 2B). The taxonomic diversity was also significantly higher during the Growth bloom phase (Wilcoxon; p < 0.001) compared with the other bloom phases, whereas no significant differences were found between the other bloom phases (Wilcoxon; p > 0.05; Supporting Information Fig. S3A). Differences in bacterial community composition were significantly correlated with differences in environmental variables, mainly with salinity (Mantel R = 0.52, p < 0.001; Fig. 3A) and temperature gradients (Mantel R = 0.37; p < 0.001), as well as with differences in the phytoplankton community composition (Mantel R = 0.51; p < 0.001). The high correlation found between differences in phytoplankton and bacterial community composition seemed to be mostly due to changes in diatom communities, as suggested by the higher Mantel R found when comparing with diatoms (Mantel R = 0.492, p < 0.001), than with dinoflagellates (Mantel R = 0.292, p < 0.001; Fig. 3A). Interesting, the correlation between changes in the bacterial assemblages and those in Chl a concentration was lower tions between differences in environmental (abiotic or plankton) variables and bacterial Bray-Curtis community structure dissimilarity matrices, pooling all samples together (N = 122). The dissimilarity matrices compared include taxonomic dissimilarities of the total phytoplankton communities (Phyto), of the diatom communities (diatom), and the dinoflagellate communities (Dinos), Chl a, all physical-chemical variables (salinity, temperature, dissolved oxygen, inorganic nutrients, dissolved inorganic and organic carbon), salinity, and temperature (tem). **p < 0.001, *p < 0.01. (B) RDA of the bacterial community including the explanatory significant environmental variables selected with the RDA forward selection: Salinity and temperature (tem), the biomass of A. taeniata, S. marinoi, T. levanderi, Chaetoceros spp., T. baltica, Melosira arctica, other diatoms, Cryptophyceae, Gymnodiniales and ciliates. The r 2 value (R adj ) and the p value (ANOVA-like permutation test) represent the most parsimonious model. the color indicates the sampling years (as in Fig. 1) and the dot size is proportional to the OTU richness (Chao 1) in each community. (C) Results of partial regression analysis including the variables selected from the RDA forward selection, partitioning the variation in bacterial community composition considering all the samples together (all) and each sampling year. Four different components are shown: Abiotic variation (orange) independent of plankton community; variation due to changes in planktonic communities (green) independent of any abiotic factors; variation attributable to a combination of plankton and abiotic factors (blue); unexplained variation (gray).
(Mantel R = 0.17, p < 0.01) than the correlation found between bacterial and phytoplankton dissimilarity matrices. A more detailed analysis on the effect of the environmental variables on the bacterial community composition revealed that, together with the abiotic variables (temperature and salinity), the abundances of several phytoplankton groups significantly explained the observed variation in bacterial communities (Fig. 3B). For example, lower temperatures and the presence of several diatom species such as A. taeniata, S. marinoi, T. levanderi, and Chaetoceros spp. shaped the bacterial communities in 2013, whereas the salinity gradient seemed to differentiate communities from 2015 and 2016, together with the biomass of the diatom T. baltica. Higher abundances of organisms, such as Cryptophyceae and some heterotrophs organisms like Gymnodiniales, ciliates and HNF, were associated with the bacterial communities sampled during 2016 (Fig. 3B).
These results were supported by the variation partitioning analysis demonstrating that both, the abiotic (salinity and temperature) and biotic variables (plankton community) influenced the variance in taxonomic composition when pooling all the data together, explaining 19% and 24% of the variance, respectively (Fig. 3C). However, when the different sampling years were considered separately, we observed a stronger effect of the plankton community than the abiotic variables on the bacterial community structure in all the sampling years except in 2014 (Fig. 3C).

Environmental drivers of major bacterioplankton groups
We then explored the individual relationships between different abiotic and biological factors and the abundances of the main bacterial groups. The relative abundances of Flavobacteriia, Cytophagia, Sphingobacteriia, Gamma-and Betaproteobacteria were positively correlated with bacterial heterotrophic production, with all the diatoms species, except M. artica, and also with the dinoflagellate complex (Fig. 4). Conversely, they were negatively correlated with the carbon biomass of Cryptophyceae, heterotrophic organisms (Gymnodiniales and ciliates), temperature, and salinity (Fig. 4). Only Flavobacteriia correlated positively with salinity. SAR11 and Cyanobacteria correlated positively with temperature and with bacterial abundance in the case of SAR11, whereas the dominance of certain diatoms or autotrophic dinoflagellates caused decreases in their relative abundances (Fig. 4). SAR11 also correlated negatively with salinity, whereas Rhodobacterales showed positive correlations. Phycisphaerae correlated positively with the biomass of the diatom T. baltica.
When exploring the relationship between the bacterial groups and all diatoms or dinoflagellates together, Gammaproteobacteria was the group showing the strongest positive correlation (R = 0.751, p = 0.0001) with the total biomass of diatoms (Fig. 5A), being also positive with dinoflagellates (Fig. 5B). In contrast, SAR11 was the group showing the strongest negative correlation with the total biomass of diatoms (R = −0.364, p = 0.0001; Fig. 5C) and dinoflagellates (R = −0.456, p = 0.0001; Fig. 5D). Betaproteobacteria, Cytophagia, and Sphingobacteria also correlated positively and greatly with diatoms compared to dinoflagellates (Supporting Information Table S2), and more specifically with the diatom groups S. marinoi, A. taeniata, Chaetoceros spp., and T. levanderi (Supporting Information Table S3).

Links between bacterial community structure and prokaryotic heterotrophic activity
Overall, the bacterial production rates measured with leucine incorporation (BPL) were significantly higher in 2013 and 2014 compared with 2015 and 2016 (Wilcoxon: p < 0.0001), similarly to the Shannon diversity values ( Table 2). The greatest value of BPL was measured in 2013 (0.62 μg C L −1 h −1 ; Table 2), following the tendency of the PP (Table 2). When bacterial heterotrophic production was measured as thymidine incorporation (BPT), the rates were significantly higher in 2014 (Wilcoxon: p < 0.0001) than in the other years, reaching values of 1.67 μg C L −1 h −1 (Table 2). No significant differences were found in BPL (Kruskal-Wallis; chi-squared = 1.1544, df = 3, p = 0.764) and BPT (Kruskal-Wallis; chisquared = 7.3596, df = 3, p = 0.06128) between bloom phases (Supporting Information Fig. S3B,C). In contrast, the highest bacterial abundances (BA) were recorded at the stations sampled in 2015 and the lowest in 2013 (Wilcoxon: p < 0.0001; Table 2). The bacterial abundance was also significantly lower during the Growth phase compared with the other bloom phases (Wilcoxon: p < 0.001; Supporting Information Fig. S3D), opposite to the pattern shown by the Shannon diversity index.
Interestingly, the differences in bacterial community composition found across our data set seemed to have an impact on community functioning. For example, most of the bacterial groups that were prevalent in the high diversity assemblages in 2013 and 2014, that is, Flavobacteriia, and Actinobacteria, showed strong positive correlations with BPL (R = 0.656, p = 0.0001; Fig. 6A) and BPT (R = 0.662, p = 0.0001; Fig. 6B), respectively, coinciding with the highest bacterial production rates (Table 2). Also, Cytophagia, Sphingobacteria, Beta-and Gammaproteobacteria showed positive correlations with bacterial production rates (Supporting Information  Table S2), whereas the abundance of SAR11 correlated negatively (BPL: R = −0.634, p = 0.0001; Fig. 6C; Supporting Information Table S2).

Discussion
Our results demonstrate that the structure of the studied bacterioplankton communities from the Baltic Sea, besides being affected by the interannual and spatial variation in salinity and temperature, is also strongly influenced by the nature of phytoplankton communities, and in particular, by the abundance of specific diatom species. As we captured much higher biomass of diatoms in 2013 than in the other years, the spatiotemporal heterogeneity covered by the sampling allowed us to explore how changes in the dominant phytoplankton groups, mainly diatoms and dinoflagellates, may be translated into changes in the structure and functioning of the associated bacterioplankton communities.

Spatiotemporal variability in environmental variables and plankton communities
Our large-scale sampling allowed covering highly heterogeneous phytoplankton communities spanning most bloom phases, ranging from those with a codominance of diatoms and dinoflagellates, to those largely dominated by dinoflagellates. During the early bloom development, particularly in the GoF, the water was cold with high inorganic nutrient concentrations, and the phytoplankton community was dominated by diatoms that typically occur during the cold-water season (Lignell 1990;Spilling 2007). The dominance of diatoms was also associated with high primary productivity suggesting an actively growing community. The increased availability of inorganic nutrients and light after the ice cover disappears are known to boost and shape the bloom dynamics in terms of Chl a concentration, primary productivity, and dominant phytoplankton species in the Baltic Sea (Tamminen and Andersen 2007). In the Bothnian Sea, which we sampled during the Peak growth phase, the phytoplankton community was dominated by the diatom species T. baltica along with high biomass of the mixotrophic ciliate M. rubrum, both of which seem to be common in this sub-basin (Andersson et al. 1996). The diatom-dominated communities co-occurred to a varying extent with several dinoflagellates, most likely G. corollarium and B. baltica (Sundström et al. 2009;Olli and Trunov 2010) in the GoF and the BP, similarly to earlier observations (Klais et al. 2011). The heterotrophic E. tripartita was most abundant during diatom dominance, likely because it grazed on diatom species such as Skeletonema sp. and Thalassiosira spp., as previously reported (Hargraves 2002).
During the declining phytoplankton bloom, mainly sampled in 2014 and 2016, the water was warmer, and dinoflagellates contributed more to the biomass than diatoms. This is the typical situation in the GoF with diatoms being more abundant in early spring and autotrophic dinoflagellates gradually becoming more dominant as the nutrients are depleted (Heiskanen 1998;Spilling et al. 2018). For instance, the stations sampled in 2014 had low phytoplankton biomass and low DSi concentration in the GoF, suggesting that a diatomdominated bloom had settled out of the surface layer before the time of sampling. In the BoB, the low phytoplankton biomass was limited by phosphate with excess NO 2 + NO 3 and DSi. This sub-basin is P-limited (Tamminen and Andersen 2007) and is more turbid due to terrestrial DOM inputs limiting light penetration and thus, resulting in light-limited PP .
The abundance of nano-and microzooplankton such as heterotrophic nanoflagellages, ciliates, and heterotrophic dinoflagellates increased toward the decline phase, pointing to a more heterotrophic community. A decrease in phytoplankton biomass and an increase in the proportion of dinoflagellates and heterotrophic organisms have been projected for a future, warmer Baltic Sea , in addition to the earlier onset of the spring bloom . Remarkably, our results reflect different temperature scenarios, as although roughly the same sub-basins were sampled in 2013 and 2016, the stations in 2016 were warmer and the bloom was clearly in a more advanced stage, in line with these predictions.

Environmental variables drive bacterioplankton dynamics in the Baltic Sea
As hypothesized, the changes in the phytoplankton community composition between the different bloom phases and subbasins were accompanied by pronounced shifts in the structure of the associated bacterial communities. This is in line with previous experimental (Pinhassi et al. 2004;Sarmento et al. 2013) and field studies (Teeling et al. 2012;Lindh et al. 2015) of phytoplankton blooms. The release of phytoplankton-derived DOM during these blooms can lead to changes in bacterial community dynamics (Buchan et al. 2014).
The studied bacterial communities clustered into two clearly different groups characterized by high and low taxonomic richness and diversity, which corresponded to the communities sampled during 2013-2014, and those sampled during 2015-2016, respectively. Among the environmental variables tested, salinity and temperature were the only abiotic variables significantly explaining the observed structural shifts in bacterial assemblages, but only the temperature gradient seemed to explain the observed partition between high and low richness communities. Salinity and temperature are known to be relevant factors in shaping the bacterial communities in the Baltic Sea (Hugerth et al. 2015;Herlemann et al. 2016). Thus, the predicted increases in water temperature and decrease in salinity due to an increase in water run-off along the Baltic Sea may lead to significant changes in the planktonic assemblages Legrand et al. 2015).
Besides abiotic factors, the biomass of specific diatom species such as A. taeniata, Chaetoceros spp., S. marinoi, and T. levanderi emerged as some of the most important environmental drivers explaining the segregation of high-and low-richness bacterial communities. More precisely, changes in diatom communities were more strongly related to taxonomic dissimilarities in bacterial assemblages than temperature variations, as shown by the Mantel tests. The important role of the diatom community on explaining the variation in the bacterial community composition was backed-up also by the variation partitioning analysis, showing that diatoms explained 24% of the variation in bacterial assemblages. It is thus possible that the diatom bloom captured in 2013 and the diatom bloom that had likely recently settled out from the surface layer in 2014, favored a higher diversity of associated bacterial communities with high abundances of groups such as Gamma-and Betaproteobacteria, Flavobacteriia, Sphingobacteriia, and Cytophagia. This was further supported by the strong and positive correlations between the different diatom groups and the abundances of these copiotrophic bacteria. These bacteria are known to be associated with the diatom blooms (Amin et al. 2012;Camarena-Gómez et al. 2018;Liu et al. 2019), in accordance with their preference for productive or bloom-like conditions. In contrast, the lower diatom biomass observed during 2015 and 2016 was associated with less diverse communities dominated by the SAR11 in 2015 and also by Rhodobacterales in 2016. This was supported by the strong negative correlations between the relative abundances of these bacterial groups and the total biomass of phytoplankton. The conditions within the less productive stations sampled in 2015 and 2016 likely promoted the dominance of SAR11, which typically dominates nutrient-poor environments (Alonso-Sáez et al. 2007). Other groups like Phycisphaerae, and in particular CL500-3, were negatively affected by salinity, in accordance with a previous study where Phycisphaerae was associated to the lower salinity of the Gulf of Bothnia and to a higher influence of allochthonous DOM (Rieck et al. 2015).
Within copiotrophic bacteria, the contribution of Gammaproteobacteria and Flavobacteria differed between sampling years and bloom phases. Flavobacteriia were always present with fluctuating relative abundance, whereas Gammaproteobacteria, dominated by the genus Crenothrix, peaked locally and decreased toward the collapse of the phytoplankton bloom in stations with complex diatom communities. The presence of Crenothrix was surprising in the surface waters since this taxon has been observed in anoxic bottoms in the Baltic Sea (Koskinen et al. 2011) and also in stratified lakes (Oswald et al. 2017). In addition, Flavobacteriia and Gammaproteobacteria have been reported to occupy different niches in terms of their metabolic capacities: Flavobacteriia are capable of degrading high-molecular weight compounds, such as proteins, chitin, and polysaccharides (Teeling et al. 2012), whereas Gammaproteobacteria are more specialized on breaking down low-molecular weight substrates during increased organic and inorganic substrate availability (Gómez-Consarnau et al. 2012). This may explain the occurrence of Crenothrix in our samples since this taxon is a methane-oxidizer bacteria specialized in the uptake of simple carbon compounds such as methane or acetate (Stoecker et al. 2006). This taxon may have benefited from the increase of the suitable substrates, either by fresh and labile DOM released or as a product from a metabolic cascade after the degradation process carried out by other bacterial groups.
We observed lower net PER measured during the cruise in 2013, coinciding with a dominance of diatoms and higher bacterial production rates, compared with the higher net PER detected during the 2016 cruise dominated by dinoflagellates and with lower bacterial production rates. This may suggest that diatoms release more labile DOC relative to particulate organic carbon (POC) production than dinoflagellates and that this diatom released-DOM was rapidly consumed by bacteria. This type of response agrees with previous studies conducted in the Baltic Sea showing that some diatoms species release more DOC than flagellates (Wolter 1982;Lignell 1990). However, PER might vary largely and be affected to some extent by bacterial consumption, but also by the physiology of the phytoplankton cells and/or by the different phases of the bloom (Myklestad 1995;Morán and Estrada 2002), making difficult to conclude whether lower PER means higher heterotrophic activity. Besides the stimulation of certain bacteria by the release of DOM, phytoplankton and prokaryotic organisms are capable of interacting in many ways and can establish species-specific relationships that can be beneficial or harmful (Amin et al. 2012). Indeed, we observed negative correlations between major bacterial groups and the abundance of specific phytoplankton groups. For instance, the production of certain aldehydes by Thalassiosira spp. (Wichard et al. 2005) negatively affect the protein production of Bacteroidetes (Balestra et al. 2011). This could explain the low relative abundance of Flavobacteriia observed in the presence of the diatom T. baltica in the BS. Antibacterial effects of phytoplankton are not well studied, but there are indications that some species inhibit bacterial growth (Amin et al. 2012).

Links between bacterial community composition and ecosystem functioning
The observed changes in bacterial community structure resulted in pronounced changes in the bacterial heterotrophic production of the studied assemblages. Most stations sampled in 2013 and 2014 showed much higher rates of bacterial production based on the incorporation of leucine (BPL) compared to the stations sampled in 2015 and 2016. The BPL rates were similar or, in some cases, higher than those reported by previous studies conducted at similar temperatures during the spring in the Baltic Sea Bunse et al. 2019). The variation in BPL between studies conducted at similar temperatures might be due to the different dominating phytoplankton groups that produce DOM of diverse lability, but also to differences in the frequency of sampling at the studied sub-basins. For instance, Bunse et al. (2019) found large BPL values during a spring dinoflagellate-dominated bloom in the BP, whereas the BPL values we observed during a dinoflagellate bloom in 2016 dominated by P. catenata in the same sub-basin were half of the values they reported. This emphasizes the importance of the phytoplankton community composition as a driver of the bacterial carbon cycling mediated by the presence of specific bacterial taxa.
The higher values of BPL in our study coincided with the high richness of the bacterial communities in 2013 and 2014 and with the high abundances of copiotrophic groups such as Betaproteobacteria, Gammaproteobacteria, and Flavobacteriia associated with diatom dominance. In addition, large peaks in bacterial production based on the incorporation of thymidine (BPT) were observed in 2014 indicating an actively dividing community. This was mostly related to the increases in the relative abundance of Actinobacteria, which is known to occur after phytoplankton blooms (Hugerth et al. 2015) and to correlate positively with thymidine incorporation (Pérez et al. 2010). Gammaproteobacteria and Bacteroidetes are quickly grazed by heterotrophic nanoflagellates due to their preference to graze on actively growing bacteria (Alonso-Sáez et al. 2007). This may explain the low bacterial abundance of the active community observed in 2013. In contrast, SAR11 correlated negatively with BPL and positively with bacterial abundance which agrees with the small cell size and slowgrowing pattern of this clade, typically occurring in oligotrophic ecosystems, and which seem capable to avoid grazing (Alonso-Sáez et al. 2007). These bacterial responses agree with those observed in a mesocosm study performed with water sampled in the Baltic Sea, where the growth of diatoms promoted higher bacterial production and higher dominance of Flavobacteriia, Gammaproteobacteria, and Cytophagia than the presence of dinoflagellates (Camarena-Gómez et al. 2018).
It is interesting to note that the highest diversity of the bacterial assemblages was not only found during the period dominated by copiotrophic bacteria (Gammaproteobacteria, Betaproteobacteria, and Flavobaceriia), during early bloom phases or in the presence of diatoms, but also during postbloom phases, dominated by Actinobacteria, and in both cases associated with higher bacterial production rates. This suggests the occurrence of a bacterioplankton succession that responds to the increase of diatom-derived DOM by the increase in the bacterial production rates, but also responds to a niche partitioning of the highly diverse bacterial community that exploits different resources. Other studies have also reported higher bacterial diversity during the spring bloom dominated by diatoms in the Baltic Sea and Southern Ocean Liu et al. 2019). However, there are also contradictory findings with, for example, lower diversity with an increase in Chl a during a Phaeocystis antarctica bloom (Richert et al. 2019) and also lower diversity with an increase in PP (Raes et al. 2011). These contrasting results when comparing species diversity and ecosystem functioning might be linked to the different environmental circumstances in each study. For instance, differences in physico-chemical conditions, different spatial-temporal scales, nature of the phytoplankton community, or bloom conditions and differences in the sampling effort may lead to specific diversity-functioning relationships that cannot be easily predicted.
The bacterial remineralization of elements is critical for marine ecosystem functioning. The primary producers benefit from recycled nutrients and bacterial fixed carbon can be channeled to higher trophic levels through the microbial loop. Structural changes to bacterial communities and metabolism could consequently have profound effects on the fate and cycling of carbon within planktonic food webs. Our results link microbial community structure and function in the Baltic Sea to ongoing shifts in the phytoplankton community structure (Lindh and Pinhassi 2018). The data suggest that specific diatoms release high-quality DOM that enhances the growth of copiotrophic bacteria with high production rates, in turn favoring the recycling of carbon through the microbial loop. In contrast, the predicted increase in dinoflagellate abundance associated with warmer winters and springs in the Baltic Sea, seems to shift the bacterial community towards more oligotrophic generalist and reduce the bacterial production rates. This could have direct consequences for the bacterial grazers communities, affecting the transfer of carbon to higher trophic levels, and for pelagic nutrient remineralization, processes that may in turn be directly affected by changes in other environmental variables.

Conclusions
We found pronounced differences in the bacterial community composition and functioning driven by gradients in salinity and temperature, but also by differences in phytoplankton community composition. The positive correlations found between specific diatom species and the highly diverse bacterial copiotrophic communities as well as high bacterial heterotrophic activity emphasize a major role of diatom-derived DOM for bacterial remineralization of organic carbon. As such, our results supported our hypothesis of less productive bacterial communities dominated by other than copiotrophic bacterial groups during dinoflagellate dominance. However, the bacterioplankton response was linked more strongly to changes in phytoplankton at the species than at the group level, as the largest shifts in bacterial community structure were explained by the presence of specific diatom species such as A. taeniata, S. marinoi, T. levanderi, and Chaetoceros spp. Changes in the carbon biomass of these species explained the variability in bacterioplankton communities to a larger extent than the steep gradients in salinity and temperature. This suggests that quantitative data at a fine phytoplankton taxonomic resolution are needed to understand the role of phytoplankton communities in shaping bacterial ecology. Finally, our results emphasize that the projected shifts in spring bloom dynamics could cause changes in bacterial mediated carbon fluxes in the Baltic Sea.