Stable Carbon and Nitrogen Isotope Composition in Suspended Particulate Matter Reflects Seasonal Dynamics of Phytoplankton Assemblages in the Gulf of Riga, Baltic Sea

The ratio of stable carbon and nitrogen isotopes in the suspended particulate matter has been widely used to study processes occurring in the marine ecosystem. At the same time, the signals provided by isotope ratios in coastal ecosystems can be difficult to interpret, due to several, often contradictory processes taking place simultaneously. In this study, we hypothesized that the carbon and nitrogen isotopic variation is predominantly affected by seasonally occurring phytoplankton species succession in the Gulf of Riga, Baltic Sea. Cyclical seasonal patterns were observed for carbon and nitrogen isotopic compositions of both SPM and phytoplankton data. Enrichment of heavy isotopes in the Gulf of Riga took place during spring phytoplankton bloom (from on average between + 7.1 and + 8.8 ‰, and between − 23.7 and − 21.9 ‰ for δ15N and δ13C, respectively) and pooled at significantly lower values (from + 3.1 to + 5.1 ‰ and from − 28.7 to − 25.1 ‰ for δ15N and δ13C, respectively) for the rest of the year. At the same time, the spatial gradient of isotope ratios was sporadic and inconclusive. The results showed that terrestrial and anthropogenic input to particulate matter is negligible from spring to autumn. Multivariate analysis revealed that the observed seasonal variability was indeed driven by variation in phytoplankton species composition. The diatoms, dinoflagellates, and the ciliate Mesodinium rubrum facilitated enrichment of 15N and 13C in spring. In contrast, atmospheric nitrogen fixation by cyanobacteria and the assimilation of their released nutrients by other organisms resulted in lower δ15N values during summer. This variability requires careful considerations for conducting food web studies in temperate coastal and estuarine environments during high phytoplankton biomass periods.


