Introduction

Coastal seas are highly involved in fundamental biogeochemical processes controlling nutrient and organic matter cycling (Middelburg and Herman 2007; Bouwman et al. 2013; Carstensen et al. 2020). They receive nutrients and organic matter from both marine (e.g., primary production) and terrestrial (e.g., riverine inputs of terrestrial material, anthropogenic outfalls) origins and act as a filter between both realms (Asmala et al. 2017). In addition to climate change, human activities influence nutrient cycling through eutrophication, resulting in changes in both absolute values and ratios between nitrogen (N) and phosphorus (P). Major changes in nutrient supply are likely to affect requirements of primary producers with regard to elemental building blocks and have repercussions on processes regulating elemental homeostasis of consumers (ecological stoichiometry sensu Sterner and Elser 2002), with consequences for food web functioning and biogeochemical cycling.

Archived biological samples from environmental monitoring programs can be retrospectively analyzed for elemental composition (C, N, P) and stable isotope ratios of carbon and nitrogen (δ13C, δ15N) to study nutrient cycling and reconstruct food webs in relation to a changing environment. Carbon isotopes provide information about ultimate carbon sources for primary production (Fry and Sherr 1984). Nitrogen isotopes can be used to trace specific nitrogen sources, such as anthropogenic (e.g., Connolly et al. 2013) or diazotrophic ones (Rolff 2000; Karlson et al. 2015), and to quantify trophic position in consumers, if estimates of trophic discrimination factors are available (Vander Zanden and Rasmussen 1999, 2001; Post et al. 2002). Organisms at the base of the food web such as filter-feeding bivalves or grazing snails, with low motility and long-life span, are commonly used as proxies for isotope baselines (i.e., the ultimate C and N sources) since they integrate intra-annual variability of nutrients (e.g., Vander Zanden and Rasmussen 1999; Post 2002) over longer time spans, several months or even year(s) (e.g., Gorman et al. 2017), compared to organisms with shorter life spans such as zooplankton or phytoplankton. Perennial macrophytes can be used as isotope baselines when no primary consumers are available (e.g., Haubrock et al. 2020). However, perennial macrophytes typically reflect nutrient sources in the water column during the growth period (e.g., Savage and Elmgren 2004), hence the relatively fast turnover rates of their tissue should reflect a seasonal or shorter baseline.

In the Baltic Sea, the blue mussel Mytilus edulis trossulus species complex (Kijewski et al. 2006; Stuckas et al. 2009; hereafter referred to as Mytilus or blue mussel) and ephemeral filamentous green algae from the genus Cladophora, are highly abundant. Mytilus occur in densities of up to ~ 100 000 individuals m− 2 (Westerbom et al. 2008) and constitute, through their efficient suspension-feeding, an important link between the pelagic and benthic ecosystem, promoting nutrient cycling (Kautsky and Wallentinus 1980; Kautsky and Evans 1987; Attard et al. 2020). Cladophora is widely distributed in the Baltic Sea and mainly occurs from the surface down to 1 or 2 m depths, providing food and shelter for invertebrates (Snoeijs-Leijonmalm 2017). It is perennial, but overwinters as a small tuft attached to shallow rocky substrates. During summer it reaches its full growth and benefits from nutrient enrichment (Thybo-christesen et al. 1993). Both taxa are found throughout the Baltic Sea; Cladophora occurs from the Bothnian Bay down to the southern Baltic Proper while Mytilus distribution is restricted in the northernmost basin, the Bothnian bay, due to low salinity levels (Kautsky and Kautsky 2000). On a smaller scale, they co-occur on hard substrates at shallow depths and are often the dominant taxa based on biomass (Snoeijs-Leijonmalm 2017). Isotope composition in bivalves, including Mytilus, is commonly used to study pelagic organic matter origin (Magni et al. 2013) and its variability over time and space (Briant et al. 2018; Corman et al. 2018). Bivalves are considered suitable baselines for food web studies (e.g., Abrantes and Sheaves 2009; Willis et al. 2017) and in contaminant monitoring (Karlson and Faxneld 2021), while Cladophora is mainly used as a proxy of nutrient levels in coastal waters, in the Baltic and elsewhere (Mäkinen and Aulio 1986; Planas et al. 1996).

The pronounced latitudinal gradients of temperature, salinity and nutrients combined with eutrophication and ongoing climate change in the Baltic Sea provide an ideal study system to link environmental changes to resultant food web change and ecosystem functioning. In this system, retrospective chemical analyses (i.e., elemental and isotope composition) of archived samples help to quantify such potential changes in biota as a consequence of their changing environment. Riverine inputs of organic carbon and nutrients have increased in the recent time period, especially in the Bothnian Sea (Wikner and Andersson 2012; Asmala et al. 2019), and is expected to continue (Andersson et al. 2015a). In the Bothnian Sea, the N:P ratio of the dissolved inorganic pool is similar to the Redfield molar ratio of 16 (ca. 13, although the slight N limitation has increased in recent years; Rolff and Elfwing 2015), while it is considerably lower in the Baltic Proper (ca. 4), indicating strong N limitation (Savchuk 2018). Diazotrophic primary producers, such as some species of bloom forming cyanobacteria, can bypass this N-limitation by directly fixing N2. Satellite images of surface accumulations indicate that these blooms have increased since the 1980s (Kahru and Elmgren 2014), and this internal N loading now exceeds external loadings from rivers and atmospheric deposition in the Baltic Proper (Olofsson et al. 2020). Cyanobacterial blooms benefit from denitrification and phosphate release from hypoxic sediments, which exacerbate the N:P imbalance in a ‘vicious cycle’ (Vahtera et al. 2007). The most recent decade has also seen the regular occurrence of cyanobacterial blooms in the Bothnian Sea (Olofsson et al. 2020). Salinity and temperature are both lower in the Bothnian Sea compared to the Baltic Proper and are the primary factors affecting species distribution, including that of Cladophora and Mytilus. Predicted increase in temperature and decrease in surface salinity of the Baltic Sea (Räisänen 2017 and references therein) are hence expected to affect organisms, food webs and ecosystems (Andersson et al. 2015a; Vuorinen et al. 2015). Recent studies have shown a decrease in mussel populations over recent decades in the southern Baltic Proper, linked to increasing water temperature and changes in food quantity (Franz et al. 2019; Westerbom et al. 2019) and quality (Liénart et al. 2020). A shift in dominance from the canopy-forming perennial macrophyte Fucus towards opportunistic ephemeral Cladophora has been reported since the 1980 s in different areas of the Baltic Sea (Kraufvelin and Salovius 2004), likely linked to eutrophication (Kautsky et al. 1986; Råberg 2004; Torn et al. 2006). However, Fucus recovery has been observed recently in some areas (Rinne and Salovius-Laurén 2020). Nonetheless, higher temperature and declining salinity promote filamentous green algae (Takolander et al. 2017) while also disadvantaging Mytilus (Westerbom et al. 2019), suggesting climate change will enforce the shift towards ephemeral macrophytes and fewer mussel beds in the Baltic.

