Phytoplankton Dynamics in the Middle Adriatic Estuary, with a Focus on the Potentially Toxic Genus Pseudo-nitzschia

The Krka River estuary is a karstic, permanently stratified estuary due to the strong freshwater inflow. It is a special environment in which to study the phytoplankton community, especially because this area is an important aquaculture site. Among other potentially toxic phytoplankton species, the diatom genus Pseudo-nitzschia occurs frequently and is a potential source of domoic acid (DA), causing shellfish toxicity and human intoxication. The main objective was to examine the dynamics of the phytoplankton community and, in particular, the genus Pseudo-nitzschia in the upper part of the Krka estuary, through monthly sampling over two years. The phytoplankton community was analysed using light microscopy and scanning electron microscopy to determine the diversity of Pseudo-nitzschia species and characterise the environmental parameters associated with a high abundance of Pseudo-nitzschia species. Seven Pseudo-nitzschia species were identified in the investigation, with higher frequencies and abundances in the less variable layer, at a 7 m depth. Blooms of Pseudo-nitzschia were noted in the late summer/early autumn, dominated by P. delicatissima/arenysensis. Winter assemblages were characterised by P. pseudodelicatissima/cuspidata, P. calliantha, and P. subfraudulenta, and were associated with domoic acid occurrence in shellfish.


Introduction
Estuaries form transitional zones between freshwater and marine environments and are characterised by large variability in physical, chemical, and biological properties under the dual influence of climate and anthropogenic changes [1]. The stratified estuaries could be considered as a unique model for studying the impacts of riverine and marine processes [2]. The Krka River estuary is a karstic estuary with a very low allochthonous source of organic matter [3]. Previous studies of the Krka estuary concerning the phytoplankton composition were based on episodic or seasonal sampling [3][4][5], which could not provide insight into the detailed variations of the phytoplankton community throughout the year. The estuary is one of the significant aquaculture areas along the Eastern Adriatic coast and has, since 2000, been included in the National Monitoring Programme of shellfish and seawater quality controls. Long-term research shows the occurrence of lipophilic toxins, produced in this area by dinoflagellates [6,7], as well as the occurrence of domoic acid (DA), produced by Pseudo-nitzschia species, which has been reported rarely and at low concentrations [8,9].
Pseudo-nitzschia is an important diatom genus, regularly occurring in the phytoplankton community of the Adriatic Sea and also frequently in the Krka River estuary, with an increasing trend in environment, DA production is very complex, according to the synergistic effect of the various ecological factors [16].
Considering the ability of Pseudo-nitzschia to produce domoic acid, which has been detected in this area, we were interested in investigating Pseudo-nitzschia composition and toxicity. The main goal was to examine the phytoplankton community comprehensively and, in particular, the genus Pseudonitzschia in the upper part of the Krka estuary, through monthly sampling over two years.
The specific objectives were: (i) to characterise the seasonal and vertical distribution of the phytoplankton community; (ii) to describe the Pseudo-nitzschia diversity; (iii) to determine the environmental drivers associated with the development of Pseudo-nitzschia assemblage; and (iv) to analyse domoic acid (DA) in shellfish.

Study Area and Sampling
The River Krka estuary is a 25 km long estuary, highly stratified, located at the central part of the Eastern Adriatic coast. Due to a low tidal oscillation, a sharp halocline exists throughout the whole year in the layer of 1 m to 5 m in depth. This halocline strongly depends on the River Krka's inflow. Previous studies confirmed that halocline plays an important role in biogeochemical processes [5] and represents the accumulation layer for freshwater and marine phytoplankton species [3]. In contrast to the lower part of the estuary, enriched with nutrients from natural and anthropogenic sources (sewage and domestic wastewaters) [17,18], the upper part of the estuary is characterised by low concentrations of nutrients and an extremely low input of terrigenous material from the Krka River, one of the most pristine European rivers [19]. Our sampling site, Strmica (43.776320° N, 15.848007° E, depth 28 m; Figure 1), is located in the upper part of the estuary and is an important shellfish farming area. Sampling was performed over two years, from August 2015 to October 2017. In general, samples were taken monthly, except in periods when domoic acid (DA) was found in shellfish (January-February 2016) and with high Pseudo-nitzschia spp. abundance (more than 10 6 cells L −1 ; September-October 2017). In these periods, samples were taken more frequently, on a weekly basis. During the sampling period, water temperature and salinity were measured using a YSI probes (YSI Professional Plus and YSI Pro 1030; YSI Incorporated, Yellow Springs, OH, USA) at five depths (surface, halocline, 2 m, 5 m, and 7 m).
Concurrently with temperature and salinity measurements, samples for phytoplankton analyses, chl a, and nutrient concentrations were collected with Nansen bottles at three depths (surface, halocline, and 7 m). Mussels (Mytilus galloprovincialis) were sampled monthly for toxin analyses, except during the period when DA was found in shellfish, when the frequency of sampling was weekly. Additional seawater samples were taken by a phytoplankton net (upper diameter 35 cm, mesh size 20 µm) in the 0-7 m layer, to determine the Pseudo-nitzschia species composition by SEM during the period with recorded shellfish toxicity and with a Nansen water sampler during the second year of the study (July 2016-October 2017). Vertical profiles of temperature and salinity were presented using Ocean Data View [20].  Concurrently with temperature and salinity measurements, samples for phytoplankton analyses, chl a, and nutrient concentrations were collected with Nansen bottles at three depths (surface, halocline, and 7 m). Mussels (Mytilus galloprovincialis) were sampled monthly for toxin analyses, except during the period when DA was found in shellfish, when the frequency of sampling was weekly. Additional seawater samples were taken by a phytoplankton net (upper diameter 35 cm, mesh size 3 of 22 20 µm) in the 0-7 m layer, to determine the Pseudo-nitzschia species composition by SEM during the period with recorded shellfish toxicity and with a Nansen water sampler during the second year of the study (July 2016-October 2017). Vertical profiles of temperature and salinity were presented using Ocean Data View [20].

