**Biodiversity of Marine Microbes**

Editor

**Savvas Genitsaris**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Savvas Genitsaris National and Kapodistrian University of Athens Greece

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Diversity* (ISSN 1424-2818) (available at: http://www.mdpi.com).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1052-1 (Hbk) ISBN 978-3-0365-1053-8 (PDF)**

Cover image courtesy of Maria Moustaka-Gouni.

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

**Savvas Genitsaris** has a Ph.D. in Biodiversity and Molecular Ecology from Aristotle University of Thessaloniki, Greece. He worked as a post-doctoral researcher in the University of Littoral (ULCO) in France, where he focused his research on the investigation of biotic and abiotic interactions of marine microbes with the use of metagenomics and bioinformatics tools. He was a post-doctoral researcher of the International Hellenic University and in Aristotle University of Thessaloniki, where he worked on coastal algal blooms and red tides. Currently he is an Assistant Professor in the National and Kapodistrian University of Athens, Greece, where his research focuses on the diversity and function of aquatic microbial communities of the three Domains of Life. His interests include the exploration of the biodiversity of microorganisms in freshwater, marine, coastal and brackish systems, using microscopy, metagenomics, DNA and RNA sequencing, bioinformatics and ecological tools. He is also interested in examining the biotic and abiotic relationships that shape aquatic microbial communities.

## **Preface to "Biodiversity of Marine Microbes"**

Marine microbial life is comprised of a variety of different evolutionary groups from all three domains of life, Eukaryotes, Bacteria, and Archaea. It is responsible for about half of the primary production on earth, plays irreplaceable roles in biogeochemical cycles and ecosystem functioning, and actively participates in complex processes and interactions. Marine microbes are the basis of marine trophic webs (autotrophs), and an important link between different trophic levels, as decomposers, parasites, and endosymbionts. They are used as biological indicators of water quality, eutrophication, and degraded marine environments, and are targeted in conservation and restoration plans. Our understanding of their responses to climate change is considered a key research field to comprehend the complex ongoing processes that will shape the planet's future. They consist of a vast diversity of organisms, with diverse morphological features, sizes, physiology, functions, trophic characteristics, distribution, ecology, evolutionary traits, genetic content, and responses to abiotic variability. Even though we understand their significance in numerous aspects that affect all life on earth, we still have a long way to go to answer fundamental questions driving recent research: How many marine microbes are there? Where can we find them? What do they do? What is their role and responses in light of climate change? What are their phylogenetic relationships? How do they respond to environmental pressures? This along with many more. Advances in high-throughput sequencing accompanied with the technological innovations of classical tools, such as microscopy, have given researchers the equipment and the incentive to attempt to tackle the above questions and shed light on the complex and diverse life of marine microbes. This book provides a platform to highlight new research and significant advances related to the biodiversity of marine microbes.

> **Savvas Genitsaris** *Editor*

### *Editorial* **Biodiversity of Marine Microbes**

#### **Savvas Genitsaris**

School of Science and Technology, International Hellenic University, 57001 Thermi, Greece; s.genitsaris@ihu.edu.gr

Received: 15 June 2020; Accepted: 15 June 2020; Published: 16 June 2020

**Abstract:** The Special Issue entitled "Biodiversity of Marine Microbes" aimed at highlighting the significance of marine microbes as primary producers, their participation in complex processes and interactions with both the biotic and the abiotic environment, and their important roles in biogeochemical cycles and ecosystem functioning. The issue includes five research papers, covering the diversity and composition of marine microbial communities representing all three domains of life in various marine environments, including coastal eutrophic areas, ice waters, and lagoons. One paper examined the diversity and succession of bacterial and archaeal communities from coastal waters in mesocosm experiments. The combination of classical tools with novel technological advances offers the opportunity to answer fundamental questions and shed light on the complex and diverse life of marine microbes.

**Keywords:** algae; prokaryotes; microeukaryotes; coastal eutrophic systems; high-throughput sequencing; blooms; climate change

Marine microbial communities comprise a vast diversity of different evolutionary groups from all three domains of life. Recent technological advances in high-throughput sequencing combined with innovations in classical tools, such as microscopy, have permitted the in-depth examination of these communities, and have advanced our knowledge on the key roles that these microbes may play in marine systems. It is well established that marine algae contribute half of the planet's primary production, microeukaryotes are key players as top-down heterotrophs on marine trophic webs, and bacteria and archaea are essential links in global biogeochemical cycles and marine ecosystem functioning. However, many questions remain underexplored, which are relevant, for example, to the response of coastal protistan communities in high eutrophication conditions, the diversity of picoeukaryotes in extreme environments such as ice waters, the phytoplankton composition along variable conductivity gradients, the effects of hydrography on the marine microbial communities' structure, and the changes in picoplankton under different light conditions and bloom events.

This Special Issue comprises five research articles attempting to tackle the above-mentioned subjects of the biodiversity of marine microbes. Genitsaris et al. investigated the plankton community's composition and abundance in an urban eutrophic coastal area of the Mediterranean, amid frequent and persistent phytoplankton blooms, red tide events, mucilaginous aggregates, and the proliferation of potentially harmful species. This paper provided evidence of the eutrophication effects on the response of eukaryotic plankton assemblages and their impact on water quality and ecosystem services [1]. A metagenomic study of the under-ice photosynthetic picoeukaryotes of the White Sea basin, by means of 18S rRNA gene high-throughput sequencing, indicated that environmental variability in extreme marine environments could lead to a significant shift in microbial communities' composition and structure [2]. In addition, a study of phytoplankton diversity patterns along salinity variations in Mediterranean lagoons suggested that the communities' heterogeneity was highly associated with the recorded differences in conductivity among the different sites [3]. Similarly, Gong et al. showed that spatial hydrographic variability in Taiwan Strait lead to significant variations in all the components

of the microbial food web, including viruses, picoplankton, nanoflagellates and ciliates, and showed that, during the cold months, a "viral loop" might contribute to the top-down control of bacterial populations [4]. Finally, one paper in the Special Issue concerned a mesocosm experimental set-up, aiming to investigate bacterial and archaeal diversity and succession in changing light regimes during phytoplankton growth. This study confirmed that light irradiance can affect marine bacterial communities' structure both directly and indirectly, which in turn can have significant implications for a marine ecosystem's response to environmental change [5].

In my view, and according to the published information in this Special Issue, the main challenge of the coming years and the main goal of the future marine microbiologists, is to successfully combine novel technological advances and integrate the recent sequencing breakthroughs with the classical analytical tools, in order to reveal the true vast diversity of marine microbes and to understand the structure of these communities in relation to abiotic disturbances, environmental pressures and biotic relations.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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

#### **Savvas Genitsaris 1,2, Natassa Stefanidou 1, Ulrich Sommer <sup>3</sup> and Maria Moustaka-Gouni 1,\***


Received: 5 July 2019; Accepted: 12 August 2019; Published: 14 August 2019

**Abstract:** 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.

**Keywords:** nutrients; HABs; mucilaginous aggregates; *Noctiluca scintillans*; *Dinophysis*; *Mesodinium rubrum*

#### **1. 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–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–16].

A large volume of domestic and industrial wastes from the city of Thessaloniki has been directed for decades to the Thermaikos Gulf and especially its inner part, Thessaloniki Bay. In the 20th century, these wastes were discharged in the Bay without treatment, causing the eutrophication of the system. Since 2001, wastewater treatment has been implemented, decreasing the effects of anthropogenic eutrophication [17]. Although it is generally accepted in the public that the water quality in the Thermaikos Gulf has been improved compared to 20–30 years ago [18,19], the urban front of the Thessaloniki Bay, with restricted water circulation and shallowness, still exhibits apparent red tides and algal blooms. These events are usually the cause of irritation to the public, often mentioned in the Greek media, with subsequent socio-economic consequences to the city of Thessaloniki, especially the touristic center. Despite the growing concerns of the citizens and authorities on the water quality of the Bay and particularly of the urban front, only scarce and isolated studies have been published on the abundance and dynamics of plankton community (both phyto- and protozooplankton) in the urban part of the Gulf [20,21]. On the other hand, several studies in the broader Thermaikos Gulf have focused on phytoplankton [18,22,23] and the occurrence of HABs [24,25]. However, comprehensive studies on red tides and mucilage aggregate phenomena are lacking for the Thermaikos Gulf.

According to the related legislation for the ecological water quality based on nutrient pressures and the phytoplankton quality element (Water Framework Directive; WFD, 2000/60/EC) there is a need for systematic and frequent monitoring of coastal waters. Furthermore, similar to the WFD objectives are those of Marine Strategy Framework Directive (MSFD, 2008/56/EC) for achieving good environmental status of EU marine waters by 2020. The MSFD eutrophication quality descriptor (D5) refers to the adverse effects of eutrophication including harmful algae blooms [26].

The aim of this paper was to examine the shift in the protagonists of the conspicuous and successive algal blooms, red tides and mucilage aggregations in the urban marine environment of the Thessaloniki Bay, by investigating the temporal and spatial changes of the unicellular eukaryotic plankton community attributes (species diversity, dominance, and abundance). This is the first study concerning the phyto- and protozoo-plankton species succession at an annual time scale with weekly samplings, with conspicuous phytoplankton blooms, red tides and a mucilage aggregate phenomenon, in the urban coastal front of the Thessaloniki Bay (Thermaikos Gulf). The present work focuses on the zone linking the terrestrial and the coastal environments heavily affected by, and influencing various human activities, such as the harbor, tourism, industry, mussel cultures, and sewage effluents.

#### **2. Materials and Methods**

#### *2.1. Sampling Sites and Sample Collection*

Samples were collected weekly from March 2017 to February 2018, from a coastal inshore sampling site in the White Tower (WT) in the center of the city of Thessaloniki (Table 1; Figure 1). During the same period, every month, additional samples were collected from three other inshore sites along the urban front of the Bay, namely at Aretsou Beach (AR), Music Hall coast (MH), and harbor (HB) (Table 1; Figure 1). In total, 47 samples were collected from WT and 12 from each other site (Table 1). All sampling sites had a maximal depth of 4 m.


**Table 1.** Sampling sites and total number of samples collected.

**Figure 1.** Study area in Thermaikos Bay, indicating the location of the four sampling sites (\*). WT: White Tower, AR: Aretsou, MH: Music Hall, HB: Harbor.

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.

#### *2.2. 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]. Nitrate and nitrite (NO3 <sup>−</sup> and NO2 <sup>−</sup>), ammonium (NH4 <sup>+</sup>), silicate (SiO4 <sup>−</sup>), and phosphate (PO4 −) were, also, measured according to Hansen and Koroleff [29].

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 <sup>−</sup> 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:


#### 5. Bad: >1.51

The E.I. was calculated according to the following equation:

$$\text{E.I.} = 0.279 \times \text{PO}\_4 + 0.261 \times \text{NO}\_3 + 0.296 \times \text{NO}\_2 + 0.275 \times \text{NH}\_3 + 0.214 \times \text{Chl.} \,\text{m}$$

#### *2.3. Microscopic Analysis*

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–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−<sup>1</sup> 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.

#### *2.4. 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).

#### **3. Results**

#### *3.1. 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), mean 37.2. The concentrations of inorganic nutrients (SiO4, PO4, NO3, NO2, NH4), and particulate organic phosphorus (POP) all exhibited strong fluctuations throughout the entire study period, with some extreme high values recorded in all sites, especially in WT (Supplementary Figure S1). In particular, remarkably high values were recorded for all nutrients on 22 March 2017 in the WT, (SiO4: 10.17 μmol L<sup>−</sup>1; PO4: 9.54 μmol L<sup>−</sup>1; NO2: 0.77 μmol L<sup>−</sup>1; NO3 and NO2: 7.53 μmol L<sup>−</sup>1; NH4: 160.3 μmol L<sup>−</sup>1; POP: 42.1 μmol L<sup>−</sup>1), followed by high values of most nutrients during the sampling on 28 June 2017 in WT and MH, 20 September 2017 in all sites, and 10 January 2018 in WT. Additionally, the highest value of NH4 (32.86 μmol L<sup>−</sup>1) was recorded on 18 October 2017 in HB. The annual means of Chl *a* for each station were as follows: WT: 2.62 μg L−1, AR: 3.11 μg L<sup>−</sup>1, MH: 1.43 μg L−<sup>1</sup> and HB: 3.15 μg L<sup>−</sup>1. Furthermore, Chl *a* showed marked variability, ranging from 0.27 μg L−<sup>1</sup> on 27 December 2017 in WT, to 17.28 μg L−<sup>1</sup> on 13 December 2017 in WT. The mean inorganic nutrient ratios were 25.1 for N:P, 0.70 for Si:N and 18.4 for Si:P. The maximum N:P (71.9) and Si:P (76.6) ratios in WT was measured on 28 June 2017 while the minimum N:P ratio (10.1) on 13 December 2017 and the minimum Si:P ratio (1) on 22 March 2017. In the rest of the sites, the maximum N:P and Si:P ratios were recorded in August samples and the minimum in November (HB) and December samples (AR, MH). The calculated Si:N ratios were relatively low (< 2) in all samples with maximum values coinciding with diatom blooms. All the calculated E.I. values exceeded the value 0.83 in all samples (data not shown), thus were indicative for eutrophication, reaching the highest value during the 22 March 2017 red tide event (49.5).

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 NO2 (with AR > MH), for NO3 and NO2 (WT > MH; AR > MH), and for POP (WT > MH; MH < HB) (Supplementary Table S1). Furthermore, higher ammonia (NH4) concentrations were recorded at HB (Table 2), even though no statistically significant differences between sites were found.


**Table 2.** Sample dates, sites and coding, and values of abiotic parameters (water temperature, salinity, SiO4, PO4, NO2, NO3 and NO2, NH4, POP – Particulate Organic Phosphorus, Chl *a*). All nutrient concentration values are given in μmol L<sup>−</sup>1. The sampling sites are shown in Figure 1.

**Table 2.** *Cont.*


**Table 2.** *Cont.*



**Table 2.** *Cont.*

#### *3.2. 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−<sup>1</sup> 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<sup>−</sup>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 (1-*D*) were: 22 March 2017 (0.10 in WT), 5 April 2017 (0.43 in WT), 31 May 2017 (0.22 in WT) and 18 October 2017 (0.29 in AR) (Supplementary Table S3).

**Figure 2.** Number of taxa (**a**), and total abundance (**b**) in White Tower (WT), Aretsou (AR), Music Hall (MH), and Harbor (HB). For sampling codes, see also Table 2.

#### *3.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) 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<sup>−</sup>1) of a taxonomic group or a single species, or/and red tides (Figure 3).

**Figure 3.** Cluster diagram according to Jaccard resemblance, calculated based on the non-transformed abundance (cells mL<sup>−</sup>1) of taxa during the study. Red clades in the dendrogram indicate sections of the plot where the observed profile corresponds to similarities that are larger than those expected under null conditions (>99% of the confidence envelope), suggesting the presence of true structure within the data. The nodes represent the dominant taxa blooming during the period covered by the corresponding clades.

During the period March–June 2017 a persistent diatom bloom was detected at the White Tower (WT) site, due to the high abundances recorded throughout for the taxa of *Leptocylindrus danicus* (max: >7000 cells mL−<sup>1</sup> on 24 May) and *Leptocylindrus minimus* (max abundance: >26000 cells mL−<sup>1</sup> on 31 May). In specific samplings, i.e., on 24 May 2017, the taxon *S. costatum* additionally showed high concentrations reaching > 1000 cells mL−1. The diatom bloom was accompanied by a Coccolithales bloom between 5 and 12 April 2017 (>5300 cells mL<sup>−</sup>1).

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<sup>−</sup>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−<sup>1</sup> on 19 July), and more rarely *C. closterium* (max: >1800 cells mL−<sup>1</sup> 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<sup>−</sup>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). In particular, on 8 November 2017, *D.* cf. *acuminata* reached 120 cells mL<sup>−</sup>1, in WT, while up to 350 cells mL−<sup>1</sup> of this species were recorded at HB on the 15th November 2017.