Introduction
Suspended particulate matter (SPM) is a dynamic pool of both living and non-living particles that can have a role in the functioning of food webs, nutrient and contaminant cycling, and system productivity, particularly in coastal and estuarine environments (Cresson et al. 2012;Golubkov et al. 2017;Jędruch et al. 2017;Xu et al. 2019). The amount and composition of SPM in estuarine and coastal systems are affected by various external sources like riverine inflows, coastal erosion, and atmospheric deposition, as well as internal processes like primary production and organic matter mineralization.
The opportunities presented by variable isotopic composition in different sources, like δ 13 C ratio of + 10‰ in dissolved inorganic carbonates (Wimmer et al. 2013) versus δ 13 C ratio of − 28 to − 27‰ in terrestrial SPM material (Marcelina et al. 2018;Winogradow et al. 2019), has been successfully explored to characterize extent of terrestrial influence in estuarine and continental shelf waters (McKinney et al. 2010;van de Merwe et al. 2016;Jędruch et al. 2017). Similarly, temporal 15 N depletion of SPM due to fixation of atmospheric nitrogen by diazotrophs or seasonally observable enriched SPM δ values has been used to characterize impact of internal processes (Rolff 2000;Montoya et al. 2002;Landrum et al. 2011;Marcelina et al. 2018;Winogradow et al. 2019).
At the same time, the limited number of seasonal studies of carbon and nitrogen isotope variation presents somewhat controversial evidence demonstrating from none or very limited seasonal change of δ 13 C and δ 15 N (Harmelin-Vivien et al. 2008;van de Merwe et al. 2016) to substantial seasonal fluctuations (Remeikaite-Nikiene et al. 2017;Marcelina et al. 2018). Therefore, it can be assumed that in some areas, the change in isotope ratio, caused by external sources, is offset by changes created by seasonally occurring internal processes or vice versa. However, there are only a few seasonal studies, like Remeikaite-Nikiene et al. (2017), that attempts to link seasonal changes of external sources, e.g., river runoff, with internal processes, e.g., algae blooms, to substantiate such an assumption. Therefore, the present study aims to evaluate factors and processes that control spatial and temporal changes of stable isotope ratios in the SPM pool in an estuarine like gulf, i.e., the Gulf of Riga, Baltic Sea.
In the Baltic Sea surface, SPM isotopic content is largely controlled by the presence or absence of phytoplankton that incorporates dissolved nutrients into SPM (Winogradow et al. 2019). As dissolved nutrient concentrations in the water column decrease, phytoplankton exhibits less discrimination to absorbing isotopically enriched and energetically more consuming dissolved carbon species, e.g., bicarbonate or atmospheric CO 2 (Golubkov et al. 2017).
In the case of nitrogen, the dissolved nitrogen in spring is rapidly incorporated into particulate phase by phytoplankton growth. This reveals how isotope kinetics influence the bulk particulate δ 15 N value-as the dissolved nitrogen pool is depleted, fractionation of isotopes decreases as well and the particulate δ 15 N values increase over spring bloom (Savoye et al. 2003). When nitrogen limiting conditions occur in the summer, diazotrophic cyanobacteria (most commonly Aphanizomenon sp., Nodularia spumigena, and Dolichospermum sp.) become an important dissolved nitrogen source in the surface layers of the Baltic Proper (Karlson et al. 2015). These photosynthetic bacteria can cause a significant drop in δ 15 N value by fixating isotopically depleted atmospheric nitrogen (δ 15 N = 0‰). Furthermore, other organisms can incorporate the dissolved ammonia and organic nitrogen freshly generated by cyanobacteria into SPM (Ploug et al. 2010).
This let us hypothesize that the seasonal variation of SPM carbon and nitrogen isotope ratios in the study area is predominantly affected by the seasonal succession and the biomass of phytoplankton species. We also hypothesized that the effect of freshwater runoff on isotopic composition of SPM is greatest in the area directly receiving it with fading signal strength further away in transitional waters. To test these hypotheses, we conducted a temporary and spatially resolved study of carbon and nitrogen stable isotope composition of SPM together with phytoplankton species composition, and concentrations of chlorophyll a.

Study Area
The Gulf of Riga is a relatively shallow, semi-enclosed subbasin of the Baltic Sea with the average depth of 26.2 m, water volume of 424 km 3 , and water residence time of 2 to 4 years (Yurkovskis et al. 1999;Purina et al. 2018). It is strongly influenced by freshwater runoff, since its drainage area (135,700 km 2 ) significantly exceeds the Gulf of Riga surface area (16,330 km 2 ). The average annual freshwater discharge (31 km 3 year −1 ) comprises around 7% of the volume (424 km 3 ) of the Gulf (Berzinsh 1995). The freshwater runoff is distributed unevenly with 86% of it being discharged in the south-eastern part of the Gulf of Riga, mainly from its largest tributaries, the rivers Daugava (20.4 km 3 year −1 ), Lielupe (3.54 km 3 year −1 ), and Gauja (2.33 km 3 year −1 ) (Berzinsh 1995;FAO 2016). This creates an estuary-like transitional water area of up to 15 km radius from the Daugava river mouth. Water exchange with the Baltic Proper is facilitated through the Irbe Strait in the northwest and the Suur Strait in the north. That in combination with river discharge forms a pronounced latitudinal salinity gradient, most evident in spring, when salinities of 2 g kg −1 and below can be observed in south-eastern part of the Gulf and 6-7 g kg −1 in the Irbe Strait (Stipa et al. 1999;Skudra and Lips 2017). Furthermore, the southern part of the Gulf of Riga receives municipal waste water from the Riga city wastewater treatment plant. Approximately 350,000 m 3 of treated waste water is discharged one kilometer from the shore daily (Aigars et al. 2017).

