Phytoplankton Blooms, Red Tides and Mucilaginous Aggregates in the Urban Thessaloniki Bay, Eastern Mediterranean

We investigated the plankton community composition and abundance in the urban marine environment of Thessaloniki Bay. We collected water samples weekly from March 2017 to February 2018 at the coastal front of Thessaloniki city center and monthly samples from three other inshore sites along the urban front of the bay. During the study period, conspicuous and successive phytoplankton blooms, dominated by known mucilage-producing diatoms alternated with red tide events formed by the dinoflagellates Noctiluca scintillans and Spatulodinium pseudonoctiluca, and an extensive mucilage aggregate phenomenon, which appeared in late June 2017. At least 11 known harmful algae were identified throughout the study, with the increase in the abundance of the known harmful dinoflagellate Dinophysis cf. acuminata occurring in October and November 2017. Finally, a red tide caused by the photosynthetic ciliate Mesodinium rubrum on December 2017 was conspicuous throughout the sampling sites. The above-mentioned harmful blooms and red tides were linked to high nutrient concentrations and eutrophication. This paper provides an overview of eutrophication impacts on the response of the unicellular eukaryotic plankton organisms and their impact on water quality and ecosystem services.


Introduction
On a global scale, the rate of coastal urbanization will increase rapidly in the next decades, and in combination with climate change is projected to result in an increased risk of coastal eutrophication [1,2]. Sewage inputs from coastal cities that are transported directly to coastal waters can act synergistically with land-based sources and river run-off causing high levels of nutrients [3,4]. Consequently, the global Indicator for Coastal Eutrophication Potential (ICEP) analyses indicate that the potential for coastal eutrophication continuously grows worldwide [2]. Worldwide eutrophication has led to phytoplankton abundance and biomass increase [5][6][7], while more coastal harmful algal blooms (HABs), with more toxic species, have been linked with eutrophication phenomena [8,9]. Numerous examples of linkages between nutrient loading and coastal phytoplankton blooms and mucilagine aggregate phenomena [10,11] include the involvement of harmful species, i.e., the diatom Pseudonitzschia spp. in the Gulf of Mexico [12], the dinoflagellates Prorocentrum sp., and Karenia mikimotoi along the coast of China [13], and the red tide-forming heterotrophic dinoflagellate Noctiluca scintillans [14][15][16].  During all samplings, in situ measurements of water temperature and conductivity were made with the use of the YSI Pro 1030 instrument (YSI Inc., Yellow Springs, OH, USA). Conductivity was transformed to salinity based on the equation in Weyl [27]. Water samples of 2 L were collected from the surface layer of 1 m, and separated as follows: (i) a subsample of 0.5 L was used for immediate microscopic observation of the living microbial eukaryotic community; (ii) a subsample of 0.5 L was preserved in Lugol's solution and kept in the dark in room temperature for microscopic analysis within the next few days; (iii) subsamples of 100-250 mL (depending on plankton and particulate matter density) were immediately filtered onto 0.7 μm pre-washed (in 5-10% HCl) and precombusted (6 h, 550 °C) Whatman GF/F filters, and the filters were stored in −20 °C until future particulate organic phosphorus and chlorophyll a (Chl a) measurements; (iv) a subsample of 50 mL was filtered through 0.2 μm cellulose acetate filters (Sartorius) and the filtered water aliquots were kept in −20 °C until future dissolved inorganic nutrient measurements.

Chl a and Nutrient Measurements
Chl a content was estimated according to Jeffrey and Humphrey [28]. Prior to the photochemical measurements (HITACHI, U2900) filters were put into 8 mL acetone (90%) for 24 h in the dark at 6 ℃.
Furthermore, the Eutrophication Index (E.I.) of Primpas et al. [30] was used in order to assess the eutrophication status of Thessaloniki Bay. The formula takes into consideration the NO3 − and NO2 − , ammonium, PO4 − , and Chl α concentrations resulting in three distinct ranges describing oligotrophy (0.04-0.38), mesotrophy (0.37-0.87), and eutrophication (0. 83-1.51). The ranges are further divided into a five-scale scheme according to the WFD requirements, in order to assess the water quality status, as follows: 1. High ecological water quality: < 0. During all samplings, in situ measurements of water temperature and conductivity were made with the use of the YSI Pro 1030 instrument (YSI Inc., Yellow Springs, OH, USA). Conductivity was transformed to salinity based on the equation in Weyl [27]. Water samples of 2 L were collected from the surface layer of 1 m, and separated as follows: (i) a subsample of 0.5 L was used for immediate microscopic observation of the living microbial eukaryotic community; (ii) a subsample of 0.5 L was preserved in Lugol's solution and kept in the dark in room temperature for microscopic analysis within the next few days; (iii) subsamples of 100-250 mL (depending on plankton and particulate matter density) were immediately filtered onto 0.7 µm pre-washed (in 5-10% HCl) and pre-combusted (6 h, 550 • C) Whatman GF/F filters, and the filters were stored in −20 • C until future particulate organic phosphorus and chlorophyll a (Chl a) measurements; (iv) a subsample of 50 mL was filtered through 0.2 µm cellulose acetate filters (Sartorius) and the filtered water aliquots were kept in −20 • C until future dissolved inorganic nutrient measurements.

