Effects of Nutrient Management Scenarios on Marine Food Webs: A Pan-European Assessment in Support of the Marine Strategy Framework Directive

Eutrophication is one of the most important anthropogenic pressures impacting coastal seas. In Europe, several legislations and management measures have been implemented to halt nutrient overloading in marine ecosystems. This study evaluates the impact of freshwater nutrient control measures on higher trophic levels (HTL) in European marine ecosystems following descriptors and criteria as defined by the Marine Strategy Framework Directive (MSFD). We used a novel pan-European marine modeling ensemble of fourteen HTL models, covering almost all the EU seas, under two nutrient management scenarios. Results from our projections suggest that the proposed nutrient reduction measures may not have a significant impact on the structure and function of European marine ecosystems. Among the assessed criteria, the spawning stock biomass of commercially important fish stocks and the biomass of small pelagic fishes would be the most impacted, albeit with values lower than 2.5%. For the other criteria/indicators, such as species diversity and trophic level indicators, the impact was lower. The Black Sea and the North-East Atlantic were the most negatively impacted regions, while the Baltic Sea was the only region showing signs of improvement. Coastal and shelf areas were more sensitive to environmental changes than large regional and sub-regional ecosystems that also include open seas. This is the first pan-European multi-model comparison study used to assess the impacts of land-based measures on marine and coastal European ecosystems through a set of selected ecological indicators. Since anthropogenic pressures are expanding apace in the marine environment and policy makers need to use rapid and effective policy measures for fast-changing environments, this modeling framework is an essential asset in supporting and guiding EU policy needs and decisions.


INTRODUCTION
Eutrophication is one of the most important anthropogenic pressures on coastal and estuarine waters (Cloern, 2001;Desmit et al., 2018), seriously threatening the functioning and structure of marine ecosystems (Diaz and Rosenberg, 2008;Nixon, 2009;Doney, 2010;Cai et al., 2011), due to an excessive amount of nutrients from agricultural run-off and sewage. In Europe, several pieces of legislation have been implemented to prevent negative effects of eutrophication, either directly with the Urban Wastewater Treatment (UWWTD CEC, 1991), the Nitrates (ND, CEC, 1991), and the Water Framework (WFD 2000/60/EC) Directives, or within an ecosystem context with the Marine Strategy Framework Directive (MSFD 2008/56/EC). Eutrophication is one of the eleven qualitative descriptors of the MSFD established by the European Commission (EC) to assess the environmental status of EU marine waters (European Commission, 2008;Cardoso et al., 2010). Despite the legislation to mitigate the negative impacts of nutrient discharge in European waters, countries have not been compliant and the measures have not been effective enough to achieve the ultimate goals of the regulations. One of the main obstacles is that marine and coastal ecosystems receiving these impacts are so complex that they prevent stakeholders from observing pressureimpact relationships between nutrient discharges and marine ecosystem state.
Ecological models are powerful tools to address the complexity of these systems and highlight these relationships (Hyder et al., 2015;Lynam et al., 2016;Heymans et al., 2018). Decades of experience and acquired knowledge resulted in a strong progress of ecological modeling, allowing e.g., to better simulate the different components of the marine environment and explore the ecological responses that might occur if alternative management scenarios were implemented (IPBES, 2016;Zandersen et al., 2019). Within this context, the scientific community has been working to build a robust and reliable "End to End Models" (E2EMs) framework, which simulate the main processes that influence the dynamics of marine ecosystems (Fulton, 2010). This framework includes different types of spatially temporally explicit models, as (1) Hydrological models: providing information on river discharge in terms of flow and nutrients; (2) Hydrodynamic models: simulating marine water transport; (3) Biogeochemical models (lower trophic level, LTL): assessing transported nutrients and nutrient processes within phytoplankton/zooplankton; and (4) Food-web/multispecies models (or higher trophic level, HTL): simulating biomass dynamics, the distribution of marine organisms (from phytoplankton to top predators) and fisheries.
Several studies have used the E2EM framework to assess the impact of nutrient reduction or other stressors on HTL organisms and ecosystem functioning (Libralato and Solidoro, 2009;Rose et al., 2010;Peck et al., 2018). In this study we explore the impact of nutrient management scenarios at the European Sea scale by running an ensemble of different HTL models, forced by an existing coupled hydrological and hydrodynamicbiogeochemical framework (see details below) covering specific areas of European regional seas. This was done to achieve a pan-European assessment of the impact of such measures in marine ecosystems while evaluating consistencies or divergences in predictions among different types of models and applications.
In particular we assess how changes in nutrient inputs and concentrations and consequently planktonic groups might impact the structure and function of the upper trophic levels of the food web. Classical food web theory suggests that nutrient enrichment affects the food web from the bottomup along with top-down effects, through predation, controlling the biomass of all trophic levels of a system (Oksanen et al., 1981;Borer et al., 2006). Several studies have investigated these synergetic effects on terrestrial and marine ecosystems (Worm et al., 2002;Isbell et al., 2013) highlighting that a decrease in nutrients would decrease the biomass of autotrophic organisms while the possible resulting decrease in the biomass of herbivores and carnivores would depend on the complexity (trophic linkages: Abrams, 1993;Leibold et al., 1997;Ward and McCann, 2017) and on the nature of the ecosystem (e.g., at mature or developmental stage and/or in oligotrophic or eutrophic conditions: Odum, 1969;Proulx and Mazumder, 1998;Schlenger et al., 2019). However, it is still unclear how biodiversity changes in relation to nutrients and which type of relationship this might be (Waide et al., 1999;Kondoh, 2001;Jara et al., 2006;Groendahl and Fink, 2017). Recent empirical and modeling studies have also suggested that net primary production is a key factor explaining fisheries yield and limiting fishery production potential (e.g., Stock et al., 2017;Link and Watson, 2019).
Our assessment follows MSFD descriptors (mainly the biodiversity related descriptors) and methodological standards, developed and agreed upon in the framework of European or international conventions (European Commission, 2017). We included additional indices, as suggested by previous studies (Shannon et al., 2014;Coll et al., 2016), to complement the assessment of the marine environment.
The use of an ensemble of models is crucial to increase the reliability of model predictions, account for prediction uncertainty and better inform decision-makers about the range of effects of selected pressures/measures on biodiversity, ecosystems and their services in general (Gårdmark et al., 2013;Maxwell et al., 2015;Bauer et al., 2019;Lotze et al., 2019). The overall aim of this paper is to utilize the pan-European model ensemble to address the impact of eutrophication on the European seas. Similar exercises have been done on regional scale (Bauer et al., 2018), but to our knowledge, this is the first pan-European multi-model comparison used to assess the impacts of land-based management measures (in this case the reduction of nutrients) on marine and coastal HTL ecosystems.