Tracing C and N origin in the Baltic Sea is complex, due to multiple interacting sources, especially in coastal areas. In general, eutrophication is associated with elevated δ15N values, for instance in Baltic Proper sediment (Voss et al. 2000). Also, sewage waters that are enriched in 15N are traceable in macrophytes (Savage and Elmgren 2004). However, the depleted 15N signal of synthetic N fertilizers used in agriculture (Bateman and Kelly 2007) can be confounded with the similar signal of diazotrophic cyanobacteria (Rolff 2000), the latter considered an indirect effect of eutrophication (Vahtera et al. 2007). Eutrophication can result in higher δ13C values in plankton and mussels due to increasing plankton biomass (Oczkowski et al. 2018) with the exception of cyanobacteria (generally low δ13C; Rolff 2000). In the Bothnian Sea, the naturally low δ15N-NO3 from pristine rivers (Voss et al. 2005) should equally be reflected in relatively depleted 15N baselines. Furthermore, the typically low δ13C both from particulate organic and dissolved inorganic terrestrial carbon (i.e., δ13CDIC) from extensive riverine input in the north (Rolff and Elmgren 2000), is similarly expected to be reflected in a low δ13C baseline in the Bothnian Sea. However, low temperatures and low light availability (the latter from the dark, humus rich water in the north) can also result in lower δ13C values for macrophytes (Wiencke and Fischer 1990; Hemminga and Mateo 1996). As an additional process, the global decrease in atmospheric δ13C linked with fossil fuel burning over the 20th century also has an influence on carbon isotope composition in the Baltic Sea (Gustafsson et al. 2015). An organism’s physiology (e.g., rapid growth, nutritional stress, reproductive stages) can also lead to substantial isotope variability in consumers (Doi et al. 2017; Gorokhova 2018). For instance, osmoregulation is an especially N-demanding process for Mytilus experiencing low-saline conditions in the Baltic Sea (Tedengren and Kautsky 1987) and likely influences 15N fractionation in mussels’ tissues, hence confounding the dietary origin of the δ15N signal (Liénart et al. 2020). To better trace ultimate N sources, the end-member δ15N signal can be measured in source amino acids (e.g., phenylalanine), which show negligible 15N fractionation during assimilation, physiological processes or trophic transfer (McClelland and Montoya 2002) compared to the bulk δ15N signal. Finally, regarding N and P elemental composition of algae and mussels, taxa-specific physiological requirements may override nutrient background in the water. Fast-growing algae such as Cladophora are under little homeostatic control compared to slow growing consumers like mussels (Smaal and Vonck 1997), and may therefore better reflect basin-specific nutrient conditions.

The general aim of this study is to test whether eutrophication and climate-related changes are mirrored in two widely distributed and abundant key taxa: the filamentous algae Cladophora and the blue mussel Mytilus. This main aim is divided into three objectives: to identify (I) potential region-specific differences and (II) potential consistency in temporal changes in elemental and isotope baselines of these two taxa and (III) link the observed year-to-year and long-term variability in the isotope and elemental data to environmental and oceanographic variables, and ultimately to population level data for each taxa.

We expect that (1) spatially, elemental and isotope composition of both taxa will reflect the latitudinal gradient in nutrients, with lower values of δ13C, δ15N, N% and P% in the more oligotrophic Bothnian Sea compared to the Baltic Proper. However, lower δ13C, δ15N, N% and P% values of both taxa in the Bothnian Sea may reflect increased influence of riverine (terrestrial) input in this basin; (Table S1 and Liénart et al. 2020); (2) temporally, the decrease in eutrophication in recent years should be reflected in decreases in δ13C, δ15N, N% and P% in our longest time series, though potential increase in riverine influence in the system is also expected to lead to decreasing δ13C but in increasing values of δ15N, N% and P% (Table S1 and Liénart et al. 2020); (3) N and P contents of Cladophora should reflect water nutrient concentrations in the water earlier in the year whereas Mytilus should have a greater homeostatic control over N and P content (and thus no link to water nutrient concentrations); (4) osmotic stress, especially in the low-saline Bothnian Sea, will influence blue mussel N requirements and hence confound δ15N values (Lienart et al. 2020). Table S1 describe the rationales for predictor selection that supports our hypotheses.

To explore these expectations, we take advantage of long-term oceanographic and land-based monitoring data and archived biological samples for the well-studied Baltic Sea. We retrospectively analyse elemental (C, N and P) and stable isotope (δ13C, δ15N) composition of Cladophora and Mytilus from three contrasting regions (coastal Bothnian Sea, coastal Baltic Proper and open sea Baltic Proper), in time series spanning 8 to 24 years and link these biotic proxies to potential drivers. Finally, we correlate chemical, environmental and oceanographic data to population level biomass data in both taxa.

Materials and methods

Study areas and sampling

This study used data from the following study regions: (1) the Bothnian Sea coastal region (Höga Kusten (HK)), (2) the Baltic Proper coastal archipelago area (Askö (A)), and (3) the Baltic Proper open sea (Gotland island (G)) (Fig. 1, details Table S2). Within the Swedish marine monitoring program of the phytobenthic community, both the blue mussels Mytilus and filamentous green algae Cladophora have been sampled once a year at different stations and archived over the periods 1993–2016 for A, 2000–2016 for G, and 2008–2016 for HK. These taxa were sampled in three replicates (three quadrates of 20 × 20 cm) collected every year along the same land-to-sea transects (different depths, three quadrates per depth) in late August/early September by SCUBA divers. Samples were oven-dried at 60 °C and stored in the dark, dry and at room temperature. From each year and region, Mytilus individuals, n = 3 (G and HK) n = 15 (A), from similar size class (length: 10.0 ± 1.5 mm, mean ± sd; Table S2) from 5 m depth at the most exposed station in each region (two stations at HK to allow n = 3 per year; Fig. 1), were selected for biometric and chemical analyses (details below). Cladophora material (a few grams dry weight homogenised material, coming from many specimens) from 1 m depth was selected for chemical analysis for three to six stations per region (Fig. 1). The rationale behind using several stations within each region for Cladophora was to avoid variability from local nutrient conditions (such as bird guano on a particular rock), whereas one station was deemed sufficient for mussels, which actively filter-feeds large amounts of water.

Fig. 1
figure 1

Study regions of the Bothnian Sea (HK: Höga Kusten) and the Baltic Proper (A: Archipelago close to Askö field station, and G: Gotland island). Stations within each region are shown in the three right panels: green triangles are Cladophora sampling stations only, blue circles are Mytilus and Cladophora sampling stations, black squares are pelagic environmental monitoring stations

