Phytoplankton Dynamics and Water Quality in the Venice Lagoon

: We analyzed the phytoplankton abundance and community structure monthly over a 20-year period (1998–2017) at ﬁve stations in the Venice lagoon (VL), one of the sites belonging to the Long-Term Ecological Research network of Italy (LTER-Italy). We focused on phytoplankton seasonal patterns, inter-annual variability and long-term trends in relation to water quality. Diatoms numerically dominated (ca. 60% on average), followed by nanoﬂagellates (37%), while coccolithophorids and dinoﬂagellates contributed less than 2%. We observed distinct seasonal and inter-annual changes in the abundance and ﬂoristic composition of the phytoplankton groups, whilst no clear long-term trend was statistically signiﬁcant. We also assessed the water quality changes, applying to our dataset the multimetric phytoplankton index (MPI), recently ofﬁcially adopted by Italy to accomplish the water framework directive (WFD) requirements. The index evidenced a temporal improvement of the water quality from “moderate” to “good” and allowed us to conﬁrm its reliability to address the changes in the water quality, not only spatially—as previously known—but also for following the yearly time trends. Overall, our results highlight the importance of long-term observations, for understanding the variability in the phytoplankton communities of the lagoon as well as the relevance of their use to test and apply synthetic descriptors of water quality, in compliance with the environmental directives.


Introduction
An intrinsically high seasonal and interannual variability characterize the phytoplankton communities in coastal transitional environments [1,2], such as the Venice Lagoon (VL). The availability of multidecadal observations on phytoplankton in these ecosystems is crucial to depict seasonal patterns of abundance and biodiversity, interannual variability, and temporal trends, as well as to inform environmental policy, such as the European Water Framework Directive (WFD 2000/60/EC, [3]). Actually, phytoplankton is the only planktonic element included as a water quality indicator in the WFD, which enforced an integrated and coordinated framework for the monitoring of ecological status and the management of transitional waters, based on key biological elements [3]. In particular, some phytoplankton-related variables, such as taxonomic composition, abundance, biomass, frequency and intensity of blooms, and cell size spectra have been proved relevant for the definition and classification of the water quality of coastal transitional environments, in order to implement the directive [4][5][6][7]. A reliable tool for the classification of transitional waters, the multimetric phytoplankton index (MPI), has been formulated by Facca and colleagues [8]: the index, based on phytoplankton biodiversity, chlorophyll-a and bloom frequency, was specifically developed in the VL and successfully applied to compare specific areas in the lagoon. The Italian Ministry of Ecological Transition (formerly Ministry of the Environment) proposed the index to the Italian local authorities as an effective tool for adhering to WFD requirements. In recent years, the MPI has also been successfully applied in some Mediterranean lagoons in Sardinia [9], and it is currently used to assess the status of the Italian, Greek and Croatian transitional waters [10,11].
The VL is the largest Italian lagoon and is one of the research sites belonging to the Italian Long-Term Ecological research network LTER-Italy (www.lteritalia.it (accessed on 4 October 2021)). Data on the phytoplankton communities in the lagoon have been available sporadically since the 1970s, but they have only been gathered regularly and with comparable methods of investigation since the late 1990s [12][13][14]. The water environment of the VL is turbulent, nutrient-enriched, and open to other connected systems, i.e., the land, the sea and the sediments: these broad ecological features are those mainly affecting the prevalent phytoplankton community composition, its seasonal cycles and its morphological traits [13,15]. Diatoms dominate in the VL, being mainly responsible for the blooms and comprising most of the taxa representative of the complex riverine, coastal, benthic and pelagic influences. The prevalent phytoplankton annual cycle has been described as unimodal, with maxima typically attained in summer, mainly driven by the seasonal fluctuations of temperature and irradiance [13], a pattern that seems to characterize most temperate enclosed coastal ecosystems, with shallow depths and permanently high nutrient concentrations [1,16].
In the last two decades, a general improvement of the water quality of the VL has been documented [13,14,17]. In particular, water transparency and oxygen have increased, while nutrients (ammonia and nitrate) and chlorophyll-a have decreased. These changes seem to be related to several factors, such as the application of more stringent environmental regulations concerning farm wastes in the watershed, declining river discharge, and a reduction in clam harvesting in the lagoon. At the same time, an upswing of macrophytes has been observed [14,18], evidencing the complex interrelations between the different primary producers, which can exploit the LV habitat heterogeneity, and the different trophic and water quality conditions.
In light of these environmental changes, in this paper, we analyze 20 years (1998-2017) of monthly phytoplankton data, gathered in the northern and central basin of the VL, with the main aim to identify the overall changes that might have occurred in the community through time. We assess the overall features and the main changes of the phytoplankton community and of the MPI, further developing its applicability and reliability to track the water quality changes through time. In particular, we focus on (i) identifying the prevalent seasonal cycle and the most characteristic taxa or taxa assemblages in each season, (ii) assessing the interannual variability and the long-term trends, and (iii) comparing the water quality index MPI in different years.

