Combined Effects of Environmental Drivers on Marine Trophic Groups – A Systematic Model Comparison

Ehrnsten, Eva; Bauer, Barbara; Gustafsson, Bo G.

doi:10.3389/fmars.2019.00492

ORIGINAL RESEARCH article

Front. Mar. Sci., 02 August 2019

Sec. Marine Ecosystem Ecology

Volume 6 - 2019 | https://doi.org/10.3389/fmars.2019.00492

Combined Effects of Environmental Drivers on Marine Trophic Groups – A Systematic Model Comparison

$\r\nEva Ehrnsten,*&#x;$ Eva Ehrnsten^1,2*†

Barbara Bauer^2,3,4†

Bo G. Gustafsson^1,2

¹Tvärminne Zoological Station, University of Helsinki, Hanko, Finland
²Baltic Sea Centre, Stockholm University, Stockholm, Sweden
³Institute of Ecology, Friedrich Schiller University Jena, Jena, Germany
⁴German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany

The responses of food webs to simultaneous changes in several environmental drivers are still poorly understood. As a contribution to filling this knowledge gap, we investigated the major pathways through which two interlinked environmental drivers, eutrophication and climate, affect the biomass and community composition of fish and benthic macrofauna. For this aim, we conducted a systematic sensitivity analysis using two models simulating the dynamics of benthic and pelagic food webs in the Baltic Sea. We varied environmental forcing representing primary productivity, oxygen conditions and water temperature in all possible combinations, over a range representative of expected changes during the 21st century. Both models indicated that increased primary productivity leads to biomass increase in all parts of the system, however, counteracted by expanding hypoxia. Effects of temperature were complex, but generally small compared to the other drivers. Similarities across models give confidence in the main results, but we also found differences due to different representations of the food web in the two models. While both models predicted a shift in benthic community composition toward an increased abundance of Limecola (Macoma) balthica with increasing productivity, the effects on deposit-feeding and predatory benthic groups depended on the presence of fish predators in the model. The model results indicate that nutrient loads are a stronger driver of change for ecosystem functions in the Baltic Sea than climate change, but it is important to consider the combined effects of these drivers for proper management of the marine environment.

Introduction

Marine ecosystems are threatened by multiple anthropogenic pressures (Halpern et al., 2008). Eutrophication and associated hypoxia have deteriorated coastal seas and continue to pose a major threat to marine ecosystems globally (Cloern, 2001; Gray et al., 2002; Breitburg et al., 2018). Climate change affects the life-histories and distributions of marine species, and is recognized as a major driver of food-web reorganizations (Poloczanska et al., 2013; Cloern et al., 2016). Unraveling the effects of these and other interacting drivers on the marine ecosystem is imperative for proper ecosystem management to ensure the future functioning of marine ecosystems and provision of ecosystem services (Atkins et al., 2011; Cloern et al., 2016).

Most studies so far have looked at the effects of isolated drivers (Richardson et al., 2012; Pascual et al., 2017), although Crain et al. (2008) show that synergistic and antagonistic interactions among drivers are common. It is logistically difficult to test the effects of multiple drivers across a large gradient of levels in empirical studies (Boyd et al., 2018). Process-based ecological models are powerful tools to conduct such tests (Christensen and Walters, 2004; Coll et al., 2015; Seidl, 2017). Additionally, such models can be used to study the relative importance of direct (e.g., physiological) and indirect (e.g., food-web mediated) effects on ecosystem components (Blackford, 2002; Fulton et al., 2004; Condie et al., 2012).

Here, we use two different process-based ecological models to systematically analyze the responses of benthic and pelagic food webs to multiple combinations of drivers across a large range of levels. In general, marine ecological modeling is associated with large uncertainties, especially when implemented on large temporal and spatial scales. Comparing and synthesizing the results of multiple models is one way to address that uncertainty (Fulton, 2010; Uusitalo et al., 2016). The models complement each other in terms of spatial and temporal scales and level of abstraction. The Ecopath with Ecosim (EwE) model focuses on the ecosystem scale and describes food web interactions in both benthic and pelagic compartments of the central Baltic Sea (Bauer et al., 2018). The Benthic Macrofauna (BM) model focuses on local macrobenthos in the area (at 40 m depth) and uses a physiologically more detailed, mechanistic description of the system (Ehrnsten et al., 2019).

There are a number of previous modeling studies investigating the potential future responses of marine ecosystems to various drivers with the help of future scenario simulations, including those using the EwE approach (Evans et al., 2013; Niiranen et al., 2013; Österblom et al., 2013) and a previous version of the BM model (Timmermann et al., 2012). Scenarios are valuable to gain an insight into possible ranges of future states of the ecosystem (Cheung et al., 2009; Brochier et al., 2013). However, as they typically vary several drivers (e.g., temperature, oxygen, and salinity) at the same time, they only provide a limited mechanistic understanding of the pathways by which such drivers and their interactions affect ecosystems. In their recent review of studies on climate change impacts on ecosystem services, Runting et al. (2017) pointed out that a systematic analysis of both cumulative and relative effects of multiple drivers is needed, but often overlooked. Such an analysis enables a better interpretation of the results of scenario studies, by identifying which specific drivers, or their interactions, drive scenario results. In addition, an estimation of the uncertainty associated with studies solely looking at the effects of single drivers is made possible.