The high abundances of *D.* cf. *acuminata* in November 2017 were followed by a red tide in all sampling sites during the period December 2017–February 2018. This bloom was dominated by the photosynthetic ciliate *Mesodinium rubrum*, first appearing in the Bay on 29 November 2017, and peaking from 13 December 2017 till 10 January 2018, reaching > 1000 cells mL−<sup>1</sup> on 13 December. The period January–February 2018 was characterized by diatom dominance, i.e., the taxa *Chaetoceros tenuissimus* (max abundance > 2000 cells mL−<sup>1</sup> on 6 February, WT), and *S. costatum* (reaching 1000 cells mL−<sup>1</sup> on 3 January, WT).

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−<sup>1</sup> at the White Tower (WT) site on 19 July 2017, just after the mucilage aggregate phenomenon. Relatively high abundances (270 cells mL<sup>−</sup>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<sup>−</sup>1) was observed at the same sampling point one year later (unpublished data).

#### *3.4. 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 NH4 and PO4, 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).

**Figure 4.** Network diagrams of highly significant connections (*p*-values < 0.01) based on the maximal information coefficient (MIC) scores between dominant taxa (comprising of > 10% of the total plankton abundance in at least one sample) and environmental parameters. Boxes (nodes) with indicated taxa names, represent bloom (detected with abundances > 1000 cells mL−<sup>1</sup> in at least one sample throughout the samplings), and red tide forming taxa. To facilitate reading only bloom and red tide forming taxa are indicated. Black lines (edges) depict positive connections, and red edges depict negative connections.

Dinoflagellates were found to be significantly positively correlated with E.I. (*p* < 0.001, *r* = 0.64), N:P (*p* < 0.01, *r* = 0.31) and Si:P (*p* < 0.05, *r* = 0.26), while Cryptophyceae were significantly negatively correlated with Si:N (*p* < 0.01, *r* = −0.31) and Si:P (*p* < 0.05, *r* = −0.22).

#### **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].

#### *4.1. 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–44]. In a recent comprehensive study on various coastal areas of Greece (1995–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 PO4 (max 6.50 μmol L−1) and NH4 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−<sup>1</sup> and 160 μmol L−<sup>1</sup> for PO4 and NH4, respectively) were even higher [45]. On the other hand, the highest NO3 value (21.09 μmol L<sup>−</sup>1) in the urban front of Thessaloniki Bay during our survey was slightly lower than the highest measured NO3 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 NH4 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 NH4 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.

#### *4.2. 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].

#### *4.3. 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].

#### *4.4. 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–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<sup>−</sup>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−<sup>1</sup> [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 < 3ms−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 NH4 and PO4 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]. NH4 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]. NH4 also increases during decline bloom phase indicating

release of intracellular NH4 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 NH3 [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].

#### *4.5. 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 chemical biochemical 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].

#### **5. 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<sup>−</sup>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,

**Author Contributions:** Conceptualization: S.G. and M.M.-G.; methodology: S.G., N.S. and M.M.-G.; validation: S.G., N.S. and M.M.-G.; formal analysis: S.G., and M.M.-G.; investigation: S.G. and N.S.; resources: U.S. and M.M.-G.; data curation: S.G. and N.S.; writing–original draft preparation: S.G., N.S., and M.M.-G.; writing–review and editing: S.G., U.S., and M.M.-G.; visualization: S.G.; supervision: M.M.-G.; funding acquisition: S.G., U.S., and M.M.-G.

**Funding:** This research is implemented through 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.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Photosynthetic Picoeukaryotes Diversity in the Underlying Ice Waters of the White Sea, Russia**

#### **Tatiana A. Belevich 1,2,\*, Ludmila V. Ilyash 1, Irina A. Milyutina 2, Maria D. Logacheva <sup>2</sup> and Aleksey V. Troitsky 2,\***


Received: 5 February 2020; Accepted: 3 March 2020; Published: 5 March 2020

**Abstract:** The White Sea is a unique basin combining features of temperate and arctic seas. The current state of its biocenoses can serve as a reference point in assessing the expected desalination of the ocean as a result of climate change. A metagenomic study of under-ice ice photosynthetic picoeukaryotes (PPEs) was undertaken by Illumina high-throughput sequencing of the 18S rDNA V4 region from probes collected in March 2013 and 2014. The PPE biomass in samples was 0.03–0.17 <sup>μ</sup>g C·L−<sup>1</sup> and their abundance varied from 10 cells·mL−<sup>1</sup> to 140 cells·mL<sup>−</sup>1. There were representatives of 16 algae genera from seven classes and three supergroups, but Chlorophyta, especially Mamiellophyceae, dominated. The most represented genera were *Micromonas* and *Mantoniella*. For the first time, the predominance of Mantoniella (in four samples) and Bolidophyceae (in one sample) was observed in under-ice water. It can be assumed that a change in environmental conditions will lead to a considerable change in the structure of arctic PPE communities.

**Keywords:** White Sea; under-ice water; picoeukaryotes; *Micromonas*; *Mantoniella*; high-throughput sequencing; metagenomics; 18S rDNA

#### **1. Introduction**

The picophytoplankton (cyanobacteria and photosynthetic eukaryotes with cell diameter <3 μm) make up the smallest component of phytoplankton populations [1–3]. In the Arctic region, photosynthetic picoeukaryotes (PPEs) are major contributors to picophytoplankton and the small phytoplankton (<5 μm) represent 59%–63% of all marine photosynthetic biomass [4]. The Arctic has been undergoing accelerated warming and freshening since the 1990s due to the melting of multiyear sea ice and increasing river runoff into the Arctic Basin [5,6]. Environmental changes affect the phytoplankton and lead to an increase in PPEs contribution to total primary production and phytoplankton biomass [7–9].

Correct taxonomical identification of PPE requires the use of molecular methods. 18S rRNA gene-based environmental surveys have been increasingly used to investigate the composition of small eukaryotes. Molecular environmental studies conducted in Arctic and subarctic waters reported the presence of diverse microbial communities [10–17]. To better assess the diversity of small photosynthetic eukaryotes, flow cytometry via chlorophyll fluorescence was used to sort cells successfully [16,18–20]. However, the use of flow cytometry does not prevent the detection of heterotrophic eukaryotes that ingest photosynthetic organisms in their food vacuoles and thus could be detected by flow cytometry sorting that targeted chlorophyll fluorescence as well as the detection of photosynthetic algae with cell sizes larger than those of picoforms [20,21]. PPEs in under-ice waters remain understudied, especially for the season preceding the under-ice bloom in spring [10,11,17,22].

The White Sea is a small (area of 90,000 km<sup>2</sup> and volume of 6000 km3) subarctic semi enclosed basin with an outlet to the Barents Sea. It has features similar to those of the Arctic shelf seas [23]. Usually from December to May the sea is covered with ice. The White Sea is strongly affected by continental runoff, and its waters are less saline (14–27 psu) than open ocean waters. The species composition and abundance of plankton algae have been studied for almost 80 years in the White Sea [24]. The species richness of nano- and microphytoplankton of the White Sea has been studied by microscopy and is represented by 450 taxa [24]. However, most of the studies examining the taxonomic diversity of the White Sea algae have been limited to the easily recognizable nano- and microsized algae, while the PPE composition remained understudied.

Previously, the taxonomic composition of PPEs was studied in summer plankton [25,26] and in the sea ice of the White Sea [27]. PPEs of the White Sea ice were represented by 16 algae genera belonging to eight classes and three supergroups. Chlorophyta, especially Mamiellophyceae, dominated among ice PPEs. The composition of the underlying ice waters' eukaryotic picophytoplankton in the White Sea was estimated for only one sample [28].

Considering the ongoing changes in the Arctic Ocean caused by global warming, and their implications, it is crucial to understand the PPE composition and provide detailed data on the prevalent taxa in subarctic waters. Hence, the objectives of this study were: i) to evaluate under-ice picophytoplankton abundance and biomass and the contribution of PPEs to total picophototroph abundance, ii) to reveal the taxonomic diversity of PPEs during the early spring by way of 18S rDNA sequencing, and iii) to compare the PPEs composition in under-ice water and ice. We targeted the study of the smallest size class of algae by a sample filtration approach, as they are abundant in the Arctic Ocean and difficult to identify by microscopy.

#### **2. Materials and Methods**

#### *2.1. Sampling and Study Area*

The samples were collected in Kandalaksha Bay, the White Sea on 19–23 March 2013 and 16–19 March 2014 near the White Sea Biological Station, Lomonosov Moscow State University (66◦33' N, 33◦06' E), from five different stations (Figure 1) with various degrees of under-ice water salinity and under-ice water current speeds (Table 1). The under-ice water salinity was lower due to the impact of the freshwater runoff at stations 1, 2, and 3 with the maximum freshening level at station 1. Under-ice water at stations 4 and 5 had salinity characteristics of the White Sea surface layer in the winter. The waters at stations 2, 4, and 5 are characterized by higher speeds of under-ice water currents than water areas of stations 1 and 3 [29]. Water samples at stations 1 and 2 were taken twice—in 2013 and 2014. The reference number of each sample consists of the station number, the last two digits of the year when the sample was taken, and the letter "w," which means water samples (e.g., 1/13w).

*Diversity* **2020**, *12*, 93

**Figure 1.** Location of the sampling stations in Kandalaksha Bay, White Sea in March 2013 and 2014.


*Diversity* **2020** , *12*, 93

At each station, a titanium manual ice corer (14 cm of diameter) was used to make holes in the ice to collect 5 L of the underlying water. The ice thickness and water temperature were measured. Water temperature was measured directly with a probe Testo 108 (Testo, Lenzkirch, Germany). Within 1 h, water samples were brought to the laboratory, where the salinity was measured with a conductivity probe Cond 3150i (Xylem Analytics Germany Sales GmbH & Co. KG, WTW, Weilheim, Germany). The air temperature recordings from the weather station at the White Sea Biological Station were also used.

#### *2.2. Chlorophylla*

Subsamples of underlying water (500–1000 mL) were filtered through Whatman GF/F filters and frozen (−80 ◦C) for subsequent analysis. On returning to the Moscow laboratory, extractions and calculations were made following the procedure [30].

#### *2.3. Enumeration of Picophototrophs*

