The model system used compares potential changes on microbial communities, which have undergone 50 years of increased and prolonged warming in the closest to natural possible conditions. For this, a nuclear power plant using Baltic Sea water for cooling and releasing the heated water into a semi-enclosed bay (only open to the open Baltic Sea) was chosen. The water from a nearby site is pumped through the system to cool down the reactor within a heat exchanger. During this process, the cooling water is not in direct contact with radioactive material and is during release, on average, 10°C warmer at the discharge site (Seidel et al., 2022). The heated water in the bay was, on average, 5.3–5.7°C warmer compared to the selected control bay with the largest differences in winter (December–February) and little to no difference in summer (June–September). The control bay was chosen based upon its proximity to the heated bay without being affected by the warmed water as it was not connected to the heated bay other than via ~1.5 km of open Baltic Sea. As described in Seidel et al. (2022), the two bays are natural systems, and environmental parameters may vary, such as a thin freshwater plume observed in the control bay surface waters during spring 2018. However, previous findings along with this study indicated temperature as a key driver in bacterial community selection.

This study was composed of a time series of sampling in September, December, April, and June performed between September 2017 and June 2018 at the two coastal bays near the city of Oskarshamn, Sweden (temperature affected bay, 57°25'09.7“N 16°40'19.3”E and control bay 57°25'58.7“N 16°41'17.2”E; Supplementary Table 1). The samples did not overlap with those from the previous Seidel et al. (2022) study except for the six sediment samples taken in June 2018. Eight sampling sites per bay were selected with each site sampled in triplicate (Figure 1, Supplementary Table 1). The temperature was continuously measured at three sites per bay from December 2017 until November 2019 with data loggers (HOBOware, Onset Computer Corporation, USA). From the end of August 2017, warm water discharge was reduced for a 2-month period from 55 to 0.3 m3/sec. A description of data loggers set up, placement, and an overview of temperature changes between bays is given in Seidel et al. (2022) and Figure 1.

Three replicate sediment cores were sampled with a Kajak gravity corer at each location to give 24 acrylic cores (inner diameter: 7 cm; length: 60 cm) for each bay and a total of 48 cores per sampling occasion as described in Seidel et al. (2022). Surface plus bottom water temperatures and bottom water oxygen concentrations were measured in situ (Multiline™ sensor, WTW™) during sampling. Additionally, overlying bottom water samples (BW; 50 ml) were taken from the sediment (SED) cores close to the SED surface and transferred into sterile tubes (Thermo Scientific™) for 16S rRNA gene amplicon sequencing (described below). An additional 15 ml of BW was filtered through a 0.7-μm Target2™ GMF Syringe Filter (Thermo Scientific™) into a pre-acid washed polypropylene tube (Thermo Scientific™) for chemical analyses. The BW was decanted from the core, and a 0–1-cm SED slice was collected in a sterile 50-ml polypropylene tube (Thermo Scientific™). The tube contents were well-mixed, and 15 ml was transferred to a pre-acid washed polypropylene tube (Thermo Scientific™) for subsequent pore water analysis. A further 2-ml reaction tube (Sarstedt, Inc.) was filled with SED for organic matter (OM) analysis. All samples were cooled during transport to the laboratory on the same day where they were frozen at either −80 or −20°C.

BW, pore water [wet sediment was centrifuged at 2,200 × g for 15 min and the water filtered through a 0.7- μM Target2™ GMF Syringe Filter (Thermo Scientific™)], and SED was analyzed. pH (pHenomenal, VWR pH electrode, VWR™) was measured for BW and pore water samples. Total Fe, PO43-, and SO42- concentrations for pore water and BW were measured based on spectrophotometric methods (Broman et al., 2018). Additionally, NO3- and NO2-concentrations for pore water were measured using Hach–Lange cuvette tests LCK339 and LCK341, respectively. Organic matter content (% wt) was measured using the loss on the ignition method as described in Broman et al. (2017b) except that the sample was dried before loss on ignition analysis for 6 days at 40°C instead of 3 days at 80°C.

