The Gulf of Trieste is a small (~ 500 km²) and shallow (maximum depth 25 m) basin in the northern part of the Adriatic Sea. In this area, freshwater inputs and atmospheric forcing greatly influence seawater temperature, salinity and water column stratification (Malačič and Petelin, 2001). The area is characterised by a marked seasonal cycle of seawater temperature (from winter minima of 8°C to summer maxima of 28.4°C) and strong salinity gradients (from 24.0, in spring during high riverine discharge, to 38.3) (Celio et al., 2006; Kralj et al., 2019). Typically, in winter, the water column is well-mixed, whereas during spring, freshwater input and surface heating lead to thermohaline stratification. The period between May and September is characterised by strong density gradients and the prevalence of respiration processes at the bottom layer, which determine low oxygen concentration and occasionally hypoxia events (Faganeli et al., 1985; Malej and Malačič, 1995; Kralj et al., 2019). In autumn, convective and mechanical mixing, induced by water cooling and wind, disrupt the vertical stratification, oxygenate the bottom water and distribute the re-generated nutrients to the entire water column.

The main riverine input in the Gulf of Trieste derives from Isonzo/Soča River on the north-western coast (Cozzi et al., 2012), which controls the salinity and nutrient concentration of the system with a highly variable outflow. On a seasonal scale, however, spring and autumn are generally characterised by the highest river discharges (due to snowmelt and rain, respectively), while drought periods occur during winter and summer (Comici and Bussani, 2007).

The trophic status of the gulf also depends on the prevailing circulation patterns and not only on the intensity of the Isonzo River discharge rate. Circulation in the gulf is mainly cyclonic at the transitional and lower layer (10 m – bottom), while the surface layer (0 – 4 m) is affected by wind conditions (Stravisi, 1983). There are mainly two dominant winds: the SE Scirocco and the NE Bora (Stravisi, 1977; Stravis, 1983). Bora-induced circulation is more frequent in autumn and winter, and it generates a cyclonic gyre at the surface layer, which causes a fast outflow of riverine waters from the gulf. When the Scirocco wind blows, instead, anticyclonic surface circulation favours eastward spreading of nutrient-enriched freshwater, which increases primary production (Cantoni et al., 2003; Querin et al., 2006).

The Gulf of Trieste is also influenced by the Eastern Adriatic Current (EAC), a current flowing northward along the Croatian coast and advecting warmer, saltier, and more oligotrophic waters coming from the Ionian Sea (Poulain et al., 2001). The ingression of EAC is more frequent in the cold seasons, when a cyclonic circulation is present, which can lead to oligotrophic conditions of the gulf.

Sampling was performed at the Long-Term Ecological Research (LTER) station C1 (45°42’2” N and 13°42’36” E, maximum depth 17.5 m) located in the Gulf of Trieste ( Figure 1 ). From March 2015 to March 2019, discrete seawater samples were collected on a seasonal basis with 5-L Niskin bottles at four depths (0.5, 5, 10, 15 m) for nutrient, phytoplankton and primary production analyses.

Seawater temperature and salinity were recorded monthly by a CTD probe model Sea-Bird Electronics SBE 19plus SeaCAT profiler. These data were compared to the monthly means (1998 - 2019) that were recorded with the following probes: from January 1998 to October 1999 with an Idronaut mod. 401 probe; from November 1999 to September 2003 with an Idronaut mod. 316 probe. These CTD probes were calibrated with an interval of 6 - 12 months. From October 2003 onwards, a Sea Bird Electronics SBE 19 plus SeaCAT was used that was calibrated every year.

To have an indication of the water column stability the Brunt-Väisälä (B-V) frequency was applied. This frequency is a measure of the natural oscillation resulting from vertical displacement of a neutrally buoyant body. Higher values indicate a large density gradient and a strong stratification, while near-zero and negative values indicate unstable conditions. For the calculation, the Matlab routine called sw_bfrq (Revision 1.12, dated 1994/11/15, Copyright CSIRO, Phil Morgan 1993) contained in the SEAWATER package was used. This routine requires as input temperature (ITS-90), salinity and pressure (dbar) data and returns the Brunt-Väisälä square frequency (N²) at mid-latitudes. Results were plotted using ODV (Schlitzer, 2021) and weighted-average gridding (Ocean Data View User’s Guide, Version 5.6.0, 2022).

From March 2015 to September 2016, Photosynthetic available radiation (PAR) was recorded monthly by a PNF-300 Profiling Natural Fluorometer (Biospherical Instruments Inc., San Diego, USA), whereas from October 2016 onwards, a BNF-2102P S/N 1020 was used. These data were compared to the monthly means (November 1999 – September 2019) that were recorded with the above-mentioned probes. In 2015 the probe was sent to the USA for calibration and several monthly PAR data were missing, therefore we chose to exclude the year 2015 from the graphical representation of the monthly data.

Samples for the determination of dissolved inorganic nutrient concentrations (nitrite, N-NO₂; nitrate, N-NO₃; ammonium, N-NH₄; phosphate, P-PO₄; and silicate, Si-Si(OH)₄) were pre-filtered through glass-fibre filters, pore size 0.7 µm (Whatmann GF/F), stored at -20°C and analysed on a four-channel continuous segmented flow analyser (QuAAtro, Seal Analytical), using standard colorimetric methods (Hansen and Koroleff, 1999). To highlight nutrient limitation for microalgal growth we applied the following inequalities to our ratios: Si:N < 0.80 indicates Si limitation (Brzezinski, 1985), N:P > 22 indicates P limitation and N:P < 13 indicates N limitation (Hillebrand and Sommer, 1999). We also applied the Redfield ratio N:P < 16 to highlight a slight N limitation (Redfield, 1958). With the exception of January 2017, sediment sampling was carried out within a few days (-2/+2) from water sampling. Due to bad weather conditions, in January 2016 and April 2017 up to 9 days passed between the two samplings.

Virtually undisturbed sediment cores were collected by an automatic KC Haps bottom corer (KC-Denmark) using polycarbonate sample tubes (13.3 cm i.d. with a sample area of 127 cm²); from one sediment core, in a N₂-filled chamber the overlying water was sampled with a syringe, filtered through Millipore Millex HA 0.45 µm pore size cellulose acetate filters and collected in acid-precleaned vials which were stored at -20°C until nutrient analysis, while from other 3 sediment cores, the uppermost oxic layer (<1 cm) was sampled, pooled, homogenised and subsampled for the analysis of microphytobenthos and benthic primary production. At the moment of sampling, PAR, sea water temperature and salinity were measured as already specified for the water sampling.