Full size image

Temperature (°C), salinity, dissolved inorganic nitrogen and phosphorus (DIN, DIP, µmol L− 1), total nitrogen and phosphorus (total N, total P, µmol L− 1; i.e., includes dissolved inorganic, dissolved organic, and particulate fractions), and the phytoplankton community were measured by the Swedish national monitoring program for the nearby pelagic ecosystem: the station C3 close to Höga Kusten, station B1 close to Askö island in Stockholm’s archipelago and the open sea station BY31 north of Gotland island (Fig. 1, Table S2). For more details about sampling and analytical methods, see supplemental material Table S2, Andersson et al. (2015b) and Liénart et al. (2020). All environmental and phytoplankton community data are available at http://www.smhi.se (Marine environmental monitoring data - SHARK database). Surface terrestrial total organic carbon transport (TOC, total loadings in tons) for each watershed in Sweden is measured monthly within the freshwater monitoring program “river mouths” (chemical analyses) and is available at http://webstar.vatten.slu.se/db.html (HAVET 2015/2016).

Elemental and isotope analysis

Dry ground material from both Cladophora (pool of many specimens) and Mytilus soft tissue (individual whole mussel soft tissues, shells not included but weight of shells taken for later condition index calculations) were used for elemental carbon, nitrogen and phosphorus (% of C, N, P expressed per dry weight), bulk carbon and nitrogen isotope (δ13C, δ15N) and amino-acid-δ15N measurements (AA-δ15N). Individual specimens of Mytilus soft tissue were used for C and N elemental and isotope measurements and P% was measured individually in different specimens, since the mussels were too small for all analyses in the same individual (soft tissue was on average 2.1 ± 1.5 mg in mussels in this size range; Table S2). Phosphorus in Cladophora was only measured for 4 of the 6 stations at A (not measured from Lac and StA) and 3 of the 4 stations at HK (not measured for Dön). Amino-acid-δ15N analysis was performed on pooled Mytilus material (ca. 5 individuals per region and year).

Phosphorus analyses of individual mussels were performed at the Department of Ecology, Environment and Plant Sciences (Stockholm, Sweden) using a Alpkem SFA system to analyse phosphate (SS-EN ISO 15681-2:2018), after combustion and persulfate digestion, following Larsson et al. (2001). C and N elemental and bulk stable isotopes were measured at the Center for Physical Science and Technology (Vilnius, Lithuania) with a Flash EA 1112 Series Elemental Analyzer (Thermo Finnigan) connected to a DeltaV Advantage Isotope Ratio Mass Spectrometer (Thermo Fisher).

In order to validate ultimate nitrogen sources independently of potential physiological effect, we further compared bulk δ15N composition with a smaller data set of δ15N in source amino acids for both taxa. Amino-acid-δ15N analysis of Mytilus was conducted at the Center for Physical Science and Technology (Vilnius, Lithuania) and details are presented in supplemental material. In brief, the extraction and derivatization of amino acids were performed following the method described in Ledesma et al. (2020), then analysed by GC-C-IRMS with a Trace GC Ultra Gas Chromatograph (Thermo scientific) coupled to a Delta Advantage Isotope Ratio Mass Spectrometer via GCC III combustion interface (ThermoFinnigan). Cladophora amino-acid-δ15N, which had a lower N content, was analysed at the University of California Davis Stable Isotope Facility (UC Davis) following the method described on https://stableisotopefacility.ucdavis.edu/gcamino.html (details in suppl. material).

All isotope data are expressed using the conventional delta notation: δ13Csample or δ15Nsample = [(Rsample/Rstandard)–1], where R = 13C/12 C or 15N/14 N, with per mil deviation (‰) from international standards, Vienna Pee Dee belemnite for δ13C and atmospheric N2 for δ15N. External and internal standards were analysed as references within each batch of samples (details in suppl. material). Analytical uncertainties were < 0.15‰ and < 0.20‰ for δ13C and δ15N respectively and < 0.30‰ for AA-δ15N. The overall analytical precision was 0.9% and 0.2% for elemental C and N respectively.

Data analyses

Condition index and population data

The Condition Index (CIMytilus = dry weight soft tissue (g) / dry weight shell (g) × 100) was calculated as an indicator of mussels’ health status (e.g., Filgueira et al. 2013; Irisarri et al. 2015) on the same individual Mytilus used for stable isotope analysis (n.b. each mussel’s soft tissue and shell weight were taken separately). The N:P ratio of both taxa (N:PCladophora, N:PMytilus) was calculated to be compared to Redfield ratio in the water.

For each station of the three regions (Fig. 1), and for each year, (i) Cladophora total biomass (g dry weight per m2) from ca. 1 m depth, (ii) Mytilus total biomass (g dry weight including shells per m2) and abundance (individual per m2), as well the ratio between the two (Bm:AbMytilus, mg dry weight per individual, a proxy for average mussel size in the population), were calculated for all size classes (from juveniles ≥ 1 mm to bigger mussels > 10 mm, only a few individuals had a shell length exceeding 20 mm, Åkermark et al. in prep) at ca. 5 m depth. Calculations for each region (average of stations) were based on geometric mean for the total biomass of both taxa and the Bm:AbMytilus ratio, to avoid influence from extreme values, and on the median for Mytilus abundance.

Environmental and oceanographic data

For each station, environmental data were first integrated over 0–10 m depth (average) for each sampling date, then averaged per month (linear interpolation was performed for missing sampling dates or months). The average values per period were then calculated: the annual mean (January–August 31st, i.e., approximate date of biota sampling each year) for temperature, salinity, total N and P and the winter mean (January-February, i.e., reliable period to measure nutrient concentrations and reflect the general conditions of each region) for DIN and DIP concentrations (Table S2). The DIN:DIP ratio in the water (mol:mol) was calculated for each year based on winter concentrations. The maximum summer temperature (Tmax) was considered over the period June-August for the Baltic Proper stations and July-August for the Bothnian Sea (Table S2). The total biovolume (mm3 L− 1) of phytoplankton community (total Phyto.) integrated over 0–10 m depth (multiplied by 2 at stations B1 and BY31 to adjust for the tube sampling depth of 0–20 m compared to the 0–10 at C3, Olofsson et al. 2020) was calculated over the productive period (March–August, see details Table S2). In addition, and over the same period, the biovolume of four groups representing most of the phytoplankton community biomass (diatoms, dinoflagellates, cyanobacteria and the ciliate Mesodinium rubrum) were recalculated as an integrated sum (area under the curve (AUC) linear interpolation), following Liénart et al. (2020) (see details Table S2). A lag in the productive period was considered for the northern station (April-August). Summer AUC values, for only cyanobacteria, were also calculated to be later compared with Cladophora variables. Yearly averages of total organic carbon loadings (TOCterr, tons) were calculated from the different watersheds of each region (Table S2). The average day of the year the shift in water temperature from 8 to 10–12 °C occurred (i.e., temperature rise triggering spawning, Kautsky pers. com.) was calculated for each region (Tshift, in Julian days) and used as a proxy of climate change (with this shift supposedly occurring earlier in recent years). Finally, to capture the variations in the Ocean-Atmosphere Regime of the Baltic Sea (i.e., large-scale climatic index related to the North Atlantic Oscillation), the Baltic Sea Index (BSI) was calculated according to Lehmann et al. (2002).