A total of 50 ml of the BW samples was filtered through a 0.22-μM GV Durapore^® Membrane filter (Merck) and the filter retained for DNA extractions. BW and SED DNA extractions were performed using the DNeasy^® PowerWater and PowerSoil Extraction Kits (QIAGEN), respectively, according to the manufacturer's guidelines. DNA concentrations were measured using a Qubit^® 2.0 (Invitrogen™, Life Technologies Corporation). Amplifications of partial 16S rRNA gene fragments were conducted using the PCR primers 341f and 805r according to a modified PCR program by Hugerth et al. (2014), while library preparations were conducted according to Lindh et al. (2015) .

16S rRNA gene amplicon sequencing was performed at the Science for Life laboratory (SciLifeLab) in Stockholm, Sweden. Sequencing was done on the Illumina MiSeq platform with a 2 × 301 bp paired-end setup with an average of 213,555 raw sequences obtained (Supplementary Table 2). The raw data (with adapter removal and demultiplexing performed by the sequencing facility) were filtered, trimmed, denoised, merged, and chimeras removed with the DADA2 paired-end pipeline (v. 1.16) (Callahan et al., 2016) on the UPPMAX cluster (Uppsala Multidisciplinary Center for Advanced Computational Science). Trimming of the sequences was at 290 bp forward and 230 bp reverse, and left trimmed at 21 bp to ensure primer sequences were removed from the analysis. Maximum consistency of the error model was set to 30, and a minimum overlap at 10 bp with 0 mismatches was used for merging the forward and reverse reads. After removing chimeras using the default settings, the taxonomy was assigned against the Genome Taxonomy Data Base (GTDB; version 86) set for DADA2 (Parks et al., 2021). The average sequence output count was 75,180 reads (min., 366 and max., 289,850), with a total of 65,952 amplicon sequence variants (ASVs; Supplementary Table 2). The data were analyzed using R [version 4.0.4 (2021-02-15) (R Core Team, 2018)]. Packages used within the R environment are stated below at the specific steps they have been used. As the 16S rRNA gene sequencing data are compositional (Reimann et al., 2012) and follow Aitchison geometry due to the inherent closed structure of the data (Aitchison, 1986), the dataset was transformed to relative abundance per sample. A total of 37,916 unique ASVs in BW and 50,337 in SED were detected, of which 0.79% in BW and 1.77% in SED relative abundance were doubletons within the whole community, respectively. Alpha diversity analysis with removed doubletons resulted in a similar pattern and statistical results, and, therefore, the analysis was continued on the whole dataset (Supplementary Figure 1, Table 3).

Metabolic response prediction analysis with the PICRUSt2.0 tool (v. 2.4.0) (Douglas et al., 2020) was performed on counts based on significant (p < 0.05) differentially abundant ASVs from the abundance analysis described below. PICRUSt2.0 was used with default settings according to the tool guidelines available at: github.com/picrust/picrust2/wiki. There are certain key limitations to the prediction of functional profiles based on 16S rRNA gene amplicon sequences (github.com/picrust/picrust2/wiki/Key-Limitations), and, therefore, the results were limited to an overall comparison between bays. Analysis was performed on each dataset (BW and SED) separately to produce EC (enzyme commission number), KO (KEGG Orthology) identifiers, and pathway abundances per sample with KO identifiers used for further analyses.

Rarefaction curves were calculated with the vegan (Oksanen et al., 2019) (v. 2.5-6) package in R, suggesting that the majority of microbial diversity had been covered in both bays (Supplementary Figure 1). Furthermore, unknown Cyanobacterial sequences were run through additional annotation using the database Phytoref (2017-04-04), and sequences annotated as eukaryotic diatoms were removed. For further analysis, samples under 1,000 reads were removed, which included one sample in the BW (X2061; Supplementary Table 1) and two in SED sets (X51 and X11; Supplementary Table 1). For Alpha diversity analysis (Shannon's H, Shannon's H evenness, and Chao1 richness), the data were normalized by scaling with ranked subsampling (Beule and Karlovsky, 2020). To test the robustness of the analysis of Alpha diversity, the data were tested and visually compared without rarefying, rarefying the data to the smallest samples size (n = 2,356), as well as removing doubletons. All these tests resulted in similar patterns to the full dataset, showing strong validation of the results (Supplementary Figure 1). To test if the diversity indices were significantly different between the two bays over the different sampling months, ANOVA was performed. Due to the nature of the dataset, a linear mixed-effects model was used. Tested interactions showed months and bays were dependent upon each other and, therefore, were modeled as interaction and were treated as a fixed variable within the model. The different sampling sites, with replicates nested within each site, were treated as random variables. The variable sites was usually nested within bays as well but overfitted the model. The different models were tested and compared with the best fit chosen based on the Akaike information criterion (AIC). The model was created with the “lme()” function within the “nlme” package [v. 3.1-152; (Pinheiro et al., 2021)]. Pairwise comparison of the bays at each sampling month was performed using the “emmeans” package (v. 1.5.4) in R [Supplementary Table 3; (Lenth, 2021)] to test for significant differences between the bays at each sampling month.