Sampling
Sampling addressing interannual and spatial variability was done from May 2015 to May 2019 from the research vessel "Salme," the LR navy hydrological vessel "Varonis," and a multiple purpose ship "Mare." Each year samples were collected in spring, summer, and autumn at 28 stations located in the Gulf of Riga and in the Daugava, Lielupe, and Gauja rivers' delta areas (Fig. 1). Seasonal variation was assessed during 2017, when sampling at stations 101A, D, G, and L ( Fig. 1) was done 1-3 times per month from March to November (except August, when sampling was done only at station 101A). SPM, primary production measurement, and phytoplankton samples in the Gulf of Riga were collected as an integrated (0-10 m) sample by a 10-m plastic hose (Ø 5 cm, collected sample volume 5 L) from the euphotic layer on the board of the research vessel. The method in detail, including quality assurance procedures, is described in HELCOM Guidelines for marine monitoring (HELCOM 2017). Water samples for nutrient concentrations were taken by Niskin type bathymeters at 0, 5, and 10 m depth horizons.
Surface layer SPM in the rivers was collected from the riverbank with a bucket attached to a telescopic shaft and filled in 5-L plastic cans, which were prewashed with 0.1 N HCl and rinsed with MilliQ grade deionized water, for shortterm storage and transportation.
Water temperature and salinity vertical profiles were measured using CTD (conductivity, temperature, and depth) water probe (SBE 19plus Sea-Cat, Sea-Bird Scientific, USA) with vertical resolution of 0.5 m. Water transparency was measured with a Secchi disc.

Analytical Procedures
The concentration of chlorophyll a (Chl a) was measured according to HELCOM protocols (HELCOM 2017). Chl a was collected as SPM on glass microfiber filters (Whatman GF/C, equivalent to 1.2 µm pore size), and afterwards Phytoplankton samples (300 mL) were fixed with acid Lugol's solution. Subsamples of 10 and 25 mL of fixed samples were settled in a sedimentation chamber for 12 h and counted according to the Uthermöl technique using an inverted microscope (Leica DMI 3000, Leica Microsystems GmbH, Germany) at × 200 and × 400 magnification. The number of counted cells in all subsamples exceeded 500 (Utermöhl 1958;Olenina et al. 2006;HELCOM 2017). The carbon content (µg C m −3 ) was calculated according to Menden-Deuer and Lessard (2000).
Nutrient concentrations were determined according to Grasshoff et al. (1983). The concentrations of ammonium (NH 4 + ) and phosphate (PO 4 3− ) were measured employing the indophenol blue and molybdenum blue methods, respectively. Nitrite (NO 2 − ) and nitrate (NO 3 − ), after reduction to nitrite in a copper coated cadmium column, were determined by nitrite reaction with an azo dye. Dissolved inorganic oxidized nitrogen (NO 2+3 ) is the sum of nitrite and nitrate. The total nitrogen (N tot ) and total phosphorus (P tot ) were analyzed as nitrate and phosphate after wet digestion with persulfate. The dissolved silicates (DSi) were determined photometrically after silicate reaction with ammonium molybdate. Quality of the analyses was checked by inter-laboratory comparison exercises in the Quality Assurance of Information for Marine Environmental Monitoring in Europe Programme (QUASIMEME), from which our laboratory took part, with Z-scores in the range of [− 2;2].

