Redox-Specialized Bacterioplankton Metacommunity in a Temperate Estuary

This study explored the spatiotemporal dynamics of the bacterioplankton community composition in the Gulf of Finland (easternmost sub-basin of the Baltic Sea) based on phylogenetic analysis of 16S rDNA sequences acquired from community samples via pyrosequencing. Investigations of bacterioplankton in hydrographically complex systems provide good insight into the strategies by which microbes deal with spatiotemporal hydrographic gradients, as demonstrated by our research. Many ribotypes were closely affiliated with sequences isolated from environments with similar steep physiochemical gradients and/or seasonal changes, including seasonally anoxic estuaries. Hence, one of the main conclusions of this study is that marine ecosystems where oxygen and salinity gradients co-occur can be considered a habitat for a cosmopolitan metacommunity consisting of specialized groups occupying niches universal to such environments throughout the world. These niches revolve around functional capabilities to utilize different electron receptors and donors (including trace metal and single carbon compounds). On the other hand, temporal shifts in the bacterioplankton community composition at the surface layer were mainly connected to the seasonal succession of phytoplankton and the inflow of freshwater species. We also conclude that many relatively abundant populations are indigenous and well-established in the area.


Introduction
The world's oceans are the cradle of life; hence, the evolution of aquatic microorganisms for 3.5 billion years has produced enormous diversity and functional plasticity, only recently assessed by the sequencing of metagenomic DNA (pioneered by [1]). There are many aspects of microbial life that make the ecology of microorganisms different from that of macroorganisms [2], including intercontinental dispersion by winds [3] and the capability to persist in environmentally hostile conditions over a long period of time [4]. Aquatic microbes are essential for life on Earth [5,6] and therefore unveiling the mechanisms underlying the spatiotemporal dynamics of bacterioplankton community composition (BCC) remains the one of the most important issues in aquatic microbial ecology.
Over the last few decades, advances in sequencing technologies have revolutionized the power of the identification process for microorganisms and thereby revealed tremendous microbial diversity and plasticity in aquatic environments [7]. 16S rRNA gene-based investigations have contributed a massive number of sequences to databases and have revealed a comprehensive uncultured microbial diversity [8,9]. The ease of microbial community profiling has been effectively utilized to determine the biogeographic patterns of the most numerous and cosmopolitan marine bacterioplankton clades and, ultimately, to determine the functional traits that make them so successful [10][11][12]. More recently, high-throughput sequencing technologies have allowed for increasing the depth of investigation and thereby unveiled a rare biosphere that accounts for most of the observed phylogenetic diversity of bacterioplankton community [13,14]. This acts as a "seed bank" from where new dominant species can emerge when the environmental conditions change [4,15].
Consequently, species-sorting by the local environment has been demonstrated to be one of the main driving processes behind shaping the BCC [16][17][18][19]. However, in some cases, the assembly mechanism can be well explained by neutral models [20][21][22][23], by mass effects [24,25], or by the combination of several mechanisms [26][27][28]; the relative importance of these mechanisms may change over time [29].
In addition to unique environmental conditions, similarities to other communities have to be considered in order to identify processes underlying the assembly of local microbial communities [23]. Hence, in this study, we combined environmental factors with phylogenetic affiliations of relatively abundant populations for that purpose. Hence, special attention was paid to associations within between ribotypes, because these interactions can have stronger correlative relationships compared to relationships between bacteria and eukaryotes, or bacteria and abiotic environmental factors [30]. The co-occurrence of networks of dominant bacterial ribotypes isolated from the marine oxygen minimum zone (OMZ) throughout the world has revealed a pattern of cosmopolitan key species filling redox-driven niches [31]. These niches revolve around functional capabilities to utilize different electron receptors and donors [32]. Next important step towards a better understanding of these microbial communities inhabiting OMZ is to define shared or specialized metabolic subsystems in different oceanic provinces [31]. Our results contribute to this effort.
The Baltic Sea is one of the largest brackish basins of the world, characterized by a long residence time. Therefore, it is not just a mixing zone for fresh water and marine species, but a habitat for microbes specialized for brackish water, which has been illustrated by the spread of different bacterial populations throughout the salinity gradient of the Baltic Sea [33][34][35]. Unlike the diversity of macro-organisms, the BCC does not decline within a salinity gradient [35].
The Gulf of Finland is the easternmost sub-basin of the Baltic Sea. The strong stratification in the central part of the gulf due to the seasonal thermocline and permanent halocline often hinders mixing in the water column [36]. Eutrophication-driven phytoplankton production leads to increased sedimentation of organic matter and hence increased consumption of oxygen for which atmospheric and photosynthetic re-oxygenation cannot compensate [37].
Furthermore, the Gulf of Finland is directly connected to the Baltic Proper, where the anoxic zone is permanent and therefore inhabited by well-established anaerobic ecotypes typical of the Baltic Proper [36,38,39]. During oxygen deficiency, certain microbes are capable of using terminal electron acceptors other than oxygen (e.g. NO 3 -, SO 4 2and metal oxides). Epsilonproteobacterial Sulfurimonas gotlandica clade GD1 has been demonstrated to be one of the most numerous chemolithoautotrophic bacteria present in the OMZ of the central Baltic Sea [40,41]. This clade has been shown to spread into the anoxic zone of the Gulf of Finland [42].
As a temperate estuary, the Gulf of Finland undergoes many seasonal changes in environmental parameters such as ice coverage, water temperature, solar radiance, inorganic nutrients, etc., which contribute to the recurring succession of phytoplankton. Variation in biotic and abiotic factors also leads to the recurring succession of bacterioplankton in the Baltic Sea [34,43,44].
The goals of present study were to investigate the spatial and temporal dynamics of the BCC of the Gulf of Finland during the spring to summer transition in order to determine the main factors driving the bacterioplankton community assembly. To these ends, we used pyrosequencing of 16S rRNA genes from community DNA samples, which were collected in parallel with monitoring of physicochemical parameters and phytoplankton community composition. Many annually recurring environmental shifts took place during the study period, mainly the formation of a thermocline, the depletion of oxygen in the deeper layers and the shift from the spring to summer phytoplankton community.

