Benthic‐pelagic coupling and trophic relationships in northern Baltic Sea food webs

Understanding marine ecosystem structure and functioning is crucial in supporting sustainable management of natural resources and monitoring the health of marine ecosystems. The current study utilized stable isotope (SI) mixing models and trophic position models to examine energy flow, trophic relationships, and benthic‐pelagic coupling between food web components. Roughly 1900 samples from different trophic levels in the food web, collected during 2001–2010 from four northern and central sub‐basins of the Baltic Sea, were analyzed for SI ratios of carbon and nitrogen. Trophic structure of the food webs among the sub‐basins was consistent, but there were differences between the proportions of energy in different trophic levels that had originated from the benthic habitat. Mysids and amphipods served as important links between the benthic and pelagic ecosystems. Much (35–65%) of their energy originated from the benthic zone but was transferred to higher trophic levels in the pelagic food web by consumption by herring (Clupea harengus). One percent to twenty‐four percent of the energy consumption of apex seal predators (Halichoerus grypus and Pusa hispida) and predatory fish (Salmo salar) was derived from benthic zone. Diets of mysids and amphipods differed, although some overlap in their dietary niches was observed. The food web in the Gulf of Finland was more influenced by the benthic subsystem than food webs in the other sub‐basins. The baseline levels of δ13C and δ15N differed between sub‐basins of the Baltic Sea, indicating differences in the input of organic matter and nutrients to each sub‐basin.

Traditionally, pelagic and benthic compartments in aquatic food webs have been studied largely in isolation and past research has mainly focused on pelagic communities, with the benthic compartment frequently viewed as a source or sink of pelagic nutrients or energy (Schindler and Scheuerell 2002;Vadeboncoeur et al. 2002). However, several studies have shown that biotic components of aquatic ecosystems are discontinuous and coupled by movement of different organisms (e.g., Reynolds 2008). Mobile, opportunistic species in particular can physically and functionally couple these habitats (Vander Zanden et al. 2006;Reynolds 2008;Baustian et al. 2014). Research in recent decades has challenged our view concerning the pathways of energy flow in marine pelagic systems, as it is now known that they are more complex in structure and more diverse in their energy sources than previously believed (Reynolds 2008;Baustian et al. 2014).
Large pelagic marine systems are mainly driven by pelagic primary production, although in some cases, for example, areas close to river deltas, substantial amounts of energy can be of terrestrial origin, transported into the system by rivers and through the atmosphere as terrestrial particulate organic matter (POM) (Rolff and Elmgren 2000;Vähätalo et al. 2011;Woodland and Secor 2013). Part of the overall pelagic biomass and energy (hereafter "pelagic energy") is retained by the pelagic food web, but eventually most of it settles onto the seabed as fecal material and remains of dead individuals, providing energy for the benthic community. Some of this benthic biomass and energy (hereafter "benthic energy") is recycled back in to the pelagic food web through benthic-pelagic coupling mechanisms, that is, through active (organism movement, trophic interactions) or *Correspondence: mikko.j.kiljunen@jyu.fi This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
passive transport (biogeochemical cycling) (Marcus and Boero 1998;Raffaelli et al. 2003;Baustian et al. 2014). Most studies of benthic-pelagic coupling in aquatic systems have concentrated on fluxes and energy flow of organic matter and effects of nutrients on production and community structure (e.g., Karlson et al. 2007) and many of these have been conducted in lakes (e.g., Schindler and Scheuerell 2002;Vander Zanden and Vadeboncoeur 2002;Vander Zanden et al. 2006). However, relatively little is known about the direct energy transfer from the benthic to the pelagic compartment via trophic links and how this energy propagates into the pelagic food web (Raffaelli et al. 2003;Baustian et al. 2014). The role of benthic-pelagic coupling has been proposed to be particularly important in shallow marine ecosystems where a large proportion of energy from pelagic primary production settles to the benthic layer (Marcus and Boero 1998;Garstecki et al. 2000;Woodland and Secor 2013;Baustian et al. 2014;Kopp et al. 2015;Duffill Telsnig et al. 2018). One such shallow system is the Baltic Sea, which is a semi-enclosed brackish water ecosystem with rather simple food webs, containing roughly a tenth or less of the species richness found in the NE Atlantic (Johannesson and André 2006).
Benthic organisms, particularly nektobenthos (e.g., mysids and amphipods), move between the benthic, demersal, and pelagic zones by passive resuspension caused by extreme hydrological events, active ontogenetic movements, and regular diel vertical migrations (Rudstam et al. 1989;Raffaelli et al. 2003). Nektobenthos is important prey not only for demersal and epibenthic fish, but also for pelagic fishes such as herring (Clupea harengus) (Casini et al. 2004;Woodland and Secor 2013). Therefore, nektobenthos may function as a dominant vector for integration of pelagic and benthic trophic pathways (Marcus and Boero 1998;Albertsson 2004;Woodland and Secor 2013). In the open Baltic Sea, few abundant demersal fish species exist, and two of them, cod (Gadus morhua) and European flounder (Platichthys flesus), are abundant only in the southern Baltic (Nissling et al. 2002;Nielsen et al. 2013). However, pelagic herring which are obligate zooplanktivores during early life stages that are confined to coastal areas (Aro 1989) usually switch to nektobenthic prey as they grow (Rudstam et al. 1992;Casini et al. 2004). Herring has a wide salinity tolerance and is abundant virtually everywhere in the Baltic Sea. Because adult herring are capable of opportunistically shifting from pelagic to benthic habitats and are important prey for fish-eating mammals, birds, and predatory fishes, they have the potential to serve as an important energy link between benthic and pelagic compartments.
Abundance and species composition of nektobenthos varies within the Baltic Sea, which can be attributed to the spatial and temporal variations in the oxygen deficiency in deep central Baltic sub-basins, the geological history of the sea, and adaptation of species to the brackish environment (Salemaa et al. 1990;Laine 2003). In coastal food webs, the species diversity of nektobenthos is higher (Berezina et al. 2017) than in northern open-sea areas where only Mysis mixta, Mysis relicta, and Mysis salemaai (Salemaa et al. 1986;Väinölä 1986;Audzijonyte et al. 2005) and the amphipods Monoporeia affinis and Pontoporeia femorata (Laine 2003) are dominant. Hypoxia and anoxia have severe impacts on both benthic and nektobenthic species (Diaz and Rosenberg 1995;Carstensen et al. 2014) and prolonged anoxic periods decrease the abundance of nektobenthos in large areas of the Baltic Sea (Välipakka 1990). Therefore, patchiness in the availability of nektobenthos to opportunistic feeders such as adult herring may cause spatial variability in the strength of coupling between benthic and pelagic pathways. In some areas, adult herring are forced to consume more zooplankton, which may eventually decrease their energy intake and growth due to the lower energy density of zooplankton (Arrhenius and Hansson 1993). Changes in the diet strongly contribute to the fluctuations in herring growth rates (Flinkman et al. 1998;Rönkkönen et al. 2004). As herring is an important food source for aquatic predators, prey switching may cause spatial differences in the strength of coupling between benthic and pelagic pathways, and eventually influence the energy flow dynamics and trophic relationships of the brackish water food webs.
Developments in stable isotope analysis (SIA) and in modeling tools have made the isotope approach increasingly attractive in the exploration of marine food webs (e.g., Davenport and Bax 2002;Woodland and Secor 2013;Kopp et al. 2015;Giraldo et al. 2017;Stasko et al. 2018). In particular, the development of Bayesian SI models such as MixSIR, SIAR, and MixSIAR (Moore and Semmens 2008;Parnell et al. 2010;Stock et al. 2018) has allowed many simultaneous sources of uncertainties, for example, in the dietary sources or in trophic fractionation values to be propagated through the model. These models can resolve the most likely set of dietary proportions given the isotope ratios in a set of possible food sources and a set of consumers. In this study, we aim to reveal trophic interactions within and between the Baltic Sea pelagic and benthic food webs and the magnitude of energy flow from primary energy resources across the food web up to the top predators using SI ratios of carbon (δ 13 C) and nitrogen (δ 15 N) (Fig. 1). Extensive data sets from four different sub-basins of the Northern Baltic Sea are utilized to estimate the trophic position (TP) of key species in the ecosystem across different trophic levels and the relative contributions of benthic and pelagic energy sources to the pelagic food web. Such comprehensive knowledge of largescale, trophic interactions is necessary for development of improved management of the marine living resources and in rebuilding the resilience of the Baltic Sea ecosystem.