Both models are parameterized to be representative of the Baltic Sea, a coastal sea that due to its location and hydrography is particularly exposed to multiple anthropogenic pressures, such as nutrient loading and climate change (Karlson et al., 2002; Belkin, 2009; Carstensen et al., 2014). Since the 1950s, increased nutrient loading has led to increased primary production (Gustafsson et al., 2012), providing food for zooplankton and planktivorous fish in the pelagic food web (Eero et al., 2016). It has also increased sedimentation of organic matter which acts as a food source for benthic communities (Cederwall and Elmgren, 1980). The degradation of organic matter consumes oxygen, and has led to spreading hypoxic and anoxic areas in Baltic Sea (Conley et al., 2009; Carstensen et al., 2014). Hypoxia diminishes the diversity and functionality of benthic communities (Villnäs et al., 2013; Gammal et al., 2017). The effects of decreasing oxygen concentrations on growth and survival differ among marine groups, with demersal fish and crustaceans generally showing high sensitivity while bivalves can survive extended periods of complete anoxia by closing their valves (Gray et al., 2002). In the Baltic Sea, the common deposit-feeding amphipods Monoporeia affinis and Pontoporeia femorata are highly sensitive to hypoxia while the predatory isopod Saduria entomon is among the most tolerant (Timmermann et al., 2012). Hypoxia also has negative impacts on reproductive conditions for demersally spawning fish, such as cod (MacKenzie et al., 2000; Margonski et al., 2010).

At the same time, the Baltic Sea has been one of the fastest warming marine ecosystems in the world during the last decades (Belkin, 2009), and water temperatures are projected to increase further during the 21st century (Meier et al., 2011, 2012a). Temperature directly affects most physiological rates across the whole food web, with varying effects of warming on different species and their interactions (Pörtner and Farrell, 2008; Brander, 2010). It increases growth rates of herring larvae (Hakala et al., 2003) and has a positive effect on herring and sprat productivity (Margonski et al., 2010; Bartolino et al., 2014). On the other hand, there are indications that several common marcobenthic species in the Baltic Sea have declined due to rising temperatures, including the glacial relicts M. affinis and S. entomon (Segerstråle, 1962; Rousi et al., 2013). Beukema et al. (2009) found negative effects of warming on the growth, reproduction and survival of the Baltic clam Limecola (Macoma) balthica in the Wadden Sea, suggesting that warming might affect also Baltic populations in the near future. Increasing temperature also increases the sensitivity of benthic groups to hypoxia (Vaquer-Sunyer and Duarte, 2011). Furthermore, rising temperatures increase the microbial degradation of organic matter, consuming oxygen, and reducing food availability for macrobenthos.

The aim of this study is to identify the most important environmental drivers of selected benthic and pelagic ecosystem components, potential non-linear responses in the ecosystem, and to explore possible interactive effects between drivers. We focus on the effects of two major external pressures of the marine ecosystem, nutrient loading and climate, on the biomass and community composition of macrobenthos and pelagic and demersal fish (Figure 1). More specifically, we hypothesize that: (1) increasing productivity increases biomass of planktivourous fish and detritivorous benthos, (2) hypoxia counteracts the effects of productivity and alters community composition with strongest negative effects on sensitive groups like deposit-feeders and demersally spawning and feeding fish, and (3) increasing temperature has interactive effects with productivity and oxygen conditions, improving growth conditions for small pelagic fish (herring and sprat) while diminishing those of BM.

FIGURE 1

Figure 1. Overview of interacting environmental drivers and marine ecosystem components considered in this study.

Materials and Methods

We use two ecological models to study the effects of changes in the environmental drivers: (1) productivity, (2) oxygen conditions, and (3) water temperature on the biomass of BM, demersal and pelagic fish in the Baltic Sea (Figure 1). These environmental drivers are represented in the two models as forcing affecting different variables (Table 1). Below the models are described in general terms with a focus on the implementation of forcing affecting ecosystem components in each model. References to full model descriptions including mathematical formulations are given in the respective sections.

TABLE 1

Table 1. Studied environmental drivers, corresponding forcing, and functional groups directly affected by the forcing in the two models.

Study System

The Baltic Sea is a semi-enclosed brackish-water basin with a short geological history, and, consequently, a low diversity of macroscopic species and simple food-webs characterized by a few dominating species (Snoeijs-Leijonmalm, 2017). The sea encompasses positive north to south gradients in environmental conditions, such as salinity and temperature, and extent of anthropogenic pressures, such as nutrient loads, as well as in species and functional diversity. Here, we focus on the Baltic Proper basin including offshore pelagic food-webs and benthic food-webs in the deep, aphotic (non-sunlit) zone. Soft sediments inhabited by the facultative suspension/deposit-feeding Baltic clam Limecola balthica dominate the sea (Gogina and Zettler, 2010; Gogina et al., 2016). Other common deposit-feeders include the amphipods M. affinis and P. femorata and the non-native polychaetes Marenzelleria spp. Key predators in the benthic system include invertebrates, such as the isopod S. entomon, and demersal fish, such as European flounder (Platichthys flesus) and cod (Gadus morhua), the latter linking the benthic and pelagic systems. Key fish species in the pelagic system are the planktivorous sprat (Sprattus sprattus) and herring (Clupea harengus) and their main predator, cod.