Whole (not prefiltered) seawater samples (10 mL) intended for analysis by epifluorescence microscopy were fixed with glutaraldehyde (AppliChem Panreac, Barcelona, Spain) at a final concentration of 1% (*v*/*v*). Nuclear filters (0.12 μm pore diameter) prestained with Sudan black were used for filtration. Cells with sizes <3 μm were enumerated at × 1000 magnification with a Leica DM2500 (Leica Microsystems GmbH, Wetzlar, Germany) epifluorescence microscope equipped with a 50 W mercury lamp under blue (Filter D; 355–425 nm) and green (Filter N2.1; 515–560 nm) excitation. The bright yellow fluorescence of the phycoerythrin-containing cyanobacteria could be distinguished easily from the deep red fluorescence of the chlorophyll-dominant picoeukaryotes. At least 300 cells at 30–50 microscopic fields were counted for each sample. Cell volumes were calculated as volumes of the relevant geometrical bodies [31] and then converted to their carbon content using the conversion factors of 470 ƒg C cell−<sup>1</sup> for prokaryotes and 0.433 <sup>×</sup> (V)0.863pg C cell−<sup>1</sup> for PPEs [32].

#### *2.4. DNA Isolation of Picoplanktonic Size-fraction*

Three to five liters of water samples were filtered through a 2-μm pore size polycarbonate filter and then filtered again through 0.2-μm Sterivex units (Millipore Canada Ltd, Mississauga, ON, Canada). The buffer was added to the Sterivex units (1.8 mL of 50 mM Tris-HCl, 0.75 M sucrose, and 40 mM EDTA; pH 8.3). These units were stored at −80 ◦C until DNA extraction using the NucleoSpin Plant kit (Macherey-Nagel, Düren, Germany) following the manufacturer's instructions.

#### *2.5. DNA Amplification and Sequencing*

The ~0.43-kb fragments of the hypervariable V4 region of the 18S rRNA were amplified with the primer pair EuF-V4 (5'-CCAGCASCCGCGGTAATWCC-3') and pico-R2 (5'-AKCCCCYAACTTTCGTTCTTGAT-3') [27]. For PCR, the Encyclo Plus PCR kit (Evrogen, Moscow, Russia) was used. The volume of the amplification mixture was 30 μL. This was divided into three equal parts (10 μL each), and then PCR was carried out for each sample at three annealing temperatures, 55 ◦C, 60 ◦C, and 65 ◦C [27]. Cycling conditions were as follows: an initial denaturation step for 3 min at 94 ◦C, followed by 30 cycles (denaturation at 94 ◦C for 20 s, annealing for 20 s, and extension at 72 ◦C for 40 s), followed by a final extension at 72 ◦C for 5 min.

The PCR products obtained at three annealing temperatures were combined and, after extraction by agarose gel electrophoresis, were purified with a Cleanup Mini kit (Evrogen, Moscow, Russia). The resulting amplicons were used to prepare the libraries for the sequencing on the Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA) with a TruSeq Nano DNA Kit (Illumina Inc., San Diego, CA, USA). The maximum read length of the Illumina MiSeq technology is about 500–600 bp, which matches and even exceeds the length of V4 of the SSU rRNA. The hypervariable V4 region of the 18S rRNA revealed an impressive hidden diversity in picoplanktonic communities [3,33].

The effective concentration of the libraries was tested by quantitative PCR with the primers I-qPCR-1.1 (5'-AATGATACGGCGACCACCGAGAT-3') and I-qPCR-2.1 (5'-CAAGCAGAAGACGGCATACGA-3'). The library PhiX Control v3 (Illumina) was used as a control. Then the libraries were diluted to 12 pM and sequenced with a MiSeq Reagent Kit v.2 for 500 cycles. The pair-end read length was 250 × 2 bp.

#### *2.6. Bioinformatics and Data Evaluation*

The raw sequencing data were processed by Mothur software [34] and other procedures implemented in the SOP protocol [35]. Reads shorter than 150 bp and longer than 550 bp were removed, as well as reads with ambiguous bases (Ns) or >6 repeated bases. Assembled contigs were 430 bp in length, with ~70 bp overlapped paired reads. Identical sequences were removed by the unique.seqs command. The sequences were aligned using MAFFT with FFT-NS-2 strategy, a gap-opening penalty of 1.53, and a gap extension penalty of 0.123. Putative chimeric sequences were identified by UCHIME v 4.2.40 [36] and removed. A distance matrix of the high-quality sequences was constructed, and the sequences were clustered into operational taxonomic units (OTUs) at a 97% similarity level, with average neighbor clustering using cd-hit-est v,.3.1.2 [37]. The classification was performed by a local nucleotide BLAST search against the nonredundant version of the SILVA 123 SSU RNA database [38] using blastn (version 2.2.28+) with standard settings [39]. Sequences affiliated with nonprotist phyla or bacteria were eliminated. All singletons were removed. Consensus sequences of the OTUs were generated by the script described earlier [27]. All sequence reads were submitted to the GenBank BioProject (PRJNA368621) under the accession numbers MK571487-MK571523, MN541095, and MN684208. The phylogenetic tree was inferred by maximum likelihood method using RAxML 8.2.10 program [40], with default options according to GTRGAMMA model with 400 bootstrap replications, the number of which was set by bootstrapping criterion implemented in RAxML. The secondary structures of the terminal hairpins of V4 rRNA region were constructed according small subunit RNA secondary structure model [41].

Since the purpose of our research was limited to the photosynthetic picoeukaryotes, from the complete list of taxa revealed in under-ice water samples filtered through a 2-μm pore size filter, we chose only those species that have a cell size ≤3 μm. Where OTUs were identified to the genus/order/class level, we only analyzed taxa that, according to the published data [1,42,43], have species corresponding to the pico-size fraction. Since photosynthetic pico-sized cryptophytes were not detected microscopically and photosynthetic pico-sized dinoflagellates are not currently described, these groups were excluded from the analysis. Classes of algae are given according to AlgaeBase [44].

#### *2.7. Statistical Procedures*

The similarity matrix was calculated after standardization of the abundance of PPEs reads and square-root transformation for reducing the influence of the most dominant taxonomic entries [45]. The PRIMER v6 software (Primer-E Ltd, Plymouth, UK) [46] was used to group samples with similar taxonomic compositions by a group-average linkage cluster analysis and a nonmetric multidimensional scaling (MDS) ordination of a Bray–Curtis similarity matrix [45]. A breakdown of species similarities (SIMPER) was used to determine which taxon combination leads to the resulting groups [47].

#### **3. Results**

#### *3.1. Environmental Conditions*

The year 2014 was warmer than 2013 during the sampling period and the preceding month (Figure S1). The air temperature occasionally rose to the water freezing point or above, even in the middle of winter. Kandalaksha Bay was partially ice covered from December, with more complete ice cover in late March in both 2013 and 2014. Ice thickness varied from 22 cm to 71 cm (Table 1). Under-ice water salinity was lowest at station 1 in both years and varied between 14.9 psu and 15.6 psu (Table 1). At the other stations, the salinity of the under-ice water ranged from 21.9 psu to 26.7 psu. The temperature of the under-ice water varied very little: between −0.7 (2/13w) and −1.2 ◦C (2/14, 5/14w) with an average of −1.0 ◦C.

#### *3.2. Total Chlorophyll a Biomass*

Total chlorophyll *a* biomass (Chl *a*) level varied from 0.05 <sup>μ</sup>g·L−<sup>1</sup> at sample 2/13w to the highest value of 0.73 <sup>μ</sup>g·L−<sup>1</sup> at sample 2/14w (Table 1). The average Chl *<sup>a</sup>* concentration in under-ice water was 0.27 <sup>±</sup> 0.21 <sup>μ</sup>g·L<sup>−</sup>1.

#### *3.3. The Abundance of Picophototrophs*

The PPEs abundance ranged from 10 cells·mL−<sup>1</sup> to 140 cells·mL−<sup>1</sup> with an average of 50 cells·mL<sup>−</sup>1. The biomass varied between 0.03 and 0.17 <sup>μ</sup>g C·L−<sup>1</sup> (Table 1). Among photosynthetic pico-sized organisms, cyanobacteria dominated in all samples except 3/14w where we did not reveal photosynthetic prokaryotes. The relative abundance of PPEs varied significantly between 5% and 100% of the total cell counts and carbon biomass of pico-sized photosynthetic organisms.

#### *3.4. Taxonomic Composition of Eukaryotes in Samples Filtered through a 2-*μ*m Pore Size Filter*

A total of 268,124 amplicons were sequenced from the seven samples, and 122,503 reads remained after quality filtering and preprocessing. The relative abundance of PPEs reads was 11%. The number of OTUs (at the 97% similarity level) that were clustered in individual samples varied between 609 and 3856 (Table 2).


**Table 2.** Summary of recovered reads and the number of operational taxonomic units (OTUs) in under-ice water picoplankton samples.

Different OTUs are grouped according to their taxonomic affiliations to major phylogenetic groups, such as Chloroplastida, Stramenopiles (Bacillariophyta, Bolidophyceae, Chrysophyceae, Dictyochophyceae, Raphidophyceae, Pelagophyceae, Eustigmatophyceae), Alveolata, Rhizaria, Cryptophyta, Haptophyta, Opisthokonta, and others (Centrohelida, Telonemia, Kathablepharidae, Picozoa, and Eukaryota incertae sedis) (Figure 2). Protists from the taxonomic groups Rhizaria, Opisthokonta, Centrohelida, Telonemia, Kathablepharidae, and Picozoa are nonphotosynthetic forms. Alveolata and Cryptophyta include heterotrophic species. A total of 175 taxa of protists, determined to the genus level, and 148 forms, determined to higher taxonomic ranks, were found in water samples filtered through a 2-μm pore size filter.

**Figure 2.** The relative abundance (%) of V4 rDNA reads of the major protist groups in picoplankton samples.

#### *3.5. OTU Richness and Taxonomic A*ffi*liation of the PPEs Sequences*

PPEs belong to three supergroups: Chloroplastida (Chlorophyta), Stramenopiles, and Haptophyta (Table 3). Since different samples yielded different total numbers of sequence reads, they were normalized based on the lowest sample size (sample 1/13w—7398 reads) for comparing OTUs richness. The expected OTUs richness of PPEs was calculated with a 95% probability (Table 3). The minimum expected OTUs richness of PPEs species was observed at the lowest water salinity, in sample 1/13w. The highest expected richness was found at the highest values of salinity, in samples 4/14w and 5/14w.

**Table 3.** Relative abundance (%) of PPE groups based on V4 rDNA reads and the expected OTUs richness of PPEs per sample (95% probability). The standardized number of OTUs for each group is indicated in parentheses.


#### *3.6. Chlorophyta*

Chlorophyta were represented by four classes: Mamiellophyceae, Pyramimonadophyceae, Palmophyllophyceae, and Trebouxiophyceae (Table 3). Mamiellophyceae were the predominant PPEs

in all samples of under-ice water except sample 2/13w, where Bolidophyceae dominated. Within the Mamiellophyceae, *Micromonas* was the dominant genus in three samples and the *Mantoniella* genus dominated in four samples (Table 4). The negative correlation between the relative abundance of *Micromonas* and *Mantoniella* reads was found (Rs = −1; *p* = 0.003). The genera *Ostreococcus*, *Bathycoccus*, *Crustomastix* and OTU similar to the uncultured clone DSGM81 were also detected (Figure 3, Table 5).

**Table 4.** A relative abundance (%) of reads found for different taxonomical groups of Mamiellophyceae based on the V4 region of the 18S rRNA gene sequences.


**Figure 3.** Phylogenetic tree of revealed Mamiellophyceae OTUs and Genbank reference sequences constructed from the V4 region of the 18S rRNA gene sequences by the maximum likelihood method. Bootstrap supporting values >0.5 are indicated. The scale is a number of nucleotide substitutions per site. Clades are designated according to Tragin and Vaulot (2019) [43] and Belevich et al. (2018) [27].


**Table 5.** The most abundant PPEs OTUs recovered in the under-ice water of Kandalaksha Bay, the White Sea in March 2013 and 2014 (clustering at 97% similarity threshold). The number of reads of each OTU is indicated in parentheses.


**Table 5.** *Cont.*

The *Micromonas* genus was represented by two species and a recently described clade, corresponding to clade F [27] or clade B3 [43]. *Micromonas polaris* (previously *M. pusilla* clade Ea) was revealed in all samples; its contribution to the total *Micromonas* reads varied between 90% (2/13w) and 100% (1/13w, 1/14w, 2/14w, 3/14w). *Micromonas commoda* clade A2 [43] and *Micromonas* clade F (B3) were revealed much less often—in only three samples (2/13w, 4/14w and 5/14w) and two samples (2/13w and 5/14w), respectively. The contribution of each to the total *Micromonas* reads was low (Table 4).

There were several *Mantoniella* phylotypes from the three clades Ms, Mb, and A in the under-ice water of the White Sea (Figure 3). Three phylotypes were found in all samples: the first was identical to *Mantoniella squamata* (X73999), the second was matched to *Mantoniella beaufortii* (KT860921), and the third from clade A [43] was similar (>98%) to the Uncultured Chlorophyta clone 5-D5 (FN690723) from the Baltic Sea. In general, three phylotypes assigned to clade A were discovered in our samples. Two of them (MK571493 and MK571492) were previously found in the ice of the White Sea and identical to environmental sequences *Mantoniella* MF589930 and *Mantoniella* MF589928, respectively. The third *Mantoniella* phylotype (MK571498), with 100% similarity to the Uncultured Prasinophyceae clone North Pole SI120\_29 (HQ438123) from the marine ice, was revealed in sample 2/13w. Substitutions in basal parts of helixes E23\_1 and H 25 of 18S rRNA are diagnostic for distinguishing the phylotypes *M. squamata* and *M. beaufortii* and three other phylotypes of clade A, MK571493, MK571492, and MK571498 (Figure 4).


**Figure 4.** Compensatory base changes in the helices of the 18S rRNA secondary structure of *Mantoniella* (helices E23\_1 and H25). The CBCs are shown in rectangles.

*Bathycoccus* OTUs found in samples 2/13w, 2/14w, 4/14w, and 5/14w were identical (100%) to the *Bathycoccus prasinos* strain RCC801 (KT860937). The relative read abundance of *B. prasinos* did not exceed 5% of Mamiellophyceae reads. The OTU matching *Ostreococcus tauri* (Y15814) was revealed only in sample 5/14w with a low (<1%) relative read abundance. Moreover, in all samples, we revealed OTUs that showed 99.2% similarity to the uncultured eukaryotic clone DSGM-81 (AB275081). The previous molecular phylogenetic analysis revealed that clone DSGM-81 belongs to Mamiellophyceae [27]. *Crustomastix* OTUs were found in four samples and showed 100% similarity to environmental sequences of uncultured *Crustomastix* (MF589934) previously identified in the White Sea ice.

Pyramimonadophyceae sequences were represented by Pyramimonas OTU, which is similar to Pyramimonas sp. (JF794047) from the Beaufort Sea. Among Palmophyllophyceae OTU, Prasinoderma similar to the *Prasinoderma coloniale* strain RCC854 was identified. Trebouxiophyceae from samples 1/13w, 1/14w, and 2/13w were identical to the freshwater algae *Choricystis minor* (X89012). Additionally, in all samples Trebouxiophyceae was represented by OTU identical to *Picochlorum* sp. RCC748 (KT860896) from the Atlantic Ocean.

#### *3.7. Stramenopiles*

Stramenopiles were represented by two classes: Bolidophyceae and Mediophyceae. All diatoms reads were revealed only in sample 5/14w with the highest sequencing depth. Among Mediophyceae, OTU identical to sequences of two different species, *Minutocellus polymorphus* (LC189088) and *Arcocellulus cornucervis* (JN934677), were revealed. *Skeletonema marinoi* (KR091067) and *Chaetoceros* cf. *neogracilis* (JN934684, KT860998) were also identified.

Bolidophyceae were revealed in all samples (Table 3, Figure 5). They were represented by sequences of *Triparma strigata* (KR998402) with 100% similarity, and OTUs similar to three phylotypes of uncultured bolidophytes earlier revealed in the ice and summer plankton of the White Sea [48], two phylotypes from the plankton of the Pacific Ocean (LC191051, LC190998), two uncultured stramenopiles from the Baltic Sea ice (FN690655, FN690656), and five uncultured eukaryotes from the Greenland Sea (KT818381, KT811782, KT814386, KT815972, and KT813573). Bolidophytes were the predominant PPEs in sample 2/13w.

37

**Figure 5.** The maximum likelihood Bolidophyceae phylogenetic tree constructed from the V4 region of the 18S rDNA. Bootstrap support values >50% are indicated. The scale is the number of nucleotide substitutions per site.

#### *3.8. Haptophyta*

The contribution of Haptophyta reads varied between 3% and 10% (Table 3). Among the Haptophyta, two genera of class Coccolithophyceae were found—*Phaeocystis* and *Chrysochromulina*. *Phaeocystis* OTUs were similar (>99%) to the *Phaeocystis pouchetii* isolate AJ01 (KR091066) from the North Sea. Sequences closely related to the uncultured *Chrysochromulina* clone MALINA (JF698782) from the Beaufort Sea, occurred at insignificant levels only in the 5/14w sample. *Chrysochromulina simplex* (AM491021) was only found in sample 2/13w. Moreover, two Haptophyta phylotypes were classified at the phylum level with >98% similarity: the uncultured haptophyte (FN690514) from the Baltic Sea and the uncultured haptophyte (KC488456) from the North Atlantic Ocean.

#### *3.9. PPEs Community Structure*

Contributions of different taxa OTU reads to total numbers of PPE reads resulted in grouping of the stations into two clusters: CI at 64% similarity (1/13w, 1/14w and 3/14w) and CII at 62% similarity (4/14w, 2/13w, 2/14w and 5/14w) (Figure 6). SIMPER analysis revealed that cluster C1 was characterized by a high contribution of *Mantoniella* reads to the total number of PPEs reads (53%), and cluster C2 was formed by stations with a high contribution of *Micromonas* (22%).

**Figure 6.** Community comparison of PPE assemblages at the different sampling stations using nonmetric multidimensional scaling (MDS) of a data matrix based on Bray-Curtis similarity.

#### **4. Discussion**

Our study revealed the most complete genetic diversity of under-ice PPEs in the White Sea, a unique marine environment combining features of temperate and Arctic seas. Such uniqueness of the abiotic environment was reflected in the composition of pico-sized photosynthetic organisms: besides widespread taxa (*Bathycoccus prasinos*, *Ostreococcus tauri*), the Arctic endemic *Micromonas polaris* (previously *M. pusilla* clade E2) and temperate waters *Micromonas commoda* A2 and *Mantoniella* were revealed. Temperate-water taxa survive in the White Sea despite extreme environmental conditions under the ice, i.e., near-freezing temperatures, polar night, and low irradiance, because of the snow and ice cover.

In the under-ice water of Kandalaksha Bay, Chl *a* concentrations in March 2013–2014 were double the values recorded in water samples taken directly underneath the ice of Chupa Inlet of Kandalaksha Bay in February 2002, but half the Chl *a* values in the same inlet in April 2002 [49]. This indicates that our studies of under-ice plankton algae were carried out in the prebloom period. Studies of the biomass plankton algae dynamics from January to April in Kandalaksha Bay also revealed the highest values in April [50]. The values of Chl *a* concentration obtained in the under-ice water of the White Sea are comparable to those in the under-ice water of the Baltic Sea in March (0.5–1.0 μg L<sup>−</sup>1) [51].

The abundance of photosynthetic picoeukaryotes in our study was lower than the total number of PPEs observed in ice-covered underlying waters of Kandalaksha Bay (near station 2 in this study) in April 2010 [52]. At the same time, the average PPE abundance was comparable to the numbers of small photosynthetic eukaryotes (<2 μm) found in the under-ice water of Franklin Bay in December-March, but was approximately one order of magnitude lower than the abundance observed there in April-May [53]. In the under-ice water, PPEs were most abundant in three samples, while cyanobacteria dominated in the biomass in the other four samples. Earlier, the dominance of cyanobacteria was revealed in the summer plankton of the Onega [25] and Kandalaksha [52] bays of the White Sea.

Eukaryotic picoplankton is phylogenetically very diverse and includes lineages not yet described [1,54,55]. High-throughput sequencing of 18S rDNA of pico-sized eukaryotes from under-ice water yielded a detailed view of the plankton PPE community in the ice-covered underlying waters of the subarctic sea. As is commonly found in picoplankton diversity studies based on filtration approaches, sequences from larger protists and metazoans were recovered, probably due to cell breakage and the deformation of flexible walled cells, allowing their DNA and RNA to pass through the 3-μm filters [14,17,27,56]. However, some factors should be noted that can potentially lead to distortions in the estimation of real picoplankton diversity when using filtration. Among them are a probable breakage and deformation of larger cells [14,17,27,56], the presence of extracellular DNA in filtrates [57–59], and inaccurate size-based fraction separation. The metagenomic approach also has its limitations because of possible overestimation of particular OTUs due to high rRNA gene copy numbers or artifacts of the sequencing procedure [60,61], and insufficient 18S rDNA V4 region sequence resolution for detection of all morphospecies, as has been shown for diatoms, for example [62,63]. The exclusion of dinoflagellate sequences from the analysis due to the lack of cultured representatives of this group with cell sizes of less than 3 μm could affect the accuracy of our PPE diversity estimates [1]. Considering the above, 16 algae genera from seven classes and three supergroups are detected among the White Sea under-ice PPEs. Mamiellophyceae dominated in most of the samples. Palmophyllophyceae and Mediophyceae were the minor component of PPEs, Bolidophyceae made a significant contribution and even predominated in one sample.

Most sequences were assigned to Chlorophyta OTUs. Chlorophyta reads were abundant in mid-April during the early phase of the spring bloom in Norwegian waters (Isfjorden, West Spitsbergen) [64]. The dominance of Chlorophyta sequences was repeatedly noted in summer plankton communities of temperate and arctic waters [22,65–67]. On the class level, most Chlorophyta sequences were assigned to Mamiellophyceae OTUs, among which *Micromonas* and *Mantoniella* reads dominated. *Micromonas* was represented by species *Micromonas polaris*, *M. commoda* A2, and a phylotype of a recently described clade F according to Belevich et al. [27], or B3 according to Tragin and Vaulot [43]. *Micromonas polaris* dominated among *Micromonas*. *M. polaris* is widespread and dominant in the under-ice and open Arctic waters [11,17,43,68,69] and regarded as arctic endemic [11,19]. Previously, we detected *M. polaris* (as *Micromonas* E2) in the ice of the White Sea [24]. Its presence in the under-ice water and summer plankton of the subarctic White Sea and the Gulf of Finland of the subarctic Baltic Sea [43] once again indicates that the area of distribution of this species is wider than previously thought [11]. This endemic *M. polaris* does not seem to show intraspecific variability [70].

*Micromonas commoda* was detected earlier in the White Sea ice and summer plankton as *Micromonas* clade C [25,27]. This study is the first to show *M. commoda* in under-ice water. Its relative contribution to total *Micromonas* reads varied from 0% to 6% and was lower than in summer plankton [25]. *M. commoda* has ubiquitous distribution [43,71,72] and constitutes <1%–40% of Mamiellophyceae reads in different regions [43]. Unlike *M. polaris*, which does not show any intraspecific variability, the genetic diversity within *M. commoda* was previously highlighted [43,73,74] and it was suggested that speciation events might be ongoing within this species [72].

This study is the first to discover *Micromonas* clade F (B3) in the under-ice water of the White Sea. Previously, the phylogenetic analysis of Mamiellophyceae revealed the existence of a new clade, *Micromonas* F, in the ice and summer plankton of the White Sea [25,27]. Later, the analysis of the taxonomic diversity and global distribution of *Micromonas* revealed the existence of *Micromonas* clade B3 [43]. This subarctic clade combined amplicons that are 100% similar to OTUs of *Micromonas* clade F from the White Sea ice. *Micromonas* clade B3 made a great contribution to Mamiellophyceae reads from Canada (32%) and amounted to more than 10% of Mamiellophyceae reads at four subarctic stations off Maine and Iceland, as well as at a temperate location off the U.K. coast in the North Sea [43]. The contribution of this taxon to the total Mamiellophyceae reads in under-ice plankton was <1%, which is much lower than in summer plankton [25].

The results of this work represent the first sighting of a broad diversity of *Mantoniella* phylotypes in the under-ice water (Figure 3): *Mantoniella* OTUs match *M. squamata* (clade Ms), *M. beaufortii* (clade Mb), and three other *Mantoniella* lineages from clade A [43] found earlier in the White Sea ice [27]. It has been suggested that clade A is potentially an ice alga [43]. This assumption does not agree with the fact that *Mantoniella* phylotype MK571493 dominated among the Mamiellophyceae reads in two out of seven samples of under-ice water and did not dominate in any of the ice samples [27]. The *Mantoniella* MK571493 and MK571492 are 98.4% similar; however, out of six substitutions making those two phylotypes different from each other, four are compensatory. Likewise, the *Mantoniella* MK571498 sequence is 97.3%–97.8% similar to two other sequences from the same A clade. Substitutions in basal parts of helixes E23\_1 and H 25 are differentiated in all known *Mantoniella* phylotypes (Figure 4). Further research may lead to the discovery of new species in the *Mantoniella* A clade.

OTUs of *Ostreococcus* found in the White Sea plankton were identical to the sequence of *O. tauri* Y15814, which was isolated from the Thau lagoon (Mediterranean Sea) with highly variable salinity from 24 to 38 psu [75]. It was suggested that *O. tauri* might represent several species, adapted to different degrees of salinity [43]. The relative abundance of *O. tauri* in under-ice water (<0.1%) was significantly lower than that in summer plankton (31%) [25].

The relative abundance of *Bathycoccus prasinos* reads in under-ice water of the White Sea did not exceed 5%, whereas in summer plankton it reached 33% [25]. *B. prasinos* is a widespread alga with global distribution from tropical to polar waters [43,76]. *Bathycoccus* is now known to be composed of two cryptic species with identical 18S rRNA sequences but differences in the ITS, as well as at the genomic level [76,77]. One of them could be coastal, while another might have adapted to warmer oceanic waters [76–78].

In under-ice water, Palmophyllophyceae and Pyramimonadophyceae were represented only by one taxon each; their respective contributions to total PPEs reads were low. The cell size of *Pyramimonas* sp. (JF794047) is unknown and its assignment to picoforms may be inappropriate. In our samples, several *Pyramimonas* taxa with nanosizes were found (*Pyramimonas mucifera*, *P. olivacea*, *P. tetrarhynchus*, etc.), which is consistent with the high diversity of nanoforms of this genus in arctic and subarctic plankton [19,79].

Trebouxiophyceae were represented by both freshwater *Choricystis minor* and marine *Picochlorum* sp. algae. Its contribution to the total PPEs reads was low (Table 3) and exceed 5% only at station 1 (samples 1/13w and 1/14w), which was most affected by river flows (Table 1). Earlier, the dominance of Trebouxiophyceae reads was revealed in the ice at this station [27].

Diatoms made an insignificant contribution to the PPE community and were found only in sample 5/14 with the maximum sequencing depth. Among all identified diatom taxa, we can confidently assert that only *Minutocellus polymorphus*/*Arcocellulus cornucervis* and *Skeletonema marinoi* match the picofraction. Unfortunately, the resolving power of V4 is insufficient to correctly identify *M. polymorphus* and *A. cornucervis*—their V4 regions are identical [63]. Earlier, *M. polymorphus* was not recorded in the under-ice plankton communities; like *S. marinoi*, it was registered in the sympagic communities of the White Sea [27]. The revealed OTU of *Chaetoceros* cf. *neogracilis* matches two culture representatives of this species deposited in the RCC culture collection (Roscoff, France) numbers RCC2318 and RCC2507. The cell size of both algae is more than 3 μ. However, small *Chaetoceros* are abundant in spring bloom waters, and the simple morphology could hide a high diversity of species [80,81].

Our study identified a limited variety of Haptophyta in the White Sea under-ice water. We found *Phaeocystis pouchetii* and *Chrysochromulina simplex* that are widespread in the plankton of the subarctic and Arctic seas [19,51,82,83]. The cell sizes of the identified *Chrysochromulina* sp. and two uncultured haptophytes are unknown, but they supposedly can match the size of picofraction. Egge et al. [83] discovered six OTUs assigned to *Chrysochromulina* that were only found in the picoplankton size fraction rather than the nanoplankton. Since Haptophyta DNA is known to amplify poorly, the molecular methods can underestimate these algae diversity [16,84]. In the Arctic, the number of haptophyte OTUs found in plankton fraction <3 μm varies significantly, between 10–12 in Beaufort Sea and English Channel [16,19] and 59 in the North Pacific [85].

The phylogenetic analysis of Bolidophyceae sequences from the under-ice plankton of the White Sea showed the highest diversity of these algae among the identified PPE classes. Bolidophytes are considered true picosized forms [86]. Bolidophytes were represented by *Triparma strigata* and 12 OTUs of uncultured forms. A large number of bolidophyte sequences are uncultured forms, as noted earlier for the Arctic or subarctic locations, as well as for the Baltic Sea [14,65,79,82,87]. Previously, ice only two phylotypes of bolidophytes were found in the White Sea [27]. Based on a phylogenetic analysis, the number of White Sea bolidophytes was increased by including sequences from the Genbank that had not been previously classified as Bolidophyceae and had been deposited at the NCBI Genbank as uncultured stramenopile or uncultured eukaryotes [47]. In the under-ice water of the White Sea, four OTUs of bolidophytes have 100% similarity with OTUs found earlier in the ice (1 OTU) and summer plankton (2 OTUs). The phylogenetic analysis recovered that the White Sea bolidophytes refer to three environmental clades (env. clades I, II and III [88] (Figure 4)), in addition to the group corresponding to the genus *Triparma*. Bolidophyte *T. strigata* has a worldwide distribution in plankton and is most abundant in polar waters [89–91]. The complete sequences of the 18S rRNA gene are almost identical (similarity of 99.9%–100%) for such morphologically distinct species as *T. strigata*, *T. laevis f. longispina, T.* aff. *verrucosa*, and flagellate *Triparma* sp. RCC1657; therefore, if only the 18S RNA gene sequences are considered, all those species are combined in the clade of *T. laevis* [90]. Therefore, the discovery of a phylotype similar to *T. strigata* does not necessarily mean that there is only this species in the under-ice plankton of the White Sea. Bolidophyceae sequences dominated among PPEs reads in sample 2/13. This is the first registration of Bolidophyceae domination in the subarctic plankton.

The dominance of *Micromonas polaris* reads in three out of seven samples corresponds to the fact that our studies were carried out in the prebloom period. *M. polaris* had exceptionally high relative read abundances during pre- and postbloom stages in Isfjorden, West Spitsbergen [64], the Amundsen Gulf, and the Canadian Beaufort Sea [68]. An unexpected result was the high share of different *Mantoniella* taxa in three samples and Bolidophyceae in one sample. Situations where the relative read abundance of *M. polaris* was lower than other taxa were noted earlier in different regions of the Arctic and subarctic: the unexpectedly high proportion of *Bathycoccus* was revealed in July surface samples in the Amundsen Gulf and Canadian Beaufort Sea [68]. It was suggested that this might have been associated with offshore upwelling, or, more speculatively, a viral attack on *Micromonas* triggered by specific oceanographic conditions. *Mantoniella squamata* made a great contribution to the Mamiellophyceae reads off Greenland [43]. The dominance of *Mantoniella* from clade A in PPE reads was identified for the first time. The spatial variability of relative read abundances may be controlled by the combined influence of abiotic and biotic factors.

Earlier, at the same stations as in the present work, we studied the taxonomic composition of PPE and protists in ice samples filtered through a 2-μm pore size filter [27]. Different assemblages colonized the under-ice water and the ice. In samples of under-ice water filtered through a 2-μm pore size filter, Stramenopiles made the most significant contribution to total quality reads (average: 34%), whereas Rhizaria dominated in the ice samples (average: 18%). At the same time, the number of identified protists taxa determined to genus level was comparable: 175 in water and 185 in ice. In the Baltic Sea, the ice community was more diverse than the wintertime water [79]. The contributions of Chloroplastida in ice and water were comparable; the Alveolata contribution was lower in ice than in under-ice water. It is interesting to note that Alveolates were not the dominant group in any of the samples of the White Sea plankton, while the domination of Alveolates in the reads abundance was noted in size-fractionated seawater (0.2–3.0 μm) of the Amundsen Gulf flaw lead system [13].

The similarity of PPE composition in under-ice water and ice was 0.75 (Sørensen index). The variety of phylotypes of certain genera in plankton was lower than in ice. For example, the genus *Ostreococcus*

in plankton was represented only by *O. tauri*, while *O. tauri* and *Ostreococcus* sp. were recorded in ice. Pelagophyceae was found in ice but not in under-ice water. Chlorophyta dominated in both habitats, but the contributions of specific genera and classes varied. For example, the domination of *Mantoniella* in water (three stations) was not observed in the ice, where *Micromonas* always made the most significant contribution. In addition, the domination of Bolidophyceae was observed only in under-ice water (one sample), while the domination of Trebouxiophyceae was only found in ice (one station).

The research undertaken on White Sea under-ice photosynthetic picoeukaryotes' genetic diversity is one stage in studying the dynamics of plankton communities in the subarctic.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/3/93/s1. Figure S1. Air temperature at the White Sea Biological Station (shaded area corresponds to sampling period).

**Author Contributions:** Conceptualization, T.A.B. and L.V.I.; methodology, M.D.L. and I.A.M.; validation, T.A.B., L.V.I. and A.V.T.; formal analysis, T.A.B., M.D.L. and A.V.T.; investigation, I.A.M. and T.A.B.; resources, T.A.B. and A.V.T.; data curation, I.A.M. and T.A.B.; writing—original draft preparation, T.A.B.; writing—review and editing, A.V.T. and L.V.I.; visualization, A.V.T. and T.A.B.; supervision, A.V.T.; project administration, A.V.T.; funding acquisition, T.A.B. and A.V.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was done in the framework of the Moscow University Project "Noah's Ark" and of the State Tasks of the Lomonosov Moscow State University (themes no. AAAA-A16-116021660052-0, and AAAA-A17-117120540067-0), with the financial support of the Russian Foundation for Basic Research (project No. 16-05-00502).

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Patterns in Alpha and Beta Phytoplankton Diversity along a Conductivity Gradient in Coastal Mediterranean Lagoons**

**Natassa Stefanidou 1,\*, Matina Katsiapi 1,2, Dimitris Tsianis 3, Maria Demertzioglou 4, Evangelia Michaloudi <sup>4</sup> and Maria Moustaka-Gouni 1,\***


Received: 16 December 2019; Accepted: 17 January 2020; Published: 19 January 2020

**Abstract:** Understanding the diversity patterns of phytoplankton assemblages in coastal lagoons is clearly important for water management. In this study, we explored alpha and beta diversity patterns in phytoplankton communities across five Mediterranean lagoons hydrologically connected to Vistonikos Gulf. We examined the phytoplankton community composition and biomass on a monthly basis from November 2018 to October 2019. For this, water samples were collected from seven inshore, brackish and coastal waters, sampling sites covering a wide range of conductivity. We found significant spatial and temporal differences in phytoplankton alpha diversity and in phytoplankton biomass metrics explained by the high variation of conductivity. Evenness remained low throughout the study period, reflecting significant dominance of several phytoplankton blooms. Harmful algal blooms of *Prorocentrum minimum*, *Alexandrium* sp., *Rhizosolenia setigera* and *Cylindrotheca closterium* occurred. The system's species pool was characterized by relatively high phytoplankton beta diversity (average ~0.7) resulting from high temporal species turnover (90%). Overall, alpha and beta diversity components were indicative of rather heterogeneous phytoplankton communities which were associated with the high differences in conductivity among the sampling sites.

**Keywords:** brackish; HABs; dinoflagellates; diatoms; cyanobacteria; transitional water bodies

#### **1. Introduction**

Coastal lagoons are shallow semi-enclosed dynamic ecosystems which are connected to the sea by one or more restricted inlets [1]. They are transitional zones between the terrestrial and the marine environment occupying 13% of the coastal areas worldwide [2]. Coastal lagoons provide essential ecosystem services [3]. Due to their shallow depth and the restricted water exchange with the adjacent sea, they are naturally high productive ecosystems [4,5] that are usually subjected to intensive human activities, mostly aquaculture, with several consequences [6–8]. In addition to human exploitation, coastal lagoons have been identified as one of the most vulnerable habitants to potential impacts associated with climate change [5,9]. Projections for atmospheric temperature increase up to 6 ◦C in parallel with the more often and stronger hot extremes for the Mediterranean region [10,11] would most probably increase the salinity in coastal lagoons affecting, in turn, the biological communities. Studies using microbial communities as experimental systems have pointed out that a potential salinity

change in coastal Mediterranean environment would significantly affect the microbial biodiversity, biomass and resource use efficiency [12,13].

Within the Mediterranean basin, there are numerous (>100) coastal lagoons of varying size [14,15]. Multiple line of evidence suggest that these ecosystems host an important fraction of global biodiversity including rare species [16–18]. A considerable number of studies examining the alpha diversity patterns (i.e., the local diversity within a site) of macrophytes (e.g., [19]), fish (e.g., [20]) and benthic invertebrates (e.g., [21]) in coastal lagoons relate them to environmental parameters (e.g., total phosphorus). More recently, studies of the spatial alpha diversity gradients of phytoplankton communities relate them to geographical variables (e.g., altitude in [22]) or other physical (e.g., water temperature in [23]) and chemical (e.g., nutrients in [24]) environmental aspects. Furthermore, investigations were carried out in order to determine the key environmental factors driving the frequent phytoplankton blooms in these shallow productive ecosystems [25]. The advancement of Next Generation Sequencing (NGS) tools in recent years has allowed the investigation of deep microbial diversity in coastal lagoons, revealing information for their alpha diversity patterns in relation to other biotic or abiotic parameters (e.g., [26–28]).

However, still today, very little is known about the microbial (including phytoplankton) beta diversity (i.e., the ratio of the local diversity divided by the regional diversity) and how it varies along the conductivity gradient in such ecosystems. Whittaker et al. [29] defined beta diversity as the extent of differentiation among biological communities along habitat gradients within an ecosystem. In particular, beta diversity is a measure of the difference in species composition between two or more local assemblages (e.g., the five Mediterranean lagoons under study) [30]. There are two potential ways in which two or more assemblages can be different: one is species replacement (i.e., turnover) and the second is species loss or gain, which implies that the poorest assemblage is a subset of the richest one (nestedness) [31]. The differentiation of the spatial turnover and nestedness components of beta diversity is crucial for understanding central biogeographic, ecological and conservation issues [32].

In the present study we explored the patterns of phytoplankton biodiversity across five Mediterranean lagoons which are located nearby and are hydrologically connected to Vistonikos Gulf, North Aegean Sea. To the best of our knowledge, this is the first study on these lagoons' phytoplankton. Particularly, we examined the spatial and temporal differences in phytoplankton species richness and other alpha diversity estimators (i.e., Shannon index) along the conductivity gradient that characterizes these lagoons. Furthermore, we explored the degree of homogeneity across the different waterbodies by computing beta diversity. By calculating nestedness and turnover, we examined whether beta diversity originated from differences in richness level (nestedness) or from species replacement (turnover). Finally, in order to study the effect of conductivity on phytoplankton communities we modelled their relationship using different phytoplankton attributes (e.g., phytoplankton biomass, species and group composition richness). Overall, we aimed to evaluate whether hydrological connectivity through water intrusion from Vistonikos Gulf or environmental heterogeneity in terms of conductivity stronger affects the phytoplankton community structure across the five lagoons. In the first case we expect that phytoplankton communities would be relatively homogeneous (i.e., low beta diversity and low turnover) with no significant relationships between conductivity and the different phytoplankton attributes. Alternatively, in the case where conductivity would significantly affect the phytoplankton communities across the lagoons, we expect rather heterogeneous phytoplankton communities (i.e., high beta diversity) resulting from species turnover.

#### **2. Materials and Methods**

#### *2.1. Sampling Sites and Sample Collection*

Phytoplankton samples were collected monthly on an annual basis (November 2018 to October 2019) from five Mediterranean lagoons located in Thrace, Northeastern Greece. All the lagoons are connected to Vistonikos Gulf (North Aegean Sea) by canals (Table 1; Figure 1). In total, 84 phytoplankton samples were analyzed from seven inshore sites of these lagoons. Vistonis lagoon has a surface area of approximately 45 km2 with a mean depth of 3 m and it is connected to Vistonikos Gulf through three tidal canals [33]. Being the largest lagoon among the five studied, samples were collected from three different sites in the southern part of Vistonis lagoon nearby Vistonikos Gulf. The second largest lagoon, Lagos lagoon, has a surface area of 3 km2 and mean depth of 2 m. In its northern part it communicates with the southern part of Vistonis lagoon while in its southern part with Vistonikos Gulf through tidal canals [34]. Lafra and Lafrouda are smaller lagoons than Vistonis and Lagos and have approximately 1.2 km<sup>2</sup> and 0.8 km<sup>2</sup> surface areas, respectively, and a mean depth <0.5 m [35]. The smallest lagoon is Palaia Koiti lagoon which occupies a surface area of <0.1 km2 and has a mean depth of approximately 1m[35]. Lafra, Lafrouda and Palaia Koiti lagoons are connected with Vistonikos Gulf through tidal canals.

**Figure 1.** Study area and sampling sites in Lafra, Lafrouda, Lagos (Lag), Palaia Koiti (PK) and Vistonis (Vis1, Vis2, Vis3) lagoons. The sampling sites are indicated with a red symbol.


**Table 1.** Latitude, longitude, mean annual conductivity (μS cm<sup>−</sup>1) and mean annual salinity (psu) per sampling site.

During all samplings, in situ measurements of water temperature and conductivity were made with the use of a HACH HQ40D portable multi meter. Conductivity was transformed to salinity based on the equation in Weyl [36]. Phytoplankton samples of 0.5 L were collected from the surface layer (0–1 m) and were preserved immediately with Lugol's solution. The samples were kept in the dark until they were microscopically analyzed within the next few days.

#### *2.2. Microscopy Analysis*

Phytoplankton samples were examined using an inverted epi-fluorescence microscope (Nikon Eclipse TE 2000-S, Melville, MSA) with phase contrast. Phytoplankton identification was based on taxonomic keys and papers and counting was done using the inverted microscope method following the same protocol described in Mazaris et al. [37]. For biomass estimation, the dimensions of 30 individuals of each abundant species were measured using a digital microscope camera (Nikon DS-L1). Mean cell volume estimates were calculated using the appropriate geometric formulae [38]. Phytoplankton biovolume for individuals > 180 μm3 was converted into carbon content using the formulae C = 0.288 V0.811 for diatoms and C = 0.216 V0.939 for other phytoplankton taxa (C is carbon content in pg, V is cell volume in μm3) [39]. For individuals < 180 μm3, conversion factors of 0.108 pg C μm−<sup>3</sup> for diatoms and 0.157 pg C μm−<sup>3</sup> for other phytoplankton taxa were used [40].

#### *2.3. Diversity*

Alpha diversity estimators for the phytoplankton community (the Shannon and Evenness indices) were calculated with the PAST3 software [41].

To compute beta diversity of the phytoplankton community for each sampling, we used the 'betapart' R package version 1.5.1 [31]. Baselga's [31,42] approach suggests that Sorensen multiple-site dissimilarity (bSOR) is partitioned into two components: spatial turnover in species composition, measured as Simpson dissimilarity (bSIM) and variation in species composition due to nestedness (bNES) measured as nestedness-resultant fraction of Sorensen dissimilarity. Furthermore, Sorensen pair-wise dissimilarity (bSOR) between site pairs was computed in order to use the values for the construction of non-linear regression fitting models. The above analyses were run in R environment version 3.5.3 [43].

#### *2.4. Non-Linear Regression Models*

To explore the relationship between conductivity and phytoplankton metrics (abundance, biomass, phytoplankton C, species richness, Evenness and Shannon index) we created non-linear regression models using conductivity as the predictor variable and phytoplankton metrics as the response variable. Models were fit to the data according to the polynomial equation *y* = α1*X*<sup>2</sup> + *a*2*X* + *a*3, where y is the response variable, X is the predictor variable and α1, α2 and α3 the model coefficients estimated by maximum likelihood. The statistical significance of the coefficients in the constructed models was assessed by their individual p-value; a p-value < 0.05 indicated a significant coefficient for the predictor variable. The overall significance of each constructed model was assessed by the F-statistic and its

p-value. Finally, the strength of the relationship between the predictor and the response variable was determined through the multiple R2. All the included data for the construction of the models were log(x + 1) transformed in order to produce the smallest error possible when making a prediction and improve the fit of the model [44].

To explore the relationship of beta diversity with conductivity we constructed a model using the difference in conductivity (absolute values) between pair of sites and Sorensen pair-wise dissimilarity (bSOR). The polynomial equation was also used to fit the data, the statistical significance of the coefficients was evaluated by their p-value and the overall significance of the model by the F-statistic and its p-value. Multiple R2 was used to examine the strength of the relationship between the difference in conductivity and pair-wise bSOR.

#### *2.5. Statistical Analysis*

Multidimensional scaling (MDS) based on Jaccard similarity coefficient of presence-absence data was used to depict the similarities of phytoplankton communities among the different lagoons. Significant differences in phytoplankton community structure among the samples were evaluated with one-way analysis of similarities (ANOSIM) and similarity percentage (SIMPER) analysis was run to determine the percentage contribution of phytoplankton species in dissimilarity between the site pairs. Finally, we used a one-way permutational multivariate analysis of variance (PERMANOVA) to test the significance of conductivity on the phytoplankton community composition and biomass among the different sampling sites. The above statistical analyses were run in PAST 3.

#### **3. Results**

#### *3.1. Phytoplankton Community Structure and Biomass*

A total of 77 phytoplankton species were identified in all the sampling sites during one year (November 2018 to October 2019) which were affiliated to 10 taxonomic groups (Supplementary Table S1). Bacillariophyceae (diatoms) had the highest overall species richness among the taxonomic groups with 31 representatives, followed by Dinophyceae (17 species), Chlorophyceae (10 species) and Cyanobacteria (8 species) (Supplementary Table S1). The rest of the taxonomic groups (Cryptophyceae, Haptophyceae, Prymnesiophyceae, Prasinophyceae, Euglenophyceae & Raphidophyceae) contributed with ≤4 species each to total phytoplankton richness. Diatoms were also the most diverse group in terms of species richness in all the sampling sites. Palaia Koiti lagoon with mean annual conductivity 49,468.2 μS cm−<sup>1</sup> (mean annual salinity 30.7 psu; Table 1) exhibited the highest diatom richness (28) while Vis1 in Vistonis lagoon with mean annual conductivity 11,664.6 μS cm−<sup>1</sup> (mean annual salinity 7.3 psu) the lowest (13) (Table 2). Dinoflagellates had higher richness in Lafra, Lafrouda, Lagos and Palaia Koiti lagoons with a mean annual conductivity > 40,000 μS cm−<sup>1</sup> (mean annual salinity > 28 psu). Chlorophytes and cyanobacteria appeared to be more diverse in Vistonis lagoon with mean annual conductivity <sup>≤</sup> 30,000 <sup>μ</sup>S cm−<sup>1</sup> (mean annual salinity <sup>&</sup>lt; 18.5 psu).

**Table 2.** Number of species for each phytoplankton taxonomic group in Lafra, Lafrouda, Lagos (Lag), Palaia Koiti (PK) and Vistonis (Vis1, Vis2, Vis3) lagoons.


Phytoplankton biomass varied highly throughout the year in all the sampling sites with the highest biomass values being most commonly measured during spring and summer and the lowest during winter (Figure 2). The highest variability was recorded in Vis1 where the maximum biomass (134.7 mg L−1) was measured in September 2019. This high biomass was made up by the diatom *Rhizosolenia setigera* forming a bloom (12,315 cell mL<sup>−</sup>1) at a conductivity of 16,090 μS cm−<sup>1</sup> (9.7 psu). The minimum biomass (4.4 mg L<sup>−</sup>1) in Vis1 was measured in January 2019 dominated by the diatom *Cyclostephanos* sp. at a conductivity of 2620 μS cm−<sup>1</sup> (2.6 psu). Similarly, high variability in phytoplankton biomass was recorded in Lafrouda lagoon reaching a maximum biomass of 127.8 mg L−<sup>1</sup> in April 2019, concurrent with a *Prorocentrum minimum* bloom (301 cell mL<sup>−</sup>1) at conductivity 46,400 μS cm−<sup>1</sup> (31.6 psu). The minimum biomass (0.2 mg L−1) in this lagoon was measured in December 2018 and was mainly made up by *Cyclostephanos* sp. at a conductivity of 25,700 μS cm−<sup>1</sup> (8.6 psu).

**Figure 2.** Total phytoplankton biomass in Lafra, Lafrouda, Lagos, Palaia Koiti and Vistonis lagoons during the period from November 2018 to October 2019.

A comparable high variability was noticed for the phytoplankton biomass proxy (phytoplankton C) in all the sampling sites (Supplementary Figure S1a). The widest range was observed in Lafrouda lagoon with a maximum value of 16,448.5 μgCL−<sup>1</sup> coinciding with the highest biomass (April 2019) and a minimum of 16.9 μgCL−<sup>1</sup> coinciding with the lowest biomass (December 2018). Vis3 followed with a similar wide range and particularly a maximum value of 12,851.4 μgCL−<sup>1</sup> was recorded in March 2019 at a conductivity of 23,900 μS cm−<sup>1</sup> (14.6 psu). The minimum value of phytoplankton C (417.6 μgCL<sup>−</sup>1) in Vis3 was observed in October 2019 at a conductivity of 31,600 μS cm−<sup>1</sup> (18.6 psu).

#### *3.2. Phytoplankton Diversity*

#### 3.2.1. Alpha Diversity Estimators

The number of phytoplankton taxa varied among samples between 7 in December 2018 in Lafrouda lagoon at a conductivity of 25,700 μS cm−<sup>1</sup> (8.6 psu) and 29 in June 2018 in Lagos lagoon at a conductivity of 50,600 μS cm−<sup>1</sup> (32.7 psu) (Supplementary Table S2). We found an overall increasing trend in the total species richness over spring and summer in all the sampling sites with the lowest number of taxa being most commonly recorded during winter. The highest number of phytoplankton taxa as well as the highest mean annual number of taxa were found in Lagos (29 taxa) and Palaia Koiti (27 taxa) while the lowest in Lafrouda (7 taxa) (Supplementary Figure S1b). Likewise, the Shannon diversity index showed an increasing trend over spring and summer compared to winter months. The highest values along with the highest mean annual values were computed for Lagos and Palaia Koiti (Supplementary Figure S1c). Evenness remained relatively low throughout the study period in all the sampling sites reflecting significant dominance by few taxa (Supplementary Figure S1d). Hence, the lowest values of Evenness were observed at the same time as heavy phytoplankton blooms (e.g., the case of *Prorocentrum minimum* bloom in April 2019 in Lafrouda lagoon).

#### 3.2.2. Beta Diversity

Phytoplankton beta diversity ranged from 0.6 (January and October 2019) to 0.76 (December 2018) and did not display a seasonal pattern similar to the majority of alpha diversity estimators (e.g., Species richness and Shannon index). Community compositional variation was almost entirely attributed to species turnover (0.60 ± 0.05 SE) rather than nestedness (0.07 ± 0.05 SE; Figure 3). Both turnover and nestedness did not show high fluctuation under different conductivity levels.

**Figure 3.** Phytoplankton beta diversity according to Sorensen similarity index (bSOR), turnover (bSIM) and nestedness (bNES) components values in Lafra, Lafrouda, Lagos, Palaia Koiti and Vistonis lagoons during the period from November 2018 to October 2019.

#### *3.3. Phytoplankton Grouping*

MDS showed a separation of the phytoplankton community composition based on mean annual conductivity (Figure 4). The majority of the samples from Vistonis lagoon, with a mean annual conductivity <sup>≤</sup> 30,000 <sup>μ</sup>S cm−<sup>1</sup> (mean annual salinity <sup>&</sup>lt; 18.5 psu), were placed together and separately from the samples of Lafra, Lafrouda, PK and Lag sampling sites with mean annual conductivity > 40,000 μS cm−<sup>1</sup> (mean annual salinity > 28 psu). This separation was supported by one-way PERMANOVA which detected significant effects of conductivity (F = 2.0, *p* < 0.01) on the phytoplankton community structure. One-way ANOSIM also indicated significant differences (R = 0.2, *p*same < 0.01) since higher similarities in the phytoplankton community composition were found within the same site in two different time points than between two different sites. According to SIMPER analysis, four to nine species were responsible for>70% of the compositional differences between the three sampling sites of Vistonis lagoon (Vis1, Vis2, Vis3) and Lafra, Lafrouda, PK and Lag with *Prorcentrum minimum* contributing with>20% to the compositional differences in all cases (Supplementary Table S3). The other species discriminating the phytoplankton communities among the sampling sites were the diatoms *Chaetoceros* sp., *Cyclostephanos* sp., *Rhizosolenia setigera* and *Skeletonema* cf. *costatum*, the dinoflagellates *Alexandrium* sp. and *Gonyaulax verior* as well as the coccoid cyanobacteria.

**Figure 4.** Multidimensional scaling plot of phytoplankton biomass variation in Lafra, Lafrouda, Lagos (Lag), Palaia Koiti (PK) and Vistonis (Vis1, Vis2, Vis3) lagoons according to Jaccard similarity index. Samples with mean annual conductivity lower than 30,000 μS cm−<sup>1</sup> are represented with a light blue rhomb and samples with mean annual conductivity higher than 40,000 μS cm−<sup>1</sup> are represented with a dark blue triangle. 2D Stress value: 0.16.

Regarding the biomass values of the discriminating species, *Prorocentrum minimum* reached the extremely high biomass of 124.6 mg L−<sup>1</sup> in April 2019 in Lafrouda lagoon with a conductivity of 46,400 μS cm−<sup>1</sup> (31.6 psu). Occasionally, very high biomass was measured for this species in the three sampling sites of Vistonis lagoon; 55.6 mg L−<sup>1</sup> in May 2019 in Vis1 at conductivity 15,910 μS cm−<sup>1</sup> (9.6 psu) and 26.1 mg L−<sup>1</sup> in August 2019 in Vis3 at conductivity 34,900 μS cm−<sup>1</sup> (22 psu). On the contrary, the biomass of *Prorocentrum minimum* remained rather low (<2.5 mg L<sup>−</sup>1) in Lagos and Palaia Koiti lagoons with mean annual conductivity > 49,000 μS cm−<sup>1</sup> (>29 psu). In contrast to *Prorocentrum minimum*, *Alexandrium* sp. and *Gonyaulax verior* were detected only in one or two samples of Vistonis lagoon throughout the study period. However, they were recorded more frequently and with higher biomass in Palaia Koiti, Lagos, Lafra and Lafrouda lagoons. The maximum biomass of *Alexandrium* sp. (81.4 mg L<sup>−</sup>1) was measured in July 2019 in Lafrouda at a conductivity of 47,400 μS cm−<sup>1</sup> (30.9 psu) and of *Gonyaulax verior* (3.7 mg L<sup>−</sup>1) at the same month in Lafra at a conductivity of 36,800 μS cm−<sup>1</sup> (25.5 psu). Even though the diatoms *Chaetoceros* sp., *Cyclostephanos* sp., *Rhizosolenia setigera* and *Skeletonema* cf. *costatum* were detected in the majority of the samples in the seven sampling sites, they all reached their maximum biomass in Vistonis lagoon, except for *Skeletonema* cf. *costatum*, which had higher biomass in Lagos and Palaia Koiti lagoons. In particular, the maximum biomass of *Chaetoceros* sp. (9.3 mg L<sup>−</sup>1) was measured in November 2018 in Vis3 at a conductivity of 44,850 μS cm−<sup>1</sup> (17.2 psu), of *Cyclostephanos* sp. (20.5 mg L−1) in December 2018 in Vis2 at a conductivity of 17,457 μS cm−<sup>1</sup> (1.7 psu), and of *Rhizosolenia setigera* (123.2 mg L−1) in September 2019 in Vis1 at a conductivity of 16,090 μS cm−<sup>1</sup> (9.7 psu). For *Skeletonema* cf. *costatum* the maximum biomass (1.4 mg L−1) was measured in August 2019 in PK at a conductivity of 49,300 μS cm−<sup>1</sup> (26.4 psu). Regarding the coccoid cyanobacteria, they were also detected in the majority of the samples. Nevertheless, their biomass was higher in Vistonis lagoon where they reached their maximum value (6.7 mg L<sup>−</sup>1) in May 2019 in Vis1 at a conductivity of 15,910 μS cm−<sup>1</sup> (9.6 psu). Furthermore, in Vistonis lagoon, we measured the maximum biomass for other two cyanobacterial species *Dolichospermum flos-aquae* in September 2019 at a conductivity of 16,090 μS cm−<sup>1</sup> (9.7 psu) and *Aphanizomenon favalaroi* in November 2018 at a conductivity of 11,716 μS cm−<sup>1</sup> (12.8 psu). It is worth noting that *Dolichospermum flos-aquae* and *Aphanizomenon favalaroi* were not detected in the samples from the other four lagoons (Lafra, Lafrouda, Lagos, Palaia Koiti) with the exception of two samples in Lagos lagoon where only few individuals were found.

#### *3.4. Non-linear Regression Models*

Five out of the seven constructed models demonstrated significant relationships (*p*-value < 0.05) between conductivity (predictor) and phytoplankton variables (phytoplankton abundance and biomass, phytoplankton C, species richness and pair-wise bSOR) (Table 3). A negative relationship was noticed in all five cases between the predictor and responses variables and in particular, we found that a potential increase in conductivity will lead to significant negative effects on phytoplankton abundance and biomass, phytoplankton C, species richness and pair-wise bSOR (Supplementary Figure S2). The strongest relationship (R<sup>2</sup> of 0.4) was detected between conductivity and phytoplankton abundance while the weakest (R2 = 0.1) was between conductivity and pair-wise bSOR. The lowest residual standard error (0.1) was computed for the constructed model between conductivity and species richness since the majority of our observation data were placed above or very close to the best fit (Supplementary Figure S2).

**Table 3.** Statistics of the constructed non-linear regression models. Significant relationships (*p*-value < 0.05) are in bold.


#### **4. Discussion**

The significant role of salinity/conductivity in determining the composition and diversity of microbial communities in aquatic ecosystems is supported by a series of studies carried out in coastal lagoons (e.g., [45]), in marine environments (e.g., [12]) as well as in lakes (e.g., [46]). Hydrological connectivity has also been identified as an emergent driver for species homogenization between two or more adjacent sites [47]. After discussing how conductivity affected the phytoplankton community structure and diversity across the five studied lagoons, we will discuss the degree of homogeneity (beta diversity) among the communities. Thus, we address the question of whether environmental heterogeneity in terms of conductivity or hydrological connectivity is the primary driver for shaping the less studied phytoplankton communities in transitional water bodies.

#### *4.1. Phytoplankton Composition and Alpha Diversity*

In our study, diatoms were the most diverse group with the highest species number in all five lagoons. This result corroborates previous findings on phytoplankton community structure in other Mediterranean lagoons where diatoms supported the highest taxonomic richness [22,24]. Commonly, diatoms are considered as an important component of phytoplankton communities in estuaries and coastal waters being typically the dominant taxonomic group in terms of species richness in such environments [48]. This success and dominance may be attributed to the higher inherent plasticity that several diatom genera (such as *Skeletonema*, *Chaetoceros*, *Thalassiosira*) have in comparison to other phytoplanktonic genera such as the dinoflagellates *Alexandrium* and *Gonyaulax* or the cyanobacterium *Aphanizomenon* and the chlorophyte *Monoraphidium*. Bussard et al. [49] associated the phenotypic plasticity of diatoms under different environmental conditions (including a salinity gradient) with their ability to change their gene expression pathways. They mentioned that this taxonomic group induces the expression of stress-related genes under stress environmental conditions. This characteristic justifies the higher acclimation capability of diatoms to environmental fluctuations in comparison to other phytoplankton taxonomic groups. On the other hand, dinoflagellates due to low inherent plasticity together with higher marine affinity were more diverse in the lagoons with relatively high mean annual conductivity/salinity (> 40,000 μS cm−<sup>1</sup> and > 28 psu). This is in contrast to chlorophytes and cyanobacteria which are characterized by higher freshwater affinity in brackish waters [50] and appeared to be more diverse in Vistonis lagoon (< 30,000 μS cm−<sup>1</sup> and < 18.5 psu).

Conductivity significantly determined the phytoplankton composition and biomass according to one-way PERMANOVA in line with previous studies in similar ecosystems (e.g., [51,52]). In addition, there was very high variation in the values of the two alpha diversity indices with annual mean values of the Shannon index being double in Palaia Koiti lagoon compared to Lafrouda lagoon. Conductivity did not account or only weakly (non-significantly) accounted for differences in alpha diversity indices across the sampling sites (based on the relationships as demonstrated in the constructed non-linear regression models). Thus, most probably, the different typological descriptors (size, morphometry) of the lagoons were responsible for the high variation in species richness. Smith et al. [53] analyzed data from 142 different natural ponds, lakes, and oceans and 239 experimental ecosystems that revealed a strong phytoplankton species–area relationship. Although there is broad agreement on the positive effect of diversity on community productivity [54,55], during our survey, the lowest values of alpha diversity indices were concurrent with the highest phytoplankton biomass. These findings are well in line with the humped-back biomass–species richness relationship in phytoplankton communities due to the habitat and trophic state of the ecosystem [56]. The resultant biomass–species richness relationship indicates that heavily eutrophic environments as our studied lagoons are expected to have strong dominance and not the optimum species richness and diversity. According to Dodson et al. [57] an optimum of phytoplankton species richness is found at intermediate productivity for lake areas <100 km2. In the majority of the cases where very high phytoplankton biomass values was measured, it was attributed to the very high proliferation of one or two species (i.e., *Prorocentrum*

*minimum* contributed with 97.5% to the total phytoplankton biomass of 127.8 mg L−<sup>1</sup> at Lafrouda lagoon measured on April 2019).

*Prorocentrum minimum* formed blooms in all studied lagoons. This is a harmful, common bloom-forming species in eutrophic, low to moderate salinity coastal and brackish environments [58]. In the Mediterranean region, this species can proliferate in a wide range of water temperatures (4–27 ◦C) showing great ability to adapt to widely varying natural conditions as shown in our study. *Alexandrium*, also formed blooms at high biomass at different salinities. The ability of *Alexandrium* to colonize multiple habitats and its adaptability and persistence over large region has been recognized in association with the ability of its species to produce toxins [59]. Blooms of the other euryhaline species, the mucilage producing diatoms *Cylindrotheca closterium*, *Rhizosolenia setigera* and *Chaetoceros* sp. have been observed in a eutrophic coastal sea contributing to frequent mucilage events [60].

#### *4.2. High Variation in Conductivity Promotes High Phytoplankton Beta Diversity*

Beta diversity is considered positively dependent on environmental heterogeneity since an increase in the latter incorporates an increase in the variety of environmental conditions to which different species are adapted, hence promoting greater variation in species composition among local communities within a region [61,62]. Studies on phytoplankton communities from different aquatic environments have provided consistent evidence and relate the enhanced beta diversity to environmental gradients within a region [63,64]. Under different circumstances, environmental homogeneity has been shown to generate biotic homogenization and thus, low beta diversity of biological communities, including phytoplankton [65]. In aquatic ecosystems, hydrological connectivity has been associated with environmental homogeneity and hence, with species dispersal and high shared diversity throughout a region [66–68]. Even though hydrological connectivity is supported for the five lagoons connected to Vistonikos Gulf by tidal canals, according to our analyses, instead of finding evidence of biotic homogenization (i.e., low beta diversity and low turnover), we calculated relatively high beta diversity for the entire monitoring period. Community compositional variation was almost entirely attributed to species turnover (i.e., species replacement among the different sampling sites). This result contradicts our hypothesis for biotic homogenization and decreased beta diversity, as a consequence of hydrological connectivity of the lagoons.

We assume that conductivity contributed to the species replacement (turnover) between the sampling sites since a significant effect was detected between conductivity and the pair-wise beta diversity (according to the constructed non-linear regression models). Hence, our hypothesis for heterogeneous phytoplankton communities is clearly supported, reflecting the direct effects of environmental heterogeneity, by means of conductivity in our case, on species turnover according to Soininen et al. [69]. It seems that lagoon typological factors (size, morphometry and especially the share of water volume or surface area in contact to the sediment) and local environmental variables were significant determinants of the community structure in agreement with previous studies [62,70]. Our study provides evidence that habitat fragmentation together with environmental heterogeneity constitute an important part of the theory of metacommunity organization [71].

Understanding the relationships between biological communities and environmental aspects remains a challenge for climate research [72] and at the same time, it is essential for establishing water management policies [70]. Therefore, during this study, we modeled conductivity with several phytoplankton metrics in order to determine the type of relationship (negative or positive) between them. As demonstrated by the constructed non-linear models, all phytoplankton metrics (phytoplankton C, abundance, biomass, species richness and pair-wise beta diversity) consistently exhibited a negative relationship with conductivity, implying that a potential increase in the latter will negatively affect them. Nevertheless, this is a first attempt to predict the phytoplankton response of five Mediterranean lagoons to enhanced salinity due to climate change and longer time series in the future will allow us to strengthen the accuracy of the model and of projections [18].

#### **5. Conclusions**

Conductivity/salinity was a key environmental factor in shaping phytoplankton communities across five Mediterranean lagoons connected to the Vistonikos Gulf. In particular, we showed that environmental heterogeneity, by means of conductivity, led to high phytoplankton beta diversity due to high turnover exceeding the effects of hydrological connectivity. This fact implies that local conditions in each lagoon were important in determining which species were capable not only of maintaining a viable population but also of forming a heavy bloom. Nevertheless, the exposure of phytoplankton to continuous salinity fluctuations had as a consequence the increase of euryhaline harmful species. Overall, our results provide useful information for ecological water management of transitional eutrophic water bodies in the Mediterranean region.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/1/38/s1. Supplementary Figure S1: Boxplots of phytoplankton C (a), species richness (b), Shannon index (c) and Evenness (d) in the lagoons Lafra, Lafrouda, Lagos (Lag), Palaia Koiti (PK) and Vistonis (Vis1, Vis2, Vis3). Supplementary Figure S2: Scatter plots conductivity and the different phytoplankton metrics. Supplementary Table S1: Species list of phytoplankton taxa per sampling site during the study period. Supplementary Table S2: Species number, total phytoplankton biomass and alpha diversity indices per sample. Supplementary Table S3: List of higher contributing phytoplankton species in communities' dissimilarity determined by SIMPER between site pairs.

**Author Contributions:** Conceptualization, N.S. and M.M.-G.; Sampling D.T., Microscopy Analysis, N.S., M.K.; Data analysis N.S.; Resources, M.M.-G.; Writing-Original Draft Preparation, N.S. and M.M.-G.; Writing Review and Editing N.S., M.M.-G., M.K., E.M., M.D., D.T.; Visualization, N.S., M.D.; Supervision, M.M.-G.; Project Administration, M.K., D.T. and M.M.-G.; Funding Acquisition, M.M.-G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by PREFECTURE OF EASTERN MACEDONIA AND THRACE, grant number ELKE-AUTh: 97071.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **E**ff**ect of Hydrographic Variability on the Distribution of Microbial Communities in Taiwan Strait in Winter**

**Gwo-Ching Gong 1,2, Hsin-Ming Yeh 3, Yu-Kai Chen 3, Chih-hao Hsieh 4, Pei-Chi Ho <sup>4</sup> and An-Yi Tsai 1,\***


Received: 2 September 2019; Accepted: 11 October 2019; Published: 14 October 2019

**Abstract:** This study investigated the spatial variation in the components of a microbial food web (viruses, picoplankton, nanoflagellates, and ciliates) in different hydrographic environments in the Taiwan Strait during winter. Water temperature and salinity varied spatially, with lower temperatures (15.3–22.8 ◦C) and salinities (32.2–33.4 psu) in the northern part of the Taiwan Strait, largely affected by runoff from the coast of China. Concentrations of nutrients and Chl *a* were significantly higher in the northern part than that in the southern part of the study area. *Synechococcus* spp., nanoflagellate, and ciliate abundance also varied significantly, with the northern strait having higher abundances of these communities. In contrast, a higher abundance of bacteria was found in the southern part of the Taiwan Strait. The results of this study, which describes two different ecosystems in the Taiwan Strait, suggest that during winter, a "viral loop" might play an important role in controlling bacterial production in the southern part of the Taiwan Strait, while nanofalgellate grazing of picophytoplankton may contribute mainly to the flux of energy in the northern part.

**Keywords:** microbial food web; bacteria; nanoflagellates; ciliates; Taiwan Strait

#### **1. Introduction**

The microbial community, which consists of bacteria, nanoflagellates, and microzooplankton, plays an important role in the life of aquatic environments [1]. An understanding of the relationship between the distributions of microbial organisms and their environments can help us better understand their eco-physiological requirements. Bacteria are able to use dissolved organic matter and incorporate it into their own biomass. Their biomass may be reduced in several ways, with the important one being grazing by phagotrophic protists (nano- and microzooplankton) and viral lysis [2–5]. These studies showed that both viral lysis and grazing can cause significant mortality, but that the impact of each varies by season, host organism, and environmental conditions. Seasonally, viral lysis is reported to be the main cause of bacterial mortality during the winter, and thus a "viral loop" might serve as an important control mechanism for bacterial production in colder seasons, while grazing can be the cause of most bacterial mortality during warmer seasons [5].

The growth of bacteria can also be controlled by physical factors, such as temperature and the supply of substrates [6–8]. Li [9] reported that their regulation can be seasonal and suggested that, over the course of a year, temperature is the dominant factor affecting bacterial growth in colder seasons. Other factors, such as substrate supply, may be important in warmer seasons. Ochs et al. [10] found that bacterial growth rates are unrelated to temperatures above 14 ◦C, which was further confirmed by Shiah and Ducklow [11].

Understanding the principles that govern trophic interactions among microbial planktonic organisms, those that regulate carbon and energy flows in microbial food webs, is an important goal in biological oceanography. The dominating predators of bacteria are heterotrophic nanoflagellates (HNFs) [3,12,13]. Some studies reported increases in HNFs in areas where there are increases in bacteria, suggesting that HNFs are resource controlled [12,14]. Other studies have found a weak coupling between bacteria and HNF density and recommended predation control [15,16]. Although correlative data can suggest possible governing factors, such data do not provide much insight into the complex mechanisms involved. For example, the lack of a strong relationship between HNFs and bacterial abundance might be related to the possibility that HNFs use other sources of carbon (e.g., picophytoplanktonic cyanobacteria) [13], or it may be that there are other more important predators of bacteria (e.g., ciliates) [17]. Other studies suggested that bottom-up control dominates HNF abundance [18,19], while other studies noted that top-down control is more important [16,20]. In addition, the relative importance of these controls may change along with the environmental trophic state, season, and other environmental factors [21,22].

Taiwan is located at the junction of Ryukyu and Luzon in the northwestern part of the Pacific Ocean. This large island country is bordered to the west by the shallow Taiwan Strait and to the east by the deep Philippine Sea. The broad (140 to 200 km) Taiwan Strait Shelf is located between the coasts of China and Taiwan and extends for more than 300 km. Generally, currents in the Taiwan region vary consistently from season to season [23]. Jan et al. [24] described the seasonal variations in currents in the Taiwan Strait in detail. During winter, the northward flow of South China Sea waters through the Penghu Channel is slowed down by the northeasterly winter monsoon and turns northwestward along the southern edge of the Changyun Ridge. In the northern area of the Taiwan Strait, the cold and fresh southward flow of China's coastal waters turn cyclonically at the northern edge of the Changyun Ridge, where they form a cyclonic cold eddy. These winter differences in physicochemical and biotic factors in the shallow Taiwan Strait may set the scene for differences in the development of microorganisms and provide an excellent environmental frame to analyze the potential changes in microbial food web components under very different trophic environments. We hypothesize that the spatial variation in temperature affects the microbial dynamics in the shallow Taiwan Strait in cold seasons. Therefore, this study investigated the surface water distributions of viral, heterotrophic bacterial, *Synechococcus* spp., nanoflagellate, and ciliate abundance, as well as possible HNF-bacterial coupling along different cruise transects in the Taiwan Strait during winter.

#### **2. Material and Methods**

#### *2.1. Study Sites and Samplings*

Samples were collected in January and February 2018 from the surface waters at 25 established stations (9 transects; A–I) located in the Taiwan Strait on board the Fishery Researcher I (Figure 1A). During this cruise, seawater samples were collected at depths of 2 m using a SeaBird CTD (Sea-Bird Scientific, Bellevue, WA, USA) and General Oceanic Rosette assembly with 10 L Niskin bottle (General Oceanic, Florida, USA) Temperature and salinity were measured with a SeaBird CTD. Water samples were filtered (25 mm GF/F: Whatman glass microfiber filters) for Chl *a* analysis and measured after extraction using an in vitro fluorometer (Turner Design 10-AU-005) (Turner Designs, Inc, San Jose, CA, USA) [25]. Nutrients in seawater samples were measured as previously described by Gong et al. [25]. Sample volumes of 50 mL were used to count bacteria, *Synechococcus*

spp., and nanoflagellates. They were fixed with gluteraldehyde at a final concentration of 0.5%. For ciliates, 500 mL water samples from the surface were fixed with neutralized formaldehyde (2% final concentration) [26] and preserved at 4 ◦C until analysis. Virus populations were identified and measured using flow cytometry (FCM) (BD FACSCaliburTM) (Biosciences, Macquarie Research Park, Sydney, Australia) [27].

**Figure 1.** Spatial variations of temperature (**A**), salinity (**B**), NO3 (**C**), and Chl *a* concentrations (**D**) during the study period in the Taiwan Strait.

#### *2.2. Microbial Counts*

Bacteria, *Synechococcus* spp., and nanoflagellates were counted using an epifluorescence microscope (Nikon Optiphot-2) (1000×). Subsamples of 1–2 mL, 4 mL, and 20 mL were filtered onto 0.2 μm or 0.8 μm black Nuclepore filters for bacteria, *Synechococcus* spp., and nanoflagellates, respectively. Samples were stained with DAPI (4 ,6-diamidino-2-phenylindole) at a final concentration of 1 μg mL−<sup>1</sup> [28] to count bacteria and heterotrophic nanoflagellates (HNFs). Pigmented nanoflagellates (PNFs) and HNFs were detected and counted based on the absence or presence of chlorophyll autofluorescence using a separate filter set optimized for chlorophyll or DAPI under a 1000× epifluorescence microscope (Nikon Optiphot-2). Bacteria and HNFs were identified using their blue fluorescence under UV illumination. *Synechococcus* spp. and PNFs were identified using their orange and red autofluorescences under blue excitation light. Furthermore, to obtain reliable ciliate counts, a 500-mL water sample was concentrated into a 100-mL subsample using a 20-μm mesh size net, and then the subsamples (100 mL) were allowed to settle in an Utermöhl chamber (Utermöhl). The entire area of the Utermöhl chamber was examined at 200× or 400× magnification using an inverted microscope (Nikon-TMD 300). To perform FCM analysis, we diluted them to 1:10 with a 0.2-μm filtered TE buffer (Tris-EDTA buffer) (10 mM Tris, 1 mM EDTA: Ethylenediaminetetraacetic acid), stained them with SYBR Green I solution (1:500 dilution; Molecular Probes), and incubated them at 80 ◦C in the dark for 10 min [29]. Cell side scatter (a proxy of cell size) and SYBR Green fluorescence (an indicator of nucleic acid) were used to distinguish viruses from heterotrophic bacteria [29].

All statistical analyses were performed using SPSS software (SPSS 14). A Student *t*-test was carried out to check for significant differences in microbial community (viruses, bacteria, *Synechococcus* spp., HNFs, PNFs, and ciliates) and environmental factors (temperature, salinity, NO3, and Chl *a*) between the northern and southern parts of the Taiwan Strait, with a confidence level of 95% (α = 0.05). Regression analyses were used to explain the relationship between microbial community and environmental factors. This study used redundancy analysis (RDA; rda function in R package *vegan*) to examine the environmental conditions to explain the microbial community composition in the studied parameters and the spatial variation of the microbial community. To identify the main patterns of spatial change of the relationship between microbial community and environmental conditions, we conducted a series of RDAs using (1) all datasets of the Taiwan Strait, (2) northern part, and (3) southern part of the Taiwan Strait.

#### **3. Results**

#### *3.1. Environmental Conditions*

The spatial variations in hydrographic parameters (temperature, salinity, and NO3 concentration) in the surface waters of the Taiwan Strait are shown in Figure 1. During the study period, water temperature and salinity were found to have different patterns in different areas, with the northern part of the Taiwan Strait having lower temperatures (15.3–22.8 ◦C) and salinities (32.2–33.4 psu) (transects A–D) (Figure 1A,B), as they were affected by China's coastal waters. Temperatures and salinities were higher in the southern part (transects E–I) (*t*-test, *p* < 0.05, Figure 2A). NO3 concentrations were higher (>4 μM) in the northern part compared to the southern part of the study area (Figures 1C and 2B). Chl *a* concentrations were highest (1.6 mg m−3) at the edge of China's coastal waters in the northern part (Figure 1D). In general, the northern and southern parts of the Taiwan Strait had significant differences in average Chl *a* concentrations (*t*-test, *p* < 0.05) (Figure 2B).

**Figure 2.** Mean values were determined from the transects A–D (-) and transects E–I () for temperature and salinity (**A**) and for Chl *a* and NO3 (**B**).

#### *3.2. Spatial Changes in Microbial Community*

During the study period, the spatial variations in viral abundance ranged from 0.4 <sup>×</sup> 107 to 2.4 <sup>×</sup> <sup>10</sup><sup>7</sup> viruses mL−1, was slightly lower in the northern (transects A–D) part (1.2 <sup>±</sup> 0.5) <sup>×</sup> 107 viruses mL−1) than that in the southern (transects E–I) part of the Taiwan Strait, (1.5 <sup>±</sup> 0.6) <sup>×</sup>107 viruses mL−1) (Figure 3A). Bacterial abundance varied between 0.4 <sup>×</sup> 106 and 1.1 <sup>×</sup> 106 cells mL−<sup>1</sup> and 0.6 <sup>×</sup> 106 and 2.9 <sup>×</sup> 106 cells mL−<sup>1</sup> in the northern (transects A–D) part and southern (transects E–I) part of the Taiwan Strait, respectively (Figure 3B). Bacterial abundances in the northern part (transects A–D) were significantly lower than those in the southern part (transects E–I) (Figure 4). Based on our field studies, bacterial abundance was probably related to temperature, and temperature explained about 42% of the difference in bacterial abundance (Figure 5A). Furthermore, viral abundance was significantly positively related to bacterial abundance in the northern part (r = 0.41, *p* < 0.05) and southern part (r = 0.82, *p* < 0.05) of the Taiwan Strait (Figure 5B). Abundance of *Synechococcus* spp. ranged from 0.4 <sup>×</sup> <sup>10</sup><sup>4</sup> to 4.2 <sup>×</sup> 104 cells mL−<sup>1</sup> and showed higher values in the northern part (transects A–D) (Figures 3C and 4). Nanoflagellate concentrations were higher in the northern part (transects A–D), HNF concentrations were (0.6–2.8) <sup>×</sup> 10<sup>3</sup> cells mL−1, and PNF concentrations were (0.3–2.3) <sup>×</sup> 10<sup>3</sup> cells mL−<sup>1</sup> (Figure 3D,E and Figure 4). When *Synechococcus* spp. and PNFs from all stations were pooled, they were found to be positively correlated with total Chl *a* values (*Synechococcus* spp.: r = 0.77, *p* < 0.05; PNF: r = 0.63, *p* < 0.05; Figure 6). Mean *Synechococcus* spp. and PNFs represented a significant part of the total phytoplankton biomass. The same pattern was noted for ciliate abundances and the nanoflagellate community, with the southern part of the Taiwan Strait having lower ciliate concentrations (<400 cells L<sup>−</sup>1) (Figures 3F and 4).

