Diversity of Pico- to Mesoplankton along the 2000 km Salinity Gradient of the Baltic Sea

Hu, Yue O. O.; Karlson, Bengt; Charvet, Sophie; Andersson, Anders F.

doi:10.3389/fmicb.2016.00679

ORIGINAL RESEARCH article

Front. Microbiol., 12 May 2016

Sec. Aquatic Microbiology

Volume 7 - 2016 | https://doi.org/10.3389/fmicb.2016.00679

Diversity of Pico- to Mesoplankton along the 2000 km Salinity Gradient of the Baltic Sea

$\r\nYue O. O. Hu$ Yue O. O. Hu¹

Bengt Karlson²

Sophie Charvet³

Anders F. Andersson¹^*

¹Science for Life Laboratory, Division of Gene Technology, School of Biotechnology, KTH Royal Institute of Technology, Stockholm, Sweden
²Oceanography, Research & Development, Swedish Meteorological and Hydrological Institute, Gothenburg, Sweden
³Leibniz Institute for Baltic Sea Research Warnemünde, Rostock, Germany

Microbial plankton form the productive base of both marine and freshwater ecosystems and are key drivers of global biogeochemical cycles of carbon and nutrients. Plankton diversity is immense with representations from all major phyla within the three domains of life. So far, plankton monitoring has mainly been based on microscopic identification, which has limited sensitivity and reproducibility, not least because of the numerical majority of plankton being unidentifiable under the light microscope. High-throughput sequencing of taxonomic marker genes offers a means to identify taxa inaccessible by traditional methods; thus, recent studies have unveiled an extensive previously unknown diversity of plankton. Here, we conducted ultra-deep Illumina sequencing (average 10⁵ sequences/sample) of rRNA gene amplicons of surface water eukaryotic and bacterial plankton communities sampled in summer along a 2000 km transect following the salinity gradient of the Baltic Sea. Community composition was strongly correlated with salinity for both bacterial and eukaryotic plankton assemblages, highlighting the importance of salinity for structuring the biodiversity within this ecosystem. In contrast, no clear trends in alpha-diversity for bacterial or eukaryotic communities could be detected along the transect. The distribution of major planktonic taxa followed expected patterns as observed in monitoring programs, but groups novel to the Baltic Sea were also identified, such as relatives to the coccolithophore Emiliana huxleyi detected in the northern Baltic Sea. This study provides the first ultra-deep sequencing-based survey on eukaryotic and bacterial plankton biogeography in the Baltic Sea.

Introduction

The Baltic Sea is the world's second largest brackish water body. It is connected to the North Sea by several narrow straits. A horizontal salinity gradient ranging from 2 to 30 traverses the Baltic Sea and parts of the North Sea in a Northeast–Southwest (NE-SW) direction. This gradient is maintained by sporadic inflows of saline deep water in the south, in combination with freshwater input from rivers. A strong halocline at 60–80 m depth restricts vertical mixing in the Baltic Proper, and large areas of the deeper waters are hypoxic or anoxic (Carstensen et al., 2014). The intermediate salinity and subarctic conditions constrain the macrobial biodiversity in the Baltic Sea, and its multicellular species are fewer and typically genetically impoverished compared to the populations of the NE Atlantic (Johannesson and André, 2006). In contrast, bacterioplankton diversity does not appear to decrease at intermediate salinities (Herlemann et al., 2011), and for protists it has even been suggested that species richness is actually highest in the horohalinicum region (salinity 5–8; Telesh et al., 2011a, 2013).

Unicellular plankton play key roles in global biogeochemical cycles and form the basis of marine food webs (Falkowski et al., 2008). In the Baltic Sea, diatoms and dinoflagellates (>10 μm) dominate primary production in spring, while smaller eukaryotic phytoplankton and cyanobacteria (<3 μm) dominate primary production over the rest of the year (Johansson et al., 2004). Heterotrophic bacteria play key roles in taking up dissolved organic matter stemming from primary production within the ecosystem and from terrestrial runoff-water (Azam and Malfatti, 2007). Bacteria also mediate most transformations in the cycling of nitrogen, phosphorus, trace metals, and other nutrients (Falkowski et al., 2008). Energy and nutrients are propagated to higher trophic levels through grazing. Large phytoplankton are primarily grazed upon by mesozooplankton (e.g., copepods) but also by micro-zooplankton, i.e., heterotrophic protists. The current view is that picoplankton (bacteria and picoeukaryotes; 0.2–2 μm) are predated upon by nanoplankton (nanoflagellates and small ciliates; 2–20 μm), while the nanoplankton are eaten by larger zooplankton (copepods, rotifers, larger ciliates and dinoflagellates) (Samuelsson et al., 2006). The grazing on heterotrophic bacteria brings the organic matter leaked by phytoplankton during photosynthesis back to the traditional food chain, closing what is called the microbial loop (Azam et al., 1983). Finally, the larger zooplankton are predated upon by fish; in the pelagic Baltic Sea, mainly by sprat and herring. In this way, all elements of the pelagic ecosystem are interconnected, leading to cascade effects between different levels of the food web. For instance loss of top predators could be propagated all the way down to affect picoplankton community composition and function (Casini et al., 2009).

Monitoring of phyto- and zooplankton serves to estimate the environmental state and to detect changes regarding eutrophication, biodiversity, harmful algal blooms etc. This is today largely based on microscopic techniques (Hällfors, 2004; Olenina et al., 2006), which are time-consuming (~4 h/sample) and rely heavily on the skills of taxonomists. Cross-comparison between datasets analyzed by different taxonomists may therefore be problematic. The smallest plankton organisms, i.e., the pico- and nanoplankton (0.2–2 and 2–20 μm), and larger organisms that lack distinctive morphological features are difficult to identify using light microscopy. Therefore, these most abundant organisms are most often simply termed “Unicells” or “Flagellates” since they cannot be taxonomically identified. Hence, there is a need for faster, more robust novel methods that can cover a larger fraction of planktonic diversity.

Due to the recent developments of high-throughput sequencing techniques, genetic methods are a promising avenue for monitoring biological diversity. Genetic barcoding consists of taxonomically assigning a specimen based on sequencing a short standardized DNA fragment (barcode) and by matching this to a reference library of sequences of known taxonomy. By metabarcoding, the approach is extended to a community of individuals of different species. Metabarcoding has been extensively used for assessing both prokaryotic and eukaryotic plankton diversity and biogeography (Sogin et al., 2006; Amaral-Zettler et al., 2009; de Vargas et al., 2015). A number of metabarcoding studies have investigated diversity and spatiotemporal distribution patterns of prokaryotes in the Baltic Sea (Andersson et al., 2010; Herlemann et al., 2011; Dupont et al., 2014; Lindh et al., 2015), however, the methodology has not yet been used to map eukaryotic plankton in this sea.