Ecopath With Ecosim (EwE) Model

The EwE modeling approach and software¹ is used to describe biomass flows among ecosystem components (Ecopath) and the system’s dynamic response to external forcing such as changing fishing effort or environmental conditions (Ecosim). The biomass dynamics of a functional group is dependent on the biomass production of the group (assimilated consumption minus respiratory losses), predation on the group, fishing and other sources of mortality, as well as egg production (juveniles) or the number of juveniles maturing into adults (adults) for the groups represented by multiple life stages (“stanzas”). The theory behind EwE, its implementation as a software and its limitations have been extensively described by Pauly et al. (2000), Christensen and Walters (2004), and Plagányi (2007). EwE has been used numerous times to investigate environmental and physical effects on ecosystem state, especially on fish production (e.g., Coll et al., 2009; Morato et al., 2009; Libralato et al., 2015).

The EwE model of the offshore central Baltic Sea (ICES subdivisions 25–29, excl. the Gulf of Riga) used here represents the ecosystem as 21 living functional groups and species and one detritus group. Living groups encompass all trophic levels from phytoplankton to gray seals (Halichoerus grypus), including the four commercially most important fish species (cod, flounder, herring, and sprat). The latter are each represented as two stanzas (adult and juveniles). Macrobenthos groups in EwE are S. entomon, two bivalve groups (L. balthica and Mytilus spp.), and “other macrobenthos,” consisting of oligochaetes, polychaetes (e.g., Bylgides sarsi), and amphipods. The model is parameterized to represent the ecosystem state in the early 2000s. It has been used previously to explore alternative fisheries management strategies (Bauer et al., 2019) and spatial effects of nutrient load management in the Baltic Sea (Bauer et al., 2018). The model parameterization used here is the same as in Bauer et al. (2018), except for small modifications described in Bauer et al. (2019). Fishing mortalities were kept constant at their values in Ecopath (representative of 2004) for all analyses (adult/juvenile cod: 0.84/0.05, flounder: 0.59/0.02, herring: 0.13/0.05, sprat: 0.4/0.14).

The effects of environmental drivers in EwE are implemented in three ways: primary production forcing, forcing affecting consumption and egg production forcing (Christensen et al., 2008). All three types of forcing are implemented as a time series of multipliers applied on certain model parameters (Table 1). The first value of each time series has to be 1, that is, all parameters are initially at their baseline values. Primary production forcing values are applied as a multiplier of the r parameter (maximum production to biomass ratio, P/B) in the equation relating primary producer biomass to production via a saturating relationship, of which r is the asymptote (Christensen et al., 2008). Forcing affecting consumption rate is here applied as a multiplier on the “effective search rate” of predators, which is equivalent to the attack rate parameter commonly used in functional response formulations (Holling, 1959). Egg production forcing is applied as a multiplier on the number of eggs produced by default (Walters et al., 2010). In previous studies on another central Baltic Sea EwE model, forcing time series were derived from time series of environmental drivers (Tomczak et al., 2012; Niiranen et al., 2013) and we assume similar relationships between forcing and environmental drivers here. Consumption rates of cod and several benthic groups are related to the extent of oxic areas (OxAr) (Sarvala, 1971; Eero et al., 2015; Köster et al., 2017). Cod egg production is related to cod reproductive volume (CodRV, the volume of water with salinity > 11 psu and oxygen concentration > 2 mL L^-1), and herring and sprat egg production to sea surface temperature in August (TSurf) (MacKenzie et al., 2000; Margonski et al., 2010). Forcing functions affecting several groups (OxAr, TSurf) have the same effect on all groups’ consumption or egg production rates. Representing environmental forcing as linear multipliers on certain rates is a major difference to the more mechanistic BM. In the latter, forcing is applied as time series of absolute values of a given environmental variable (e.g., temperature in °C, see below) affecting ecological rates described by different mechanistic relationships.

Benthic Macrofauna (BM) Model

The BM model simulates the flows of carbon through functional groups of macrofauna and their food banks in the sediment. The model is based on the physiological macrofauna model of Timmermann et al. (2012), but uses some novel formulations and parameterizations, including a dynamic feedback between benthic fauna and sediment organic carbon and a new functional response formulation for predation on multiple prey (Ehrnsten et al., 2019). The carbon biomass dynamics of a functional group is described by the physiological processes of ingestion, assimilation, egestion and respiration as well as the population-level processes of predation, other mortality and recruitment, similarly to the benthic biology module in the European Regional Seas Ecosystem Model (Ebenhöh et al., 1995; Butenschön et al., 2016; Lessin et al., 2019).