Chl a and Nutrient Measurements
Chl a content was estimated according to Jeffrey and Humphrey [28]. Prior to the photochemical measurements (HITACHI, U2900) filters were put into 8 mL acetone (90%) for 24 h in the dark at 6 • C.
Particulate organic phosphorus (POP) was measured colorimetrically by an element analyzer (Thermo Scientific Flash 2000) at 882 nm, following the protocol by Hansen and Koroleff [29] Planktonic unicellular eukaryotes were examined in sedimentation chambers using an inverted epi-fluorescence microscope (Nikon Eclipse TE 2000-S, Melville, USA), with phase contrast. Taxa were identified based on taxonomic keys and relevant papers [31][32][33]. Light and phase-contrast micrographs of live and Lugol-preserved cells were taken using a digital microscope camera (Nikon DS-L1, Melville, USA). Plankton counts (cells and colonies) were performed using the inverted microscope method [34]. At least 400 individuals in total, and 100 individuals of the most abundant taxa, were counted per sample in sedimentation chambers. Taxa comprising of > 10% of the total plankton abundance per sample were arbitrarily considered to be dominant for that particular sample. Population density of 1000 cells mL −1 for a particular phytoplankton taxon in a sample was considered as a baseline bloom density in this urban coastal environment. This threshold is based on the Greek eutrophication scale and the total phytoplankton abundance (960 cells mL −1 ) given as an indicator of bad water quality or eutrophic coastal waters [35]. Potentially harmful plankton taxa identified during the study, were acknowledged according to the IOC-UNESCO Taxonomic Reference List of Harmful Microalgae.

Data Analysis
Alpha-diversity estimators (the Simpson, Shannon, Evenness, Equitability, and Berger-Parker indices) were calculated with the PAST 2.17c software [36] in all samples. These indices have been reported to better describe general properties of the communities [37] and reflect anthropogenic or environmental variability effects on ecosystem functions [38]. Paired t-tests were applied in PAST 2.17c software to compare the (i) physical and chemical variables and (ii), richness, abundance, and the alpha-diversity estimators between the four sampling sites. The p-values < 0.05 indicated significant differences between pairwise comparisons. Furthermore, pairwise comparisons of sampling sites, based on the relative abundance of individual taxa, were implemented with the Kolmogorov-Smirnov test in the PAST 2.17c software.
The plankton assemblages of the different samplings were compared using the Plymouth routines in the multivariate ecological research software package PRIMER v.6 [39]. The Jaccard coefficients were calculated to develop the matrix based on taxa abundance in order to identify interrelationships between samples and construct cluster and MDS (multi-dimensional scaling) plots. The similarity profile (SIMPROF) permutation test was conducted to determine the significance of the dendrogram branches resulting from cluster analysis.
Network analysis was performed in order to explore strong relationships among plankton taxa, and between plankton taxa and environmental parameters in all samplings. The relationships were characterized through MINE (maximal information-based nonparametric exploration) statistics by computing the maximal information coefficient (MIC) between each pair of taxa, and pairs of taxa and environmental parameters [40], considering abundance values for each taxon. MIC is a non-parametric method which captures associations (linear or non-linear) between data pairs. It provides a score that represents the strength of the relationship. The matrix of MIC values corresponding to p-values < 0.01, based on pre-computed p-values of various MIC scores at different sample sizes, was used to visualize the networks of associations with Cytoscape 3.5.1 [41]. Furthermore, correlation analysis was conducted in order to investigate the relationships between the plankton taxonomic groups and the E.I., and the inorganic nutrient molar ratios (N:P, Si:N, Si:P).

Environmental Parameters
Seawater temperature recorded in the WT sampling site during the period of the study ranged from 9.60 to 29.7 • C and salinity from 32.8 to 38.8 (Table 2) Considering the common sampling dates conducted in the different sampling sites, no significant differences were found in almost all paired comparisons of environmental parameters, based on t-tests (for a visualization of mean values of environmental parameters in each sampling site, see Supplementary Figure S2). Significant differences were found between sites for NO 2 (with AR > MH), for NO 3 and NO 2 (WT > MH; AR > MH), and for POP (WT > MH; MH < HB) (Supplementary Table S1). Furthermore, higher ammonia (NH 4 ) concentrations were recorded at HB (Table 2), even though no statistically significant differences between sites were found.