For phytoplankton (2-200 µm) analysis, samples were collected in 500 mL-dark bottles and preserved with pre-filtered and neutralized 1.6% formaldehyde (Throndsen, 1978). Cell counts were performed following the Utermöhl method (Utermöhl, 1958). A variable volume of seawater (10-50 mL) was settled depending on cell concentrations. Cell counts were performed using an inverted light microscope (Olympus IX71 and LEICA BMI3000B) equipped with phase contrast. Cells (minimum 200) were counted along transects (1–2) at a magnification of 400×. In addition, one half of the Utermöhl chamber was also examined at a magnification of 200×, to obtain a more correct evaluation of less abundant phytoplankton taxa. Phytoplankton specimens were identified to the lowest possible taxonomic level referring to Tomas (1997); Bérard-Therriault et al. (1999); Horner (2002); Young et al. (2003) and Malinverno et al. (2008). Species/genus names were checked for validity against AlgaeBase (Guiry and Guiry, 2022). The biovolume of phytoplankton cells was calculated according to Edler (1979) and Hillebrand et al. (1999). Cell volumes were converted to carbon content using the formula introduced by Menden-Deuer and Lessard (2000).

Cyanobacteria were not included in this study, however on average (over a 12-month study) they accounted for 13.3% of the total phytoplankton biomass in the same site (Cibic et al., 2018a), and their contribution to PP was previously presented and thoroughly discussed (Cibic et al., 2018b). Phytoplankton community loss factors (grazing pressure, viral infection) were not considered in this study.

In this study, with the term microphytobenthos – MPB, we refer to the microscopic eukaryotic algae (diatoms, dinoflagellates, flagellates, etc.), and prokaryotic photosynthetic organisms, such as filamentous cyanobacteria. For MPB analyses, three aliquots of homogenized sediment (2 cm³) were withdrawn using a syringe and directly fixed with 10 mL of formaldehyde (4% final concentration) buffered solution CaMg(CO₃)₂, in pre-filtered bottom seawater (0.2 μm filters). After manual stirring, 20 μL aliquots of the sediment suspension were drawn off from the slurries and placed into a counting chamber. Only cells containing pigments and not empty frustules were counted under a Leitz inverted light microscope (Leica Microsystems AG, Wetzlar, Germany) using a 32× or 40× objective (320× or 400× final magnification) (Utermöhl, 1958). When possible, at least 200 cells were counted per sample to evaluate rare species, too. The microalgal taxonomy was based on AlgaeBase (Guiry and Guiry, 2022) and WoRMS (WoRMS Editorial Board, 2022) websites. The qualitative identification of MPB assemblages was carried out using floras listed in Cibic and Blasutto (2011). Quantitative data are reported as cells cm^-3 of wet sediment (cells cm^-3) and as Relative Abundance (RA). To estimate the biomass (expressed as μg cm^-3), the biovolume of MPB cells was calculated according to Hillebrand et al. (1999). Afterwards, the MPB biomass was obtained multiplying the abundance (cells cm^-3) by the carbon content of each counted cell using the formulas introduced by Menden-Deuer and Lessard (2000).

We estimated gross primary production (GPP) in water samples and in surface sediments by the ¹⁴C uptake (Steeman Nielsen, 1952). GPP represents the sum of net primary production (NPP) and community respiration (CR). The ¹⁴C technique measures values between GPP and NPP, depending on the incubation time: shorter incubation times are closer to GPP whereas incubation times ≥ 6 h are closer to NPP (Gazeau et al., 2004). In our study, the incubation time was about 2 hours for the water samples and 45 minutes for the sediment slurries, therefore in both cases a GPP rate was measured. From March 2015 to March 2019, pelagic primary production (PPw) was estimated in situ. Water samples were poured into three light and one dark 70 mL polycarbonate carboys (Nalgene) per each depth (0.5-5-10-15 m). The samples were kept in the dark for 30 minutes to stop residual photosynthetic activity. Subsequently, 0.22 MBq (6 μCi) of NaH¹⁴CO₃ (DHI, Denmark) was added to each carboy. The samples were then fixed on a rosette, lowered to the corresponding sampling depth, and incubated for 2 h around noon. At the end of the incubation, the samples were stored in dark and cold conditions, and immediately transferred to the laboratory. From each sample, an aliquot of 25 mL was filtered through 0.2 μm polycarbonate filters (Nuclepore) applying low vacuum pressure (100 mmHg) in order to avoid cell damage. The filters were placed in 6 mL plastic scintillation vials (Perkin Elmer), acidified with 200 μL of HCl 0.5 M (Cibic and Virgilio, 2011) to remove the residual ¹⁴C-bicarbonate not assimilated by the phytoplankton, and an aliquot of 5 mL of scintillation cocktail (Filter Count; Perkin-Elmer) was added to each vial.

Benthic primary production (PPs) was estimated in the laboratory from ¹⁴C-incubation of slurries. An aliquot of 10 cm³ of homogenised surface sediment was withdrawn with a syringe, resuspended in 190 mL of overlying filtered seawater (0.2 μm filter), withdrawn from undisturbed sediment cores, and inoculated with 20 μCi (0.74 MBq) of NaH¹⁴CO₃ (DHI, Denmark) (Steemann Nielsen, 1952). After stirring, the slurry was transferred into 21 glass vials containing 9 mL which were divided as follows: 3 replicates to assess the sediment matrix effect, 3 dark replicates and 3 replicates for each of the 5 light intensities used. In a thermostatic chamber the samples were incubated at the in situ temperature under a gradient of light intensities (20–50–100–200–500 μE m^-2 s^-1) and after 45 min carbon incorporation was stopped by adding 200 μL of HCl 5 N (final HCl concentration 0.11 N) (Cibic and Virgilio, 2010). Subsequently, samples were treated as described by Cibic et al. (2008). In this study only PPs rates, referred to data obtained at in situ light conditions, are presented and discussed.