Here, we use the implementation of the model representing a 40 m deep soft-sediment area in the western Baltic Proper described in Ehrnsten et al. (2019). It includes three functional groups of macrofauna and two sediment food banks (freshly deposited and older organic carbon). The freshly deposited carbon is consumed by the surface deposit-feeder functional group (hereafter referred to as “deposit-feeders”) dominated by the amphipods M. affinis and P. femorata. The second functional group only comprises the key species L. balthica, which is capable of suspension-feeding when plankton concentrations are high. This behavior is modeled as a switch to feeding on suspended matter when bottom water chlorophyll-a concentrations are > 2 mg m^-3 (see Timmermann et al., 2012). However, in the present implementation the chlorophyll-a threshold is rarely reached, and L. balthica mostly consumes organic carbon from the older sediment food bank. Invertebrate predator/scavengers (hereafter “predators”) are a taxonomically diverse group including the priapulid Halicryptus spinulosus, the isopod S. entomon and the polychaetes B. sarsi and Hediste diversicolor. Predators feed on both other macrofauna groups with a strong preference for deposit-feeders.

Forcing to the BM model consists of daily time-series of bottom water temperature (TBot), concentration of oxygen (O2Conc), concentrations of suspended plankton and detritus, as well as organic matter sedimentation rate (Sed). These are provided by the physical-biogeochemical BALTSEM model (Gustafsson et al., 2012; Savchuk et al., 2012). Sedimentation is accumulated in the food banks, which are consumed by the fauna and degraded at temperature-dependent rates. Temperature also affects the production of macrofauna through ingestion, basal respiration, and mortality rates. These direct temperature effects are modeled according to a standard, exponentially increasing Q₁₀-formulation. Low oxygen concentrations cause mortality of the fauna via a temperature-dependent negative exponential function, with different sensitivities toward hypoxia for the three groups (deposit-feeders > L. balthica > predators) (Timmermann et al., 2012).

Mean values for January 1st 2000–2005 simulated by Ehrnsten et al. (2019) were used as initial conditions for the state variables. Otherwise all model settings and parameter values were applied as in Ehrnsten et al. (2019).

Analysis of the Effects of Environmental Drivers on Ecosystem Components

To test our hypotheses on the effects of changes in productivity, oxygen conditions, and temperature compared to recent conditions on a number of ecosystem components, both models were run under systematically varied forcing (Table 1). Subsequently, we calculated biomass change of key functional groups compared to the baseline run. In EwE the baseline meant a simulation with all forcing set to the value of 1, i.e., all affected parameters were kept at their default values. In BM, baseline forcing were the daily values of all forcing time series simulated by BALTSEM averaged over the years 2000-2005 (Figure 2). In both models, productivity- and temperature-related forcing were varied in the range of 0.5–2 times the baseline. This range corresponds roughly to changes that can be realistically expected to occur in the region in the 21st century (Meier et al., 2011, 2012a; Saraiva et al., 2018), except for a decrease in temperature which is included in the test for completeness. Cod reproductive volume (CodRV) was tested in the range of 0–10 times the baseline, approximately corresponding to a range of no reproductive volume to 1500 km³ (the latter representing a tenfold increase in CodRV, which could only happen under an extremely optimistic scenario without climate change). OxAr were tested in the range of 0.75–1.5, approximately corresponding to half to the whole of the Baltic Proper seafloor being oxic. In BM, bottom water oxygen concentrations (O2Conc) of 0–1 times the baseline was tested, as the baseline corresponds to approximately fully oxygen saturated conditions, so only decreases are possible. We show biomass changes under selected forcing combinations, all combinations are shown under Supplementary Figure S2 (EwE) and Supplementary Figure S3 (BM). The sensitivity of EwE model results to important parameter settings are shown in Supplementary Figure S1.

FIGURE 2

Figure 2. Forcing time-series used in the BM model, showing the tested levels of sedimentation (A), oxygen concentration (B), and bottom temperature (C) over 1 year. Level 1 denotes the baseline.

In EwE simulations were run for 400 yearly time steps under constant forcing. That is, we applied a time series of 400 steps, the first one being 1, the next 399 being the tested value of the forcing. In BM forcing was repeated as an identical seasonal cycle every year for 100 years, using daily time steps. The range of tested values of the forcing time-series in BM are shown in Figure 2. The tested time periods in both models were chosen to be long enough for the dynamics to surpass the transition period and reach an equilibrium state, defined as the lack of interannual increasing or decreasing trends in biomass of the studied functional groups.

Results

Ecopath With Ecosim (EwE) Model

The P/B forcing, representing phytoplankton productivity, was the most influential factor affecting community biomass (Figure 3). Cod and sprat showed non-linear positive responses, adult flounder, L. balthica and “other macrobenthos” showed linear positive responses and S. entomon a saturating positive response. Herring biomass did not respond strongly to P/B or to any other forcing (Figure 3C and Supplementary Figure S1.2).

FIGURE 3

Figure 3. Change in biomass of cod (A), flounder (B), herring (C), sprat (D), Limecola balthica (E), “other macrobenthos” (F), and Saduria entomon (G) as a function of relative changes in P/B and OxAr, where 1 denotes baseline forcing. Biomass change is relative to baseline equilibrium value. Error bars show the variability in biomass change due to varying TSurf and CodRV. Horizontal line indicates no change. Results from simulations with CodRV = 0, resulting in 0 cod biomass, are excluded from the figure. Forcing abbreviations are explained in Table 1.