Hydrological and Hydrodynamic-Biogeochemical Models
The hydrodynamic-biogeochemical dynamics of the European seas resulting from different river scenarios have been simulated using an End to End Model called Modeling Framework (MF), developed at the European Commission by its science Directorate-General (DG), the Joint Research Centre (JRC). The MF consists of coupled (either offline or online) hydrological, hydrodynamic-biogeochemical, and food-web models (Supplementary Figure S1). These models have been implemented and validated at different spatial (regional and sub-regional) and temporal (past and future) scales (Garcia-Gorriz et al., 2016;Macias et al., 2018b) across Europe.
The MF has been designed to simulate changes in the state of European marine ecosystems and derived services in response to different pressures and management scenarios with the overall goal of providing explicit support to the decisionmaking process. In particular, in relation to eutrophication, the MF has investigated two realistic nutrient management scenarios, following measures reported and suggested by Member States within the WFD implementation plans. The two scenarios covering inland water quantity and quality (nutrients) in Europe include: (1) actual nutrient loads from river discharge (reference scenario, REF) and (2) maximum technically feasible reduction (MTFR scenario) of nutrient input to surface water. The nutrient reductions can be achieved by, e.g., keeping the nutrient surplus in agricultural areas to a minimum, optimizing mineral fertilizer applications and upgrading waste water treatments to the highest level of nutrient removal (Grizzetti et al., 2021).
The water flow and the effectiveness of measures for preventing water scarcity were simulated by the LISFLOOD model (De Roo et al., 2020); a GIS-based hydrological rainfallrunoff-routing model used to simulate the hydrological processes that occur in a river basin. Annual total nitrogen and total phosphorus loads reaching the sea were estimated by the GREEN model (Geospatial Regression Equation for European Nutrient losses; Grizzetti et al., 2012); a statistical regression model that represents the processes of nutrient transport and retention in the river basin as well as the nutrient sources and physical characteristics that influence nutrient processes. GREEN was coupled to LISFLOOD by incorporating the modeled water flow, to simulate annual nutrient loads for the period 2005-2012, corresponding to the most recent spatial data homogeneously available at the European scale (Grizzetti et al., 2021).
The MTFR scenario comprised increased water use efficiency in irrigation and domestic usage, changes in cooling water requirements and the implementation of wastewater re-use for irrigation (De Roo et al., 2020). Measures to reduce nutrient pollution in the water consisted of upgrading all waste water treatment plans in the EU to a high level of nutrient removal (i.e., a tertiary treatment with an enhanced reduction of phosphorus) and lowering the mineral fertilization in agricultural fields by setting maximum nitrogen surplus in agricultural areas to 10%, but without changing the current level of livestock and manure production (Grizzetti et al., 2021). The simulation of nutrient inputs and measures was implemented in a high spatial resolution grid (catchments of 7 km 2 average size). The simulated annual river flow and nutrient concentrations estimated by GREEN, for both scenarios, were then coupled to the LTL module of the MF to assess how changes in nutrient load might impact nutrient and plankton concentrations at sea.
The LTL module consisted of 3D hydrodynamicbiogeochemical models representing the four main MSFD regions (Mediterranean Sea, Black Sea, Baltic Sea and North East Atlantic). Details of the hydrodynamic-biogeochemical models can be found in Garcia-Gorriz et al. (2016) and Macias et al. (2018b). The mean difference in nutrient concentrations between the two scenarios (MTFR and REF) estimated by the hydrological model for the land-sea interface, and integrated for the assessed period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), ranged from −2 to −31.8% for riverine total nitrogen and from −4 to −46.3% for total phosphorus, depending on the assessed ecosystem ( Table 1). Relative changes in nutrient concentrations and primary production at sea estimated by the hydrodynamic-biochemical models were smaller but varied according to the HTL marine ecosystem considered.
In addition, an example of the spatial output scenarios produced by the hydrological and hydrodynamicbiogeochemical models for the Mediterranean Sea, and integrated for the assessed period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), is presented in Supplementary Figure S2. The full Pan-European assessment for the hydrological and LTL modules can be found in Grizzetti et al. (2021) and Friedland et al. (2021), respectively.