Statistical analysis

Data treatment and statistical analysis were performed using R free software (R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/, version 4.0.3 (2020-10-10)), except for the DistLM models, which were performed in PERMANOVA + of PRIMER version 6 (Anderson et al. 2008). For all statistical univariate tests, critical p-value was set to 0.05.

Pre-analysis of Cladophora data

Since Cladophora were sampled at 3–6 stations within each region, local variability in elemental and isotope composition among these stations within each region was investigated prior to subsequent analyses. Unidirectional trends in time series using Mann-Kendall tests (details below) were explored for the individual stations within a region. Potential station effect over the entire time period was tested with a Kruskal-Wallis test followed by Bonferroni corrected Mann-Whitney pairwise comparison tests. Since little difference among stations was found (see results 3.1.1.), the annual mean of the 3–6 stations was computed for each variable and region for use in further analyses.

Region-specific differences and temporal trends in all data sets

To test for region-specific differences in elemental and isotope baselines (objective I) the last 9 years of each time period (only 9 years long for Höga Kusten) was compared using a Kruskal-Wallis test followed by Bonferroni corrected Mann-Whitney pairwise comparison tests.

To test for unidirectional trends in elemental and isotope baselines (objective II) as well as in the other biological and environmental variables for each region, Mann-Kendall tests were used. Autocorrelation was checked (acf(), {stats}), and if required, a modified version of the test for autocorrelated time series was applied (mmkh3lag(), {modifiedmk} package, version 1.5.0, allowing autocorrelations). To further test similarities in temporal development among regions (objective II), correlations were calculated for each unique variable of the same taxa over the same time series length. Correlations were evaluated, for % C, N, P, N:PCladophora with Spearman’s rank correlation while Pearson’s correlation were used for stable isotope and biometric data on Mytilus (δ13C, δ15N, CIMytilus). To test similarities between the two taxa within each region, correlations (as described above) were calculated for each unique variable between Mytilus and Cladophora.

Finally, in order to summarize the major patterns of variation among the different regions, principal component analyses (PCA, scaled data) were performed for (i) elemental and isotope composition of Cladophora and Mytilus (separate analyses) and (ii) environmental data. The results of these analyses are presented in Supplemental Material.

Linking biological and environmental and oceanographic data

Several statistical approaches were used in order to explore linkages between the observed year-to-year and long-term variability from chemical biomarkers to population level data with environmental and oceanographic variables (objective III).

First, to identify which of the environmental variables (11 for Cladophora, 15 for Mytilus, Table S1) are correlated to the combined dataset of elemental and isotope composition of Cladophora or Mytilus, distance-based linear models were used (DistLM, resemblance matrix created using Euclidean distances, evaluated using both forward and stepwise selection of explanatory variables based on Akaike information criterion). The final DistLM model was plotted as a redundancy analysis.

Second, to model the relationship between a single response variable in each taxon (δ15N, δ13C, N%, P%) or the population biomass (and for Mytilus, the abundance and Bm:AbMytilus) with the environmental and phytoplankton explanatory variables (11–17 predictors, the rationales for including predictors in the different models are presented Table S1, see also our general expectations in the end of introduction), partial least squares regression (PLSR) were used. PLSR analyses were performed on mean-centred and variance-standardized data, with models optimized to two components (as detailed in Lienart et al. 2020). Selection and removal of predictors were performed stepwise (i) according to the variable of importance for projection (VIP) scores values (with a cut-off of 0.8 Wold et al. (2001) extended to 0.9 when needed for more parsimonious models i.e., ≤ 7 predictors) and (ii) trying to maximize both R2Q (indicates model predictive capacity) and R2Y (explanatory capacity, analogous to coefficient of determination in regression analysis, R2X is the explained variance). Quality of the model was evaluated based on both explanatory R2Y > 0.6 and predictive R2Q > 0.4 capacities. Standardized regression coefficients were calculated for each predictor from the two-component model to ascertain the significance of environmental factors explaining and predicting the variance in response variables. PLSR analyses were performed using a modified version of the {pls} R package (v. 2.7-0, August 2018) as detailed in Liénart et al. (2020).

Finally, in order to test the direct relationship between the N:P ratio of the biota and the DIN:DIP in the water of each region over time (expectation 3), linear models were applied for each taxon (lm(), allowing the interaction region x N:P water). Variability in N:PCladophora was also investigated locally (4 stations, Fur, Str, Jus, Iss) for our longest time series only (A), using distance to the shore (in km from an arbitrary runoff point selected with regard to the agricultural landscape occupying the bay and its position relative to the stations), and compared with the local DIN:DIP in the water (inner station, sö28 and outer station sö25; 2004–2016 August only; Fig. 1).

Results

Region specific differences in elemental and isotope baselines

Comparison of Cladophora within regions

There were no significant differences in Cladophora δ13C and δ15N between stations in any of the regions with the exception of δ13C in the coastal Bothnian Sea (p = 0.05, between Bön and Sjä at Höga Kusten; Table S3, Fig. S1). There were no significant trends in any of Cladophora elemental and isotope variables over time (Table S4). Elemental composition differed only between some stations of the coastal archipelago (Askö) and the open Baltic Proper (Gotland, Table S3). Thus, mean values for Cladophora δ13C and δ15N and C, N and P % were calculated within regions, and used in subsequent analyses.

Comparisons among regions (expectation 1)

Both taxa had significantly higher δ15N values from the Baltic Proper coastal archipelago station (Askö) compared to the other regions (Figs. 2, S2, Table S5). Cladophora δ13C was lower in the coastal Bothnian Sea station (Höga Kusten) and Mytilus δ13C was higher in the open Baltic Proper station (Gotland) compared to the other regions (Figs. 2, S2, Table S5). For both taxa, N% and N:P ratio did not differ between the regions and P% was generally lower at Gotland compared to Höga Kusten (Figs. 2, S2, Table S5). The condition index of Mytilus (CIMytilus) was significantly higher for the mussels of Höga Kusten than for the Baltic Proper, and was significantly lower for Askö archipelago than for Gotland (Figs. 2, S2, Table S5). Within each region, δ13C and δ15N were never significantly correlated for the mussels nor algae (Table S6).

Fig. 2
figure 2