Here, we have conducted ultra-deep Illumina sequencing of surface-water bacterial and eukaryotic plankton communities sampled in the summer of 2013 along a 2000 km transect following the salinity gradient of the Baltic Sea. The aim of the study is to give a comprehensive view of plankton diversity across the salinity gradient, and to investigate how microbial alpha- and beta-diversity is structured by salinity. To our knowledge, this is the first high-throughput sequencing-based investigation of eukaryotic plankton biogeography in the Baltic Sea and we combined it with a renewed description of the bacterial diversity highlighted by previous studies (Herlemann et al., 2011; Dupont et al., 2014).

Materials and Methods

Sampling

Twenty-one water samples were collected in the Kattegat, the Baltic Proper and the Gulf of Bothnia using a FerryBox system installed in the ship TransPaper during 13th–19th of July 2013. The ship followed the route: Gothenburg (Sweden)—Kemi (Finland)—Oulu (Finland)—Lübeck (Germany)—Gothenburg. The FerryBox system consists of a pump with a water inlet at 3 m depth, a circuit of multiple sensors for temperature, conductivity, chlorophyll and phycocyanin fluorescence, turbidity, and oxygen as well as automated water sampling devices. A detailed description of the FerryBox system is found in Karlson et al. (in press). Manual water sampling for DNA analysis was carried out both on the Northward and Southward legs. Approximately, 10 L of seawater were collected in a polycarbonate carboy. Subsamples of 200–500 mL were filtered onto 0.22 μm pore-size mixed cellulose ester membrane filters (Merck Millipore co., Cat. No. GSWP04700) to capture plankton. The filters were frozen in liquid nitrogen on board and kept at −20 to −80°C until DNA extraction.

DNA Extraction, PCR Amplification and Sequencing

Genomic DNA was extracted using the PowerWater^® DNA isolation kit (MO-BIO Laboratories Inc, Carlsbad CA, USA) following the instructions provided by the manufacturer. The V3-V4 regions of bacterial 16S rRNA genes were PCR amplified with primers 341F (CCTACGGGNGGCWGCAG) and 805R (GACTACHVGGGTATCTAATCC) (Herlemann et al., 2011), and the V4–V5 regions of eukaryotic 18S rRNA genes were amplified with primers 574*F (CGGTAAYTCCAGCTCYV and 1132R, CCGTCAATTHCTTYAAR; Hugerth et al., 2014a), using KAPA HiFi HotStart ReadyMix (2X) (KAPA Biosystems, Kit Code KK2602). A two step PCR procedure was applied (Hugerth et al., 2014a), with 35 (25 + 10) and 38 (28 + 10) PCR cycles for 16 and 18S, respectively. Between the first and second PCR, and prior to pooling libraries, amplicons were purified with 8.8% PEG 6000 (Polyethylene Glycol 6000) (Merck Millipore co., Cat. No. 528877-1KG) precipitation buffer and CA beads (carboxylic acid-coated superparamagnetic beads) (Dynabeads^® MyOne™ Carboxylic Acid, Cat. No. 65012; Lundin et al., 2010). Agilent 2100 Bioanalyzer (Agilent, Technologies, DNA 1000 LabChip kit) and Qubit^® 2.0 Fluorometer (Invitrogen, Qubit-IT™ dsDNA HS Assay kit) were used for checking the amplicon fragment sizes and quantification. Equimolar amounts of indexed samples were mixed and sequenced with Illumina MiSeq (Illumina Inc, USA) at NGI/Scilifelab Stockholm. The sequencing reads have been submitted to the European Nucleotide Archive (ENA) under accession numbers PRJEB12362.

Sequence Processing

The OTU table building procedures followed https://github.com/EnvGen/Tutorials, using USEARCH for quality trimming and OTU-clustering (Edgar, 2010), but with some optimizations. The primers were trimmed from forward and reverse reads with FastX (Pearson et al., 1997), as well as 40 bases in the 3′ end of the reverse reads (these displayed lower quality than the forward reads). This resulted in average Phred scores (Q-values) of every position in the trimmed reads to be >20. In order to avoid generating artificial OTUs while assigning as many reads as possible, we used more stringent criteria when generating OTU-centroid sequences than in the mapping (assigning step). For the centroid generation step, 3′-ends were further trimmed such that all remaining bases had Q ≥ 25. Forward and reverse reads were merged, and merged pairs that were shorter than 300 bp, had >3 mismatches, or overlapped by <100 bp were discarded. Remaining merged reads were clustered into OTUs using 99% identity level. For processing reads for the assigning step, the same procedure was applied, but with Q20 as the quality trimming parameter. The trimmed reads were mapped to the OTU centroids to build the 16S OTU table. 16S OTUs (their centroid sequences) were classified with SINA-1.2.13 using the SILVA 119 SSU database (Altschul et al., 1997; Quast et al., 2013; Decelle et al., 2015).

The 18S OTU table was built similarly to the 16S OTU table, but the overlap between read pairs was not sufficient for merging, so each pair was concatenated rather than merged. Primers were trimmed from forward and reverse reads, as well as 44 bases in the 3′ end of the reverse reads, resulting in average Phred scores of every position to be >20. After trimming, forward and reverse reads were concatenated prior to OTU clustering. For the taxonomic classification, OTU centroid sequences were split at the concatenation position, and the forward and reverse counterparts were separately BLAST (Altschul et al., 1997) searched against the PR2 database (based on GenBank 203–October 2014; Guillou et al., 2013). BLAST results were filtered such that only alignments of at least 250 and 210 bp in length (90% of the length of reads), for the forward and reverse read, respectively, were kept. Then, the highest taxonomic level was found where the forward and reverse read aligned to the same database sequence, with alignment identities exceeding the taxonomic level-specific cutoff. If just one database sequence was hit, the OTU inherited the taxonomy of this sequence. If multiple database sequences were hit, the bit-scores of the forward and reverse reads to each hit were summed, and the hits having at least 95% of the best sum of bit scores were selected, and the most detailed consensus taxonomy was built from these hits and assigned to the OTU. The taxonomic-level specific identity cutoffs were as follows: kingdom to family level (PR2 level 1–6): 90%; genus level (PR2 level 7): 97%, species level (PR2 level 8): 99%. These cutoffs were determined by running an in-silico experiment where artificial reads were generated from 1000 randomly selected PR2 sequences, and classified against the remaining sequences in PR2, using a range of identity cutoffs at each taxonomic level. The selected cutoffs were good compromises between classification accuracy (precision) and percentage of sequences classified (sensitivity), and these were for the simulation as follows: kingdom to family level: precision >95%, sensitivity >85%; genus level: precision 78%, sensitivity 42%; species level: precision 42%, sensitivity 38%. In addition to this automatic taxonomic annotation procedure, for a subset of the individual OTUs (n = 30) that were of special interest (e.g., being of high relative abundance) but could not be automatically annotated to a satisfactory detailed level, manual annotations by online BLAST searches against the NCBI nt database were performed, in which case only hits with coverage and identity >99% where considered, and among these the most detailed consensus taxonomy was assigned to the OTU.