Environmental variables of BW and SED (both n = 191) at the different sampling months between the different bays were shown within a radar plot from the ggplot2 package (v. 3.3.5) (Wickham, 2016) using the rescale function and the discrete scale option. A principal component analysis (PCA) of the environmental variables (n = 191 for both BW and SED) showed an association with season, bays, and sites for BW and SED (Supplementary Figure 2). This result indicated that the same model, which has been used for the diversity analysis, could also be used for testing statistical differences between the bays at the different sampling months within the environmental variables. Testing of the best fit using the AIC was performed and resulted in “bay” and “month” as fixed variables and modeled as interaction, while replicates were nested within a site as a random variable [“lme()” function, “nlme” package, v. 3.1-152, Supplementary Table 3]. For the higher resolution of the results of the environmental variables, boxplots, including all replicates and sites within each bay and each month (n = 24) plus the individual concentrations of each sample stated as points, were plotted with the ggplot2 package (Supplementary Figures 3, 4). To compare for significant differences between the surface and bottom water temperatures in each bay, a mixed-effects linear model and ANOVA were used. The temperature of surface and bottom waters was modeled with interaction to months and treated as fixed variables, while sampling site (n = 8) was treated as a random variable (Supplementary Table 3).

Statistical differences between both bays on bacterial communities were evaluated based on 999 permutations in the perMANOVA (permutational multivariate analysis of variance) based on the Bray-Curtis dissimilarities “adonis()” function within the “vegan” package (Oksanen et al., 2019; Supplementary Table 3). Within the perMANOVA, the variables bays and months were stated as interactions (bay × month), while the variable site was controlled by the “strata” option. Samples “X69,” “X78,” “X93,” “X1141,” “X1045,” “X2013,” “X2014,” and “X2015” (further details of the samples within Supplementary Table 1) were excluded from further distance-based redundancy analysis (db-RDA) of the bacterial communities with environmental variables as well as Pearson correlation. This was due to missing geochemical variables such that they did not increase the knowledge of interactions between the communities and environmental variables. Db-RDA (n = 187 BW, n = 185 SED) based on the relative abundance of ASVs and environmental parameters as explanatory variables was performed within R using the “vegan” package. ANOVA (n = 999 permutations) was used to test for statistically significant environmental variables. Results based on environmental variables (PCA in Supplementary Figure 2) and ordination plots of the microbial communities (Supplementary Figure 5) showed that seasons were the main driver of changes in bottom water communities, while surface sediment communities were less affected. Therefore, the permutations were performed within bays in BW to test for seasonal influence, while it was performed within months in SED to test for differences between bays. To test for microbial community dissimilarities between samples outside of a constrained ordination environment, the unconstrained ortholog to a db-RDA non-metric dimensional scaling (nMDS) was performed for BW and SED (Supplementary Figure 5).

SIMPER (similarity of percentages) analysis was performed to discriminate taxa (on an order level) responsible for the dissimilarities (based on Bray-Curtis dissimilarities) between the heated and control bay for BW and SED samples, using the “vegan” package (Supplementary Table 4).