P/B and extent of OxAr had indirect effects on several groups due to trophic relationships. OxAr directly increased consumption rates of cod and some other benthic groups on their prey (Table 1). Thus, OxAr had a positive effect on cod (Figure 3A) and L. balthica (Figure 3E), due to their increased consumption of benthic food. In contrast, OxAr had a weak negative effect on sprat, probably as an indirect effect of increased cod biomass under higher OxAr, as sprat is a main non-benthic prey of cod (Figure 3D and Supplementary Figure S1.1). Even though OxAr also increased consumption rates of the “other macrobenthos” group, OxAr ≥ 1 resulted in lower biomasses than OxAr = 0.75 (Figure 3F), presumably due to higher cod biomass causing higher predation mortality under higher OxAr values. For the same reason, the highest tested level (1.5) of OxAr, resulted in lower S. entomon biomasses than at OxAr < 1.5, although only in combination with high (> 1) P/B (Figure 3G).

Although surface temperature (TSurf) and cod reproductive volume (CodRV) only directly affected fish groups, they also impacted some benthic groups (Figure 4) due to a chain of indirect effects. Temperature had a positive effect on juvenile sprat (Figure 5D) which feeds on zooplankton, thereby freeing phytoplankton from predation pressure, and increased phytoplankton increased L. balthica biomass (Figures 3E, 4A). Increasing CodRV had a negative effect on S. entomon (Figure 4B) due to increased predation especially by juvenile cod (Figure 5C).

FIGURE 4

Figure 4. Change in biomass of L. balthica as a function of relative changes in TSurf (A) and S. entomon as a function of relative changes in CodRV (B) at four different levels of P/B, where 1 denotes baseline forcing. Biomass change is relative to baseline equilibrium value. Error bars show the variability in biomass change due to varying OxAr and CodRV (A) and OxAr and TSurf (B). Horizontal line indicates no change. Results from simulations with CodRV = 0, resulting in 0 cod biomass, are excluded from the figure. Forcing abbreviations are explained in Table 1.

FIGURE 5

Figure 5. Change in biomass of adult (A) and juvenile (C) cod and adult (B), and juvenile (D) sprat as a function of their respective egg production forcing (Table 1) at four different levels of P/B, where 1 denotes baseline forcing. Biomass change is relative to baseline equilibrium value. Error bars show the variability in biomass change due to varying forcing factors other than P/B and egg production. Horizontal line indicates no change. Forcing abbreviations are explained in Table 1.

We found that the effects of egg production forcing were modified based on productivity (Figure 5). Except for P/B = 0.5, where adult cod biomass was mostly insensitive to changes in CodRV, adult cod biomass increased up to a point with increasing CodRV but saturated (at P/B = 1) and even decreased (at P/B > 1) when CodRV increased above ∼1.75, hinting at a negative effect of a large number of juveniles produced on the adult population, especially at high primary productivity. The relationships between TSurf, P/B and sprat were very similar: at P/B = 0.5 increasing TSurf had no effect, at P/B ∼1 it had a positive effect on both adult and juvenile sprat and at P/B > 1 it had a positive effect on juvenile, but not on adult sprat. Thus, competition between juveniles and adults of both species was higher at high P/B. That seems counterintuitive, as at high P/B there is more primary production to support higher trophic levels. However, the production of juveniles is, without any change in egg forcing functions, already relatively high at high P/B because of the high adult fish biomass implying a large number of mature fish. A further increase when applying an increase in egg production forcing thus means a disproportionately large (compared to the increased P/B) increase in juvenile fish and their predation pressure on groups that also serve as prey for the adults. High juvenile sprat biomass had a stronger negative effect on adult sprat biomass than it was the case for cod, due to a very high diet overlap between juveniles and adults of sprat.

Benthic Macrofauna (BM) Model

Organic matter sedimentation (Sed) as an indicator of productivity and food availability was the most influential forcing affecting benthic community biomass. It had a mostly linear positive effect on L. balthica and predator biomasses (Figures 6, 7). Temperature (TBot) counteracted the positive effects of sedimentation on these two groups: increased bottom water temperature decreased the slope of the biomass response to sedimentation (Figure 6). For L. balthica, this was due to both to increased degradation rate of the food source counteracting the potential increase in ingestion rate, and to increased metabolic losses with increasing temperature. Oxygen concentrations (O2Conc) showed synergistic effects with sedimentation, i.e., increasing oxygen concentrations increased the positive response to sedimentation for the two groups (Figure 7).

FIGURE 6

Figure 6. Change in biomass of deposit-feeders (A), L. balthica (B), and predators (C) as a function of relative changes in Sed and TBot, where 1 denotes baseline forcing. Biomass change is relative to baseline equilibrium value. Error bars show the variability in biomass change due to varying O2Conc. Horizontal line indicates no change. Forcing abbreviations are explained in Table 1.

FIGURE 7

