Water samples from 5 m depth were collected from 9 stations in the Gulf of Finland (GoF), the Baltic Sea, from winter to late spring (2013 Dec – 2014 May; 4 stations), and from summer to autumn (2013 May-Oct; 5 stations; Fig. 1, Suppl. material 2, Table S1). In general, the stations were sampled once per month, however stations A12, AP2, and AP5 were sampled on multiple occasions in May, June, and July (Suppl. material 2, Table S1). The data from those occasions was pooled and displayed based on the month in the manuscript. A rosette sampler (M1018, General Oceanics, U.S.A.) equipped with 8 Niskin water samplers (volume 1.7 l, model 1010, General Oceanics, U.S.A.) was used for sampling. Salinity and water temperature were measured with an OS320plus CTD probe (Idronaut s.r.l, Italy). The samples were collected into sterile bottles (342020-100, Nalgene, Thermo Scientific Nalgene Inc., U.S.A.) and filtered immediately onboard the RV through 5.0 μm and 0.2 μm filters (Whatman 10462300, Puradisc FP 30, Merck, Germany) as described by Laas and others (2014, 2015). The sample volume varied between 0.5 and 1.0 liters depending on the amount of particulate matter in the samples. Filters were kept frozen at -20°C until the DNA extraction by CTAB (Lian et al. 2001) with the following modifications: 300 μl of CTAB solution was pipetted into the syringe filters and the luer inlets and outlets were capped with Combi-Stopper Closing Cones (B. Braun Melsungen AG, Germany) to prevent the CTAB buffer leaking during the incubation. After incubation at 65 °C for 1 hour, 400 μl of TE buffer was added to the DNA extract taken from the filters to facilitate purification by chloroform-isoamyl alcohol (24:1) following the protocol by Lian and others (2001).

Figure 1.

Sampling locations in the Gulf of Finland (GoF), the Baltic Sea. Stations labeled with “A” were sampled from December to May in 2013-2014; stations labeled with “AP” were sampled monthly from May to August-October in 2013.

DNA extracted from 0.2 and 5 µm filters was mixed in equal volume and used as a template for paired-end library preparation. Eukaryotic species were targeted by universal primers for the 18S (V7-V9 region; SSR-F1289-sn, F: TGGAGYGATHTGTCTGGTTDATTCCG; SSR-R1772-sn, R: TCACCTACGGAWACCTTGTTACG, Dzhembekova et al. 2018) and 28S (D1 region; D1R_M, F: ACCCGCYGAAYTTAAGCATA, modified by S. Nagai after Scholin et al. (1994); D1_YK_R: AGWCCGATAGCRMACAAGTA, Sildever et al. 2019) rRNA genes. In addition, universal primers mainly targeting dinoflagellates (28S rRNA gene, Dino28S_325F-L: CAAGTACCATGAGGGAAAG, Dino28S_759R_L: GACTCCTTGGTCCGTGTTTC, Nagai et al. unpublished data) were also used. For prokaryotic species, universal 16S rRNA gene primers (V3-V4 region; Bakt_341F: CCTACGGGNGGCWGCAG, Bakt_805R: GACTACHVGGGTATCTAATCC, Sinclair et al. 2015) were employed. Two-step PCR for the construction of paired-end libraries and HTS on Illumina Miseq 300 PE platform (Illumina Inc., USA) followed the protocol described by Dzhembekova and others (2017). Briefly, in the first-round PCR, the target regions were amplified and in the second PCR, the index primers were attached to each sample. For the first round PCR, the reaction mixture (25 µl) contained 0.2 mM of each dNTP; 1 x PCR buffer; 1.5 mM Mg2+; 1.0 U KOD-Plus-ver.2 (TOYOBO, Osaka, Japan), 0.4 µM of each primer, and 1 µl of template DNA. The amplification was done using the Bio-Rad T100 Thermal Cycler with the initial denaturation at 94˚C for 3 min, 30–32 cycles at 94˚C for 15 s, 56˚C for 30 s, and 68˚C for 40 s. PCR amplification was verified by 1.5% agarose gel electrophoresis. The PCR products were purified using an Agencourt AMPure XP (BECKMAN COULTER, Life Sciences, Brea, California, USA) and eluted in 25 mL of TE buffer following the manufacturer protocol. When unexpected PCR bands appeared, only the target amplicon was excised from the gel and purified using a High-Pure PCR Product Purification Kit (Roche, Basel, Switzerland) and eluted in 25 mL of TE buffer following the manufacturer protocol.

The purified PCR products were used as a template in the second-round PCR. The reaction mixture composition was the same as in the first PCR, however, the total volume of the mixture was 50 µl, 2 µl of DNA template, and the following primers were used: 5’AATGATACGGCGACCACCGAGATCTACAC-8 bp index-ACACTCTTTCCCTACACGACGC (forward) and 5’CAAGCA GAAGACGGCATACGAGAT- 8 bp index GTGACTGGAGTTCAGACGTGTG (reverse). The PCR conditions were as follows: initial denaturation at 94˚C for 3 min, 8–12 cycles at 94˚C for 15 s, 56˚C for 30 s, and 68˚C for 40 s. PCR amplification was verified by agarose gel electrophoresis, and the PCR products were purified using an Agencourt AMPure XP. The amplified PCR products were quantified, pooled in equal concentrations, and stored at -30 ˚C until sequencing. Sequencing on Illumina Miseq 300 PE platform was ordered from the FASMAC, Co. Ltd., Japan.