Disintegrations per minute (DPM) were measured by a QuantaSmart TRI-CARB 2900 TR Liquid Scintillation Analyzer (Packard BioScience, USA) including quenching correction, obtained using internal standards. Assimilation of carbon was calculated as described by Gargas (1975), assuming 5% isotope discrimination. Activities of the added NaH¹⁴CO₃ and inorganic carbon concentration (tCO₂) were calculated on the basis of total alkalinity measured in the same samples. Standard deviation (SD) of three replicate values was below 25% except for rates close to zero for which SD was over 50%.

Volumetric PPw data were converted to areal data using standard trapezoidal integration over the top 15 m of the water column. Total incident photosynthetically active radiation was continuously recorded at the OGS laboratory at S. Croce using the scalar PAR surface reference sensor QSR2100 (Biospherical Instruments). PP hourly rates were then converted into daily rates considering the amount of light during incubation in situ, and the total amount of light during the sampling day computed by integrating the PAR data over the hours of daylight.

Physical data (sea water temperature, salinity and PAR) of the sampling period were compared to the monthly mean climatology. The latter was computed on seawater temperature and salinity data recorded between 1998 and 2019, and on PAR data measured from November 1999 to September 2019. The monthly mean climatology was computed for each parameter at four sampling depths (0.5, 5, 10 and 15 m), except for surface PAR measurements that were carried out at 1 m instead of 0.5 m in order to reduce light scattering in the coastal waters. Whenever CTD and PAR measurements were performed with different sampling frequency (weekly and biweekly), data were averaged to a monthly value.

All statistical analyses, except for the correlations, were performed using PRIMER 7.0.21 (Clarke et al., 2014). To infer the temporal changes, the samples were grouped based on the sampling season considering the meteorological calendar. To highlight seasonal variability in the distribution of the phytoplankton and MPB biomass, a PERMANOVA test was used. “Season” was applied as a fixed factor in a one-way analysis. Unrestricted permutation of raw data, 9999 permutations and Monte Carlo pairwise comparisons were applied. To highlight relationships between abiotic and biotic variables, a Spearman rank correlation analysis (r) was performed using the software Past 4 (Hammer et al., 2001). Only statistically significant (p < 0.05) results are presented and discussed. For each water depth and for the sediment, the biological data were divided into two separated matrices, with the abundances and biomasses of phytoplankton and MPB taxa respectively. For each sampling layer, an additional matrix with the abiotic parameters (temperature, salinity, PAR, inorganic nutrients) was constructed. Univariate diversity analysis was applied to phytoplankton and benthic diatom abundances, after removing undetermined forms, considering richness (d, Margalef, 1986), equitability (J′, Pielou, 1966), diversity (H’(log₂), Shannon and Weaver, 1949) and dominance (λ, Simpson, 1949). Before the multivariate analysis, each biotic matrix was log (x + 1) transformed. The data were then analyzed using cluster analysis (performed with the complete linkage clustering algorithm) and for the biotic matrices, a Bray-Curtis similarity was applied.

A Principal Component Analysis (PCA) was carried out on environmental data of each sampling layer (or sediment) in order to visualize the temporal distribution of main abiotic variables (salinity, temperature and nutrients).

To visualize differences in taxa assemblages among the different samplings, a non-metric multidimensional scaling ordination (nMDS) (Kruskal and Wish, 1978) was performed separately on the Bray-Curtis similarity matrix of phytoplankton at each sampled depth (0.5, 5, 10 and 15 m) and of MPB, after log (x + 1) transformation of the data. To highlight which taxa mainly contributed to the temporal variation of the assemblages, the taxa with the highest relative biomass (average ≥ 1% for the phytoplankton and ≥ 2% for the MPB, over the study period), were overlaid on the nMDS plot. The normalized (z-standardization) environmental variables (i.e. temperature, salinity, PAR and dissolved inorganic nutrients: N-NO₂, N-NO₃, N-NH₄, P-PO₄ and Si-Si(OH)₄; for sediments also TOC and TN data presented in Franzo et al. (2019) were used) were fitted as supplementary variables (vectors) onto ordination spaces to investigate their effects on community structure, using an Euclidean distance matrix for physical-chemical data. The analyses were carried out on biomass values, since is the biomass, over the abundance that drives the primary production rate. On all the communities, we performed multivariate analyses to examine seasonal composition changes in relation to environmental factors. The following abiotic variables were considered for each sampling layer: i) physical parameters (temperature, salinity and light intensity); ii) inorganic nutrients (N-NH₄, N-NO₂, N-NO₃, P-PO₄, Si-Si(OH)₄); and TOC and TN for sediments only. The significancy of differences among a priori-fixed seasons (winter: December, January, February; spring: March, April, May; summer: June, July, August; autumn: September, October, November) was tested by ANOSIM (ANalysis Of SIMilarity) (Clarke et al., 2014). The resulting pair-wise R-values give an absolute measure of separation among groups with zero (0) indicating no difference among groups, while one (1) indicating that all the samples within groups are more similar to one another than any samples from different groups.

Comparing the seawater temperature to climatological (1998 - 2019) data, much lower values than the monthly means were recorded particularly in August 2015 (up to -4.73°C at the surface) and July 2018 (up to -2.16°C at the surface), and much higher values in May 2018 (up to +6.16°C at the bottom) and July 2015 (up to +4.30°C at the 5-m layer). ( Figures S1A–D ) In the sediment overlying water, the temperature, measured during seasonal samplings, ranged between a minimum of 9.45°C in January 2016 and a maximum of 22.16°C in July 2016 ( Table S1 ).

The minimum of salinity recorded in May 2016 at the surface was -6.76 lower than the climatological (1998-2019) monthly mean and lower salinity persisted in the water column also in the following month. Higher values than the climatological monthly means were recorded at the surface particularly in spring-summer of 2018 ( Figures S1E–H ). In the sediment overlying water, the salinity, measured during seasonal samplings, ranged between a minimum of 36.99 in October 2015 and a maximum of 38.06 recorded in March 2019 ( Table S1 ).

The degree of the water column stability was obtained by the Brunt-Väisälä (B-V) frequency that allows to identify the position of the seasonal pycnocline and further determine periods of substantial vertical gradients in water properties. The maximum B-V frequency along the water column was generally observed in late spring-summer, indicating the predominance of the thermal effect ( Figure 2 ). However, marked differences among the investigated years were noticed. In particular, in 2016 a strong water stratification established at the surface in May, while a minimum B-V frequency was recorded in July. In 2017 maxima interested almost the whole water column both in June and July, with a minimum in July at -5m. The year 2018 was characterized by peaks at mid-water column from April to September, while in November 2018 a pycnocline appeared at about 2.5-m depth.