Stable Isotope Analysis
For stable isotope analysis, the water sub-samples were vacuum filtered for at least 30 min on pre-combusted (at 450 °C for 2 h) 24-mm diameter glass microfiber filters (Whatman GF/F, equivalent to 0.7 µm pore size) to collect sufficient amount of material (160-2540 mL of sampled water until the filtering rate was 1 mL s −1 ). Once the samples were filtered, 5 mL of MilliQ grade water was added to wash off residual soluble salts. All filters were air dried at room temperature and thereafter lightly folded in aluminum foil cups for storage in desiccator until further preparation and analysis.
For gravimetric SPM mass determination, pre-weighted nitrocellulose membrane filters (Millipore, 45 mm diameter, 0.45 μm pore size) were used. After filtration (330-2080 mL of sampled water until the filtering rate was 1 mL s −1 ) and drying, the nitrocellulose filter weights were measured and mass concentration γ f (mg L −1 ) and consequently mass of SPM on GF/F filter was calculated as follows: where m f+s is filter + SPM mass, m f filter mass, and V f filtered water volume through nitrocellulose filter. Due to restricted sampling (1 sample per station, no replicates) and non-homogenous distribution of SPM on GF/F filters (e.g., randomly distributed visible cyanobacteria filaments), acidification to remove inorganic carbonates was not performed as it has been reported to affect δ 15 N values (Kolasinski et al. 2008;Brodie et al. 2011). The carbonates in the Gulf of Riga form a very minor proportion of total carbon as demonstrated by Carman et al. (1996); therefore, it has been expected that carbonates would not have detectable impact on carbon isotopic ratio.
Prior to the stable isotope analyses, dried GF/F filters holding SPM were cut with a cross-section into 4 equal pieces (pseudo-subsamples). Thereafter, each pseudo-subsample was wrapped in a separate tin cup and analyzed in the Laboratory of Analytical Chemistry at Faculty of Chemistry, University of Latvia, using an elemental analyzer (EuroEA-3024, EuroVector S.p.A, Italy) coupled with a continuous flow stable isotope ratio mass spectrometer (Nu-HORIZON, Nu Instruments Ltd., UK). Isotope ratios were reported relative to Vienna Pee Dee Belemnite with a lithium carbonate anchor (VPDB-LSVEC) for δ 13 C and to atmospheric nitrogen (AIR) for δ 15 N as parts per thousands (‰): where R = 13 C/ 12 C, 15 N/ 14 N. The average value of the 4 pseudo-subsamples was used to represent sample isotopic ratio. For stable isotope ratio measurement quality control, an internal standard sample of glutamic acid and international reference material L-glutamic acid USGS-40 (Reston Stable Isotope Laboratory of the US Geological Survey, Reston, Virginia, NIST®RM 8573) were used. Repeated measurements of 121 internal standard exhibited reproducibility of 0.14‰ for δ 13 C and 0.21‰ for δ 15 N. Reference material USG-40, where stable carbon isotopic and nitrogen isotopic compositions with combined uncertainties are δ 13 C VPDB-LSVEC = − 26.39 ± 0.04‰ and δ 15 N AIR = − 4.52 ± 0.06‰ (Qi et al. 2003), was used to check accuracy of the stable isotope ratio determination. Our results for reference standard USGS-40 were δ 13 C = − 26.38 (SD = ± 0.03‰, n = 20) and δ 15 N = − 4.54 (SD = ± 0.08‰, n = 20).

Primary Production Measurements
In order to evaluate the study area productivity, the light and dark bottle oxygen technique was used (Olesen et al. 1999).
Water samples were filled in 18 transparent, calibrated glass bottles for oxygen measurements. Bottles were divided into 6 groups and wrapped in the plastic optical filters (GAM-PRODUCTS) with 3 replicates to imitate the light conditions at specific depth of euphotic layer: 100%, 66%, 53%, 23%, 13%, and aluminum folium for 0% light transmittance. Initial oxygen concentrations were fixed with Winkler reagents (1 mL manganese chloride and 1 mL alkaline iodide) before incubation. The remaining samples were placed on a rotating wheel in the incubator with an electric thermometer that controls the water temperature close to natural conditions. At the end of incubation, samples were fixed with Winkler reagents. Oxygen concentrations were determined using the iodometric method-titration with standardized sodium thiosulfate in acid solution according to ISO 5813:1983. Oxygen consumption in the dark bottles (0% light transmittance) was used as a proxy of the phytoplankton community respiration, while the other five groups were used to evaluate the daily primary production rates of the water column. Oxygen concentrations were converted to carbon units according to stoichiometry of the photosynthesis equation. Water elimination coefficient (k) was calculated for each sample from measured Secchi depth as follows: where D s is Secchi depth (m). The depth of specific light conditions (z; units expressed as meters below water surface) was calculated from: where I z is light intensity at specific depth (66%, 53%, 23%, 13%) and I o the light intensity below surface (100%). Daily water column net primary production (NPP, g C m −2 day −1 ) rates were estimated by trapezoidal integration of the data from the various light conditions. Gross primary production (GPP, g C m −2 day −1 ) was calculated by summing the NPP and the community respiration (R).