Figure 7. Change in biomass of deposit-feeders (A), L. balthica (B), and predators (C) as a function of relative changes in Sed and O2Conc, where 1 denotes baseline forcing. Biomass change is relative to baseline equilibrium value. Error bars show the variability in biomass change due to varying TBot. Horizontal line indicates no change. Forcing abbreviations are explained in Table 1.

The effects of O2Conc were sigmoidal and saturating for all groups: oxygen levels below a certain threshold led to extinction, with a strong increase in biomass above the threshold (Figure 8). The form of the response was group-specific with L. balthica being least sensitive to hypoxia. Reduced O2Conc affected predators both directly through mortality and indirectly through reduction in prey biomass, leading to a higher sensitivity than would be predicted from hypoxia mortality rates alone.

FIGURE 8

Figure 8. Change in biomass of deposit-feeders (A), L. balthica (B), and predators (C) as a function of relative changes in O2Conc and TBot, where 1 denotes baseline forcing. Biomass change is relative to baseline equilibrium value. Error bars show the variability in biomass change due to varying Sed. Horizontal line indicates no change. Forcing abbreviations are explained in Table 1.

The effects of forcing on deposit-feeders were complex with interacting effects of all three forcing types. Low sedimentation rate, indicating low food availability, in combination with high temperature and/or low oxygen concentration led to severely reduced biomass or even extinction. This means that the ability to survive in increased temperatures and/or reduced oxygen concentrations was greatly reduced by food shortage. However, at forcing combinations outside of this “lethal space,” i.e., when assuming a higher sedimentation rate and/or more suitable abiotic conditions, the biomass response to increasing sedimentation was slightly negative (Figures 6, 7). The effects of temperature were slightly positive except at hypoxia or very low sedimentation rate (Figures 6, 8).

The effects of sedimentation and temperature on deposit-feeders was opposite to those on L. balthica and predators. This pattern was due to food web interactions that became an important driver of deposit-feeder biomass responses under non-lethal environmental conditions, as can be seen from the relationship between benthic predators and their prey (Figure 9). At severe hypoxia there was no correlation between predators and prey. In normoxic conditions, there was a positive relationship between L. balthica and predators, but a negative relationship between deposit-feeders and predators. This indicates that L. balthica biomass drives predator biomass, but there is little reciprocal effect as L. balthica biomasses are generally one order of magnitude larger than predator biomasses (bottom-up control). For deposit-feeders, the negative relationship indicates top-down control, i.e., that predation is a significant driver of deposit-feeder biomass.

FIGURE 9

Figure 9. Relationship between biomass of L. balthica and predators (A) and surface deposit-feeders and predators (B) at hypoxia (O2Conc < 0.1, blue) and normoxia (baseline O2Conc, yellow). Biomass values are in mg C m^-2. Forcing abbreviations are explained in Table 1.

Discussion

We have used two distinct food web models to study the effects of environmental drivers on marine ecosystem components in the Baltic Sea. Both models showed the same main responses to environmental drivers. Primary productivity as a basis of the food web was the strongest driver of community biomass, having a positive effect on most groups as hypothesized. In EwE, all groups but herring showed a clear positive response to increases in phytoplankton production. In BM, two out of three groups (L. balthica and predators) showed a positive response to increased sedimentation. The other two groups (herring in EwE and deposit-feeders in BM) showed weak responses to productivity. The effects of oxygen conditions, measured as extent of oxic areas (EwE), cod reproductive volume (EwE) or oxygen concentration (BM) forcing, were mostly non-linear. A strong negative effect of low oxygen concentrations on all benthic groups was evident in BM, but biomasses saturated quickly as concentrations increased. Similarly, cod reproductive volume had a strong effect on cod biomass at low values, but the biomass response saturated quickly. Temperature effects were more variable than initially hypothesized (e.g., contrasting effects between adult and juvenile sprat and between L. balthica and deposit-feeders). Nevertheless, they were generally small compared to the other two environmental drivers within the tested range, representative of anticipated change in the 21st century.

The models showed several interactive effects of the studied drivers, both synergistic and antagonistic. In many cases these effects were indirect and modified by food web interactions. For example, the result that the relationship between sedimentation and deposit-feeders was slightly negative in BM seems counterintuitive, but can be explained by indirect effects through the food web. Sedimentation had a positive effect on the food availability for L. balthica, which represents the majority of benthic biomass. Increased biomass of L. balthica was one of the main drivers of predator biomass, which in turn asserted an increased predation pressure on deposit-feeders. This is a phenomenon known as “apparent competition” when two prey biomasses are negatively related due to a shared predator (Holt, 1977; Holt and Lawton, 1994). Thus, the indirect effect of increased sedimentation via intensifying apparent competition between L. balthica and deposit-feeders had a stronger effect on deposit-feeder biomass than the direct effect via increased food availability. Tight predator-prey coupling between S. entomon and the deposit-feeder M. affinis has been recorded in the northern Baltic Sea where L. balthica was absent (Sparrevik and Leonardsson, 1999; Englund et al., 2008; Englund and Leonardsson, 2008). In an area with few amphipods present, S. entomon predated actively L. balthica (Ejdung and Bonsdorff, 1992), but to our knowledge there are no studies that included both prey types that could have revealed apparent competition in the lab or field.