Differential abundance analysis was performed on each dataset from the 16S rRNA gene ASVs within BW and SED samples (both n = 191). Taxa not identified in at least 50% of the total samples were filtered out (thus removing rare abundant ASVs), which resulted in n = 476 taxa retained for BW and n = 536 taxa for SED. A zero-inflated negative binominal model was used to account for a large number of zeros in the dataset. The analysis was carried out within the “DESeq2” package in R (version 1.30.1) based on counts to calculate possible statistical differences between ASVs of the two bays at the different sampling months (Supplementary Table 5; Love et al., 2014). To account for the nature of the dataset, the model used within the differential abundance analysis included bays, months, and sites as fixed factors (as DeSeq2 does not support random effects), while replicates were nested within sites. Significant (p < 0.05) differentially abundant ASVs over 1% from the previous step and environmental variables were used to calculate Pearson correlations (n = 191 BWand n = 163 SED). The counts were summarized on the order taxonomic level and clr (the centered log ratio) transformed using the packages “zComposition” (version 1.3.4) (Palarea-Albaladejo and Martín-Fernández, 2015) and “CoDaSeq” (version.99.6) (Gloor and Reid, 2016). Correlation was calculated using the “ggcormat” function from the “ggstatsplot” package (Patil, 2018) with Benjamin–Hochberg correction for p-values (Supplementary Table 6).

To test for potential differences between the two bays within bottom water and surface sediments, a differential abundance analysis on predicted abundances of the KO identifiers based on PICRUSt2.0 predictions was performed using the ALDEx2 (version 1.22.0) package in R (Gloor et al., 2016). The differential abundance analysis tool was used as it was suggested by the creator of the PICRUSt analysis tool. The analysis was performed for each dataset (BW and SED) on default settings with “bay” as condition setting. Only significant (p < 0.05) results from Welch's t-test and the Wilcoxon rank test were retained (Supplementary Table 8).

Both the bottom and surface water temperatures were significantly higher in the heated bay compared to the control bay [mixed linear model, ANOVA, F_{(1, 14)} = 331.8, p < 0.0001 surface water and F_{(1, 14)} = 227.7, p < 0.0001 BW; Supplementary Table 3]. This difference was, on average, 5.7 and 5.3°C higher in the BW and surface waters, respectively (Figure 1, Supplementary Table 1). Furthermore, temperature fluctuations in the control bay were greater due to warm water discharge, limiting heated bay winter lows. The water from the cooling outlet was, on average, 10°C warmer compared to the ambient water temperature except during the period of 18 August−22 October 2017 when the reactor was turned off and similar water temperatures were observed between the bays. However, microbial community and chemistry comparisons were not strongly influenced during this period (Supplementary Figures 2, 5), suggesting changes did not rapidly occur, and the reactor shutdown did not notably influence the data.

Salinity concentrations were significantly higher in the surface [mean + s.d., a 6.6 ± 0.2% heated bay vs. a 4.8 ± 1.4% control bay, a linear mixed model, ANOVA, F_{(1, 14)} = 43.4, p < 0.0001] and bottom waters [6.6 ± 0.2% vs. 6.1 ± 0.6%, F_{(1, 14)} = 8.6, p = 0.01] of the heated bay compared to the control bay (Table 1, Supplementary Table 3). Oxygen concentrations were similar in the two bays [9.9 ± 5 vs. 10.9 ± 4.9 mg/L, F_{(1, 14)} = 0.1, p = 0.712]. Similarly, there was no significant difference in sediment organic matter (%) between the two bays (32 ± 13.3 vs. 34.6 ± 13.3%, F_{(1, 14)} = 0.1, p = 0.81; Figure 2, Table 2, Supplementary Figure 4, Supplementary Table 3].

Bottom water NO3-/NO2- concentrations (Figure 2, Table 2) were significantly higher in the heated bay during September (a 2.02 ± 2.2-μM heated bay vs. a 67 ± 1.2-μM control bay, pairwise comparison, p = 0.005; Supplementary Table 3) and had significantly lower concentrations in April (0.97 ± 0.6 vs. 2.35 ± 1.3 μM, p = 0.044). In addition, significantly higher heated bay SO42- concentrations were measured during April (2.75 ± 0.7 vs. 2.3 ± 0.4 mM, p = 0.013) and December (3.1 ± 0.8 vs. 2.7 ± 0.6 mM, p = 0.028). Furthermore, significantly lower concentrations of heated-bay totally dissolved Fe were measured during April (2.9 ± 0.23 vs. 3.9 ± 1.3 μM, p = 0.0077; Figure 2, Table 2 plus Supplementary Figure 3, Supplementary Table 3). The SO42- concentrations in both bays ranged between 2 and 5 mM (Figure 2, Table 2, Supplementary Figure 3). BW PO43- was not significantly different between the bays [10.6 ± 2.5 vs. 11.1 ± 3.5 μM, ANOVA, F_{(1, 14)} = 0.21, p = 0.65].