Environmental Parameters and Shellfish Toxin Analyses
The chlorophyll a concentration was analysed according to the method of Strickland and Parsons [21]. Seawater samples (500 mL) were filtered through a glass microfiber filter (Whatman GF/F; pore size 0.7 µm) and frozen at −20 • C until the analyses. Chl a was extracted in 90% acetone for 2 h. The fluorescence was measured using a Trilogy fluorimeter (Turner Designs, Sunnyvale, CA, USA), calibrated with Chl a (Sigma, Taufkirchen, Germany) solutions prepared at six concentrations.
Water samples collected for nutrient analyses were stored in the freezer immediately upon sampling and transferred to the laboratory for further analyses. Domoic acid (DA) and epi-DA in shellfish samples were analysed by a High-performance Liquid Chromatography system coupled with Photodiode Array and Ultraviolet detector (HPLC-DAD-UV; Varian ProSTAR; Varian Inc., Walnut Creek, CA, USA) according to Quilliam et al. [23]. Approximately 100 g of soft shellfish tissue was homogenised; afterwards, a subsample of 8 g was prepared following the protocol proposed by Quilliam et al. (1995), including strong anion exchange (SAX) to avoid interference with tryptophan. A volume of 20 µL of the sample was injected into the HPLC-DAD-UV system. The retention time of DA was close to 13 min. The limit of detection (LOD) was 0.1025 µg DA g −1 and the limit of quantification was determined as 3 × LOD and was found to be 0.3415 µg DA g −1 .
A detailed description of the HPLC method and composition of the mobile phase was previously given in [8,9]. A total of 33 shellfish samples were analysed during the investigated period.

Phytoplankton Analyses (LM and SEM)
Phytoplankton community composition was analysed using a Leica inverted light microscope (LM) DMI4000B (Leica Microsystems CMS, Wetzlar, Germany), according to the Utermöhl method [24]. After collecting, seawater samples (250 mL) were preserved with Lugol's solution and stored in the dark until the analyses. Subsamples of 25 mL were settled in counting chambers for at least 24 h. Afterwards, cells were counted along one transect at magnification 400×. Additionally, half of the chamber was examined at a magnification of 200× for more precise identification of less abundant microphytoplankton (>20 µm) taxa. Phytoplankton specimens were identified to the lowest possible taxonomic level. The phytoplankton community was analysed based on four of the most important phytoplankton groups, diatoms, dinoflagellates, coccolithophores, and nanoflagellates, which together accounted for, on average, 96.4% of the total phytoplankton community. A total of 102 seawater samples were analysed by LM.
The composition of Pseudo-nitzschia species was analysed using scanning electron microscopy (SEM, Tescan MIRA3, Brno, Czech Republic), at the depth with the highest Pseudo-nitzschia spp. abundance determined by LM. Organic material from the valves was removed by adding 10% hydrochloric acid (HCl), 30% sulphuric acid (H 2 SO 4 ), concentrated potassium permanganate (KMnO 4 ) and, after 24 h, concentrated oxalic acid (COOH) 2 , after which the sample was rinsed three times in distilled water, as described in detail in [25]. Afterwards, samples were filtered on polycarbonate membrane filters (pore size 2 µm, Nucleopore, Whatman) dried for at least 24 h in a desiccator, coated with gold, and examined with SEM. In each sample, at least the first 30 valves were measured: the length, width, number of fibulae and interstriae in 10 µm, as well as the number of poroids in 1 µm and number of sectors observed within the poroid. The observed morphological characteristics were compared to the available literature data to determine the Pseudo-nitzschia species. In total, 21 samples were analysed by SEM.

Statistical Analyses
To assess the relationship between salinity, as a proxy for river inflow, and nutrient concentrations, as well as the relationship between environmental parameters and phytoplankton groups, a nonparametric Spearman rank correlation was performed in Statistica 13. Additionally, to estimate the phytoplankton community's response to the set of environmental data, distance-based redundant analysis (dbRDA) was performed in PRIMER with an add-on PERMANOVA+ software package. Before running the dbRDA analyses, the phytoplankton species abundance data were log(x + 1) transformed, and the Bray-Curtis similarity matrix was calculated. In contrast, the environmental dataset was normalised, and a Euclidian distance similarity matrix was used.