Statistical Analysis
The relationship between δ 13 C and δ 15 N and the distance from the mouth of the river Daugava (a proxy of runoff impacts) as well as salinity were analyzed by Pearson's correlations. The analysis was done per sampling period (months) when data from at least 10 different stations were obtained.
K-means clustering algorithm was used to identify seasonal differences in δ 13 C and δ 15 N values (Hartigan and Wong 1979). The k-means algorithm establishes k centroids and clusters points by assigning them to the nearest centroid so that the total intra-cluster variation is minimized. The value of k was prespecified based on an exploratory analysis and, ultimately, it was set to 2. Principal component analysis (PCA) was performed to classify the main sources of data spatial variability. Each cluster, defined by k-means algorithm, was analyzed independently to identify environmental pressures specific to the season. PCA was conducted on the correlation between the isotope values and environmental parameters (temperature, salinity, Chl a, concentration of nutrients ammonium, phosphate, nitrite + nitrate) and also abundant phytoplankton groups, separately. Statistical analysis and visualization of the results were done using software R 3.6.1 (Wickham 2009; R Core Team 2019).

Physicochemical Parameters
The dynamics of temperature and nutrient concentrations in the Gulf of Riga (Online Resource 1) followed a typical seasonal pattern for temperate coastal waters. Low temperatures and high nutrient concentrations were observed during the autumn and winter months, whereas the summer months were characterized by high temperature and low nutrient concentrations. Salinity also followed a seasonal pattern with a slight decrease during the winter. Moreover, it significantly correlated (Table 1) with the distance from the Daugava river mouth in all studied periods (r = 0.61-0.87; p < 0.01), expressing persistent spatial gradient within the Gulf of Riga.

Spatial and Temporal Variation of δ 13 C and δ 15 N in SPM
δ 13 C and δ 15 N values varied from − 30.6 to − 18.5‰ and from − 0.6 to + 12.9‰, respectively, over the entire 2015-2019 study period in the Gulf of Riga. δ 13 C showed sporadic correlation with both salinity and the distance from the region of joint inflow of the main rivers (Table 1). δ 15 N demonstrated consistently negative correlation to the distance and salinity, though with varying significance levels ( Table 1).
The carbon isotopes in SPM at station 101A appeared to reach an equilibrium with river SPM only during the convective water mixing periods from autumn to early spring coinciding with larger river runoffs and low phytoplankton GPP ( Fig. 2; Fig. 3). In correspondence with the increase of phytoplankton biomass (Fig. 2), the values of δ 13 C in the SPM of the Gulf of Riga rapidly increased during April 2017 reaching maximum by April 25 (Fig. 3). Thereafter, the slow decline of δ 13 C was observed throughout May and the summer months. The δ 13 C values in riverine SPM followed an opposite pattern that coincided with runoff intensity by demonstrating decrease of values during spring and summer. δ 15 N values of SPM in the Gulf of Riga exhibited a bimodal distribution with two maxima in spring and autumn. The spring (Fig. 3) maxima was reached on May 12 after a rapid 15 N enrichment in April. Thereafter, unlike the δ 13 C, the δ 15 N values quickly decreased, reaching a summer minimum on June 28. The broader second δ 15 N maxima was observed during August and September (Fig. 3). Although the δ 15 N in riverine SPM exhibited general increase during spring and summer, the seasonal dynamic substantially differed from that observed in the Gulf of Riga. Nitrogen isotope ratio in 101A appeared to even out with river δ 15 N values over winter; however, equilibrium was reached a month later than for δ 13 C.

