Temporal and spatial changes in benthic invertebrate trophic networks along a taxonomic richness gradient

Abstract Species interactions underlie most ecosystem functions and are important for understanding ecosystem changes. Representing one type of species interaction, trophic networks were constructed from biodiversity monitoring data and known trophic links to assess how ecosystems have changed over time. The Baltic Sea is subject to many anthropogenic pressures, and low species diversity makes it an ideal candidate for determining how pressures change food webs. In this study, we used benthic monitoring data for 20 years (1980–1989 and 2010–2019) from the Swedish coast of the Baltic Sea and Skagerrak to investigate changes in benthic invertebrate trophic interactions. We constructed food webs and calculated fundamental food web metrics evaluating network horizontal and vertical diversity, as well as stability that were compared over space and time. Our results show that the west coast of Sweden (Skagerrak) suffered a reduction in benthic invertebrate biodiversity by 32% between the 1980s and 2010s, and that the number of links, generality of predators, and vulnerability of prey have been significantly reduced. The other basins (Bothnian Sea, Baltic Proper, and Bornholm Basin) do not show any significant changes in species richness or consistent significant trends in any food web metrics investigated, demonstrating resilience at a lower species diversity. The decreased complexity of the Skagerrak food webs indicates vulnerability to further perturbations and pressures should be limited as much as possible to ensure continued ecosystem functions.