A significant decrease in concentration was measured in heated-bay sediment of NO3-in December (20.5 ± 14.2 vs. 37.45 ± 14.8 μM, pairwise comparison, p = 0.0182), and SO42- concentrations [3.70 ± 1.2 vs. 3.89 ± 1.3 mM, ANOVA, F_{(1, 14)} = 23.3, p < 0.0001; Figure 2, plus Supplementary Figure 4, Supplementary Table 3]. The lower concentration of sediment NO3- in the heated bay indicated increased microbial activity, resulting in lower (but not significantly) total dissolved Fe concentrations (3.38 ± 1.5 vs. 5.03 ± 2.4 μM, p = 0.1493) and increased use of SO42- as an alternative electron acceptor in December (Figure 2, Table 2, plus Supplementary Figure 4, Supplementary Table 3). Heated-bay sediment pore water PO43- was lower (but not significantly) compared to the control bay [a 63.1 ± 74.7 μM heated bay, a 125.9 ± 104.9 μM control bay, ANOVA, F_{(1, 14)} = 1.8, p = 0.21; Figure 2, Supplementary Table 3] during the whole year and showed the highest difference in September and December (Table 2, Supplementary Figure 4).

In summary, the heated bay bottom waters showed overall higher temperatures with less fluctuations over the sampling months. Additionally, higher salinity concentrations were observed within the heated bay bottom waters, with similar oxygen concentrations compared to the control bay, due to continuous mixing of the water column. Furthermore, the heated bay showed decreased NO3-and SO42- concentrations during winter (December), indicating the use of alternative electron acceptors plus lower PO43- concentrations compared to the control bay.

The Alpha diversity and evenness indices in the BW communities showed a significantly higher Shannon's H diversity in the control bay compared to the heated bay [mean ± s.d., a 5.2 ± 0.9 heated bay vs. a 5.7 ± 0.6 control bay, ANOVA, F_{(1, 14)} = 7.3, p = 0.0173] and evenness [0.8 ± 0.1 vs.0.9 ± 0.1, F_{(1, 14)} = 16.7, p = 0.001; Figure 3, Supplementary Figure 1, Supplementary Table 3]. Overall, more ASVs were found within the control bay (Chao1, mean, 720 ± 402.7 vs. 861.5 ± 321.2; Supplementary Table 1). The opposite pattern was observed in the 0–1-cm SED with an insignificant increase in the heated bay Shannon's H diversity (6.3 ± 0.5 vs. 5.7 ± 0.8) compared to the control bay, significantly higher number of ASVs [761.1 ± 162.8 vs. 630.8 ± 193.9, ANOVA, F_{(1, 14)} = 4.5, p = 0.05], and more equal distribution (0.9 ± 0.04 vs. 0.8 ± 0.1; Figure 3, Supplementary Figure 1, Supplementary Tables 1, 3).

16S rRNA gene amplicon sequencing of the BW microbial communities showed that the most abundant phyla in both bays were Actinobacteria [mean of relative abundance of total reads (%) of all sampling occasions: 29.9 ± 7.9 heated bay vs. 28.4 ± 6.3% control bay], Proteobacteria (21.8 ± 4. vs. 23.2 ± 3.2%), Bacteriodota (19.8 ± 6.3 vs. 13.2 ± 3.4%), and Cyanobacteria (11.3 ± 4.9 vs. 13. ± 5.2%; Supplementary Figure 6). The most abundant SED phyla in both bays were similar as found in the BW with Proteobacteria (27.7 ± 4.1 vs. 25.9 ± 3.1%), Cyanobacteria (5.2 ± 3.8 vs. 11.1 ± 8.2%), Bacteriodota (13.5 ± 3.5% vs. 16.5 ± 4.4%), and Desulfobacteriodota (9.8 ± 3.7 vs. 11.2 ± 2.1%; Supplementary Figure 6). Similarity of percentages (SIMPER) analysis showed 56.1% (BW) and 53.3% (SED) of overall dissimilarity between communities of both bays on the order level. The most likely taxa responsible for the differences between the BW communities were, e.g., Flavobacteriales (8.1% average contribution of dissimilarity), Synechococcales A. (6.1%), Microtrichales (3.9%), Betaproteobacteriales (3.6%), and Nanopelagicales (3.4%) (Table 3, Supplementary Table 4). Those taxa make up 39.2 ± 11.4% within the heated (37.7 ± 12.1% of the community significantly differentially abundant ASVs) and 36.2 ± 7.8% (33.7 ± 8.1% significant ASVs) within the control bay of the relative abundance of the whole bacterial community dataset. SIMPER analysis of the SED communities on the order level showed that Cyanobacteriales were most likely one of the main contributors to differences with 3.85%, followed by Betaproteobacteriales (3.5%), Desulfobacteriales (3.3%), Anaerolineales (2.8%), and Flavobacteriales (2.7%; Table 3, Supplementary Table 4). These top five orders make up 16.9 ± 7.8% (9.2 ± 3.9% significant ASVs) and 21.9 ± 5.8% (13.8 ± 5.3% significant ASVs) of the relative abundance in the heated and control bays, respectively.