Our study shows that the choice of system boundaries are important to consider when estimating impacts of environmental drivers on ecosystems. In EwE, the benthic predator S. entomon did not show a similar strong response to increased environmental productivity as in BM. This is at least partly due to a stabilizing effect of predation by cod. Cod strongly benefited from increased primary productivity through increased availability of one of its main prey, sprat. Thus, under high productivity cod caused high predation mortality on S. entomon, counteracting the positive effect of increased prey availability for S. entomon. In BM, predation on benthic predators is included as a closure term, and is therefore less dynamic than in EwE that explicitly models the dynamics of the top predator (cod) feeding on both predatory benthos and pelagic prey. Thus, EwE showed an apparent competition effect between sprat and S. entomon, mediated by cod predation, that was missing from BM. Real world marine ecosystems are extremely complex, and models are by definition a simplification of that immense complexity. There is no “best model” for any given system, as the choice of elements to be included in any model depends on the question to be addressed, Data Availability and planned uses of the model (Fulton et al., 2003; Fulton, 2010). There are a number of modeling approaches describing marine ecosystems and populations that differ in their representations of main processes, but they can all be used the generate projections of how the system reacts to external pressures. As such projections are sensitive to model formulation, as shown by this as well as other studies (e.g., Fulton et al., 2003; Gårdmark et al., 2013), they should be viewed as alternative hypotheses and as a potential basis for empirical studies.

Besides model scale, another difference between the models is the implementation of forcing. In EwE, it is assumed that the effects of environmental drivers on the affected parameters, such as consumption or egg production rates, are linear (Christensen et al., 2008). Additionally, due to the complexity and large number of groups included in our EwE model, each environmental forcing only affected a selected number of groups to keep model dynamics more tractable. In BM, environmental variables were used as forcing in a more realistic manner. For example, temperature effects were only included in EwE on herring and sprat egg production, while in BM all state variables were affected by temperature through more mechanistic Q₁₀-relationships for several process rates. We note that there are recent approaches to incorporate more mechanistic forcing in EwE models as well (de Mutsert et al., 2012; Libralato et al., 2015). Despite these very different approaches, both models showed relatively small effects of temperature on the studied food-webs. Overall temperature responses were largely dampened by compensatory effects in different parts of the system, e.g., in BM the increased potential ingestion rates of L. balthica and deposit-feeders were counteracted by faster degradation of the food sources under increasing temperatures. It has been argued that estuarine areas with strong seasonality and low diversity like the Baltic Sea should harbor organisms with large tolerance ranges to environmental drivers such as temperature (Jernelöv and Rosenberg, 1976). Nevertheless, some of the dominant species of the deep zone included in our models (M. affinis, Mysis relicta, and S. entomon) are glacial relicts that are believed to be cold-adapted (Segerstråle, 1962). The reason for a lack of a larger temperature response from these groups may partly be that some important temperature effects are not captured by either model, e.g., on species distributions, phenology or body sizes (Sommer et al., 2012; Poloczanska et al., 2013; Birchenough et al., 2015).

There are additional relevant environmental drivers with potential effects on both the benthic and pelagic food web in the Baltic Sea that were not considered here. Salinity in this study was only represented as a component of cod reproductive volume, as in an earlier EwE model of the Baltic other salinity effects had little impact on model dynamics (Niiranen et al., 2013). However, salinity is considered a major determinant of present and future species distributions in the Baltic Sea (Gogina and Zettler, 2010; Ojaveer et al., 2010; Jokinen et al., 2015; Vuorinen et al., 2015). Oxygen conditions on the seafloor do not only affect the benthic communities through permanent low oxygen concentrations, but also through occasional extreme hypoxic events. These occur with increasing frequency under elevated eutrophication (Conley et al., 2011). From our two models, only the BM is suitable to study the effects of short extreme events, due to its daily timescale. Here, we focused on effects of permanent changes in environmental drivers, such as elevated water temperatures over the whole year, which could be studied by both models.

In general, both models consider a fixed number of functional groups in a non-spatial framework representing the central Baltic ecosystem. Thus, potential effects of migration, both small-scale movement of species to avoid adverse conditions (e.g., local hypoxia), and invasion of completely new species, are beyond the scope of this study. Small-scale movement of species could be addressed using spatial models like in Bauer et al. (2018). The impacts of specific novel species or potential invaders could be addressed in a straightforward way in both modeling frameworks by including them as additional functional groups. In this case, the biomass responses of current functional groups presented in this study would serve as a baseline. Changes in those responses compared to the baseline due to the inclusion of invasive groups are a quantitative measure of how much a novel ecological configuration affects the system’s sensitivity to environmental pressures (Radeloff et al., 2015). It should also be noted that the results may not be directly generalizable to shallow coastal systems or other basins of the Baltic Sea, as environmental conditions and food-web structures in those habitats are quite different from the ones considered here. For example, suspension-feeding had a minor role in the current study, but may be important for pelagic-benthic organic matter fluxes and benthic dynamics in other areas (Kautsky and Evans, 1987; Lessin et al., 2019).