PAR reached the maximum (1654.8 µE m^-2 s^-1) of the study period at the surface in May 2017, while the minimum (2.5 µE m^-2 s^-1) was recorded at the bottom in January 2019 ( Figures S2A–D ). Except for the surface, PAR was, on average, lower than the monthly mean at the other sampling depths. In the sediment overlying water, PAR, registered during seasonal samplings, was at its lowest (6.0 µE m^-2 s^-1) in October 2016, while it reached the maximum (126.3 µE m^-2 s^-1) in June 2015 ( Table S1 ).

Nutrient concentrations analysed at the four water depths and in the overlying water, and N:P and Si:N ratios are shown in Table S2 . N-NH₄ displayed two peaks (5.0 and 2.6 µM) in October 2015 at -5 m and in November 2018 at -15 m, respectively, and three peaks (5.6, 4.9 and 5.9 µM) in the overlying water in October 2017, April 2018 and November 2018, respectively ( Figure S3A ). Also N-NO₂ displayed high concentrations along the water column and in the sediment overlying water in January 2016, 2017, and 2018 ( Figure S3B ) whereas N-NO₃ in the overlying water in January 2016, and in November 2018 ( Figure S3C ). Similarly to N-NH₄, P-PO₄ reached its maximum (1.1 µM) in October 2015 at -5 m and in July 2017 (0.5 µM) in the overlying water ( Figure S3D ). Si-Si(OH)₄ reached the highest concentrations in November 2018 along the entire water column, while the minimum was recorded at -5 m in October 2015. In the overlying water, Si-Si(OH)₄ displayed three major peaks in July and October 2017, and in November 2018 ( Figure S3E ).

Applying the ratios (Redfield, 1958) as proposed by Hillebrand and Sommer (1999), and the Si:N ratio reported by Brzezinski (1985) to our data, we found Si limitation at the top 10 m-layer in October 2015, in April 2017, at the surface layer in January 2016 and in the overlying water in April 2018. N-depleted conditions were observed particularly in June 2015, April 2017, July 2018 and March 2019, whereas a slight N-limitation occurred in the overlying water in October 2016 and July 2017 ( Table S2 ). When N was not depleted, P became the limiting nutrient for the phytoplankton growth in most of the other samplings or water layers.

Considering the whole study period, the phytoplankton community was dominated by flagellates in terms of absolute abundance, and by diatoms in terms of biomass. The average abundance of the whole study period was equal to 1.0·10⁶ ± 8.3·10⁵ cells L^-1 for flagellates and to 3.0·10⁵ ± 3.8·10⁵ cells L^-1 for diatoms. Three peaks in abundance were reached in June 2015 at the surface, when the flagellates accounted for 3.4·10⁶ cells L^-1, and in July 2017 and 2018, when this group reached 4.3·10⁶ and 2.9·10⁶ cells L^-1, respectively ( Figure S4A ). In July 2018, a peak was observed also at -10 m when flagellates and diatoms were co-dominant displaying similar abundances ( Figure S4C ). In addition, diatoms reached other two peaks in July 2018 at -5 m ( Figure S4B ) and April 2016 at -10 m ( Figure S4C ). Compared to the other two groups, dinoflagellates and coccolithophores represented a negligible portion of the phytoplankton community ( Figures S4A–D ).

In terms of biomass, the composition of the phytoplankton community differed substantially compared to what emerged from the abundance data. The diatoms dominated the community with an average biomass (considering the whole study period and all sampling depths) equal to 31.2 ± 45.4 µg C L^-1 while that of flagellates was 19.5 ± 14.2 µg C L^-1. Diatom maxima were observed over the whole water column in April 2017 and January 2018 ( Figures S4E-H ), while the maxima of flagellates and dinoflagellates were recorded at the surface in June 2015 and November 2018, respectively ( Figure S4E ). A significant seasonal variability in the distribution of the phytoplankton biomass was highlighted by PERMANOVA at all depths, except for the 15-m layer ( Table 1 ).

To have an indication of the community composition over the study period, the phytoplankton abundance and structure at the four depths were integrated over the water column and the relative abundance of taxa was calculated. On average, small flagellates (<10 µm) of uncertain taxonomic identification dominated the phytoplankton community (60.60%), followed by small (< 20 µm) diatom species belonging to the genus Chaetoceros (9.87%), undetermined cryptophyceae (5.89%) and the coccolithophore Emiliania huxleyi (3.47%) ( Figure 3A ). A very different picture was obtained considering the phytoplankton biomass. When large diatoms occurred in the water column, particularly Lauderia annulata and Leptocylindrus spp., but also Chaetoceros, they prevailed over the other groups in terms of carbon content ( Figure 3B ).

Univariate diversity indices applied to the phytoplankton community, integrated over the water column, revealed that the highest value of richness (d = 4.441) was observed in November 2018 when the highest number of species was identified in the water column (S = 78). The second maximum of richness was obtained in October 2016 (d = 4.389, S = 73), while the highest values of diversity and equitability were calculated in January 2018 (H’ = 3.452, J’ = 0.561) in correspondence to the lowest value of dominance (λ = 0.182) ( Table 2 ).

Over the study period, MPB displayed two abundance peaks (1.6·10⁵ and 1.4·10⁵ cells cm^-3) in March and June 2015, respectively, while two minima (3.0·10⁴ and 2.9·10⁴ cells cm^-3) were recorded in October 2015 and December 2016, respectively. In contrast to the phytoplankton community, the MPB biomass dynamics were consistent with those of the abundances ( Figure S5A ). Two biomass peaks (50.8 and 48.9 µg C cm^-3) were reached in March and June 2015, while the minima (9.3 and 8.4 µg C cm^-3) were detected in October 2015 and December 2016 ( Figure S5B ). The MPB biomass significantly differed among seasons but not when the pairwise comparisons were considered (Monte Carlo test, Table 1 ).