Temporal trends in elemental and isotope composition of Cladophora (upper panels, triangles) and Mytilus (lower panels, circles) from the 3 regions (Archipelago: black, Gotland: orange, Höga Kusten: blue). For a multivariate summary of region-specific pattern see Fig. S2

Full size image

Ultimate nitrogen sources: δ15N in amino acids for both taxa (expectation 1)

There was a significant correlation between bulk δ15N and δ15N-Phe for both taxa (Mytilus: p = 0.028, rho = 0.61; Cladophora: p = 0.008, rho = 0.83, Fig. 3). Between regions, Mytilus δ15N-Phe was lowest for Gotland in the open Baltic Proper (– 2.9 ± 1.6‰) and for Höga Kusten in the coastal Bothnian Sea (– 1.9 ± 0.4‰) and higher for the coastal archipelago of Askö (– 0.7 ± 2.0‰). Cladophora δ15N-Phe was lowest at Höga Kusten (– 0.7 ± 0.03‰), intermediate in Gotland (1.1 ± 3.5‰), and highest for Askö (4.6 ± 3.4‰, Fig. S3). Regardless of region, Mytilus had generally negative δ15N values of the source amino acid Phenylalanine (δ15N-Phe) compared to more positive values for Cladophora (Fig. 3).

Fig. 3
figure 3

Nitrogen isotope ratio of the source amino-acid phenylalanine (δ15N-Phe) for Mytilus (blue) and Cladophora (green) of the 3 regions

Full size image

Comparison between the two taxa within regions (expectation 1)

Regardless of region, the %P and N were higher in the mussels (N: 9.5 ± 1.0%, P: 0.8 ± 0.1%) than in the green algae (N: 1.9 ± 0.4%, P: 0.1 ± 0.04%), and δ13C was higher for Cladophora (-15.9 ± 1.8‰) than Mytilus (– 22.3 ± 0.8‰). δ15N was generally higher for Mytilus (5.1 ± 1.3‰), except for Cladophora from Askö, which showed high values comparable to the mussels at this site (5.8 ± 1.2‰), but was lower in the other regions (ca. 2.0 ± 1.5‰) (Fig. 2). Cladophora and Mytilus δ13C were positively correlated at Gotland (r = 0.65, p = 0.01), and for their δ15N values at Askö (r = 0.41, p = 0.04), while no significant correlation was observed at Höga Kusten (please note n = 9; Table S7).

Temporal changes in elemental and isotope baselines (expectation 2)

Over time, Mytilus δ13C and δ15N at Askö (tau = – 0.37, p = 0.01, tau = – 0.41, p = 0.005 respectively) and N% at Gotland (tau = – 0.50, p = 0.01) in the Baltic Proper were significantly decreasing, while no patterns in isotopes could be identified for the shorter time series of the Bothnian Sea station Höga Kusten (Table S8, Fig. 2). There was no significant trend for any Cladophora variables at the regional level (Table S8). Mytilus δ13C was positively correlated between the two coastal stations from different basins: Askö and Höga Kusten (r = 0.82, p = 0.01), and also between the two stations of the Baltic Proper (Gotland and Askö, r = 0.59, p = 0.02; Table S9). Mytilus N % was positively correlated between Gotland and Höga Kusten (rho = 0.78, p = 0.02; Table S9, no correlations between regions for Cladophora variables).

Environmental, oceanographic and population data

The biomass of Cladophora was on average similar at all regions (A: 7.7 ± 7.6 g m− 2, G: 6.6 ± 5.4 g m− 2, HK: 5.6 ± 4.7 g m− 2; Fig. 4, Table S8), and there was no temporal trend over time nor correlation between regions (Tables S8 and S9). Mytilus biomass was high at Askö in the coastal Baltic Proper archipelago (453 ± 163 g m− 2), about four times lower at Gotland in the open Baltic Proper (129 ± 43 g m− 2), and 100 times lower at Höga Kusten in the coastal Bothnian Sea (4 ± 3 g m− 2, Fig. 4, Table S8). Further details are presented in supplemental material (Fig. S4, Tables S5 and S8).

The station C3 in the Bothnian Sea was characterized by high total organic carbon loadings (TOCterr), low salinity and water temperature, low total nitrogen (N) and phosphorus (P), and low dissolved inorganic phosphorus (DIP) concentrations (Figs. 5, S5). Generally, the stations B1 and BY31, in the coastal and open Baltic Proper respectively, had higher salinity, total and dissolved N and P concentrations and total phytoplankton biovolumes (total Phyto.) compared to C3 (Figs. 5, S5). The station B1 was characterized by higher concentrations of dissolved inorganic N (DIN), while BY31 had more dinoflagellates, N2-fixing cyanobacteria and larger phytoplankton biovolumes (Figs. 5, S5). Over time (Fig. S5, Table S10), C3 showed a significant decrease in DIN:DIP. Both B1 and BY31 showed a significant increase in total P and significant decrease in temperature shift (Tshift), indicating that the rapid temperature rise from 8 to 10–12 °C is now occurring earlier in the year in the Baltic Proper. Total N and phytoplankton biomass, especially dinoflagellates, significantly increased at BY31. DIN, DIN:DIP and the ciliate Mesodinium rubrum significantly decreased and TOCterr significantly increased at B1. Environmental variables were rarely correlated (Table S11).

Fig. 4
figure 4

Temporal trends in Cladophora (triangles) and Mytilus (circles) population biomass (g dry weight per m2) from the 3 regions (Archipelago: black, Gotland: orange, Höga Kusten: blue). For Mytilus biomass, the left axis indicates Askö archipelago and Gotland and the axis on the right indicates Höga Kusten

Full size image
Fig. 5
figure 5

Principal component analysis (PCA) of environmental variables for B1 (black), BY31 (orange) and C3 (blue). Scaled data. Ellipses represent 50% of the data. See Fig. S5 for the raw dataset over time. Tmax: water maximum temperature in summer, DIN or DIP: dissolved inorganic nitrogen/phosphorus, TOCterr: total organic carbon from land, Tshift: water temperature shift, Diat.: diatoms, Dino.: dinoflagellates, Meso.: Mesodinium rubrum, NfixCyb.: N2-fixing cyanobacteria

Full size image

Links between elemental, isotope, population, and environmental data

Generally, elemental and isotope composition of Cladophora (R2adj. = 0.35) and Mytilus (R2adj. = 0.31) were associated in the DistLM analyses with 5 environmental variables, of which three were the same (DIN, salinity, Tmax, Fig. 6, Table S12). Higher δ13C values in Cladophora were linked to high salinities and maximum water temperatures in summer (Tmax). For Mytilus, high δ13C was linked to high dinoflagellate biovolumes. Higher δ15N values were linked to high dissolved inorganic nitrogen (DIN) in both Cladophora and Mytilus, driven by the high values in samples from the coastal station of Askö archipelago.

Fig. 6
figure 6