Significantly differentially abundant ASVs summarized and annotated on the order level (with >1% relative abundance per sample) within BW samples of the heated and control bays showed a strong seasonal influence on bacterial community composition (Figure 4, Supplementary Table 5). Between both analyzed bays, BW samples comprised a slightly higher proportion of ASVs that differed significantly in abundance compared to SED communities. On average, 80.8% of the ASVs were significantly and differently abundant in the heated bay, while 79.4% were significant in the control bay (Figure 4). Within the high variation of bacterial community composition, the order Nanopelagicales (mean, a 6.6 ± 2.4% heated bay vs. a 7.8 ± 2.6% control bay) were more differentially abundant in the control bay compared to Flavobacteriales (a 12.1 ± 5.8% heated bay vs. a 6.4 ± 2.9% control bay) that were higher in the heated bay (Figure 4, Supplementary Table 6).

In contrast, surface SED communities did not show strong seasonal influence on significant differentially abundant ASVs (mean summarized on the order level >1% relative abundance; Figure 4), and the number of non-significant differentially abundant ASVs was, on average, 23.3% higher within sediment communities (19.9% BW vs. 43.3% SED). The main abundant taxa were Desulfobacteriales (1.82 ± 1.00% and 3.10 ± 1.07%), Flavobacteriales (2.3 ±0.9 vs. 3.7 ± 2.1%), and Betaproteobacteriales (2.4 ± 1.5 vs. 4.7 ± 2.6%) that had higher relative abundances in the control vs. the heated bay (Figure 4). Betaproteobacteriales were significantly positively correlated with phosphate (Pearson correlation, p = 0.02, r = 0.26), sulfate (p = 0.002, r = 0.29), and organic matter (p = 0.0002, r = 0.35) while decreasing with increasing temperature (p = 0.04, r = 0.22; Figure 4, Supplementary Figure 7, Supplementary Table 6). Furthermore, Flavobacteriales were positively correlated with sulfate (p = 0.03, r = 0.21), while Desulfobacteriales positively correlated to phosphate (p < 0.001, r = 0.32), organic matter (p = 0.02, r = 0.23), and temperature (p = 0.003, r = 0.28; Figure 4, Supplementary Figure 7, Supplementary Table 6).

Bacterial communities in BW were clustered by the different sampling months. The fluctuations in geochemical parameters included temperature (permutation test, ANOVA, F = 14.12, p = 0.001) and temperature-related parameters, including O₂ (F = 17.75, p = 0.001), NO3-/NO2- (F = 5.06, p = 0.001), total Fe (F = 6.88, p = 0.001), and SO42- concentrations (3.39, p = 0.003; Supplementary Table 3), associated with the contrasting communities (Figure 5). Microorganisms in BW in both bays clustered together significantly on the first axis (permutation test, ANOVA, F = 27.09, p = 0.001) according to the sampling time (season) and on the second axis according to the different bays (F = 13.62, p = 0.001; Figure 5, Supplementary Table 3). This indicated that short-term changes, such as seasonal variation, were the dominant drivers for community composition compared to long-term changes of increased temperature and prolonged warming (Figure 5). The communities of both bays were positively related with high O₂ concentrations in April and showed lower O₂ during the rest of the year, indicating high-water circulation-replenished O₂ during spring (Figure 5). Potential influences on microbial community composition in the heated bay could also be the higher SO42-concentrations in June in the heated bay and lower NO3-/NO2- (in combination) in April plus high concentrations in September. Finally, temperature was positively related with the heated bay communities (Figure 5).