Statistical Analysis

All statistical analysis and plotting were conducted in R (www.r-project.org) using the R libraries: vegan (α-diversity; β-diversity; subsampling; Mantel test), rworldmap (geographic plotting), ape (principal coordinates analysis), cluster (hierarchical clustering), and pheatmap (heatmap). For calculating the p-value of the observed level of matching between the taxonomies of correlated 16 and 18S OTUs, we randomly shuffled the taxonomic annotations among the 18S OTUs that displayed significant correlations to 16S OTUs, and compared the taxonomic annotations between the chloroplast 16S OTUs and the 18S OTUs. The shuffling was redone 10⁶ (N) times and we counted how many times (n) 5 or more of the 6 chloroplast 16S OTUs were annotated to the same taxonomic class as the paired 18S OTU (i.e., the same or higher level of matching as the real data). The p-value was calculated as n/N.

Results

Trends in Beta-Diversity

On average 32,948 amplicons of bacterial (16S) and eukaryotic (18S) rRNA genes were sequenced from surface water samples from 21 stations within the Baltic Sea and Kattegat, collected in July 2013 (Figure 1A; stations are described in Supplementary Table 1 and seasonal variation in temperature and salinity of the sampled water bodies presented in Supplementary Figure 1). The sequences were clustered into 1243 and 1860 16 and 18S operational taxonomic units (OTUs), respectively (Supplementary Tables 2, 3). To investigate overall trends in community composition, Principal Coordinates Analysis (PCoA) was conducted on the bacterial and eukaryotic communities (Figure 1). Similar to what has been shown before in the Baltic Sea (Herlemann et al., 2011; Dupont et al., 2014), bacterial community composition changed gradually along the salinity gradient, with the first principal coordinate being highly correlated with salinity (Spearman rho = 0.99, P < 10⁻⁵; Figure 1B). Here we could also, for the first time, show a similar pattern for eukaryotic plankton in the Baltic (Spearman rho = 0.88, P < 10⁻⁵; Figure 1C). Likewise, hierarchical clustering grouped the samples according to their salinity, based on both bacterial and eukaryotic community composition (Supplementary Figure 2). Difference in salinity between samples was also directly correlated with difference in community composition (beta-diversity), for both bacterial and eukaryotic communities (Supplementary Figure 3). Difference in temperature, that to some extent covaries with salinity (Supplementary Table 1), was also correlated with beta-diversity, but the correlation was weaker (Supplementary Figure 3). Since both eukaryotic and bacterial communities correlated with salinity, it followed that beta-diversity of the two community types were correlated (Mantel test, Spearman rho = 0.86, P < 10⁻⁶; Figure 1D).

FIGURE 1

Figure 1. Trends in bacterial and eukaryotic beta-diversity. (A) Map showing sampling sites within Kattegat and the Baltic Sea, with salinity indicated by color. (B–C) Ordination of bacterial (B) and eukaryotic (C) communities using Principal Coordinate Analysis (PCoA) based on beta-diversity calculated using Spearman correlation of OTU frequencies. Samples are colored and sized according to (A). Variation explained by the principal components (PCs) are given within parenthesis. PC1 was highly correlated with salinity for both bacterial and eukaryotic communities (Spearman ρ > 0.97 and P < 10⁻⁵ for both). (D) Pairwise dissimilarity score (beta-diversity) between all samples in the dataset, where the eukaryotic scores (y-axis) are plotted against the bacterial scores (x-axis) from the same samples.

Bacterial Community Composition

Our analysis confirmed the overall trends in bacterial taxonomic composition observed in previous studies, with an increase in Alpha- and Gammaproteobacteria and a decrease in Actinobacteria and Betaproteobacteria with increasing salinity levels (Herlemann et al., 2011; Dupont et al., 2014). Similar to what was found by Herlemann et al. (2011), Verrucomicrobia peaked in intermediate salinity levels, which here was also the case for Planctomycetes and Cyanobacteria (Figure 3A).

Of the 1243 16S OTUs observed in the 21 samples, 1187 could be taxonomically classified at least to the class level. The majority (59%) of the OTU centroid sequences matched SILVA database reference sequences with >99% identity. Identity levels however varied substantially among taxa (Figure 2). Chlamydiae and the newly discovered non-photosynthetic sister clade to Cyanobacteria, the Melainabacteria (Di Rienzi et al., 2013), displayed notably low identity levels (Figure 2A).

FIGURE 2

Figure 2. Identity levels (%) between OTU centroid sequences and their best matching database sequence. (A) Bacterial OTUs. (B) Eukaryotic OTUs. Number of OTUs included in each taxonomic group is given within parenthesis; only groups with ≥10 OTUs are shown.

The gradual shift in community composition observed at the phylum/class level was also evident within specific taxonomic groups, as exemplified by the Cyanobacteria and the Alphaproteobacteria cluster SAR11 in Figures 3B,C. The Cyanobacteria were represented by 26 OTUs, of which 12 had a mean abundance of >0.1% of reads, most of which were related to the picocyanobacteria Synechococcus (Figure 3B). OTUs associated to filamentous cyanobacteria (Aphanizomenon, Nodularia, and Dolichospermum (syn. Anabaena)) were also detected from the brackish region but in much lower abundances than the dominant Synechococcus, consistent with cyanobacterial cell numbers and biomass being dominated by picocyanobacteria in the Baltic Sea (Stal et al., 2003). Phycocyanin fluorescence, as recorded by the FerryBox system, correlated well with relative abundance of sequenced filamentous cyanobacteria (Spearman rho = 0.69, P < 10⁻³), in agreement with Seppälä et al. that observed a linear relation between biomass of filamentous cyanobacteria and phycocyanin fluorescence (Seppälä et al., 2007). The different Synechococcus populations showed clear differences in salinity preference, with OTU_25 and OTU_34 being dominant in Kattegat at more than 99.7% of cyanobacterial reads, while OTU_3 and OTU_4 dominated the Baltic Proper with on average 88% of cyanobacterial reads. OTU_3 was identical to the 16S sequence of a metagenome-assembled genome reconstructed from the Baltic Proper (BACL30; Hugerth et al., 2015) that possess the type IIB phycobilisome gene organization recently identified as being dominant in Baltic Sea picocyanobacteria (Larsson et al., 2014).

FIGURE 3