Research area and species composition
The Baltic Sea is a large brackish-water body comprising several distinct sub-basins. This study focused on food webs of the four northernmost sub-basins ( Fig. 2): the Bothnian Bay (mean depth = 42 m), the Bothnian Sea (69 m), the Gulf of Finland (38 m), and the Northern Baltic Proper (71 m). The large inflow of freshwater from the drainage area and irregular saline water inflow through the Danish Straits create a progressive south to north gradient of decreasing salinity (Lass and Matthäus 2008). Together these water masses induce strong vertical salinity stratification hence limiting the vertical mixing of the water body. These characteristics, together with anthropogenic eutrophication, have induced deep-water oxygen deficiency in large areas of the Baltic proper and at times in the Gulf of Finland. In the Bothnian Sea and the Bothnian Bay, the deep-water oxygen concentrations are higher because shallow sills limit the inflow of saline, oxygen-depleted water, and the anthropogenic nutrient loading is lower (Lass and Matthäus 2008).
The number of marine species decreases and of freshwater species increases toward the northern and eastern parts along with decreasing salinity. The zooplankton communities in the study areas are dominated by calanoid copepods, but also cladocerans and rotifers are abundant (Suikkanen et al. 2013;Kuosa et al. 2017). Three species of mysids (Crustacea) exist in the study area (M. mixta, M. relicta, and M. salemaai). However, in this study, the ecologically very similar M. relicta and M. salemaai were grouped as M. relicta (Audzijonyte et al. 2005). The most common species of soft-bottom benthic macrofauna in the area are Saduria entomon (Isopoda), Limecola balthica (Bivalvia; previously called Macoma balthica), the invasive Marenzelleria spp. (Polychaeta), and the amphipods Monoporeia affinis and Pontoporeia femorata, the latter being scarce in the Gulf of Bothnia (Laine 2003).
Marine fish species, herring and sprat (Sprattus sprattus) in particular, are the most abundant planktivorous species in the Baltic Sea pelagic ecosystem. Three-spined stickleback (Gasterosteus aculeatus) occupies both coastal and open sea areas of the Baltic sub-basins. The relative abundance of freshwater planktivorous species (vendace [Coregonus albula] and smelt [Osmerus eperlanus]) increases near estuarine areas and toward the northern and northeastern parts of the Baltic (HELCOM 2012). Riverspawning anadromous species, including Atlantic salmon (Salmo salar), support commercial and especially recreational fisheries around the Baltic Sea.
Low biodiversity of the brackish sea is also reflected in the marine mammal diversity. Only three species of pinnipeds are found: harbor seal (Phoca vitulina), gray seal (Halichoerus grypus), and ringed seal (Pusa hispida). Harbor seals occupy only the southernmost part of the Baltic Sea. Gray seals are distributed throughout the region. Ringed seals inhabit mostly the northern regions. The ringed seal population (around 20,000 individuals) inhabits three distinct breeding areas: the Gulf of Riga, the Gulf of Finland, and the Bothnian Bay in the central, eastern, and northern Baltic Sea, respectively. Most of the ringed seal population is located in the Bothnian Bay (Sundqvist et al. 2012). The gray seals (population size over 30,000 individuals) are highly mobile and exhibit long-range movements encompassing a major proportion of the Baltic Sea, although they are most abundant in the northernmost parts of the Baltic Proper (Oksanen et al. 2014).