Seasonal fluctuations likely did not substantially affect SED communities as there was a strong clustering on the first axis due to the different bays (permutation test, ANOVA, F = 9.51, p = 0.001; Figure 5, Supplementary Table 3), with no strong seasonal distribution between the bays, indicating that long-term changes, such as prolonged increased temperature, were the dominant drivers. The main influencing factors shaping the microbial communities within the different sampling months were temperature (permutation test, ANOVA, F = 4.59, p = 0.001), O₂ (F = 3.43, p = 0.002), SO42- (F = 3.24, p = 0.004), PO43- (F = 1.83, p = 0.043), and total Fe (F = 1.91, p = 0.052: Figure 5, Supplementary Table 3). The data indicated that long-term increased temperatures of spatial factors played an important role in shaping the microbial communities in the SED surface (Figure 5).

The BW community Shannon's H Alpha diversities between bays were significantly different in June (ANOVA, pairwise comparison, p < 0.01) and September (p < 0.01; Figure 3, Supplementary Table 3). Within the BW, the order Actinomycetales peaked in September (mean relative abundance, a 16.2 ± 5.4% heated bay vs. a 5.5 ± 1.9% control bay) with the highest abundance within the heated bay positively correlated with temperature (Pearson correlation, p = 0.0001, r = 0.36) and negatively with oxygen (p = 0.001, r = −0.26). In addition, Nanopelagicales-relative abundance increased in both bays in September (a 12.4 ± 3.4 heated bay vs. a 15.4 ± 3.5% control bay) and was positively correlated with higher temperature (p < 0.0001; Supplementary Table 6). The main contributor to this increased relative abundance in September was the genus Planktophila (mean ± s.d. of significant ASVs summarized on the genus level, 5.9 ± 1.5 vs. 7.6 ± 2.6%) that had higher abundances in the control bay in December (1.8 ± 0.9%) compared to the heated bay (0.9 ± 0.5%; Figure 6). Heated-bay BW showed a rapid change in Cyanobacteria compared to the control bay, while, in June, the cyanobacterial abundances were higher in the heated bay (23 ± 12.6%); they dropped rapidly in September (1.8 ± 0.5%), and stayed low in December (1.9 ± 0.8%). The Cyanobacteria order Synechococcales A peaked in June with positive correlation with temperature (p < 0.0001, r = 0.36) plus sulfate (p < 0.0001, r = 0.39) and higher relative abundance in the heated compared to control bay with mean relative abundances of 23 ± 12.6 vs. 15.7 ± 9% (Figure 4, Supplementary Figure 7). These high abundances were mainly due to the genus Cyanobium (16.6 ± 13.5%). Notably, the dominant bloom in June in the control bay was RCC307 (13.1 ± 9.3%) within the Synechococcales A order (Figure 6). The second bloom in September was dominated by Cyanobium in both bays but already decreased in the heated bay (1.3 ± 0.4%) compared to the control bay (5.9 ± 1.9%), indicating that peak blooms may have occurred earlier within the heated bay. With oxygen-rich spring BW, the peak in abundance of, e.g., Rhodobacteriales (positively correlated with oxygen, p < 0.0001, r = 0.41) and Flavobacteriales (p < 0.0001, r = 0.42; Figure 4, Supplementary Figure 7), was observed in both bays in April. The main contributors to the increased abundance were the genus Cellulophaga within the order Flavobacteriales, with higher abundances within the heated bay (14.4 ± 6.3%; Figure 6). On the order level, Betaproteobacteriales showed similar peak abundances in both bays within April (7.6 ± 2.4 vs. 5.2 ± 1%), while the genus Thiobacillus showed high abundance within the control bay in December (2. ± 1.3%) compared to the heated bay (0.8 ± 0.7%; Figure 6).