Phytoplankton Community Structure
Total phytoplankton abundances were generally low until the autumn and winter (November-December) 2016, when a large peak of small diatoms occurred in all depths of the investigated water column. Abundances increased again in late summer/early autumn 2017, when coccolithophores and diatoms started to bloom. During the investigated period, the phytoplankton community was characterised by the lowest abundances in spring ( Figure 2).
Diatoms dominated in abundance and contribution in all investigated layers (0 m, halocline, and 7 m), ranging from 2.6 × 10 3 to 3.5 × 10 6 cells L −1 and contributing 6-100% in samples. The abundance of nanoflagellates ranged from 6.5 × 10 3 to 7.2 × 10 5 cells L −1 and continuously occurring. Coccolithophores reached a very high density of 2.3 × 10 6 cell L −1 , but only occasionally occurred in the phytoplankton community. The abundance of dinoflagellates was low, as well as their contribution to the total phytoplankton density.
The vertical distribution of phytoplankton groups showed a similar pattern in the water column. Diatoms were the most represented group in all depths, with an average contribution of 48-63%. Dinoflagellates and coccolithophores contributed less, not exceeding 20%. Nanoflagellates contributed the most to the surface layer (28%), while in the other two layers their numbers slightly decreased (23% and 22%). Furthermore, the vertical biodiversity was very similar. The largest number of species was recorded at 7 m, where the Shannon diversity index was slightly higher (H' = 1.47) in comparison to the other two layers (H' = 1.38 and 1.42, respectively).
The calculated diversity indices showed slightly higher average values in the spring compared to the other seasons. So, the average Shannon index was H' = 1.64 as opposed to the winter, summer, and autumn seasons (H' = 1.38, 1.51, and 1.26, respectively). A similar pattern was observed with the other analysed indices.
Considering the whole study, 140 phytoplankton species were recorded at Strmica station (presented in Supplementary Table S1). Dinoflagellates were the most diverse group, with 70 recorded species. Diatoms numbered 59 and coccolithophores 10 species. Despite the high diversity, dinoflagellates were very low in abundance, with only two slight enhancements in December 2015 (halocline) and June 2016 (surface). In those samples, Prorocentrum cordatum was the prevailing species, reaching 9.0 × 10 4 cells L −1 (December 2015) and 3.3 × 10 5 cells L −1 (June 2016). The most frequently occurring dinoflagellate was Gymnodinium spp. (77%), with a maximum abundance of 6.6 × 10 4 cells L −1 . The increase in coccolithophore abundance occurred in the summer period in both analysed years. In August 2016, the abundance of Syracosphaera pulchra HOL oblonga type reached a density of 2.8 × 10 5 cells L −1 at the surface. The next year, 2017, the coccolithophore bloom started in August but reached its maximum of 2.3 × 10 6 cells L −1 in October 2017. Unexpectedly, the blooming species was Syracosphaera halldalii [26]; it was the first time this species was documented to have such a high proliferation. The high proliferation of diatoms in winter 2016 was due to an increase in the species Chaetoceros subtilis, which reached a maximum of 3.5 × 10 6 cells L −1 (December 2016), and in February 2017 when the abundance of Bleakeleya notata reached 2.5 × 10 6 cells L −1 at halocline.
The most frequently occurring taxa in the entire phytoplankton community were nanoflagellates recorded in all samples. Moreover, Pseudo-nitzschia and the Chaetoceros genus (87% and 70%, respectively) were widespread, followed by Navicula sp. (75%), Gymnodinium spp. (74%), and Cylindrotheca closterium (65%). Coccolithophores occurred in 53% of the analysed samples. Regarding diatoms, four genera presented more than 71% of the diatom community: Pseudo-nitzschia, (29.  The high proliferation of diatoms in winter 2016 was due to an increase in the species Chaetoceros subtilis, which reached a maximum of 3.5 × 10 6 cells L −1 (December 2016), and in February 2017 when the abundance of Bleakeleya notata reached 2.5 × 10 6 cells L −1 at halocline.

Pseudo-nitzschia Assemblage Seasonality and Diversity
Pseudo-nitzschia species were present in the phytoplankton community throughout the year, with a frequency of appearance of 85%, 82%, and 94% at the surface, halocline, and 7 m depth, respectively. The seasonal distribution of Pseudo-nitzschia spp. showed that, on average, the highest abundance and contribution to the diatom community were recorded during late summer. The annual peak of Pseudo-nitzschia spp. was in September in both analysed years ( Figure 4). In September 2016, the highest cell abundance was determined at halocline (0.91 × 10 6 cells L −1 ), while in 2017 the maximal abundance was 1.84 × 10 6 cells L −1 at the 7 m depth. An increase in abundance was also recorded in the winter, but the abundance was one order of magnitude lower than in the summer, at 3.69 × 10 5 cells L −1 (January 2016) ( Figure 4). The lowest abundance was found in the late winter to spring period (March-April).
The highest contribution of Pseudo-nitzschia spp. to the diatom assemblage was during the bloom event in 2017, as it accounted for 86-99% of the total diatom abundance within four weeks, at 7 m (September-October). In the same period of the previous year, during a high Pseudo-nitzschia spp. abundance, its contribution was significantly lower, not exceeding 67%. Moreover, in February and June 2016, Pseudo-nitzschia spp. dominated the diatom community, contributing up to 86% and 98%, respectively, while diatom abundance was in general low.
The vertical distribution of Pseudo-nitzschia spp. mostly showed both higher abundance (Supplementary Figure S1A) and a more substantial contribution to the diatom community at 7 m depth (Supplementary Figure S1B) in comparison to the surface and halocline depths. The only exceptions, when higher abundance was found at the halocline, were August 2015 and September 2016 ( Figure 4).

Pseudo-nitzschia Assemblage Seasonality and Diversity
Pseudo-nitzschia species were present in the phytoplankton community throughout the year, with a frequency of appearance of 85%, 82%, and 94% at the surface, halocline, and 7 m depth, respectively. The seasonal distribution of Pseudo-nitzschia spp. showed that, on average, the highest abundance and contribution to the diatom community were recorded during late summer. The annual peak of Pseudo-nitzschia spp. was in September in both analysed years ( Figure 4). In September 2016, the highest cell abundance was determined at halocline (0.91 × 10 6 cells L −1 ), while in 2017 the maximal abundance was 1.84 × 10 6 cells L −1 at the 7 m depth. An increase in abundance was also recorded in the winter, but the abundance was one order of magnitude lower than in the summer, at 3.69 × 10 5 cells L −1 (January 2016) ( Figure 4). The lowest abundance was found in the late winter to spring period (March-April).
The highest contribution of Pseudo-nitzschia spp. to the diatom assemblage was during the bloom event in 2017, as it accounted for 86-99% of the total diatom abundance within four weeks, at 7 m (September-October). In the same period of the previous year, during a high Pseudo-nitzschia spp. abundance, its contribution was significantly lower, not exceeding 67%. Moreover, in February and June 2016, Pseudo-nitzschia spp. dominated the diatom community, contributing up to 86% and 98%, respectively, while diatom abundance was in general low.
The vertical distribution of Pseudo-nitzschia spp. mostly showed both higher abundance (Supplementary Figure S1A) and a more substantial contribution to the diatom community at 7 m depth (Supplementary Figure S1B) in comparison to the surface and halocline depths. The only exceptions, when higher abundance was found at the halocline, were August 2015 and September 2016 ( Figure 4).
During the summer bloom in September 2017, Pseudo-nitzschia spp. dominated in the phytoplankton community, reaching 91.8% of total phytoplankton abundance.
Seven Pseudo-nitzschia species were identified by SEM in the period from July 2016 to October 2017. The most frequently occurring were species from the P. delicatisisma/arenysensis complex, which was found in all analysed samples and prevailed in Pseudo-nitzschia assemblage during summer blooms in September 2016 (100%) and 2017 (94%). On average, a high percentage of P. delicatissima/arenysensis was also recorded in other seasons (Figures 5 and 6E,F). Molecular analyses confirmed the dominant species during the bloom event in 2017 as P. cf. arenysensis [27]. Still, as these species are impossible to distinguish based solely on morphological characteristics, throughout the study, we designate them as complex P. delicatissima/arenysensis.
Another commonly occurring species was P. pseudodelicatissima/cuspidata, which significantly contributed to the Pseudo-nitzschia assemblage in winter (January 2016) and summer (August 2017) ( Figures 5 and 6C,D). Pseudo-nitzschia subfraudulenta was recorded in higher percentages during the winter (Figures 5 and 6I,J), while in the summer, it was found to make a small contribution. Similarly, P. calliantha was also determined in higher percentages during the winter and lower in summer ( Figures 5 and 6A,B). Pseudo-nitzschia multistriata was found in all seasons but at a very low contribution, not exceeding 6% ( Figures 5 and 6K,L). Species P. fraudulenta ( Figure 6G,H) and P. lundholmiae were found only on two occasions, contributing 2-7%, and 2-4%, respectively ( Figure 5). In general, high species diversity was recorded throughout all seasons, with five species found in winter, summer, and autumn and four in spring. The species' morphological characteristics are presented in Table 1 and SEM micrographs of the most representative species are given in Figure 6.  Seven Pseudo-nitzschia species were identified by SEM in the period from July 2016 to October 2017. The most frequently occurring were species from the P. delicatisisma/arenysensis complex, which was found in all analysed samples and prevailed in Pseudo-nitzschia assemblage during summer blooms in September 2016 (100%) and 2017 (94%). On average, a high percentage of P. delicatissima/arenysensis was also recorded in other seasons (Figures 5 and 6E,F). Molecular analyses confirmed the dominant species during the bloom event in 2017 as P. cf. arenysensis [27]. Still, as these species are impossible to distinguish based solely on morphological characteristics, throughout the study, we designate them as complex P. delicatissima/arenysensis.
Another commonly occurring species was P. pseudodelicatissima/cuspidata, which significantly contributed to the Pseudo-nitzschia assemblage in winter (January 2016) and summer (August 2017) ( Figures 5 and 6C,D). Pseudo-nitzschia subfraudulenta was recorded in higher percentages during the winter (Figures 5 and 6I,J), while in the summer, it was found to make a small contribution. Similarly, P. calliantha was also determined in higher percentages during the winter and lower in summer ( Figures 5 and 6A,B). Pseudo-nitzschia multistriata was found in all seasons but at a very low contribution, not exceeding 6% (Figures 5 and 6K,L). Species P. fraudulenta ( Figure 6G,H) and P. lundholmiae were found only on two occasions, contributing 2-7%, and 2-4%, respectively ( Figure 5). In general, high species diversity was recorded throughout all seasons, with five species found in winter, summer, and autumn and four in spring. The species' morphological characteristics are presented in Table 1 and SEM micrographs of the most representative species are given in Figure 6. Seven Pseudo-nitzschia species were identified by SEM in the period from July 2016 to October 2017. The most frequently occurring were species from the P. delicatisisma/arenysensis complex, which was found in all analysed samples and prevailed in Pseudo-nitzschia assemblage during summer blooms in September 2016 (100%) and 2017 (94%). On average, a high percentage of P. delicatissima/arenysensis was also recorded in other seasons ( Figures 5 and 6E,F). Molecular analyses confirmed the dominant species during the bloom event in 2017 as P. cf. arenysensis [27]. Still, as these species are impossible to distinguish based solely on morphological characteristics, throughout the study, we designate them as complex P. delicatissima/arenysensis.
Another commonly occurring species was P. pseudodelicatissima/cuspidata, which significantly contributed to the Pseudo-nitzschia assemblage in winter (January 2016) and summer (August 2017) ( Figures 5 and 6C,D). Pseudo-nitzschia subfraudulenta was recorded in higher percentages during the winter (Figures 5 and 6I,J), while in the summer, it was found to make a small contribution. Similarly, P. calliantha was also determined in higher percentages during the winter and lower in summer ( Figures 5 and 6A,B). Pseudo-nitzschia multistriata was found in all seasons but at a very low contribution, not exceeding 6% (Figures 5 and 6K,L). Species P. fraudulenta ( Figure 6G,H) and P. lundholmiae were found only on two occasions, contributing 2-7%, and 2-4%, respectively ( Figure 5). In general, high species diversity was recorded throughout all seasons, with five species found in winter, summer, and autumn and four in spring. The species' morphological characteristics are presented in Table 1 and SEM micrographs of the most representative species are given in Figure 6.

Shellfish Toxicity
During the study, the domoic acid concentration in shellfish samples was below the limit of detection except in the winter of 2015-2016. Low concentrations of DA were measured in December 2015; afterwards, the DA concentration increased and reached the highest value of 1.17 µg g −1 at the beginning of February 2016 ( Figure 7A). During this period, Pseudo-nitzschia species represented up to 72% of the total phytoplankton abundance in the layer at 7 m depth, in mid-February 2016 ( Figure 7B). Additionally, five Pseudo-nitzschia species were identified in the phytoplankton community ( Figure 7B). The most common was P. pseudodelicatissima/cuspidata, found in all samples and representing 50-93% of the Pseudo-nitzschia assemblage. Additionally, P. calliantha was recorded in all samples but with the highest contribution up to 36% in February 2016. Other species, such as P. delicatissima/arenysensis and P. subfraudulenta, occasionally contributed up to 30% while P. lundholmiae was found only once at a low percentage. The increase in DA content in shellfish occurred after the observed increase in Pseudo-nitzschia spp. abundance and its contribution to total phytoplankton abundance. Moreover, changes in species composition were noticed as the percentage of P. subfraudulenta, and P. calliantha rose, but P. pseudodelicatissima/cuspidata remained dominant.

Shellfish Toxicity
During the study, the domoic acid concentration in shellfish samples was below the limit of detection except in the winter of 2015-2016. Low concentrations of DA were measured in December 2015; afterwards, the DA concentration increased and reached the highest value of 1.17 µg g −1 at the beginning of February 2016 ( Figure 7A). During this period, Pseudo-nitzschia species represented up to 72% of the total phytoplankton abundance in the layer at 7 m depth, in mid-February 2016 ( Figure  7B). Additionally, five Pseudo-nitzschia species were identified in the phytoplankton community ( Figure 7B). The most common was P. pseudodelicatissima/cuspidata, found in all samples and representing 50-93% of the Pseudo-nitzschia assemblage. Additionally, P. calliantha was recorded in all samples but with the highest contribution up to 36% in February 2016. Other species, such as P. delicatissima/arenysensis and P. subfraudulenta, occasionally contributed up to 30% while P. lundholmiae was found only once at a low percentage. The increase in DA content in shellfish occurred after the observed increase in Pseudo-nitzschia spp. abundance and its contribution to total phytoplankton abundance. Moreover, changes in species composition were noticed as the percentage of P. subfraudulenta, and P. calliantha rose, but P. pseudodelicatissima/cuspidata remained dominant.

Environmental Factors Associated with Phytoplankton Community
During the study period, temperature and salinity were highly variable. Surface temperature and salinity values were within the ranges of 4.9-26.4 °C and 0.15-31.8, respectively. Halocline was regularly present at a depth from 1 to 4 m depending on the River Krka inflow (https://hidro.dhz.hr). The period with low salinity was November to March, with the lowest values in February in both analysed years. The lowest temperature was measured in January, while the highest temperatures

Environmental Factors Associated with Phytoplankton Community
During the study period, temperature and salinity were highly variable. Surface temperature and salinity values were within the ranges of 4.9-26.4 • C and 0.15-31.8, respectively. Halocline was regularly present at a depth from 1 to 4 m depending on the River Krka inflow (https://hidro.dhz.hr). The period with low salinity was November to March, with the lowest values in February in both analysed years. The lowest temperature was measured in January, while the highest temperatures  Throughout the investigated period, phytoplankton biomass was in general low, with average values recorded from 0.23 ± 0.14 µgL −1 at 7 m depth to 0.37 ± 0.22 µgL −1 at the surface. Regarding the seasonal distribution, slightly higher biomass was found in the autumn to winter, with the peak (1.14 µg L −1 ) recorded in October 2015 in the surface layer (data not shown).
The distribution of nutrients showed elevated concentrations of nitrates and orthosilicates in the winter to autumn, in surface and halocline layers ( Figure 9A,D). The highest values of nitrates (33.82 µM) and orthosilicates (63.23 µM) were recorded in February 2016 during a high river discharge. There was negative correlation of nitrites (r = −0.21, p < 0.05), nitrates (r = −0.56, p < 0.001) and orthosilicates (r = −0.59, p < 0.001) with salinity pointing to the high river influence on environmental conditions, as was previously reported [4,19,28] (Table 2). Orthophosphates, nitrites, and ammonia did not follow the same pattern. A higher concentration of orthophosphates was determined in the spring and summer at a 7 m depth, with the highest value (0.416 µM) in May 2017 ( Figure 9C). The nitrite concentration showed an increase in the autumn in all sampling layers, with the highest value (1.51 µM) recorded in October 2015. The distribution of ammonia did not show any regular pattern as increased concentrations were found occasionally ( Figure 9B). Nutrient ratios were very variable during the investigated period. The N:P ratio was in the range 0.49-1174.27; Si:N was 0.14-79.98 and Si:P was 5.46-1966.10 (data not shown). Throughout the investigated period, phytoplankton biomass was in general low, with average values recorded from 0.23 ± 0.14 µgL −1 at 7 m depth to 0.37 ± 0.22 µgL −1 at the surface. Regarding the seasonal distribution, slightly higher biomass was found in the autumn to winter, with the peak (1.14 µg L −1 ) recorded in October 2015 in the surface layer (data not shown).
The distribution of nutrients showed elevated concentrations of nitrates and orthosilicates in the winter to autumn, in surface and halocline layers ( Figure 9A,D). The highest values of nitrates (33.82 µM) and orthosilicates (63.23 µM) were recorded in February 2016 during a high river discharge. There was negative correlation of nitrites (r = −0.21, p < 0.05), nitrates (r = −0.56, p < 0.001) and orthosilicates (r = −0.59, p < 0.001) with salinity pointing to the high river influence on environmental conditions, as was previously reported [4,19,28] (Table 2). Orthophosphates, nitrites, and ammonia did not follow the same pattern. A higher concentration of orthophosphates was determined in the spring and summer at a 7 m depth, with the highest value (0.416 µM) in May 2017 ( Figure 9C). The nitrite concentration showed an increase in the autumn in all sampling layers, with the highest value (1.51 µM) recorded in October 2015. The distribution of ammonia did not show any regular pattern as increased concentrations were found occasionally ( Figure 9B). Nutrient ratios were very variable during the investigated period. The N:P ratio was in the range 0. 49-1174.27; Si:N was 0.14-79.98 and Si:P was 5.46-1966.10 (data not shown). Contrary to the main phytoplankton groups, which showed weak or no correlations with the environmental parameters, Pseudo-nitzschia showed a significant association with most of them, except orthophosphates ( Table 2). The Pseudo-nitzschia spp. abundance was positively correlated to salinity and temperature, highlighting the development of late summer blooms. Furthermore, a positive correlation was recorded for Pseudo-nitzschia spp. abundance and the Si:N ratio. The negative correlation of Pseudo-nitzschia spp. with nitrates, orthosilicates, and the nutrient ratios of N:P and Si:P, demonstrates the influence of River Krka runoff as the primary source of these nutrients (Table 2). In the layer of 7 m depth, Pseudo-nitzschia spp. were positively correlated with temperature (r = 0.379, p < 0.05, n = 34).  Contrary to the main phytoplankton groups, which showed weak or no correlations with the environmental parameters, Pseudo-nitzschia showed a significant association with most of them, except orthophosphates ( Table 2). The Pseudo-nitzschia spp. abundance was positively correlated to salinity and temperature, highlighting the development of late summer blooms. Furthermore, a positive correlation was recorded for Pseudo-nitzschia spp. abundance and the Si:N ratio. The negative correlation of Pseudo-nitzschia spp. with nitrates, orthosilicates, and the nutrient ratios of N:P and Si:P, demonstrates the influence of River Krka runoff as the primary source of these nutrients ( Table 2). In the layer of 7 m depth, Pseudo-nitzschia spp. were positively correlated with temperature (r = 0.379, p < 0.05, n = 34).
The distance-based redundant analysis (dbRDA) revealed a seasonal difference in phytoplankton assemblage distribution, as summer and winter samples are in general well separated ( Figure 10). The first dbRDA axis was correlated with nitrates (r = 0.772), orthosilicates (r = 0.500), and temperature (r = −0.350), while dbRDA2 was correlated to nitrites (r = 0.695), temperature (r = 0.528), and salinity (r = −0.343). A superimposed bubble plot of the Pseudo-nitzschia spp. community ( Figure 10) showed higher abundances during the summer, indicating that salinity and temperature influenced it most, while during the winter Pseudo-nitzschia spp. were correlated to higher nitrate and orthosilicate concentrations. The distance-based redundant analysis (dbRDA) revealed a seasonal difference in phytoplankton assemblage distribution, as summer and winter samples are in general well separated ( Figure 10). The first dbRDA axis was correlated with nitrates (r = 0.772), orthosilicates (r = 0.500), and temperature (r = −0.350), while dbRDA2 was correlated to nitrites (r = 0.695), temperature (r = 0.528), and salinity (r = −0.343). A superimposed bubble plot of the Pseudo-nitzschia spp. community ( Figure  10) showed higher abundances during the summer, indicating that salinity and temperature influenced it most, while during the winter Pseudo-nitzschia spp. were correlated to higher nitrate and orthosilicate concentrations.

Phytoplankton Community
In this study, the phytoplankton community was analysed, with a particular focus on the potentially toxic diatom genus Pseudo-nitzschia, over two years. Phytoplankton displayed an unusual seasonal cycle, with the spring characterised by the lowest abundances. Two abundance peaks occurred: in the winter, dominated by small diatoms of genus Cheatoceros, and in the late summer and early autumn, when Pseudo-nitzschia and coccolithophores increased. This was not in accordance with the results of Cetinić et al. [4], who found the spring characterised by the highest phytoplankton abundance. In the Northern Adriatic, the spring peak was dominated by nanoflagellates and small diatoms, and a second slight increase in autumn was dominated by nanoflagellates, larger diatoms, and coccolithophores, as observed by Cerino et al. [29].
In our study, diatoms and nanoflagellates co-occurred, while diatoms dominated in abundance and contribution at all investigated depths. This is also typical of other Adriatic areas, where diatoms often dominate the phytoplankton maxima [30][31][32][33]. Turbidity in estuaries could regulate Bubbles present Pseudo-nitzschia spp. abundance in different seasons.

Phytoplankton Community
In this study, the phytoplankton community was analysed, with a particular focus on the potentially toxic diatom genus Pseudo-nitzschia, over two years. Phytoplankton displayed an unusual seasonal cycle, with the spring characterised by the lowest abundances. Two abundance peaks occurred: in the winter, dominated by small diatoms of genus Cheatoceros, and in the late summer and early autumn, when Pseudo-nitzschia and coccolithophores increased. This was not in accordance with the results of Cetinić et al. [4], who found the spring characterised by the highest phytoplankton abundance. In the Northern Adriatic, the spring peak was dominated by nanoflagellates and small diatoms, and a second slight increase in autumn was dominated by nanoflagellates, larger diatoms, and coccolithophores, as observed by Cerino et al. [29].
In our study, diatoms and nanoflagellates co-occurred, while diatoms dominated in abundance and contribution at all investigated depths. This is also typical of other Adriatic areas, where diatoms often dominate the phytoplankton maxima [30][31][32][33]. Turbidity in estuaries could regulate phytoplankton growth and also the community composition, with small diatoms better adapted to the turbid conditions caused by river inflow [1].
The phytoplankton community was mainly composed of fast-growing diatoms such as Chaetoceros spp., Pseudo-nitzschia spp., Bleakeleya notata, Cyclotella spp., and Cylindrotheca closterium. These taxa seem to respond quickly to nutrients introduced by rivers or precipitation [34]. Nanoflagellates were less pronounced; they consistently contributed to the community, but their densities did not exceed diatom abundances. Cetinić et al. [4] reported the maximum and minimum nanophytoplankton abundance in the same period. The highest abundance of nanophytoplankton was associated with the halocline, as was also noticed in previous studies [2]. The dominance of small nanoflagellates, particularly in spring and summer, has been reported in coastal areas of the southern Adriatic Sea [35]. Coccolithophores occasionally occurred in higher abundances. An extreme increase occurred in October 2017, when for the first time a bloom of S. halldalii was recorded [26]. Interestingly, coccolithophores appeared on the surface and at the halocline, as the most variable layers. Regarding the coccolithophore community in this area, Šupraha et al. [28] found a high number of heterococcolith and holococcolith forms in two contrasting seasons. The distinct ecological preferences and seasonality of coccolithophore life phases in the Mediterranean Sea suggest that environmental parameters may regulate the life-cycle transitions [28,36].
The vertical distribution of the phytoplankton groups' abundance showed no distinction between the investigated depths. Moreover, vertical biodiversity was very similar, even though the largest species number was recorded at 7 m depth. Interestingly, despite the low abundance, the diversity indices show higher average values in spring in comparison with the other seasons.
The phytoplankton composition observed in this study area, with 140 taxa, was quite diverse. It was similar to the previous studies carried out in adjacent areas [4,11]. The most diverse phytoplankton group was dinoflagellates. This group was dominated by the potentially toxic species Prorocentrum cordatum (previously named P. minimum), which first appeared in the Adriatic Sea in 1983 [37] and caused blooms in the surface layer due to nutrient availability and favourable light intensity [38]. Other toxic species occurring during the study period were Dinophysis acuminata, D. acuta, D. sacculus, and Phalacroma rotundatum, associated with DSP toxicity. Other dinoflagellates such as Alexandrium cf. minutum, A. cf. tamarense, and Alexandrium sp. are causative organisms of PSP toxicity. The occurrence of toxic dinoflagellates, such as Dinophysis and Alexandrium species was sporadic. Ninčević et al. [13], in their research dealing with a 14-year dataset, observed a correlation between Alexandium species and substantial rainfall, river inflow, and stratification. The same authors confirmed that this area is characterised by the more frequent occurrence of species A. cf. tamarense, D. acuta, and P. cordatum. A spring Alexandrium bloom with increased rainfall, freshwater inflow, and water column stabilisation was reported for the Mediterranean lagoon [39].

Pseudo-nitzschia Assemblage
During our study, Pseudo-nitzschia spp. were present in the phytoplankton community throughout the whole year, with a high frequency of occurrence (87%). Seasonal and vertical distribution showed that, on average, the highest abundance of Pseudo-nitzschia spp. and its contribution to the diatom community was recorded during the late summer/ early autumn and at a depth of 7 m. An increase in abundance was also noted during the winter. Different seasonal patterns of Pseudo-nitzschia spp. have een observed in different seasons: winter [8,11,40], late summer/early autumn [14,30], or both [12,13,41,42], and also in spring [14,[42][43][44][45][46][47]. Throughout our study, the abundance of Pseudo-nitzschia spp. in spring was very low in both years.
Given the high frequency of occurrence and the high contribution to the total phytoplankton and diatom assemblage, Pseudo-nitzschia is distinguished as an important microplankton component of the Krka River ecosystem. Such a high rate (87%) is in accordance with those previously reported for the Central Adriatic Sea [11], but in general higher than those reported for the Northern Adriatic [10,30,41,43]. In accordance with our observations, Pseudo-nitzschia spp. contributed up to 97% to total diatom abundance during blooms in the Northern Adriatic Sea [30], but were less represented during other blooming events [14,40,48,49]. Pseudo-nitzschia has also been reported at other aquaculture sites in the central and southern Adriatic as a significant portion of the microphytoplankton assemblages, although the contribution is lower than observed in this study [50,51].
Throughout the study, the Pseudo-nitzschia assemblage consisted of seven species, morphologically determined by SEM. The majority of them have already been found along the eastern Adriatic coast [9,25,30,46], except P. multistriata and P. lundholmiae, as we reported before in [27]. The morphological characteristics of the identified species in general correspond to those previously reported from the Adriatic Sea [9,14,25,30,46]. The most commonly occurring were species complexes P. delicatissima/arenysensis and P. pseudodelicatissima/cuspidata, which were found in the majority of samples and even in high percentages in different seasons (P. delicatissima/arenysensis: spring to autumn and P. pseudodelicatissima/cuspidata: winter and early summer). Due to crypticism, there is a possibility that these species complexes consisted of at least two species having different environmental preferences. However, P. delicatissima/arenysensis was associated with warmer seasons (summer and autumn), as in the spring the Pseudo-nitzschia spp. abundance was very low. The seasonal peak of cryptic species P. delicatissima and P. arenysensis was usually found in spring [14,45,52], while in our study it was recorded during the late summer/early autumn.
Pseudo-nitzschia pseudodelicatissima/cuspidata were highly present in assemblage throughout the year, but prevailed during the winter. Previously, P. pseudodelicatissima/cuspidata was commonly found in Pseudo-nitzschia assemblage of the Central Adriatic Sea [9,25], but, contrary to our findings, in some Mediterranean areas it was scarce or not detected at all [14,[53][54][55].
Species P. subfraudulenta and P. fraudulenta were more represented in the period from December to February, as reported in our previous studies [9,25]. This is commonly observed for the Mediterranean Sea, as these species usually contributed more in the period from autumn to spring [14,30,45,56].
In the Krka River, the estuary species Pseudo-nitzschia calliantha was found in winter and summer, but in higher percentages during the winter. However, Pseudo-nitzschia multistriata was found in all seasons in low rates. Reoccurrence of P. multistriata during the summer/autumn in the Gulf of Naples was reported by Zingone et al. [57]. Afterwards, it was confirmed in all seasons with the abundance peak also in the winter [56]. In the Adriatic Sea, P. multistriata occurrence was limited to the autumn [33] and winter [14].
Pseudo-nitzschia calliantha is widespread in the Mediterranean Sea and shows broad seasonal distribution. In the Adriatic Sea, it was associated with low temperatures [58], but blooms were recorded in the winter and also in late summer/early autumn [40,48]. Similar seasonality was observed in the NW Mediterranean, where a persistent bloom was present from August to November and from January to March [59]. In other Mediterranean areas, a seasonal peak was observed in spring and autumn [57], and also in the summer [48,50]. The bloom of P. calliantha that occurred from August to October 2007 in the Northern Adriatic was associated with a low concentration of domoic acid in shellfish [48].
During the two-year study period, domoic acid was determined in shellfish only during the winter of 2015-2016. The measured concentrations of DA were very low, far below the regulatory limit of 20 µg g −1 . Concurrently, five Pseudo-nitzschia species were identified: P. pseudodelicatissima/cuspidata, P. calliantha, P. delicatissima/arenysensis, P. subfraudulenta, and P. lundholmiae. Given that P. pseudodelicatissima/cuspidata, followed by P. calliantha and then P. subfraudulenta, were the most frequent and contributed the most, it can be assumed that these species caused the shellfish toxicity. In the investigated area, similar DA values were previously determined in shellfish during the winter of 2006 [8] and 2011 [9]. In 2006, DA was associated with Pseudo-nitzschia bloom. Still, the species composition was unknown [8]. In 2011, particulate DA (pDA) was also measured in plankton net samples, in which P. pseudodelicatissima and P. calliantha were the dominating species. In autumn 2015-winter 2016, low DA was also found in shellfish from Kaštela Bay [25]. Moreover, in that study pDA was confirmed in plankton net samples, and six Pseudo-nitzschia species were identified in the water column. The species composition was similar to what we found during this study, except for P. fraudulenta, which was found only in Kaštela Bay. All of the mentioned species that were found in this study, during shellfish toxicity, are considered to be toxic. However, for the majority of them, the toxicity of the local population was not analysed. To our knowledge, toxicity has only been confirmed for Adriatic strains of P. delicatissima [44], P. multistriata [60], and P. calliantha [61]. The increase in DA concentrations in shellfish tissue (M. galloprovincialis) was recorded 5 days after the increase in Pseudo-nitzschia abundance ( Figure 7A). Although, it is possible the increase in DA occurred even earlier, which we could not detect due to our sampling frequency.
A slight increase in DA concentration occurred after the increase in Pseudo-nitzschia spp. abundance, and its contribution in total phytoplankton. Additionally, changes in species composition were noticed during that event. A similar Pseudo-nitzschia species composition was recorded in a previous study of DA occurrence in the same area, as P. pseudodelicatissima and P. calliantha prevailed over P. subfraudulenta and P. pungens [9].
Thorel et al. [62] observed that Pseudo-nitzschia species composition and toxicity was associated with a low Si:N ratio (<1), while nitrate was still replete. In our study, the occurrence of DA coincided with a period of low Si:N ratio. After the end of January, when the very low Si:N ratio (0.82) was recorded, the highest content of DA was measured in shellfish. An increase in Pseudo-nitzschia toxicity in low-silicate conditions was also observed by [63,64].
Pseudo-nitzschia assemblage was significantly correlated with salinity and temperature (positive correlation) and also showed a negative correlation with nutrients, except orthophosphates. The investigated area is under the considerable influence of the Krka River inflow, causing decreased salinity in the upper layers and permanent halocline throughout the whole year. Nitrates and orthosilicates are mainly originated by river inflow and, thus, are negatively correlated with salinity, as was noticed in this study and previously reported by [4,18,65]. The influence of salinity was evident in the vertical distribution of Pseudo-nitzschia spp., which were more frequent and abundant at 7 m depth than at the surface and halocline. In this layer, Pseudo-nitzschia spp. were associated with higher temperatures and the low turbulence caused by river inflow. In the upper Krka River estuary, domination of Pseudo-nitzschia spp. in the phytoplankton community of deeper layers during summer-autumn was found by Cetinić et al. [4]. An investigation of the time series of toxic phytoplankton species' occurrence with environmental parameters showed that, in this area and the adjacent area, Pseudo-nitzschia spp. were the most abundant during the weaker halocline, and during the summer they were associated with the surface sea temperature [13].
The dbRDA analyses showed a seasonal separation of Pseudo-nitzschia spp. as the winter population was more related to high nitrates, orthophosphate concentrations, and lower salinity, while the summer population was influenced by temperature. However, statistical analyses did not show a significant correlation of Pseudo-nitzschia spp. with phosphates; it was noted that the increase in abundance was preceded by an increase in orthophosphates. Similarly, in the Northern Adriatic, the winter/spring Pseudo-nitzschia population was associated with low temperatures and increased orthosilicate concentrations, while the summer/autumn population was related to higher temperatures [46]. Pseudo-nitzschia blooms were positively correlated to phosphate concentration [40] and the increased nutrients caused by high freshwater inflow [64,66].
Our nutrient data (total inorganic nitrogen, orthophosphates, and total phosphorus), and the data of phytoplankton biomass (chl a) are within the nutrient range gained from National monitoring of transitional waters. Moreover, the entire water body (where Strmica station is located) was characterised by a very good environmental status [67]. Additionally, considering dissolved organic matter (DOC), this area is categorised as oligotrophic [67]. Regarding the observed phytoplankton diversity and entire phytoplankton community structure in this study, we presume that this area is not influenced by intense anthropogenic pressure and can be characterised as oligotrophic. Although previous investigations stated differently [3,4,11,18], all these studies were focused on the lower part of the estuary, which is under the anthropogenic influence of the city of Šibenik.
To summarise; the phytoplankton community was quite diverse, with 140 identified taxa. Phytoplankton abundance showed two major peaks. The winter peak was more pronounced than that observed in the summer/autumn. Diatoms were the dominant group throughout the investigated period. Vertical distribution did not show a significant difference in the phytoplankton community structure. Pseudo-nitzschia spp. was present in the phytoplankton community throughout the year, with a higher frequency and abundance in the less variable layer at a 7 m depth. An increase in abundance occurred during the winter, while the development of bloom was noted in the late summer/early autumn. Seven Pseudo-nitzschia species were determined by SEM throughout the investigated period. Winter assemblage was characterised by P. pseudodelicatissima/cuspidata, P. calliantha, and P. subfraudulenta and was associated with domoic acid occurrence in shellfish. Although the DA production was related to the winter species composition, DA production could be triggered by low Si:N values. The late summer/early autumn bloom was characterised by the dominance of P. delicatissima/arenysensis. Concerning environmental parameters, Pseudo-nitzschia spp. was associated with higher salinity and temperature. The development of the winter Pseudo-nitzschia population was associated with high nitrate and orthosilicate concentrations, while in the summer/autumn it was affected more by temperature. According to the development of winter and late summer/autumn blooms, special attention should be paid to monitor Pseudo-nitzschia species and domoic acid in these seasons. Further morphological, molecular, and toxicological studies of Pseudo-nitzschia genus are needed in this area. Funding: This study was fully funded by Croatian Science Foundation under the project IP-2014-09-3606 "Marine plankton as a tool for assessment of climate and anthropogenic influence on the marine ecosystem, MARIPLAN".