Variability and Community Composition of Marine Unicellular Eukaryote Assemblages in a Eutrophic Mediterranean Urban Coastal Area with Marked Plankton Blooms and Red Tides

: The Thessaloniki Bay is a eutrophic coastal area which has been characterized in recent years by frequent and intense phytoplankton blooms and red tides. The aim of the study was to investigate the underexplored diversity of marine unicellular eukaryotes in four different sampling sites in Thessaloniki Bay during a year of plankton blooms, red tides, and mucilage aggregates. High-Throughput Sequencing (HTS) was applied in extracted DNA from weekly water samples targeting the 18S rRNA gene. In almost all samples, phytoplankton blooms and/or red tides and mucilage aggregates were observed. The metabarcoding analysis has detected the known unicellular eukaryotic groups frequently observed in the Bay, dominated by Bacillariophyta and Dinoflagellata, and revealed taxonomic groups previously undetected in the study area (MALVs, MAST, and Cercozoa). The dominant OTUs were closely related to species known to participate in red tides, harmful blooms, and mucilage aggregates. Other OTUs, present also during the blooms in low abundance (number of reads), were closely related to known harmful species, suggesting the occurrence of rare taxa with potential negative impacts on human health not detectable with classical microscopy. Overall, the unicellular eukaryote assemblages showed temporal patterns rather than small-scale spatial separation responding to the variability of physical and chemical factors.


Introduction
During the past decade, the advancement of High-Throughput Sequencing (HTS) technologies have provided unprecedented insights into the diversity [1][2][3][4][5], spatial community composition [6][7][8], temporal patterns [9][10][11][12], and functional and ecological processes affected by marine unicellular eukaryotes [9,13,14]. In a recent study, Ramond et al. [14] have attempted to link the taxonomic diversity of protistan taxa detected with HTS with 30 different traits in samples from ecosystems across the coasts of North Atlantic Ocean, suggesting that small taxa are characterized by broader taxonomic and functional diversity. Furthermore, previously underestimated trophic strategies, such as parasitism and decomposition, have been recognized as central to marine planktonic assemblage structure and ecosystem function [9,15]. In addition, the environmental effects [9,16] and biotic relationships in unicellular eukaryotic communities have been explored, using classical and modern ecological tools [17,18]. It is evident that the plethora of the genetic and ecological information produced provide a better understanding of marine global ecosystem processes [19,20].
In particular, coastal areas are systems with high risk of pronounced modification in view of climate change and coastal urbanization leading to increased eutrophication [21,22]. Thus, there is an increase in the attention of marine microbiologists thanks to the development of HTS tools. Either by focusing on specific taxonomic groups, such as the parasitoids Perkinsozoa [23] and Labyrinthulomycetes [24,25], the cosmopolitan and conspicuous tintinnid ciliates [26], the newly described Vampyrellida (Rhizaria) [27], or by investigating the overall marine protistan diversity [28,29], temporal dynamics [11], function [14], biotic interactions [7], relationships with abiotic variables [9], effects of climate change [13] and potential activity [5], marine microbiologists have attempted to answer the classical questions of who, when, why, and where in coastal systems. However, few HTS-based studies have investigated marine unicellular eukaryotic assemblages in degraded eutrophic coastal areas, with high anthropogenic impacts. These studies, albeit scarce, have revealed high and novel diversity [3,9,30], and have attempted to provide insights of the trophic relationships between known groups and previously overlooked links [7,31].
Here, we use HTS to investigate the marine unicellular eukaryotic diversity and spatial distribution throughout a year of frequent samplings in a coastal eutrophic Mediterranean urban area (Thessaloniki Bay). [32,33]. This approach was used as a follow-up to a previous publication examining via classical microscopy tools the microbial plankton communities in the same area and the same sampling period, which was characterized by marked harmful algal blooms, red tides, and mucilage aggregate incidents [33]. The Thessaloniki Bay is the inner part of Thermaikos Gulflocated in northern Greece-and has been accepting for decades a large volume of domestic and industrial wastes from the one million residents of the city of Thessaloniki. Although wastewater treatment has been implemented, decreasing the effects of anthropogenic eutrophication [34], the restricted water circulation and shallowness of the Bay in combination with high nutrient inputs, lead to harmful algal blooms, apparent red tides, and mucilage aggregates [33], leading to growing concerns of the citizens and authorities on the water quality of the Bay and particularly of the urban front. Data on the abundance and dynamics of plankton community (both phyto-and protozooplankton) in the urban part of the Gulf obtained by classical microscopy [32,33,35], have attempted to identify and characterize the marine protist abettors and perpetrators responsible for these phenomena. These studies have shown the contribution of taxa known to thrive due to eutrophication and to cause harmful events, such as the diatoms Cylindrotheca closterium, Skeletonema costatum, Chaetoceros spp., and the dinoflagellates Noctiluca scintillans and Dinophysis acuminata. Taking into consideration the continuous increase of eutrophication in coastal areas worldwide, there is a need to identify the causative organisms of these events and implement strategies to mitigate their severe impacts to ecosystem and public health. HTS-based studies investigating the unicellular eukaryotic diversity, dynamics, spatial distributions, and environmental effects in coastal eutrophic areas could contribute towards the early and accurate characterization of the microbes responsible for harmful bloom events and red tides.
The aim of this study was to investigate and describe the overall marine unicellular eukaryotic assemblages of different sites of the urban front of Thessaloniki throughout a year in short sampling intervals via HTS tools, following a previous publication concerning the same area and sampling period, where different periods of plankton blooms were identified and analysis of the environmental variables was presented. The sites are located in close proximity, thus allowing us to examine the variability of unicellular eukaryotic diversity and dynamics in small temporal and spatial scales. The main questions were: (i) Do the assemblages include previously undetected taxa in the area? (ii) How do these taxa contribute to the phenomena of blooms, red tides, and mucilage aggregates? (iii) What are the temporal and spatial distributions of the unicellular eukaryotic assemblages? (iv) Are they affected by nutrients and how?