HTL Models
Fourteen HTL models were used to run the nutrient management scenarios, all covering either full MSFD regions, sub-regions, or smaller zones within single MSFD areas (Figure 1 and Supplementary Table S1). The criteria used to select these models were that they were validated, fitted to observed data, and published in peer-reviewed scientific literature. The model types included in this assessment were: food-webs (79%, Ecopath with Ecosim or EwE: Christensen and Walters, 2004), multispecies/individual-based (14%, Osmose: Shin et al., 2005) and end-to-end (7%, Atlantis: Fulton et al., 2011) models, which corresponded to the main modeling tools used in Europe to assess HTL compartments of marine ecosystems (Piroddi et al., 2015). A detailed description of these models is given in Supplementary  Table S1. As shown in Supplementary Table S1, the HTL models were all different in complexity and structure depending on their specific goals e.g., number of functional groups, temporal/spatial scale, and environmental/anthropogenic drivers (e.g., fishing pressure, primary production, temperature). Yet, a commonality across these tools was the ability to capture historic (hindcast) ecosystem dynamics using environmental variables (e.g., nutrients concentrations, phytoplankton biomass or primary production) as drivers. This was an essential prerequisite for exploring the effect of nutrient scenarios on these ecosystems.
To run the nutrient simulations, each HTL model was forced with specific environmental output from the REF and MTFR scenarios of the MF hydrodynamic-biogeochemical module covering the 2005-2012 period. Each scenario was run separately by the individual HTL models using the 8-years (2005-2012) simulations from the LTL models as "forecast" scenarios (Supplementary Table S1), keeping the hindcast from the original settings, as the main goal of this study was to detect changes between the two scenarios. The types of drivers provided correspond to the drivers utilized by the specific HTL models to hindcast the historical observations. For example, simulated phytoplankton biomass or primary production were used to simulate nutrients in the EwE models while simulated nutrients at sea were used in the Atlantis model (Supplementary Table S1). The other forcing (e.g., fishing mortality/effort, temperature) already incorporated in each HTL model were kept unchanged as the purpose of the exercise was to assess changes in the HTL ecosystems impacted solely by changes in nutrient management (Supplementary Table S1).
TABLE 1 | The HTL European marine ecosystems, HTL model type, acronym and HTL spatial extent (ranges of latitude and longitude), which was used to extract LTL models outputs (details in Supplementary Table S1). Relative changes (%) between the two scenarios (MTFR and REF) of riverine total nitrogen (TN) and total phosphorus (TP) loads estimated by the hydrological model (Grizzetti et al., 2021), and total nitrate (DIN), phosphate (DIP) and primary production (PP) at sea estimated by the hydrodynamic-biogeochemical models (Friedland et al., 2021)

MSFD Criteria
To assess the impact of the nutrient management scenarios, a set of criteria from MSFD Descriptor 3 (D3: Commercially exploited species) and Descriptor 4 (D4: Food webs) were used. In particular, we chose criteria and species that would likely have a direct response to these scenarios, such as small pelagic fishes (e.g., herrings, sardines, anchovies) bottom-up controlled by primary production, and criteria that would be able to capture changes within the food web (e.g., species diversity). To be able to compare output across regions and models, we selected criteria common to most of the models ( Table 2). MSFD Descriptor 1 (D1: biodiversity) was not used because the majority of the models were not able to capture dynamics of non-commercial iconic species such as marine mammals and seabirds. The full list of criteria is shown in Table 2.

Trophic-Based Indicators
In addition to the MSFD criteria, two other indicators were assessed to test if nutrient scenarios would have an impact on the trophic structure of the ecosystems. These include the trophic level (TL) of the community, excluding TL < 2 (mTLco) (Shannon et al., 2014), calculated as: where TB is total biomass, B i is the biomass of species i and TL i is the trophic level of species i (note: TB, B i and TL i vary in time); and the TL of the landed catches (TL c ) (Christensen, 1996;Pauly et al., 1998), calculated as: where Y L is total landings, Y i is the landing of species i and TL i is the trophic level of species i (note: Y L , Y i and TL i vary in time).
Modeled criteria and indicators (I) were extracted annually (for the 8 years of simulations) for each HTL model and for both REF and MTFR scenarios. The relative mean change between these scenarios was calculated as: and presented per descriptor-criterion/TL indicator and per model. For D3 the scale of assessment required by the GES (Good Environmental Status) Decision (European Commission, 2017) is at stock level (typically at one or more geographical subareas [GSAs 1 /ICES 2 ] as defined by FAO) and for D4 it is at MSFD regional/sub-regional scale. As the models in this assessment had different spatial scales (Figure 1 and Table 1) depending on their final goals, the estimation of the mean change of the selected criteria was done using the original scale of each model, not following, thus, GES requirements. The same scale was also applied for trophic level indicators.

Metrics Over Space
Only six of the HTL models used in this assessment had a spatial component; of these, three were built for the Mediterranean Sea, two covering the entire basin and one set up for the Western Mediterranean Sea (Supplementary Table S1). For the purpose of comparing the spatial predictions from the models and assessing their consistencies/divergences, the Western Mediterranean Sea was used as a reference study. In particular, the relative difference between the two scenarios for each model was calculated as the relative mean change per grid cell. Food-web criteria and their relative changes were presented per model for the whole sub-basin, for shelf (<200 m) and open water (>200 m) areas; the same scale was also applied for trophic level indicators. For criterion D3C2 (the spawning stock biomass of the European pilchard, Sardina pilchardus), the relative mean change was presented at geographic subarea (GSAs) level. D4C1 and D3C2 were not available for the Osmose model.
Finally, coherence maps were created for the three indicators (the small pelagic fishes [D4C2] and the TLs) common to the models to evaluate the coherence of the projections. Trends of relative changes were compared per grid cell and per model, looking at the signs, indicating whether an increase (or decrease) in the selected indicator occurred under the MTFR scenario. The percentage of coherence was calculated for the whole sub-basin, the shelves (<200 m) and open waters (>200 m).