May 2015
May 2017  The results of PCA showed that δ 13 C and δ 15 N values were affected by different abiotic parameters and gradients (Fig. 5). Salinity and nitrates + nitrites were the main environmental drivers (PC1; 34.4%) on overall data variability during the spring season (cluster I; Fig. 5; Online resource 2). Chl a, δ 13 C, and phosphates together contributed considerably to PC2 (cluster I; Fig. 5) implying δ 13 C relation to spring phytoplankton dynamics. The impact of seasonal changes on δ 13 C became more pronounced within cluster II data. They revealed temperature, phosphates, and δ 13 C as the main contributors to PC1 during the summer-winter period (cluster II; Fig. 5). δ 15 N had no influence on cluster I data variability, but together with the spatial (salinity) and temporal (chl a dynamics, nitrates + nitrites, ammonium) environmental gradients contributed to PC2 during the summer-winter period (cluster II; Fig. 5; Online resource 2).

Impact of Seasonal Succession of Phytoplankton on Isotopic Composition of SPM
The spring bloom in 2017 had started by late March-early April reaching maximum biomass in late April ( Fig. 2; Online Resource 3). The early spring bloom (March and April) was dominated by arctic and arctic-boreal diatoms (mainly Pauliella taeniata, Thalassiosira baltica, and Chaetoceros wighamii) while the late spring bloom (May, early June) was dominated by dinoflagellate Peridiniella catenata and mixotrophic ciliate Mesodinium rubrum. The gross primary production (GPP) (Fig. 2; Fig. S2 in Online Resource 3), mainly manifested as planktonic community respiration, started to increase with the development of phytoplankton spring bloom (April-May) reaching the first maximum in late May-early June.
The second maximum of GPP was observed in July-August, indicating intensive development of summer species. Summer phytoplankton from late June to September was characterized by a mix of different functional groups-cyanobacterium Aphanizomenon flosaquae, ciliates, chlorophytes, and cryptophytes ( Fig. 2; Online Resource 3). During the autumn after convective mixing in October-November, diatoms became more abundant again, though decreasing GPP values indicated low activity of the phytoplankton community.
PCA conducted to identify impacts of biotic (phytoplankton succession) parameters on δ 13 C and δ 15 N affirms their relation ( Fig. 6; Online resource 2). The δ 13 C, as one of the main contributors to PC1 in both seasonally distinct data pools (cluster I; cluster II), evidently is directly linked to processes in the phytoplankton community. This demonstrates a clear negative correlation with diatom biomass in the spring (cluster I) and shows a direct link to carbon mass of dinoflagellates and small-sized M. rubrum (cluster II). Although to a lesser degree, the δ 15 N shows clear correlation to specific taxa as well, contributing to PC2 for both data pools, diatoms during spring (cluster I), and large-sized M. rubrum and cyanobacteria during the summer-winter period (cluster II).