Figure 3. Abundance distribution of major bacterial groups across the Baltic Sea salinity gradient. (A) Major bacterial taxonomic groups (>0.1% mean abundance). (B) Dominant Cyanobacteria OTUs (>0.1% mean abundance). In Subsection I (picocyanobacteria), OTU_56 is a Microcystis (a freshwater species), remaining OTUs belong to Synechococcus. In Subsection IV (filamentous cyanobacteria), OTU_777, OTU_49 and OTU_19 belong to Aphanizomenon, Nodularia, and Dolichospermum (syn. Anabaena), respectively. The black curve indicates phycocyanin fluorescence (in arbitrary units). (C) Dominant SAR11 OTUs (>0.1% mean abundance).

SAR11 is the most abundant bacterioplankton clade of the world's oceans, and is also found in brackish and fresh waters. It consists of several subclades that together span the candidate order Pelagibacterales (Grote et al., 2012). Eight OTUs were classified as SAR11 in our dataset, of which six had >10 reads (Figure 3C), representing five different subclades of SAR11 (Supplementary Figure 4 for a phylogenetic tree). An OTU (OTU_16) identical to Pelagibacter ubique (open ocean subclade Ia) dominated the marine samples but was also present in the brackish samples, while two OTUs belonging to subclade IIb (OTU_5) and IIIa (OTU_1) dominated the brackish region. Subclade IIIa is frequent in brackish waters and has previously been found in the Baltic Sea (Herlemann et al., 2011; Bergen et al., 2014; Dupont et al., 2014), while subclade II is mainly associated with coastal and mesopelagic waters (Carlson et al., 2009; Beman et al., 2011). Interestingly, contrary to previous studies (Bergen et al., 2014; Dupont et al., 2014), we also found an OTU representing the freshwater subclade IIIb (LD12) in the low-salinity Bothnian Bay.

The deep sequencing provided by the Illumina platform allowed clustering of reads with higher phylogenetic resolution than with earlier 454-data (99% identity; compared to 97% in the study by Herlemann et al., 2011) while still maintaining enough reads per OTU to accurately quantify them. This uncovered highly similar OTUs with distinct spatial distributions. Herlemann et al. (2011) found that the most abundant OTU was a Verrucomicrobia belonging to the class Spartobacteria, for which the genome was later reconstructed with metagenomics (Herlemann et al., 2013), and the name Candidatus Spartobacteria baltica was proposed. In the current dataset we obtained an OTU identical to S. baltica (OTU_2) that was among the most abundant OTUs in the Baltic Proper, but we also obtained closely related OTUs (>98% identity to OTU_2) with contrasting distribution patterns over the salinity gradient (Supplementary Figure 5).

Eukaryotic Community Composition

We retrieved 18S rRNA gene reads corresponding to nine eukaryotic super-phyla, 24 phyla, 77 classes and 256 genera (see Figure 4 for details). Since no prefiltration was applied before the plankton were captured on filters, the eukaryotic organisms span a large size range, and include unicellular organisms (protists) as well as multicellular organisms (metazoa). The phyla Metazoa, Ciliophora, Dinophyta, Cercozoa, Chlorophyta, Katablepharidophyta, Stramenopiles_X, Ochrophyta, Telonemia, Haptophyta, and Cryptophyta each makes up >1%, and together comprise 94% of the total 18S sequence reads. Stramenopiles_X is a division used in the PR2 database and includes heterotrophic unicellular heterokonts, e.g., the Oomycota, MAST and Labyrinthula. Some samples were highly dominated by metazoan sequences (mainly copepods), likely due to the large number of rRNA genes in their genomes and the fact that these are multicellular organisms. Nevertheless, they never exceeded 70% of the sequences in a sample, leaving at least 17,226 sequences for description of the remaining community of unicellular protists.

FIGURE 4

Figure 4. Summary of detected eukaryotic taxonomic groups. The number of detected OTUs and the mean abundance is given for each taxonomic group at the superphylum, phylum, class and genus level (corresponding to PR2 levels 2, 3, 4, and 7). Only groups displaying >0.1% mean abundance are shown. The dot plot illustrates the distribution of the annotated genera along the salinity gradient, with circle area corresponding to relative abundance.

Of the total 1860 18S OTUs, 283 could be attributed to a known genus or species, 309 could be assigned to PR2 phylogenetic clusters corresponding to the genera level, and 740 could only be classified to higher taxonomic levels. Remaining OTUs could not be classified using the PR2 database, but these corresponded to only 5.63% of the reads in the samples on average. As was observed for the bacteria, identities of closest-related database sequences varied substantially between clades (Figure 2B). Groups with unusually low similarity levels (median < 98%) were Amoebozoa, Choanoflagellida, and Mesomycetozoa.

Eukaryotic Phytoplankton

Recognizing that some protistan groups include phototrophic, mixotrophic and heterotrophic organisms, a division of the organisms in phytoplankton and heterotrophic plankton has been made. Among the dinoflagellates many are mixotrophic, i.e., they carry out photosynthesis but also feed on other organisms. Here mixotrophs are treated as the part of the phototrophic phytoplankton.

Microphytoplankton, made up of diatoms and dinoflagellates, are characteristic of the saline southern and western parts of the Baltic Sea and Kattegat, while smaller eukaryotic phytoplankton groups are more dominant in the northern Baltic Sea (Ojaveer et al., 2010). Pronounced seasonal fluctuations also prevail, with diatoms and dinoflagellates blooming in spring, and smaller eukaryotic phytoplankton and cyanobacteria dominating in summer (Johansson et al., 2004). Diatoms (Bacillariophyceae) and dinoflagellates (Dinophyceae) had higher relative abundances in the more saline regions of our survey (Figure 5A). However, although the samples were taken after the spring bloom, a substantial fraction of the reads from the Baltic Proper also broadly classified as dinoflagellates (Figure 5A). On average, dinoflagellates contributed to 10% of the reads. The relative proportion of dinoflagellates decreased from 50% to nearly 0% of the protist fraction, going from highest to lowest salinity (Figure 5A).

FIGURE 5

Figure 5. Abundance distribution of major eukaryotic taxonomic classes of plankton across the Baltic Sea salinity gradient. Classes displaying >1% mean abundance are shown and colored according to the superphylum they belong to. Green, purple, orange, blue, gray, and red represent Archaeplastida, Hacrobia, Stramenopiles (or Heterokonta), Alveolata, Rhizaria, and Opisthokonta, respectively. (A) Phytoplankton, here broadly defined as groups of organisms that include mainly phototrophic and mixotrophic representatives, although some groups, such as Dinophyta, also include organisms that are heterotrophic. (B) Heterotrophic plankton excluding metazoa. (C) Metazoa. Relative abundances in (A,B) indicate percentages of protist (i.e., non-metazoan) sequence reads; in (C) percentages of all eukaryotic sequence reads.