MSFD Criteria-Mean Change
The scenario outputs showed differences between and commonalities among models, criteria, size and locations (Figure 2). In particular, at MSFD regional scale, the Spawning Stock Biomass [SSB] of commercial small pelagic fishes (D3C2) showed a slight decrease in all the areas, with the exception of the English Channel (Osmose; +0.7%), with average values ranging between −1.3 and −2.0% (Figure 2A). The Baltic Proper (the only model available for the Baltic Sea used to represent the entire MSFD region) showed a decrease but also a high level of variability in the modeled SSB values. The Mediterranean Sea models highlighted a reduction in the SSB in all the models/areas considered (from sub-regional to localized coastal/shelf level). Overall, the most impacted MSFD sub-region was the Adriatic Sea (−3.0%), and among smaller areas within an MSFD region/sub-region, the North Sea (−5.5%) and the West Coast of Scotland (−3.3%).
Species diversity (D4C1) did not show a significant change ( Figure 2B) among the MSFD regions except for the Baltic Sea where a slight increase of +1.3% was observed, and for the North-East Adriatic Sea where a decrease (−3.1%) was found. Small pelagic fish biomass (D4C2) had the highest variability among the areas/models (Figure 2C). At MSFD regional scale, this criterion decreased (−2.3%) in the Black Sea followed by the North East Atlantic and the Mediterranean Sea with slight decreases between −1.8 and −0.9%. In contrast, the The diversity (species composition and their relative abundance) of the trophic guild is not adversely affected due to anthropogenic pressures Species diversity index Expresses species diversity by considering the biomass of those organisms with trophic levels 3 or higher (Kempton and Taylor, 1976;Ainsworth and Pitcher, 2006). This index is calculated as follow (Ainsworth and Pitcher, 2006): where Fg is the total number of functional groups in the model, R 1 and R 2 are the representative biomass values of the 10th and 90th percentiles in the cumulative abundance distribution.

C2
The balance of total abundance (number of individuals or biomass in tons) between the trophic guilds is not adversely affected due to anthropogenic pressures   Table 1). Yearly modeled data points are plotted as colored circles. D3C2 is not displayed for Med-Osmose and NE_Adri as not available, same for D4C1 and Med-Osmose.
biomass of small pelagic fish in the Baltic Sea increased slightly (+2.0%). At sub-regional level, the Adriatic Sea from the MF-Mediterranean model (Adri_JRC) was the area with the highest reduction (−2.0%). Conversely, the Osmose model for the Adriatic Sea (Adri_OSM) projected an increase (+1.3%). For the other Mediterranean sub-regions (the Western, the Ionian and the Eastern) the available models agreed on the marginal impact of the nutrient scenarios for these regions, with an agreement in the mean change for the Western Mediterranean Sea [JRC_EwE (−0.015%); ICM_EwE (−0.011%) and Osmose (−0.017%)] and a disagreement in trend for the other two sub-regions [Ionian: JRC_EwE (−0.1%), and Osmose (0.5%); Eastern: JRC_EwE (−0.09%), and Osmose (0.86%)]. Regarding ecosystems at smaller scales within an MSFD region/subregion, in the Mediterranean Sea, the North-East Adriatic Sea showed the highest reduction (−8.4%) in small pelagic fish biomass (D4C2) followed by the Thermaikos Gulf (−1.7%) and the Inner Ionian Sea (−1.2%). In the North-East Atlantic, the North Sea was the area most impacted by the nutrient reduction with a −7.5% decrease followed by the West Coast of Scotland (−2.5%), the Celtic Sea (−1.05%) and the Irish Sea (−0.67%), while for the English Channel, the two available models highlighted that the change was minimal (Atlantis: −0.08%; Osmose: +0.04%).

TL Indicators-Mean Change
Overall, the TL indicators did not show clear changes between the two nutrient scenarios (Figure 3). For the MSFD regions, both indicators showed a mean change very close to zero. At sub-regional level, some small increases (<0.5%) were observed through the Mediterranean Osmose model and, for the mTLco, also through the Western Mediterranean Sea from ICM (West_ICM) (Figure 3A). Looking at mTLco, the highest negative changes were observed at smaller scale with a reduction in the North-East Adriatic (−0.7%) followed by the North East Atlantic models (West Coast of Scotland (−0.41%), Celtic Sea (−0.36%) and English Channel in case of the Osmose model (English Channel_OSM; −0.3%) (Figure 3A). Similarly, small reductions in TLc were observed in the North-East Adriatic (−0.2%) and  Table 1). Yearly modeled data points are plotted as colored circles.
two areas of the North-East Atlantic Sea, the West Coast of Scotland (−0.29%) and the English Channel_OSM (−0.24%) ( Figure 3B).