Sampling Sites and Sample Collection
The sampling procedure was described analytically elsewhere [33]. Briefly, samples were collected weekly from March 2017 to February 2018 at an urban site in Thessaloniki Bay (Thermaikos Gulf), namely the White Tower (WT) station. Simultaneously, samples were taken from WT and three other adjacent sites (Aretsou Beach, AR; Music Hall coast, MH; and Harbour, HB) every four weeks. The stations were located along the coastline with a distance of ~ 3 Km between the closest consecutive ones and 10 Km between the most distant ones, which were located at the edges of the sampling area ( Figure 1, Table 1). All sampling sites had a maximal depth of 4 m. In total, 83 water samples were collected. In order to better illustrate temporal vs. spatial unicellular eukaryotic diversity and community structure, the results are presented separately for the 47 samples of the WT site (representing the temporal aspect), and for the 48 common sampling dates of all four sites (representing the temporal vs. spatial aspect). For all samples, in situ measurements of water temperature and salinity were made with the use of the YSI Pro 1030 instrument (YSI Inc., Ohio, USA). Water samples of 2 L were collected from the surface layer of 1 m, and subsamples of 100-250 mL (depending on plankton and particulate matter density) were immediately filtered onto 0.7 μm prewashed (in 5% to 10% HCl) and precombusted (6 h, 550 °C) Whatman GF/F filters, and the filters were stored in −20 °C for particulate organic nutrient and Chlorophyll a (Chl a) measurements. Moreover, subsamples of 50 mL were filtered through 0.2 μm cellulose acetate filters (Sartorius) and were kept in −20 °C for dissolved inorganic nutrient measurements. Finally, the rest of the water volume was prefiltered through a 200 μm sterile mesh, and subsamples of 100-200 mL from each water sample were subsequently filtered through 0.2 μm nucleopore filters and similarly were kept in −20 °C for molecular analysis.

