Reducing eutrophication increases spatial extent of communities supporting commercial fisheries: a model case study

Reducing eutrophication increases spatial extent of communities supporting commercial ﬁsheries: a model case study. In this study we investigate if eutrophication management has the potential to substantially affect which areas are going to be most suitable for commercial ﬁshing in the future. We use a spatial ecosystem model, forced by a coupled physical-biogeochemical model, to simulate the spatial distribution of functional groups within a marine ecosystem, which depends on their respective tolerances to abiotic factors, trophic interactions, and ﬁshing. We simulate the future long-term spatial developments of the community composition and their potential implications for ﬁsheries under three different nutrient management scenarios and changing climate. The three nutrient management scenarios result in contrasting developments of bottom oxygen concentrations and phytoplankton abundance, with substantial effects on ﬁsh production. Nutrient load reduction increases the spatial extent of the areas suitable for the commercially most valuable demersal ﬁsh predator and all types of ﬁsheries. This suggests that strategic planning of ﬁshery management strategies could beneﬁt from considering future changes in species distributions due to changes in eutrophication. We show that combining approaches from climate research, physical oceanography, biogeochemistry, biogeography, and trophic ecology with economical information provides a strong foundation to produce scientiﬁc knowledge that can support a multisectoral management of ecosystems.


Introduction
Eutrophication-induced habitat degradation directly affects demersal and demerso-pelagic fish, and may affect the fisheries exploiting them as well (Stortini et al., 2017;Townhill et al., 2017). Such fish commonly function as key predators in aquatic ecosystems. Therefore, changes in their spatial distribution as a V C International Council for the Exploration of the Sea 2018. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. result of management actions modifying underwater habitat quality can have large effects on the spatial distribution of their prey and the whole community. These predators and their prey are in some cases targeted by different segments of the fishery. Thus, eutrophication reduction actions may actually have different effects across the fisheries sectors. To predict such effects it is important to reliably estimate species distribution changes, which necessitates to consider not only direct effects of changes in habitat quality on commercial fish but also their biotic interactions (Godsoe et al., 2017).
Here, we use a modelling framework to investigate the causal chain between nutrient load management and the spatial distribution of fishing efforts: changing abiotic conditions affecting species distributions and fish production across space, and the latter influencing relative suitability of fishing grounds. The framework consists of a regional climate model, a coupled physical-biogeochemical model and an ecosystem model incorporating economic information, parametrized to describe the central Baltic Sea ecosystem. Eutrophication is one of the main pressures on the Baltic Sea ecosystem and the extent of hypoxic areas increased 10-fold during the past 100 years (Carstensen et al., 2014). Increased nutrient loading is proposed to have increased production of forage fish (Eero et al., 2016), but reduced the suitable habitat of the eastern Baltic cod (Gadus morhua) (Casini et al., 2016) causing a mismatch in the spatial overlap of cod and its main forage fish species, which might be one of the reasons of the failed recovery of this cod stock from overfishing (Eero et al., 2012).
Even though there is a number of models focusing on different aspects of the Baltic Sea ecosystem, there is a lack of processbased understanding of the spatial effects of environmental drivers on the whole food web. Previous studies on species distributions and pressures in the Baltic Sea (Gogina and Zettler, 2010;Casini et al., 2011Casini et al., , 2014Voss et al., 2012;Uusitalo et al., 2016;Bartolino et al., 2017) ignore dynamic feedbacks among ecosystem components. Similarly, spatial process-based models of eastern Baltic cod stock and fisheries (Röckmann et al., 2007(Röckmann et al., , 2008Kraus et al., 2008;Bastardie et al., 2010aBastardie et al., , b, 2017 have not taken trophic interactions into account so far. Radtke et al. (2013) model spatial distributions of Baltic fish based on plankton food availability, omitting direct effects of environmental drivers on fish and the benthic part of the food web. Models looking at combined effects of environmental drivers and fisheries while representing food web interactions Ö sterblom et al., 2007;Niiranen et al., 2013) lack a spatial component, with the exception of the model developed by Lindegren et al. (2014), which modelled the central Baltic Sea as three interlinked subbasins. Previous studies generally showed a link between high nutrient loads, pronounced eutrophication and an increase of sprat abundance, whereas low nutrient loads are generally thought to lead to decreased eutrophication and an increase in cod abundance. However, it is an open question how these effects are going to be realized in space and if there are areas within the Baltic Sea that are going to especially benefit from the positive effects of reduced eutrophication.
To answer this question, we use a modelling approach that goes beyond previous studies by incorporating both information on abiotic drivers of species distributions, trophic interactions, and fisheries effects on the food web in space. Ecospace is the spatialtemporal module of the commonly used Ecopath with Ecosim (EwE) suite of models (Walters et al., 1999;Pauly et al., 2000).
The newest addition to Ecospace, the habitat capacity model, combines the strength of Species Distribution Models (Peterson et al., 2011) with dynamics approaches by incorporating a dynamic niche model that considers the responses of functional groups to any number of (changing) environmental conditions (Christensen et al., 2014). In the present study we use the habitat capacity model of Ecospace to identify potential shifts in distributions of functional groups as a result of changing environmental conditions under three different nutrient management scenarios and changing climate. In addition, we are going to investigate to what extent the suitability of different areas for fishing may change under these scenarios.

Study system
The area represented in our model is the central Baltic Sea, a large brackish water body in northern Europe. Weather-driven inflows from the North Sea and anthropogenic nutrient loads from land determine oxygen concentrations (Meier et al., 2006;Matthäus et al., 2008). During the last decades, hypoxic conditions on the sea bottom have become more widespread (Figure 1), with adverse effects on the reproductive potential and stock production of demersal spawning fish and on benthic macroinvertebrates (Karlson et al., 2002;Meier et al., 2012a;Carstensen et al., 2014;Casini et al., 2016).
The offshore central Baltic Sea contains a highly productive but low diversity ecosystem with three main commercially important fish stocks, the Eastern Baltic cod, and two clupeid stocks, sprat (Sprattus sprattus) and central Baltic herring (Clupea harengus) (ICES, 2016a). Flounder (Platichthys flesus) is also a relatively abundant species and caught commercially as well. Even though the grey seal (Halychoerus grypus) population is steadily increasing, the number of grey seals is still low (Härkönen et al., 2013), thus, cod is the main piscivore. Cod, flounder, and to some extent herring, consume benthic preys while herring and sprat are the main planktivores. Mysids (mainly Mysis mixta, M. relicta, and Neomysis integer) consume both phyto-and zooplankton as well as benthic material, thus, they provide an important trophic link between the benthic and pelagic parts of the food web.

Regional ocean climate model
We use scenario simulation results of the regional ocean climate model RCO-SCOBI which consists of the physical Rossby Centre Ocean (RCO) model (Meier et al., 2003) and the Swedish Coastal and Ocean Biogeochemical (SCOBI) model (Eilola et al., 2009) performed within the project ECOSUPPORT 2009-2011 (Advanced modeling tool for scenarios of the Baltic Sea ECOsystem to SUPPORT decision making, see Meier et al., 2014).
The ocean model is coupled to a Hibler-type sea ice model and the subgrid-scale mixing in the ocean is parametrized using a k-e turbulence closure scheme with flux boundary conditions (Meier et al., 2003). A flux-corrected, monotonicity-preserving transport scheme is embedded without explicit horizontal diffusion. In the northern Kattegat open lateral boundary conditions are used, where in case of inflow temperature, salinity, and nutrient values are nudged toward observed climatological profiles. Horizontal and vertical resolutions amount to 3.7 km and 3 m, respectively.
SCOBI describes the dynamics of nitrate, ammonium, phosphate, oxygen, and hydrogen sulphide concentrations (the latter as negative oxygen), three phytoplankton species, zooplankton and Spatial effects of eutrophication on fisheries detritus (Eilola et al., 2009). The sediment contains nutrients in the form of benthic nitrogen and benthic phosphorus. Processes like assimilation, remineralization, nitrogen fixation, nitrification, denitrification, grazing, mortality, excretion, sedimentation, resuspension, and burial are considered. Resuspension of organic matter is calculated using a simplified wave model (Almroth-Rosell et al., 2011).

Downscaling of projected climate change
Atmospheric forcing fields of RCO-SCOBI were calculated applying a dynamical downscaling approach using a regional coupled atmosphere-ice-ocean model (Meier et al., 2012b) with lateral boundary data from a global climate model HadCM3 (Gordon et al., 2000). For the projections 2001-2098 the greenhouse gas emission scenario A1B was selected (Naki cenovi c et al., 2000). Bias correction of atmospheric forcing data for the ocean model was not applied, except that wind speed extremes were improved using simulated gustiness . River runoff was calculated from the net water budget over land (precipitation minus evaporation) using a statistical model (Meier et al., 2012b). Finally, nutrient loads were derived from the product of river flow and riverine nutrient concentrations. For details of the modeling approach and climate model results, the reader is referred to Meier et al. (2012b, c).

Food web model
We constructed a food web model describing the environmental drivers of the functional groups and their trophic interactions in the offshore central Baltic Sea using the EwE food web modelling approach (Walters et al., 1997;Christensen and Walters, 2004). The first component of the suite, Ecopath, describes the average trophic flows in an ecosystem during one year in our case. The Ecosim model is a set of differential equations describing the temporal behaviour of the ecosystem, using the Ecopath model as initial condition. More details on the EwE approach are included in the Supplementary material. The capabilities and limitations of the approach have been described by Christensen and Walters (2004), Plagányi and Butterworth (2004), and Plagányi (2007).
Ecospace is the spatially explicit component of EwE (Pauly et al., 2000;Christensen et al., 2014;Romagnoni et al., 2015). Ecospace is represented by a set of water and land grid cells. Functional groups and fisheries interact with each other within the water cells according to modified versions of Ecosim equations (see Supplementary Appendix S3). The representation of life histories in Ecospace compared to Ecosim is modified (Walters et al., 2010) and an effect of habitat capacity on predator-prey interactions is introduced. Low habitat capacity for a consumer species is modelled as decreased vulnerability of its prey to predation (Christensen et al., 2014). Habitat capacity in a cell for a functional group depends on the values of environmental drivers in the cell and the group's response function to these (Supplementary Appendix S3.1).
To initialize Ecospace simulations, biomasses of functional groups are distributed based on their respective overall relative habitat capacity values. These biomass distributions change in the following time-steps due to food web interactions. These biomass distributions change in the following time-steps due to the interplay of food web interactions, fishing, and species dispersal until Ecospace reaches spatial equilibrium. Therefore it is necessary to have a spin-up period under stable conditions in Ecospace, before introducing spatio-temporal forcing.
Spatial migration among cells is represented by redistributing the functional groups' biomasses among cells with a speed depending on their basal migration rate. Overall relative habitat capacity is inversely proportional to the rate of migration out of grid cells, as organisms are assumed to be less likely to leave habitats with higher capacity and more likely to migrate out of habitats with lower capacity (Christensen et al., 2014). Fishing efforts of fleets are distributed among cells based on the attractiveness of each cell for the fleet (eq. 7, Supplementary Appendix S3). Fishing mortality caused by each fleet on its target species in each cell is proportional to its fishing effort in that cell.

Model parameterization and calibration
Our Ecopath model describes annual trophic flows in the Baltic Proper during the early 2000s between 21 functional groups (composed of developmental stages of fish groups, species or groups of species) and detritus ( Figure 2). Consistency of Ecopath input parameters with basic ecological principles was checked using the Prebal procedure (Link, 2010), described in detail in ICES (2016b, Annex 3).
The Ecopath model includes the effects of fisheries on the food web by defining 10 fishing fleets operating in the region and the fishing mortality caused by them, calculated based on their landings and discards ( Figure 3). We implemented three types of fleets in the model: (1) active demersal (ACT; mostly otter trawls and demersal seine) in three size categories: <18 m, 18-24 m, 24-40 m; (2) passive demersal (PAS; gillnets, trammel nets, longlines, and pots) in three size categories: For this study, the Ecosim model described in ICES (2016b, Annex 3) was refitted to a number of reference time series using environmental forcing functions derived from RCO-SCOBI outputs, corresponding to the time period 2004-2013 (please see Supplementary Appendix S2.2 for details of the fitting procedure). The period 2004-2013 was chosen for fitting as 2004 was the first year when fishing effort (kW days at sea) data from STECF became available and 2013 the last year when an analytical stock assessment for the Eastern Baltic cod was performed (at the time of this study). Both types of information were used during the model fitting procedure. The procedure was the same as described in ICES (2016b). During the model fitting process, first we assessed the sensitivity of the sum of squared deviations measure (SS) to the number of "vulnerability blocks" (v-s) fitted using the "Stepwise fitting" plug-in of Ecosim (Christensen et al., 2008). Second, we set the v values to those maximizing model fit to time-series using the "Fit to time series" plug-in (Supplementary Appendix S2.1). As suggested by Heymans et al. (2016), we did not simply use the v-s resulting in the best fit to observed time series data, but applied additional tests on stockrecruitment and fishing mortality-catch relationships (Heymans et al., 2016;Stäbler et al., 2016) and model stability (Mackinson and Daskalov, 2007) to test for ecologically credible model behaviour and modified a few v-s accordingly (Supplementary Appendix S2.1).
To set up the Ecospace model, driver maps were generated for each environmental driver (Supplementary Table S5). All environmental driver maps we used are derived from the outputs of the RCO-SCOBI model, with the exception of the depth map. The latter is based on the Depth Relief Map published by the Temora sp.
Other ZP Spatial effects of eutrophication on fisheries HELCOM Map and Data Service (www.helcom.fi). We use yearly average phytoplankton biomass as relative primary production map, similarly to Coll et al., (2016). To parametrize environmental response functions (ERF) in the Ecospace model, we collected information about the responses of functional groups and species biomasses to abiotic factors from the species distribution modelling literature (Supplementary Table S5). We assumed three types of ERFs, "left-shoulder" (Supplementary Figure S8a), "trapezoid" (Supplementary Figure S8b), and "right-shoulder" (Supplementary Figure S8c) shapes. The choice of shape for a particular group-environmental driver pair does not reflect some general ecological characteristic of that group, rather it shows whether the environmental driver in the Baltic Sea have been described to encompass the entirety of the groups' preferred range and values above and below that (trapezoid shape) or whether the group is only possibly limited by that driver because of too high (left-shoulder) or too low (right-shoulder) values in that ecosystem.

Meiozoobenthos
In contrast to the Ecosim module, which we fit to non-spatial time series data, we assessed a fit of the Ecospace output to spatially explicit but temporally static empirical data (maps). There is no automated fitting procedure available for Ecospace. In the lack of temporal forcing, our model describes ecosystem behavior approximately of the year 2004. However, to make the model validation less sensitive to potential noise in the data and inherent natural variability in the system, we compared averaged observed stock biomass, catch and fleet effort distributions from the period 2004-2008 to model outputs. The Ecospace model validation process is described in more detail in Supplementary Appendix S3.2.

Sensitivity analysis
We tested the sensitivity of our biomass simulations to key ecological assumptions. First, we iteratively tested how excluding ERFs from the model influenced the correlation with data. This way we could identify those ERFs that were crucial to reproduce key patterns in observational data (Supplementary Appendix S3.3.1).
We also investigated the sensitivity of model fit to two parameters related to fisheries (Supplementary Appendix S3.3.2): port placement, which influences spatial distribution of fleets via the fishing cost map (Supplementary Figure S9), and Effective Power (1=r in eq. 7, Supplementary Appendix S3). We reran the model using the same settings as for the validation run, with five variations of randomly placed ports and with values for Effective Power ¼ 0.5, 1, 5 and 10.

Scenario simulations
First, we simulated three scenarios driven by differing nutrient loads using the RCO-SCOBI model. We then used environmental driver maps and temporal forcing derived from that model to drive distributions of functional groups, and, consequently, fishing efforts in Ecospace. The three scenarios of nutrient concentrations were selected to reflect rather contrasting socio-economic developments in the Baltic Sea catchment area (Meier et al., 2012b): (1)  We used annually averaged maps in EwE as drivers as we focus on the effects of long-term changes in environmental conditions and not on the seasonal cycle or extreme events like salt water inflows. The driver maps were inserted into the running Ecospace model through the spatial-temporal data framework (Steenbeek et al, 2013). We considered the same warming climate and increasing seal population in all scenarios, to be able to compare eutrophication effects in a realistic environmental context. We kept the total level of fishing efforts per fleet over the whole modelled area constant at 2013 levels. However, the spatial distribution of efforts within the area was changing every time-step as a consequence of changes in spatial distributions of the targeted fish. This means that total fishing mortality caused by each fleet on the species they catch remained constant over time, but varied in space according to the simulated effort distribution. Temporal forcing used in the scenarios is described in Supplementary Appendix S2.3.

Results
First we compare the ecosystem response among the three modelled nutrient scenarios BSAP, REF, and BAU. Second, we present the main outcomes of the sensitivity analysis.

Spatial ecosystem structure
In our EwE projections, species or groups sensitive to O 2 concentrations close to the seafloor generally benefit from reduced nutrient loads. Cod, flounder (Figure 4 deep basins (Supplementary Figure S26). Cod and flounder biomass density is low in the direct proximity of the coast and south-east from the island of Ö land due to low bottom salinity in all scenarios, and west of Gotland due to low oxygen content especially in the REF and BAU scenarios. In the BAU and REF scenarios both species are mostly concentrated in the southern, and, in the case of flounder, eastern parts of the Baltic. They (especially flounder) reach high densities along the coasts, just beyond the shallowest areas, in these two scenarios. Changes in demersal fish distributions substantially affect the spatial distributions of some of their prey and top predator species. Besides clupeids, both juvenile and adult cod and adult flounder are important prey for grey seals, and therefore seal concentration is predicted to shift southwards under both REF and BAU scenarios compared to BSAP (Supplementary Figure S25). Sprat is present in all of the modelled area in all scenarios, but under BSAP it is rather concentrated toward shallower areas (Figure 4). Under REF and especially BAU sprat has a very high density across the whole area although it is relatively more concentrated in deep areas. In both cases, the distribution of sprat is negatively related to the distribution of cod, most probably due to strong cod predation on sprat. Compared to other fish, the spatial distribution of herring is less affected by the nutrient load scenarios, even though also for this species there is a general increase in density across the whole area in the REF and BAU scenarios. This is probably due to various factors affecting its distribution simultaneously (predation by cod and seal, benthic food availability, competition with sprat for zooplankton). Spatial distributions of the intermediate trophic level predators (the clupeids) affects the distributions of lower trophic level groups. Sprat and herring are the most important predators of the Pseudocalanus spp. and "other zooplankton"' functional groups in the model, which therefore benefit from the relatively low densities of clupeids in the deep sea east of Gotland under the BSAP scenario ( Figure 5). Even though the smaller Acartia spp. and Temora sp. are also consumed by clupeids, they are significantly predated upon by mysids as well. This is probably the reason why they do not show substantial differences in their distributions among the scenarios (Supplementary Figure S25). Differences in the spatial distribution of the primary producer group among the scenarios are the result of differences in zooplankton predation and nutrient loads. Phytoplankton density overall is increasingly higher when comparing BSAP, REF, and BAU scenarios due to an increasing level of nutrients available for primary production ( Figure 5). While in the BSAP scenario phytoplankton in the deep offshore area east of Gotland is consumed by zooplankton, the low densities of Pseudocalanus spp. and the "other zooplankton" groups under REF and especially BAU result in an accumulation of phytoplankton biomass in the area. Figure 6 shows the distribution of fishing efforts of three selected fleet segments (one vessel size per each gear type) under three scenarios. Note that socio-economic drivers, such as port placement, fleet composition, and structure are assumed to be constant in time. Thus, modelled differences in fishing effort distributions across scenarios reflect differences in their target species distributions, higher priced fish having a larger influence. Thus, effort distributions indicate the relative suitability of fishing grounds under the three scenarios. Effort distributions of fleets using active and passive gears strongly reflect the biomass distribution of cod. Consequently, under the BSAP scenario their efforts are more evenly distributed over a larger area than in BAU and REF. This means that while under BSAP there are many similarly suitable fishing grounds in the model, increasing nutrient loads lead to intense fishing in small areas. Comparison of weighted center points of fishing effort distributions in 2004 to those from 2088 to 2098 shows that under BSAP fishing efforts of the demersal active and passive fisheries shift in a north-east direction, especially in the case of passive fleets. Under the REF and BAU scenarios weighted center points do not shift in space compared to 2004. Fleet effort distributions are projected to be very similar among fleet segments using differently sized vessels. Thus, the effort distributions shown in Figure 6 for mid-sized demersal trawlers and small vessels using passive gears are representative for all modelled vessel size categories of demersal trawlers and vessels using passive gears, respectively.

Distribution of fishing effort
The fishing effort distribution of the pelagic fleet segment ( Figure 6) reflects herring and sprat distribution in the BSAP scenario ( Figure 4). Although this fleet segment mostly targets clupeids, it catches cod as well. This explains our projections which indicate that under the REF and BAU scenarios the location of the most suitable fishing grounds mirror the changes in clupeids' distribution at the broad scale and the cod distribution at a finer scale. For all fleet segments, but especially for the pelagic fleets, the weighted center points during 2088-2098 are concentrated in a small area in the BSAP scenario relative to the other two scenarios, where they are more scattered. This means that under BSAP the year-to-year variability in effort distributions is smaller, indicating less change in the location of the most suitable fishing grounds between subsequent years. In contrast to the demersal fleets, effort distribution varied with vessel size in case of the pelagic fleet. In our model, vessels <24 m have a higher share of cod in their landings and therefore their distributions mostly reflect that of cod in all scenarios, similarly to demersal fleet segments. In contrast, landings of vessels >40 m consist almost entirely of clupeids and therefore their distributions follow that of the clupeids in all scenarios (Supplementary Figure S28).

Sensitivity analysis
The correlation between the modelled functional groups and the fleet effort distributions to empirical data and its sensitivity to model assumptions are described in Appendices S3.2 and S3.3, respectively. In general, model fit to observations measured by correlation is similar among biomasses, catches, and efforts (Supplementary Figure S10

Discussion and conclusions
In this study we present a mechanistic framework to assess how future nutrient management measures potentially alter the capacity of different areas of the central Baltic Sea to support commercial fisheries under climate change. We show that the Spatial effects of eutrophication on fisheries implementation of a strong nutrient reduction policy, such as the BSAP, would strongly increase the spatial extent of areas that can support all types of fisheries. On the other hand, a smaller part of the Baltic Sea may experience increased densities of fish under scenarios assuming constant or further increasing nutrient loads. Such increased densities may cause population pressures and responses that are not included in our modelling framework, such as increased parasite infection rates and decreased individual growth (Eero et al., 2015;Casini et al., 2016).
We found large differences among three modelled nutrient management scenarios in terms of spatial community composition and, consequently, fishing effort distributions across the whole modelled area. Although one region, the southeastern Baltic Sea, remained an important fishing ground in all scenarios, its relative importance compared to other areas changed dramatically. While in the highest nutrient load scenario it was the only area which could sustain both demersal and pelagic fisheries, in the lowest nutrient load scenario other, more northern areas also became suitable. Therefore, the relative location of most suitable fishing grounds for demersal fisheries shifted northwards, especially for the segments using passive gears. An extended potential range of operations may be particularly important for this segment as it is considered to be the most vulnerable within the Baltic fishery (Strehlow, 2010). Not only the spatial distribution of suitable fishing grounds, but also their interannual variability, differed among the scenarios. Under the low nutrient loading scenario, larger areas were suitable for fisheries and their location tended also to be more stable among years. This sort of spatial reliability of fish production may facilitate the inclusion of fisheries into marine spatial planning in the future.
One of the most important outcomes of our study is that the differences in species distributions among modelled scenarios were the result of cumulative impacts of several environmental factors, in agreement with Stortini et al. (2017). While for individual groups one or two factors could be pinpointed as important drivers, changes in the spatial structure of the community as a whole were the result of the combined effects of changes in oxygen, salinity, primary productivity, and food web interactions. In the Baltic Sea, currently cod is the most important top predator and changes in its abundance potentially cause multilevel trophic cascades (Casini et al., 2008). Hypoxia-induced habitat compression of cod and its consequences for the spatial distribution of intermediate trophic level forage fish in the Baltic Sea are well documented (Casini et al., 2011(Casini et al., , 2016. Our model results indicate that the habitat compression of cod may be reversed if nutrient load reduction policies are implemented. While constant and high nutrient load scenarios had adverse effects on benthic and demersal groups, abundances of phytoplankton and pelagic fish were predicted to increase. Similarly to other seas, there has been a positive link between increased nutrient loads and (especially forage) fish production in the Baltic (Chassot et al., 2007;Eero et al., 2016) also supported by our model results. However, it is questionable whether this relationship will hold in the future. Some evidence suggests that further increases in the eutrophication levels compared to today, especially under higher temperatures and lower abundances of higher trophic level predators, could contribute to shifts in primary producer composition to an unfavourable state for consumers. Such shifts include an increased proportion of smaller-sized organisms (Suikkanen et al., 2013), a more frequent occurrence of toxic cyanobacterial blooms (Lehtiniemi et al., 2002;Neumann et al., 2012), and an increased dominance of filamentous algae in coastal habitats (Borg et al., 1997).
Our results also point out the environmental dependency of suitable areas for fisheries and possibly all human activities based on ecosystem functioning. This means that long-term, adaptive marine spatial planning needs to take into account changing abiotic conditions (Miller et al., 2013). Our modelling study suggests that the provision of wild-captured fish food, one of the important ecosystem services, may have a more even spatial distribution across the central Baltic Sea when nutrient loads are reduced. This could have important economic consequences for the fishing industry as spatial relation to the most productive fishing grounds is an important determinant of fleet efficiency (Hutniczak et al., 2015;Bastardie et al., 2017). When fish distribution consists of small pockets of high densities in space, as predicted under increasing nutrient loads, the risk for overexploitation is higher. Discard issues may also increase if species which are targeted and those that are caught as bycatch have similar requirements and their distributions become restricted to overlapping areas, such as cod and flounder in our model. Additionally, fisheries have to share the marine space with other human activities (Tidd et al., 2015;Yates et al., 2015). For example, in one area within the Sound (part of the Baltic Sea) a trawling ban has been in place since 1932 due to intense shipping traffic in the area (Lindegren et al., 2013). When the extent of areas suitable for Spatial effects of eutrophication on fisheries fishing operations is decreased, together with the extent of areas supporting ecosystem functions, managers may face more difficult trade-offs in allocating marine areas for exploitation, conservation and other uses.
In their recent study, Zurell et al. (2016) have shown that mechanistic modelling approaches, such as dispersal or population dynamics models and Bayesian process-based dynamic range models, outperform correlative species distribution models in predicting species range dynamics under climate change. We argue that the approach presented here is a useful complement to those evaluated by Zurell et al. (2016), as it simultaneously provides projections of all functional groups in an ecosystem without necessarily needing spatio-temporal data on abundances of all groups. For our ecosystem, the model was also not sensitive to the number of ERFs included and major patterns in the data could be reproduced by including a few key functions only (see Supplementary Appendix 3.3.1). Still, there is a need for the development of a consistent methodology for the parameterization of ERFs that express the responses of functional groups to abiotic factors. The empirical measurement of such responses is a highly active research area (e.g. Birnie-Gauvin et al., 2017). In addition, developing standard methodology to reliably assess the skill of spatial ecosystem models such as Ecospace is important to have an insight about the uncertainty of their predictions. Ideally, such a methodology would be based on a combination of metrics including correlation as used here, but also neighbourhood-based methods as described by Rose et al., (2009) and Stow et al., (2009) for oceanographic models.
As Ecospace model parameterization is not based on automated statistical fitting but on expert judgement and literature values, it is especially important to explore the sensitivity of the results to assumptions made during model parameterization. Romagnoni et al. (2015) conducted an extensive sensitivity analysis of their Ecospace model for the North Sea. We have tested our model's sensitivity to some of the same parameters they have found to be important. Both studies found a reasonably good agreement between modelled population distributions and spatial data from scientific surveys and a large effect of the parameter "Effective Power," which affects the level of dispersion of modelled fleet efforts around profitable fishing areas. The agreement between modelled fishing efforts and spatial data from commercial fisheries was better in our Baltic model. In contrast to Romagnoni et al. (2015), we found an effect of port placement on fishing fleet distributions. The placement of fishing ports affects the calculation of a fishing cost map that reflects distance from ports. The fishing cost map is then used to distribute fishing effort, evaluating fleet-and cell-specific fishing costs based on the fleet-specific ratio of sailing-related costs to fixed fishing costs. The latter ratios were much higher in the case of the Baltic model which explains the higher sensitivity of our modelled fleets' to port placement. Notably, some species distributions were also sensitive to port placement (see Supplementary Figure S21). The reason for this is that port placement influenced how fishing mortality was distributed in space via making areas far away from ports relatively less attractive for fishing fleets. This underlines the importance of considering both economic and environmental factors when making predictions about future species distributions.
Compared to other modelled populations, our approach proved to be less successful in reproducing the distribution of Figure 6. Projected fishing effort (average values 2088-2098) of selected fishing fleets, up to down: mid-sized demersal trawlers, small vessels using passive gears, and mid-sized pelagic trawlers, under three nutrient management scenarios: Baltic Sea Action Plan (left), Reference (middle) and Business-As-Usual (right column), in the modelled area (see Figure 1). Values express fishing efforts relative to each fleet's average effort over the entire modelled area in the initial year, 2004. Darker shades represent higher values. Brown triangles indicate the locations of the modelled weighted center points of the effort distributions in each of the last 11 simulated years (2088)(2089)(2090)(2091)(2092)(2093)(2094)(2095)(2096)(2097)(2098). Orange circles show the same from 2004 (initial model state after spin-up period).
one commercially important group, herring. Some earlier studies have shown that modelled distributions of species with distinct environmental preferences, such as cod, generally fit better to data than those of species with wide tolerances, such as herring (Somodi et al., 2017). In addition, pelagic species have more variable distributions than demersal ones which is harder to reproduce by models (Thorson et al., 2016). These results suggest that spatial management of such groups inevitably involves more uncertainty.
Changes in habitat quality due to human impacts are increasingly common across the globe. As species shift their distributions to adapt to altered environmental conditions, the spatial provision of ecosystem services changes as well. Here we presented the projected effects of various nutrient management policies on various environmental variables and the cumulative effects of those factors across the marine food web and on commercial fisheries in the example of the Baltic Sea. Where data are available, the same approach could be used to evaluate the potential consequences of various environmental policies in other systems. In the Baltic Sea, it may provide inspiration for studies more focused on certain functional groups or areas. Our results indicate the effectivity of nutrient load reduction policies in recovering ecosystem function across large areas of the Baltic Sea, which may motivate environmental managers to further pursue such policies.

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.