Considering the whole study period, the MPB community was dominated by diatoms (94.78%), whereas diatom spores (2.27%), cyanobacteria (2.02%), undetermined flagellates (0.63%), dinoflagellates and their resting stages (0.30%) represented minor fractions of the benthic microalgal community. Focusing on diatoms only, Nitzschia was the most abundant genus over the study period (29.34%), followed by Paralia (15.49%) and Gyrosigma (11.28%). The tychopelagic diatom species Paralia sulcata was the most represented over the study (15.41%), Nitzschia fasciculata (5.73%) was particularly abundant in summer, other well represented species were Gyrosigma acuminatum (5.40%), G. spencerii (3.32%) and Nitzschia sigma var. sigmatella (3.01%) ( Figure 4 ).

Univariate diversity indices applied to the benthic diatom community revealed that the highest value of richness (d= 3.368) and diversity (H’=4.217) were observed in July 2018 in correspondence to the highest number of species (S= 36). In October 2016 the highest value of dominance (λ = 0.222) and the lowest equitability (J’ = 0.692) were obtained in correspondence to a major density of Nitzschia sigma var. sigmatella. Similarly, the second maximum of dominance (λ = 0.191) and the second minimum of equitability (J’ = 0.701) were obtained in January 2018 due to a high relative abundance of Paralia sulcata ( Table 2 ).

The absolute PPw maximum of the study period was estimated at the surface in November 2018 (6.71 ± 0.82 µg C L^-1 h^-1) ( Figure 5A ) in low light conditions, while a second maximum was obtained in July 2017 at -5 m (6.46 ± 0.52 µg C L^-1 h^-1) ( Figure 5B ). Other three relative maxima were estimated at the surface in January 2018 (5.13 ± 0.07 µg C L^-1 h^-1) and July 2017 (4.76 ± 0.74 µg C L^-1 h^-1), and at -5 m in October 2016 (4.74 ± 0.85 µg C L^-1 h^-1). In general, at the surface and -5 m PP dynamics were uncoupled with the light intensity, whereas at -10 m and -15 m the relative PP peaks were more in line with the PAR availability at these depths ( Figures 5C, D , S2 ). PP minima (<0.2 µg C L^-1 h^-1) were reached in April 2017 at the surface and in October 2017 at the bottom ( Figures 5A, D ).

The primary production depth-integrated rates (PPi) displayed the absolute maximum in July 2017 (69.30 mg C m^-2 h^-1) and the second highest value in July 2018 (55.89 mg C m^-2 h^-1), while the minima were obtained in March 2019 and April 2018 (11.47 and 13.17 mg C m^-2 h^-1, respectively) ( Figure 6 ).

The highest PPs rate was estimated in January 2018 (28.50 ± 2.08 mg C m^-2 h^-1), a second maximum was reached in October 2016 (17.55 ± 1.17 mg C m^-2 h-¹) while negative values were obtained in October 2015 (-20.88 ± 3.93 mg C m^-2 h^-1), December and January 2016 (-7.58 ± 1.67 and -1.20 ± 0.30 mg C m^-2 h^-1 respectively) in low light conditions ( Figure 6 ).

We made an estimate of the total primary production (PPt) rate per m². Therefore, PPs rates were added to PPi values. Even when PPs values were highly negative, the sum of PPs + PPi was consistently > 12 mg C m^-2 h^-1, except for winter 2016 - 17 when the water sampling was much delayed compared to the sediment sampling ( Table 3 ). To obtain an estimate of PPs to PPt rate, we first excluded the 3 negative PPs values then calculated the percentage of the benthic contribution to the overall primary production. The highest benthic contribution was obtained in January 2018 (43.2%) and October 2016 (28.3%), while the mean PPs value over the study period was 11.3% ( Table 3 ).

To highlight how the physical and chemical features of the water column and surface sediments affected the seasonal development of pelagic and benthic phototrophs and the overall primary production, a Principal Component Analysis (PCA) was first carried out on the five abiotic data matrices, separately. The total variance explained by the first two PC axis displayed increasing values from the surface to the 15-m depth layer, accounting for 64.2%, 69.6%, 70.0% and 73.8% of total variance for the 0.5 m, 5 m, 10 m and 15 m layer, respectively, whereas for the sediment PC1 + PC2 explained 63.3% of the total variance (data not graphically shown). For the PCA performed on the 0.5-m data, temperature (0.60) and salinity (-0.50) were the predominant elements of the first factor, while the major contributors to the second one were N-NO₃ (-0.66) and silicates (-0.64). For the PCA performed on the 5-m data, N-NO₃ (-0.47) and silicates (-0.46) were the predominant elements of the first factor, while the major contributors to the second one were P-PO₄ (-0.53) and N-NH₄ (-0.51). For the PCA performed on the 10-m data, N-NO₃ (0.53) and N-NO₂ (0.48) were the predominant elements of the first factor, while the major contributors to the second one were temperature (-0.56) and salinity (0.54). For the PCA performed on the 15-m data, N-NH₄ (0.53) and salinity (-0.43) were the predominant elements of the first factor, while the major contributors to the second one were N-NO₂ (0.59) and N-NO₃ (0.51). Finally, for the PCA performed on sediment data, silicates (0.47) and P-PO₄ (0.44) were the predominant elements of the first factor, while the major contributors to the second one were N-NO₂ (0.57) and N-NO₃ (0.45). On all PCA plots, the four groups of seasonal samplings significantly differed (R_ANOSIM varying from 0.366; p = 0.2% for the 5-m layer to R = 0.611; p = 0.1% for the sediment).

The nMDS analysis based on the biomass of either the phytoplankton (separately for each layer) or MPB revealed temporal differences among samplings ( Figures 7A–E ). The superimposed PP and abiotic variables indicated the main discriminating factors responsible for the separation of samplings. At the surface layer, the phytoplankton assemblages significantly differed among seasons (R_ANOSIM = 0.571; p = 0.1%) ( Figure 7A ). During spring, the diatom biomass prevailed in April 2017 but mostly it was co-dominant with that of flagellates, and the measured photosynthetic rates were quite low (1.13 ± 0.99 µg C L^-1 h^-1). In summer, flagellates were consistently dominant and higher mean PP rates were estimated (3.30 ± 1.50 µg C L^-1 h^-1). In autumn, diatoms contributed to higher PP (October16 and October17), whereas the absolute maximum PP (6.71 µg C L^-1 h^-1) obtained in November 2018 was due to the highest biomass of small dinoflagellates which co-occurred with flagellates. In January 2018, the second maximum PP (5.13 µg C L^-1 h^-1) was estimated on the account of the diatoms Lauderia annulata and several Chaetoceros species. In the other two winter samplings, PP rates were low due to a minor overall phytoplankton biomass and low light availability.