It is uncertain how primary production in the Baltic Sea is going to develop during the next century. Climate change alone would increase primary production due to increased water temperatures and nutrient runoff from land (Meier et al., 2011, 2012b; Saraiva et al., 2018). Our simulations indicate that this would increase biomass and alter community composition in both benthic and pelagic parts of the system. Regarding the fish community, EwE predicts that further increases in primary productivity favor sprat and cod equally. However, when coupled with hypoxia, positive effects may only be seen in sprat. Negative effects of hypoxia on cod have been previously proposed based on empirical data (Eero et al., 2015; Casini et al., 2016), as well as positive effects of increased productivity on sprat production (Chassot et al., 2007; Eero et al., 2016). The model indicated that positive effects of increased productivity could be limited by density-dependence within the cod and especially sprat population (Casini et al., 2011).

Regarding the macrobenthos, the simulated increases in community biomass with increased organic matter supply are supported by studies from coastal OxAr of the Baltic Sea, but the balance with biomass losses due to simultaneously spreading hypoxia are not well quantified on the ecosystem scale (Cederwall and Elmgren, 1980, 1990; Karlson et al., 2002; Timmermann et al., 2012). Both models predict a change in community composition toward a stronger dominance of L. balthica with increased productivity. An uneven community composition could decrease the functional diversity contributing to efficient nutrient recycling as well as decrease the resilience of the benthic community to further environmental perturbations (Villnäs et al., 2013). Changes in macrofaunal and fish community composition could also have large effects on the magnitude of fish predation on benthos and therefore on the biomass flows from the benthic to the pelagic food web, with important consequences on marine ecosystem functioning (Griffiths et al., 2017).

On the other hand, continued efforts to reduce nutrient loads could lead to decreased primary production (Saraiva et al., 2019). This could mean less production of certain commercial fish, as was the case during the first half of the 20th century, when low nutrient availability limited the growth of the cod population (Eero et al., 2011). Additionally, it may have negative consequences for benthivorous birds (Laursen and Møller, 2014).

Our models assume a simple relationship between productivity and food availability for the studied groups, as the measures of productivity used here only represent the production of edible plankton, such as diatoms. In the last decades, increases in primary production have mostly been due to increases in cyanobacteria and late summer phytoplankton species in the pelagic, which can be of poor nutritional quality or even toxic (Raateoja et al., 2005; Suikkanen et al., 2007; Kahru et al., 2016), and in filamentous algae in coastal areas that reduce habitat quality for commercial fish and benthic infauna (Norkko and Bonsdorff, 1996; Troell et al., 2005). Additionally, higher water temperatures may increase respiratory losses during pelagic heterotrophic processing of organic matter, decreasing sedimentation (Griffiths et al., 2017; Tamelander et al., 2017). Thus, it is possible that both fish and BM community biomasses would decrease as response to lower food and/or habitat availability despite continued increases in primary production.

In conclusion, this study demonstrates the insights that can be gained through a combined use of complementary modeling approaches. Both models generally agree that nutrient loads are a stronger driver of pelagic and benthic community composition than warming in the study system. However, the effects on individual functional groups were sensitive to the number of groups and trophic interactions included in the model. Ecosystem models like EwE consider a larger number of such interactions and also include socio-ecological components. Thus, EwE is more capable to generate broad predictions on the state of ecosystem services under environmental change. However, flexible parameterization and more physiologically explicit forcing make BM more readily validated against experimental or field studies. In both models the strengths of predator-prey interactions were allowed to vary in response to the change in biomasses of the interacting groups. As we observed extreme biomass changes across the simulations ranging from almost extinction to a tenfold increase, we consider this an important added value compared to statistical models where no or only limited change is possible in interaction strengths (e.g., Lindegren et al., 2014; Blenckner et al., 2015). The models also showed multiple interactive effects of environmental drivers on different ecosystem components, adding to the evidence that multiple stressors should be considered in combination in empirical studies and environmental management.

Data Availability

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

BB, EE, and BG conceived and designed the study. EE and BB performed the model simulations, analyzed the results, and wrote the manuscript with considerable support from BG.

Funding

This study received funding from the EU, the Academy of Finland, and the Swedish Research Council Formas through BONUS (Art 185) as part of the BONUS BALTICAPP project, and from the Swedish Agency for Marine and Water Management through their grant 1.11 – Measure for Marine and Water Environment.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Prof. Alf Norkko for the fruitful discussions and comments related to this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2019.00492/full#supplementary-material

Footnotes

^www.ecopath.org

References

Atkins, J. P., Burdon, D., Elliott, M., and Gregory, A. J. (2011). Management of the marine environment: Integrating ecosystem services and societal benefits with the DPSIR framework in a systems approach. Mar. Pollut. Bull. 62, 215–226. doi: 10.1016/j.marpolbul.2010.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartolino, V., Margonski, P., Lindegren, M., Linderholm, H. W., Cardinale, M., Rayner, D., et al. (2014). Forecasting fish stock dynamics under climate change: baltic herring (Clupea harengus) as a case study. Fish. Oceanogr. 23, 258–269. doi: 10.1111/fog.12060