In the high salinity Kattegat, dominant dinoflagellate genera were Alexandrium and Ceratium of the order Gonyaulacales (Supplementary Figure 6B), while in the Baltic Proper, typical brackish dinoflagellates from the order Suessiales, such as Biecheleria baltica (OTU_25), were most abundant. Major OTUs of diatoms belonged to species often found in Baltic monitoring programs, such as the polar centric diatoms Attheya longicornis (OTU_144) and Chaetoceros neogracile (OTU_136), and the raphid diatom Cylindrotheca closterium (OTU_94). Interestingly, at the station south of Öland in the Baltic Proper (salinity 6.97), the Bacillariophyta and Dinophyceae were both dominated by OTUs (OTU_125 and OTU_32) that displayed highest similarity (99.7 and 98.1% identity) to sequences of uncultured diatoms and dinoflagellates from anoxic layers of the Framvaran Fjord in Norway (Behnke et al., 2010; Supplementary Figures 6A,B). Why these potentially anaerobic diatoms and dinoflagellates flourished in the surface water at this particular station remains to be explained.

Unlike the distribution pattern of dinoflagellates and diatoms, other phytoplankton, such as Chrysophyta, Chlorophyta, and Cryptophyta, were observed at lower salinities. The typical freshwater groups, Chrysophyceae-Synurophyceae (these two classes are combined in the PR2 database), were abundant in the surface water in the Gulf of Bothnia, peaking with 25% of protist reads at salinity 5.4 in the Bothnian Sea (Figure 5A). The dominant OTU (OTU_4) in this area belonged to Uroglena (100% identity), a common colony-forming genus, that made up nearly half of the phototrophic protist community reads. The green algae group, Chlorophyta, was found prevalent in the Baltic Proper and Gulf of Bothnia (Figures 4, 5A). In the Baltic Proper it was dominated by the genus Choricystis (class Trebouxiophyceae) (OTU_9), while in the Gulf of Bothnia two Chlorophyta freshwater OTUs (OTU_103: order Chlamydomonadales, OTU_65: Monoraphidium convolutum), together with one Cryptophyta freshwater species (OTU_56: Rhodomonas minuta) made up one third of the phytoplankton reads (Supplementary Figure 6C).

The ultra-deep sequencing also revealed organisms not easily identifiable with a microscope, such as picoeukaryotes of the class Mamiellophyceae within the Chlorophyta (Figure 4), from which we detected OTUs of the genera Ostreococcus, Micromonas and Bathycoccus (Supplementary Figure 6E). Though they never contributed with more than 2.5% of the protist reads, they displayed distinct distribution patterns along the salinity gradient. For example, Ostreococcus, the smallest known eukaryote, was associated with two OTUs that were separately abundant in Bothnian Sea (OTU_315) and Arkona Basin (OTU_464). Notably, the OTU that peaked in the Bothnian Sea was absent from the high salinity Kattegat region. This OTU belonged to Ostreococcus clade D according to its high similarity (1 deletion) to clade D strains isolated from the Mediterranean (Subirana et al., 2013), while OTU_464 belonged to clade C with its representative sequence being identical to 18S of the widespread Ostreococcus tauri (Keeling, 2007). Our results were consistent with a study on global distribution of Ostreococcus (Demir-Hilton et al., 2011), which reported that Clade D was only found in estuarine or brackish waters, while Clade C was also found in marine waters.

The most abundant Haptophyta OTUs were most closely related to Coccolithophorids (class Prymnesiophyceae), which are globally important phytoplankton with characteristic exoskeletons of calcium carbonate plates (coccoliths). They are found in large numbers throughout the ocean's euphotic zone but only one species, Balaniger balticus, has been described in the Baltic Sea (Thomsen and Oates, 1978). The sparsity of coccolithophorids in the Baltic has been attributed to undersaturation of calcium carbonate in winter, potentially resulting in dissolution of the coccoliths (Tyrrell et al., 2008). We found two coccolithophore-related OTUs (OTU_127 and OTU_2007) that reached a maximum of 4.2% of the protist reads in the high salinity samples in Kattegat. The representative sequences of these two OTUs showed 100% identity to both coccolithophore species Emiliania huxleyi and Gephyrocapsa muellerae from the family Noelaerhabdaceae, which, incidentally, have identical sequences (Bendif et al., 2015). Intriguingly, OTU_127 was also detected, though at much lower abundances, in the Gulf of Bothnia (Supplementary Figure 7). Here we also detected the non-calcifier Isochrysis galbana (OTU_585) that belong to the same order (Isochrysidales) as the coccolithophores. As far as we know, neither Emiliania/Gephyrocapsa or Isochrysis have been observed in the Baltic Sea before. Although OTU_127 dominated the prymnesiophytes in Kattegat, coccolithophores were not the predominant organisms from this group in the Baltic regions. Instead, order Prymnesiales was dominant with OTUs widely spread according to their salinity preferences (Figure 4). Prymnesium, a genus that include the potentially toxic spring-bloom forming species Prymnesium polylepis (Gorokhova et al., 2014), was present in all samples except those collected in the Bothnian Bay (Figure 4).

While the above data on eukaryotic phytoplankton are based on 18S rRNA gene amplicons, phytoplankton can potentially also be monitored using 16S sequencing of chloroplast sequences (Eiler et al., 2013; Lemieux et al., 2014). In line with the fact that our 16S primer pair matches a fair amount of chloroplast sequences (65% of chloroplast sequences in RDP v.11.4, Cole et al., 2014), 120 of the 16S OTUs, corresponding to on average 0.3% of the reads, were classified as chloroplasts by SILVA, and a subset of these (98 OTUs) could be further classified by PhytoREF (Decelle et al., 2015). Ideally, the 16 and 18S data should reveal the same trends across samples for the same clade or species. We tested this by first correlating relative counts at the class level. For the seven classes that had at least 1000 reads in total in each dataset (16 and 18S), six displayed significant correlations (P < 0.05; average Spearman correlation coefficient = 0.75). Comparing at a more detailed taxonomic level was complicated by the fact that the 16S OTUs lacked a detailed taxonomic annotation. Instead we correlated all 16S OTUs to all 18S OTUs and extracted the pairs displaying a Spearman rank order correlation coefficient of at least 0.8. Of the 304 pairs found, six included a 16S OTU classified as chloroplast. For five of these six pairs, the 16 and 18S OTU were classified to the same taxonomic class (four different classes; Supplementary Figure 8). The probability that this level of matching would occur by chance is extremely low (P < 10⁻⁶; see Methods), and hence these five pairs of OTUs likely represent the same species or genus. This demonstrates the robustness of the metabarcoding approach and indicates that this type of correlation analysis can be fruitful for linking chloroplast 16S sequences to their nuclear 18S counterparts.

Eukaryotic Heterotrophs