Sampling and sample preparation
Invertebrate samples for SIA were collected between 2003 and 2010 (May-August), as part of regular monitoring cruises of the research vessels Muikku and Aranda. Zooplankton samples (herbivorous cladocera, predatory cladocera [Cercopagis pengoi] and copepods) were collected by vertical net hauls (mesh size = 100 μm) from the bottom to the surface. Benthic macrofauna (isopods, polychaetes, and bivalves) were collected with a van Veen grab (0.1 m 2 ), filtered through a 1 mm sieve, and immediately dried at 60 C for 48 h. Nektobenthos samples (mysids and amphipods) were collected by vertical net hauls (mesh size = 500 μm) from the bottom to the surface or with a benthic sledge sampler (mesh size = 500 μm). Immediately after collection, all zooplankton and nektobenthic samples were frozen at −20 C. In the laboratory, they were carefully thawed, hand-sorted, and identified into previously defined groups/species under a stereomicroscope.
In 2002 (Table 1). Fish were stored in ice at sea and later frozen (−20 C). Samples of dorsal muscle were dissected in the laboratory for SIA.
Seal muscle tissue sampling for SIA was coordinated by the Finnish Game and Fisheries Research Institute (now the Natural Resources Institute Finland) between the years 2000 and 2009. Samples from gray seals (n = 248) were collected from Finnish hunters and occasionally from fishermen's incidental by-catch hauls from the Bothnian Bay (ICES SD 31) (n = 176), the Bothnian Sea (SD 30) (n = 18), Gulf of Finland (SD 32) (n = 22), and the Northern Baltic Proper (SD 29) (n = 9) between April and November, annually. Ringed seals (n = 82) were sampled from the Bothnian Bay annually in April-May under a special permit for scientific sampling (approximately 10 individuals per year) granted by the Ministry of Forestry and Agriculture in Finland. In the laboratory, muscle samples were cut from the scapular or pelvic regions of the seals. Seal muscle tissue samples were stored in a −20 C freezer for later SIA.

Stable isotope analyses
All samples were either freeze-dried or oven-dried at 60 C to constant weight. Dried samples were ground to a fine powder using either a freezer mill (6750 SPEX; CertiPrep, Metuchen, NJ) or mortar and pestle. SIA for carbon and nitrogen was conducted by continuous-flow stable isotope-ratio mass spectrometry (CF-IRMS). All the samples, except the 26 L. balthica samples, were analyzed at the University of Jyväskylä, Finland, using a Thermo Finnigan DELTA plus Advantage continuous-flow stable isotope-ratio mass spectrometer (CF-SIRMS) coupled with a FlashEA 1112 elemental analyzer. L. balthica were analyzed in a similar manner using a PDZ Europa 20-20 CF-IRMS connected Table 1. Stable isotope ratios (SD in parentheses), trophic position (SD), and number of analyzed individuals of organisms (n) from four Baltic Sea sub-basins. to a PDZ Europa ANCA-GSL elemental analyzer at the UC Davis Stable Isotope Facility, U.S.A. Results are expressed using the standard δ notation as parts per thousand (‰) differences from the international standard. The reference materials used were internal standards of known relation to the international standards of Vienna Pee Dee belemnite (for carbon) and atmospheric N 2 (for nitrogen). Precision was better than 0.2‰, based on the standard deviation of replicates of the internal standards. Sample analysis also yielded percentage of carbon and nitrogen from which C : N ratios (by weight) were derived. Since lipids are known to be 13 C-depleted (lower δ 13 C) relative to other major tissue constituents (DeNiro and Epstein 1978;Focken and Becker 1998), carbon SI ratios were lipidnormalized (indicated here as δ 13 C 0 ) using the C : N ratio. Seals were corrected according to Ehrich et al. (2011) (Eq. 3), fish according to Kiljunen et al. (2006), and invertebrates according to Logan et al. (2008) (Eq. 1a).

