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

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.


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][7][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][17][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][27][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.

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 Diversity 2020, 12, 38 3 of 16 samples were analyzed from seven inshore sites of these lagoons. Vistonis lagoon has a surface area of approximately 45 km 2 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 km 2 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 2 and 0.8 km 2 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 km 2 and has a mean depth of approximately 1 m [35]. Lafra, Lafrouda and Palaia Koiti lagoons are connected with Vistonikos Gulf through tidal canals.

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 km 2 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 km 2 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 2 and 0.8 km 2 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 km 2 and has a mean depth of approximately 1 m [35]. Lafra, Lafrouda and Palaia Koiti lagoons are connected with Vistonikos Gulf through tidal canals.   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.

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 µm 3 was converted into carbon content using the formulae C = 0.288 V 0.811 for diatoms and C = 0.216 V 0.939 for other phytoplankton taxa (C is carbon content in pg, V is cell volume in µm 3 ) [39]. For individuals < 180 µm 3 , conversion factors of 0.108 pg C µm −3 for diatoms and 0.157 pg C µm −3 for other phytoplankton taxa were used [40].

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

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 = α1X 2 + a2X + a3, 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 R 2 . 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 R 2 was used to examine the strength of the relationship between the difference in conductivity and pair-wise bSOR.

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.

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 −1 (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 −1 (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 −1 (mean annual salinity > 28 psu). Chlorophytes and cyanobacteria appeared to be more diverse in Vistonis lagoon with mean annual conductivity ≤ 30,000 µS cm −1 (mean annual salinity < 18.5 psu). 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  (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 −1 ) at a conductivity of 16,090 µS cm −1 (9.7 psu). The minimum biomass (4.4 mg L −1 ) in Vis1 was measured in January 2019 dominated by the diatom Cyclostephanos sp. at a conductivity of 2620 µS cm −1 (2.6 psu). Similarly, high variability in phytoplankton biomass was recorded in Lafrouda lagoon reaching a maximum biomass of 127.8 mg L −1 in April 2019, concurrent with a Prorocentrum minimum bloom (301 cell mL −1 ) at conductivity 46,400 µS cm −1 (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 −1 (8.6 psu).    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).

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.  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).

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.

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 ≤ 30,000 μS cm −1 (mean annual salinity < 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 −1 (mean annual salinity > 28 psu). This separation was supported by one-way

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 ≤ 30,000 µS cm −1 (mean annual salinity < 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 −1 (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. 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, psame < 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.

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 2 of 0.4) was detected between conductivity and phytoplankton abundance while the weakest (R 2 = 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).

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.

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 −1 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 −1 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 km 2 . 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 −1 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].

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][67][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].

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.

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