**Figure 3.** Spatial variations of viral (**A**), bacterial (**B**), *Synechococcus* spp. (**C**), heterotrophic nanoflagellate (**D**), pigmented nanoflagellate (**E**), and ciliate (**F**) abundance during the study period in Taiwan Strait.

**Figure 4.** Mean values were determined from the transects A–D (-) and transects E–I () for viral, bacterial, *Synechococcus* spp., HNF, PNF, and ciliate abundances.

**Figure 5.** Relationship between temperature and log(bacterial abundance) (**A**) and between viral and bacterial abundance (**B**) during the study period. The solid line is the regression line for all data (-); (-): the data for transects A–D; (): the data for transects E–I.

**Figure 6.** Relationship between *Synechococcus* spp. (**A**), PNF abundance (**B**), and Chl *a* concentration during the study period. The solid line is the regression line for all data.

Based on analysis of all the data, environmental conditions explained 43.7% of variance in the data for microbial community composition in the Taiwan Strait. We observed that HNFs, PNFs, ciliates, and *Synechococcus* spp. showed similar patterns of spatial fluctuations (Figure 7A). On the other hand, temperature and salinity were associated with bacteria (Figure 7A). This result is supported by the correlation analysis between both of them (Figure 5A). Interestingly, a different pattern was observed between the northern and the southern parts of the Taiwan Strait (Figure 7B,C). In the northern part of the Taiwan Strait, we noticed the marked impact of the bacteria characterized by temperature, salinity, and viruses (Figure 7B). However, in the southern part of the Taiwan Strait, the impact parameters appeared in a different way, as bacteria variations were associated with an abundance of viruses and HNFs, highlighting the importance of top-down control as a limiting factor on bacteria (Figure 7C).