Modeling energy pathways
Using SI data, we estimated the proportion of benthic energy relative to pelagic energy in each of the four Baltic Sea sub-basin pelagic food webs. Because the SI values of basal trophic level energy sources are often spatially and temporally variable in aquatic systems (Rolff 2000;Rolff and Elmgren 2000;Syväranta et al. 2006), baseline values are preferably taken from primary consumers (Vander Zanden et al. 1997), which integrate temporal variability in primary producer SI values. In this study, baseline SI values of pelagic and benthic POM were back-calculated from SI ratios of herbivorous cladocerans and L. balthica, respectively, assuming trophic enrichment of 0.5 AE 1.54‰ for δ 13 C 0 and 3.15 AE 1.31‰ for δ 15 N (Gorokhova and Hansson 1999;Sweeting et al. 2007). SIAR 4.2 (Stable Isotope Analysis in R, Parnell et al. 2010), solving linear mixing models in the R statistical computing environment, was used to determine the contribution of different prey items for each consumer in the system. This package uses data on consumer SI ratios and fits a Bayesian model to their diet proportions based on Gaussian likelihoods with Dirichlet-distributed priors of the mean. Prey δ 13 C 0 and δ 15 N ratios were corrected for trophic enrichment (Δ) using the respective fractionation factors (see above).
The number of prey types included in the mixing models should be kept reasonable (Fry 2013). Following these tenets, a maximum of four of the most important prey sources per consumer were selected for the models based on an extensive literature search (see following section "Diet of consumers", and Table 2). Averages and standard deviations of lipidnormalized carbon (δ 13 C 0 ) and nitrogen (δ 15 N) values of the prey were used in the models. A separate model was constructed for each individual consumer. SIAR reports the dietary contributions as the mean of all feasible solutions with 5 th -95 th percentiles of the distribution ranges.
We estimated in a stepwise manner, starting from primary consumers, the proportions of benthic and pelagic energy for different consumers in the food webs. The proportion of the energy originating from the benthic compartment (BC) in the diet of consumer j was estimated as where BP i is the proportional contribution of benthic energy of the i th prey component and PP i is the proportional contribution of i th prey component consumed by consumer j.