Plankton Diversity and Abundance
A total of 117 plankton morphospecies were identified in all four sampling sites during the study period (Supplementary Table S2). The taxonomic group of Bacillariophyta (diatoms) had the highest overall species richness as 44% of the total number of taxa belonged to this group, and was followed by Dinophyta (including also mixotrophic and heterotrophic dinoflagellates) (37% of the total number of taxa), Cryptophyta (5%), Haptophyta (3%), Chlorophyta (2%), Dictyochophyta (2%), and Euglenozoa (2%), while other groups (Cercozoa, Chrysophyceae, Telonemida, Xanthophyceae, and Ciliophora) contributed with < 2% of the total number of taxa (Supplementary Figure S3). In the four sampling sites, dinoflagellates were more diverse in terms of species richness during March 2017-November 2017, while diatoms appeared to be more diverse, in all sites, during December 2017-February 2018 (Supplementary Figure S3). The other taxonomic groups had a more or less consistent representation that altogether did not exceed in any case 40% of the total number of taxa in a sample.
The number of identified taxa varied among samples between 24 (22 March and 19 April in WT, and 9 May 2017 in MH) and 57 taxa on 8 November 2017 in WT (Supplementary Table S3), with the highest values in December-February when the measured water temperature was lower than 15 • C (Figure 2a). High variability was recorded in total cell abundance of phytoplankton reaching a maximum of 42,000 cells mL −1 on 19 July 2017 in WT, dominated by the diatom Skeletonema costatum (see Figure 2b). Heterotrophic dinoflagellates dominated by red tide forming N. scintillans exhibited highest values on 22 March 2017 in WT (>3250 cells mL −1 ). The mean total taxa number and abundance were the only a-diversity estimators that were found significantly different in some paired comparisons between the different sites, based on t-tests; in particular: WT > AR, WT > MH and MH < HB (Supplementary  Table S4; for a visualization of mean values of taxa number and abundance values in each sampling site, see Supplementary Figure S4). However, no significant differences in the distribution of the taxa relative abundances between sites were detected according to the Kolmogorov-Smirnov test (p > 0.05). The other a-diversity estimators calculated (Simpson, Shannon, Equitability, Evenness, and Berger-Parker), fluctuated during the study, and showed sometimes relatively low values, reflecting high dominance by one (or few) taxa, and high variation between taxa abundances within the community. In particular, the sampling dates with the low Simpson index  Table S3).  Table S4; for a visualization of mean values of taxa number and abundance values in each sampling site, see Supplementary Figure S4). However, no significant differences in the distribution of the taxa relative abundances between sites were detected according to the Kolmogorov-Smirnov test (p > 0.05). The other a-diversity estimators calculated (Simpson, Shannon, Equitability, Evenness, and Berger-Parker), fluctuated during the study, and showed sometimes relatively low values, reflecting high dominance by one (or few) taxa, and high variation between taxa abundances within the community. In particular, the sampling dates with the low Simpson index  Table S3).

Phytoplankton Blooms, Red Tides, and a Mucilage Aggregate Phenomenon
Based on the plankton community composition and abundance during the study period, four major clusters were identified at a similarity level ~ 35% (Figure 3) grouping together the samplings irrespectively of the sample collection site, according to the sampling dates: March-June 2017 (Cluster I); July-October 2017 (Cluster II); November 2017 (Cluster III); and December 2017-February 2018 (Cluster IV). This is in accordance with the results of the t-test paired comparisons of a-diversity indices between sites that showed no significant differences in most occasions. The samples in each cluster were further grouped together based on higher similarity levels (> 40% similarity). These groupings included a small number of samples taken in close dates and were characterized by phytoplankton blooms (> 1000 cells mL −1 ) of a taxonomic group or a single species, or/and red tides ( Figure 3).