Ethics statement
No specific permits were required for the described field studies. Our study area is not privately owned or protected in any way. The study did not involve endangered or protected species.

Sample collection and extraction of community DNA
Water sampling was performed aboard the RV Salme in the spring and summer of 2011 on two transects in the central part of the Gulf of Finland (Fig 1). The coordinates of sampling stations are listed in Table 1. Samples were collected from three depths: at 5 m, 40 m and about 5 m above the seafloor; detailed information about each sample is provided in Table 2. A rosette sampler (M1018, General Oceanics) equipped with Niskin water samplers (volume 1.7 l) was used for sampling. For background information, the depth profiles of conductivity and temperature were obtained with a SBE19plus CTD probe (conductivity-temperature-depth probe, Sea-Bird Electronics), chlorophyll a fluorescence with WETStar fluorometer (WETLabs) and dissolved oxygen with SBE 43 probe (Sea-Bird Electronics).
The samples for nutrient (NO 2 − + NO 3 − and PO 4 3− ) analysis were deep-frozen at −20°C after collection and analyzed at a shore-based laboratory using a Lachat QuikChem 8500 Series 2 automatic nutrient analyzer (Lachat Instruments, Hach Company). The nutrient analyses were performed according to the ISO 15681-1 method for PO 4 3− and ISO 13395 method for The lower detection ranges for PO 4 3− and NO 2 − + NO 3 − were 0.02 and 0.03 μmol l -1 , respectively. Water samples for microbial community extraction were collected into sterile bottles (Nalgene) and immediately filtered through 0.2 μm filters (Whatman, Puradisc FP 30) after preliminary filtration through 5.0 μm prefilters (Whatman, Puradisc FP 30). The scheme of the filtration system was described by Laas et al. (2014) [42]. The sample volume varied between 0.5 and 1.0 liters. Filters were kept frozen at -20°C until community DNA was extracted with a PowerSoil DNA Isolation Kit (MO BIO Laboratories, Inc.). A few modifications were made to the protocol: syringe filters were incubated with the lysis buffer in the casing at 60°C for 30 min and then the eluate was removed.