**Figure 7.** Redundancy analysis (RDA)-based microbial community composition (bac = bacteria, vir = viruses, HNF = heterotrophic nanoflagellate, PNF = pigmented nanoflagellate, cil = ciliate, and Syn = *Synechococcus* spp.), and environmental conditions (chl a = chlorophyll a, T = temperature, S = salinity) as the explanatory variables, for all data (**A**), northern part (**B**), and southern part (**C**) of the Taiwan Strait. In (A), samplings from northern part and southern part of the Taiwan Strait are labeled using magenta and black, respectively. The percentages of the community composition variation explained by the first and second RDA are labeled on the two axes of each RDA plot. The explanatory environmental conditions and microbial compositions are presented as blue arrows and red lines on RDA plots, respectively.

#### **4. Discussion**

This study explored spatial variations in microbial communities across the Taiwan Strait during the winter months (January and February) of 2018. We used temperature and salinity in our study period to distinguish two hydrographic regions in the Taiwan Strait. Transects A–D and transects E–I were assumed in this study to be representative of the northern and southern parts of the strait, respectively. To the best of our knowledge, there are no previous studies of aquatic–microbial community dynamics in the Taiwan Strait, and so this study represents the first to investigate aquatic microbial ecology in this region. We found viral and bacterial concentrations to be higher at higher temperatures and in oligotrophic oceans, which have higher salinities in the southern part of the Taiwan Strait. *Synechococcus* spp., nanoflagellates, and ciliates had different spatial variations, with higher abundances of both found in the northern part (transects A–D) of the Taiwan Strait containing cold China coastal waters.