Redundancy analysis of (a) Cladophora and (b) Mytilus elemental and isotope data (black arrows) for the 3 studied regions (Archipelago: black, Gotland: orange, Höga Kusten: blue) constrained by environmental variables (red vectors) selected from the DistLM models based on Akaike information criterion. DIN: dissolved inorganic nitrogen, Total Phyto.: total phytoplankton bloom, Tmax: water maximum temperature in summer, Dino.: dinoflagellates, Tshift: water temperature shift. See Table S12 for model selection

Full size image

When isotope and elemental response variables were tested separately using PLSR across regions (Table 1), the δ15N models of both taxa met the Lundstedt criteria (predictive capacity R2Q > 0.4; Lundstedt et al. 1998). In line with DistLM models, high values of δ15N for both taxa were explained by high DIN. High δ15N in Cladophora was generally linked to high nutrient levels (DIP, total N) and to an earlier warming of waters (Tshift). Large blooms of N2-fixing cyanobacteria were linked to low δ15N of Mytilus. The condition index (CIMytilus) was also included as a predictor in the mussel’s δ15N model, with better condition (high CIMytilus) linked to lower δ15N. In turn, Mytilus condition was best explained by δ15N (negative relationship). It was furthermore linked to the Baltic Sea Index, a proxy for global climate change, and to earlier warming of water (Tshift), which negatively affected the mussels’ condition. Higher TOCterr was associated with higher CIMytilus. High δ13C for Cladophora was linked to high water temperature (annual and summer maximum Tmax), high phosphorus concentrations (dissolved inorganic), and high total organic carbon (TOCterr), the latter a proxy for low light penetration (Table S1). Low δ13C in Mytilus was linked to high DIN, high TOCterr, a proxy for food quality (i.e., degraded food (Attermeyer et al. 2018), bearing a low δ13C signal), and to lower salinities as well as lower total N and P. The models for elemental content for both taxa generally had low predictive and explanatory performance. High N% in Cladophora was linked to high nutrients (DIN, DIP), and to low phytoplankton bloom (more light availability, a proxy for lower competition for nutrients, Table S1). In the mussels, high N% was linked to higher CIMytilus and lower water temperatures, as well as earlier warming of water. High P% in Cladophora was linked to high summer water maximum temperature, low total phytoplankton bloom and low TOCterr loadings (higher light availability). Mytilus P % was positively linked to DIN and to proxies for diet quality (diatoms and dinoflagellates).

Table 1 Results of the PLSR models for the 6 Cladophora and 8 Mytilus response variables tested (in white: individual level, in dark grey: population level)
Full size table

The simpler linear models testing for stoichiometric relationship between the elemental ratio N:P of the biota and water between different regions did not show any effect for Mytilus. However, Cladophora N:P was explained by region-specific differences in water nutrient DIN:DIP ratio (Fig. S6). These results are supported by a local dataset from Askö archipelago, where N:PCladophora is low in the inner station (Fur) compared to high values for the outer archipelago (Iss, significantly higher than the inland stations), with this pattern also observed for the DIN:DIP in the water between these two locations (Fig. S7). Regarding population models; the Mytilus biomass model met the Lundstedt criteria (predictive capacity R2Q > 0.4) and had high explanatory power (R2Y ≥ 50%; Table 1). High Mytilus biomass was linked to high δ15N, high nutrient concentrations (DIN, DIP, total N) and to low TOCterr.

Discussion

Spatial δ13C pattern differs among taxa

The δ13C values for Cladophora followed a latitudinal pattern, with most depleted values in the north, in support of expectation 1. This pattern was mainly explained by physical variables: the low δ13C in the algae from Höga Kusten, compared with algae from the Baltic Proper, was linked to low water temperature and terrestrial loading, which was our proxy for light penetration but might also indicate dissolved inorganic carbon with a depleted endmember 13C value. Light and temperature are known to drive growth and C uptake in macrophytes, thereby influencing their δ13C values. Enriched 13C (less negative δ13C) is associated with enhanced photosynthetic activity in warm areas like the tropics, experiencing high irradiance, compared to temperate and polar regions (Wiencke and Fischer 1990; Hemminga and Mateo 1996; Stepien 2015). Cladophora δ13C was further linked to phosphorus levels in the water (DIP), which likely influence its growth condition and hence δ13C. In limnic systems, P is the limiting factor for Cladophora growth (Howell 2018 and references therein) but, to our knowledge, no studies have tested whether N or P limit Cladophora growth in the Baltic Sea.

The δ13C of Mytilus showed no latitudinal pattern, but instead differed between coastal areas and open sea, in contrast to expectation 1. Similarly depleted 13C values were observed for the two coastal stations (Askö archipelago and Höga Kusten), despite them being situated in different basins. In the statistical models, low δ13C was linked to higher eutrophication (dissolved inorganic nitrogen: DIN) and total organic carbon loadings from land (TOCterr), and to lower salinity. Terrestrial or organic particles from resuspended old sediment, both bearing a low δ13C signal (Rolff and Elmgren 2000; Voss et al. 2000), generally have a higher contribution to the pool of suspended particulate organic material in shallow areas (i.e., close to the coast) or at the vicinity of river mouths (Liénart et al. 2017). This signal can be mirrored in the tissues of filter-feeders (Lefebvre et al. 2009; Pernet et al. 2012; Briant et al. 2018). In a recent study comparing the isotope signal of benthic and pelagic food web components in different sub-basins of the Baltic Sea, Kiljunen et al. (2020) found a progressive north to south 13C enrichment pattern in both benthic and pelagic baselines, supposedly influenced by the amount of allochthonous organic material from freshwater inflows. However, it is less likely that the low δ13C signal of Mytilus from the coastal station of Askö archipelago, nearly identically to Höga Kusten, can be explained by temporally similar TOCterr inputs, since total organic carbon loadings are 10-fold lower in the archipelago area compared to Höga Kusten. Other processes are likely involved, and the significant correlation over time between δ13C in mussels from Askö archipelago and Gotland suggests a global driver influencing the pelagic ecosystem, such as the Suess effect (Gruber et al. 1999; Quay et al. 2007; Gustafsson et al. 2015, see below). Despite correlation over time, higher δ13C for Mytilus from the open sea region Gotland reflects a more marine signal, less influenced by terrestrial runoff (Rolff and Elmgren 2000). An explanation for the higher δ13C in Mytilus from Gotland could be greater eutrophication of the open Baltic Proper, since 13C enriched values in mussels may be the result of high plankton biomass (Oczkowski et al. 2018). However, since phytoplankton biomass was never selected as an important predictor in our models, this explanation is less likely. In addition, P content was also lowest in Gotland mussels (and algae from this region)).

Common δ13C temporal trends may reflect global environmental changes