Discussion
The spatial gradients of δ 13 C and δ 15 N have been previously documented in estuarine and continental shelf waters (McKinney et al. 2010, van de Merwe et al. 2016, Jędruch et al. 2017. Furthermore, it has been argued that observed seasonal variability of δ 13 C values of SPM in coastal lagoons is driven by the seasonal dynamic of riverine inputs (Remeikaite-Nikiene et al. 2017;Marcelina et al. 2018). However, in contrast with findings in these studies, we did not establish a clear spatial gradient of δ 13 C and δ 15 N values in the Gulf of Riga. The results strongly indicate that the cause of the variability in C and N isotope fractionation are biological processes. The co-variation of δ 13 C and δ 15 N with seasonally changing abiotic factors, like temperature and nutrient concentrations, establishes the seasonal nature of δ 13 C and δ 15 N values as observed by Savoye et al. (2003), reflecting the seasonal succession of phytoplankton species.
The dynamic of phytoplankton biomass and species succession in the Gulf of Riga followed the general pattern for cold temperate coastal waters (Olli and Heiskanen 1999;Yurkovskis et al. 1999;Jurgensone et al. 2011;Gustafsson et al. 2013;Purina et al. 2018). The pronounced phytoplankton bloom resulted in enrichment of 13 C in SPM, similarly to that observed by Remeikaite-Nikiene et al. (2017), as well as in enrichment of 15 N. The δ 15 N curve observed in the river mouth stations (Fig. 3c) supports this statement as the phytoplankton biomass (expressed as chl a) in the river mouth reaches maximal values during summer months, as demonstrated by earlier study of Aigars et al. (2014), whereas the chl a values in the open part of the Gulf of Riga peak during the spring. The decrease of δ 13 C values in the riverine water during the summer were in clear contrast to the increase of δ 13 C values observed in the Gulf of Riga. This demonstrates the insignificance of allochthonous terrestrial material to phytoplankton in transitional waters of the Gulf during summer in contrast to other studies conducted in estuaries or river mouths of the Baltic Sea (Golubkov et al. 2017;Remeikaite-Nikiene et al. 2017;Marcelina et al. 2018). On the other hand, an increased proportion of allochthonous SPM in transitional waters was indicated during periods with low phytoplankton activity by analogous δ 13 C and δ 15 N values in pooled river SPM and transitional water SPM values in station 101A during winter (Fig. 3), and elevated C:N ratios before the spring bloom and in late autumn (Table S2 in Online Resource 1).
The PCA clearly indicates that capacity to affect seasonal fractionation of 13 C varies among the phytoplankton groups during active biomass growth periods. Variation in δ 13 C values is mostly driven by diatoms in early spring and dinoflagellates and 33-65 μm sized M. rubrum in late spring, while other phytoplankton groups have minor or no effect. δ 13 C values in non-vernal periods are driven by dinoflagellate and 16-32 μm sized M. rubrum blooms. The relationship between δ 15 N and phytoplankton species is more complex than that of δ 13 C. The rapid increase in δ 15 N values in March-April can be explained by the development of the early spring bloom. The dominating phytoplankton species during this period were diatoms, known to prefer the more 15 N-enriched NO 3 − rather than NH 4 + for nitrogen uptake during active growth phase (Domingues et al. 2011). After exhaustion of DSi pool (Table S2 in Online resource 1), diatoms were replaced by dinoflagellates and large-sized M. rubrum in late spring (May) further depleting the inorganic nitrogen winter pool by transferring the reminder of it from dissolved phase (majority of it in NO 3 − form) to particles (Savoye et al. 2003). The relations between δ 15 N, chl a, and salinity in the spring indicate enhanced growth in transitional waters. This may be due to increased agricultural runoffs being rapidly taken up by the fast-growing phytoplankton species before their dispersion across the Gulf of Riga. The observed decrease of δ 15 N values in the summer was most likely caused by the development of the N 2 -fixing cyanobacterium A. flosaquae. In the Baltic Sea, this species has more N 2 -fixing heterocysts in June than in late summer (Klawonn et al. 2016). The observed difference in the number of heterocysts can explain 15 N enrichment in SPM after June, while A. flosaquae still maintains a significant proportion in phytoplankton biomass.
However, the negative correlation between cyanobacteria biomass and δ 15 N values is not obvious in summer according to PCA (Fig. 6, cluster II) most likely because the cyanobacteria were never the overwhelmingly dominating taxa in the non-vernal phytoplankton community during our study. δ 15 N values by the end of spring bloom in April only slightly exceeded the range of values that are typical for deep ocean inorganic nitrogen (Jędruch et al. 2017;Pantoja et al. 2002) or Baltic Proper surface SPM δ 15 N values (Winogradow et al. 2019) suggesting very limited impact of the river runoff, which delivers inorganic nitrogen of terrestrial and anthropogenic origin (see isoscapes of DIN, total nitrogen, and δ 15 N over the study period in Online resource 4, visualized with Ocean Data View v5.5.1. (Schlitzer 2021)). This is in line with the budget calculations presented by Müller-Karulis and Aigars (2011), which demonstrate that on average, 81% of annual nitrogen input was denitrified and 3% was buried. Therefore, only a minor fraction of annual input with clear anthropogenic δ 15 N signal is being retained in the water column.
Consequently, it could be expected that in the absence of a significant external nitrogen pool by the end of spring bloom, the δ 15 N value in SPM would be determined by the recycling of an already assimilated nitrogen pool. Contrary to this assumption, the highest δ 15 N value was observed in May, when no further phytoplankton carbon biomass increase could be detected (due to lack of dissolved phosphates and silicates; Table S2 in Online Resource 1). The most plausible explanation is the shift in phytoplankton species composition, e.g., successive diatom replacement by ciliates M. rubrum in May. It is highly possible that the fast vertical migration ability of this ciliate manifests δ 15 N value increase from utilizing the bottom layer nitrogen pool , as can be observed in the δ 15 N values transitioning from April 25 to May 23. Furthermore, in early June, M. rubrum did not increase the bulk δ 15 N values despite its large abundance in the phytoplankton biomass, likely due to depletion of deeper dissolved nitrogen pools. Additionally, it is possible that unfiltered water, as sampled in our study, contains other taxa, e.g., small-sized zooplankton that could increase the overall δ 15 N values of bulk SPM as shown by Rolff (2000).

Conclusions
The Gulf of Riga receives substantial riverine discharges, which are enriched by nitrogen originating from agricultural lands and municipal wastewater plants located in the drainage basin, as well as direct inputs from municipal wastewaters from Riga and Jurmala cities. The combined discharges form a strong and geographically distinctly located source of nitrogen and carbon of anthropogenic and terrestrial origin. Nevertheless, the δ 13 C or δ 15 N values of combined river SPM input signal can only be detected in the stations under direct impact of the river inflow over winter. Results of the present study demonstrate that neither terrestrial or anthropogenic sources or even a combination of them has sufficient strength of its δ 13 C or δ 15 N signal to create distinct spatial gradient of δ 13 C or δ 15 N values of SPM in the Gulf of Riga during the productive period (spring-autumn). We established that the main driver for δ 13 C or δ 15 N variability in SPM during this period is the succession of phytoplankton species. It evidently defines carbon and nitrogen isotopic ratios of SPM during the spring bloom and substantially affects isotope ratios until mid-autumn. Diatoms, dinoflagellates, and M. rubrum show the strongest positive relation to isotopic changes in the Gulf of Riga, whereas the diazotrophic cyanobacteria have an observable but a statistically insignificant negative effect on δ 15 N values. The results of this study emphasize the need to consider seasonal and spatial variations of the phytoplankton community, while assessing the spatial or seasonal variability of carbon and nitrogen isotope ratios of the SPM. In addition, the indirect evidence provided by this study let us hypothesize that the ciliate M. rubrum facilitates δ 15 N increase by incorporating near-bottom nitrogen before vertical migration to euphotic zone. Fig. 6 The results of principal component analysis for cluster I (a, b) and cluster II (c, d) in combination with carbon mass of abundant phytoplankton groups: diatoms (Diat), dinoflagellates (Dino), ciliates Mesodinium rubrum with cell size between 33 and 65 μm (Cili_33-65) and cell size between 16 and 33 μm (Cili_16-33), cyanobacteria (Cyano), and other taxa (other). Variation explained (var.expl.) represented by eigenvalues of the axes (a, c) and biplot (b, d) of individual data scores-points-and loadings of explanatory variables-arrows