Similarly, at the 5-m layer, the phytoplankton assemblages significantly differed among seasons (R_ANOSIM=0.443; p=0.2%) ( Figure 7B ). Except for March 2015, in the other spring samplings the diatom biomass dominated over the other phytoplankton groups, particularly in April 2017 when a remarkable diatom biomass was obtained on the account of Leptocylindrus spp. However, not particularly high PP rates (1.66 ± 1.08 µg C L^-1 h^-1) were estimated in this season at the 5-m depth layer. With the exception of July 2018, when the diatom biomass prevailed at this depth, in the other summer samplings the flagellates dominated the phytoplankton biomass as for instance in July 2017, when the second maximum PP rate of the study period was estimated (6.46 µg C L^-1 h^-1). In autumn, diatoms and flagellates equally contributed to the phytoplankton biomass (October15 and October17), whereas in October 2016, when a relative PP maximum was reached (4.74 µg C L^-1 h^-1) several diatom species co-occurred (Pseudo-nitzschia “seriata complex”, Chaetoceros socialis, Guinardia flaccida and others). Overall, higher PP rates were estimated during the autumn samplings (3.47 ± 0.95 µg C L^-1 h^-1). In January 2018, the phytoplankton biomass was dominated mostly by Lauderia annulata, but also Chaetoceros curvisetus and other diatom species contributed to PP, while in the other winter samplings a scarce phytoplankton biomass was present at this water depth.

Also at the 10-m layer, the phytoplankton assemblages significantly differed among seasons (R_ANOSIM = 0.432; p = 0.2%) ( Figure 7C ). The highest PP values estimated in spring at 10 m depth were observed in April 2016 (3.43 µg C L^-1 h^-1) and April 2017 (1.93 µg C L^-1 h^-1) when the phototrophic biomass was dominated by diatoms. In the former sampling, small Chaetoceros forms prevailed whereas in the latter a large Leptocylindrus species contributed to the maximum biomass obtained at this water layer. In summer, the highest PP rates were reached in July 2017 (3.87 µg C L^-1 h^-1) and July 2018 (3.37 µg C L^-1 h^-1), in the first case with a dominance of flagellates, in the second case with a major diatom biomass (Chaetoceros anastomosans, Chaetoceros spp. < 20 µm). In autumn, flagellates, diatoms, dinoflagellates and coccolithophores co-occurred at this layer and an average PP rate of 1.96 ± 0.74 µg C L^-1 h^-1 was obtained. Lauderia annulata was present also at this layer and contributed to PP in January 2018.

At the 15-m layer, the phytoplankton assemblages significantly differed among seasons (R_ANOSIM = 0.414; p = 0.1%) ( Figure 7D ). In spring 2015 and 2016, higher PP rates were estimated at this layer compared to spring 2017, although a major diatom biomass, mostly on the account of Leptocylindrus spp. persisted at the bottom layer in April 2017. A remarkable PP rate (4.46 µg C L^-1 h^-1) was measured at the bottom layer in July 2018, in correspondence to a photosynthetically active phytoplankton assemblage (mostly Leptocylindrus spp., Chaetoceros anastomosans, small flagellates). In the other summer samplings, lower PP rates were obtained (1.91 to 2.50 µg C L^-1 h^-1) when flagellates were the major contributors to the phototrophic biomass. Low PP rates, often close to zero (0.56 ± 0.38 µg C L^-1 h^-1) were estimated in autumn and winter due to light depletion at 15 m depth.

At the surface sediment, the MPB biomass did not significantly differ among seasons (R_ANOSIM = 0.208; p = 5.5%). The maximum PP rate (28.50 µg C L^-1 h^-1) was estimated in January 2018, in correspondence to a major contribution of Paralia sulcata (16.6%) and Pinnularia spp. (16.4%) to the benthic phototrophic biomass ( Figure 7E ). The second maximum was obtained in October 2016 mostly on the account of Nitzschia sigma var. sigmatella (21.4%). High PP values were further measured in April 2016 (9.42 mg C m^-2 h^-1) and April 2017 (4.12 mg C m^-2 h^-1). The highest benthic PP rates were consistently obtained in late winter-early spring or early autumn, i.e. in not excessively high temperature and light conditions. In contrast, the lowest PP rates were estimated in correspondence to the minima of MPB biomass, in October 2015 and December 2016, and in low light conditions. In each summer sampling, although a major MPB biomass was reached, the photosynthetic rates did not display high values. In June 2015, the benthic community was dominated by several Nitzschia species (N. fasciculata, N. sigma, N. lorenziana) and the large Gyrosigma acuminatum, but an uncoupling between the autotrophic biomass and PP occurred. Some of these taxa, namely N. fasciculata, Nitzschia spp., but also larger Naviculacea, seemed to be driven by higher concentrations of TOC and TN (Spearman R: Nitzschia spp. vs TN, p = 0.001; N. fasciculata vs TN, p = 0.013, and vs TOC, p = 0.007; Navicula sp.2 vs TN, p = 0.010, and vs TOC, p = 0.012). Similarly, in July 2017, the presence of very large diatom species (Tropidoneis longa, Entomoneis alata, Nitzschia dubia) and their high biomass was not mirrored in an equally high PP rate. The same pattern was observed on 30 March 2015, when the highest MPB biomass of the study period, due to the co-occurrence of taxa belonging to both winter (e.g. Pinnularia spp.) and spring (e.g. Gyrosigma acuminatum) assemblages, was uncoupled with the benthic photosynthetic rate.