MSFD Criteria Over Space
The spatial outputs produced by the HTL models for the two scenarios in the Western Mediterranean Sea highlighted differences and commonalities among models, criteria/indicators, and between the whole subregion, shelf areas and open waters (Figures 4, 5 and Supplementary Figure S3). SSB (D3C2) results showed a slight reduction for European pilchard with values that, depending on the model considered, ranged between −0.07 and −6% (Figure 4). Looking at GSAs (refer to Figure 4 for the names of the different GSA regions), the highest reductions were found in the Balearic Islands (GSA5: −3.1%) and Algeria (GSA4: −2.2%) using the Western Mediterranean Sea from the JRC (West_JRC) model and the Gulf of Lion (GSA7: −0.39%) from the West_ICM model (Supplementary Table S2). Figure S3) showed good agreement in the direction of change but differences in the amplitude. Only a marginal reduction (between −0.1 and −2%) was depicted by the available models around the continental shelf of the Gulf of Lion, Northern Spain and the Balearic islands and a slight increase along the coasts of the Gulf of Lion and the North Tyrrhenian Sea. For the small pelagic fish biomass (D4C2), two out of the three available models showed a decrease along European coastlines with the highest negative values (around −3%) concentrated around the Gulf of Lion, Northern Spain and the North Tyrrhenian region (Figure 5). In these models the continental shelves were the most negatively impacted areas. A different pattern was found using the Western Mediterranean Sea from the Osmose (West_OSM) model, where high variability was observed in the whole sub-region with open waters more impacted than shelf areas (Figure 5).