The detected heterotrophic protists mainly belonged to the phyla Ciliophora, Cercozoa, Dinophyta, Katablepharidophyta, and Telonemia (Figure 5B). Ciliates (Ciliophora) that carry out photosynthesis by using klepto-chloroplasts belong to the functional group phytoplankton but are treated among the heterotrophs is in this article. Heterotroph dominance varied between class Spirotrichea from the phylum Ciliophora, in the Gulf of Bothnia and the Southern Baltic, and class Filosa-Thecofilosea from the phylum Cercozoa, in the Northern Baltic Proper (Figure 5B). The abundance of both these classes dropped in the Kattegat, as well as in the “intermediate” Arkona sample (salinity 9.7).

Of the 173 OTUs classified as Ciliophora, 32 had a mean relative abundance >10⁻³. The most prominent OTUs belonged to the genus Vorticella (OTU_2) from class Oligohymenophorea and to Strombidium (OTU_5) from class Spirotrichea (Figure 5, Supplementary Figure 6F). Only two OTUs (OTU_425, OTU_666), that displayed rather low abundance, belonged to the genus Mesodinium (Mironova et al., 2009), normally abundant in the Baltic Sea. The low yields could be due to diurnal vertical migration, or that the Mesodinium population was present in thin layers that were not sampled (Sjöqvist and Lindholm, 2011). Cercozoa is a highly cryptic phylum that contains widespread protozoan omnivores and parasites with vague morphological characters (Cavalier-Smith, 1998; Bass et al., 2009; Weber et al., 2012). Of the 170 Cercozoa OTUs, 25 had a mean relative abundance >10⁻³ (Figure 5, Supplementary Figure 6G).

Other heterotrophic strategies, such as parasitism, also play essential roles in the microbial food web, but monitoring parasites or symbionts by microscopy can be difficult. Sequencing of marker genes can greatly facilitate depiction of the distribution pattern of these organisms. For instance, Syndiniales, a class of Dinophyta that are parasites on other dinoflagellates, were found where other dinoflagellates were abundant (Figures 5A,B). In total, 93 OTUs were found from three Syndiniales groups (group I, II, III) (Guillou et al., 2008). The dominant OTUs were mainly from group I (Figure 4), while group II was the most diverse with 56 OTUs.

In the Hacrobia division, three major heterotrophic classes, Picobiliphyta (Picozoa), Katablepharidaceae, and Telonemia, showed contrasting distribution patterns in the Baltic. The Picobiliphyta, initially described from cold polar waters (Lovejoy et al., 2006), and recently renamed as Picozoa (Seenivasan et al., 2013), is a highly diverse deep-branching clade, with representatives found in a variety of marine ecosystems (Not et al., 2007, 2009). In the present study, the dominant Picozoa OTU (OTU_44), possessing almost 80% of the Picozoa reads, was identical to the only cultured member of the Picozoa, Picomonas judraskeda, that feed on small organic particles (Seenivasan et al., 2013). These picoeukaryotes (2–5 μm length) have previously been enriched from brackish Baltic Sea water under dark incubations (Weber et al., 2012). Abundance of Picozoa peaked in the Bothnian Sea, while the larger Katablepharydophyta and Telonemia accumulated mainly in the Baltic Proper and the Arkona basin.

The superphylum Stramenopiles were represented by the Ochrophyta and Stramenopiles_X. The 70 Stramenopiles_X OTUs were mainly annotated to the classes Labyrinthulea, Oomycetes, and MAST. One of the samples from Bothnian Bay (salinity 2.4) displayed markedly higher abundance of Labyrinthulea than the others (Figure 5B). Class Labyrinthulea was highly dominated by an OTU (OTU_54) that was found in all Baltic regions with no apparent salinity preference. The representative sequence of OTU_54 was highly similar (99–100% identity) to endophytic Labyrinthula strains isolated from either mangrove leaves or eelgrass, which indicates that OTU_54 is an endophyte associated with plants. Unlike other sampling sites, the depth at this station is rather shallow (17-m), and directly to the east a region of shallow waters (10-m depth on average) extends to the Finnish coast. Hence, the high abundance of Labyrinthula here could be due to debris from benthic macro vegetation captured on the filter.

Of the 70 OTUs from Stramenopiles_X, 30 were annotated to the MAST lineage, a diverse group of pico- to micro-sized heterotrophic marine stramenopiles that was discovered by environmental sequencing only 10 years ago (Massana et al., 2004) and has not yet been cultured (Massana et al., 2006). Many studies have since then investigated the links between environmental factors and MAST lineages (Massana et al., 2004, 2014; Piwosz and Pernthaler, 2010; Lin et al., 2012; Thaler and Lovejoy, 2013). Here, seven MAST clades (“ribogroups”; Massana et al., 2014) were detected. The most prevalent were MAST-1 and -2 (Figure 4, Supplementary Figure 6H), which are commonly found in marine and coastal regions (Massana et al., 2004, 2006). To date, four subclades of MAST-1 are known (Massana et al., 2014), however, only MAST-1C was observed in our dataset, with five associated OTUs (Supplementary Figure 6H). Besides the most prevalent MAST clades (MAST-1C, -2D), MAST-4D and -6 were abundant in the brackish region (salinity 4–10). For MAST-6, all OTUs (OTU_145, 245, and 644) were missing in Bothnian Bay and Kattegat, which indicates these are brackish specialists. In the Gulf of Bothnia, MAST-2A, -3J, -12C, and -12_X displayed pronounced preference of low salinity (<4) (Supplementary Figure 6H). So far, only sequences from MAST-2 and -12 have been found in freshwater (Massana et al., 2014). The OTU (OTU_337) from MAST-3J that only appeared in the Bothnian Bay indicates that also this subclade contains freshwater members. In contrast, other MAST-3 subclades only appeared in the high salinity Kattegat samples.

The detected unikonts (superphylum Opisthokonta) in this study included Choanoflagellates, Mesomycetozoa, Fungi, and Metazoa, with 195 OTUs annotated in total. The most divergent and abundant group was Metazoa (Figure 5C) with 105 OTUs of which 15 had a mean relative abundance >10⁻³. This can be compared with 217 species of metazooplankton in the Baltic Sea based on zooplankton inventories (Telesh et al., 2009). A recent review of plankton diversity in the Baltic ranked 23 abundant mesozooplankton species in different subbasins based on monitoring data (Ojaveer et al., 2010). OTUs corresponding to these Metazoa species displayed remarkably similar distribution patterns (Figure 6). Unicellular unikonts (Choanoflagellates, Mesomycetozoa, Fungi) made up a limited fraction of protist reads (Figure 5B). The low abundance of choanoflagellates may be attributed to seasonal dynamics, since unlike other heterotrophic nanoflagellates, they have been reported as being sparse during summer in the Baltic Sea (Samuelsson et al., 2006).