In the oligotrophic Gulf of Trieste, the phytoplankton density and structure are mostly driven by the nutrient availability (Mozetič et al., 2010; Cibic et al., 2018b; Cozzi et al., 2020). This was confirmed by our results that highlighted a persistent P-limitation for the phytoplankton development regardless of the degree of the water column stability. Indeed, the nutrient ratios applied to our data revealed P-limitation indiscriminately during both stratification and mixing of the water column. This was replaced by Si-limitation during a diatom bloom, as occurred in April 2017 ( Figure S4 ), and coupled with N-limitation. N-depleted conditions were further observed in July 2018 and March 2019, when the diatom biomass prevailed over the other phytoplankton groups in the water column ( Figures 3 , S4 ). The N/P ratios higher than the Redfield one are a common feature not only for the gulf of Trieste (Lipizer et al., 2012) but also for the northern Adriatic Sea (Giani et al., 2012). However, the anomalous thermohaline conditions of the water column (extreme deviations in temperature and salinity compared to the climatological monthly means, Figure S1 ), particularly evident at the surface layer, seemed to influence greatly the phytoplankton development, as already observed in the area (Cerino et al., 2019) and in other shallow coastal waters (e.g. Trombetta et al., 2019). This clearly emerged in the PCA built on surface abiotic data, in which temperature and salinity were the predominant elements of the first factor, while in the other PCAs, dissolved inorganic nutrients were the major contributors to the total variance.

Focusing on the photosynthetic rates, we did not find a clear relationship between PPw and the stability of the water column (Spearman test p > 0.05, data not shown), neither considering the integrated rates. Indeed, only the highest PPi values of July 2017 and 2018 were obtained with a strong stratification of the water column, in typical summer conditions, whereas the other relative PPi peaks (November 2018, April and October 2016) were estimated in unstable or weak stratification conditions. During the other two summer samplings (June 2015 and July 2016), abiotic factors other than PAR availability ( Figure S2 ) likely affected the phytoplankton productivity. As a general rule, major autotrophic densities observed in the summer period were not mirrored in exceptionally high PPw rates. Indeed, phytoplankton constantly adjust their photosynthetic output to match environmental constrains and optimize their growth (Talaber et al., 2014). Moreover, the maximum PPw of the study occurred in November 2018, when a pycnocline was observed at -2.5 m depth, established after a freshwater input from the Isonzo River following the autumnal precipitations (flow rate obtained from the hydrometric height provided by: Regione Autonoma FVG, Water resources management service). This freshwater-induced pycnocline, likely triggered the not exceptionally abundant phototrophic community ( Figure S4 ), leading to this major PPw. Our PPw rates were in line with previous results from the same site (Cibic et al., 2018 b). This could suggest that the trophic conditions at this coastal site have not changed after more than a decade. Interestingly, the maximum of the study period (November 2018) was consistent with that observed in November 2006 at the surface, and the high values recorded below 5 m in summer. Our PPi rates also displayed a very similar pattern to those calculated in 2006 - 07, with the highest values recorded in July 2017 and July 2018. Further, they were similar to the values estimated in the middle of the gulf from March 2011 to March 2013 (Ingrosso et al., 2016) and the south-eastern part of the gulf in the period 2010-2011 (Talaber et al., 2018). However, our mean PPi rate was about half of that estimated by Fonda Umani et al. (2007) from January 1999 to December 2001 at the same station (C1). Our daily rates ( Table 3 ) were comparable or slightly lower than those obtained by Pugnetti et al. (2005) using the ¹⁴C uptake in the northern Adriatic Sea; lower than those estimated by Decembrini et al. (2009) in the southern Tyrrhenian Sea in December (429 mg C m^-2 d^-1) but higher than their rates in July (273 mg C m^-2 d^-1); and much higher than those measured in the Mar Piccolo of Taranto (Ionian Sea) by the same ¹⁴C technique (Cibic et al., 2016).

The maximum MPB abundances were observed in spring/early summer 2015, in correspondence to the highest light availability at the bottom ( Table S1 ), in accordance with previous results (Cibic et al., 2012). However, the highest PPs rates were obtained in late winter - early spring or early autumn, with MPB abundances close to the average value of the study period. It is widely accepted that not all the viable cells present in sediments are likely to be active at once (Cahoon and Cooke, 1992). In contrast, the lowest or negative PP rates were estimated in correspondence to the minima of MPB density, and in low light conditions. The productivity of the microalgal community inhabiting surface sediments is known to be strongly dependent upon the light availability at the bottom (Miles and Sundbäck, 2000), especially in subtidal areas. Negative PPs values were previously reported, particularly from sublittoral sites worldwide (Meyercordt and Meyer-Reil, 1999; Gillespie et al., 2000; Forster et al., 2006). Overall, our PPs values are comparable to monthly data obtained from the same site in 2003-04 (Cibic et al., 2008), and to seasonal data reported from July 2010 to July 2012 (Franzo et al., 2016) as well as to those estimated in the port of Trieste at similar depths in June 2013 (Rogelja et al., 2018). In all these studies, with a few exceptions, the same pattern in primary production was observed, with higher rates from late winter to early summer, in correspondence with low temperature and increasing light conditions, and low or negative rates from late summer to early winter. As for the water column, the results in this study suggest that the trophic state of the benthic ecosystem has remained unaltered after more than 12 years.

The nutrient ratios applied to our data revealed that the MPB community was mostly P-depleted. When minima of all nutrients occurred in the overlying water (April 2016, 2017, 2018 and March 2019), the MPB abundances were consistently slightly lower, ranging from 4.5·10⁴ to 5.7·10⁴ cells cm^-3, than the average value (6.1·10⁴ ± 3.7·10⁴ cells cm^-3). Notwithstanding, in those months the MPB community was photosynthetically active (PPs ranging from ~ 1 to 9.42 mg C m^-2h^-1). Very large diatom cells (length >100 µm: Entomoneis alata, Gyrosigma acuminatum, G. attenuatum, Pinnularia spp., Tropidoneis longa) represented a considerable fraction of the MPB biomass in these sampling months. Large diatoms have relatively larger nutrient storage vacuoles that allow them to achieve a slow but steady growth rate in pulsed-nutrient or nutrient-depleted environments and use the stored nutrients when needed (Finkel et al., 2005).

At this sublittoral site, the benthic-pelagic coupling is not so tight as it might be in an estuarine or lagoon environment. The benthic microalgae occasionally enter the pelagic domain, mostly due to sediment resuspension and, on the other hand, the settlement of planktonic forms at the seafloor can increase the microphytobenthic biodiversity and the overall benthic photosynthetic activity. The mean annual contribution of benthic diatoms to the microphytoplankton (cells > 20 µm) was estimated to be 9.9% (Cibic et al., 2018b). Pennate diatoms (mostly Nitzschia, Navicula, Pleurosigma and Diploneis) were temporarily present in the phytoplankton assemblage following strong wind events that resuspended them from the bottom. However, in October 2006 when the microphytoplankton was quite scarce in the water column, Nitzschia alone reached up to 36% of this assemblage. We further calculated the mean benthic contribution to the total eukaryotic phytoplankton (2-200 µm), presented in this study, and in this case, it accounted for less than 1%, therefore it could be considered negligible.