DNA Extraction and Sequencing
DNA was extracted from the 83 samples using a Macherey-Nagel NucleoSpin ® Soil, Genomic DNA Isolation Kit, according to the manufacturer's instructions. The concentration and quality of recovered DNA was confirmed using the Thermo Scientific™ NanoDrop™ spectrophotometer.
The extracted DNA was subjected to PCR using specific primers targeting the V4 hyper variable region of the 18S rRNA gene (E572F = CYGCGGTAATTCCAGCTC; E1009R = AYGGTATCTRATCRTCTTYG) [38]. These primers have been found to successfully amplify approximately 440 bp in the V4 hypervariable region of all the major high-level eukaryotic taxonomic groups up to 92% coverage, after in silico analysis via the SILVA TestPrime 1.0, done at the Integrated Microbiome Resource (IMR) at Dalhousie University (Halifax, Canada). The PCR products were verified visually by running on a high-throughput Hamilton Nimbus Select robot (Hamilton Company, Reno, USA) using Coastal Genomics Analytical Gels (Hamilton Company, Reno, USA) and were cleaned-up and normalized using the high-throughput Charm Biotech Just-a-Plate 96-well Normalization Kit (Charm Biotech, USA). The amplicon samples were sequenced on Illumina MiSeq using 300 + 300 bp paired-end chemistry which allows for overlap and stitching together of paired amplicon reads into one full-length read (http://cgeb-imr.ca/protocols.html). The PCR amplification step and sequencing steps were performed at the Integrated Microbiome Resource (IMR) at Dalhousie University, Halifax, Canada.

Read Processing
The produced reads were subjected to downstream processing using the mothur v.1.34.0 software [39], following the proposed standard operating procedure [40]. Briefly, forward and reverse reads were joined, and contigs below 200 bp, with >8 bp homopolymers and ambiguous base calls were removed from further analysis. The remaining reads were dereplicated to the unique sequences and aligned independently against the SILVA 132 database, containing 163,916 eukaryotic SSU rRNA gene sequences [41]. Then, the reads suspected of being chimeras were removed using the UCHIME software [42]. The remaining reads were clustered into Operational Taxonomic Units (OTUs) at 97% similarity level [43], using the average neighbor method in mothur. In order to obtain a rigorous dataset, OTUs with a single read in the entire dataset were removed from the analysis, as they were suspected of being erroneous sequences [7,44,45]. The resulting dataset was rarefied with the subsample command in mothur v.1.34.0 to 14,258 reads, while seven samples with lower total number of reads (between 9186 and 11,851 total reads per sample) were not removed from the dataset. We chose this course of action as the best compromise in order to retrieve the majority of the biodiversity detected by sequencing, as rigorous subsampling may result in some OTUs loss, but also to include all samples in the analysis, even those with lower number of reads. Taxonomic classification was assigned using BLASTN searches against the "Protist Ribosomal Reference database" (PR2 database), with curated protistan taxonomy [46], by applying a lowest accepted level of > 80% similarity with a closest relative. The reads belonging to OTUs related to metazoa and Embryophyceae (land plants) were removed from the dataset, while reads affiliated to unicellular eukaryotes were included in the final data matrix. The raw reads were submitted to GenBank-SRA under the accession number PRJNA552665.

Data Analysis
All data analyses were performed on the rarefied dataset. Α-diversity estimators (the richness estimator SChao1; the Simpson, Shannon, and Evenness indexes) were calculated for all samples with the PAST 3.16 software [47]. Paired t-tests were applied in PAST 3.16 software to compare the richness, and the α-diversity estimators between the four sampling sites. The eukaryotic assemblages of the different samples at the four sites during the common sampling dates were compared using the Plymouth routines in the multivariate ecological research software package PRIMER v.6 [48]. The Bray-Curtis similarity coefficients were calculated based on OTUs abundances with log(x + 1) transformation, in order to identify interrelationships between samples and construct a cluster dendrogram. Finally, the relationship among all samples (by using the OTUs abundance matrix) and environmental variables in each of the four sites was explored with the Canonical Correspondence Analysis (CCA; CANOCO) [49]. Briefly, the Monte Carlo permutation test (n = 999) was implemented in order to identify the set of parameters that best signify/explain the observed clustering of samples. The CCA was run using only the environmental parameters with P-value < 0.05, as indicated by the Monte Carlo test. Furthermore, the relationship between the dominant OTUs individually, with relative number of reads > 1% of the total number of reads in each site (indicated as instantly abundant by Logares et al. [6]) and environmental variables was also investigated via CCA. Firstly, the Monte Carlo test was run in order to select the abiotic variables that significantly affected the samples or OTUs, respectively (P-value < 0.05, Inflation factor < 20), and then CCA was performed by CANOCO [49]. The ordination plots, including samples / OTUs and environmental variables, were created in CanoDraw. The significance of the axes obtained by the CCA was determined based on the Monte Carlo permutation test [49].

Environmental Variables and Bloom Periods
Environmental data are presented elsewhere [33]. Briefly, the seawater temperature during the year of the study ranged from 9.6 °C on 24 January 2018 at the WT site, to 29.7 °C on 09 August 2017 at the same site (Supplementary Figure Figure S1).
Based on microscopy data from Genitsaris et al. [33], conspicuous plankton blooms were detected and analyzed. The periods of the blooms are presented briefly in Table 2.

Temporal Dynamics of the Unicellular Eukaryotic Community
Rarefaction curves of all samples reached a plateau in the majority of cases (Supplementary Figure S2). The ratio of observed to expected (SChao1) number of OTUs in all samples was > 0.73, and in > 90% of the samples, the ratio was > 0.8 (Supplementary Tables S1 and S2). A total of 2821 OTUs, affiliated to 19 high-level taxonomic groups, were detected in all sites (Supplementary Figure S3). Dinoflagellata was the most diverse high-level group comprising of 35% of the total number of OTUs, followed by Cercozoa (13%), Ciliophora (9%), and Bacillariophyta (5%). All other high-level taxonomic groups consisted of < 3% of the total number of OTUs each (Supplementary Figure S3).
In the 47 samples from the WT site, a total of 2256 OTUs were recovered. The highest number of OTUs (378 OTUs) was detected at 23 August and the lowest (44 OTUs) at 22 March (Supplementary  Table S1). The α-diversity estimators fluctuated during the study, with extremely low values on 22 March 2017 at the WT site (Simpson: 0.03, Shannon: 0.12, and Evenness: 0.03), indicating high dominance by one (or few) taxa, and high variation between taxa abundances within the community. Indeed, OTU_2 closely related to the dinoflagellate Noctiluca scintillans (Table 3) Table S2). Overall, 560 OTUs were found to be common in all four sampling sites (Supplementary Figure S4). Concerning each site individually, 172, 231, 351, and 178 unique OTUs were detected at WT, AR, MH, and HB, respectively (Supplementary Figure S4). The average values of the a-diversity estimators were generally similar between all four sites, i.e., the Simpson estimator between 0.73 at WT and 0.86 at MH, the Shannon estimator between 2.55 at WT and 2.99 at AR, and the Evenness estimator between 0.07 at HB and WT, and 0.1 at the other two sites (Supplementary Table S2). The plot of the cumulative number of the OTUs (Figure 2b) revealed a continuous influx of newly introduced OTUs per new sample in all sampling sites, even higher than in WT during the weekly samples, due to the largest interval between samplings (four weeks vs. one week). The average number of newly introduced OTUs ± the Standard Deviation were 90 ± 59.5, 96.2 ± 49.2, 105.3 ± 44.9, and 80.7 ± 49.4 at WT, AR, MH, and HB, respectively. Dinoflagellata in all cases comprised > 30% of the number of reads per sample, while during the warm months of summer, Bacillariophyta showed higher relative abundance reaching > 50% of the total number of reads at all sites except HB (Figure 3). Furthermore, > 40% of the total number of reads were affiliated to Bacillariophyta in November 2017 and February 2018 at WT and in November 2017 at MH, while also Ciliophora-related OTUs were recorded in high relative abundances (> 40% of the total number of reads) in January 2018, again at MH (Figure 3b). It is noteworthy, that even though no significant differences were found in almost all paired comparisons of OTUs numbers and α-diversity estimators of common sampling dates (Supplementary Table S3), significant differences were only found between WT-MH (with WT < MH), and MH-HB (with MH > HB) concerning the Simpson estimator, and between MH-HB (with MH > HB) concerning the Shannon estimator (Supplementary Table S3) which take into account OTUs abundances. Based on the Bray-Curtis similarities of the unicellular eukaryotic community composition and abundance of the four sites during the common sampling dates three large groupings at a similarity level > 30% were detected (Figure 4). These large groups comprised of seasonal samples, which were further broken into smaller temporal clusters, comprising of the samples of the different sites at the same date, at a similarity level of > 40% (Figure 4). These clusters grouped the samples of the different sites based on the sampling dates, supporting the idea that the temporal scale (rather than the spatial) is the main signal for shaping eukaryotic assemblages. This is in accordance with the results of the ttest paired comparisons of a-diversity indexes that generally did not reveal significant differences between sites (Supplementary Table S3). Two samples (18OctHB and 15NovHB) were clustered together with low similarity in contrast to the general trend of samples of the same date clustering together rather than samples of the same site. The sample of 18OctHB was one of the seven samples with lower number of reads than the rarefication level, and it might have skewed this sample's clustering.