Diet of consumers
Dietary source groups in the SIAR model were benthic POM, pelagic POM, zooplankton (cladocerans, copepods/C. pengoi) isopods, amphipods, mysids, planktivorous fish, and, in some cases, other species (see Table 2). The dominant calanoid copepods and cladocerans are mainly herbivorous, feeding on different phytoplankton species (Bleiwas and Stokes 1985;Meyer-Harms and von Bodungen 1997), that is, pelagic POM in the SIAR model. C. pengoi is a carnivorous zooplankter feeding on other cladocerans and copepods (zooplankton in the SIAR model) Lehtiniemi and Gorokhova 2008). Nektobenthic mysids (M. mixta, M. relicta) are omnivorous, feeding on detritus, phyto-and zooplankton (SIAR: benthic and pelagic POM, zooplankton). Their diet switches from containing more phytoplankton to zooplankton dominating as they grow (Viherluoto et al. 2000;Lehtiniemi et al. 2002Lehtiniemi et al. , 2009. Amphipods (M. affinis and P. femorata) feed mainly on surface sediment and utilize sunken phytoplankton (SIAR: pelagic POM) and detrital organic matter, but their diet also includes bacteria, protozoa, and meiofauna (SIAR: benthic POM) (Elmgren 1978). They also eat postlarval bivalve L. balthica (Ejdung et al. 2000), when larval bivalve end their pelagic larval phase (considered as zooplankton signal in the SIAR model). Amphipods are known to rise from the bottom to swim around in the overlying water and perform vertical migrations (Donner et al. 1987). The periods of free swimming expose amphipods to fish predation (Elmgren 1978;Donner et al. 1987).
Zooplankton (SIAR: zooplankton) are abundant in the diets of herring (e.g., Peltonen et al. 2004;Ojaveer et al. 2018), but at times herring also feed on nektobenthos (SIAR: amphipods, mysids) (Casini et al. 2004). In pelagic areas, sprat and three-spined stickleback mainly utilize small zooplankton (Peltonen et al. 2004;Lankov et al. 2010;Jakubaviči ut _ e et al. 2017a,b;Ojaveer et al. 2018), that is, cladocerans and copepods, C. pengoi in the SIAR model. The diet of smelt, living in low saline areas of the Baltic Sea, resembles the diet of stickleback (Lankov et al. 2010). The diet of vendace in the Baltic is poorly known, although in freshwater its diet consists almost entirely of crustacean zoo plankton (Northgote and Hammar 2006). Herring, sprat, and three-spined stickleback dominate the salmon diet. Sprat are rare in the diet of salmon in the Bothnian Bay and  (Hansson et al. 2001). Gray seals in the Baltic Sea feed exclusively on fish (Lundström et al. 2007(Lundström et al. , 2010Kauhala et al. 2011;Scharff-Olsen et al. 2019). In total, 23 fish species have been identified in the diet of the gray seal but only a few species typically contribute substantially to the diet. Herring (SIAR: herring) always dominates the diet of gray seals in all studies, both by numbers and biomass. In addition to herring, sprat (SIAR: planktivorous fish) and Atlantic salmon (SIAR: salmon) are important prey species. Some 18 prey taxa have been identified in the diet of the ringed seal in the Baltic Sea, but again only a few species typically contribute substantially. Herring (SIAR: herring) and three-spined stickleback (SIAR: planktivorous fish) are the most frequent prey items, both by numbers and biomass (Tormosov and Rezvov 1978;Sinisalo et al. 2006;Scharff-Olsen et al. 2019). In the northern part of the Baltic Sea, vendace (SIAR: planktivorous fish) are also important prey (Suuronen and Lehtonen 2012). In addition, the ringed seal diet included a substantial crustacean component (SIAR: Isopods) (Tormosov and Rezvov 1978; Sinisalo et al. 2006).

Isotope baseline corrections
Because of area-specific differences in the baseline SI ratios, corrections are needed for direct spatial comparisons of isotope data. Therefore, we applied corrections for SI ratios for both carbon and nitrogen before conducting statistical analyses. For δ 15 N, the TP was estimated for each individual food web component using a two-source model (Post 2002): where λ is the TP of organisms used for benthic (δ 15 N benthic ) and pelagic (δ 15 N pelagic ) nitrogen baselines (in this case λ = 2 for L. balthica and cladocerans), δ 15 N consumer is nitrogen isotope ratio of the consumer organism, and Δ 15 N is enrichment of δ 15 N for one TP, which was here assigned to be 3.15 (see section "Modeling energy pathways"). δ 15 N benthic and δ 15 N pelagic were assigned based on respective mean values of L. balthica and cladocerans for each area separately. α is the proportion of nitrogen in a consumer derived from the base of the benthic compartment and was estimated using Eq. 1.
To standardize the carbon stable isotope values of nektobenthos and herring and to be able to compare them on the same scale, we calculated an index of benthic-pelagic coupling (BPCI) following the equations in Olsson et al. (2009): where δ 13 C org is the carbon SI ratio of the organism to be corrected, δ 13 C meanbase is the mean δ 13 C 0 of all baseline organisms, and δ 13 C rangebase is the range of δ 13 C 0 (δ 13 C 0 max-δ 13 C 0 min) for the same baseline organisms used as a baseline when calculating the TP (Eq. 2). The higher the BPCI value, the larger is the contribution of benthic energy sources relative to other standardized organisms in the system.

Statistical tests
A nonparametric Kruskal-Wallis H test was applied to investigate differences between sub-basins in the index of benthic-pelagic coupling (BPCI, standardized carbon stable SI ratios) and in the TP of species/groups taking part in benthicpelagic coupling (nektobenthos and herring). Sub-basinspecific differences in nektobenthos were tested in two groups (mysids and amphipods), since only one of the analyzed species occurred in all the sub-basins. Pairwise comparisons within groups were done using Dunn's post hoc tests with Bonferroni adjustment. The statistical tests were conducted with IBM SPSS Statistical package (Released version 24.0.0.1).

Carbon and nitrogen isotope ratios (δ 13 C´and δ 15 N)
There was notable variability in the baseline (Cladocera and L. balthica) carbon and nitrogen SI ratios (δ 13 C´and δ 15 N) among the different basins of the Baltic Sea. The pelagic baseline δ 13 C´(the isotopic ratio in cladocerans, i.e., in herbivorous zooplankton) constituted an increasing sequence from the Bothnian Bay, to the Bothnian Sea, to the Gulf of Finland and finally, to the Northern Baltic Proper (Table 1). A similar spatial pattern was observed in δ 13 C 0 values of the benthic baseline organism L. balthica. In all study areas, δ 13 C values tended to be higher in L. balthica than in cladocerans (Table 1). Differences in the baseline levels induced differences in the isotopic signatures in the other levels of the food webs in the different sub-basins. The benthic baseline organisms (L. balthica) also exhibited higher δ 15 N values (up to~4‰) than the pelagic baseline (cladocerans) in all study areas. A common trend was that the δ 15 N values of both baseline organisms increased progressively from north to south.
In each sub-basin, the pelagic food web clearly reflected those isotope values measured from the baseline organisms. In most cases, the δ 13 C 0 values of the studied organisms fell between those measured from two baseline organisms. A general increase δ 15 N values was observed through the food web (e.g., from planktivorous fish to top predators) ( Table 1).

Trophic position
TP estimates for different species/groups were similar among basins. As expected, there was a clear increase in TP through the pelagic food web from primary consumers, through planktivorous fish, and finally to top predators (Figs. 3 and 4a, Table 1). There was some variability between basins in TP of the most important groups involved in benthic-pelagic coupling. TP in both groups of nektobenthic invertebrates, amphipods, and mysids differed statistically between sub-basins (amphipods: H = 31.985, df = 2, p < 0.001; mysids: H = 121.427, df = 3, p < 0.001). Pairwise comparison revealed that average TP of both groups in the Bothnian Bay was significantly lower than in other areas (p ≤ 0.001). TP of mysids in the Bothnian Sea also differed significantly from those in the other basins, being slightly lower than in the Gulf of Finland and in the Northern Baltic Proper (p ≤ 0.001). We found no statistically significant differences in herring TPs between sub-basins (H = 6.379, df = 3, p = 0.095), although average TP was slightly higher in the Northern Baltic Proper.

Utilization of prey and basal sources
The dietary proportions and utilization of benthic energy in the food web differed between sub-basins ( Table 2, and Figs. 4b).
The two groups of nektobenthic invertebrates (amphipods and mysids) used a substantial amount of energy from benthic sources (Fig. 5). Statistically significant differences in BPCI of amphipods (H = 24.36, df = 2, p < 0.001) and mysids (H = 122.43, df = 3, p < 0.001) between sub-basins indicate spatial variability in basal resource utilization (Fig. 4b). In the Gulf of Finland, amphipods were significantly more dependent on energy from benthic sources than in the Bothnian Sea (p < 0.001) or the Bothnian Bay (p = 0.039, Fig. 4b). Mysids in the Gulf of Finland were also more reliant on benthic sources than those in the other sub-basins (p < 0.001, Fig. 4b).
In the Bothnian Bay, the pelagic energy pathway clearly dominated the diet of herring, while this energy pathway decreased in importance in the sequence of the Northern Baltic Proper, the Bothnian Sea, and the Gulf of Finland (Table 2, Fig. 5). Within all sub-basins, 10-50% of the total energy intake of herring originated from the benthic compartment. Mysids were the most important prey items for herring in the Northern Baltic Proper and in the Gulf of Finland, while in the Bothnian Sea amphipods and in the Bothnian Bay zooplankton were the dominant prey items (Table 2, Fig. 5). Significant differences between sub-basins in BPCI of herring (H = 189.06, df = 3, p < 0.001) also point to the existence of spatial differences in herring diet (Fig. 4b). Significantly higher herring BPCI in the Gulf of Finland and in the Bothnian Sea (p < 0.001) compared to the two other sub-basins indicate proportionally higher benthic source utilization in these areas (Fig. 4b). These results are comparable with the results obtained from the SIAR model (Fig. 5). Based on the SIAR, other pelagic fish (sprat, stickleback, smelt, and vendace) had a similar diet composition, with copepods together with the predatory cladoceran C. pengoi (if present) being somewhat more important dietary items than herbivorous cladocerans (Table 2, Fig. 5).
In all sub-basins, up to~24% of the energy consumption of the pelagic top predators (seals and salmon) originated from the benthic compartment. Sprat and three-spined stickleback dominated the salmon diet, whereas the average proportion of herring in the diet varied between 11% and 27%. In the Bothnian Bay, vendace and smelt also had a minor importance in salmon diet. The results of the mixing model suggested that the diet of ringed seals consisted not only of fish, but that benthic fauna (S. entomon) was also an important source of energy. In the Bothnian Bay, the diet of gray seals comprised mainly herring and salmon. However, in the three more southern sub-basins, one third of the gray seal diet consisted of planktivorous fish species, for example, sprat and sticklebacks (Table 2, Fig. 5).

Discussion
Our analysis showed substantial differences in the structure and functioning of the food webs between the studied subbasins in the central and northern Baltic. Mysids and amphipods were dominant vectors to coupling of benthic and pelagic systems in all the studied sub-basins. Depending of the sub-basin, 30-70% of their energy originated from benthic sources and was transferred to higher trophic levels and to the pelagic subsystems, in particular, through consumption by adult herring. Up to~24% of the energy of the top predators, gray seals, ringed seals, and salmon originated from benthic sources. The most striking difference between the sub-basins was the greater importance of benthic energy sources in the Trophic position (TP)    diets of the key organisms in the Gulf of Finland than in the other sub-basins in the northern Baltic Sea.
Although it is widely accepted that benthic-pelagic coupling often plays an important role in the functioning of aquatic ecosystems (Schindler and Scheuerell 2002;Raffaelli et al. 2003;Baustian et al. 2014), its quantification has proved difficult. Recent developments in SIA and modeling techniques have opened up new possibilities to study the pathways in energy transfer between these two habitats. However, obtaining temporally and spatially representative stable isotope data from a dynamic ecosystem such as the Baltic Sea is challenging. Individual organisms may constantly move, change their diets (ontogenetically and seasonally) and vary in their physiologies, and therefore have different isotope turnover rates. Moreover, the specification of dominant prey items for each consumers in this study is based on information from the literature and plays a central role in the model outputs. Based on literature information, we included only the most abundant prey species. In some cases, different taxa were combined to reduce the number of sources in the analyses. Thus, diets of omnivorous species (e.g. amphipods, mysids, and seals) were unavoidably simplified. To obtain data from all trophic levels and both benthic and pelagic systems, we pooled data collected over a period of 10 yr. There were also seasonal differences when benthos, zooplankton (May-August), fish (February-December), and seals (April-November) were sampled. These limitations may introduce uncertainty and noise to the data and affect interpretation of the diet and benthic-pelagic coupling estimates. Regardless of these shortcomings, we conclude that our models captured   Table 2).
the magnitude of benthic-pelagic coupling among the subbasins and provide estimates that are in good agreement with earlier stable isotope studies conducted in shallow marine ecosystems (Woodland and Secor 2013;Kopp et al. 2015;Giraldo et al. 2017;Duffill Telsnig et al. 2018). For example, Kopp et al. (2015) andDuffill Telsnig et al. (2018) showed that in the English Channel and in the northern North Sea herring obtained on average 14-35% of its energy through benthic pathways, while our results suggest that herring obtain 29% (average of all sub-basins) of energy from benthic sources.
Considering the high abundance of herring in all the studied basins (ICES 2019), adult herring is likely to be an important link between benthic and pelagic habitats in northern Baltic Sea ecosystems. Variability in reliance on benthic resources has been investigated by a few other studies (Woodland and Secor 2013;Kopp et al. 2015;Giraldo et al. 2017). In our study, the key taxa in the benthic-pelagic coupling, nektobenthos, and herring exhibited relatively large spatial variability in benthic reliance, which could be a result of variability in basin morphometry and prey availability. Depth is considered as one of the key regulators of benthic-pelagic coupling (Baustian et al. 2014) as it determines the proximity between the two compartments (Schindler and Scheuerell 2002). The basins in the current study are generally shallow, but depth varies between basins (average depth = 38-71 m). It has been shown that depth influences the resource use of many pelagic fish species (Kopp et al. 2015;Giraldo et al. 2017). Our results suggest that this is partly true also for the herring, but prey availability in general is likely to play a more important role in the Baltic Sea. Herring were more reliant on benthic energy in the Bothnian Sea and the Gulf of Finland, where 37-45% of their overall assimilated energy originated from the benthic sources. In the Bothnian Bay, however, herring rely mainly on pelagic sources probably due to less abundant nektobenthos. In the Baltic Proper, low oxygen concentration near the bottom has decreased the area suitable for nektobenthic species, and thus may limit the possibility for growing herring to shift from zooplankton to nektobenthic prey. Indeed, three decades ago Elmgren (1989) stated that oxygen deficiency had led to loss of formerly important benthic food webs in the Baltic Proper, representing an area of 100,000 km 2 .
Particularly high benthic reliance of nektobenthos in the Gulf of Finland could be a result of anthropogenic factors and general characteristics of this basin. Based on the HELCOM Eutrophication Assessment Tool, the Gulf of Finland is one of the most eutrophicated areas in the Baltic Sea (Ranft et al. 2011), mainly due to the large input of nutrients from the River Neva (Viktorsson et al. 2012). Indeed, in our study, higher δ 15 N values of the entire pelagic ecosystem in the Gulf of Finland compared to other areas point to substantial runoff from agricultural areas (McClelland et al. 1997;Voss et al. 2006). Corresponding results have been shown from rivers draining from densely populated catchments to the Southern Baltic Sea, which have δ 15 N ratio in particulate organic nitrogen twice those of a river draining to the northern part of the Baltic Sea with more pristine vegetation in the catchment area (Voss et al. 2006). On the other hand, intense blooms of N 2 -fixing cyanobacteria during the summer and autumn occur in the Gulf of Finland, which can decrease the δ 15 N in organic particles in the surface layer (Voss et al. 1997). Nevertheless, our study shows that the higher δ 15 N values of terrestrial runoff are sufficient to yield identifiable spatial patterns in the δ 15 N of the Gulf of Finland samples, even given the occurrence of cyanobacteria blooms. It is likely that in this highly eutrophicated shallow system excess nutrients of allochthonous origin end up on the bottom through sedimentation, providing a more important source for benthic secondary production relative to more oligotrophic systems.
Although nektobenthos has been shown to act as an important vector between pelagic and benthic food webs through daily migrations in the water column followed by a return to bottom refugia (Jumars 2007;Woodland and Secor 2013), these movements can exhibit considerable plasticity (Jumars 2007). Simply because of the smaller distance from surface to bottom, nektobenthos has a higher probability to be eaten by pelagic predators while not in the bottom refugia. Therefore, it is plausible that in the Gulf of Finland mysids will spend more time feeding in the bottom refugia than in the deeper sub-basins.
The Baltic Sea has experienced several regime shifts during the 20 th century, with the latest in the late 1980s (Österblom et al. 2007) after which the zooplankton community composition changed due to decreasing salinity (Viitasalo et al. 1995). This resulted in the domination of smaller freshwater copepods and cladocerans over larger-bodied marine copepods, further changing the herring diet (Flinkman et al. 1998;Rönkkönen et al. 2004). Indeed, a shift in herring diet has taken place in the Southern Baltic Sea where herring with length > 15 cm feed on nektobenthos particularly during the autumn and winter (Casini et al. 2010). Our results from the Northern Baltic Sea are in very good agreement with this observation, except in the Bothnian Bay where the large contribution of zooplankton in the diet of herring (75%) may indicate less severe competition for food resources with other fish species or shortage of nektobenthos prey.
Based on our results, benthic energy is transferred through the trophic steps all the way to the top predators in all the analyzed food webs. The studied food webs consisted of about five trophic levels with gray seals being the top consumer (see also Sinisalo et al. 2006). This is in accordance with earlier studies from brackish and marine ecosystems where marine mammals have been estimated to reach TP from 4 to 5 (Hobson and Welch 1992;Pauly et al. 1998a,b). The TP of ringed seal is lower than that of the more piscivorous and pelagic energy dependent gray seal. Although the herring typically forms the main bulk of the diet of seals in all sub-basins, gray seals also utilize fish species from higher trophic levels such as salmonids (Sinisalo et al. 2008;Lundström et al. 2010;Suuronen and Lehtonen 2012;Kauhala et al. 2016). In the Bothnian Bay, a substantial fraction of the energy intake of ringed seals is derived from the benthos. This is in accordance with earlier results from the Bothnian Bay, where Sinisalo et al. (2008) observed similar dietary differences between gray and ringed seals. It has been observed that ringed seals can specialize either on pelagic or benthic prey (Thiemann et al. 2007), for example, large benthic invertebrates such as the isopod Saduria entomon (Helle 1983;Sinisalo et al. 2006). However, the results regarding highly mobile species such as gray seals (e.g., Oksanen et al. 2014) or salmon (e.g., Torniainen et al. 2014) must be considered with caution as they may have been recently feeding in other areas in the Baltic Sea other than where the individuals were sampled. In particular, the salmon caught from the Bothnian Bay typically have been feeding in the Baltic Proper and are on their way back to their natal rivers when caught (Kallio-Nyberg et al. 1999).
The raw SI data provide intriguing clues, for example, of the contribution of allochthonous energy sources to the marine food webs in the Baltic Sea. We observed a clear progressive north to south enrichment in δ 13 C values of the organisms used to represent both the pelagic and benthic baselines, which can be explained by the amount of allochthonous organic matter in freshwater inflow in relation to the total water volume that differs greatly between basin (Lass and Matthäus 2008). This causes decreased terrestrial influence from the Bothnian Bay to the Baltic Proper (Bianchi et al. 1997;Rolff and Elmgren 2000;Hoikkala et al. 2015). This area-specific freshwater-seawater mixture creates an isotope imprint on primary consumers, which is further transferred via trophic linkages to higher consumers (Rolff and Elmgren 2000). The differences in the pelagic δ 13 C´baseline among the sub-basins indicate apparent differences in the input of organic carbon from the catchment area (Rolff and Elmgren 2000). Interestingly, the δ 13 C values were consistently around 2-3‰ higher in zooplankton in the current study compared to the estimates of Rolff and Elmgren (2000). For mysids the difference in carbon values was smaller, but still from 1 to 3‰ depending on the sub-basin. Likewise, the difference in δ 13 C values of herring in this study was about 2‰ higher in the Bothnian Bay but negligible in the Bothnian Sea and the Northern Baltic Proper compared to the estimates by Rolff and Elmgren (2000). These differences can result from the lipid correction made in this study, which tends to increase the δ 13 C values (e.g., Kiljunen et al. 2006;Logan et al. 2008;Ehrich et al. 2011). The observation that L. balthica had markedly higher δ 13 C values than pelagic organisms has also been made in the Arctic Chukchi Sea (Iken et al. 2010) where δ 13 C values of pelagic POM (containing phytoplankton, bacteria, other particulate matter) were roughly 7.4‰ lower than L. balthica.

Conclusions
Our results emphasize that benthic and pelagic food webs cannot be treated separately (e.g., Marcus and Boero 1998), and that benthic energy sources have to be accounted for when reconstructing food web and energy flow models in shallow seas. It is reasonable to assume that benthic-pelagic coupling is a highly variable phenomenon across spatial and temporal scales and between ecosystems, and benthic-pelagic coupling can be sensitive to human influence. The magnitude and variability of these biological processes are rarely assessed (Griffiths et al. 2017), mainly because quantification of complex interactions have been extremely difficult. SIA and other emerging analytical and statistical techniques can help in making quantitative comparisons of different systems in the future. The Baltic Sea is an excellent model system to test such techniques, since it has relatively simple food webs, but large spatial differences in environmental characteristics. The Baltic Sea would also benefit greatly from such information about trophic relationships in particular, as there is a plethora of environmental problems in the area to be solved (e.g., eutrophication, climate change, environmental toxins/pollutants, overfishing, and invasive species). Understanding of aquatic ecosystem structure and functioning are key elements to support sustainable management of aquatic resources, ecosystem restoration, and the monitoring of the status of the Baltic Sea and other multistressed ecosystems.