Over time, our longest dataset (1993–2016) from Askö archipelago in the coastal Baltic Proper revealed a significant decrease in Mytilus δ13C, as reported in Lienart et al. (2020). A similar decrease in Mytilus δ13C was also demonstrated from the same time period in Kvädöfjärden, a coastal station 100 km south of Askö archipelago (Ek et al. 2021). However, the Kvädöfjärden time series was longer (starting in 1981), and there was no evident decrease when considering the entire series of nearly 4 decades. Still, the nearly identical temporal variation (and absolute values) in δ13C for mussels from Askö archipelago and Höga Kusten, and the positive correlation over time in δ13C for Mytilus from Kvädöfjärden, Askö archipelago and Gotland, suggest a large-scale effect on δ13C values in recent decades (in contrast to our more local expectation 2). Furthermore, even though not significant (p < 0.08 in both cases), the results suggest a positive correlation trend between mussels and algae δ13C at both Gotland and Askö archipelago. This large-scale similarity in δ13C over time may hence reflect common global or at least northern hemisphere changes not included in our statistical approach. The worldwide decrease in δ13C of atmospheric CO2 as a result of global anthropogenic activities (e.g., fossil fuels use and deforestation; referred to as the Suess effect; Gruber et al. 1999; Quay et al. 2007 and references therein), has significantly affected marine biota over recent decades (Schloesser et al. 2009). This global process is reflected in ultimate source of carbon in the Baltic Sea (Gustafsson et al. 2015) and could underlie the more recent (since 1993 and onwards) decrease in the δ13C of phytoplankton mirrored in our mussels. We found a significant 1‰ depletion in δ13C values of the mussels over the 24 year time period (see calculations in suppl. material), which is consistent with the decreases reported for other marine organisms (Druffel and Benavides 1986; Bauch et al. 2000; Schloesser et al. 2009). Similar multi-decadal decreasing trends for δ13C have been observed in bivalves (including Mytilus) from the English Channel and the Mediterranean Sea (Briant et al. 2018). A general decrease since the early 90 s has also been demonstrated from various organisms and for different trophic levels in coastal and pelagic food webs of the in the North Sea and Baltic Sea, including Herring gull (Corman et al. 2018) and Atlantic salmon (Torniainen et al. 2014).

δ15N pattern is region-specific for both taxa

The high and similar δ15N values of Cladophora and Mytilus from the coastal Baltic Proper (Askö archipelago) were both linked to the high DIN concentration in this region. Savage and Elmgren (2004) described high δ15N values in Fucus vesiculosus close to a sewage outfall north of Askö archipelago but reported that values were back to background level (δ15N of 4‰) towards the Askö island, 24 km away from the sewage outlet. The δ15N of Cladophora in our study, however, remained high around Askö island (ca. 6.5‰). This could be explained by the higher nutrient uptake and growth rate of the seasonal filamentous algae, with a high surface:volume ratio (Snoeijs-Leijonmalm 2017) compared to the slow growing perennial Fucus spp. (Wallentinus 1984), which makes Cladophora a more effective N-sink. The significant decrease in δ15N of both taxa in the Askö archipelago over time is statistically linked to the significant decrease in DIN concentration in this area, probably reflecting the general decrease in land-based nutrient loadings over recent decades as a result of eutrophication mitigation (Savage et al. 2010; Elmgren et al. 2015). An alternative explanation is that it mirrors changes in the internal denitrification potential from increased hypoxic bottom waters since the 90s resulting in less DIN (Carstensen et al. 2014). Natural features such as coastal upwellings of bottom water, which are frequent on the Swedish east coast of the Baltic Proper (Lehmann et al. 2012), bear an enriched 15N signal that can contribute to high δ15N values in macrophytes (Marconi et al. 2011; Viana and Bode 2013). Bivalves including Mytilus can also capture the high δ15N signal from eutrophication (Fukumori et al. 2008; Carmichael et al. 2012; Thibault et al. 2020; Liénart et al. 2020) or from upwellings indirectly via plankton ingestion, as observed for invertebrates in other systems (e.g., Hill and McQuaid 2008).

In the Askö archipelago, both bulk δ15N and δ15N of phenylalanine (δ15N-Phe) for Cladophora where generally high and correlated, supporting the assimilation of 15 N enriched nutrients in this area. For Mytilus, there was also a positive correlation between δ15N-bulk and δ15N-Phe, but δ15N-Phe values were much lower than for Cladophora, demonstrating differences in utilization of ultimate N sources, with negative values indicating diazotrophic N (Rolff 2000; Eglite et al. 2018). Alternatively, this may indicate different turnover time of N sources between these taxa. Indeed, low δ15N in Mytilus was explained by high amounts of N2-fixing cyanobacteria and also linked to the mussel’s condition index, indicating the role of physiology on bulk δ15N variability in the mussels, with lower δ15N indicative of better condition (Liénart et al. 2020). For Mytilus, slow N turnover (Smaal and Vonck 1997), and hence likely large 15 N fractionation, may confound this signal in bulk δ15N values. In addition, the low salinity (expectation 4) of the Baltic Sea reduces the scope for growth and affects respiration rates of the mussels (Tedengren and Kautsky 1987), thus influencing N requirements and possibly confounding δ15N signal.

Cyanobacterial N supports zooplankton and deposit-feeders during summer in the Baltic Proper (Karlson et al. 2015), but there is currently no published study reporting active feeding of Mytilus on cyanobacteria. However, a recent experimental study suggests it is a relevant food source for the mussels during summer (Liénart et al. in prep). At the open sea station of the Baltic Proper, cyanobacterial N fixation is higher than in the coast (Olofsson et al. 2020). Mytilus bulk δ15N in this region was accordingly lower than in the Askö archipelago, better resembling the low δ15N-Phe in the mussels, and also lower in algae, possibly indirectly utilizing leaked N of cyanobacterial origin. Negative δ15N-Phe for particulate organic matter in surface water of Gotland has been confirmed previously (– 1.7 to – 6‰, Eglite et al. 2018). The large year-to-year variability in Cladophora bulk δ15N, and possibly in δ15N-Phe, could also reflect the occasional influence of local upwellings at the southern tip of Gotland (Lehmann et al. 2012) bringing up more enriched 15N waters.

Finally, in the Bothnian Sea Höga Kusten, the low δ15N of Cladophora is likely related to the low δ15N signal of NO3 from pristine Nordic rivers (δ15N-NO3 of 0.6 ± 1.1‰, δ15N-PON of 2.9 ± 2.1‰; Voss et al. 2005). Similarly, Kiljunen et al. (2020) show increasing δ15N values from north to south of both pelagic and benthic baselines in the Baltic Sea, likely related to eutrophication patterns. However, increased cyanobacteria blooms in the coastal zone of the Bothnian Sea over the past decade (Andersson et al. 2015b; Olofsson et al. 2020) could additionally contribute to the low bulk δ15N and negative δ15N-Phe observed in both taxa here, especially in Mytilus.