#### *4.1. Spatial Variations in Bacterial and Viral Abundance*

In marine environments, substrate availability and water temperature are important factors regulating bacterial production [9,30–33]. Some previous studies have suggested that, over the course of a year, temperature is the dominant factor affecting bacterial growth in colder seasons. Other factors, such as substrate supply, may be more important in warmer seasons [9,34]. In a study of a Delaware estuary, Hock and Kirchman [34] showed that its temperature varied over an annual cycle and that bacterial growth rate was directly related to temperatures below, but not above, 12 ◦C. In the present study, we found a clear relationship between temperature and bacterial abundance (ANOVA, *p* < 0.05) (Figure 5A), but not between Chl *a* and bacteria concentrations.

Viruses control bacterial mortality and may help to maintain bacterial community diversity [35]. In most aquatic environments, viral abundance was positively correlated with bacteria [27,35,36]. Interestingly, we found a closer linear relationship between viruses and bacteria in oligotrophic oceans in the southern part than those in the northern part of the Taiwan Strait (Figure 5B). One previous study of two lakes in the French Massif Central found that the relationship between viral and bacterial abundances to be weaker in the more productive lake, suggesting that there is an increase in relative abundance of non-bacteriophage viruses, such as cyanophages, in more productive environments [37]. In the present study, *Synechococcus* spp. abundance was significantly higher in the northern part of the study area (transects A–D) (Figure 4). This might explain the weaker linear relationship between the abundance of viruses and bacteria in the northern part of the Taiwan Strait (Figure 5B). Aside from this, differences in these relationships may be related to the burst size and infection rate of host cells [36]. Thus, future studies may want to include viral lysis and viral production when assessing the ecological role of phages and their potential for controlling bacterial dynamics.