Amplification of bacterial 16S rRNA gene sequences
The bacterial 16S rRNA gene V1-V2 hypervariable regions were amplified in two polymerase chain reactions (PCR). For the first reaction, universal bacterial primers BSF8 and BSR357 were complemented with 8 nt barcode and partial adapter sequences (Table 3) [46]. PCR was performed with Smart-Taq Hot Red 2X PCR Mix (Naxo, Estonia), 1 μl of extracted DNA and 0.2 μM each primer, using the following cycling parameters: 15 min denaturation followed by three cycles (30 sec at 95°C, 30 sec at 50°C, 60 sec at 72°C), 28 cycles (30 sec at 95°C, 30 sec at 65°C, 60 sec at 72°C) and a final extension at 72°C for 7 min. To achieve full length sequencing adapters, second PCR amplification was performed.
The second reaction was run with Smart-Taq Hot Red 2X PCR Mix (Naxo, Estonia), 1 μl of 10 X diluted amplicon, 0.2 μM each primer; using the following cycling parameters: 15 min denaturation followed by five cycles (30 sec at 95°C, 30 sec at 62°C, 60 sec at 72°C), 20 cycles (30 sec at 95°C, 60 sec at 72°C) and a final extension at 72°C for 10 min. PCR reactions were run on a thermal cycler (model 2720, Applied Biosystems). Each PCR product was gel purified on a 1.5% agarose gel. DNA was isolated using the QIAquick Gel extraction kit (Qiagen, Inc.). DNA concentrations were measured with a Qubit fluorometer (Invitrogen). Sequencing was performed on a Roche GS FLX next generation sequencing platform (IMGM Laboratories).