Different nutrient turnover rates in taxa reflect system nutrient dynamics

There were no significant differences in Cladophora and Mytilus N% or N:P across regions, with only P% significantly lower for both taxa in Gotland compared to the other regions (Table S5). This is surprising since the overall Baltic Sea is classified as a ‘problem area’, with eutrophication status rated as mostly poor for the Central Baltic/Baltic Proper compared with moderate to poor for the Bothnian sea according to the HELCOM Eutrophication Assessment Tool. However, the N:P ratio of Cladophora was statistically linked to the water N:P ratio at a regional scale (simple linear model, Fig. S6), and this was supported at a local scale within the Askö archipelago (inner versus outer stations, Fig. S7), in accordance with expectation 3. In our PLSR models, elemental composition for Cladophora was mainly linked to nutrient levels and, for P%, to summer maximum water temperature and phytoplankton bloom intensity, with a larger bloom driving higher competition for nutrients and resulting in less light available. Nutrient concentrations in the water were generally higher for the Baltic Proper stations, and the DIN:DIP ratio was close to Redfield values in the Bothnian Sea (ca. 13), but N limited in the Baltic Proper (ca. 6–7, Fig. S5). The high nutrient absorption rates of Cladophora to sustain rapid growth during summer (Wallentinus 1984) can thus explain the correlation between the algal N:P ratio and nutrient levels in the water.

N:P ratio in Mytilus did not differ among regions, despite differences in nutrients and salinity (expectation 3 and 4), reflecting the relatively higher degree of homeostatic control consumers have over their elemental ratios. Mytilus has a low N turnover (Smaal and Vonck 1997), which likely explains the general absence of significant temporal trends and spatial pattern in elemental content of the mussels. A high N% in the mussels, used here as a proxy for protein content (Sterner and Elser 2002), was linked to higher condition index and lower temperature, with negative effects of higher temperatures on protein content, as expected (see Liénart et al. 2020 and references therein). Resembling the δ15N model, this highlights the importance of physiological processes for interpreting consumers’ N elemental and isotope values. Nonetheless, Mytilus P%, used here as a proxy for growth (Elser et al. 2003), was significantly lower at the open sea station (Gotland), despite the potential for highest P exposure here. Still, P content was best explained by eutrophication and diet related predictors (dinoflagellates, diatoms). This emphasizes the importance of food quality for the growth of the mussels (e.g., Bracken 2017), with implications for their role as an ecosystem P sink.

From nutrients to populations: understanding population level effects

Our data showed a rather similar biomass per m2 of Cladophora across the Baltic Sea, even though a strong year to year variability was observed within regions. There are no recent studies reporting Cladophora distribution, biomass or chemical composition in the Baltic Sea. Although Cladophora has been reported as increasing in the Baltic since the 1980s (Kraufvelin and Salovius 2004), our data did not show any significant temporal trends in Cladophora over the past two decades. That Cladophora thrive in shallow areas has been mainly discussed in relation to the eutrophication-induced decline of perennial brown macroalgae Fucus (Kautsky et al. 1986; Torn et al. 2006). The two macrophytes compete for space, light and nutrients (Kautsky et al. 1986), which could explain the rather high biomass of Cladophora in the Bothnian Sea where Fucus is absent due to low salinity. In our analyses, the biomass of Cladophora was linked to large-scale environmental conditions (i.e., Baltic Sea Index), and stable isotopes were also included as main predictors, with high δ15N reflecting eutrophication (positive correlation between δ15N and DIN in the model, Table 1). However, the biomass model had low explanatory and predictive capacities, and should be interpreted with caution. Adding predictors related to top-down control, such as grazing, and biomass of competitors for nutrients, like Fucus, and ice cover and duration would likely improve the model.

The biomass and abundance of Mytilus populations were low in the Bothnian Sea, as expected due to the low salinity of this area (i.e., affect growth, Kautsky et al. 1990; Westerbom et al. 2002). High mussel biomass was positively linked to higher nutrient concentration (DIN, DIP), with both decreasing over time as a result of eutrophication mitigation, although no corresponding decrease in total phytoplankton biovolume has been shown for this time period (opposite trend apparent at all regions, see Fig. S5 and Ek et al. 2021). However, compositional changes of phytoplankton have occurred (Fig. S5; Hjerne et al. 2019), which influence diet quality and hence the growth or survival of mussels as discussed by Liénart et al. (2020). For instance, the decrease in average mussel size (Bm:AbMytilus) in the Baltic Proper was linked to decrease in DIN but also to increase in N2-fixing cyanobacteria blooms, a food source for Mytilus whose quality is still debated. Finally, high biomass and large average mussel size were also linked to high δ15N, although δ15N was simultaneously linked to a lower condition index in individual mussels, which is contradictory. In general, the condition index was highest for the mussels at Höga Kusten. This is surprising considering the low salinity of the Baltic Proper, which is already at the Mytilus distribution area salinity limit of 4 (Snoeijs-Leijonmalm et al. 2017). One ecological explanation could be related to spawning, which is characterised by dramatic weight loss during early summer (up to 50%, Kautsky 1982). It is possible that spawning may not occur for mussels living at the salinity margin at Höga Kusten, hence resulting in a high condition index contrary to Askö archipelago mussels, which are likely now undergoing two spawning events per year (Westerbom pers. com. and field obs.). Experimental studies testing the effects of the predictors discussed here on Mytilus individual growth and condition are hence needed to provide a mechanistic understanding for observed population biomass declines.

Conclusions

Our study documents large-scale and long-term patterns in elemental and isotope composition in the ephemeral filamentous macroalgae Cladophora spp. and the suspension-feeding blue mussel Mytilus sp., two key-taxa from the Baltic Sea. We statistically link the observed patterns to environmental and oceanographical monitoring data. The δ15N of both taxa responded to regional and local drivers (mainly water nutrient concentrations), yet δ15N in the source amino acid Phenylalanine (δ15N-Phe) revealed clear differences between algae and mussels regarding ultimate N source. Diazotrophic N signal was detected in Mytilus, which was masked in bulk δ15N. Clear differences in carbon isotope baselines between the two major Baltic Sea basins were found: the δ13C followed a latitudinal gradient for the algae, while it instead differed between the coastal and open sea environment for mussels. Nonetheless, mussels from the different regions had similar temporal development in δ13C, suggesting a global driver influencing δ13C dynamics in the pelagic ecosystem. Elemental composition only differed among regions for P%, which was lower at the open sea station for both taxa. Cladophora N:P reflected regional nutrient levels in the water, while Mytilus P% was linked to food quality. This study highlights the use of a multi-taxa and combined elemental and isotope approach to quantify the effects of eutrophication and climate-related environmental changes in food webs and ecosystem functioning.