Phytoplankton Blooms, Red Tides, and a Mucilage Aggregate Phenomenon
Based on the plankton community composition and abundance during the study period, four major clusters were identified at a similarity level~35% (Figure 3 IV). This is in accordance with the results of the t-test paired comparisons of a-diversity indices between sites that showed no significant differences in most occasions. The samples in each cluster were further grouped together based on higher similarity levels (>40% similarity). These groupings included a small number of samples taken in close dates and were characterized by phytoplankton blooms (>1000 cells mL −1 ) of a taxonomic group or a single species, or/and red tides (Figure 3). On the other hand, conspicuous red tides, macroscopically visible, appeared in the front of the Bay at three occasions, during this period, making the water viscous. The red tides were detected at 22 March, 12 April, and 14-21 June 2017, and mainly consisted of the known red tide forming dinoflagellates N. scintillans and its close relative Spatulodinium pseudonoctiluca. Especially on 22 March 2017, the event was so intense that the sample consisted entirely of N. scintillans cells, reaching 3250 cells mL −1 , comprising > 99% of the total abundance. The co-occurrence of these species with bloom-forming, mucilage-producing diatoms, e.g., Cylindrotheca closterium, Chaetoceros spp., L. minimus, L. danicus, S. costatum, the haptophyte Phaeocystis sp. and the dinoflagellate Gonyaulax cf. fragilis, were observed. They were producing or being embedded in mucilage, before and during the development of an extreme aggregation of mucilage, between 28 June and 4 July 2017. N. scintillans was observed to feed on diatoms, and most commonly on Chaetoceros spp.
During the period July-October, diatoms remained in high numbers, dominated by the taxa Chaetoceros spp. (max: > 6000 cells mL −1 on 19 July), and more rarely C. closterium (max: > 1800 cells On the other hand, conspicuous red tides, macroscopically visible, appeared in the front of the Bay at three occasions, during this period, making the water viscous. The red tides were detected at 22 March, 12 April, and 14-21 June 2017, and mainly consisted of the known red tide forming dinoflagellates N. scintillans and its close relative Spatulodinium pseudonoctiluca. Especially on 22 March 2017, the event was so intense that the sample consisted entirely of N. scintillans cells, reaching 3250 cells mL −1 , comprising > 99% of the total abundance. The co-occurrence of these species with bloom-forming, mucilage-producing diatoms, e.g., Cylindrotheca closterium, Chaetoceros spp., L. minimus, L. danicus, S. costatum, the haptophyte Phaeocystis sp. and the dinoflagellate Gonyaulax cf. fragilis, were observed. They were producing or being embedded in mucilage, before and during the development of an extreme aggregation of mucilage, between 28 June and 4 July 2017. N. scintillans was observed to feed on diatoms, and most commonly on Chaetoceros spp. During the period July-October, diatoms remained in high numbers, dominated by the taxa Chaetoceros spp. (max: >6000 cells mL −1 on 19 July), and more rarely C. closterium (max: >1800 cells mL −1 on 20 September). On 19 July and 13 September 2017, the taxon S. costatum, additionally showed high abundances reaching > 25,000 and > 1500 cells mL −1 , respectively.
The diatom bloom was followed by an increase in the abundance of the harmful dinoflagellate Dinophysis cf. acuminata, during November 2017 (Figure 3 Based on the IOC-UNESCO Taxonomic Reference List of Harmful Microalgae it can be stated that at least 11 out of 117 plankton taxa found in the present study have been reported as harmful. These taxa are the diatoms Pseudonitzschia cf. delicatissima, Pseudonitzschia cf. multistriata, Pseudonitzschia cf. pseudodelicatissima, and Pseudonitzschia cf. pungens, the dictyochophycean Vicicitus globosus, the haptophyte Phaeocystis sp. and the dinoflagellates D. cf. acuminata, Dinophysis caudata, Karenia brevis, Karlodinium spp., and the epiphytic Prorocentrum cf. lima. In particular, the diatoms P. cf. delicatissima, P. cf. pseudodelicatissima, and P. cf. pungens were detected in concentrations > 500 cells mL −1 at the White Tower (WT) site on 19 July 2017, just after the mucilage aggregate phenomenon. Relatively high abundances (270 cells mL −1 ) were recorded for K. brevis at WT on the 28 June, right in the middle of the mucilage aggregate phenomenon. An extremely high bloom (>10,500 cells mL −1 ) was observed at the same sampling point one year later (unpublished data).

Links of Environmental Parameters and Plankton Bloom, Red Tide, and Mucilaginous Aggregate-Forming Taxa
Connections between all detected taxa, including all phytoplankters and red tide/bloom-forming taxa, were investigated according to the MIC correlation coefficient. Only the strong connections between phytoplankters/red tide forming species and environmental parameters were visualized in network analysis ( Figure 4). The strong connections represented MIC values corresponding to pre-calculated p-values (with p < 0.01), based on the total number of samples. Network analysis showed negative connections between salinity/water temperature and the majority of diatom taxa included in the network, and all Cryptophyta and Dictyochophyceae. However, the diatoms Chaetoceros spp., S. costatum, and L. minimus, all mucilage producers, were positively connected with salinity/water temperature. Dinophyta, on the other hand, showed mostly positive strong connections with salinity/water temperature and negative with POP ( Figure 4). To note, the red tide forming N. scintillans showed strong positive connections with NH 4 and PO 4 , while the other red tide forming species, M. rubrum, exhibited negative connections with salinity and water temperature. Finally, the harmful alga D. cf. acuminata displayed negative connection with water temperature (Figure 4).

Discussion
The reason for our focus on plankton community weekly dynamics of the urban marine system in the Thessaloniki Bay's front was motivated by the lack of relevant data in this particular system in eutrophication studies of Thermaikos Gulf, despite the recurrent phenomena of harmful algal blooms (HABs) and conspicuous mucilage phenomena. These phenomena are of great ecological importance for the coastal system and have significant socio-economic impact to the city's residents. After discussing nutrient pressure in the system, species diversity will be discussed, dominance of blooms and red tide forming species, and the key species which were the cause of the mucilage phenomenon verifying results of the marine eutrophication research and the related eutrophication symptoms [26].

Environmental Conditions
In the study area, a heavily modified marine water body according to WFD, annual mean salinity was 37.2 and close to the highest threshold value (37.5) for type IIA of the Mediterranean coastal water types that have been intercalibrated (applicable for phytoplankton) according to Commission Decision 2018/229/UE. This type of coastal water is considered moderately influenced by freshwater inputs, while the annual salinity average value is close to the boundary value (37.5) for type IIIE. Phytoplankton metrics have been intercalibrated only for type IIIE in Greece [35].
Throughout the study, nutrients (N and P) which are indication of eutrophication exhibited high values and were among the highest values reported for nutrient-rich coastal areas of the

Discussion
The reason for our focus on plankton community weekly dynamics of the urban marine system in the Thessaloniki Bay's front was motivated by the lack of relevant data in this particular system in eutrophication studies of Thermaikos Gulf, despite the recurrent phenomena of harmful algal blooms (HABs) and conspicuous mucilage phenomena. These phenomena are of great ecological importance for the coastal system and have significant socio-economic impact to the city's residents. After discussing nutrient pressure in the system, species diversity will be discussed, dominance of blooms and red tide forming species, and the key species which were the cause of the mucilage phenomenon verifying results of the marine eutrophication research and the related eutrophication symptoms [26].

Environmental Conditions
In the study area, a heavily modified marine water body according to WFD, annual mean salinity was 37.2 and close to the highest threshold value (37.5) for type IIA of the Mediterranean coastal water types that have been intercalibrated (applicable for phytoplankton) according to Commission Decision 2018/229/UE. This type of coastal water is considered moderately influenced by freshwater inputs, while the annual salinity average value is close to the boundary value (37.5) for type IIIE. Phytoplankton metrics have been intercalibrated only for type IIIE in Greece [35].
Throughout the study, nutrients (N and P) which are indication of eutrophication exhibited high values and were among the highest values reported for nutrient-rich coastal areas of the Mediterranean Sea [42][43][44]. In a recent comprehensive study on various coastal areas of Greece (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) influenced by anthropogenic activities, mainly by sewage and riverine outflows [45], the Thessaloniki Bay which is one of the most polluted coastal areas of Greece, exhibited the highest PO 4 (max 6.50 µmol L −1 ) and NH 4 concentrations (max 15 µmol L −1 ). Comparing to the highest values of nutrients during 1995-2007 [45], the values in this study (max 9.50 µmol L −1 and 160 µmol L −1 for PO 4 and NH 4 , respectively) were even higher [45]. On the other hand, the highest NO 3 value (21.09 µmol L −1 ) in the urban front of Thessaloniki Bay during our survey was slightly lower than the highest measured NO 3 value (23.5 µmol L −1 ) during the period April-May 2012 [46]. Even so, both of them are extremely high for coastal sites, and are indicative of nitrogen pollution due to anthropogenic activities [47]. These outliers can be used as a sensitive tool for assessing water quality in coastal management studies according to Karydis [48], who showed that outliers are more sensitive in characterizing pollution/eutrophication levels than whole datasets, which usually include a large number of low values. The average N:P ratio in this study was higher (mean 25) compared to the N:P ratio (6.40) of the period 1995-2007 [45]. The high N:P ratio during the present study, in combination with the extreme NH 4 concentrations, may be linked to the relatively high contribution of dinoflagellates, such as N. scintillans and S. pseudonoctiluca, to plankton community biomass [49]. The much higher NH 4 values in 2017-2018 in combination with the high N:P ratio indicated nitrogen pollution, an increasing global problem [50]. Agriculture is the largest source of nitrogen pollution to many of the planet's coastal marine ecosystems [51].
In addition to high nutrient concentrations and ratios, the annual mean values of Chl a for all stations were higher than the values measured in the same area in 2012 [47]. Based on our data, ecological water quality can be classified as bad according to Simboura et al. [47]. Low Chl a values coincided both with low and high cell abundances in WT and HB. In most samples of WT the low Chl a coincided with high cell densities in spring and summer under high irradiance/day-length and temperature when L. danicus and L. minimus were the dominant species. These results may reflect the physiological state of these diatoms on their Chl a content due to the effect of temperature and irradiance [52,53]. Similarly, low Chl a was measured simultaneously with high Leptocylindrus densities in HB.
In addition to the evaluation of the eutrophication and the nitrogen pollution of the study area, based on the individual nutrient variations and their extreme values (outliers), the multimetric eutrophication index E.I. [30] is of great interest for coastal management. In our samples, the E.I. values were always > 0.83 (mean 2.56 after excluding outliers) indicating a heavily eutrophic system reflecting bad environmental status according to Pavlidou et al. [46]. The poor to bad water quality of Thessaloniki Bay according to the phytoplankton-based indices and the E.I. index used, is indicative of both nitrogen and phosphorus enrichment. There is evidence for Greek coastal waters that phytoplankton-based indices are highly sensitive to nitrogen enrichment while the E.I. index is highly sensitive to phosphorus enrichment [43]. It is noteworthy that according to Pavlidou et al. [46] the E.I. reflected the integral eutrophication status of a water body as a whole and has been proposed as a reliable tool regarding the assessment of eutrophication status, and the implementation of nutrient management strategies under the EU WFD and the EU MSFD.

Diversity and Composition of the Plankton Community
Various a-diversity indices have been used (Shannon, Simpson, Equitability, Evenness, Berger-Parker) to describe the structure of the community in terms of its species diversity, dominance and evenness. The species pool of the unicellular eukaryotic plankton community reflected by the a-diversity indices [54] was found similar in the four sampling sites, according to pairwise comparisons with t-test (see Table S4 for a-diversity pairwise comparisons). Additionally, the Kolmogorov-Smirnov test showed no significant differences on the distribution of taxa between sampling sites. A seed bank of the local pool, persistent or transient according to Partel et al. [54] and the plankton life history traits [55] contributed to maintain relatively high biodiversity in this urban degraded marine environment. Even though a-diversity indices showed no significant differences between sites, when considering only the mean number of species identified in this study, significantly lower mean values were found in HB, a site with the highest ammonia and nitrite nitrogen pollution (Table 2), and a generally stressed area because of the harbor daily activity. The lower species diversity might be explained by environmental change consisting of several stressors, which can cause stress-induced community sensitivity [56]. The impacts of environmental stress on biodiversity are well known [57].
Diatoms were the most diverse taxonomic group with the highest species numbers in all sites during the period of December 2017-February 2018 in contrast to dinoflagellates that were more diverse during the period of March 2017-November 2017 (Supplementary Figure S3). The different temporal pattern of diatom and dinoflagellate species richness, as also shown by their contrasting relationships to water temperature and salinity, might be explained by their different response to vertical water mixing versus stratification conditions [58]. However, very few diatoms are strictly restricted to the periods of deep mixing, while sinking is an important factor for the growth and bloom formation only for large, sinking diatoms and large, non-sinking dinoflagellates [55].

Phytoplankton Blooms
During the study, at least one taxon per sample was recorded in bloom abundances, dominating the plankton community. The taxa that were detected in bloom abundances mainly belonged to diatoms, mainly Chaetoceros spp., C. closterium, L. minimus, L. danicus, and S. costatum. The persistent growth and bloom formations by these diatoms under various turbulent and stratification conditions can be explained by their small cell size in combination with mucilage production [55,58,59].
These species with extended blooms were the most important constituents of Thermaikos Gulf phytoplankton 30 years ago, when untreated sewage entered the Gulf [18]. Even though no connections were found between them and nutrient concentrations according to network analysis in the present study, it is well known that under P-limited conditions, certain diatoms become increasingly dominant with increasing Si:P ratios [60]. The most persistent diatoms blooms during our survey coincided with the highest Si:P ratios (>20). Dense diatom blooms in marine ecosystems suffering from eutrophication can generate highly dominant diatom communities within phytoplankton assemblages [61]. Nevertheless, apart from the proved impacts of nutrient concentration and ratios on the occurrence of algal blooms [62], in many cases it seems that algal bloom proliferation is more complicated, and the quantity and the ratio of inorganic nutrients alone cannot sufficiently explain high abundance blooms of extended duration [8].
The known harmful species D. cf. acuminata, P. cf. delicatissima, P. cf. multistriata, P. cf. pseudodelicatissima, P. cf. pungens, D. caudata, K. brevis, Karlodinium spp., and the epiphytic P. cf. lima, with worldwide distribution, were detected in relatively high abundances, but not exceeding bloom densities during individual sampling dates. Furthermore, the known harmful alga V. globosus was recorded occasionally in live water samples taken during the sampling period, but its cells usually could not be preserved with Lugol's solution. These plankton species have been previously reported in the Thermaikos Gulf [18,21,63] indicating a persistent seed bank of the local pool [64]. Previous studies reported evidence for a diverse cyst bank, with high recorded abundance of cysts even in periods when the corresponding species were absent from the water column [65]. These cysts were associated with the formation of dense algal blooms in the water column and a high risk of HABs, as could be the case of the Dinophysis bloom observed in the present study during the period October-November 2017, and a short-term excessive Karenia bloom of extreme densities that was observed in spring 2018 (unpublished data). The urban Thessaloniki Bay exhibits the sustained increases in algal blooms and in HABs in accordance with high nutrient levels, similar to reports in other coastal areas of Mediterranean Sea and the Black Sea [23,66].

Red Tides
Several occasions of macroscopically visible red tides were documented over a temperature range of 10 to 25 • C, and a salinity range of 36 to 38.5. They were attributed to the known red tide forming dinoflagellates of N. scintillans together with its close taxonomic relative S. pseudonoctiluca, and the photosynthetic ciliate M. rubrum.
Noctiluca scintillans is one of the most important red tide forming dinoflagellate worldwide in the water temperature range of 10-25 • C, and salinity range of 28 to 36 in eutrophic areas dominated by diatoms [15,67], similar to our study area. Noctiluca red tides have been linked to eutrophication in several areas of the world and especially in the Black Sea, the Sea of Marmara [68][69][70], the Aegean Sea [71], and Adriatic Sea [72]. In contrast, S. pseudonoctiluca has been reported rarely, although new records from many areas suggest a cosmopolitan distribution. Its distribution has been underestimated due to its complex life cycle, morphological variability, and taxonomic issues in its identification [73]. In the Mediterranean Sea among these red tide forming species, another species of Spatulodinium has been found based on DNA analysis [73]. In the Mediterranean Sea, Spatulodinium showes a wide range of temperature preference, similar to the temperature range (19-30 • C) reported in the Mexican Pacific [74]. Both N. scintillans and S. pseudonoctiluca have been considered exclusively heterotrophic and inclusion of diatoms, dinoflagellates, and dictyochophytes have been observed in their cells [74].
The red tide forming N. scintillans terminated its growth in our study area by the increase of water temperature above 25 • C, as in many other studies, but temperature did not correlate with the start of its growth ( [67] and references therein). A rich food supply of a broad spectrum of food items (from bacteria to fish eggs) is needed to start massive growth and formation of red tides, while availability of phytoplankton as a prey is a key factor [67,75]. Particularly, Noctiluca red tides are known to coincide or follow diatom blooms [16,67,76]. A strong temporal overlapping of N. scintillans and diatoms blooms has been also observed in the present study. High numbers of N. scintillans (>400 cells mL −1 ) coincided with high numbers of Chaetoceros spp., L. minimus, S. costatum, and C. closterium cells. Different species of diatoms (mostly Chaetoceros spp.) have been observed in food vacuoles of N. scintillans in agreement with other studies [77,78]. Additionally, N. scintillans was feeding on harmful Dinophysis spp. Noctiluca scintillans containing toxigenic Dinophysis and Pseudonitzschia species may act as a vector of toxigenic algae to higher trophic levels or transport to shellfish aquaculture [79]. On the other hand, grazing pressure by N. scintillans on the growth of other toxigenic dinoflagellates should be considered as a potential regulator of phytoplankton toxins production [80]. This is of particularly interest in our study area, due to its close vicinity to the biggest mussel culture of Greece, where a harvest ban is often implemented due to Dinophysis spp. abundance > 1 cell mL −1 [81].
Accumulation of N. scintillans cells in the surface water, forming a red tide, was observed under calm weather days (generally with daily mean wind speed < 3 m s −1 ) in this urban front of the bay protected from intense water circulation [82]. It is established that meteorological conditions and topography are crucial factors for red tide formation [75]. In Thessaloniki Bay, N. scintillans appeared to prefer higher salinity (>36) relative to those found in other studies [16,67,70]. Based on our results and N. scintillans abundance dynamics during spring-early summer in the Black Sea and the Northern Adriatic Sea, it is suggested that weather forecasts, and in particular wind speed projections, can be used for medium-term prediction of red tides [83].
A strong positive connection between N. scintillans cell abundance and NH 4 and PO 4 in our study area might indicate nutrient regeneration by this heterotrophic dinoflagellate and contribution to the local nutrient pool. The significant role of N. scintillans as a nutrient regenerator and an efficient recycler of nitrogen has been linked to extremely high concentrations of nitrogen in its cells and excretion regulated by nutrient quality of its food items. Nutrient liberation of senescent cells would stimulate the phytoplankton growth near the red tide patches while improving the food quality for N. scintillans [78]. NH 4 regeneration by N. scintillans in coastal seas has been reported by Montani et al. [84] whereas high ammonia concentrations released from Noctiluca cells during the decay process of the red tide were also shown by Schaumann et al. [85]. NH 4 also increases during decline bloom phase indicating release of intracellular NH 4 accumulated through Noctiluca grazing according to Baliarsingh et al. [86]. Direct toxicity to fish by ammonium/ammonia is possible although at seawater pH, approximately 5% of total ammonia is unionized NH 3 [87].
Mesodinium rubrum is a globally distributed photosynthetic ciliate that sometimes causes red tides in coastal waters [88]. M. rubrum is a marine plankton of great cytological, physiological, and evolutionary interest, which has an exceptional type of cellular organization not realized by other species, supported by organelle robbery [89]. M. rubrum and its accompanying cryptophytes showed a strong positive correlation (p < 0.001, r = 0.59), according to correlation analysis in the present study. M. rubrum reached high numbers (>700 cells mL −1 ) in December (temperature range from 11.1 to 13.3 • C, and salinity range from 37.3 to 38.0), just after the drop of the D. cf. acuminata maxima at all sites. The red tide formed was spatially extended and the abundance of M. rubrum was found negatively correlated with both temperature and salinity. An important factor for M. rubrum seasonal dynamics and its short-lived bloom in the study area seems to be the persistent occurrence of Dinophysis spp. and several common heterotrophic dinoflagellates, which it is known to feed [90,91].

Mucilage Aggregates
Noctiluca red tides have been linked to mucilage phenomena, either as a shift from red tides to these events [92] or an overlapping that could be observed in Lapseki coastal area of the Dardanelles in early summer where gelatinous surface layers were recorded [16]. In Thessaloniki Bay, mucilage aggregates appeared on 22 June 2017, characterized by creamy whitish-brownish and gelatinous surface layers [93], which became progressively darker with age. Before the appearance of the phenomenon in Thessaloniki Bay, the plankton community consisted of known mucilage producing species such as the common diatoms in the bay C. closterium, L. minimus, L. danicus, S. costatum, the dinoflagellate G. cf. fragilis and the slime producing red tide dinoflagellates N. scintillans, and S. pseudonoctiluca and the foam forming Phaeocystis sp. Many studies have been published that relate diatom extracellular polymer production with the well-known phenomenon of marine mucilage in the Adriatic Sea, in particular the diatom species C. closterium [94]. This species has been regularly observed as dominant in the mucilage macroaggregates and it has been demonstrated experimentally that polysaccharides from it can form a gel network similar to the macroscopic gel phase occurring in the northern Adriatic Sea irrespective of any bacterial mediation or interaction with inorganic particles [95]. In our study, C. closterium has been observed abundant before, during and after the mucilage phenomenon, within abundant transparent freshly formed mucilage, whereas Chaetoceros spp. and S. costatum chains were also embedded in the mucilage. The dinoflagellate G. cf. fragilis was observed actively producing mucilage in the samples from June 2017 similarly as in the Emilia-Romagna coast (Northern Adriatic Sea) by Pompei et al. [96]. The same phytoplankton mucilage producers, i.e., G. cf. fragilis, S. costatum, and C. closterium were identified as abundant species also in a mucilage phenomenon in the Sea of Marmara [97]. The well-known foam-forming Phaeocystis pouchetii caused mucilage problems in the Evoikos Gulf [98].
Small mucilage aggregates were observed both by active and decaying N. scintillans cells during red tide formation and termination. Decaying N. scintillans cells contribute high amounts of organic matter to the local pool while active cells excrete mucus for trapping food items [78]. The aggregates formed by decaying N. scintillans sampled in the Northern Adriatic Sea presented a similar chemicalbiochemical composition to the different typologies of mucilage aggregates in the same area [94]. This showed that the organic matter of N. scintillans could form a part of the mucilage organic matter in the Adriatic Sea. The accumulation of excess autochthonous organic material (dead and alive material) from the preceding red tides and phytoplankton blooms producing mucilage (from late March to June) and the mucilage they produce, in combination with the hydrodynamic conditions in the Bay (initiation of thermal stratification in May) [82] are suggested as main factors for the formation of the creamy and gelatinous surface layer in the urban Thessaloniki Bay. Our results agree with Umani et al. [10] study on the microbial community of a coastal area in the northern Adriatic Sea with frequent reports of mucilage aggregates, suggesting that mucilage is derived from accumulated slow-to-degrade dissolved organic matter. The months preceding the mucilage events (March-May) in the Northern Adriatic Sea were assumed to be an 'incubation' period. Mucilage was the consequence of a coupling between the accumulation of organic matter and the temporal pattern of meteorological and oceanographic conditions [99] similar to our observations. Strong north winds (>10 m s −1 ) in the beginning of July 2017 were successful to degrade the gelatinous surface layer suddenly and disperse it as small mucilage aggregates. After a week, microscopic aggregates were observed concentrated above the pycnocline (5 m) in deeper areas in the Thermaikos Gulf [82].

Conclusions
During the study period, analysis of the weekly water samples from the urban coastal frontal zone of Thessaloniki Bay (Thermaikos Gulf) provided an outlook of the effects of eutrophication in this Mediterranean urban environment with further implications on marine eutrophication research and coastal management. In the majority of the samples, phytoplankton abundance, and nutrient concentrations indicated high eutrophic conditions and bad environmental status according to the implementation of the EU WFD and the EU MSFD. In addition, in all samples, the Eutrophication Index (E.I.) indicated a heavily eutrophic system, which was characterized by persistent phytoplankton blooms and conspicuous red tides. The phytoplankton blooms were dominated by the diatom species Cylindrotheca closterium, Chaetoceros spp., Leptocylindrus minimus, Leptocylindrus danicus, and Skeletonema costatum reaching high abundances during the spring-summer 2017, while the species Chaetoceros tenuissimus and S. costatum formed blooms during January-February 2018. Red tides of the species Noctiluca scintillans accompanied with Spatulodinium pseudonoctiluca in March 2017, and the species Mesodinium rubrum in December 2017 were observed in the Bay, while a mucilage aggregate phenomenon formed by the mucilage-producers C. closterium, Chaetoceros spp., L. minimus, L. danicus, S. costatum, Phaeocystis sp., and Gonyaulax cf. fragilis was observed in June 2017. These mucilage producers were linked to high temperatures/low salinity, while, on the other hand, red tide forming N. scintillans was linked to high nitrogen and phosphorus concentrations and higher salinities. These harmful events, along with the occurrence of several harmful algae, such as the known toxin-producer Dinophysis cf. acuminata, illustrate the need for continuous monitoring of target indicators of nutrient pollution, ecological water quality and environmental status in the Bay. In the prism of climate change and the increase of eutrophication conditions in coastal areas, this study sounds the alarm and highlights the need to reduce the causes contributing to the bad environmental status and the development of the described phenomena causing severe socio-economic impacts to the public.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/11/8/136/s1; Figure S1. Environmental parameters examined in the four sampling sites; Figure S2. Annual mean values of the environmental parameters examined in the four sampling sites; Figure S3. Pie chart: Relative number of taxa belonging to major high-level unicellular eukaryotic taxonomic groups in all samplings; Figure S4. Bars showing the mean values of the number of taxa and abundance (cells mL −1 ) in the four sampling sites; Table S1. Differences in physical and chemical parameters among the different sites; Table S2. Species list of unicellular eukaryotes found in Thessaloniki's Bay during the study period; Table S3. Species number, total cell abundance, and alpha diversity measurements per sample; Table S4. Differences in taxa number, cell abundance, and the diversity indices among the different sites,

Conflicts of Interest:
The authors declare no conflict of interest.