#### *4.2. E*ff*ects of Bottom-Up and Top-Down Controls on HNFs*

According to one study, HNF abundance should be higher in warm environments than in colder ones, owing to consistent differences in the food web structure along the latitudinal gradient [38]. Some models of the predator–prey relationship have predicted that an increase in trophy of the freshwater environment should be reflected in both higher prey and higher predator abundances [39]. For example, a positive relationship has been found between bacterial and HNF abundances in marine systems [12,40–42]. Therefore, we also expected to find a relationship between the abundances of bacteria and HNFs in our study, but after considering all the data we collected, we found no consistent positive relationship between bacterial and HNF abundances (Figure 8), which is a correlation reported previously by many studies [12,41–43]. Šoli´c et al. [44] also reported a very low association between bacteria and HNF concentrations during the colder part of the year, suggesting that HNFs did not graze as much on bacteria during the colder months in the middle of the Adriatic Sea. In addition, top-down control on HNFs weakens the relationship between the abundances of bacteria and HNFs [15,16]. Therefore, a negative relationship was found between bacteria and HNFs in the northern part of the

Taiwan Strait (r = −0.57) (Figure 8), showing a top-down control on bacteria by HNFs in this study area. The significant and positive relationship between water temperature and bacterial abundance (r = 0.65) suggests a stronger bottom-up control of bacteria in the northern part of the Taiwan Strait (Figure 5A). For another possibility of both bacteria and HNFs showing a negative relationship, the high correlation between HNFs and ciliates (strong top-down control of HNFs) could explain their relationship through the "trophic cascade effect." This predation effect is known as the trophic cascade and is based on predation limitation at several trophic levels [45].