From the benthic point of view, a two-year study focused on the MPB biodiversity carried out in the same site, identified 8 planktonic species out of 103 diatom taxa (Cibic et al., 2007 b) and found that on average 9% of the species in the sediment were planktonic. However, they were mostly in appalling condition and only species with a hardy frustule, e.g. Pseudo-nitzschia seriata, could be recognized (Cibic et al., 2007a), and were likely still photosynthetically active.

Our PPt rate per m², obtained adding PPs to PPi rates, displayed values varying from 12.44 to 70.52 mg C m^-2 h^-1, indicating that the Gulf of Trieste, although oligotrophic (Mozetič et al., 2012; Cozzi et al., 2020) is a productive ecosystem throughout the year. When the photosynthetic rates are lower along the water column, i.e. in late winter/early spring, they are compensated, to some extent, by higher rates of an active MPB community, once sufficient light reaches the surface sediments. In contrast, in summer and autumn much higher rates occur along the water column, when the benthic primary production is inhibited by low light or supra-optimal temperatures (Guarini et al., 1997; Kirk, 2000). Of course, to some degree there is an overlapping, particularly when a late winter phytoplankton bloom is triggered by favourable conditions and an abundant and photosynthetically active MPB community has already developed at the seafloor; or in summer when the water column is very clear and high PAR irradiance is available at the sea bottom. In any case, the pelagic contribution to total PP is consistently predominant in this coastal site. Indeed, excluding negative PPs values, the mean PPs contribution to PPt was 11.3%, while it reached up to 43% in January 2018. Although these results were obtained at a 17-m deep sublittoral site, they are in line with the literature, since the microphytobenthos have been reported to contribute up to 50% of the total PP in shallow coastal systems (Perissinotto et al., 2002; Montani et al., 2003). Only a few studies, mostly carried out in estuaries, investigated simultaneously the PPw and PPs (using the ¹⁴C method), and calculated the PPs contribution to PPt. Perissinotto et al. (2003) and Anandraj et al. (2007) reported that the pelagic production was much higher than the microphytobenthic one in a river-dominated temporarily open/closed estuary in South Africa (up to 64x in the closed phase) as a direct consequence of the light availability at the bottom. van der Molen and Perissionotto (2011) found that in St. Lucia Estuary (South Africa, average depth about 2 m) the pelagic contribution to the total production varied greatly, from 2% to 100% and was below 15% in 4 out of 20 sampling dates/points. In the shallow (<12 m) semi-enclosed Mar Piccolo of Taranto (Ionian Sea) the PPs contribution to PPt was very low as it varied from 1% to 10% due to high concentration of contaminants (e.g. PCBs, PAHs, Hg) accumulated in the surface sediments (Cibic et al., 2016).

To detect which abiotic factors were most related with the occurrence and density of phytoplankton groups and the overall photosynthetic rate, an exploratory approach based on the Spearman rank correlation analysis was applied separately to samples from the four water depths and sediment. Only the significant relationships were extrapolated from the five matrices and shown in Table 4 . Salinity and temperature highly affected the abundance of the main phytoplankton groups, i.e. undetermined flagellates and dinoflagellates, at the top 5 m of the water column, whereas their influence was less marked below the 10-m layer and not relevant at the sediment surface. This was also confirmed by the PCA. In contrast, light availability, which was not a limiting factor at the water surface, did not affect the phytoplankton density and photosynthetic rate but became an important driver of species richness below the 5-m layer, and particularly at the sediment surface, where a clear relation between PAR and biodiversity emerged. Among the considered inorganic nutrients, N-NO₂ highly influenced the abundance of flagellates along the water column, and of dinoflagellates below the 5-m layer. However, neither N-NO₂ nor other nutrients showed to directly affect the MPB total numbers and overall biomass that seemed to be rather influenced by TOC and TC concentrations. Yet, the nutrient availability showed to affect the number of species at the sediment surface. Indeed, changes in the relative availability of N and Si may influence the relative as well as absolute abundance of various diatom species (Gilpin et al., 2004). Interestingly, the significance of N-NH₄ on shaping the structure of the phytoplankton community increased particularly at the lower layers of the water column (<10 m) where this nutrient is regenerated from the microbial activity, particularly in summer (Souza et al., 2011; Cossarini et al., 2012), or become available after sediment resuspension. Primary production seemed to be directly influenced by salinity and temperature at the upper 10-m layer of the water column, only. PP was also dependent upon the community richness at the first 5-m layer, but a direct relation between PP and the overall abundance emerged only for the benthic community.

In conclusion, in this study, for the first time, pelagic and benthic primary production were estimated quasi-synchronously in the northern Adriatic Sea. Although the Gulf of Trieste is persistently P-limited, it is still a productive ecosystem throughout the year. We did not find a clear relationship between PP and the water column stability; however, the absolute maximum was obtained in a weak pycnocline at the surface in autumn whereas the highest integrated rates in stratified waters in typical summer conditions. We found higher benthic PP rates from late winter to early summer, in correspondence with low temperature and increasing light conditions, and low or negative rates from late summer, when PP was inhibited by supra-optimal temperature, to early winter when PP was light limited. The benthic contribution to total PP was generally low but it increased in clear water periods and/or when pelagic PP was strongly nutrient-depleted. We found PP to be influenced by salinity and temperature at the upper 10-m layer of the water column, and mostly by the N-NO₂-availability at the lower water layers. Light availability influenced the structure of the community at the sediment surface. N-NO₂ also emerged as the driving factor of the abundance of flagellates and dinoflagellates along the water column.

On a global scale, this study represents one of the few attempts of integrating pelagic and benthic primary production using the ¹⁴C technique to quantify the ecosystem productivity. The main findings, obtained in the Gulf of Trieste, could be extended beyond the geographical limits of this particular ecosystem. Our results add on the limited database on primary production from sublittoral areas and could be used as proxy of ecosystem functioning in future studies on climate change impacts.