For the sediment, significantly higher Shannon's H Alpha diversities could be observed in April (ANOVA, pairwise comparison, p < 0.05) and June (p < 0.05; Figure 3, Supplementary Table 3) in the heated bay. Overall, the proportion of abundant microorganisms was more stable with the changing season but showed differences between the bays. For example, the order Anaerolineales showed overall higher concentrations in the heated (mean ± s.d. 1.2 ± 0.4%) compared to the control bay (0.9 ± 0.5%), while their abundance peaks in September (1.7 ± 0.8%) and June (1 ± 0.5%). The main contributor of the peak in September was the genus UBA6092 (0.81 ± 0.42%), which did not occur in higher abundance in the control bay (0.4 ± 0.3%). A similar pattern was observed within Rhodobacterales that had an overall low abundance in both bays but was, on average, 37.5% higher in the heated bay (Figure 4). On the other hand, the control bay showed a higher significant differential abundance of Betaproteobacteriales over all seasons compared to the heated bay. The main contributor was the genus Thiobacillus, which was the dominant taxon in all sampling sites and seasons within the control bay (4.8 ± 2.6%). Within the control bay, the order Desulfobacteriales, which was, on average, 28.6% higher abundant compared to the heated bay with peaks in September (3.5 ± 1.3%) and June (3.2 ± 1.1%), as well as Flavobacteriales, with an average of 32.3% higher abundance, was found.

In summary, cyanobacteria RCC307 and Cyanobium increased in relative abundance within the control bay bottom waters in the June and September samples, respectively. Within the heated bay, there was an increase of bottom water cyanobacteria-relative abundance in June due to the genus Cyanobium that decreased to a low-relative abundance in September. Furthermore, surface sediment microbial communities were more stable in composition over the year, showing overall higher relative abundances of Anaerolineales and Rhodobacteriales in the heated bay compared to higher relative abundances of Thiobacillus in the control bay.

PICRUSt2.0 predictions of BW metabolic responses (based on KO-identifiers) of the differentially abundant ASVs between the bays showed a significant increase in sulfur and nitrogen metabolism-associated genes within the heated bay, while a higher number of significantly different genes associated within photosynthesis pathways were found in the control bay (Figure 7, Supplementary Table 7). In addition, a significant increase of predicted genes associated with sulfur oxidation (SOX-system) and nitrogen metabolism-associated genes were found in the heated bay (Figure 7). Furthermore, nitrate and potentially nitrite uptake genes were significantly higher abundant in the heated bay and likely associated with Cyanobacteria (Flores et al., 2005). Genes associated with the dissimilatory nitrate reduction to ammonium (DNRA, genes nrfA and nrfH) were significantly more abundant in the heated bay (Figure 7, Supplementary Table 7). The control bay BW showed a different picture with the photosynthesis associated genes psb28 and psb27 significantly more abundant than the heated bay, which were likely associated with Cyanobacteria blooms (Figure 7). Further oxygenic photosynthetic light harvesting system-associated genes, including cpeZ, cpeY, and cpeU, had significantly higher relative abundance in the control bay (Figure 7, Supplementary Table 7).

Prediction of genes involved in SED community metabolism (Supplementary Table 7) showed a slightly shifted pattern within energy metabolism. While nitrogen metabolism genes were significantly differently abundant in the heated bay BW, there was no significant difference found within the heated bay sediment (Figure 7). Instead, genes related to the sulfur and methane metabolism were the main abundant predicted genes in the heated bay. In contrast, in total, more genes were significantly differently abundant in the control bay compared to the heated bay, including those coding for oxidative phosphorylation, photosynthesis, and nitrogen metabolism (Figure 7, Supplementary Table 7). Within the heated bay, sulfur oxidation soxAC plus conversion of methanol to formaldehyde within methane oxidation (mxaLKJGDCA) genes were significantly more abundant (Figure 7). In comparison, photosynthesis genes (cpeZYUTSREB) plus gapA, gap2, and tktB carbon fixation genes from photosynthetic organisms were significantly higher in the control bay (Figure 7). The higher abundance of Cyanobacteria within the control bay SED surface (Supplementary Figure 6) was further strengthened by the predicted significantly differentially abundant nrtC and nrtB genes for nitrate/nitrite uptake (Figure 7). Finally, sulfur metabolism played a minor role in the control bay differentially abundant genes, while methane oxidation gpmLB was relatively more abundant (Figure 7).