Environmental Effects on Unicellular Eukaryotic Communities
Concerning the relationship between all samples in WT and the corresponding environmental variables, overall, five environmental variables were included in the CCA analysis according to their P-values. The ordination plot produced explained 54.9% of the variance of samples-environmental variables. The first axis (x-axis) was positively related to water temperature and SiO4 and negatively related to salinity and the second axis (y-axis) was positively related to all of the included environmental variables (mainly NO3 − + NO2 − ) except for water temperature. According to the ordination results almost all samples of June-August 2017 were ordinated with a strong positive relation to the water temperature and SiO4 ( Figure 5). The majority of the samples of April-May 2017 were positively related to the included dissolved nutrients and in particular NO3 − + NO2 − and PO4. Regarding the samples of autumn and winter (September-November 2017 and December 2017-February 2018), these were positively related to salinity but negatively to water temperature ( Figure  5). Apart from the separate ordination plots for each sampling site, an overall ordination plot was produced including all samples, which corroborates these results (see Supplementary Figure S5). The nutrients included in all plots did not show any strong and consistent relations among sites. Regarding the relationship of individual OTUs and environmental variables, the OTUs with an overall relative abundance of > 1% of the total number of reads per site were included in the CCA analysis. In all four sites the OTUs 4, 5, 18, 21, and 28 (Table 3), whenever they were included in the plots, were mainly related to the water temperature, the OTUs 1 and 13 (Table 3) were mainly related to salinity, and different OTUs were related to nutrient concentrations in each site (Supplementary Figure S6). It is worth mentioning that OTU_2, being closely affiliated to the red tide forming dinoflagellate Noctiluca scintillans (Table 3), was positively related to NH4 at the WT site (Supplementary Figure S6).
Among the above taxa, OTU_1, OTU_2, OTU_4, OTU_5, OTU_6, and OTU_8 (Table 3) exhibited relative abundances > 10% of a sample's total number of reads, representing arbitrarily bloom forming OTUs, in at least one sample in all four sites (Supplementary Figure S7). In particular, OTU_1, with closest relative the dinoflagellate Scrippsiella trochoidea, followed by OTU_2 (the dinoflagellate Noctiluca scintillans), were detected in relative abundances > 10% of a sample's total number of reads for long periods during the study, including before and during the mucilage aggregate event, especially in WT for OTU_2. It is noteworthy that OTU_4, closely affiliated to Chaetoceros tenuuissimus was frequently detected during the warm months and was strongly positively related to the water temperature, while the OTU_5 (Gonyaulax fragilis) had high abundances during the mucilaginous aggregate phenomenon in June 2017 (Supplementary Figure  S7). Other dominant OTUs were observed in peak abundances in individual samples throughout the study. On the other hand, OTUs closely related to known harmful species (such as the dinoflagellate Dinophysis cf. acuminata, which was found as a rare overall OTU) were found in low abundance (number of reads), suggesting the presence of taxa with potential negative impacts on human health.