Bioinformatics and statistics
Reads with low quality and those shorter than 150 bp (basepairs) were removed from the dataset. The PyroNoise algorithm was used to discard homopolymer-derived errors [47] and UCHIME to remove chimeric DNA sequences caused by PCR errors [48]. OTUs (Operational Taxonomic Units) were defined using the average neighbor-clustering algorithm of MOTHUR 1.19.1 [49] with a 97% similarity threshold. Reference sequences were selected from the SILVA ribosomal RNA database [50]. Taxonomic assignments were processed by the Ribosomal Database Project (RDP) naïve Bayesian Classifier [51]. For database affiliations RDP Seqmatch [52] and BLAST [53] search algorithms were used against RDP and a NCBI (National Center for Biotechnology Information) nucleotide databases, respectively. Chloroplast and mitochondrial sequences were discarded from the dataset. Statistical analyses were carried out with the R program version 2.14.0 (http://www.r-project.org), ACE (Abundance-based Coverage Estimation) [54] and Chao1 [55] richness estimates; multivariate statistics were calculated using the VEGAN package [56]. The similarity matrices and clustering were generated using the gplots package [57]. The sequences have been deposited into GenBank (accession numbers from KM489611 to KM491167).

Environmental parameters structuring the bacterial community
The sampling took place from April to July in 2011 and the bacterial 16S rRNA gene libraries were generated from three sampling depths: (i) surface water at 5 m, (ii) the intermediate layer  Table 2). These depths were chosen because, in summer, the gulf is stratified into three layers. Four hydrographic profiles are visualized in Fig 2 to provide overview of changes of environmental parameters in the water column throughout the sampling period. The 5 m and 40 m horizons always remained above the permanent halocline and at those depths the salinity varied in the range of 4.8-6.6 g kg -1 and 6.7-7.2 g kg -1 , respectively. Salinity decreased from west to east, with the lowest value recorded in the easternmost surface sample. The near-bottom layer samples were collected at the halocline or below it and the salinity ranged from 8.5-10.4 g kg -1 . The temperature varied little in the near-bottom layer (3.7-5.0°C in April and 0.1-4.0°C in July), but at the 5 m depth, warming was observed from 0.8-0.9°C in April to 17.1-18.7°C in July. The most remarkable change in temperature occurred between May and June, when a well-pronounced thermocline started to form at 10-25 m.
Oxygen conditions were hypoxic or anoxic in deeper sampling layer (> 63 m) throughout the sampling period. Oxygen concentrations at the 5 m depth were significantly (p < 0.01, Student's t-test) lower during the summer (10.0-12.3 mg l -1 ) compared with the spring (13.9-16.0 mg l -1 ). At 40 m depth, the overall tendency remained similar to the upper layer, but on two occasions in May, the oxygen concentrations were below 9.5 mg l -1 .
Chlorophyll a concentrations were highest in April (4.2-12.5 mg m -3 , Table 2), when phytoplankton was dominated by Diatomphyceae (data not shown). Also, concentrations of inorganic nutrients at 5 m depth we higher during the spring months compared with June or July. Nitrites and nitrates remained at concentration levels from the detection limit to 1.2 μmol l -1 and phosphates varied between 0.6-0.8 μmol l -1 . In June-July, the concentrations of nitrites and nitrates remained below the detection limit and phytoplankton was dominated by diazotrophic Nostocophyceae (data not shown). The concentration of phosphates was also lower in the surface layer during the summer months (Table 2).

Bacterioplankton community diversity
A total of 73,494 partial rRNA gene sequences and 48 different samples were used in this analysis, on average 2,161 sequences per sample. The minimum number of sequences per sample was set to 350. In Table 2, the number of observed OTUs is accompanied with the Chao1 and ACE species richness estimates for each sample. Overall, the number of OTUs was significantly higher in the near-bottom layer (on average 144, SD = 34) than at 5 m and 40 m (98, SD = 44 and 90, SD = 43; respectively). Based on the rarefaction curves outlined from samples (S1 Fig), we assume that deeper sequencing (more sequences per sample) would have resulted in significantly higher estimates. Therefore, in this study, we could not study rare species and therefore concentrated on abundant members of the bacterioplankton community.
A descriptive multivariate statistical method, i.e. detrended correspondence analysis (DCA), was used to describe the overall OTU abundance and thereafter a vector-fitting procedure was applied to determine which environmental variables significantly related to BCC patterns (Fig 3, Table 5). As a result, oxygen, depth and salinity were distinguished as the most important co-varying environmental factors (r 2 and p values are presented in Table 5). Temperature, sampling time and Chl α could be identified as being less significant to the DCA space (in declining order). The geographic location (longitude and latitude) was rendered non-significant. The DCA also indicates how measured environmental parameters co-varied. To further investigate the dynamics of BCC, Pearson correlation analysis was carried out based on the relative abundance of OTUs, which was used to construct a similarity matrix (Fig 4). Organized according to clustering of the similarity matrix, all community profiles are visualized on the class level in Fig 5 for an overview and the OTU level (97% similarity) in Fig 6 to provide detailed insight into the relative abundance of dominant ribotypes. On a broader scale, the communities were divided into three groups: (i) surface communities during the summer months that were dominated by the unicellular cyanobacterium Synechococcus (OTU BSNS2840); (ii) hypoxic/anoxic near-bottom communities that contained a large fraction of chemolithotrophic Sulfurimonas (OTU BSNS3177); and (iii) a larger group with all remaining communities with notable subdivisions among them.
To examine how the BCC was related to environmental factors, each measured environmental parameter was correlated to the relative abundance of the dominant bacterial classes (Fig 7) of OTUs (Fig 8). These heatmaps illustrate how occurrence patterns of classes differ in relation to environmental conditions and that on OTU level there are subdivisions within these groups. Interactions between microbes were as important as other environmental factors. Therefore, similar correlation analyses were carried out between dominant ribotypes (Fig 9). The clustering achieved in this way is slightly different and provides additional information. The efficiency of this approach is demonstrated by the fact that OTUs with related database matches (isolation source, etc.) clustered together. Therefore, relatively abundant ribotypes with the closest affiliations are given in the same order in Table 4.

Connection between locally established bacterioplankton community and global redox-specific metacommunity
We used 16S rRNA gene-based community profiling to identify spatiotemporal patterns of bacterial picoplankton (including 0.22 to 5 μm fraction) in the Gulf of Finland. As discussed in detail below, the collected data demonstrate that the dynamics of the BCC in the Gulf of Finland exhibits striking parallels with other OMZ that are also characterized by salinity gradients and a temperate climate, like Chesapeake Bay and the Saanich Inlet [25,58]. This suggests that, in addition to the dynamics of different electron acceptors and donors, salinity also plays an important role. In community ecology, there is an increasingly popular concept of metacommunities, defined as a set of local communities that are linked by the dispersal of multiple interacting species [59]. Considering the interconnectedness of aquatic ecosystems, we propose that bacterioplankton communities in the OMZ can be considered as a globally distributed redoxspecialized metacommunity. In such a framework, OMZs occurring in freshwater and marine water mixing zones should form a salinity-dependent subsystem of this metacommunity. Due high selectivity by salinity-and redox-driven niche partitioning in these systems, species-sorting becomes the driving force behind community assembly in these systems.
The overall phylogenetic makeup of the bacterioplankton community observed on the bacterial class level resembled those routinely found in Baltic Sea pelagic waters using culture     independent methods, although an underrepresented fraction of Gammaproteobacteria and Verrucomicrobia was noted [34,39,42,44,60]. Moreover, many relatively abundant OTUs are closely affiliated with sequences previously isolated from the Baltic Sea (Table 4), as expected due to the long residence time of the Baltic Sea and annually reoccurring patterns of bacterioplankton succession [44]. Two OTUs that managed in some cases to contribute to over 50% of the BCC (Fig 6) can be considered permanent and well-adapted local populations. One of these OTUs, BSNS2840, had an identical match to a sequence isolated from laminae of the Gulf of Finland sediments dating back to the Late Litorina Sea [61]. Hence, this population may have been present in the area for over 3000 years. It was classified as a member of the cyanobacterial clade GpIIa (Synechococcus by database affiliation) and it contributed up to 88.7% of the BCC in the surface layer in June, making it the most abundant picocyanobacterium of the dataset (Fig 6). OTU BSNS2840 was also numerous in the winter picoplankton community in the same sampling area [42].
The prevalence of Cyanobacteria during summer is typical of the seasonal succession of phytoplankton in the area and is caused by multiple environmental conditions, most notably the increase in temperature and the depletion of nitrates in the surface layer after the spring bloom (Table 2), which gives a distinct advantage to diazotrophic cyanobacteria [62]. The  variability of the fraction contributed by Cyanobacteria within samples that were collected during same cruise indicates patchy nature of the bloom. In addition, it is important to point out that Synechococcus has been shown to occupy the anoxic and dark zone in the central Baltic Sea [38,63] and Chesapeake Bay [25]. The underlying mechanisms behind this phenomenon remain unexplored to our knowledge. The most abundant OTU in the hypoxic/suboxic layer, BSNS3177, was classified as Sulfurimonas (genus of Epsilonproteobacteria). One characteristic feature of chemolithoautotrophic bacteria in the Baltic Sea is that the majority of cells belong to the Sulfurimonas GD17 group [64], which oxidizes H 2 S with NO 3-. The prevalence of single domant strain of epsilonproteobacterium is also reflected in our results. Ribotype BSNS3177 contributed to a large fraction of the BCC in the hypoxic/suboxic near-bottom layer at the westernmost stations and at AP2 and AP5 in mid-summer (Fig 6). This suggests dispersion from the Baltic Proper to the Gulf of Finland, as proposed by our previous study [42]. Moreover, the Sulfurimonas strain co-occurred with OTU BSNS3126 (Fig 9), classified as Bacteroidetes, and both were affiliated with sequences isolated earlier from the anoxic zone of central Baltic Sea (Table 4) [65]. Linkage to the same isolation source was also demonstrated by the co-occurring pair of ribotypes BSNS3168 and BSNS3178, classified as representatives of Gammaproteobacteria (affiliated with the SUP05 clade) and Flavobacteriaceae, respectively. These correlations probably represent cooperation activities. The Epsilonand Gammaproteobacteria have been identified as major chemolithoautotrophic groups in the central Baltic Sea [40,41,[64][65][66][67] and also globally in other marine OMZs [31,68]. The GD17 group is capable of a chemotactic response to nitrate [69]. The chemotactic swimming strategies of marine bacteria are especially significant in patchy nutrient seascapes [70], and in OMZs, the corresponding genes have been shown to be over-represented [39,71].
Heterotrophic Alphaproteobacteria dominated in the upper oxygenated layer of water. In addition, Candidatus Pelagibacter (ribotypes BSNS3076, BSNS3084 and BSNS3171 accounted for 15.9%) was the most abundant genus in the whole dataset. Candidatus Pelagibacter belongs to the SAR11 clade, which is the most abundant type of organism in the world's oceans [12,72,73] and has been found to be numerous in the Baltic Sea as well [35,42]. Another prevalent alphaproteobacterial OTU in oxygenated waters, BSNS3149, was classified as Rhodobacteraceae (order Rhodobacterales). It had an identical match with a sequence isolated from a coastal North Sea diatom bloom [74]. Both groups, Rhodobacterales and SAR11, contribute to the conversion of dimethylsulfoniopropionate to dimethylsulfide [75]. Interestingly, both groups have been found to be numerous in OMZs in both oxic and anoxic layers, but their metabolic adaptations for anaerobic growth remain to be uncharacterized [31].

Importance of methane and trace metal gradients
Methane is mainly generated in the sediment from which it transfers into the water column. A methane gradient with decreasing concentrations toward the water surface is usually observed in the Baltic Sea [76]. Microbial oxidation of methane in the water column represents an important sink of methane. Different electron acceptors can be used for methane oxidation (sulfate, nitrate, nitrite, iron and manganese), but using oxygen produces a considerably larger energy yield [77].
As a consequence, the hotspot for microbial methane oxidation activity in central Baltic Sea is at the upper part of the redoxcline [78]. Interestingly, the corresponding group consists of only type I methanotrophic bacteria [79,80]. Our results suggest that the same group is also represented in the Gulf of Finland with gammaproteobacterial OTU BSNS3164 (Methylobacter) as the most abundant representative. It is important to consider that our samples were collected above and under the redoxcline. In addition, OTU BSNS3143 corresponded via databases to Methylosinus trichosporium, a type II methane-oxidizing alphaproteobacterium previously found in surface sediments of the Gotland Deep [81]. Considerably more relevant was betaproteobacterial OTU BSNS3170, classified as Methylophilus, an obligate methylotroph common to the Baltic Sea [38,42,82]. All these three groups were more abundant at the 40 m depth (Fig 6).
Redox-sensitive trace metals (e.g. iron and manganese) provide another energy source for microbes and can serve as electron donors above the redoxcline. When oxygenated, insoluble compounds form aggregates that will sink below the redoxcline and can be used as electron acceptors [83][84][85]. Our results suggest that members of Actinomycetales (order of Actinobacteria) may inhabit this niche, because several relatively abundant OTUs (BSNS2659, BSNS3154 and BSNS3169) were classified as Ilumatobacter. All known and described members of Acidimicrobiaceae are capable of iron oxidation [86]. The found OTUs exhibited different patterns of distribution indicating separate strategies. For example, OTU BSNS2659 was more abundant at 5 m in July. However, two other Ilumatobacteria and three unclassified Actinobacteria (OTUs BSNS3032, BSNS3120, BSNS3133) occurred mainly at 40 m and in the near-bottom layers, supporting the iron-oxidation hypothesis. Interestingly, almost all of these OTUs (except OTU BSNS3169) had nearly identical matches to brackish and temporarily anoxic estuaries (Table 4).

Spring bloom associated community influenced by freshwater clades
The Gulf of Finland has a relatively large input of freshwater compared to the Baltic Proper and southern basins. This impacts the BCC through the salinity distribution of the gulf and the influx of populations originating from freshwater sources [33][34][35]44]. Consequently, some relatively abundant OTUs were highly similar to sequences obtained from rivers and lakes (Table 4). Overall, the prevalence of phylogenetic groups considered to be characteristic of freshwater ecosystems, like Actinobacteria and Betaproteobacteria, exhibited a positive correlation with longitude (Fig 7), i.e. increasing in abundance towards less saline parts of the gulf. This concurs with previous reports of BCC dynamics along the salinity gradient of the Baltic Sea [35].
Annually recurring phytoplankton spring blooms in the Gulf of Finland are co-dominated by diatoms and dinoflagellates and exhibit high spatiotemporal variability [87]. The alternating dominance of bottom-up and top-down interactions result in a succession of heterotrophic bacterial groups with different growth strategies throughout the phytoplankton bloom [74,88]. Although our sampling period covered only part of the spring bloom, groups representing different strategies could be distinguished.
Bacterial lineages considered fast growing "opportunistic" types that utilize dissolved organic matter (DOM) like Betaproteobacteria and Flavobacteria [88][89][90] were accompanied by more slow-growing and grazing-resistant members of the AcI lineage of Actinobacteria [88,91,92]. Most of the OTUs connected to the spring bloom were clustered together by the correlation analysis (upper 12 OTUs in Fig 7), including the discussed unclassified Rhodobacteraceae that were affiliated with the diatom bloom.
Members of Flavobacteria are often overrepresented in the spring BCC of the Baltic Sea and are major contributors to the degradation of high molecular weight carbon [39,82,93]. OTUs BSNS2870, BSNS3078 (both classified as Flavobacteriaceae) and BSNS3163 (unclassified Bacteroidetes) were closely related to sequences previously isolated from the Baltic Sea (Table 4). We identified several OTUs classified as Polaribacter (member of Flavobacteria), with OTU BSNS3072 as the only relatively abundant representative. This group contributes to the degradation of phytoplankton-derived organic matter via high expression of sulfatases [74].
Overall, Betaproteobacteria had the strongest correlation with Chl a concentrations, especially unclassified Comamonadaceae (BSNS3044) which stood out on the OTU level (Figs 6 and 7, respectively). Members of Comamonadaceae have been previously identified in spring bacterioplankton communities on several occasions [39,82]. Moreover, mesocosm experiments have shown a negative impact of higher temperature on this psychrotolerant lineage (Lindh et al., 2012).
We identified three relatively abundant members of Actinobacteria (OTUs BSNS3107, BSNS3110 and BSNS3157) that were closely related to the freshwater AcI lineage characterized by small cell sizes. The presence of small bacteria can be considered a clear indication of grazing pressure by bacterivorous nanoflagellates. Metagenomic profiling of pelagic and benthic bacteria during the spring (Landsort Deep) reveled that genes for the degradation of polyaromatic hydrocarbons (like cellulose and chitin) belong mainly to Actinobacteria and more specifically to Mycobacterium [39]. These genes were mainly found in the sediments, and because they are affiliated with aerobic groups, the authors considered sedimentation most likely. In our dataset, OTU BSNS3016, classified as Corynebacterineae and closely affiliated with Mycobacterium, comprised 8.7% of the total sequences isolated at a depth of 40 m in the beginning of June; therefore, our results support the sedimentation hypothesis.

Conclusions
The availability of electron acceptors is a critical determinant of the marine ecosystem structure, and we conclude that oxygen concentration is a major environmental factor impacting the BCC in the near-bottom layer of the Gulf of Finland. We conclude that chemolithotrophic groups dispersing from the central Baltic Sea become dominant members of the BCC in the suboxic/anoxic layer when it is formed. Some members of Actinobacteria inhabit the layer above the redoxcline and most probably contributing to the oxidation of ferrous iron. Our results also led to the conclusion that the Gulf of Finland has likely a more diverse composition of methanotrophic bacteria than the central Baltic Sea. The BCC at the surface layer was strongly impacted by phytoplankton seasonal succession, as many abundant OTUs in April and May could be associated with the spring phytoplankton bloom. In addition, in the beginning of summer, inorganic nutrient depletion and rising temperatures led to the proliferation of picocyanobacteria. We determined one dominant lineage of Synechococcus, which due affiliation contributes to the conclusion that some well-established phylogenetic lineages have persisted in the area for over 3000 years. Furthermore, our results support the emerging pattern of related microbes occupying the OMZs throughout the world and suggests that core of the bacterioplankton community of the Gulf of Finland is part of that a redox-specialized bacterial network, which by definition can be considered a metacommunity.