Study Site
The VL ( Figure 1) is a large (550 km 2 ) Mediterranean coastal transitional environment located in Italy, in the Northeastern part of the Adriatic Sea. It has an average depth of 1 m, and it is morphologically characterized by the presence of canals (5-10 m deep) and of a multiplicity of aquatic (e.g., salt marshes, shoals, and mud flats) and terrestrial (e.g., islands, coastal strips) habitats. The VL is microtidal, poly-and euhaline, and, in terms of confinement, it has both choked and restricted waters [19]. Two long barrier islands separate the lagoon from the Adriatic Sea to the east, allowing for water exchange through three large inlets (Lido, Malamocco and Chioggia), each corresponding to a specific basin inside the lagoon (north, central and south). Water exchange in VL is governed by the Adriatic tides, whose average amplitude ranges between 20 cm at neap tides and 100 cm at spring tides. The average residence time of VL waters is around 10 days, varying from a few days close to the inlets to 30 days for the internal areas [20,21]. The yearly freshwater discharge in the VL, coming from twelve main tributaries, is about 35 m 3 s −1 on average, with seasonal peaks in spring and autumn. The VL is subjected to high levels of various anthropogenic pressures, due to the presence of urban areas, tourism, fishing, The northern one is characterized by the presence of tributaries draining the waters from an intensively cultivated hinterland; the central one is subject to strong nutrient enrichment from the city of Venice and from the major islands, and it receives treated urban waste waters from the urban centre of Mestre and the discharge of wastewaters and cooling waters from the nearby industrial areas of Fusina and Marghera (Figure 1). The northern one is characterized by the presence of tributaries draining the waters from an intensively cultivated hinterland; the central one is subject to strong nutrient enrichment from the city of Venice and from the major islands, and it receives treated urban waste waters from the urban centre of Mestre and the discharge of wastewaters and cooling waters from the nearby industrial areas of Fusina and Marghera ( Figure 1).