Discussion
The study area is considered a heavily modified marine water body according to the Water Framework Directive (WFD, 2000/60/EC), one of the most polluted coastal areas of Greece [50], which is influenced by both anthropogenic activity and freshwater inputs. During the study period N and P concentrations were among the highest recorded for eutrophic coastal areas of the Mediterranean Sea [51,52]. Moreover, environmental variability between the closely located sampling stations was reported as nonsignificant for the majority of the parameters examined [33].
In this nutrient-rich coastal environment, high biodiversity and abundance of marine unicellular eukaryotes has been previously reported based on classical microscopy observations [32,53,54]. The most abundant and diverse planktonic groups in these observations comprised of Bacillariophyta (diatoms) and Dinoflagellata (dinophytes), while less diverse groups included Haptophyta, Cryptophyta, Chlorophyta, Euglenozoa, etc. Evidenced by High-Throughput Sequencing (HTS) applied in the present study, the dominant taxonomic groups of unicellular eukaryotes belonged to the "usual suspects", frequently found in the area, i.e., Dinoflagella, Bacillariophyta, and Ciliophora, but also previously undetected groups were revealed, belonging to Cercozoa, the Marine Alveolates group (MALVs), the Marine Stramenopiles group (MAST), Labyrinthulea, Picozoa, Oomycota, Fungi, and other members of marine unicellular eukaryotes. In WT site, with higher sampling frequency (weekly samplings) a higher diversity was detected, in comparison to the monthly samplings. It is reasonable that with higher sampling efforts, diversity estimates will increase, as rarely occurring taxa have better chances to be detected [55]. Moreover, because of the short generation times and high turnover rates of marine microbes, rare taxa may have low residence and detection times in the marine environment. Accordingly, frequent temporal (e.g., one week sampling interval) and highly resolved spatial samplings can provide with higher diversity coverage for planktonic organisms [56]. Table 3. The OTUs with relative abundance of > 1% of the total number of reads in each sampling site, their putative higher taxonomic affiliation, their closest relative based on BLAST searches against the PR2 database and confirmed against GenBank, the isolation source of the closest relative, and sampling site in which they were detected in > 1% of the total number of reads. Moreover, OTUs with relative abundances > 10% of the total number of reads in at least one sample are indicated with an asterisk (*) and are arbitrarily defined as blooming OTUs.