FIGURE 6

Figure 6. Abundance distribution of dominant mesozooplankton (Metazoa) across the salinity gradient in comparison to historical monitoring data. Data from Figure 1 and Table 3 of Ojaveer et al. (2010) is indicated with triangles connected by gray lines (right y-axis). Relative abundances of the OTUs for the corresponding species (when a species corresponds to multiple OTUs these were summed) are indicated with circles connected by blue lines (left y-axis).

Trends in Alpha-Diversity along the Salinity Gradient

The deep sequencing allowed us to investigate trends in alpha-diversity over the salinity gradient without the limitations associated with microscopic detections. Since alpha-diversity estimates are biased by sequencing depth (Lundin et al., 2012), the same number of reads (n) were subsampled from each sample, where n equals the lowest number of reads available among samples for the community or phylum analyzed. No clear trends in species richness (number of observed OTUs) or combined richness-evenness (Shannon-Wiener index) could be observed for either bacteria or protists at the total community level (Figure 7). When zooming-in on specific eukaryotic phyla (here n reads were sampled from within each phylum) no consistent trends were observed, but one of the phyla, Ochrophyta, appeared to have lower diversity in the intermediate salinity range. Comparing the species richness and Shannon index between the samples in the horohaliculum (salinity 5–8) with the other samples gave significant results only for Ochrophyta, that had significantly lower richness in the horohalinicum (Wilcoxon rank sum test, False Discovery Rate-adjusted p = 0.03; Figure 7).

FIGURE 7

Figure 7. Trends in planktonic alpha-diversity along the salinity gradient. Circles connected by colored lines indicate observed OTU richness (left y-axes) and triangles connected by gray lines represent Shannon-Wiener diversity index (right y-axes). In order to remove biases due to variation in sequence depth, the same number of reads (indicated with n in each panel) were subsampled from each sample. The subsampling procedure was repeated 1000 times and the averages and standard deviations (error bars) were calculated. For Shannon-Wiener standard deviations were negligible and are not shown. (A) Total bacterial community. (B) Total protist community (non-metazoan eukaryotes). (C–L) Eukaryotic phyla with mean relative abundance >1%. For each phylum (n) reads were subsampled from reads belonging to the respective phylum, to avoid that diversity indices were biased by the relative abundance of the phylum within the sample. A minimum of 1000 reads were subsampled per phylum/sample; samples containing less than 1000 reads from the phylum were excluded. Phyla were colored according to superphylum in accordance with Figures 4, 5. Q-values indicate False Discovery Rate-adjusted p-values of Wilcoxon rank sum tests comparing Shannon-Wiener index and OTU richness in the horohalinicum samples (salinity 5–8; n ≤ 11) and the other samples (n ≤ 10).

Discussion

The Baltic Sea has an exceptional long record of plankton monitoring, based primarily on microscopy, from the HELCOM-COMBINE monitoring program (HELCOM, 2015), but to our knowledge this is the first investigation of eukaryotic plankton diversity along the full salinity gradient using high-throughput sequencing. Compared to traditional microscopy based methods, high-throughput sequencing of taxonomic marker genes (metabarcoding) potentially offers a more comprehensive view of microbial diversity. While not yielding absolute cell counts, sequencing data presents integral information of microbial communities, and reveals trends in community changes. While morphological features in some cases can be used to classify plankton with high taxonomic resolution, this is not the case for the majority of species that are too small to allow classification by light microscopy. This is probably the reason why, with only 1 week's sampling, the present metabarcoding-based study revealed a microbial diversity on par with what has been recorded in the Baltic Sea with microscopy over a century (Hällfors, 2004; Telesh et al., 2009). In total we observed 1860 OTUs (548 broadly classified as phytoplankton and 774 as heterotrophs) which can be compared with 2075 known species of phytoplankton (Hällfors, 2004) and 1030 known species of micro-, meso-, and macrozooplankton (Telesh et al., 2009).

When performing metabarcoding, the taxonomic accuracy is limited primarily by the resolution of the marker gene, the coverage of the primers, and the quality and completeness of the reference database. Here we used the SSU rRNA gene, the most widely used taxonomic marker for bacterioplankton and protists (Amann et al., 1995; Pace, 1997). While it provides a lower taxonomic resolution than e.g., the internal transcribed spacer (ITS) region, its highly conserved regions allow the use of primers with domain-wide taxonomic range. Databases are also more complete than for other markers. We amplified and sequenced 450–600 bp of the hypervariable regions V3–V4 in bacteria and V4–V5 in eukaryotes, using in-house designed primer pairs that cover a large proportion of taxa in each domain (Hugerth et al., 2014a,b). Of the generated 18S OTUs, 31.8% could be classified to the genus level using our automatic classification scheme (see Methods), while 58.3% did not have sufficiently close matches (>98% identity) in the PR2 database (i.e., insufficient database completeness) and 9.8% matched multiple genera (i.e., insufficient taxonomic resolution of the marker). An example of the former condition is present for sequences classified to Amoebozoa, which displayed unusually low identity levels to reference sequences (Figure 2). Interestingly, the bacterial phylum Chlamydiae, that are symbionts on e.g., fish and amoeba in aquatic environments (Collingro et al., 2005; Stride et al., 2014) also displayed unusually low identity to reference sequences. Thus, the Baltic Sea appears to be inhabited by a set of poorly characterized Amoebozoa with poorly characterized bacterial symbionts. Examples of the latter were E. Huxleyi and G. muellerae (Bendif et al., 2015) that have 100% identical sequences in the V4–V5 region, as was also the case for species of Labyrinthula isolated from different hosts. Ideally, one would like to sequence both a part of the SSU rRNA gene and a more rapidly evolving marker, such as the ITS. With recently developed emulsion-based protocols, it will be possible to connect different markers from the same cell (Borgström et al., 2015; Spencer et al., 2016).

OTU relative abundances need to be carefully interpreted due to variations in the respective copy numbers of rRNA genes among phylogenetic groups, especially for eukaryotes (Prokopowich et al., 2003; Větrovský and Baldrian, 2013). Due to large cell volumes and vast genomes, dinoflagellates, diatoms, and ciliates in general possess more copies of the 18S rRNA gene than smaller flagellates, such as prasinophytes (Zhu et al., 2005; Godhe et al., 2008; Gong et al., 2013). Although, previous studies have shown a correlation between cell size (volume or length) and rRNA gene copy number (Godhe et al., 2008), copy numbers can also vary among species of the same genus. Thus, comparing the relative abundance among different taxonomic groups may be misleading. However, the variation in relative abundance of the same group across samples is still valid, allowing to test for correlations with variations in e.g., environmental parameters. Most taxonomic groups, at all taxonomic levels, displayed strikingly smooth abundance changes along the salinity gradient sampled (Figures 3–5, Supplementary Figure 6) in contrast to the spiky abundance distributions to be expected if the methodology had introduced substantial random noise. This was also true for metazoa (Figure 6), which is perhaps unexpected given that the low volumes of water sampled here (<1 l) are suboptimal for quantifying metazoa with their relatively low abundances in the water column. Systematic biases still likely exist, caused by e.g., primer mismatches and differences in rRNA gene copy numbers, but these should affect all samples equally.