**Figure 8.** Simultaneous observations of bacterial and HNF abundances plotted following Gasol's model [16]. MAA is the maximum attainable abundance line, and MRA is the mean realized abundance line. (-): the data for transects A–D; (): the data for transects E–I.

We looked to interpret our data using Gasol's [16] theoretical framework to better understand the HNFs' limitation in our study (Figure 8). According to Gasol's model, points located below the MRA line (mean realized abundance) suggest top-down control of HNFs and points above it suggest low top-down control of HNFs. Points close to the MAA (maximal attainable abundance) line indicate clear bottom-up control of HNFs. The points were all below the MRA line in the southern part of Taiwan Strait in this study, denoting that HNFs were controlled more by predation (top-down) than available resources (bottom-up) (Figure 8). However, we did not observe a higher abundance of ciliates, top-down predators of HNFs, in the southern part of the Taiwan Strait (Figures 3D and 4). However, viral lysis is a major source of bacterial loss [46,47], and its effect on bacterial mortality has been compared with the effect of nanoflagellate grazing in many oceanic systems [47–49].

This study found a closer linear relationship between viruses and bacteria in the southern part than those in the northern part of the Taiwan Strait (Figure 5B), thus suggesting that a significant fraction of bacterial carbon and energy is probably not transferred to higher trophic levels, but instead is cycled in a bacteria–virus–dissolved organic matter (DOM) loop. This "viral loop" [50] might serve as an important control mechanism for bacterial production in the southern part of the Taiwan Strait. Therefore, future studies may want to take viral lysis into account when studying mechanisms controlling bacterial abundance. Aside from this, mixotrophic nanoflagellates (MNFs) and HNFs have been described as important grazers on bacteria; in particular, the relevance of MNFs as grazers in oligotrophic environments was recognized [8]. Furthermore, a group of grazers (MNFs) may be related to bacteria, but not to HNFs, which probably affects the bacteria–HNF relationship.

We applied our field data of bacteria and HNF abundances in the northern part of the Taiwan Strait to the model by Gasol [16] (Figure 8). Because most points lay between the MRA and MAA lines in this part of the strait in this study, top-down control on HNFs should be low and bottom-up control should be high (Figure 8). This part of the Taiwan Strait is relatively high in nutrients. Thus, our results were not in agreement with Gasol regarding dominant bottom-up control of HNFs in oligotrophic sites [16]. At the stations in the northern part of the strait, bacterial abundance was relatively low, because low temperature reduced bacterial growth. In this situation, HNFs could not satisfy their carbon demand through only grazing on bacteria in the northern part of the strait.

The increase in HNF abundance in the northern part of Taiwan Strait may have to do with the availability of other types of prey, such as *Synechococcus* spp. or picoeukaryotes, as well as low top-down control by predators such as ciliates. Some studies have found variations in the abundance of ciliates to be important in explaining the variations in abundance of nanoflagellates. Weisse [51] and Nakano et al. [52] reported a negative relationship between abundances of the two. According to earlier studies, predation control on HNFs should increase with productivity [12,16]. However, in the current study, we found a positive relationship between ciliates and nanoflagellates (ANOVA, *p* < 0.05). As for the ciliate community, ciliates of <30 μm in equivalent spherical diameter (ESD), including *Strombidium* spp., *Strobilidium* spp., and *Tontonia* spp., were the most abundant group (65–90%) of the ciliate community in the northern part of the strait in this study. These small ciliates probably control the production of picoplankton [53]. This positive relationship between ciliates and nanoflagellates may suggest that the same factors control the abundance of both communities, since they basically have the same resource requirements (picophytoplankton) and potential predators (zooplankton) [21]. Thus, the high abundance of HNFs obtained in relation to the low abundance of bacteria we found in the northern part of the strait in this study was probably due to the availability of other prey sources. HNFs are known to have other food resources (e.g., picophytoplankton) besides bacteria [54]. However, this does not change the fact that we found that HNFs were primarily resource limited in the northern part of the study area. The spatial variations in *Synechococcus* spp. abundance in the present study showed higher values in the northern part (transects A–D) (Figure 4). Thus, we believe that the major flux in energy in the northern part of the strait during winter was nanofalelltes grazing on picophytoplankton.

#### **5. Conclusions**

We found higher abundances of bacteria in waters with warmer temperatures and lower salinities, while nanoflagellate and ciliate abundances were found to be higher in the northern part of the Taiwan Strait due to cold China coastal waters. Our findings characterize two different ecosystems in the Taiwan Strait. We suspect that a "viral loop" plays a major role in controlling bacterial production in the southern part of the strait and that picophytoplankton plays a major role in the flux of energy in the northern part of the strait during winter, though further studies are needed to investigate these possibilities.

**Author Contributions:** Conceptualization: G.-C.G. and A.-Y.T.; methodology: A.-Y.T.; validation: A.-Y.T.; formal analysis: G.-C.G., C.-h.H., P.-C.H., and A.-Y.T.; investigation: H.-M.Y. and Y.-K.C.; resources: G.-C.G. and A.-Y.T.; data curation: G.-C.G. and A.-Y.T.; writing–original draft preparation: G.-C.G. and A.-Y.T.; writing–review and editing: G.-C.G. and A.-Y.T.; supervision: G.-C.G.; funding acquisition: G.-C.G. and A.-Y.T.

**Funding:** This study was supported by a grant (NSC 106-2611-M-019-013) from the Ministry of Science and Technology, Republic of China.

**Acknowledgments:** We appreciate the language editing and helpful comments from James Steed on this manuscript.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