Sampling Strategy, Hydrochemical Parameters and Laboratory Methods
Monthly samplings were carried out from 1998 to 2017 at the near-surface layers of each station, at neap tide and prevalently in the morning, maintaining the duration of the whole operations as short as possible (around 4 h in total). For each parameter, we collected and analyzed one sample at each station. The parameters measured and the relative methods are reported in Table 1. Analytical quality was assessed via participation in the Quality Assurance of Information for Marine Environmental Monitoring in Europe international laboratory proficiency-testing programme (QUASIMEME; http://www.quasimeme.org (accessed on 4 October 2021)).

Sampling Strategy, Hydrochemical Parameters and Laboratory Methods
Monthly samplings were carried out from 1998 to 2017 at the near-surface layers of each station, at neap tide and prevalently in the morning, maintaining the duration of the whole operations as short as possible (around 4 h in total). For each parameter, we collected and analyzed one sample at each station. The parameters measured and the relative methods are reported in Table 1. Analytical quality was assessed via participation in the Quality Assurance of Information for Marine Environmental Monitoring in Europe international laboratory proficiency-testing programme (QUASIMEME; http://www.quasimeme.org (accessed on 4 October 2021)). Chlorophyll a Perkin Elmer LS5B spectrofluorometer, [24]. From 2010 Turner Trilogy Laboratory Fluorometer. Phytoplankton was fixed in neutralized formalin with hexamethylentetramine, recognized and counted with an invertoscope after samples had settled in 2-25-mL chambers, for 12-24 h. Cell counting was carried out at 400× magnification along transects, whose number varied depending on cell abundance, to count a minimum of 200 cells (but often more than 500) in each sample [25][26][27]. Over the study period, two different operators analyzed the phytoplankton samples using the same technique and with accurate intercalibration. Taxa composition was mainly defined in accordance with [28,29]. All the species were reviewed and conformed to any synonyms through the World Register of Marine Species WoRMS (https://www.marinespecies.org (accessed on 4 October 2021)) and the World's Algae databases Algaebase, (http://www.algaebase.org/ (accessed on 4 October 2021)). The analysis was confined only to those forms that were detectable by light microscopy, i.e., up to 3 µm, thus not considering the picophytoplankton fraction. For a detailed analysis of interannual variations and trends, we focused only on the four major groups of phytoplankton: diatoms, nanoflagellates, dinoflagellates and coccolithophorids. For the dinoflagellates, no distinction has been made between autotrophic and heterotrophic species, which were therefore all included in the countings. The nanoflagellate group contained all the undetermined organisms, whose size was always less than 10 µm, mostly around 3-4 µm, mainly consisting of cryptophyceans, chrysophyceans, prymnesiophyceans (except coccolithophorids), chlorodendrophyceans and undetermined flagellates.

Phytoplankton Data Processing and Statistical Analyses
Statistical analyses were carried out on a matrix composed of 422 taxa and 1202 samples, after log-transformation of biological data [30] to account for non-normal data distributions [31].
Since the main focus of this study was the individuation of the overall features of the central and northern basins of the VL, we did not address any comparison among the single stations, which we pooled together to obtain an overall view of the area. Indeed, previous studies [13] evidenced that all stations shared the most abundant phytoplankton species and that taxonomic composition did not exhibit significant spatial differences.
We performed different types of analysis in order to: (i) identify the prevalent seasonal cycle and the most characteristic taxa or taxa assemblages in each season, (ii) assess the interannual variability and the long-term trends, and (iii) compare the water quality index MPI in the different years.

Annual and Seasonal Variations
The average phytoplankton seasonal cycle was assessed calculating the averages abundance of each month in the 20 years.
In order to identify the most characteristic species or species assemblages of each season, the Indicator Value (IndVal [32]), which combines the relative abundance of a species with its relative frequency of occurrence in a given seasonal period, was calculated on the seasonal average abundance. The IndVal is among the indicators most suitable to describe the average seasonal pattern of the plankton community, when a pluriannual set of data is available [13,[33][34][35]. The higher the mean abundance and the relative frequency of occurrence of one species in a season, the higher the IndVal. These two terms are multiplied and then scaled to 100 to express the indicator value of a species with respect to the cluster in terms of percentages. The significance of the indicator values for each species was tested using a Monte Carlo permutation test (999 random permutations).
We adopted a conventional division of the seasons, which we defined as follows: January-March, winter; April-June, spring; July-September, summer; October-December, autumn. This division was used to calculate the IndVal and to individuate the taxa most typical in each season. In presenting the results, we evidenced whenever a taxa characterized different part of a season, with an early or late presence.

Interannual Variability and Pluriannual Trends
For providing an overview of phytoplankton interannual variability, we used a 12-month moving average. The pluriannual trends were analyzed by applying to the whole dataset a seasonal Kendall-tau nonparametric analysis following [36]. For the MPI, we performed the principal component analysis (PCA, R-mode) on the annual MPI values and on the yearly average values of each hydrochemical variable in order to identify their relationships and to display groups of observations in a 2-dimensional space. To assess the evolution of MPI, estimated every year, a simple Mann-Kendall test was used according to [37].

Multiparametric Phytoplankton Index
The MPI is an annual assessment of the water quality, which combines a diversity index (Menchinick's Diversity, [38]) and a dominance index (Hulburt's Dominance) with bloom frequency and chlorophyll-a concentration. The rationale behind the index and the detailed procedure for its application are reported in [8]; here, we give just the main information regarding the components of the index: (i) Hulburt's index, the first metric, is a species dominance index, calculated as δ 2 = 100 (n 1 − n 2 )/N t , where n 1 is the abundance of the dominant species, n 2 is abundance of the second most abundant species, and N t is the total abundance. In order to relate Hulburt's index to water quality, the "100-δ 2 ' value is used. (ii) Bloom frequency, the second metric, is the number of times in each year in which the abundance contribution of a single species exceeds the 50% of the total. The numerical value obtained is a percentage, and it is used as "100-Bloom Frequency". (iii) Menhinick's index, the third metric, is a diversity index calculated as D = S/ √ N, where S is the number of species and N is total abundance. To reduce the error caused by deletion of multiple indeterminate taxa, a correction factor was introduced, estimated as the ratio of the sum of the determinate taxa to the original total abundance (determinate + indeterminate taxa). (iv) Chlorophyll-a concentrations, the fourth metric, were obtained after log 10 -transformation, removal of outliers (values outside the range |average ± 2.5 × std. dev.|) and calculation of antilogs.
Following the indications provided by the WFD, the range of the index was divided into five classes of ecological status, which were already adopted in [8], spanning from "bad" (eutrophic conditions with very low biodiversity, opportunistic phytoplankton species dominating the community, leading to very high chlorophyll-a concentrations) to high (communities characterized by high biodiversity and low chlorophyll-a): Bad: 0-0.2, Poor 0.21-0.4, Moderate 0.41-0.6, Good 0.61-0.8 and High 0.81-1.

Results
The range of the values of the physical and chemical variables in the 20 years (Table 2), which have been thoroughly analyzed in [17], clearly testifies to the marked environmental heterogeneity of the features of the lagoon waters, as indicated by the variability of temperature (0.6 • C to 33.  Phytoplankton abundance was overall mainly made up of diatoms (60% on the average), followed by nanoflagellates (37%), dinoflagellates (2%) and coccolithophorids (0.2%). The relevance of the other taxa was negligible (far below 0.1%).

Seasonal Pattern
The average annual pattern of the total phytoplankton abundance (Figure 2a They were characterized by minor peaks in the late-winter period, by yearly maxima attained in summer, and then by a gradual decline in autumn, until reaching the lowest abundance typical of the winter conditions. As regards the species composition (Table 3), resuspended species (Gyrosigma fasciola and Melosira spp.) characterized the beginning of the winter, whilst Skeletonema marinoi dominated the first phase of the growing season in late-winter, followed by resuspended diatoms (Navicula spp., Halamphora hyalina, Licmophora gracilis) and small (<10 µm) Cyclotella spp. in spring. The highest summer diatom peaks were mainly due to Cylindrotheca closterium, small (<10 µm) Thalassiosira spp. (undetermined species and T. pseudonana), small (<15 µm) pennate diatoms (Nitzschia frustulum and Phaeodactylum tricornutum) and small (<10 µm) forms of Chaetoceros spp., together with Pseudonitzschia spp. (P. galaxiae and P. multistriata), Leptocylindrus minimus and Skeletonema tropicum. Resuspended and benthic diatoms (Placoneis elginensis, Navicula spp. and Cocconeis spp.) typified the low autumn abundance. Some samples were also characterized by the presence, with low numbers, of rare species of allochthonous origin ( Table 3): marine species such as Proboscia alata, Cerataulina pelagica, Leptocylindrus minimus and Pseudo-nitzschia spp. in summer and freshwater species in autumn (Craticula cuspidata and Sellaphora pupula) and winter (Navicula tripunctata and Amphora lineolata).
together with Pseudonitzschia spp. (P. galaxiae and P. multistriata), Leptocylindrus minimus and Skeletonema tropicum. Resuspended and benthic diatoms (Placoneis elginensis, Navicula spp. and Cocconeis spp.) typified the low autumn abundance. Some samples were also characterized by the presence, with low numbers, of rare species of allochthonous origin ( Table 3): marine species such as Proboscia alata, Cerataulina pelagica, Leptocylindrus minimus and Pseudo-nitzschia spp. in summer and freshwater species in autumn (Craticula cuspidata and Sellaphora pupula) and winter (Navicula tripunctata and Amphora lineolata). Nanoflagellates (Figure 2c) were second to diatoms as abundance throughout the time-series, averaging 1,333,818 cells L −1 (SD 32,091,079 cells L −1 ), and they were generally most abundant between April and September.
Dinoflagellates (Figure 2d) represented less than 2% of the total abundance, with an average of 72,496 cells L −1 (SD 317,697 cells L −1 ). Overall, they were generally most abundant in spring (April) and summer (August), although rarely exceeding 1 × 10 6 cells L −1 , while during autumn and winter the abundances were the lowest, around 1-2 × 10 4 cells L −1 . The peaks were due to Prorocentrum spp. (P. minimum, P. gracile and P. triestinum) in spring and naked species (small (<10 µm) undetermined species, Gymnodinium spp.), Blixaea quinquecornis, Protoperidinium bipes and Prorocentrum rhathymum in summer (Table 3). Table 3. List of main phytoplankton taxa characterized by the highest and significant IndVal for each season in the period 1998-2017. For each taxa, we also indicate the divisions to which they belong. IndVal values indicated with ** are significant at p < 0.01, and those with *** are significant at p < 0.001. Taxa are listed in order of decreasing importance in each season. Coccolithophorids, represented almost exclusively by Emiliania huxleyi (71% on average), were recorded in low numbers throughout the study (on average 9405 L −1 , SD 51,957 cells L −1 ), with an average contribution to the total abundance of 0.2%. Emiliania huxleyi was responsible for both the absolute peak, recorded in September, and of the highest relative abundance, attained in winter (Figure 2e and Table 3). The summer and autumn samples were characterized also by the presence of Michaelsarsia adriatica and Calciosolenia murrayi, respectively (Table 3).

Division
Other taxa reported in Table 3 are mostly rare and of allochthonous origin, their presence being a result of freshwater inputs: the euglenophyceans Eutreptia globulifera and the chrysophycean Dinobryon coalescens marked the winter period; Ankistrodesmus spp. From early spring to late summer, we occasionally observed some species that may potentially trigger Harmful Algal Blooms (HAB), such as the dinoflagellates Alexandrium minutum, Dinophysis caudata, D. fortii, D. sacculus, Gonyaulax polygramma, G. spinifera and Prorocentrum rhathymum and the diatoms Halamphora coffeiformis, Pseudo-nitzschia delicatissima, P. pseudodelicatissima complex, P. fraudulenta, P. galaxiae and P. pungens. These species were recorded on average in 4% of the samples (varying between 0 and 15%) and with values rarely exceeding 1 × 10 6 cells L −1 , contributing to a total abundance on average of less than 4%. Their highest values of abundance were found in summer, in accordance with the general pattern of phytoplankton growing season. Detrimental effects of these species were never recorded in the VL.

Interannual Variability and Long-Term Trends
The yearly total phytoplankton abundance varied from 1,943,797 to 10,119,165 cell L −1, and the changes were mainly determined by the yearly diatom variations (Figure 3). The diatoms peaks were always due to blooms, made by a few species: in 2001, Thalassiosira spp. and Cyclotella spp. in spring and Cylindrotheca closterium, Skeletonema tropicum and Nitzschia frustulum in summer. The year 2003 was characterized by summer blooms of Cyclotella spp., Thalassiosira spp. and Nitzschia frustulum; the year 2006 by an intense summer bloom of Thalassiosira spp. In 2009 and 2013, late winter blooms of Skeletonema marinoi occurred, which were followed in 2013 by summer blooms of a mixed community made by Chaetoceros spp., Thalassiosira spp. and Nitzschia frustulum. In the years with the lowest diatom abundance, this same pool of species was always responsible for blooms, even though they were far less intense.
The variation of the other phytoplankton groups was less marked ( For what concerns the long-term trend, the seasonal Kendall test evidenced negative slopes for the total phytoplankton abundance, for diatoms, nanoflagellates and dinoflagellates, and a positive one for coccolithophorids. All the trends, however, were not statistically significant. On the contrary, the MPI and each of its single four metrics showed statistically significant increasing trends: metric 1 (Mann-Kendall trend test: τ = 0.58; p < 0.05) and metric 3 (τ = 0.45; p < 0.05) highlighted an increase of biodiversity, metrics 2 (τ = 0.33; p < 0.05) a lower frequency of bloom and metric 4 a decreasing trend of phytoplankton chlorophyll-a (τ = 0.67; p < 0.05).
Water 2021, 13, x FOR PEER REVIEW 10 of 18 community made by Chaetoceros spp., Thalassiosira spp. and Nitzschia frustulum. In the years with the lowest diatom abundance, this same pool of species was always responsible for blooms, even though they were far less intense. For what concerns the long-term trend, the seasonal Kendall test evidenced negative slopes for the total phytoplankton abundance, for diatoms, nanoflagellates and dinoflagellates, and a positive one for coccolithophorids. All the trends, however, were not statistically significant. On the contrary, the MPI and each of its single four metrics showed statistically significant increasing trends: metric 1 (Mann-Kendall trend test: τ = 0.58; p < 0.05) and metric 3 (τ = 0.45; p < 0.05) highlighted an increase of biodiversity, metrics 2 (τ = 0.33; p < 0.05) a lower frequency of bloom and metric 4 a decreasing trend of phytoplankton chlorophyll-a (τ = 0.67; p < 0.05).
Indeed, the annual values of the MPI (Figure 4) showed a prevalence (12 years) of moderate conditions, with the frequency of good water quality levels increasing in the most recent years (2014-2017), while poor conditions were detected only twice (2002)(2003) in the 20 years. To investigate the relations of MPI with hydrochemical parameters, we performed a PCA ( Figure 5). The bi-plot of the first two PCA components explained about 62% of the total variance. The squared cosine values identify the variables best linked to each of the two axes: the greater the squared cosine, the greater the link with the corresponding axis and, consequently, the influence of the variable along the axis. The first axis mainly explains the variability related, on one side, to salinity and, on the other, to nitrogen nutrients and N/P ratio. The second one explains the variability related to O2 percentage, water transparency, water temperature and MPI, which are inversely correlated with N-NH3 and N/P ratio. These variables allow evidencing in the PCA two clusters of consecutive years (Table 4 and Figure 5): cluster A and cluster B that group, respectively, the first 14 years (1998-2011) and the most recent ones (2012-2017). The division of these two groups occurs at the value 0.5 of the second component. The ordering of the years in relation to the first two main components appears to be fairly consistent with the MPI values ( Figure 4) that appear to be the highest in well-oxygenated, transparent and less nutrient-rich waters. Indeed, the comparison of the two periods ( Figure 5 and Table 4) To investigate the relations of MPI with hydrochemical parameters, we performed a PCA ( Figure 5). The bi-plot of the first two PCA components explained about 62% of the total variance. The squared cosine values identify the variables best linked to each of the two axes: the greater the squared cosine, the greater the link with the corresponding axis and, consequently, the influence of the variable along the axis. The first axis mainly explains the variability related, on one side, to salinity and, on the other, to nitrogen nutrients and N/P ratio. The second one explains the variability related to O 2 percentage, water transparency, water temperature and MPI, which are inversely correlated with N-NH 3 and N/P ratio. These variables allow evidencing in the PCA two clusters of consecutive years (Table 4 and Figure 5): cluster A and cluster B that group, respectively, the first 14 years (1998-2011) and the most recent ones (2012-2017). The division of these two groups occurs at the value 0.5 of the second component. The ordering of the years in relation to the first two main components appears to be fairly consistent with the MPI values ( Figure 4) that appear to be the highest in well-oxygenated, transparent and less nutrient-rich waters. Indeed, the comparison of the two periods ( Figure 5 and Table 4) indicated statistically significant differences for transparency, temperature and relative oxygen, which were highest in the second period. As regards nutrients, ammonium, nitrites and nitrates decreased significantly in recent years, resulting as well in a lower N/P ratio [17]. These changes in environmental variables were accompanied by an increase of the water quality assessed with MPI, which was significantly higher in the most recent years (2012)(2013)(2014)(2015)(2016)(2017). For what concern the bloom species, three of them decreased (Skeletonema spp., small Nitzschia spp., Cylindrotheca closterium), while two (small Chaetoceros spp. and Thalassiosira spp.) increased.
indicated statistically significant differences for transparency, temperature and relative oxygen, which were highest in the second period. As regards nutrients, ammonium, nitrites and nitrates decreased significantly in recent years, resulting as well in a lower N/P ratio [17]. These changes in environmental variables were accompanied by an increase of the water quality assessed with MPI, which was significantly higher in the most recent years (2012)(2013)(2014)(2015)(2016)(2017). For what concern the bloom species, three of them decreased (Skeletonema spp., small Nitzschia spp., Cylindrotheca closterium), while two (small Chaetoceros spp. and Thalassiosira spp.) increased.

Discussion
This paper could be considered as a "third episode of a series" addressing the changes of the water quality in the VL within the same 20-year period (1998-2017), from different points of view. Actually, two previous papers analyzed the VL water quality considering, respectively, (i) the abiotic parameters and the chlorophyll-a [17], and (ii) the relations between phytoplankton and macrophytes [14]. They both evidenced an overall improvement of the lagoon water quality, which seemed to enter a new phase in the most recent years. This phase is characterized, in particular, by a reduction in the trophic state [17], mainly related to lower inputs of inorganic dissolved nitrogen nutrients from the basin-which followed the Council Directive 91/676/EC introducing good agricultural practices concerning the prohibition of spreading nitrogen rich farm wastes-and by an increase in water transparency and dissolved oxygen, also related to the augmented biomass of both macroalgae and seagrasses, which started recolonizing the lagoon, boosted by the clam fishing regulation [14].
In this paper, following and grounded in the information coming from the former two, we focused only on the phytoplankton communities, examining their composition, seasonal patterns and trends, and addressing their use for the assessment of the changes in the water quality through time. This last issue, in particular, was tackled through the application of the multiparametric index MPI [8], based on phytoplankton biodiversity, chlorophyll-a and bloom frequency, which was specifically developed in the VL and, then, became an effective tool adopted by the Italian Ministry, local authorities, Greece, and Croatia to adhere to WFD requirements [11].
It is well known that long-term research is a unique tool for obtaining a consistent picture of the annual cycle of phytoplankton communities and for assessing trends and shifts [39][40][41][42][43]. Long-term series are crucial for revealing and interpreting the pattern of changes in phytoplankton biomass and composition, in particular in transitional environments, where an intrinsically high spatial and temporal heterogeneity characterizes phytoplankton, making its characterization and the consequent development of adequate indicators quite difficult [6,[44][45][46]. In this respect, long-term series also appear essential to provide a robust base for assessing water quality changes as required by the WFD, which explicitly mentions the composition, abundance, frequency and intensity of algal blooms and biomass as phytoplankton quality elements in coastal and transitional waters [3].
Our study confirmed some features of the LV phytoplankton communities, which were already previously evidenced by a study based on a 10-year period [13]. First of all, the prevalent unimodal annual pattern of phytoplankton abundance, with a typical summer peak, was observed also in the present study. This pattern characterizes temperate enclosed coastal ecosystems, with shallow depths and permanently high nutrient concentrations [16]: the observed decrease in nutrient availability in the VL did not seem to affect the phytoplankton unimodal cycle, so that the seasonal climate cycle still remains the most recognizable driver of phytoplankton. This seasonal pattern is mainly driven by the diatoms, which are so far the most prevalent phytoplankton group in the lagoon, as typical of permanently fertile and turbulent environments [2,47]). Actually, diatoms are-from one perspective-the main responsible for the annual blooms in the VL, which occur each year mainly in late winter/spring and summer; from another, they are inclusive of taxa that represent the different systems (rivers, sea and sediments) connected with the lagoon waters. Blooms, which are essential features of many transitional ecosystems, are prevalently associated with diatoms due to their capacity to grow rapidly in turbulent high-nutrient environments [48]. Diatoms also dominate the phytoplankton in many other Mediterranean lagoons [49][50][51][52][53], and some of the genera responsible for blooms in the VL (i.e., Skeletonema, Thalassiosira and Chaetoceros) are key species known to produce blooms in transitional waters globally [48].
The second most abundant phytoplankton group in the VL is a taxonomically heterogeneous one, composed of undetermined flagellates with a size lower than 10 µm, which constitute on average almost the 40% of the total abundance in the lagoon. Unde-termined nanoflagellates are also a typical and relevant component of other transitional environments [54][55][56][57][58]. In the VL, they contribute to the spring and summer peaks, and they represent a quite stable and important component throughout the years. Nanoflagellates are difficult or impossible to identify with light microscopy and in routine phytoplankton long-term monitoring they will very likely always be left undetermined, so that an important part of the phytoplankton biodiversity will remain hidden. Actually, recent studies on the VL protists through metabarcoding [59,60] revealed that nanoflagellates are represented by several distinct species, some of which (e.g., Picochlorum, Micromonas) were reported for the first in the lagoon in that study.
The average composition of the phytoplankton seasonal cycle in the VL (Table 3) is still quite comparable with that described for the previous ten-year period 1998-2007 [13]: diatom species imported from the benthic community (resuspended species, in particular Navicula spp., Amphora spp. and Cocconeis spp.) and coccolithophorids from the coastal sea (Calciosolenia murrayi, Emiliana huxleyi) typify the periods of lowest abundance (autumn and early winter). Planktonic diatoms (Skeletonema marinoi), together with Cyclotella spp. and some resuspended species (Amphora spp., Navicula spp., Lichmophora spp.), prevail in late winter and spring, together with a few dinoflagellates species (Prorocentrum spp.). In summer the species composition become characterized by a mixed assemblage of diatoms (Chaetoceros spp., Thalassiosira spp., Cylindrotheca Closterium and Nitzschia frustulum), which remain dominant, together with dinoflagellates (Gymnodinium spp., Protoperidinium bipes).
The trends of the total abundance and of each of the four main taxa appeared to decrease in the 20 years, although not in a statistically significant way, thus providing a weak signal that they are related to water quality improvement. A quite univocal sign came, instead, from the temporal variations of the MPI, which showed a statistically significant increasing trend, as a whole and in each single component. The value of the index shifted from the prevalent "moderate" quality class, which characterized the former years (until 2010), to the "good" class, which appeared quite stable since 2011. Blooms, chlorophyll-a, and biodiversity, when pooled together, seem to represent a most powerful tool to address the changes in water quality in the lagoon. The possibility to successfully apply the MPI has been tested, until now, only to compare spatially different areas of the same ecosystem or different lagoons [8,9], in particular due to its sensitivity to the water trophic conditions. In this work, we also provide evidence of the possibility to use it to compare the water quality changes over time, a necessary requisite to evaluate the effectiveness of restoration measures, in particular in relation to the WFD. It is worth noticing that the index appears to work effectively, even in the presence of many undetermined taxa, as is the case of the nanoflagellate component in the VL.
The MPI actually proved to be well correlated with the abiotic parameters in the VL, showing the highest value in well-oxygenated and less nutrient-rich waters. This correlation also allowed one, through the use of PCA, to separate two clusters of years (1998-2011; 2012-2017), the most recent ones characterized by improved water quality conditions. A statistically significant decrease in the frequency of the blooms was recorded comparing the two groups of years, even though there is not a univocal tendency among the taxa that typically cause them: actually, some, such as Cylindrotheca, Nitzschia and Skeletonema, decreased, while others (i.e., Thalassiosira and Chaetoceros) increased and, overall, the very same taxa were always mainly responsible for the blooms, whatever their intensity.
In conclusion, we wish to highlight that a general water quality improvement has taken place in the VL in the 20-year period 1998-2017, in particular for what concerns the trophic conditions and the recovery of the macrophytes. The single signals of changes in the phytoplankton communities were somewhat weak: as already evidenced in the past, the prevalent unimodal pattern of seasonal variation, the diatoms as dominant group and the wide interannual variability remained the typical unchanged features of the communities, which did not show a statistically significant decreasing trend of abundance over the years. The MPI, which pools together different features of the phytoplankton, appears instead to robustly confirm its reliability to address the changes in the water quality, not only spatially-as previously known-but also for following the yearly time trends. A multimetric index that considers as many fundamental attributes of phytoplankton as possible (e.g., biomass, community structure and frequency of blooms) is indeed considered more sensitive and robust than the use of single metrics in assessing the status of transitional water ecosystems [7,46,58] and its spatial and temporal variations. This is particularly relevant considering that the trajectories of coastal ecosystems during oligotrophication may be more complex than expected as different controlling factors may change at the same time [61,62].

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.