Extracellular DNA released from dead (Kloos et al., 1994; Choi et al., 2004; Rice et al., 2007) or live (Böckelmann et al., 2006) cells may also lead to skewed estimates of abundance and diversity. Furthermore, metabarcoding detects all life forms that carry the marker gene, muddling the distinction between vegetative cells and other life stages. Some species of dinoflagellates, chrysophytes, cryptophytes, and haptophytes are known to form resting cysts after blooms (Lichtlé, 1980; Adam and Mahood, 1981; Jordan and Chamberlain, 1997; Lotter et al., 1998; Warns et al., 2012). One example is Alexandrium, a toxic genus of dinoflagellates (Teegarden and Cembella, 1996) that sporadically blooms in the Baltic (Witek, 2004; Kremp et al., 2009), for which we observed high abundance (9.1% of protist reads) in one of the Kattegat samples (Supplementary Figure 6B). Whether this reflects a bloom is hard to tell, since the sequences may be derived from cysts (Oshima et al., 1992). However, cysts are also toxic (Oshima et al., 1992), thus in some cases detecting both cysts and vegetative cells can be useful to track toxic outbreaks.

Similar to what was observed in previous studies (Herlemann et al., 2011; Dupont et al., 2014), bacterial community composition changed gradually along the salinity gradient, and difference in community composition (beta-diversity) was significantly correlated with difference in salinity. Here we showed for the first time a similar pattern for eukaryotic plankton composition using molecular data. The observed pattern can probably also be driven by a combination of other factors such as availability of micronutrients and composition of dissolved organic matter that covary with salinity across the gradient, and may be different or weaker during other seasons than summer, but it is well established that salinity difference per se is a major barrier for species to cross (Logares et al., 2009). A recent study shows that the bacterioplankton inhabiting the Baltic Proper are not locally adapted freshwater or marine populations, but are rather members of a global brackish metacommunity that most likely adapted to brackish conditions before the Baltic Sea was formed (Hugerth et al., 2015). Whether this also holds true for protists inhabiting this sea remains to be investigated. For endosymbiotic protists, isolated from the surrounding water and protected from osmotic stress, the distribution patterns are likely mainly determined by the salinity preference of the host.

The “Remane curve” (Remane, 1934) illustrates how diversity of freshwater, marine, and brackish specialist species vary along a salinity gradient. Total diversity is lowest in the horohaliculum (salinity 5–8) because few freshwater and marine species can tolerate this salinity level and relatively few true brackish specialists exist that occupy this niche. Remane based his “Artenminimum” (species-minimum) concept (Remane and Schlieper, 1971) on macrozoobenthos, but since then the generality of the pattern has been debated (Reid, 1961; McLusky, 1989; Attrill, 2002; Portillo et al., 2012). Recently, Telesh et al. suggested the opposite pattern for protists, with a species-maximum in the horohaliculum (Telesh et al., 2011b). The analysis was based on pooling species observations within salinity ranges from multiple studies, but has been questioned, since the intermediate salinity where the highest diversity was observed was also the most frequently sampled (Ptacnik et al., 2011; Telesh et al., 2011a), hence the pattern could stem from sampling biases. Here we took the opportunity to measure alpha-diversity with our metabarcoding dataset. This type of data has the potential to give an estimate of diversity that is not biased by morphological characters of the specimens, and the deep sequencing allows subsampling of each sample to the same depth, improving comparisons between samples (Lundin et al., 2012). Similar to what was previously observed (Herlemann et al., 2011), total bacterial alpha-diversity did not display any clear trends over the freshwater-marine continuum. Likewise, the total protist community (eukaryotic data excluding metazoa) displayed neither a species-minimum nor a species-maximum in the horohaliculum. When zooming in on specific eukaryotic groups, only one phyla, Ochrophyta, displayed significantly different diversity in the horohaliculum. Ochrophyta includes among other groups the diatoms, and since the sampling was conducted well after the spring bloom when diatoms are abundant in the Baltic Proper, the pattern may have been different if seasonal sampling had been conducted. Overall, there is considerable noise in the data, and the limited number of samples were collected during a single cruise in the summer, thus more data is obviously needed to derive stronger conclusions. However, if pronounced differences in alpha-diversity did exist between salinity regimes, this should have been visible. Telesh attributed the lack of an “Artenminimum” for protists to their fast evolutionary rate and hence ability to adapt to brackish conditions (as opposed to macroorganisms; Telesh et al., 2011b), which is also what was suggested for bacteria (Herlemann et al., 2011). However, our recent study showing that many of the bacterioplankton inhabiting the Baltic Sea are brackish specialists adapted to these conditions before the Baltic Sea was formed (Hugerth et al., 2015), and likely migrated from other brackish waters, opens the possibility that the same applies to some of the protist species. Future population genomic analysis may be able to reveal such patterns.

Author Contributions

YH, BK, and AA conceived and designed the study. BK performed sampling. YH performed molecular work. YH, AA analyzed the data. All authors wrote the paper, and read and approved the final version of the manuscript.

Funding

This work was supported by BONUS BLUEPRINT project, supported by BONUS (Art 185), funded jointly by the EU and the Swedish Research Council FORMAS. It is also funded by the Swedish Research Council VR (grant 2011-5689) through a grant to AA. The work was in part supported by the JERICO project from the European Union Seventh Framework Programme (FP7/2007–2013) under grant agreement no. 262584, and through the Swedish Research Council VR through the Swedish Lifewatch project. YH was supported by a scholarship from the China Scholarship Council (CSC#201206950024).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to TransAtlantic AB and the captains and crew of the cargo vessel TransPaper for providing space for the Ferrybox-system and technical assistance during the installation of the system and valuable help during sampling. We are also grateful to Jürg Brendan Logue and Conny Sjöqkvist for insightful comments and suggestions on the manuscript. Sequencing was conducted at the Swedish National Genomics Infrastructure (NGI) at SciLifeLab in Stockholm. Computations were performed on resources provided by the Swedish National Infrastructure (SNIC) through the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX).

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb.2016.00679

References

Adam, D. P., and Mahood, A. D. (1981). Chrysophyte cysts as potential environmental indicators. Geol. Soc. Am. Bull. 92, 839–844. doi: 10.1130/0016-7606(1981)92<839:CCAPEI>2.0.CO;2