vulnerability, and generality, have increased the knowledge of how different pressures affect different parts of the food web disproportionately (Gibert, 2019;.

Anthropogenic pressures result in changes in species richness
and composition through not only environmental filtering but also through biological interactions. For example, eutrophication not only directly influences the primary producers but also indirectly other trophic levels through bottom-up trophic cascades that often result in simplified communities (Jochum et al., 2012;Shurin et al., 2012), which can be visualized through increased trophic network connectance, or proportion of possible links realized. On the other hand, overfishing reduces top predators, resulting in a top-down trophic cascade (Gilarranz et al., 2016;Myers & Worm, 2003;Pace et al., 1999;Pauly, 1998), which can increase middle trophic level biomass through predation pressure release, and overgraze primary producers (Elmgren et al., 2015;Pace et al., 1999).
Such trophic network changes can be visualized through decreased vulnerability, or mean number of predators per prey, and shorter trophic chains. Climate warming has been correlated with increased connectance and increased variability in generality, or mean number of prey species per predator . In addition, pressures can act disproportionately on taxa, such as habitat destruction impacting sessile and long-lived species more than mobile shortlived ones (Bradshaw et al., 2012). Thus, when external pressures act on one organism, indirect changes to the trophic network can be predicted through trophic cascades (Shurin & Seabloom, 2005). For these reasons, trophic networks are ideal candidates for monitoring and management (Gray et al., 2014).
However, there are challenges to integrating trophic networks into management applications. Long-term data observations with high taxonomic resolution at multiple trophic levels are commonly lacking over larger spatial scales to characterize trends in trophic networks. In addition, food webs are dynamic, and links vary over spatial and temporal scales, which often results in a necessary tradeoff between resolution (i.e., all species and links are known) and scale (McMeans et al., 2015;Thompson et al., 2012). Importantly, food web structure is not always easily inferred based on species richness or community composition alone (Frelat et al., 2022), advocating for detailed assessments of trophic interaction structure to complement those of taxonomic community structure that is commonly the focus in monitoring.
Benthic invertebrate species composition and trophic interactions vary with environmental conditions and provide a range of ecosystem functions (Johannesson & André, 2006;van der Putten et al., 2004;Winder & Schindler, 2004). Detritus dominates the diet of many benthic invertebrates, and this processing of detritus contributes to nutrient recycling and mineralization of organic matter in interactions with microbes (Dossena et al., 2012;Lauringson et al., 2014;Moore et al., 2004). Detritivores dominate soft sediments, one of the largest habitats in the world, covering 80% of the ocean floor, but are one of the least studied ecosystems in terms of trophic networks (Lenihan & Micheli, 2001;Nybakken & Bertness, 2005). For these important ecosystem components, much remains unclear about how large-scale or longterm changes in benthic communities are manifested in ecological interaction structure and what the potential consequences are for functioning at the base of the food web.
The Baltic Sea has strong salinity, temperature, and species richness gradients over a relatively small geographic distance (Leppäranta & Myrberg, 2009;Snoeijs-Leijonmalm et al., 2017), and anthropogenic pressures vary as well, both by spatial and temporal scales. The effects of the interactions between these environmental gradients and pressures can be well studied through existing research, monitoring, and governmental structures in the Baltic Sea (Reusch et al., 2018). The physical and biological gradients found in the Baltic (e.g., salinity, temperature, and species richness) allow for a better understanding of impacts of various anthropogenic pressures such as eutrophication and climate change on aquatic ecosystems (Reusch et al., 2018). Pressures like climate change are particularly relevant in the Baltic Sea, where sea surface temperatures are rising three times faster than the global mean sea temperature increase (>1°C per decade; Belkin, 2009;Reusch et al., 2018). As a result, we can expect direct and indirect impacts to food webs, shown through smaller organism body size, behavioral changes such as shifted habitat preference (poleward (Kortsch et al., 2015) and vertical (John & Post, 2021) migrations) and decreased vulnerability, i.e. the number of predator species per prey (Snickars et al., 2015). Additionally, the Baltic Sea is particularly prone to eutrophication due to its long water retention time and large historical inputs of nitrogen and phosphorus (Andersen et al., 2017), which have been found to simplify benthic trophic networks by reducing number of taxa, links, and chain length, while increasing connectance O'Gorman et al., 2012).
While anthropogenic pressures are known to affect biodiversity and ecosystem functions in the Baltic Sea, it is unknown how these changes are cascading through the ecosystem and which functions could be altered by additional pressures. Anthropogenic pressures are likely to affect the Baltic Sea basins disproportionately, affecting benthic communities, trophic networks, and ecosystem functions in varying magnitudes (Griffiths et al., 2017). Benthic trophic networks, with important ecosystem-based management implications, have not been continuously monitored and we do not know how their structural changes have led to functional changes in the Baltic Sea (EU MSFD, 2008;Olivier et al., 2019;Rogers et al., 2010;Tam et al., 2017). The low species richness of the Baltic Sea have been well studied (Snoeijs-Leijonmalm et al., 2017), and feeding interactions thoroughly documented (see Janas et al., 2017, and references therein), with the exception of the fully marine Skagerrak (where the Baltic Sea meets the North Sea), where biodiversity and feeding links among the taxa there have not been investigated to the same degree.
In this study, we aim to evaluate changes in the benthic invertebrate food webs of the Swedish coast, utilizing long-term monitoring data and literature dietary links to model highly resolved food webs in a stressed ecosystem. We explore how environmental gradients and combined anthropogenic pressures have affected the benthic trophic network architecture, assessing changes in a broad suite of well-established, complementary metrics, as well as through changes in species richness and community composition during the study period. In particular, we expect food web structure to vary along the salinity gradient of the Baltic Sea, with higher salinities correlating with higher trophic network complexity due to the cooccurring species richness gradient. We also predict that changes in food web structure over time will be more evident on the Swedish west coast than for other basins due to large species loss during the study period (Obst et al., 2018).

| Study area
The Swedish national benthic invertebrate monitoring program has been continuously sampling Baltic Sea benthic invertebrates since the 1980s once per year in May. We chose to investigate food webs in 1980-1989, when benthic monitoring was conducted over a large geographical region using standardized methods for the first time, and 2010-2019, as the most recent decade of gathered data. The 26 included stations cover most of the Swedish coast of the Baltic Sea, from the almost freshwater conditions in the north of the Bothnian Sea (average salinity 5.8) to the marine conditions of the Skagerrak (average salinity 34), where the Baltic Sea meets the North Sea ( Figure 1). Benthic organisms were sampled with a van Veen sediment grab (0.1 m 2 ) and sieved through 1 mm mesh (Broman et al., 2019). When more than one replicate was taken, organism abundance was averaged over the replicates to reduce heterogeneity between sampling effort in the basins. Organisms were preserved, sorted, counted, and weighted in the lab according to the European standard, and identified to the lowest possible taxonomic classification (82% to species level, 11% to genus level, and 7% to family or order level; hereafter referred to as "species" for simplicity ;EN 16665:2014EN 16665: , 1992. Taxa composition and abundance data were retrieved from the Swedish Meteorological and Hydrological Institute's database for marine monitoring data, Sharkweb (sharkweb.smhi.se). Data from Sharkweb (the Baltic Proper, Bornholm Basin, and Skagerrak) were complemented with additional regional monitoring data from the Bothnian Sea collected by Umeå Marine Sciences Centre using the same methods previously described.  Table S1.

| Literature review, food web construction, and metrics evaluated
A literature review of feeding ecology of all taxa found in our dataset was conducted to determine their respective trophic links. Feeding links were only used when peer-reviewed sources showed feeding either from direct observation, feeding trials, or gut content analysis (Nordström et al., 2015;. References span the years 1937 to 2020, and Baltic Sea area studies were preferentially selected when possible, but when those were not available, references from Europe were taken. We did not impose a minimum number of individuals used to determine the diet of a given predator. Feeding link matrices and literature sources can be found in the Dryad Digital Repository: https://doi.org/10.5061/dryad.bk3j9 kdfk. We included two basal resources (phytoplankton and detritus) as nodes, although we recognize that selection within these resources is occurring. The benthic invertebrate monitoring program unfortunately does not record information about composition or biomass of phytoplankton and detritus. Five rare taxa were excluded from the analysis as no feeding information could be found or inferred from the available literature (the gastropod family Diaphanidae, gastropod species Mangelia attenuate and Taranis moerchii, and the annelid species Drilonereis filum and Glyphonesione klatti). Additionally, if there was uncertainty that the taxonomic resolution would interfere with the integrity of the trophic network construction, then the node was removed (e.g., all recordings of the class "Hydrozoa" were removed from Skagerrak databases due to the large diet diversity of Hydrozoa, but species-level observations were kept).
A metaweb was assembled for each basin of the study, consisting of all taxa present in all years and all possible feeding links ( Figure S2). From each basin metaweb, individual food webs were constructed for each station and year. While each basin has its own metaweb, we did not adjust the metawebs for spatial or temporal variation, meaning that links were assumed to be consistent across basins and decades when species co-occurred. This approach is conservative in estimating species interactions based on species composition changes and removes allowances for rewiring, as no in situ documentation could be provided. The package "igraph" in R (Csardi & Nepusz, 2006;R Core Team, 2021) was used for food web assembly and metric calculations, along with additional functions (trophic level and omnivory) from Kortsch et al. (2021). To identify changes in food web structure, we selected 12 complementary and widely used food web metrics in order to describe whole-network properties, including vertical and horizontal dimensions of the food web across strong gradients (taxonomic richness and salinity) and time.

| Statistical analysis
Data visualization was conducted in R using the "ggplot2" and "ggmap" packages (Kahle & Wickham, 2013;Wickham, 2016). All statistical tests were run in R (v4.1.2; R Core Team, 2021). Permutational multivariate analysis of variance (PERMANOVA) tests to detect differences in food web metrics were run using basin, decade, and station nested within basin as independent variables and each food web metric as a response variable separately. The adonis2 function in the F I G U R E 1 Map of paired benthic monitoring stations used in this study. Stations in blue were present in both decades, in pink only in the 1980s, and in green only in the 2010s. Note that the stations present in only the 1980s (pink) and only 2010s (green) were then paired for the analysis based on similar physical characteristics (see Methods) "vegan" package was used, with 999 permutations utilizing Euclidean distance, and pairwise t-tests were used to determine differences between basins and decade using the Holm adjustment method for multiple comparisons (Oksanen et al., 2019). Unpaired Welch ttests were conducted between the decades for the paired stations after checking for assumptions. P values were adjusted for multiple comparisons by the Holm adjustment method. Non-metric multidimensional scaling (NMDS) was run using Bray-Curtis similarity on community composition for the different basins and decades in the "vegan" package. An additional PERMANOVA using the same function to 999 permutations, but Bray-Curtis distance, was run to test significant differences in abundances of community composition with decade, basin, and their interaction as independent variables.
A SIMPER (similarity percentages) analysis was run using the "vegan" package within each basin to detect significant drivers of community differences. Finally, in order to visualize correlation between the TA B L E 1 Food web metrics utilized in this study, formulas, definitions, and ecological implications Gives an estimate of how connected species are, on average, within a food web Proportion of realized interactions from the total possible interactions  Indicates robustness or resistance to change Mean number of prey per predator (Schoener, 1989) Relative bias toward generalist or specialist predators; lower values indicate more specialist predators Mean number of predators per prey (Schoener, 1989) Relative risk of becoming prey; lower values indicate less predators and h is the intermediate species between species i and j Average of the number of links between any two given species in the network  Higher values indicate longer paths and less efficient energy transfer Mean shortest path (P)

the intermediate species between species i and j
Average of the shortest path between any two given species in the network (Kortsch et al., 2021) Indicative of stability, shorter paths are more stable than long ones Mean trophic level (T) Short-weighted trophic levels, combination between shortest trophic level and prey-averaged trophic level (Kortsch et al., 2021) Related to complexity, but also more trophic levels result in less efficient energy transfer Proportion of species in the network that feed on multiple trophic levels (Kortsch et al., 2021) Omnivores are stabilizing to food webs due to their high linkage and ability to prey switch

Number of motifs (M)
Total number of threespecies connected subgraphs in the network (Bascompte & Melián, 2005) More motifs indicate a more complex network; structure underling basic interactions Number of more interconnected groups of species in the network (Clauset et al., 2004;Grilli et al., 2016) Compares if groups ("modules") are more or less connected than random aggregations of species; indicates structure of network, with more modules being more stable food web metrics investigated here, a correlation analysis was run with "corrplot" ( Figure S3; Wei & Simko, 2021). A principal component analysis (PCA) was also run to identify the main spatial (basin) and temporal (decade) dynamics in food web structure, as well as assess complementarity among metrics. The package "ggbiplot" was utilized to visualize the PCA (Vu, 2011).

| Basin metawebs
The Bothnian Sea dataset had 19 species total (all years) and increased in species richness from 12 in the 1980s to 16 in the 2010s.
There were three unique species in the 1980s (16%) and seven unique

| Basin comparison in food web metrics
There were significant differences between the four basins for all  Table 2). Additional metrics can be found in Figure S4 and Table S3.

| Temporal changes
Together with basins, decade was also a significant explanatory factor for all food web metrics, with the exception of trophic level, indicating that a significant shift in food web complexity has oc-  Table 2), which indicates that the effect of decade was not consistent throughout all the basins. Additional metrics can be found in Figure S4 and Table S3.
The PCA also revealed a different clustering between the two decades for the Skagerrak region, suggesting that this basin changed the most between 1980s and 2010s, mainly driven by reduction in species richness, number of links, linkage density, generality, vulnerability, and number of motifs (Figure 3).

| Station-level changes
Results from the unpaired Welch t-tests on the station level between the decades are listed in Table S3 and visualized in Figure S5. In the TA B L E 2 Permutational multivariate analysis of variance (PERMANOVA) model results with basin, decade, station nested within basin, and an interaction effect between basin and decade as independent variables and food web metrics for the response variable, tested individually

| Macrofauna community composition
PERMANOVA results showed significant differences in macrofauna community composition between basins (F 3,0.31 = 70.8, p = .001), with the NMDS indicating clustering based on our defined basins (NMDS, stress = 0.102; Figure S6). In particular, the Skagerrak com-  Our results indicate that food webs in the higher salinity basin Skagerrak were more affected than the less saline basins, and it is likely that species adapted to the physiological stresses of the Baltic Sea brackish water environment are also resilient to additional correlated pressures, such as pH fluctuations and high seasonality (Rousi et al., 2019). We recognize that salinity is correlated with a number of other factors in the Baltic Sea (species richness and temperature), and the positive network complexity-salinity relationship seen here will not be necessarily true in all systems. For example, in the Mediterranean Sea, pelagic species richness is inversely correlated with salinity, with consequences for pelagic food web

| Anthropogenic pressures and food web changes
The species richness loss found in the Skagerrak is comparable to previously reported 51.9% reduction in species richness over 88 years from the same area (Obst et al., 2018). Obst et al. (2018) found similar numbers of species (their 254 to our 244) in the 2000s, but found a total of 607 species were present between 1920 and 1940, considerably higher than the 336 species from the 1980s found here. As sampling effort rarefaction curves indicate undersaturation in species richness in the earlier samplings of 1920-1940 (Obst et al., 2018), it is possible that macrozoobenthos decrease in species richness and impact in food web metrics is in fact even higher than what is indicated by our data. The current study differs from previous studies in that we investigate how species loss impacts trophic networks.
Previous studies proposed bottom trawling fishing and increase in turbidity as main drivers of macrozoobenthos biodiversity loss in the Skagerrak (Obst et al., 2018). Most bottom trawling occurs in the Skagerrak and Bornholm Basin for demersal fish species (mostly cod Gadus morhua; ICES, 2019), although there is a prawn fishery (Pandalus borealis) in the Skagerrak, but this species was not present in our dataset (Linders et al., 2018). Around 80 to 100% of the seabed is effected by bottom trawling in the Skagerrak and Bornholm Basin (HELCOM, 2018a), despite declines in cod catch per unit effort since the peak in the 1980s (Möllmann et al., 2021). In addition to physical disturbance, bottom trawling can also resuspend contaminants, and large stores of persistent organic pollutants exist in the soft sediments of the Baltic Sea (Jonsson, 2000). Such contaminants disproportionately affect larger and long-lived species (Bradshaw et al., 2012). Predators are generally larger and longer lived than their prey (Nordström et al., 2015), and we would expect to see re-  (Bock & Miller, 1995;Sañé et al., 2013), but there was an increase in proportion of suspension feeders in the Skagerrak from the 1980s to the 2010s (from 13.1% to 14.3%, but overall decrease in species richness (from 43 to 34 species)), demonstrating that suspension feeders were not more impacted than other trophic groups during the study period. As such, bottom trawling is not likely to be the sole explanation for reduced species richness in the Skagerrak basin.
In addition to decreased species richness, we observed an in- declines at all stations except one ( Figure S5i). Vertical network compression has been reported in Baltic Sea (Kortsch et al., 2021) and North Sea (Frelat et al., 2022) trophic networks previously.
Vulnerability decreased in the Bornholm Basin and Skagerrak, indicating prey are less vulnerable to predation by other benthic invertebrates. One explanation is an increase in specialist predators; another is that loss of prey items is driving decreased generality in the Bornholm Basin and Skagerrak. Conversely, the Bothnian Sea showed increased vulnerability, contrary to other basins, while generality decreased. We hypothesize that this increased vulnerability while generality decreased is due to increased importance of specialist predators that interact with a subset of prey, while increasing the overall number of predators that each prey is consumed by. Indeed, the significant decline in proportion of omnivory at most stations in the Bothnian Sea ( Figure S5j and Table S3) supports the specialization of predators, and increased proportion of predators has been reported in the Bothnian Sea before (Törnroos et al., 2019). However, only one guild (benthic invertebrates) was investigated in this study, and predation from other guilds should be considered in future research.
In the Bothnian Sea and Baltic Proper, increases in taxa richness indicate local species invasions during the study period, but these invasions, while increasing the dispersion, did not significantly change the means of most trophic network metrics, which suggests that the invasions had little effect on the benthic invertebrate trophic network architecture here investigated. However, rates of species introductions in marine ecosystems are increasing (Hulme, 2009), and monitoring for consequences on trophic networks should continue (Ojaveer et al., 2017).
Eutrophication is known to simplify food webs through fewer nodes, links, shorter chains, and increased connectance O'Gorman et al., 2012) (Karlson et al., 2002;Pearson & Rosenberg, 1978), and earlier reference periods would be needed to evaluate food web metrics from communities unaffected by eutrophication.
However, the benthic monitoring program did not begin to regularly sample until the 1980s, and thus, evaluating changes before eutrophication began to impact the Baltic Sea is difficult.
Climate change and eutrophication interact to produce largescale impacts on trophic networks (Gilarranz et al., 2016). Climate change may nullify nutrient reductions and increase hypoxic zones in the Baltic due to temperature rise (Reusch et al., 2018 The Baltic Proper and Bothnian Sea are more vulnerable to hypoxia than other basins due to longer water residence times and more pronounced stratification (Conley et al., 2009), and the low species richness indicates that any local extinctions will bring severe impacts for trophic network functioning. While we did not observe the potential network effects of temperature and eutrophication in these basins, the interaction of local, regional, and global pressures should be further investigated, as these have been shown to cause both food web simplification and functional reorganization in the Baltic (Pecuchet et al., 2020;Törnroos et al., 2019). Climate change has been connected to increased connectance, which could be behind the changes seen in the Skagerrak, and increased variability of generality, which we could see evidence of in Figure 2e  , but more studies focusing on these variables with more complete abiotic data will be required to make definite conclusions.
Of course, there are limitations with this study. Long-term monitoring datasets, while valuable, come with difficulties of standardization and shifting methods (Magurran et al., 2010). The monitoring program samples once per year in May and could have missed species present at different times of the year. Thus, we were not able to consider the seasonal differences in food webs. In addition, only considering two decades of monitoring data could also skew our findings, as we may miss shifts in community composition. We also focus only on benthic invertebrate macrofauna, and it would be beneficial to expand the study to include lower trophic levels, such as meiofauna, as well as higher trophic levels. We assume feeding links from literature, and this information is incomplete due to the stochastic nature of species interactions, and will likely overestimate the number of links (Roslin & Majaneva, 2016). In addition, we were not able to quantify the strength of the feeding links (Kortsch et al., 2021) or present high resolution of basal resources. We hope that future studies will consider these limitations and aid in better trophic network resolution. Nevertheless, we were able to show changed trophic networks utilizing existing longterm monitoring data, which has important management implications.

| Management suggestions and conclusions
We found a simplification of modeled Baltic Sea macrozoobenthic trophic networks during the period 1980 to 2019, and the Skagerrak region was most affected, in agreement with our predictions. We also detected changes in benthic invertebrate food web structures that were not related to taxonomic richness, underlying the importance of not only monitoring benthic communities but also their interactions through metrics to understand and predict changes in ecosystem functions. Commercially important fish are likely to depend on benthic invertebrates in aphotic habitats, and decreasing benthic invertebrate biomass is a concern that there may be trophic cascades resulting in decreased fish stocks in the future (Snickars et al., 2015).
Conservation of benthic invertebrate trophic networks is thus essential to maintain valuable ecosystem functions. Other studies of anthropogenic pressures on whole ecosystem marine food webs have found that marine-protected areas improve food web resilience (Gilarranz et al., 2016). We recommend further studies investigating whether expanding and adding Swedish marine protected areas, particularly in the Skagerrak, would reduce anthropogenic pressures and conserve benthic invertebrate trophic networks.

ACK N OWLED G EM ENTS
We are grateful to all those involved in the Swedish national benthic monitoring program for collecting samples throughout the

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data is available at the Dryad Digital Repository: https://doi.