Treatment of obtained sequences, selection of operational taxonomic units (OTUs), and taxonomic characterization of OTUs was done according to the workflow described by Dzhembekova et al. (2017) with the exception that sequences longer than 300 bp were truncated to 300 bp by trimming the 3′ tails. The trimmed sequences shorter than 250 bp were filtered out. Demultiplexing and trimming were performed using Trimmomatic version 0.35 (http://www.usadellab.org/cms/?page=trimmomatic). The remaining sequences were merged into paired reads using Usearch version 8.0.1517 (http://www.drive5.com/usearch/) and singletons were removed. Sequences were aligned using Clustal Omega v 1.2.0. (http://www.clustal.org/omega/) and only sequences that were contained in more than 75% of the read positions were extracted. Filtering and a part of the multiple alignment process were performed using the screen.seqs and filter. seqs commands in Mothur, as described in the Miseq SOP (http://www.mothur.org./wiki/MiSeq_SOP) (Schloss et al., 2011). Erroneous and chimeric sequences were detected and removed using the pre.cluster (diffs = 4) and chimera.uchime (minh = 0.1; http://drive5.com/usearch/manual/uchime_algo.html) (Edgar et al., 2011) commands in Mothur, respectively. Using the unique.seqs command of Mothur, the same sequences were collected into operational taxonomic units (OTUs). The contig sequences were counted as OTUs by count.seqs and used for the subsequent taxonomic identification analysis. Demultiplexed, filtered, but untrimmed sequence data were deposited in the DDBJ Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra/index-e.html; BioProject (PRJDB12325) and BioSamples (SAMD00406059-SAMD00406250), DRA accession number DRA012806).

The selected OTUs were taxonomically identified as follows: for the BLAST search, a subset of NCBI non-redundant nucleotides (NT) consisting of sequences that satisfied the below-mentioned conditions was prepared. The sequences used to identify OTUs were downloaded from the GenBank based on the following search criteria: “ribosomal”, “rrna” or “rdna”, and the search excluded keywords associated with “metagenome”, “uncultured” and “environmental”. The sequences of retrieved GenBank IDs were downloaded from the nucleotide database at the NCBI FTP server on March 23, 2019, and were used to construct a template sequence database. Subsequently, the taxonomic characterization for each OTU was performed using a BLAST search (Cheung et al., 2010) in NCBI BLAST+2.2.26+ (Camacho et al., 2009) with default parameters. The nucleotide subset described above was treated as a database, and all OTU representative sequences were treated as queries. Subsequently, taxonomic information was obtained from the BLAST top bit-score hit for each query sequence.

For eukaryotes, only OTUs with blast top hit similarity of ≥ 0.99% and a minimum query cover of 70% were included in the further analysis. OTUs associated with multiple records from the same genus were merged if the multiple records consisted of sp. of the same genus or when a single species was clustered together with other records identified as sp. from the same genus. For prokaryotes, only OTUs with > 50 total sequence reads were included in the further analysis. This rule was not applied for cyanobacteria that were identified to species level. The phylum and class/order level taxonomic groupings were also confirmed based on the World Register of Marine Species (www.marinespecies.org) or the UniProt taxonomy database (www.uniprot.org). For OTUs associated with Dinophyceae and Bacillariophyceae, the class level grouping is shown instead of phylum level as both represent important and dominant phytoplankton groups. The autotrophic ciliate, Mesodinium rubrum, is included both in the phyto- and zooplankton results to show its dynamics concerning both groups as it is usually analyzed together with phytoplankton.

The presence of toxin-producing HAB species among the detected OTUs was investigated by comparing the list of OTUs associated with phytoplankton species against the IOC-UNESCO Taxonomic Reference List of Harmful Micro Algae (Moestrup et al. 2021). Only OTUs associated with a single toxin-producing HAB species were included. To avoid the inclusion of OTUs with ambiguous identities, representative sequences of all OTUs associated with the toxic HAB species were also manually BLAST searched from the GenBank online database (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The validity of the taxonomic names was checked against the AlgaeBase (http://www.algaebase.org/) and the World Register of Marine Species (http://www.marinespecies.org/). Only OTUs that had the ≥ 99.0 % BLAST top hit similarity with the associated single species were included in the further analysis.

The numbers of shared OTUs among different sampling years were counted and depicted as a Venn diagram by using the “vegan” package (Oksanen et al. 2020) in R version 4.0.2 (R Core Team, 2021). Presence-absence OTU data was used for all statistical analysis, and non-rarefied relative abundance data was used for visualizations in Fig. 3–7. Non-metric multidimensional scaling (NMDS) based on the Jaccard similarity index was applied to investigate the biodiversity among stations and sampling months by using “vegan”. Shepard diagrams were employed to assess the goodness of fit of the resulting NMDS analysis (Suppl. material 1, Fig. S1). Differences in biodiversity between the stations and sampling months were investigated by PERMANOVA using the Jaccard similarity index and 999 permutations and the assumption for homogeneity of multivariate dispersions was tested by ANOVA (betadisper function in “vegan”). Spearman correlations were applied to determine associations between common, dominant, or toxin-producing phytoplankton species and bacteria/zooplankton based on presence/absence data from 18S and 16S datasets. P-value was adjusted for multiple comparisons by Benjamini-Hochberg correction (Benjamini and Hochberg 1995).

During winter-spring, salinity ranged from 5.11 g kg^-1 to 6.56 g kg^-1 with higher salinities at stations A1 and A4 located in the western GoF (Fig. 2; Suppl. material 1, Fig. S2). The physical-chemical background during winter-spring (2013–2014) is described in more detail by Lips and others (2017). In summer-autumn, salinity ranged from 4.58 g kg^-1 to 5.70 g kg^-1 with no notable differences between the stations. In winter-spring, the water temperature ranged from 0.71 ˚C to 5.97 ˚C with the lowest temperatures in February-March. In summer-autumn, the temperature ranged from 4.69 ˚C to 18.47 ˚C, with the highest temperatures in July-August (Fig. 2; Suppl. material 1, Fig. S2).

Figure 2.

Seasonal changes in salinity and temperature during the sampling months (May-December 2013; January-May 2014). Data until 40 m depth is shown based on stations AP5 and A8.

The number of eukaryotic OTUs detected was 1755, 551, and 229 based on 18S, 28S_E (universal eukaryote), and 28S_D (universal dinoflagellates) markers, respectively. At ≥ 0.99% BLAST top hit similarity level and a minimum of 70% query cover length, the highest number of OTUs (n=314) was identified by 18S followed by 28S_E (95) and 28S_D (76) primers (Table 1, Suppl. material 2, Tables S2–S4). The eukaryotic OTUs at ≥ 0.99 % BLAST top hit similarity level originated from 9 unique supergroups and 34 phyla (Suppl. material 1, Figs S3, S4). The highest number of OTUs detected were associated with Dinophyceae, Bacillariophyta, and Chlorophyta based on 18S universal primers (Table 1). Based on 28S_E, the phyla with the highest number of OTUs were Ascomycota and Basidiomycota, followed by Dinophyceae. In the case of 28S dinoflagellate primers (28S_D), the highest number of OTUs was expectedly associated with dinoflagellates, followed by Bacillariophyta and Haptophyta. From the OTUs detected by 18S universal primers, 29 and 20 OTUs were also found by using 28S_D or 28S_E primers, respectively, and 8 OTUs were detected by both genes targeted and all eukaryote primers employed (Suppl. material 1, Fig. S5).

In total, 24 and 29 unique OTUs associated with micro-and mesozooplankton were detected. The majority of those were identified only based on the 18S rRNA gene (18 and 25 OTUs, respectively), whereas the 28S_D marker facilitated the detection of 2 unique mesozooplankton OTUs (Suppl. material 2, Table S5). The majority of microzooplankton OTUs were associated with Ciliophora, whereas one of the detected OTUs represented Picozoa. The most diverse genera were Strombidium and Mesodinium that included 5 and 3 species, respectively. From mesozooplankton, OTUs were associated with 9 phyla: Annelida, Arthropoda, Bryozoa, Cnidaria, Ctenophora, Mollusca, Nemertea, Priapulida, and Rotifera (Suppl. material 2, Table S5). The majority of the detected OTUs were associated with Arthropoda (14 OTUs), followed by Annelida, Mollusca, and Rotifera (3 OTUs per phyla). In addition, HTS-approach facilitated the detection of 8 and 7 micro-and mesozooplankton species/genera previously not reported from the GoF (Suppl. material 2, Table S5).

From phytoplankton species, 15 OTUs associated with common or dominant species (Suikkanen et al. 2007; Toming and Jaanus 2007; Hällfors et al. 2013; Lips et al. 2014; Jaanus et al. 2017) were detected, with 7 of those only by 18S universal primers at > 0.99 similarity level (Suppl. material 2, Table S6). In addition, three common or dominant species (Chaetoceros wighamii, Dinophysis acuta, Eutreptiella gymnastica) were detected at < 0.99 similarity level. Furthermore, 14 OTUs associated with toxin-producing HAB species (Moestrup et al. 2021) were recorded from the Gulf of Finland (Suppl. material 2, Table S6). OTUs associated with Alexandrium ostenfeldii and Prymnesium polylepis were detected based on all the gene regions targeted (18S and 28S rRNA genes) and primers used, whereas species associated with 7 OTUs (Aphanizomenon flos-aquae; not toxic in the Baltic Sea; Sivonen et al. 1989; Šulčius et al. 2015), Aureococcus anophagefferens, Azadinium dexteroporum, Dinophysis norvegica, Dolichospermum planctonicum, Microcystis aeruginosa, Nodularia spumigena, Prymnesium faveolatum, and Pseudo-nitzschia pungens) were only registered by one of the markers used. The spring-bloom dinoflagellate species complex (DinoComplex), not distinguishable in Lugol-preserved samples by light microscopy, consisting of Apocalathium malmogiense, Biecheleria baltica, and Gymnodinium corollarium (Sundström et al. 2010) were also recorded. However, only G. corollarium sequence was ≥ 0.99% similar to the sequences in the database. Apocalathium malmogiense and B. baltica were detected as a multihit by 18S universal primers, whereas G. corollarium could not be identified based on the 18S rRNA gene. In the case of 28S_E, B. baltica was not recorded and A. malmogiense was detected at a < 0.99 similarity level, whereas sequences associated with G. corollarium were present in several samples. The universal primers mainly targeting dinoflagellates (28S_D) could detect G. corollarium at ≥ 0.99 % similarity level, whereas B. baltica was present at < 0.99 similarity level, and A. malmogiense was not recorded. Another important spring bloom dinoflagellate Peridiniella catenata, was identified both based on 18S and 28S_D primers, but not based on the 28S_E. Skeletonema marinoi, a common diatom species (Hällfors et al. 2013; there as S. costatum), was only detected at ≥ 0.99 similarity level by 18S, but at < 0.99 similarity level by 28S_E and 28S_D primers.

Based on 16S universal primers, the presence of 893 prokaryotic OTUs (> 50 sequences) from 22 phyla and 63 order/class was detected (Table 2). Around 90% of the bacterial sequences could be assigned to a specific phylum and around 82% to a class/order (Suppl. material 2, Table S7). On the phylum level, the highest number of OTUs was associated with Proteobacteria as well as with the Fibrobacteres, Chlorobi, and Bacteroidetes (FCB) group (Fig. 6). On the class/order level, the highest number of OTUs belonged to Alphaproteobacteria, Bacteroidetes/Chlorobi group, and Gammaproteobacteria (Table 2). Among the stations, the number of OTUs detected ranged from 23 to 247 for eukaryotes (depending on the marker) and from 326 to 851 for prokaryotes (Suppl. material 2, Table S8). The highest number of eukaryotic and prokaryotic OTUs was detected from station AP5 based on 28S_E, 28S_D, whereas in the case of 18S and 16S, the highest number of OTUs was present at station A8. The lowest number of OTUs was present at station AP8 based on 28S_E, and 16S, but at AP13, based on 18S and 28S_D (Suppl. material 2, Table S8).

3.3.1 Eukaryotes

From all the class/phyla detected, Bacillariophyta, Ciliophora, Chlorophyta, Dinophyceae, and Picozoa were present during all the sampling months based on different eukaryotic markers. Based on 18S and 28S_E, the relative sequence abundances of Dinophyceae generally increased from February with the highest values from March to May (Figs 3, S6). Bacillariophyta had higher relative sequence abundances in April, whereas Ciliophora had higher relative sequence abundances from December to February/March as well as in May and June (Figs 5, S6). Chlorophyta displayed high relative sequence abundances in spring (March, April), but also in June, August, and October (18S, 28SE; Figs 3, S6). The highest relative sequence abundances for Picozoa were detected in May or in June based on both markers. The appearance patterns of phyla not detected from all sampling occasions also varied depending on the marker used. From all the OTUs detected by 18S and 28S_D <10 % were unique among the stations (winter-spring and summer-autumn stations compared separately), whereas from the OTUs detected based on 28S_E up to 34 % could be unique, depending on the station (Suppl. material 2, Table S8). Among the winter-spring and summer-autumn stations >50% of the OTUs were shared in the case of all 3 markers (Suppl. material 1, Fig. S4).

Figure 3.

Relative sequence abundances of different phyla detected by 18S in different sampling months. TOP 5 most abundant phyla/groups are shown for each station and sampling occasion (X-axis). Y-axis displays the relative sequence abundances for each phyla/group. Sequences belonging to the phyla/group that were not among the TOP 5 phyla/group for the particular sampling occasion are displayed as “Others”.

From 53 unique micro-and mesozooplankton OTUs detected by different markers, five were present in all sampling occasions and stations. Those OTUs were associated with ciliates (Askenasia sp., M. major, M. rubrum, S. lemnae) and a Picozoa, P. judraskeda (Fig. 5). Six unique OTUs were only detected once throughout the sampling season (Arthropoda: Mesocyclops pehpeiensis, Metridia gerlachei; Bryozoa: Einhornia crustulenta; Ciliophora: Epicarchesium pectinatum, Linostomella sp., Vorticella aequilata). On average, the highest number of unique OTUs were detected in July and August based on the 18S and 28S_D markers, but in January and October based on the 28S_E marker (data not shown). From the OTUs identified by multiple markers: Askenasia sp., Mesodinium rubrum, and Picomonas judraskeda were detected from more samples based on the universal 18S marker compared to 28S_E or 28S_D. An exception to this was Gonothyraea loveni, which was detected from more samples by 28S_E (n=3) than by the 18S (n=1) marker. A ciliate, Stylonychia lemnae, was detected only based on markers targeting the 28S rRNA gene (28S_E: n=39, 28S_D: n=38). In general, the detected zooplankton species did not display distinct appearance patterns, however, Acartia bifilosa and Mertensia ovum had higher relative sequence abundances from December to early February (Fig. 5).

From the 31 common/dominant/toxin-producing phytoplankton species detected, about 6 displayed a distinct seasonal pattern by being present mainly in spring or summer, e.g. Chaetoceros holsaticus, Heterocapsa triquetra, G. corollarium, N. spumigena, Pauliella taeniata, and Protoceratium reticulatum (Fig. 4). Some OTUs displayed higher relative sequence abundances in both spring and late summer, e.g. Dinophysis acuminata, D. norvegica, Karlodinium veneficum, and Phalacroma rotundatum. At the same time OTUs associated with Azadinium dexteroporum, and Plagioselmis prolonga were present during all the sampling occasions and at all stations without notable changes in their relative sequence abundances. A photosynthetic ciliate, Mesodinium rubrum, was detected from all samples, however, its relative sequence abundances were higher in spring, but occasionally also during summer (Fig. 4). The OTUs associated with the species reported for the first time from the Baltic Sea or GoF: A. dexteroporum, D. planctonicum, K. veneficum, and P. faveolatum were detected during several months, whereas A. anophagefferens and P. pungens were detected only in a few months throughout the year (Fig. 4).

Figure 4.

Relative sequence abundances of common/dominant/toxin-producing species detected based on different markers. The X-axis represents stations and sampling months, Y-axis displays species detected. Species marked with bold are known to produce toxins (Moestrup, et al. 2020). Station and month labels marked with “_p” indicate months with pooled samples. Different colors indicate associations with different phyla. The size of the bubbles indicates relative sequence abundances. Vertical lines differentiate between the sampling months. *no toxicity confirmed from the Baltic Sea (Sivonen et al. 1989; Šulčius et al. 2015)

Figure 5.

Relative sequence abundances of zooplankton species detected based on different markers. The X-axis represents stations and sampling months, Y-axis displays species detected. Station and month labels marked with “_p” indicate months with pooled samples. Different colors indicate associations with different phyla. The size of the bubbles indicates relative sequence abundances. Vertical lines differentiate between the sampling months.

3.3.2 Prokaryotes

OTUs associated with the FCB group and Proteobacteria had high relative sequence abundances in the majority of the stations and sampling occasions (Fig. 6). In May, the FCB group became more dominant at several stations, whereas Actinobacteria relative sequence abundances were higher in June (Fig. 6). Cyanobacteria became more abundant from April with higher relative sequence abundances in June and July (Fig. 6). The majority of the bacteria that contributed more than 0.5% of the entire bacterial dataset (top 39 OTUs) based on the relative sequence abundances were detected throughout the sampling months (Fig. 7). The most abundant OTUs belonged to bacterial classes Actinobacteria, Alphaproteobacteria, Bacteroidetes, Gammaproteobacteria, Verrucumicrobia (Fig. 7). Three members of Actinobacteria were prevalent throughout the year. Some phyla displayed general appearance patterns. For example, Actinobacteria had the highest relative sequence abundances in February, April as well as from June to August. The relative sequence abundances of unicellular picocyanobacteria (Synechococcus and Synechocystis) were highest during the summer months, from the beginning of June to the end of August, with a small peak in April. Some OTUs were mainly present during winter/cold months, like Roseovarius, Nitrospina, and unclassified Bacteroidetes. Several phyla had the highest relative abundances after phytoplankton spring bloom (April, May), e.g. Alphaproteobacteria (Rhodobacter), Bacteroidetes, and Gammaproteobacteria (Fig. 7).

Figure 6.

Relative sequence abundances of different phyla detected by 16S in different sampling months. TOP 5 most abundant phyla/groups are shown for each station and sampling occasion. (X-axis). Y-axis displays the relative sequence abundances for each phyla/group. Sequences belonging to the phyla/group that did not belong to the TOP 5 phyla/group for the particular sampling occasion are displayed as “Others”.

Figure 7.

Relative abundances of bacteria-associated OTUs (TOP 39). The X-axis represents stations and sampling months, Y-axis displays selected taxa. Order of taxa is derived from pair-wise similarity matrix based on r-values (Pearson). Station and month labels marked with “_p” indicate months with pooled samples. Different colors indicate associations with different phyla. The size of the bubbles indicates relative sequence abundances. Vertical lines differentiate between the sampling months.

Based on 18S, 28S_D and 16S markers, samples from winter-spring and autumn (December-April and October) clustered together on the NMDS plots (Fig. 8). Samples from May/June had a wider distribution, whereas samples from July and August formed a separate group. No clear pattern in sample distribution could be detected based on the 28S_E marker. The biodiversity in sampled communities was significantly different between months (P < 0.01) based on all the markers (Table S9). However, the assumption on the homogeneity of multivariate dispersions was not met for any of the datasets. Thus, it is not certain if the significant difference in biodiversity between the months reflects true dissimilarity between the months or differences in biodiversity within the months compared (Anderson et al. 2006).

Figure 8.

A plot of NMDS analysis of diversity patterns in 38 samples of plankton collected from 9 stations in the Gulf of Finland, the Baltic Sea, A: 18S, B: 28S_E, C: 28S_D: 16S (k=3, stress 18S: 0.07; 28S_E: 0.15; 28S_D: 0.08; 16S: 0.06).

Common/dominant or toxin-producing species displayed statistically significant (P < 0.05) correlations with 24 zooplankton species (Suppl. material 1, Fig. S7). Aureococcus anophagefferens, A. dexteroporum, A. flos-aquae, and D. acuminata had the highest number of significant correlations (n=8) with zooplankton species, with the majority of those being ciliates. An exception to this was A. anophagefferens that had significant correlations with zooplankton species belonging to various phyla. The majority of the significant correlations between phyto-and zooplankton were positive. Two phytoplankton species, Pauliella taeniata and Protoperidinium bipes, did not have any statistically significant correlations with zooplankton. From zooplankton, ciliates Strombidium conicum, S. biarmatum, Mesodinium sp., and a ciliate resting stage Hexasterias problematica together with a comb jelly, Mertensia ovum, had the highest number of statistically significant associations (n= 5 to 8) with common, dominant, or toxin-producing phytoplankton species (Suppl. material 1, Fig. S7). All common/dominant or toxin-producing species also displayed statistically significant correlations with various bacteria (Suppl. material 1, Fig. S9). For phytoplankton, the number of statistically significant correlations ranged from 5 in the case of Coscinodiscus granii to 208 in Plagioselmis prolonga. In the case of bacteria, the number of significant correlations ranged from 1 to 6. Five bacterial OTUs displaying the highest number of statistically significant correlations were associated with Alphaproteobacteria, Micrococcales, and Cyanobacteria. On the phylum level, Proteobacteria and FCB group had the highest number of correlations with common/dominant or toxin-producing species (data not shown). From the common, dominant, or toxin-producing species the highest number of statistically significant correlations were detected for A. anophagefferens (169), Plagioselmis prolonga (208), and Protoperidinium pellucidum (106).

4.1.1 Eukaryotes

In this study, the highest number of eukaryotic OTUs detected belonged to Dinophyceae (phytoplankton). An exception to this was the 28S_E marker, which detected the highest number of OTUs associated with Ascomycota (fungi). A high number of OTUs associated with Dinophyceae has also been reported in previous studies from the Baltic Sea as well as from other localities around the world (de Vargas et al. 2015; Nagai, et al. 2016b; Dzhembekova et al. 2017; Gran-Stadniczeñko et al. 2018; Moreno-Pino et al. 2018; Zhang et al. 2020). This has been explained by the higher rRNA gene copy numbers in Dinophyceae compared to some other groups, which can increase the number of sequences detected from this group (Zhu et al. 2005; Gran-Stadniczeñko et al. 2018). Based on the 28S_E primers, two fungal phyla (Ascomycota and Basidiomycota), were associated with the TOP2 highest numbers of OTUs detected. Those phyla were also detected by the 18S rRNA gene. The detection of more fungi-associated OTUs by the 28S_E marker can be explained by its high species detection power in fungi (Schoch et al. 2012). In a previous study on planktonic fungi in the Baltic Sea, a higher dominance of Chytridomycota was reported at salinities < 8, whereas Ascomycota and Basidiomycota were more dominant at higher salinities (Rojas-Jimenez et al. 2019). Despite this general pattern, other phyla (e.g. Ascomycota) dominated at some stations even at salinities of < 8 (Rojas-Jimenez et al. 2019). In addition to salinity, the dominance of different phyla might be further influenced by other parameters, e.g. nutrient concentrations, and phytoplankton biomass (Wang et al. 2018), which might explain the differences between the results of the present and previous studies.

4.1.1.1 Zooplankton

Based on the HTS-approach, 8 species of micro- and 7 species of mesozooplankton previously not recorded from the Baltic Sea were detected. Based on the zooplankton monitoring data from 2013/2014 (Estonian Env. Agency., 2021a), only 7 mesozooplankton species were detected by the traditional morphology-based approach compared to the 29 species detected by the HTS-approach. Microzooplankton is not a part of routine zooplankton monitoring, except for M. rubrum that is included in the phytoplankton analysis. The importance of including microzooplankton in the monitoring program has been recently emphasized due to their notable contribution (> 50 %) to the nano- and microplankton biomass (Lipsewers and Spilling 2018). To estimate biomass, microscopy-based analysis is necessary, however, to improve species-level detection, molecular methods, such as HTS-approach, have been suggested (Lipsewers and Spilling 2018).

In the case of mesozooplankton, improved detection of zooplankton species based on HTS-approach compared to the morphology-based detection has been reported before (Lindeque, et al. 2013; Hirai, et al. 2015; 2017, Djurhuus, et al. 2018; Meredith et al. 2021). Especially, eDNA samples in comparison to bulk community samples have shown increased diversity detection, as those are not influenced by the mesh size, net avoidance, and destruction of fragile organisms, e.g. jellyfish (Djurhuus et al. 2018). In this study, two of the five newly detected species from the Baltic Sea were benthic species: Gonothyraea loveni and Halicryptus spinulosus. Increased detection of larvae of benthic organisms compared to morphology-based observations has also been reported in other studies utilizing HTS (Lindeque et al. 2013; Hirai et al. 2017; Djurhuus et al. 2018). A potential explanation for this may be the difficulty of identifying the larvae of benthic organisms based on morphology (Lindeque et al. 2013; Hirai et al. 2017). In this study, larvae of two non-indigenous Marenzelleria species were also detected. Interestingly, only one of them is reported from the GoF (M. arctia; Blank et al. 2008; Maximov 2015). However, the other species M. viridis has been reported to be present until 2005 (Estonian Env. Agency 2021 b), after which it was reconsidered to be M. neglecta due to new information on salinity preferences (Bastrop and Blank 2006). However, based on the results obtained in this study, M. viridis might still be the species present in the GoF instead of M. neglecta. To confirm the presence of those non-indigenous species, further molecular study with increased sampling effort in the GoF is needed as already conducted for several other basins of the Baltic Sea (Bastrop and Blank 2006; Blank et al. 2008).

The HTS-approach also detected three of the most abundant mesozooplankton species in 2013/2014 based on the monitoring reports (Eurytemora affinis, Keratella quadrata, Limnocalanus macrurus; Estonian Env. Agency 2021a). However, some dominant species (Eubosmina maritima, Synchaeta baltica, Pseudocalanus acuspes, and Pleopsis polyphenoides) were not detected by any of the three eukaryotic markers used. This might be due to the lack (S. baltica) or low number (max. 6) of 18S rRNA gene sequences available in the database (NCBI 2021). At the same time, more sequences were available for E. maritima, P. acuspes, and P. polyphenoides based on the mitochondrially encoded cytochrome c oxidase I (COI) gene (NCBI 2021). Thus, future studies targeting zooplankton by molecular tools should also include this marker (Machida et al. 2009; Bucklin et al. 2010; Zaiko et al. 2015). However, it should be noted that COI may detect fewer OTUs than the 18S rRNA gene (Djurhuus et al. 2018; Sawaya et al. 2019). Also, to obtain species-level identification group-specific primers may be needed when targeting COI, potentially resulting in inconsistent amplification success (Bucklin, et al. 2016). Interestingly, no representatives of Cladocera were detected by the three eukaryotic markers used in this study, although according to the morphology-based monitoring data from 2013/2014 cladocerans were present (Estonian Env. Agency 2021b). However, this order has been previously detected based on the 18S rRNA marker (Hirai, et al. 2017). Thus, the lack of detection may be influenced by the sampling method as surface water samples were used in this study compared to the net samples taken from the bottom of the water column to the surface by Hirai and others (2017). As cladocerans perform daily vertical migration (Boch and Taylor 1973; Wong et al. 2008), the usage of surface water samples might have hindered their detection in the present study.

4.1.1.2 Common/dominant/toxin-producing species

From the OTUs associated with toxin-producing species (Moestrup et al. 2021), two species, Azadinium dexteroporum, and Prymnesium faveolatum, have previously not been reported from the Baltic Sea (Hällfors 2004). Recently P. faveolatum was also detected from Skagerrak by the HTS-approach, but not by light microscopy (Gran-Stadniczeñko et al. 2018). In the routine phytoplankton monitoring, Prymnesium is only identified to a genus level (Olenina et al. 2006), which may explain the lack of previous records from the Baltic Sea. For A. dexteroporum, potential misidentification as Heterocapsa rotundata that has a similar shape and size (A. dexteroporum: average length 8.5 µm, width 6.2 µm; Percopo et al. 2013; H. rotundata: length 10–15 µm, width 5–10 µm; Olenina et al. 2006; Iwataki 2008) might explain why it has not been recorded before. Aurococcus anophagefferns, another toxin-producing species recorded in this study, is not reliably identifiable by light microscopy due to its small size (around 3 µm) and lack of morphological characteristics (Wang and Lei 2016). However, the species has been previously reported from the GoF and Kattegat based on HTS-approach (Hu et al. 2016; Majaneva et al. 2017), but not by light microscopy (Hällfors 2004).

Other new records of toxin-producing phytoplankton species were Dolichospermum planctonicum (previously Anabaena planctonica) and Pseudo-nitzschia pungens that have been previously reported from the southern Baltic, but not from the GoF (Hällfors 2004). Sequences of D. planctonicum were present in almost all sampling occasions and stations, potentially indicating that the species is common. The lack of detection of this species by light microscopy might be due to the absence of notable morphological differentiation from other Dolichospermum species present in the GoF (D. circinale, D. flosaquae, D. spiroides; Hällfors 2004; Rajaniemi et al. 2005) due to the high proportion of short filament fragments present in the phytoplankton samples instead of longer filaments. At the same time, Pseudo-nitzschia species have not been previously detected from the GoF or adjacent basins (Hällfors 2004). As Pseudo-nitzschia is identifiable to genus level by light microscopy (Hasle et al. 1996; Lelong et al. 2012), its presence in the GoF should be detectable. Pseudo-nitzschia pungens can grow even at a salinity of 5 (Pednekar et al. 2018), which is similar to the salinity in May when OTUs associated with P. pungens were present. However, based on the available information, the species does not grow at temperatures below 10 ˚C (Korean strain; Cho et al. 2001). The water temperature in May at both stations was around 5 ˚C. In the Kiel Bight (southern Baltic Sea), Pseudo-nitzschia spp. bloom mainly in autumn (Wasmund et al. 2008), whereas at the entrance to the Baltic Sea (Skagerrak) P. pungens is present throughout the year (Hasle et al. 1996). As the sequences were only detected in May, it might also reflect introduction by ballast water. To confirm the presence of P. pungens or Pseudo-nitzschia spp. in the northern Baltic Sea further investigation is needed, involving morphology-based identification by electron microscope and molecular characterization by the internal transcribed spacer (ITS) region, which is a more informative marker than the 28S rRNA gene for this genus (Lim et al. 2018).

From the spring bloom dinoflagellate species, G. corollarium from the DinoComplex and P. catenata, could be detected by the HTS-approach. However, two other members of the DinoComplex, A. malmogiense, and B. baltica were not detected or present at a low similarity level (<0.99%). In the case of the 18S rRNA gene, both species were classified together with other species, e.g. A. malmogiense was identified together with A. aciculiferum and B. baltica together with B. cincta and B. brevisulcata at ≥ 0.99 % similarity level. In the case of A. malmogiense, this can be explained by the lack of differences in the 18S rRNA gene sequences between the two species (Annenkova et al. 2015). However, B. baltica and B. cincta belong to closely related, but different clades based on the 18S rRNA gene (Balzano et al. 2012). Although this indicates the potential for distinguishing between the species based on the 18S rRNA gene, the detection may be limited due to the shorter sequences obtained from HTS and following bioinformatic treatment (369 to 489 bp) compared to around 1500 bp used for the phylogenetic analysis (Balzano et al. 2012). As a future perspective, usage of longer sequences obtained by the third-generation sequencing or synthetic long-read sequencing approaches has the potential to improve species detection from the eDNA samples (Schlaeppi et al. 2016; Singer et al. 2016; Heeger et al. 2018; Tedersoo et al. 2018).

4.1.2 Bacterioplankton

Overall, class-level taxa prominent within the bacterioplankton community composition (BCC) and their seasonal succession, although providing a very general overview, complemented near-perfectly the previous findings from the GoF that have used the same size fractioned filtration system (Laas et al. 2014; 2015; 2016). Together with these previous studies from the GoF, the dynamics of BCC at the surface layer are now well-documented from December of 2010 to May of 2014 and demonstrate repeating interannual patterns of changes in BCC at the surface layer. The dominance of various bacterial groups has also shown positive and negative associations with other parameters, e.g. nutrient concentrations, temperature, Chl α, salinity, turbidity, phytoplankton biomass, diatom: dinoflagellate ratio (Andersson et al. 2010; Laas et al. 2014; 2015; Bunse et al. 2016), but especially with phytoplankton community composition (Camarena-Gómez et al. 2020). This was further demonstrated by the data from this study as BCC followed seasonal changes in phytoplankton community composition.

The highest number of OTUs were associated with Alphaproteobacteria and FCB group (Bacterioidetes/Chlorobi group), followed by Actinobacteria. In general, Cyanobacteria became more prevalent during summer and Flavobacteria during the spring bloom, although some members of the FCB group occur exclusively during summer bloom. Interestingly, the unicellular cyanobacteria, Synechocystis and Synechococcus, increased in relative sequence abundance also in April. Alphaproteobacteria were dominant at the surface layer throughout different seasons, compared to previous studies, SAR11 bacteria (Pelagibacterales) were less prevalent in their relative abundances, but they were still present throughout the year (Laas et al. 2014; 2015; 2016). Rhodobacter followed a similar occurrence pattern as observed during previous years (Laas et al. 2014; 2015; 2016). Studies encompassing other basins of the Baltic Sea have found salinity as a major factor controlling BCC (Herlemann et al. 2011). Actinobacteria and Betaproteobacteria (now reclassified as Betaproteobacterales) are generally more abundant in oligohaline areas (Herlemann et al. 2011; Laas et al. 2015), in the present study many OTUs belonging to these taxa remained even during the months with higher salinity (December-February). However, the elevated occurrence of Betaproteobacterales in May and July was probably caused by the wider spread of freshwater in the surface layer of the GoF (Lips, et al. 2017).

4.2.1 Zooplankton

From the microzooplankton species registered in this and a previous study from the northern Baltic Proper (Johansson et al. 2004), Askenasia sp. was detected from a wider range of months, whereas Strobilidium sp. was present during a lesser number of months based on HTS. This can be explained by the yearly variability, as well as by the differences in the months covered: December to August and October in this study, vs. April to December in the previous study. (Johansson et al. 2004). Detection in a wider time range may also be explained by the more sensitive detection based on DNA, than by morphology (Abad et al. 2016). Interestingly, in the case of mesozooplankton, the OTUs were generally detected within the same months as based on monitoring data from the GoF in 2013/2014, but on a lower number of occasions (Estonian Env. Agency 2021a). Since it was a general pattern, differences in sample collection may serve as a potential explanation. In the case of zooplankton monitoring data, the samples are collected by using a 100 µm net throughout the water column (HELCOM 2021), which may help to detect the organisms present also below the 5 m depth used in this study. Further, the sequences present in low relative abundance may also be masked by sequences from other organisms present in higher total abundance or with higher rRNA gene copy numbers (Schneider et al. 2016; Mäki et al. 2017; Deagle et al. 2018). Thus, to obtain a comprehensive overview of both micro- and mesozooplankton appearance patterns in the GoF, further study applying both HTS and morphology-based approaches and taking into account different sampling strategies (water vs. zooplankton net) is needed.

From the Mesodinium species detected, M. rubrum and M. major displayed almost identical appearance patterns. These two species have been identified based on their morphology and division into distinct clades based on rRNA genes and internal transcribed spacer region (Garcia-Cuetos et al. 2012). Differences in their seasonal abundance have also been reported from the Öresund, Baltic Sea, where M. rubrum is more abundant during summer and early autumn and M. major in winter and early spring (Garcia-Cuetos et al. 2012). In this study, no such differences were observed. As both M. rubrum and M. major have been detected from the GoF previously (Johnson, et al. 2016), further investigation is needed, whether the two species display temporal differences in their relative abundances similar to the Öresund. This information can be useful for providing more detailed information on the appearance of those important primary producers in post-bloom conditions (Lips and Lips 2017). In addition, M. rubrum appearance patterns can also reflect the presence of toxin-producing Dinophysis species (Sjöqvist and Lindholm 2011) as it is the main prey item for D. acuminata and several other Dinophysis species (Park et al. 2006; Nagai et al. 2020b).

4.2.2 Common/dominant/toxin-producing species appearance patterns

Investigation of common/dominant/toxin-producing species appearances patterns based on their relative sequence abundances displayed similar patterns with 20 years of phytoplankton monitoring data (for the species identifiable under light microscopy; Estonian Env. Agency 2021c; Suppl. material 2, Table S6). It is recognized that the relative sequence abundances do not reflect the number of cells/organisms or biomass directly (Egge et al. 2013; Weber and Pawlowski 2013; Lamb et al. 2019). However, it may provide an overall impression of species/classes relative abundances (Eiler et al. 2013; Giner et al. 2016; Gran-Stadniczeñko et al. 2017; Schenk et al. 2019). For example, the relative sequence abundances obtained by HTS-approach have been reported to coincide with the yearly total abundance patterns detected by using light microscopy for the toxin-producing dinoflagellate Alexandrium catenella (Nagai et al. 2019). Furthermore, in combination with light microscopy-based counts, the relative sequence abundances may provide a semi-quantitative estimate of cell abundances, which has high relevance for morphologically similar species, e.g. Pseudo-nitzschia spp. (Nagai et al. 2017).

In this study, the majority of the common/dominant/toxin-producing species detected appeared in various months throughout the year, and 2/3 of those species had the highest relative sequence abundances within the months with the highest biomass reported (Estonian Env. Agency, 2021c). For example, all detected cyanobacterial species had higher relative sequence abundances in summer, which is consistent with the previous knowledge from published studies (Kanoshina et al. 2003; Laamanen and Kuosa 2005) and long-term phytoplankton monitoring dataset (Estonian Env. Agency 2021c). However, 1/3 of the highest relative sequence abundances were recorded in months different from the months with the highest biomass reported in the phytoplankton dataset. As the high relative sequence abundances were within the months of the known appearance for the species and the difference was usually one month earlier or later than previously recorded, it probably reflects variation in the conditions favoring or hindering their growth e.g. nutrient depletion, upwelling (Lips and Lips 2010; Lips et al. 2014). Thus, the detected species generally followed the previously known patterns in their appearance and abundance.

For species that were detected in this study, but not present in the long-term phytoplankton dataset, the appearance patterns were investigated based on literature or based on the genus level information from the phytoplankton dataset. For example, the toxin-producing dinoflagellate A. ostenfeldii was present during almost all the sampling months with the highest relative sequence abundances in July. However, in the Åland archipelago, the cells are present from May to September with the highest abundances in July, August, and September (Kremp et al. 2009; Hakanen et al. 2012; Le Tortorec et al. 2014). The detection of this species in winter and early spring is surprising as the A. ostenfeldii cells have not been detected from the water column in the Åland archipelago at 10˚C (Jerney et al. 2019). Furthermore, the germination of A. ostenfeldii resting cysts is inhibited if the temperature is <10˚C (Jerney et al. 2019). Thus, the early detection in this study might result from the presence of resuspended cysts in the water column. However, further investigation with higher frequency and morphology-based observations is needed to explain the detection of A. ostenfeldii associated sequences during winter and early spring. For another toxin-producing species, Aureococcus anophagefferens, the sequences were detected from July and August, whereas previously it has been recorded in March from the GoF (Majaneva et al. 2017) and in July from Kattegat (Hu et al. 2016). The species can grow from 0 ˚C to 25 ˚C and its blooms have been detected at various temperatures (Gobler et al. 2005 and references therein). For example, in the Sea of Okhotsk, the species is present throughout the year with higher relative sequence abundances during winter (Sildever et al. 2019). Its appearance in various environmental conditions is further supported by the capability to utilize organic carbon and nitrogen (Mulholland et al. 2002). This may explain its presence in the GoF during summer when inorganic nutrients, especially nitrogen, are limiting phytoplankton growth (Graneli et al. 1990; Moisander et al. 2003; Tamminen and Andersen 2007).

Contrary to A. anophagefferens, A. dexeteroporum was detected throughout the year with the highest relative sequence abundances in January, February, and June. Based on the previous studies, the species can also grow in a wide temperature range (0.6 ˚C to 28˚C), depending on the geographic origin of the strain (Percopo et al. 2013; Tillmann et al. 2015, 2020; Cavalcante et al. 2018). Furthermore, based on six years of weekly data analyzed by metabarcoding and HTS-approach, the species was detected throughout the year in the Sea of Okhotsk, from water temperatures of -1.7 ˚C to around 20˚C (Sildever et al. 2019). Thus, the high relative sequence abundances detected in this study during winter and summer, probably reflect the capability of the species to utilize suitable conditions for growth in a broad temperature range. Another species detected throughout the year was Karlodinium veneficum, displaying peaks in relative sequence abundances also in various seasons, e.g. winter: from December to February, and summer: May to July. A similar pattern of continuous presence and peaks in biomass/relative sequence abundances at different months have been also shown in the east coast of the U.S.A. and the Sea of Okhotsk (Zhang et al. 2008; Sildever et al. 2019). At both localities, the species was detected from a wide salinity and temperature range: salinity of 0 to 32 g kg^-1 and temperature of 4 ˚C to 30 ˚C in the east coast of US; salinity of 20 to 33 g kg^-1 and temperature to -1.7 ˚C to 22 ˚C in the Sea of Okhotsk (Zhang et al. 2008; Sildever et al. 2019). Since K. veneficum is mixotrophic, its appearance may reflect the availability of its prey, as its growth rates can be 2 to 3 times higher when grown mixotrophically (Li et al. 1999). However, no pattern between the relative sequence abundances of K. veneficum and its prey (Cryptophyta) was detected in this study (data not shown).

Prymnesium parvum f. patelliferum is a form of P. parvum (Larsen 1999), a cosmopolitan euryhaline (salinities of 3 to 30 g kg^-1), and eurythermal (5 ˚C to 30 ˚C) species (Edvarsen and Paasche 1998; Larsen and Bryant 1998). It has formed toxic blooms around the world (Edvarsen and Paasche 1998; Johnsen et al. 2010; Roelke et al. 2016), as well as in the coastal waters of SW Finland and a brackish water lake on the Åland islands (Lindholm and Virtanen, 1992; Lindholm et al. 1999). Prymnesium parvum is present in several basins of the Baltic Sea, including the GoF (Hällfors 2004). The species has been detected throughout the year with the highest abundances recorded at the end of June and July in the Baltic Sea (Fistarol et al. 2003). The species has not been reported in the long-term monitoring data (Estonian Env. Agency, 2021c), potentially explained by the difficulties with the morphology-based identification due to small size (Green et al. 1982) and the negative influence of sample fixation on cell morphology (Galluzzi et al. 2008). However, the long-term dataset provided information on the appearance of Prymnesiales and Prymnesium: presence throughout the year with the highest biomass detected from May to August (Estonian Env. Agency 2021c). In this study, P. parvum f. patelliferum was recorded from May to August, but also in December, with the highest relative sequence abundances in May, and July-August. The lack of detection before May by the HTS approach might be potentially due to its presence in very low abundances, further influenced by the amount of water filtered (up to 1 L) as the species was present in low abundances also at times when it was detected by the HTS-approach (relative sequence abundances < 0.3%). Another Prymnesium species detected in this study, P. faveolatum, has only been previously recorded from the Mediterranean Sea and the northern Atlantic (Fresnel et al. 2001; Edvardsen et al. 2011) and no detailed information is available on its ecological preferences. Thus, further investigations are needed to confirm its appearance in different geographic locations and environmental conditions.

Based on the correlation analysis between the common, dominant, or toxin-producing phytoplankton and zooplankton species, the statistically significant (P < 0.05) associations probably reflect co-appearance and not interactions as the majority of the species from both groups were detected throughout the sampling period. This is further illustrated by the significant negative correlations between dinoflagellate, A. anophagefferens, and a ctenophore, Mertensia ovum, as the dinoflagellate was only present in the samples from July and August, whereas the ctenophore species was present in the majority of the samples, except during those months. A similar, but opposite pattern was displayed between a chrysophyte Plagioselmis prolonga and several zooplankton species with the cryptophyte missing from the July-August samples. The significant positive correlations may also reflect the zooplankton presence following the presence of its prey, e.g. copepod Limnocalanus macrurus and diatoms Melosira arctica and Skeletonema marinoi (Mäkinen et al. 2017) or ciliates and Heterocapsa triquetra (Stoecker and Evans 1985; Hansen, 1995). However, zooplankton generally avoids grazing on toxin-producing phytoplankton species, e.g. Alexandrium ostenfeldii, A. anophagefferens, Prymnesium faveolatum (Graneli et al. 2003; Gobler et al. 2004; Sopanen et al. 2011), which should lead to negative associations between the taxa if the detected patterns reflect feeding.

In the case of correlations between common, dominant, or toxin-producing phytoplankton species and bacteria, the majority of the significant correlations detected were positive, except for A. anophagefferens. However, as the negative correlations were detected between A. anophagefferens and bacteria belonging to various phyla and classes, the significant negative correlations might reflect differences in appearance patterns. On average, the highest number of bacterial OTUs from all stations were present from December to April (average number of OTUs: 506–603), whereas during the detection of A. anophagefferens in July-August, the average number of bacterial OTUs ranged around 380. The interactions between phytoplankton and bacteria can range from mutualistic to competitive or antagonistic (Cole, 1982; Seymour, et al. 2017) with potential shifts in the type of the relationship (Wang et al. 2010; Seyedsayamdost et al. 2011). The associations detected between phytoplankton and bacteria in this study might thus provide a broad snapshot of the diversity present in the phycosphere of the phytoplankton cells in different species (Johansson et al. 2019; Blifernez-Klassen et al. 2021). This information might provide useful insight for more detailed studies on the bacterial diversity associated with various phytoplankton species.