OTUs
Putative High-Level Taxonomic  HTS tools have been increasingly used in investigations of marine microbial diversity and have similarly revealed high novel and undetected diversity in coastal areas including the eukaryotic taxonomic groups detected in our dataset [5,[28][29][30]. Christaki et al. [3] reported "new players in the succession scene" of the coast of eastern English Channel, and were associated with MALVs, MAST, and Cercozoa, which were also detected in our study. These new players were suggested to participate in complex interactions [7,15], shaping the spatial [7] and temporal unicellular eukaryote assemblages in the area [9]. Especially the rare microbial fraction of these assemblages appeared to be particularly active and play significant ecological roles [57].
During the period of the study, conspicuous, frequent, persistent, and successive phytoplankton blooms were alternated with red tides and mucilage aggregates, comprehensively described in [33]. These phenomena were detected by means of microscopy almost continuously throughout the study period and were of great ecological importance for the coastal system. They, also, had significant socio-economic impact to the residents of the city of Thessaloniki, as in several occasions the red tides and HABs were macroscopically visible impacting the touristic urban front of the city center [33]. The key plankton abettors and perpetrators contributing to the formation of the phytoplankton blooms were the known mucilage-producing diatoms Cylindrotheca closterium, Leptocylindrus spp., Skeletonema costatum, Chaetoceros spp., and the dinoflagellate Gonyaulax cf. fragilis, while the species responsible for the development of the red tides were the dinoflagellates Noctiluca scintillans and its close relative Spatulodinium pseudonoctiluca, and the photosynthetic ciliate Mesodinium rubrum. The characterization of the unicellular eukaryotic communities with HTS tools showed that the dominant biosphere, taking over a high percentage of the reads, was closely related to the abovementioned species. In particular, OTUs closely related to the red tide forming Noctiluca scintillans, the bloom forming Chaetoceros spp., Skeletonema pseudocostatum, Gymnodinium aureolum and Peridinium quinquecorne, and the mucilage aggregates abettor Gonyaulax fragilis were among the most abundant OTUs in all samples (Table 3). Indeed, these OTUs were detected in high relative number of reads during the same periods that their closest relatives were identified by microscopy in [33].
Furthermore, during these phenomena the system was being taken over by the dominant taxa, as evidenced by the hampered rate of newly introduced OTUs in the system during these periods, compared to the nonbloom periods. Among these OTUs, previously undetected taxa in the area were observed, namely the dinoflagellates Adenoides eludens and Islandinium tricingulatum, the haptophyte Haptolina sp., taxa belonging to MAST and the ciliate Parundella sp., all of them with unusual and complex life cycles and trophic preferences [58][59][60] influential in global processes [61], and suggesting alternating trophic states and carbon paths [3,14]. In addition to microscopic data, OTUs closely related to known harmful species were opportunistically found with high read numbers during bloom / red tide phenomena, such as parasitic MALVs [62][63][64] and Cercozoa [65,66], revealing a (previously undetected) significant biodiversity with potentially negative impacts for the unicellular eukaryotic community and the ecosystem. Especially MALVs have been detected in high abundances in virtually all studies of marine unicellular eukaryotes, suggesting their ubiquitous role as parasites of various marine organisms and their significant part as a trophic link in marine systems [63].
The taxonomic spatial variability between the four sampling sites appeared similar by Bray-Curtis similarities, possibly reflecting the temporally homogenous environmental pressures in the general sampling area. To highlight this, our results suggested that temporal succession represented a stronger force rather than sampling location (i.e., small spatial scale), as seasonal samples grouped together by b-diversity, irrespective of the sampling site. Numerous studies have suggested relations between the marine microbial community composition and environmental parameters variability [13,[67][68][69]. Thus, distance between different sampling sites is only important in shaping the microbial assemblages, only when it coincides with variations in environmental pressures. Strong spatial variability has been reported in distant and coastal areas [28,70,71], but also in proximate sampling sites [7], consistently associated in part with different environmental stresses. Here, the environmental variability explained > 54% of the overall community structure for each sampling site during the study period. In particular, temperature and salinity were associated with distinct seasonal clusters, i.e., temperature was positively related to the warm months cluster (June-August) and salinity was mostly related to the colder months, mainly September-February. Genitsaris et al. [9] using canonical analysis showed that environmental variability explained only 30% of the seasonal succession of marine protists in the eastern English Channel, while similarly Chow et al. [72] reported that only 10% of the variability of protistan assemblages at the San Pedro Oceanic site could be predicted according to environmental values. Therefore, biotic relations were suggested to play important roles in shaping marine protistan communities additionally to the environment [7,73]. In the present study, Maximal Information-Based Nonparametric Exploration (MINE) by computing the Maximal Information Coefficient (MIC) among the dominant OTUs, and between dominant OTUs and environmental parameters [74] in each sampling site, was performed in order to explore potential significant correlations among OTUs and between OTUs and environmental variables. However, the analysis showed no indications of significant correlations between these dominant OTUs, and between OTUs and the environmental factors, according to Weiss et al. [75], who proposed to use a P-value of < 0.001 to identify strong relations with the MIC (data not shown). This could suggest that biotic interactions might not contribute highly to the shaping of the marine unicellular eukaryotic assemblages.
Temperature and salinity were the two environmental variables that seemed to influence most the overall unicellular eukaryotic communities, but also individually the majority of the dominant OTUs, as expected. Both water temperature and salinity are among the key components of climate change, which are projected to directly affect marine microbial diversity and functions [76]. Temperature has been found to affect the size [77], metabolic rate [78], and activity [13] of microbes, thus influencing marine protistan community and food web structure [79,80]. Stefanidou et al. [13], using HTS tools in mesocosm salinity manipulations of Thessaloniki Bay unicellular eukaryotic assemblages detected a number of OTUs associated with euryhaline taxa, also found in our dataset (e.g., Skeletonema, Chaetoceros, Prorocentrum, Gyrodinium, and Gymnodinium-related OTUs). They concluded that in this Mediterranean environment, the standing biodiversity can act as a buffer against environmental fluctuations. Beyond the water temperature and salinity nutrient concentrations were only partly associated with specific OTUs. In particular, nitrogen was positively correlated with the OTU closely related to the dinoflagellate Noctiluca scintillans, which was identified as the chief responsible species for the observed red tides during the study. N. scintillans is among the most frequently observed red tide forming dinoflagellates worldwide in eutrophic systems [81,82] and water temperatures ranging from 10 to 25 °C, and salinity range from 28 to 36 [83,84], similar to our study area. The positive connection between N. scintillans read abundance and NH4 in our study area agrees with the reported correlation of its cell abundance and NH4 with microscopy data in the same samples [33] and might indicate nutrient regeneration by this heterotrophic dinoflagellate and contribution to the local nutrient pool. In addition, the role of N. scintillans as a recycler of nitrogen has been previously linked to high concentrations of nitrogen in its cells [85]. Even though microscopy data of the same samples showed strong connections of most nutrients measured in the present study [33], metabarcoding data did not reveal similar trends with the exception of the positive correlation of NH4 and the OTU closely related to N. scintillans.

Conclusions
The high-throughput sequencing tools implemented in this study showed a high diversity of known, frequently occurring protists, also detected with microscopy in previous studies of the area. It further revealed several taxa previously undetected in the area, with seemingly important roles in the structure and functioning of the entire unicellular eukaryotic communities, and key contribution to the phenomena of red tides, harmful algal blooms, and mucilage aggregates observed during the study. The unicellular eukaryotic assemblages showed temporal patterns rather than small-scale spatial separation, possibly responding to the water temperature and salinity seasonal variations. The rich species pool along with the eutrophic status of the area facilitates the unicellular eukaryotic succession, mainly favoring the taxa affected by the established cycles of temperature and salinity. Overall, environmental pressures (eutrophication, nitrogen pollution) under seasonal changes of water temperature and salinity (circulation pattern), are suggested to be the major drivers of the composition changes and the appearance of blooms and red tides in the study area as expected.
Supplementary Materials: The following are available online at www.mdpi.com/1424-2818/12/3/114/s1, Figure  S1. Temporal variation of physical and chemical variables, i.e., temperature (°C), salinity, silicates (SiO4), nitrates and nitrites (NO3 − + NO2 − ), ammonium (NH4), particulate organic phosphates (POP), and Chl a (μg L −1 ) during the study, Figure S2. Rarefaction curves representing the number of OTUs of unicellular eukaryotes against the number of high-quality reads after rarefication and removal of multicellular metazoan and streptophyta affiliated reads, Figure S3. Pie chart showing the relative number of OTUs belonging to major high-level unicellular eukaryotic taxonomic groups, based on searches against the PR2 database; to facilitate reading, the groups with ≥ 1% of the total number of OTUs are shown, Figure S4. Venn diagram, illustrating the number of unique and shared OTUs among the four sampling sites during the common sampling dates, Figure S5. Biplot of environmental parameters (see Figure S1) and all samples, based on OTUs number of reads. Different colors represent different seasons, according to the coloring of groups in Figure 4, Figure S6. Biplots of environmental parameters (see Figure S1) and the dominant OTUs with relative abundance of > 1% of the number of reads in each sampling site per sampling site (instantly abundant), based on OTUs number of reads, Figure S7. Heatmap of the OTUs with relative abundance of > 10% of the sample's total number of reads in at least one sample in the four sites. Columns are mean-centered, with relative number of reads per OTU represented by colour (pink: Lower numbers; green: Higher numbers), Table S1: Number of OTUs, number of reads, the richness estimator (SChao1), and the heterogeneity of the α-diversity estimators (Simpson, Shannon, and Evenness) of all samples at the White Tower sampling site. The total and average number of OTUs, number of reads, and average values of the α-diversity estimators are also shown. N/A: Not applicable, Table S2: Number of OTUs, number of reads, the richness estimator (SChao1), and the heterogeneity of the α-diversity estimators (Simpson, Shannon, and Evenness) of the common sampling dates at White Tower, Aretsou, Music Hall, and Harbour sampling sites. The total and average number of OTUs, number of reads, and average values of the α-diversity estimators are also shown. N/A: Not applicable, Table S3: Differences in OTUs number, and the diversity indexes among the different sites (White Tower (WT); Aretsou Beach (AR); Music Hall coast (MH) and Harbour (HB)) as revealed by paired t-test. (df: Degrees of freedom. (*: P-value < 0.05, **: P-value < 0.01, ***: P-value < 0.001, ns: Nonsignificant; for all cases n = 12). Funding: This research is implemented through the IKY scholarships programme and co-financed by the European Union (European Social Fund-ESF) and Greek national funds through the action entitled "Reinforcement of Postdoctoral Researchers", in the framework of the Operational Programme "Human Resources Development Program, Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) 2014-2020.