The species diversity index (D4C1) (Supplementary
The spatial coherence for the small pelagic fish biomass was heterogeneous for the whole sub-region with an approximate equal percentage of decrease (45%) and increase (55%) (Figure 6 and Supplementary Figure S6). When looking at shelves and open waters, all or most models showed that 69% of the grid cells were associated with a decrease (shelves) and 61% associated to an increase (open waters) (Figure 6 and Supplementary Figure S6).

TL Indicators-Spatial Scale
The West_OSM and West_ICM models suggested a slight increase of both trophic level indicators (Figure 7 and Supplementary Figure S4) in all the assessed areas (whole region, shelf, open waters) with high variability in both models for the TL of the community (mTLco; between −4 and +10%) (Figure 7) and for TL of the catches in the Osmose model (mTLc; between −6% and +4%) (Supplementary Figure S4). Conversely, the West_JRC model showed a slight decrease in mTLco and mTLc particularly in the shelf areas (mTLco: ∼−0.011%; mTLc: ∼ −0.016%) and around the coastal areas of the Gulf of Lion, Northern Spain and the North Tyrrhenian Sea (mTLco: ∼ −0.3%; mTLc: ∼ −0.2%).
The spatial coherence for the mTLco was extremely heterogeneous for all the areas (the whole sub-region, shelves and open waters). Here, 55% of the grid cells were associated with a decrease in mTLco while 45% were associated with an increase (Supplementary Figures S5, S6). In all models, ca. 70% of the cell were associated with a decrease in TLc (Supplementary  Figures S5, S6).

Mean Change in MSFD Criteria and TL Indicators
This study provides a first pan-European assessment of the impact of nutrient management scenarios on marine ecosystems FIGURE 4 | Maps representing the mean change (%) for the Spawning stock biomass of a commercially important small pelagic fish (D3C2) per GSA (# 1. Northern Alboran Sea; 2. Alboran Island; 3. Southern Alboran Sea; 4. Algeria; 5. Balearic Islands; 6. Northern Spain; 7. Gulf of Lion; 8. Corsica; 9. Ligurian and Tyrrhenian Seas; 10. South and Central Tyrrhenian; 11.1. Sardinia West; 11.2. Sardinia East; 12. Northern Tunisia) and per model type.
and related marine resources. The reduction of nutrients from river run-off showed no substantial changes in the structure and function of the HTL ecosystem models included. From a regional MSFD perspective, the mean change of SSB of commercially important small pelagic fish species (D3C2) and small pelagic fish biomass (D4C2) showed the highest decrease if only with values below 2.5%. Interestingly, all the available models confirmed a decline in the biomass of commercial small pelagic fish species among the main MSFD regions, suggesting that a reduced primary production as a result of nutrient reductions (Table 1), might negatively influence the dynamics of these fish stocks, although only modestly. This phenomenon has already been observed in other systems as shown by Breitburg et al. (2009) andde Mutsert et al. (2016). Population dynamics of small pelagic fishes are tightly coupled to the dynamics of their food; i.e., plankton. If nutrient levels do not lead to strong adverse environmental impacts such as anoxia, harmful algal blooms and shifts in the zooplankton community, then a reduction in nutrients eventually leads to a reduction in food availability for small pelagic fishes. When fishing pressure (top down effect) on small pelagics is combined with this reduced resource availability (bottom up effect), then a decline of small pelagic fishes can be expected (Ramírez et al., 2018). Biomass variability of small pelagic species may be linked to overexploitation, climate change and other environmental factors. Because of their rapid growth and short lifespan, small pelagic fishes, and especially their recruitment success (Brosset et al., 2017), are vulnerable to climate and environmental forcing Saraux et al., 2019;Tsikliras et al., 2019). However, the potential decrease in SSB of commercially important stocks might be largely compensated by a decrease in fishing mortality (Froese et al., 2018), assuming low predator pressure, and improved size selectivity of the fisheries, notably in the Mediterranean Sea (e.g., Colloca et al., 2013).
When looking at the biomass of commercial stocks and non-commercial small pelagic fish species, the pattern was  similar to the mean change of SSB of commercially important stocks (D3C2), but with larger variability across regions/subregions/small areas within an MSFD region. Variability across models might be related to the structure of the available models (e.g., some models might have more commercially important small pelagics than non-commercial species, or vice versa). Within each model, the higher variability for small pelagic fish biomass could be driven either by the different strength of species responses to lower productivity, and/or by the temporal variability of the responses. For example, in the West Coast of Scotland, this study showed that the biomass of small pelagic fishes (herring, sprat and horse mackerel) declined with lower ecosystem productivity (bottom-up control). However, sprat, which has a relatively healthy stock and a high turnover rate in this ecosystem, quickly increased again after the initial decline due to an overall reduction in predation pressure (topdown controls). Despite differences in the way models were constructed, which could cause fundamental differences when comparing ecosystems Heymans et al., 2014), the majority of the ecosystem models projected decreases in small pelagic fish biomass. Hence, the differences in the structures of the models utilized in this study did not prevent diagnosing the most likely direction of change in European marine ecosystems under the nutrient reduction scenario.
For species diversity (D4C1), no clear responses were observed at a regional or sub-regional scale, except in the Baltic Sea with a slight increase and in the North-East Adriatic Sea with a decrease in diversity. The Baltic Proper was the only area that showed a slight increase in the two food-web criteria compared to the reference scenario; this might be related to the highly eutrophic nature of this ecosystem and the fact that a reduction of nutrient inputs might lead to an improvement of the marine environment e.g., better bottom oxygen levels, as observed in Saraiva et al. (2019) and Friedland et al. (2021), and thus better spawning conditions, leading to an increase in the Eastern Baltic cod stock. This response was not observed in other ecosystems that are also considered to be eutrophic such as the North-East Adriatic Sea, the Black Sea and the whole Adriatic Sea (from the JRC model) where nutrients reduction, mainly from the Po and Danube rivers, would reduce the assessed D3 and D4 criteria (Figure 2). The reduction in dissolved nutrients and primary production indicated by the biogeochemical models (Table 1) reflect these differences. In these systems the decrease in nutrient loads resulted in an average reduction in the primary production of 0.3% for the Baltic Sea, 3.4% for the Black Sea and 4.8% for the North-East Adriatic Sea. The differential changes observed in the marine environment of these food webs are reflected in the assessed HTL criteria.
The reason why these systems responded differently to a reduction in nutrient load is probably because the level of eutrophication observed in the Baltic Sea is worse than in other regional seas, as shown also by McQuatters-Gollop et al. (2009), and oxygen is one of the main ecosystem drivers (Ehrnsten et al., 2019) thus, the system is more prone to improve, even if little as in this case, than the others. Theoretical studies suggest that an increase in species diversity might occur when productivity shifts from high (eutrophic) to intermediate levels and predator and/or fishing pressure stays low to medium (Kondoh, 2001;Worm et al., 2002). In the case of the Baltic Proper, the increase was due to a reduction of nutrient input and partially to a lower predator pressure (seals) on cod, which increased due to improved oxygen conditions. By contrast, species diversity might decrease if the reduction of productivity is combined with high level of predation (predators and/or fishing), as in the case of the Black Sea and the Adriatic Sea.
The size and location of the studied ecosystems are also important. Coastal and shelf areas are more sensitive to environmental changes than larger sub-regional or regional ecosystems, since the former are at the interface between land and sea, and are subjected to a variety of anthropogenic pressures (e.g., eutrophication, fishing pressure, pollutants; Halpern et al., 2019;Duarte et al., 2020). The North-East Adriatic Sea, a small shelf and coastal ecosystem where bottom up processes are fundamental in the structure and function of the food web (Celić et al., 2018), is one of the most impacted marine system in the overall assessment (Table 1 and Figures 2, 3). The apparent susceptibility to the variation of nutrient input of this small shelf/coastal ecosystem might be also a consequence Whole Western: -0.005% Shelf: -0.011% Open waters: 0.003% Whole Western: +0.29% Shelf: +0.082% Open waters: +0.35% Whole Western: +0.11% Shelf: +0.077% Open waters: +0.12% of the difficulties in compromising the spatial resolution of the LTL model with the small area represented. Adequate spatial resolution for both LTL and HTL would improve the representation of ecosystem dynamics for small shelves and coastal areas (Solidoro et al., 2010).
The relative contribution of river input in the total provision of nutrients, and hence primary production in the marine environment, may also control to a large extent the intensity of impacts affecting the pelagic fish community. For instance, the larger decline in small pelagic fishes observed in the North Sea compared to the Celtic Sea might be the result of the North Sea having proportionally larger riverine discharges and greater levels of mixing (Heath and Beare, 2008;Holt et al., 2012) stimulating pelagic production and the pelagic pathway in the ecosystem (Heath, 2005).
However, Pérez-Ruzafa et al. (2019 showed that in coastal lagoons and coastal areas the pelagic productivity might not reflect changes in nutrient input at sea. This is because the system could channel the production and main fluxes toward the benthic system (Agnetta et al., 2019 andCresson et al., 2020), retaining excessive production in the sediment or exporting it outside the system (e.g., through species migration). These mechanisms might impede observing clear changes in the ecosystem (mainly in the pelagic indicators) despite fluctuations in nutrient input. Furthermore, although the surface of the shelf in each model domain was considered to account for the observed changes, the different proportion of coastal environments and coastal areas included in each model might have inevitable effects on coastal biogeochemical processes ( Table 1) cascading up in the food webs (Figures 2, 3) and partially explaining the observed differential changes.
In addition, the stronger reduction in total nitrogen and phosphorus in the land-sea interface rather than in the dissolved nutrients at sea, [as shown in Table 1 and in more details in Friedland et al. (2021)], indicate that marine biogeochemical systems depend not only on total nitrogen and phosphorus from riverine inputs, but also on the internal nutrient dynamics, e.g., mixing and stratification of DIN and DIP from deeper layers. Therefore, it is expected that reducing one input source (river) does not relate to a 1:1 reduction of marine (dissolved) nutrients at sea. This could explain the little reduction of primary production and the little impact of these changes on the HTL criteria/indicators.
The TL indices assessed in this study, the trophic level of the community (mTLco) and the trophic level of landed catches (TLc), showed little variation (depending on the scale considered) when applying these nutrient management measures. Shannon et al. (2014) already reported no clear pattern or response of TL-based indicators to changes in Chl-a, suggesting also that these indicators might not necessarily reflect changes at the bottom of the food web. According to Heymans et al. (2014), the trophic level of the catch is highly influenced by ecosystem traits such as latitude, basins and depth, which should be taken into account when evaluating these indicators as proxy of food web dynamics. Coll et al. (2016) highlighted the usefulness of the trophic level of the community indicator in assessing the impact of fishing on the whole ecosystem. However, our study shows that this index is not sensitive enough to capture changes in the food web from environmental drivers, such as nutrients and associated ecosystem productivity, as is also observed by Fu et al. (2019). This result might be related to the fact that bottomup modifications, such as those induced by changes in nutrient inputs, are not as evident in the food chain as top-down forces, which generally resonate beyond the planktivore level, causing trophic cascades (Borer et al., 2006).
In addition, the weak responses of some of the criteria/indicators to changes in nutrients might be due to the short time series of forcing data utilized in this assessment. It is well known that ecosystems that accumulated nutrients during eutrophication require long recovery times to see large effects of load changes on ecosystem dynamics (Moloney et al., 2010;Murray et al., 2019). Thus, limiting the simulations to 8 years could have impeded a clear cause-effect relationship. This limitation was also highlighted in the results of LTL modules (Friedland et al., 2021;Grizzetti et al., 2021) which showed that 8 years of simulations were not enough to reach a new equilibrium. The internal nutrient dynamics and long residence times hampered the effect of the nutrient input reductions in the assessed ecosystems. Further efforts should be made to assess the potential impact of longer time series of changed nutrient loads in the various EU ecosystems.
It is also important to acknowledge that coastal processes are not well represented by the spatial models available here (which include both HTL and LTL models), e.g., the responses of species fished near the coast such as European seabass, Dicentrarchus labrax, or Black seabream, Spondyliosoma cantharus, are not well captured in these models. Similarly, potential improvements due to reductions in eutrophication in the pelagic habitat of coastal spawning species such as pikeperch, Stizostedion lucioperca, are not captured by the whole-sea models. The impact of bottom up forces on the marine food web might also be difficult to predict if inedible autotroph species are not properly modeled. These, in fact, might act as trophic dead-end species (Akoglu et al., 2014) or nutrient "sponges" (Murdoch et al., 1998), reducing the carrying capacity of edible species and, consequently, their trophic role (Murdoch et al., 1998;Borer et al., 2006;Akoglu et al., 2014). Finally, most of the models used in this exercise did not account for changes in dissolved oxygen (DO) impacting benthic and pelagic organisms/functional groups (e.g., impact on metabolism and recruitment success), which is another important component of eutrophication. This requires further research in the future.

Mean Change in MSFD Criteria and TL Indicators at Spatial Scale
This study assessed the impact of changes in nutrient concentrations on the spatially explicit ecosystem model of the Western Mediterranean sub-region. Our analysis confirmed that coastal and shelf ecosystems will be the most impacted when nutrients are reduced. Two out of three models suggested a slight decrease, more or less pronounced depending on the model/area considered, along the coasts/shelves of the Western Mediterranean Sea. These results are in line with previous studies that showed how freshwater pollution control measures will not impact the NW Mediterranean marine ecosystem at large, given the relatively smaller importance of river-borne nutrients for the marine productivity in the area (Macias et al., 2018a).
When inspecting change in the spawning stock biomass of European pilchard (D3C2) or the biomass of small pelagic fishes (D4C2), the two available spatial EwE models confirmed a slight decline along the continental shelf for both criteria and a marginal increase in open waters for small pelagic fishes. Yet, spatial differences were detected between the two models which could be related to the model structure and/or to the drivers used to spatially distribute the marine species and condition growth and consumption (such as changes in sea temperature). Ecosystem productivity, together with other environmental drivers such as temperature and salinity, are important factors affecting the distribution of small pelagic fishes (Bonanno et al., 2014;Quattrocchi et al., 2016;Quattrocchi and Maynou, 2017) so, including or excluding these factors may produce different outputs. The Osmose model, on the other hand, showed a different picture from the two EwE models, highlighting no clear spatial patterns among the criterion and TL indicators assessed. There are three main explanations for these differences. First, Osmose and EwE models differed in terms of model structure, process formulations (e.g., trophic assumptions) and parameterization (e.g., spatial distribution of species), which can lead to marked differences in the projections (Travers et al., 2010;Smith et al., 2015). Second, in the Osmose model, predation is a size-based and opportunistic process, which tends to buffer and dilute, via the numerous trophic links, the direct effects of changes in primary and secondary productions on the biomass of predators such as small pelagic fishes (Travers et al., 2010;Smith et al., 2015). As suggested by Travers et al. (2010), this predation formulation makes the system modeled by Osmose more variable and resilient to perturbations (nutrients reduction measures here). Third, the Osmose model was forced by the biomass of six phyto-and zoo-planktonic groups. Each group responded differently in time and space to changes in nutrient concentrations, which may have led to more heterogeneous spatio-temporal changes in the criteria considered. The spatial analysis confirmed the results obtained from the mean change for the TL indices and the species diversity index (see section "Mean Change in MSFD Criteria and TL Indicators"), with weak signals in all the assessed spatial compartments.
Overall our study confirmed that spatial modeling is still a challenging component of HTL ecosystem approaches as previously shown (Piroddi et al., 2015). Yet, because it is a fundamental aspect for guiding policy decisions (Liquete et al., 2016), its importance has increased considerably in recent years. Despite their ascertained importance in supporting policy and policy makers, spatial ecosystem models covering the entire food web from nutrients, phytoplankton to top predators, particularly at large scale, have hardly ever been utilized in the policy making process because of their high level of uncertainty (Fulton, 2010). Several studies have suggested possible ways (e.g., Bayesian network; Paradinas et al., 2015;Coll et al., 2019) to reduce specific aspects of such uncertainty. One of these approaches is the use of model ensembles to increase the reliability of model predictions, estimate the associated uncertainty (e.g., lack of spatial distribution data for many fish species and many marine ecosystems) and better assist our policies (Boyce et al., 2020). Our study, while going in this direction, highlights current challenges to fully implement such an ensemble framework in the context of the European Regional seas, because of the limited number of available spatial HTL models. The same area is rarely covered by more than one model, and existing models are not structured following a standardized approach for inter-comparison of results. Future work should address these shortcomings, benefiting from recent novel modeling developments (Spence et al., 2018), protocols for comparative modeling ensembles (Tittensor et al., 2018), and future opportunities for modeling development under the new Decade of Ocean Science for Sustainable Development (2021-2030) (Heymans et al., 2020).

Using Ecosystem Models to Support EU Marine Policies
This study investigated the effects of applied inland management measures that aim to reduce nutrient pollution in the marine environment. Using a broad set of HTL marine ecosystem models covering most of the European seas, this study was able to assess the response of these marine ecosystems to land-based measures using the criteria defined by the GES Decision (European Commission, 2017). The MSFD Descriptors and criteria have been designed to directly support and facilitate Member States in the assessment and reporting of GES. When trying to align the available modeling tools to such regulation, some difficulties were encountered. For example, not all models included species diversity explicitly (such as single species of dolphins or seabirds) and for this reason, D1 criteria could not be assessed in this study. However, the species diversity indices, which in the legislation are associated to the integrative Food Web Descriptor (D4), could also be used to evaluate the state of D1. Other modeled indicators have the potential to be useful for assessing GES (Tam et al., 2017), if easily interpretable, capable of describing food web changes and sensitivity to pressures and should be provided to policy makers to complement the assessments of the marine environment. In addition, the scale of assessments defined by the different policies (e.g., Common Fisheries Policy [CFP], MSFD, WFD) is specific for specific indicators (e.g., for eutrophication, one of the required scales is the coastal zone defined as one nautical mile from the coast; for fisheries, the scale is the FAO divisions/subdivisions) and should be assessed, if possible, using the existing modeling tools but it is not a binding requirement. Such scales are, in fact, used by Member States in their monitoring programs and reporting of GES. Overall, the role and strength of marine ecosystem models in support policies should be to: (1) highlight issues that cut across criteria/indicators and descriptors; (2) evaluate the responses of indicators to single or multiple stressors; (3) provide decision makers with a tool that can assess compartments of the marine ecosystems not directly measurable by current monitoring programs; and (4) evaluate the response of ecosystems to potential management measures.
Another important aspect where ecosystem models can support policy is in the setting of meaningful threshold/target values. As highlighted by Piroddi et al. (2015), there is a lack of tested and validated threshold/target values compliant with specific legal requirements. This is mainly due to disagreements among stakeholders in the settings of targets, which currently limits the full assessment of specific measures on marine ecosystems using the modeling tools.

CONCLUSION AND PERSPECTIVES
This study suggests that improved nutrient management, in line within European directives to preserve and/or recover the status of coastal and marine water status, will have little impact on the assessed HTL marine ecosystems. Riverine nutrient discharge, though, is just one of many stressors impacting our seas and further modeling studies should investigate the impact from synergistic and antagonistic stressors. Pressures like climate change, overfishing, chemical pollutants and plastics are expanding rapidly throughout the world (Halpern et al., 2019;Duarte et al., 2020). Thus, a holistic approach to marine management, such as that provided by the MSFD, is essential. Mechanistic modeling tools, like the ones used in this study, have the capabilities of assessing and providing useful information on the impacts of cumulative anthropogenic pressures on every trophic level of a marine ecosystem, and evaluating short/long term forecast of selected policy measures. This pan-European ensemble modeling study represents the first exercise for large scale policy evaluation within the EU framework. This work gives indications for further improvement of modeling tools, such as standardization of model evaluation, their sensibility and creation of model ensembles that will provide reasonable confidence intervals for policy making decisions.
In 2019, the United Nations (UN) declared the Decade on Ecosystem Restoration with the purpose of "recognizing the need to massively accelerate global restoration of degraded ecosystems, to fight the climate heating crisis, enhance food security, provide clean water and protect biodiversity on the planet" and in 2021, the UN Decade of Ocean Science for Sustainable Development will begin, aiming to "developing scientific knowledge, building infrastructure and fostering relationships for a sustainable and healthy ocean." It is now time to utilize these modeling tools to better guide and support decisions making by managers and policy makers.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
CP developed the protocol, led the data analysis, plotting, and writing of the manuscript. BG provided riverine nutrient inputs. DM, EG-G, RF, SM, and OP provided nutrients and phytoplankton data at sea. All authors contributed to the writing of the article, provided model results and, discussions on the analysis.