The Celtic Sea Through Time and Space: Ecosystem Modeling to Unravel Fishing and Climate Change Impacts on Food-Web Structure and Dynamics

Both trophic structure and biomass flow within marine food webs are influenced by the abiotic environment and anthropogenic stressors such as fishing. The abiotic environment has a large effect on species spatial distribution patterns and productivity and, consequently, spatial co-occurrence between predators and prey, while fishing alters species abundances and food-web structure. In order to disentangle the impacts of the abiotic environment and fishing in the Celtic Sea ecosystem, we developed a spatio-temporal trophic model, specifically an Ecopath with Ecosim with Ecospace model, for the period 1985–2016. In this model, particular attention was paid to the parameterization of the responses of all trophic levels to abiotic environmental changes. Satellite remote sensing data were employed to determine the spatial distribution and annual fluctuations of primary production (PP). Spatial and temporal changes in the habitat favorable for zooplankton were predicted with a novel ecological-niche approach using daily detection of productivity fronts from satellite ocean color. Finally, functional responses characterizing the effect of several abiotic environmental variables (including, among others, temperature, salinity and dissolved oxygen concentration, both at the surface and at the bottom) on fish species groups’ habitat suitability were produced from the predictions of statistical habitat models fitted to presence-absence data collected by multiple fisheries-independent surveys. The dynamic component of our model (Ecosim) was driven by time-series of fishing effort, PP, zooplankton habitat suitability and abiotic environmental variables, and was fitted to abundance and fisheries catch data. The spatial component of our model (Ecospace) was constructed, for specific years of the period 1985–2016 with contrasted abiotic environmental conditions, to predict the variable distribution of the biomass of all functional groups. We found that fishing was the main driver of observed ecosystem changes in the Celtic Sea over the period 1985–2016. However, the integration of the environmental variability into the model and the subsequent improvement of the fit of the dynamic Ecosim component highlighted (i) the control of the overall pelagic production by PP and (ii) the influence of temperature on the productivity of several trophic levels in the Celtic Sea, especially on trophic groups with warm and cold water affinities. In addition, Ecospace predictions indicated that the spatial distributions of commercial fish species may have substantially changed over the studied period. These spatial changes mainly appeared to be driven by temperature and may, therefore, largely impact future fisheries given the continuity of climatic changes.

Both trophic structure and biomass flow within marine food webs are influenced by the abiotic environment and anthropogenic stressors such as fishing. The abiotic environment has a large effect on species spatial distribution patterns and productivity and, consequently, spatial co-occurrence between predators and prey, while fishing alters species abundances and food-web structure. In order to disentangle the impacts of the abiotic environment and fishing in the Celtic Sea ecosystem, we developed a spatio-temporal trophic model, specifically an Ecopath with Ecosim with Ecospace model, for the period 1985-2016. In this model, particular attention was paid to the parameterization of the responses of all trophic levels to abiotic environmental changes. Satellite remote sensing data were employed to determine the spatial distribution and annual fluctuations of primary production (PP). Spatial and temporal changes in the habitat favorable for zooplankton were predicted with a novel ecological-niche approach using daily detection of productivity fronts from satellite ocean color. Finally, functional responses characterizing the effect of several abiotic environmental variables (including, among others, temperature, salinity and dissolved oxygen concentration, both at the surface and at the bottom) on fish species groups' habitat suitability were produced from the predictions of statistical habitat models fitted to presence-absence data collected by multiple fisheries-independent surveys. The dynamic component of our model (Ecosim) was driven by time-series of fishing effort, PP, zooplankton habitat suitability and abiotic environmental variables, and was fitted to abundance and fisheries catch data. The spatial component of our model (Ecospace) was constructed, for specific years of the period 1985-2016 with contrasted abiotic environmental conditions, to predict the variable distribution of the biomass of all functional groups. We found that fishing was the main driver of observed ecosystem changes in the Celtic Sea over the period . However, the integration of the environmental variability into the model and the subsequent improvement of the fit of the dynamic Ecosim component highlighted (i) the control of the overall pelagic production by PP and (ii) the influence of temperature on the productivity of several trophic levels in the Celtic Sea, especially on trophic groups with warm and cold water affinities. In addition, Ecospace predictions indicated that the spatial distributions of commercial fish species may have substantially changed over the studied period. These spatial changes mainly appeared to be driven by temperature and may, therefore, largely impact future fisheries given the continuity of climatic changes.

INTRODUCTION
Food webs govern the structure and functioning of marine ecosystems (Libralato et al., 2014) and are, therefore, typically modeled using trophodynamic models. However, marine ecosystems are far from being spatially homogeneous and the productivity of species at the base of the trophic network and prey-predator interactions can vary largely over space (Polis et al., 1997). This spatial variability determines the characteristics of food webs, including their stability and persistence (Huxel and McCann, 1998;Melian and Bascompte, 2002;Gravel et al., 2011).
The main driver of the spatial heterogeneity in preypredator interactions in marine ecosystems is the abiotic environment (hereafter simply referred to as the "environment"). Environmental variables such as salinity, temperature and chlorophyll-a govern species distributions by defining their "environmental niche" (Hutchinson, 1957;Chase and Leibold, 2003). Hence, by regulating the degree of overlap between species niches, the environment largely influences species interactions and shapes communities and food-web structure along gradients (Pellissier et al., 2013;Kortsch et al., 2019). In addition, topography and the structural habitat play an essential role in the spatial variability of marine food webs. In particular, water depth greatly influences exchanges between the pelagic and benthic compartments of ecosystems through food accessibility (Giraldo et al., 2017;Agnetta et al., 2019). Seabed features (including substrate composition) and water column characteristics also mediate species interactions within marine ecosystems by influencing the accessibility of potential prey items to predators (Gardner, 1981), by providing refugia to the prey (e.g., Menge and Lubchenco, 1981) or by facilitating predation through aggregation (Le Fèvre, 1987).
The complexity of marine food webs is enhanced by the responses of predators to their prey, and by competition (Hunsicker et al., 2011). The feeding rate of any predator on a given prey depends on the abundance of that prey (Holling, 1959) and the relative preference of the predator for other cooccurring prey items (Murdoch, 1969). The spatial distribution of a predator can be substantially affected by the spatial distribution of its potential prey (Logerwell and Hargreaves, 1996;Torres et al., 2008), and vice versa (Barnett and Semmens, 2012), as well as by the spatial distribution of the other predators that compete for those prey (Hansson, 1984;Ward et al., 2001). Therefore, identifying the relative role of the drivers of the structure and functioning of marine ecosystems and their changes requires spatially-explicit ecosystem models where the spatial patterns of species biomass depend on both a specific environment and trophic interactions.
Spatially-explicit ecosystem models such as Atlantis (Fulton et al., 2004;Audzijonyte et al., 2019), OSMOSE (Shin and Cury, 2001) and Ecopath with Ecosim and Ecospace (Pauly et al., 2000;Christensen et al., 2014) simulate food-web functioning while representing the spatial distribution patterns of all ecosystem components. These spatially-explicit ecosystem models are powerful tools to disentangle the intricate effects of trophic relationships, the environment and anthropogenic stressors, e.g., fishing, on species distribution patterns and biomass flows within marine ecosystems. In the context of climate change, spatially-explicit food-web models that integrate predictions from statistical habitat or niche modeling approaches are particularly useful to predict how ecosystems may respond to future modifications in environmental conditions (Pellissier et al., 2013;Moullec et al., 2019a). However, such spatial foodweb models are currently relatively scarce and these models have often difficulty to comprehensively represent heterogeneous spatial processes.
The Celtic Sea represents a large portion of the Western European continental shelf, and extends from the Western English Channel to the Celtic break delimiting Porcupine Sea Bight and the South-West of Ireland (Figure 1). This ecosystem has a complex spatial organization since it encompasses a broad range of depths and a diversity of substrate types, resulting in a variety of benthic habitats that host a wide diversity of species assemblages (Le Danois, 1948;Ellis et al., 2013Ellis et al., , 2002Martinez et al., 2013). This diversity of species assemblages is conducive to complex trophic interactions (Trenkel et al., 2005). In addition, the Celtic Sea is located at the interface of the Lusitanian and Boreal provinces, therefore defining the northern or southern limit of the distribution range of many species (Dinter, 2001). The presence of pelagic mesoscale structures such as productivity fronts (Pingree, 1980;Pingree et al., 1982) is another source of spatial heterogeneity in the Celtic Sea ecosystem. Ultimately, spatial heterogeneity in Celtic Sea fish communities is reflected by the highly mixed fisheries that operate in the Celtic Sea region (Sharples et al., 2013;Mateo et al., 2017;Dolder et al., 2018;Moore et al., 2019). Fishing has been identified as the primary driver of changes in the Celtic Sea ecosystem since 1950 and this pressure has culminated in the late 1990s  FIGURE 1 | General location of the Celtic Sea (top left corner) and delimitation of the study area regarding physical and management boundaries. The Celtic Sea is located at the border between Lusitanian and Boreal provinces. The Celtic Sea EwE and Ecospace model (red dashed contour) is developed for the shelf area shallower than 200 m depth matching ICES divisions (polygons) 7.e (Western English Channel), 7.f (Bristol Channel), 7.g (Celtic Sea North), 7.h (Celtic South) and 7.j.2 (Southwest of Ireland-East belonging to the NE Atlantic Fisheries Commission regulatory area). Hernvann and Gascuel, 2020;Mérillet et al., 2020). Nonetheless, long-term variability in hydro-climatic conditions has also largely affected the Celtic Sea ecosystem (Beaugrand et al., 2000;Hernvann and Gascuel, 2020) and the ongoing climate change may greatly impact the basal ecosystem production (McGinty et al., 2011), fish populations (Brunel and Boucher, 2007), and entire communities (ter Hofstede et al., 2010;Simpson et al., 2011) of the Celtic Sea.
In this study, we analyzed the interdependent impacts of fishing and climate change on the structure and functioning of marine ecosystems using the Celtic Sea as a regional case study. We developed a spatially-explicit trophic model for the Celtic Sea to represent the effects of the environment on the distribution and productivity of functional groups, from phytoplankton to top predators. Multiple data types were considered for this aim. Spatio-temporal dynamics of plankton groups were represented in the spatially-explicit trophic model using the predictions of remote sensing based models. The responses of higher trophic level functional groups to changes in the environment were represented using the predictions of statistical habitat models, which were fitted to data from large scale survey and environmental databases. This overall information was used to calibrate the trophic model and determine how changes in primary production, zooplankton habitat and temperature may have interacted with fishing to shape the Celtic Sea food-web dynamics over the period 1985-2016. Then, the spatial structure and functioning of the Celtic Sea food web were analyzed for contrasted years in terms of environmental conditions to better understand the processes of changes in the Celtic Sea ecosystem over the studied decades.

MATERIALS AND METHODS
The EwE Framework in a Nutshell and General Approach to Model the Trophic Functioning of the Celtic Sea

The EwE Modeling Framework in a Nutshell
Our approach relies on the Ecopath with Ecosim (EwE) framework (Pauly et al., 2000;Christensen and Walters, 2004). Hereafter, we only briefly summarize the EwE modeling steps since more details about the EwE framework and its main equations can be found in Christensen and Walters (2004) and are summed up in Supplementary Appendix SA.
Ecopath is the basal component of EwE, and provides a snapshot of the ecosystem of interest for a given period, i.e., generally a specific year. Ecopath represents the ecosystem of interest through functional groups (i.e., groups of species/taxa that have similar trophic and ecological niches and life-history characteristics) that exchange matter and energy through trophic relationships while the mass balance between the production of each functional group and all consumptions is assumed.
Ecosim is the dynamic module of EwE that allows for hindcast or forecast of changes in functional group biomasses. Ecosim predator-prey dynamics (under top-down or bottom-up control) are determined based on the foraging arena theory (Ahrens et al., 2012) whereby vulnerability coefficients modulate the proportion of prey biomass accessible to predators by altering predator consumption rates. Vulnerability coefficients can be estimated by fitting Ecosim to observed biomass, abundance or catch time-series. The response of the ecosystem to stressors FIGURE 2 | Schematic of the data integration framework used to build and calibrate the Ecopath with Ecosim (EwE) with Ecospace model. Gray rectangles represent data or EwE input. Blue rectangles correspond to statistical model runs to estimate EwE inputs from data.
can be simulated in Ecosim: (1) time-series of fishing effort or fishing mortalities can, respectively, drive exploited groups dynamics through top-down processes; (2) time-series of primary production can drive phytoplankton dynamics controlling bottom-up the food-web; and (3) environmental time-series can affect the consumption rates of predators according to linear response of following specific user-defined functional responses.
Ecospace is the spatial component of EwE, which predicts the spatial distribution of functional group biomasses over a horizontal 2-D grid (Walters, 1999;Pauly et al., 2000). The most recent version of Ecospace includes an habitat foraging capacity model that determines the suitability of each model grid cell for a functional group as the relative proportion of the cell over which the group can feed (Christensen et al., 2014). Each grid cell's capacity is calculated as the product of functions that reflect the preferences of the group of interest for specific environmental parameter values.

General Approach of Integrated Trophic Modeling in the Celtic Sea
The Celtic Sea ecosystem was previously modeled by Moullec et al. (2017) using the Ecopath and Ecosim modules of EwE to represent its dynamics and fishing impact between 1980 and 2013. However, this model relied on a limited amount of fieldwork data and did not investigate the effect of any environmental driver and their spatial heterogeneity on ecosystem structure and functioning over this period. To improve the ecological realism of the Celtic Sea food-web model, we develop an integrated trophic modeling approach consisting in successively (i) updating the Ecopath model and integrating local trophic data to define trophic relationships, (ii) exploring spatio-temporal changes in the dynamic of plankton groups, (iii) investigating the environmental drivers of high trophic levels, including temperature, and studying their spatio-temporal changes, and (iv) integrating fishing and changes in temporal variation in environmental drivers of functional groups in the fit of the Ecosim model over the period 1985-2016 and using their spatial distribution to build the Ecospace model. The different steps of this integrated trophic modeling framework are displayed in Figure 2.

Updating and Improving the Pre-existing Celtic Sea Ecopath Model
We updated the pre-existing Ecopath models for 2013 and 1980 (Moullec et al., 2017) and built a new Celtic Sea Ecopath model for the year 1985. The species composition of a few functional groups was refined for a greater coherence in the diet of the species they contain (see details in Supplementary Appendix SB). A major upgrade of the Celtic Sea model was the provision of a new diet matrix for the Ecopath model, estimated using the Bayesian integrated model "EcoDiet" and its R package implementation (Hernvann et al., in press;Théro et al., 2020) 1,2 . EcoDiet allows to simultaneously estimate diet proportions of all functional groups by integrating multiple data types, i.e., stomach content analyses, biotracer analyses and literature knowledge.
Here, EcoDiet was applied to trophic analyses conducted on fish collected in the Celtic Sea during a fisheries-independent survey as part of the EATME project (Issac et al., 2017;Rault et al., 2017;Day et al., 2019), and information from an updated bibliographic search. The new model comprised a total of 50 single and multispecies functional groups, including exploited fishes and invertebrates, non-commercial benthic invertebrates (five groups), zoo-and phyto-plankton (four and two groups, respectively), and mammals and seabirds (four groups). Some of these functional groups were further split into younger and older life-stage groups (cod Gadus morhua, hake Merluccius merluccius, and anglerfish Lophius sp.) based on onthogenetic changes in diet identified by stomach content analyses (Issac et al., 2017;Day et al., 2019). Details on the Celtic Sea Ecopath model, its parameters, and the species it represents can be found in Supplementary Appendix SB. The validity of the 1985 Ecopath model was analyzed using the standard PREBAL diagnostics (Link, 2010; Supplementary Appendix SB).

Representing the Dynamics of Plankton Groups
Spatial-temporal changes in the production of phytoplankton groups (primary production, PP) were taken from the Vertically Generalized Production Model (VGPM, Behrenfeld and Falkowski, 1997) for the period 1998-2016, and from a regional application of hydrodynamic-biogeochemical NEMO-ERSEM to the European shelf (Baretta et al., 1995;Madec, 2008;Edwards et al., 2012; UK Met Office 3 ) for the period 1985-1997. The VGPM model uses remote-sensed chlorophyll-a data collected by SeaWiFS and MODIS-Aqua to quantify PP over the entire water column (Morel and Berthon, 1989). PP is one of the direct outputs of the NEMO-ERSEM model and was extracted over the entire water column to align with the VGPM outputs. Data products were obtained from the Ocean Productivity website 4 and the E.U. Copernicus Marine Environment Monitoring Service 5 (eu), respectively. Time-series of PP in the Celtic Sea were used to force the productivity of phytoplankton groups in Ecosim (Christensen et al., 2009;Mackinson et al., 2009) while annual maps of primary productivity were used in Ecospace (Steenbeek et al., 2013). PP maps averaged by 5-year periods are provided in Supplementary Appendix SC.
The spatio-temporal changes of the suitable feeding habitat of zooplankton (small and large mesozooplankton, macrozooplankton) were estimated for the 1998-2016 period using a niche model of mesozooplankton in the North Atlantic (Druon et al., 2019). This zooplankton suitable habitat (ZSH) model relies on the identification of productivity fronts at mesoscale. Because these fronts persist long enough in the pelagic realm (at least 3 to 4 weeks), they allow the production of elevated biomass of mesozooplankton which, in turn, can sustain higher trophic level organisms. The mesozooplankton habitat model is fitted to data of the Continuous Plankton Recorder (CPR) and relies on remotely-sensed horizontal gradients of chlorophyll-a (SeaWiFS and MODIS-Aqua sensors). This habitat model only results of the daily observation of productivity fronts by satellite sensors, therefore the ZSH model is an algorithm that translates chlorophyll-a distribution directly into a daily habitat value. Seasonal maps of ZSH in the Celtic Sea are provided in Supplementary Appendix SC.

Linking the Dynamics of Higher Trophic Level Organisms to Environment
To characterize the response of higher trophic level organisms to environment, we determined functional relationships for 38 Ecopath groups (over 45 high trophic level groups, excluding mammals and birds, and benthic meiofauna) using the statistical habitat modeling approach developed by Grüss et al. (2020Grüss et al. ( , 2018. This approach involved the following steps: (1) the construction of fisheries-independent survey and environmental databases; (2) the development of generalized additive models (GAMs); and (3) the use of the GAMs to estimate functional relationships.

Construction of Survey and Environmental Databases
For each functional group, presence/absence (encounter/nonencounter) data were collected by 13 International Bottom Trawl Surveys (IBTS; Cronne, 2016) conducted between 1987 and 2017 in European waters (Supplementary Appendix SD; ICES database) 6 . We allocated presence/absence data, collected by the surveys (only 3 surveys for 4 low TL benthic invertebrates groups; Supplementary Appendix SD), to a 0.125 • spatial grid to address the issue of heterogeneity in sampling schemes (Grüss et al., 2020). A set of environmental variables were compiled for all of these grid cells. These environmental variables included time-invariant topographic and habitat variables: bathymetry, distance to the coast, and percentages of seabed substrate types. The environmental variables also included hydro-climatic time-varying variables: surface (first 10 m) and bottom sea temperature, salinity and dissolved oxygen concentration, and surface phytoplankton concentration. Spatial-temporal estimates of these variables for each grid cell at a monthly resolution were obtained from the General Bathymetric Chart of the Oceans (GEBCO) 7 , the European Marine Observation and Data Network (EMODnet) 8 , and the Copernicus platform. More details are provided in Supplementary Appendix SD.

GAM Fitting
The literature was screened to determine which of the variables included in the environmental database may influence the spatial distribution patterns of each functional group. This list was reduced through a collinearity analysis to avoid the selection of variables that are strongly correlated with one another together or with eastings and northings (longitude and latitude expressed in UTM coordinates) (Grüss et al., 2018). Then, for each functional group, we fitted binomial GAMs of the following form to the presence/absence data from the survey database ("mgcv" package in the R environment Wood, 2017): where η is the probability of presence of the functional group; g is the logit link function between η and each predictor; te(X,Y) is a tensor product smooth fitted to eastings and northings that accounts for spatial autocorrelation at broad scale; x 1 , x 2 ,. . ., x n are the predictors from the environmental database; s is a thin plate regression spline fitted to a given predictor and limited to four degrees of freedom for interpretability; and Survey and Year are treated as fixed-effect factors (Farmer and Karnauskas, 2013;Grüss et al., 2018). Only predictors with a significant effect (statistical criteria detailed in Supplementary Appendix SE) were kept in the GAMs, and the models were evaluated using the interactive cross-validation procedure described in Grüss et al. (2018). More details can be found in Supplementary Appendix SE.

Production of Preference Functions
Fitted GAMs were utilized to produce preference functions (also called functional reponses below) expressing for the preferences of each functional group for specific environmental parameter values. For each environmental variable i integrated in a GAM, predictions were made with the GAM over a vector of values ranging between the minimum and maximum values of environmental variable i in the Celtic Sea, while: (1) setting eastings X and northings Y at their values at the barycenter of the Celtic Sea; (2) keeping the other variables at their mean value from the GAM modeled dataset; and (3) setting the Survey and Year factors at their mode from the GAM modeled dataset. GAM predictions were then rescaled to provide preferences functions ranging between 0 and 1, which account for the relative effects of each environmental variable on the probability of presence of the functional group of interest. The procedure for rescaling follows Grüss et al., 2018 and is described in Supplementary Appendix SE.

Ecosim
The Ecosim model was first fitted to observed biomasses from available ICES stock assessments, abundance indices from scientific surveys for fish (EVHOE; Duhamel et al., 2018; UK-WCGFS, CEFAS), the CPR survey 9 for plankton (Johns, 2019), 9 https://www.cprsurvey.org/ and catch time-series (Statlant) 10 . Time-series of fishing mortality (for assessed stocks) or catch (for non-assessed stocks) were used as forcing functions. A summary of the time-series used is provided in Supplementary Appendix SF. Ecosim calibration was conducted using the stepwise fitting procedure implemented into EwE by Scott et al. (2016) and following the protocol found in Bentley et al. (2019) to determine the most suitable number of predator-prey dynamics parameters (vulnerabilities Vs) to estimate. To assess the relative impacts of fishing and environment on the ecosystem, we adopted an approach similar to the one developed by Mackinson (2014). We applied the fitting procedure, integrating fishing forcing functions only, then integrating the environmental drivers of plankton and fish dynamics, first separately and then simultaneously. PP and zooplankton suitable habitat (ZSH) time-series drove the production of phytoplankton groups (see Christensen et al., 2009) and the consumption of zooplankton (small and large mesozooplankton, macrozooplankton), respectively. More specifically, ZSH time-series were used as multiplier of the arena area for zooplankton-phytoplankton interactions, of zooplankton vulnerability for fish-zooplankton interactions, and of both parameters for zooplankton-zooplankton interactions (see Supplementary Appendix SF; Dahood et al., 2019). Temperature variations and corresponding preference functions to the model drove the consumption of fish groups (foraging environmental response; see Bentley et al., 2017;Serpetti et al., 2017;Corrales et al., 2018). In order to avoid any potential confounding effects with temperature-induced changes, that we precisely wanted to quantify in this study, we chose to not use the available information on recruitment anomalies (e.g., from stock assessments) as forcing function of the dynamics of functional groups represented by multiple life-stages. To quantify the relative contribution of each driver, the various fits of the models were compared using the sum of squares between observed and predicted time-series and the bias-corrected Akaike information criterion (AICc) of the model (Akaike, 1974). We addressed the uncertainty in Ecopath input parameters in Ecosim outputs through the Monte Carlo uncertainty routine based on confidence intervals determined from Pedigrees Steenbeek et al., 2018). More details on the fitting procedure can be found in Supplementary Appendix SF.

Ecospace
The fitted Ecosim scenario was run over the period 1985-2016, and the "Ecopath model from Ecosim" plug-in  was employed to generate 6 Ecopath models corresponding to mean conditions over successive 5-year periods (from 1985-1989 to 2005-2010, and 2011-2016 for the last period). Ecospace models were built for each of the 5-year periods to facilitate the comparison of spatial distributions through time and run over a spin-up period of 5 years. The Ecospace models integrated the average maps of PP to force its distribution, those of ZSH to feed the distribution of habitat foraging capacity for zooplankton, and those of environmental variables employed in the GAMs and the associated preference functions of higher TL groups to spatially distribute fish habitat foraging capacities Grüss et al., 2018Grüss et al., , 2016. Trophic and diversity indicators ( Table 1) were projected based on functional groups' distributions, partly using the ECOIND plug-in (Coll and Steenbeek, 2017). This approach decoupling temporal and spatial changes will be further discussed in the Discussion section. Parameters of the Ecospace model are given in Supplementary Appendix SG.

Improvement of Ecosim Model Fit and Relative Importance of Ecosystem Drivers
The statistically best model was obtained when trophic interactions, fishing, and environment were included together in the model fitting (SS = 340 and AICc = −3365, Table 2), hence improving the model fit by 86%. The proportions of each driver in the Sum of Square SS reduction suggest that fishing was the main driver of ecosystem changes, followed by trophic interactions and, with much lower but significant contributions, the environment ( Table 2 and Figure 3). Temperature effects on higher trophic levels and overall changes in low trophic levels productivity had similar weight in the fit. Integrated independently from plankton data, temperature forcing especially increased the hindcast quality for species with thermal preferenda at outside the temperature ranged covered in the Celtic Sea (warmer or colder, Figure 3B). Conversely, PP and ZHS integration essentially led to a better fit for functional groups belonging to the pelagic pathway.

Low Trophic Levels
With regard to the basal productivity in the Celtic Sea ecosystem, no clear trend was apparent in the Ecosim model over the period 1985-2016. Despite this relative stability, some differences can be observed in CPR time-series within the phytoplankton community. Most significant changes in observed phytoplankton data occurred during the late 1980s when the abundance of large phytoplankton (diatoms and dinoflagellates) dropped whereas the one of small phytoplankton increased (Figure 4). However, the trends for both large phytoplankton and small phytoplankton were not correctly captured by Ecosim and remained overall constant over the entire period 1985-2016 in the model. Similarly to PP, the ZSH remained stable over the study period, yet, looking at the data, substantial changes were observed within the zooplankton community, mostly during the late 1980s. When PP and ZSH were integrated, Ecosim reproduced the overall observed trend in small mesozooplankton abundance, which substantially declined before 1990 and was still decreasing until the mid-2000s. By contrast, large mesozooplankton abundance tended to increase over the period 1985-2016. Again for this zooplankton group, the trend predicted by Ecosim correctly matched observations both in terms of trend and of amplitude of variations. Both trends and interannual variability in zooplankton abundance were not captured when PP and ZSH were not included in the fit (Figure 4). The stability in the predicted phytoplankton timeseries precludes us from linking changes in mesozooplankton abundance with changes in primary productivity in the Celtic Vs columns indicate the number of vulnerability parameters that have been estimated and whether they have been estimated by predator or by couple predator-prey. K is the total number of parameters estimated by the model, which is equal to the number of Vs estimated as no PP anomaly is estimated. The maximum number for K to avoid overfitting is 56 (Supplementary Appendix SE). The best model in terms of minSS, the sum of squares between observed and predicted values reached at the end of the fit, and AICc is the one including both trophic interactions, fishing, plankton-related data and temperature effect on high trophic level groups. The best model statistically is written in bold. Sea between 1985 and 2016. Although Ecosim was able to correctly predict the overall zooplankton abundance level over the study period, the predictions did not perfectly match the observed patterns, especially at the beginning of the timeseries. FIGURE 4 | temperature forcing; orange line: with integration of temperature forcing but without PP and ZSH. The gray shaded area is delimited by the 95% intervals of the Monte Carlo runs that propagate the uncertainty on Ecopath input parameters to Ecosim. The points correspond to observed time-series used during Ecosim calibration.

Higher Trophic Levels
The observed time-series of biomass, abundance and catch of most of high trophic level groups were satisfactorily reproduced by Ecosim. Overall, Ecosim predictions were better for those single-species groups for which stock assessments are carried out than for multispecific fish groups, and they were less consistent for exploited invertebrate groups.
In 1985 were already subject to high fishing mortality that far exceeded natural losses. Driven by an increase in fishing mortality, the biomass of most of these groups declined between 1985 and the mid-2000s while their catches simultaneously stabilized or decreased. In 2005, their fishing mortality decreased and Ecosim predicted a stabilization of their abundance that preceded a progressive recovery from 2010 onward (Figure 5).
The pressure exerted on species that were initially less exploited such as smaller flatfish (i.e., endobenthivorous), medium demersal fish (i.e., epibenthivorous), or cephalopods increased after 1985. Concomitantly, the natural mortality exerted on these species decreased as a result of the depletion of highly commercial groups such as large predatory fish (e.g., hake, cod, whiting, megrim). The predicted biomass trends of lower TL fish and invertebrates varied according to their exploitation patterns ( Figure 5). Dramatic changes in the composition of the pelagic community were also predicted by Ecosim. In parallel to horse mackerel depletion, the biomass of pelagic fishes that feed upon similar ressources such as suprabenthivorous fish, sprat (Sprattus sprattus) or boarfish (Capros aper) increased, while these pelagic fishes were initially almost absent from the Celtic Sea ecosystem.
The recovery of commercial species in the 2010s is noticeable, particularly the large increases in the biomasses of large predators, including hake (whose biomass increased 5-fold) and a few other species. According to the model, these species recoveries translated into an increase in predation pressure from large predators and a decline in the abundance of cephalopods, small pelagic and small demersal fish.
The rising temperatures resulting from climate forcing had various effects on the different species represented in the model, and they mainly impacted fish productivity during the mid-1990, when a substantial increase in both surface and bottom temperatures was observed. As highlighted by the comparison of model's fit on observed data according to scenarios (Figures 3B, 5), the integration of temperature forcing allowed to better capture the overall trends in the abundance of Boreal and Lusitanian species. The foraging ability of Boreal species, including cod, plaice, whiting, herring and endobenthivorous demersal fish, dropped at the same time that fishing mortality increased, and this accentuated the depletion of the above-mentioned species. Conversely, for Lusitanian species that have warmer thermal preferences, such as sardine, boarfish, hake and cephalopods, increased productivity in the Celtic Sea ecosystem in the 1990s acted antagonistically to the increasing fishing pressure and synergistically with decreasing competition pressure. The differential trends in the productivity of warm and cold water affinity species between the start and the end of the study period were responsible for differences in the recovery capacities of commercial groups in the 2010s. For instance, cod biomass remains quite stable when temperature is not driving productivity, and plaice recovers much slower from fishing pressure release. Conversely, sole depletion in the first part of the time-series is stronger as its recovery is strengthened (Figure 5). These variations triggered changes in the fish assemblages. Interannual variability in high trophic level abundance was better reproduced when incorporating PP and ZSH data, especially for pelagic species (Figure 5).

Predicted Spatial Distribution of Food-Web Components in the Celtic Sea
Habitat Projection for the Celtic Sea Evaluation tests suggested that GAMs (Supplementary Appendix SE) were good enough to produce yearly foraging habitat maps for the 38 higher trophic level functional groups (year 2016, Supplementary Appendix SE). The tests also indicated that GAM fits and predictive capacities were better for monospecific groups than for heterogeneous multispecies groups. The fit was also lower when the sampling gear employed during the surveys was not adapted to a given functional group (e.g., benthic invertebrates). The type of environmental predictors retained in final GAMs varied across functional groups, as their number, ranging from 3 for the non-commercial invertebrate groups to 9 for hake. The preference functions estimated from fitted GAMs exhibited diverse shapes. The retained predictors often included bathymetry (32 groups) and temperature (surface or bottom; 35 groups), two variables of importance in shaping the Celtic Sea environment. Depth in the Celtic Sea drove marked northeast-southwest gradients, while surface temperature led to pronounced latitudinal gradients and bottom temperature resulted in two opposed pools of warm and cold waters in the English and Bristol Channel and in southern Ireland.
Salinity functional responses quantify habitat suitability in the riverine-influenced Bristol Channel. When they were retained in final GAMs, seabed substrate variables resulted in patchy predicted habitats for demersal species. In those instances, it was possible to make a clear distinction between restricted areas of rocky or mixed substrate and large irregular strands of muddy,

Predicted Biomass Distributions
The Ecospace model predicted the biomass distribution of all components of the ecosystem, from the phytoplankton to top predators. The biomass distributions of functional groups reflected the overall habitat foraging capacities of these groups, hence their ecological niche, but were also influenced by their trophic interactions with prey, predators and competitors (Figure 6).
The distributions of those functional groups for which no GAM was developed were estimated by Ecospace only according to their feeding relationships (birds, mammals, bacteria, benthic meiofauna, large pelagic fish). A wide variety of spatial patterns were exhibited by the functional groups represented in Ecospace. Due to their high mobility in the water column and their eurythermic characteristics, horse mackerel and mackerel were FIGURE 6 | Predicted biomass distribution in 2016 for functional groups for which statistical habitat models (generalized additive models) were developed. Biomass is represented in relative units, i.e., scaled on the highest density predicted in the Celtic Sea. The color scale represents the relative abundance of the considered group, from zero (blue) to the maximum (red).
Frontiers in Marine Science | www.frontiersin.org predicted to be the most widely distributed monospecific groups in the Celtic Sea ecosystem (Figure 6). Other fish groups that were predicted to occupy a large fraction of the Celtic Sea were multispecific groups such as benthivorous elasmobranchs, epibenthivorous demersal fish, pouts (Trisopterus minutus and T. luscus) and benthic cephalopods, which are all made of several species and, therefore, display a broader, composite ecological niche. By contrast, some species were predicted to occupy only a relatively small fraction of the Celtic Sea ecosystem. This was particularly the case for benthic groups, which are highly dependent on the seabed substrate. For instance, the distribution of Norway lobster (Nephrops norvegicus) was restricted to muddy patches in the western part of the central Celtic Sea where the species can more easily bury in mud. Ecospace modeling allowed us to distinguish between groups that have similar trophic niches but different environmental niches. Thus, sole and plaice have similar diets but showed limited spatial overlap because of their, respectively, warm and cold water affinities; the same was observed for sardine (Sardina pilchardus) and herring. Species that are taxonomically close such as the gadoids cod, haddock and whiting were shown to regularly coexist in the same areas but could also have substantially different distributions: whiting inhabited shallower areas than cod and haddock; both haddock and cod had affinities for muddy areas, but haddock also frequented sandy substrates whereas cod preferred coarse sediments. The distributions of adults and juveniles for hake, anglerfish and cod were comparable yet exhibited some differences (i.e., juveniles showed a greater affinity for muddy substrates), which highlighted the relevance of representing multiple stanzas for some species in Ecospace to better simulate changes in spatial distribution and predation pressure over fish life cycle and the potential for cannibalism.

Species Assemblages
The distribution of functional groups provided a view on fish community assemblages in the Celtic Sea ecosystem. Very coastal areas and bays around Ireland, United Kingdom and France and the inner Bristol Channel were found to be inhabited by functional groups tolerant to low salinities, such as large concentrations of benthic invertebrates, small benthivorous organisms and, upper in the surface, suprabenthivorous and small epipelagic fish (i.e., herring, sardine, sprat). Higher trophic levels in these areas were mainly represented by whiting and sea bass. Also, intense primary and secondary pelagic production took place in these areas. Areas next to the shelf edge, located from the tip of Brittany to southwestern Ireland, were characterized by deeper waters with a steeper slope. Here, abundant populations of megrim, hake, anglerfish and other piscivorous demersal fish were found together with widely distributed pelagic fish species such as boarfish, blue whiting (Micromesistius poutassou), mackerel and horse mackerel. Benthic invertebrates in these areas were mainly deposit and suspension feeders. Primary production was observed limited in these areas, but the pelagic habitat was favorable to zooplankton development. In the cold central and northwestern Celtic Sea, the mix of sediment types and the moderate depth was favorable to diverse demersal species, including epibenthivorous demersal fish, hake, anglerfish and three gadoid species, cod, haddock and whiting. Mackerel and horse mackerel were the dominant pelagic species in this area. The density of benthic carnivores and necrophageous invertebrates was especially high in this area. Finally, the Western English Channel and the outer Bristol Channel differed from the other regions of the Celtic Sea by their relatively warm waters and coarse substrate and their higher plankton productivity sustained by mesoscale structures. The Western English Channel and the outer Bristol Channel were characterized by the presence of small benthivorous fish, epibenthivorous fish, pouts, flatfish (plaice and sole), commercial invertebrates (crustaceans and bivalves), cephalopods species, some gadoid species, and a mix of pelagic species including mainly horse mackerel and mackerel but also small epipelagic species.

Temporal Changes in Functional Groups' Distributions
Overall, relative habitat capacity was constant for most functional groups over the study period. Most changes in relative habitat capacity occurred between the late 1980s and the mid-1990s (Figure 7). The intensity of changes was variable across functional groups and can be related to the water temperature preferences of functional groups. The fitted GAMs indeed highlighted that temperature, which is the variable most affected by climatic changes in the area, was an ecological niche descriptor common to almost all components of the ecosystem. Thus, changes in the modeled habitats over the studied period were observed for those functional groups with most extreme thermal preferences in the Celtic Sea such as pilchard (Sardina pilchardus), sea bass, medium pelagic fish and squids (mainly or exclusively composed of Lusitanian species), as well as cod, haddock, sprat, herring, endobenthivorous fish and whiting (mainly or exclusively composed of Boreal species).
Functional groups dominated by Boreal species illustrate biomass distribution changes related to environment (Figure 8). Most changes in the distribution of species that prefer cold waters occurred early in the study period, with a slight contraction of the area occupied by these species in the eastern and southwestern Celtic Sea in response to the more intense warming of already warm bottom waters in those areas. For the remainder of the study period, the relative distribution of biomass of Boreal species remained relatively stable.

Ecosystem Indicators
Primary production in the Celtic Sea ( Figure 9A) was essentially located in coastal areas, from the southern coast of Ireland to the Western English Channel and western Brittany, and along the southern part of the Celtic Sea shelf break. Differences in primary production between the lowest and the highest productive areas never exceed a factor of four. The biomass of secondary producers (i.e., zooplankton and benthic invertebrates that not carnivorous; Figure 9B) was more spread over the shelf and did not display similarly high densities in the center of the Western English Channel. The total biomass of high trophic level groups (i.e., TL ≥ 3; Figure 9C) was less variable over space and more patchy. The biomass of primary and secondary producers in the Celtic FIGURE 7 | Index of change in the spatial distribution of functional groups' habitat in relation to their thermal preferences between the 1985-1989 and 1995-1999 periods and between the 1985-1989 and 2010-2016 periods: Left: Values by functional group and averages over functional groups according to their temperature preferenda (cold, warm and intermediary corresponding to groups with thermal preferendum lower than, higher than, or within the range of temperatures in the Celtic Sea, respectively). The index is the absolute difference between 1 and the slope of the linear regression between Ecospace model cells' values of relative habitat for the 1985-1989 period and the two other periods. Right: mean index per group of thermal affinity, "cold," "intermediate," and "warm" for species with thermal preferendum, respectively, inferior, within, and superior to the temperature range covered in the Celtic Sea over 1985-2016. calculated as the ratio between these selected periods and the initial biomass. In (B-D) the colors correspond to the ratio of absolute biomass between the predicted biomass and biomass for 1985-1989. The black lines delimit areas for which the fraction in the total Celtic Sea Boreal demersal species biomass increased ("+" area) and decrease ("-" area) compared to 1985-1989. Sea was concentrated around particular areas, mainly around the coasts and in the English Channel, and this available production is efficiently transferred to higher trophic levels and then likely spread horizontally.
The spatialized mean trophic level (MTL) and high trophic index (HTI) highlighted that higher trophic levels were mainly located in the western Celtic Sea. In particular, high values of the MTL were concentrated in the northwestern Celtic Sea, whereas the HTI is higher in deeper areas along the shelf-break (Figure 10, left panel). Areas associated with a high MTL were those areas that were occupied by many medium to large demersal functional group including epibenthivorous and piscivorous fish, FIGURE 9 | Primary production in the Celtic Sea and biomass for trophic levels 2-3.5 and greater than 3.5 in 1985. (A) Primary production; (B) total biomass for trophic levels 2-3.5; (C) total biomass for tropic levels greater than 3. elasmobranchs, gadoids, anglerfish and hake. Areas characterized by a high HTI were those where the biomasses of hake, megrim and anglerfish were high. In the eastern part of the Celtic Sea, the MTL was lower, due the large abundance of small demersal fish species feeding on benthic invertebrates (e.g., sole, plaice, small benthivorous fish) in this area, particularly medium and epipelagic species, which dominate in the Western English Channel. In the eastern part of the Celtic Sea, top-predators were mainly represented by large pelagic sharks and sea bass.
The Kempton's index highlighted spatial patterns in functional group diversity. Spatial patterns in functional group diversity were independent from that of biomass levels. For instance, mean values of Kempton's index were found in areas with either low total biomass (shelf-break) or high total biomass (fronts and muddy areas) densities. Functional group diversity was found to be particularly high at the border between the main Celtic Sea regions identified previously. The SOI was similar but not identical to the Kempton's diversity index. It was highest in the northwestern part of the central Celtic Sea, which was inhabited by generalist species and species that can feed on both benthic invertebrates and fish, including gadoids and anglerfish, but not hake which primarily feeds on medium-sized pelagic fish.
All the indicators displayed measurable variations over time and space (Figure 10, middle panels). The trends in the HTI and MTL were in opposed directions, which means that the ratio of low trophic level biomass to total biomass increased overall but that, in parallel, the proportion of species with a TL greater than 4 increased. These results were mainly driven by the dramatic changes in horse mackerel and hake abundances that occurred over the study period, while, after a deep depletion, the abundance of most mesopredator in the 2000s returned to levels similar to the 1980s. For the HTI and MTL, the direction of the change was the same among Ecospace cells, partly because the species driving observed changes were widely distributed in the Celtic Sea region. Trophic diversity declined overall over the whole Celtic Sea, especially in coastal areas. While functional groups diversity also declined in the southern Celtic Sea and the English Channel, it increased in the northwestern part of the central Celtic Sea. The standard deviations of the four indices (Figure 10, right panels) highlighted a measurable modification of the whole structure of the Celtic Sea ecosystem over the period 1985-2016, yet some specific areas were more affected than others. In particular, the food-web structure of the central Celtic Sea displayed substantial changes between 1985 and 2016, which were reflected by both an increase in MTL and a decrease in HTI and SOI. The standard deviation of spatial indices relative to their mean value also suggests a greater amplitude of changes occurred in the Kempton index.

Integrating Multiple Data Types for More Ecological Realism in Trophic Models
The data types integrated in the Celtic Sea EwE model significantly improved its ability to represent the functioning of the Celtic Sea ecosystem by providing a more relevant picture of the food-web structure and by depicting more realistically the drivers that rule its dynamics, and notably the environmental drivers.

Improving the Description of Food-Web Components and Interactions
To ensure that a food-web model correctly represents the temporal and spatial dynamics of an ecosystem, trophic relationships have to be defined as accurately as possible, through the combined use of classical and novel analyses (Pethybridge et al., 2018;McCormack et al., 2019). In this study, the use of outputs of the EcoDiet model (Hernvann et al., in press) applied to the stomach contents and biotracers compiled within the recent and extensive EATME project was advantageous to provide a reliable diet matrix to Ecopath, in better agreement with local data than by the use of literature alone. Additionally, the substantial collection of data compiled for the present study insights on the abundance trends of some food-web components, hence allowing the Celtic Sea EwE model to better capture ecosystem dynamics. Though the integration of abundance timeseries for multispecific groups in EwE was informative, the dynamic of these groups remained poorly understood compared to the dynamic groups under stock assessment. Gathering relevant time-series of fishing effort for the period 1985-2016 for multispecific groups would greatly improve the hindcast. However, the large diversity of gears and métiers from many countries operating in the Celtic Sea makes this endeavor particularly challenging.

Primary and Secondary Production
The integration of primary and secondary production related data in EwE only slightly improved the fit of the Ecosim model. The reason for such a low improvement was the absence of clear variation of primary and secondary production in the time-series. This suggested a relative stability of the hydro-climatic conditions in the Celtic Sea between 1985 and 2016 as the main driver of plankton development. PP and ZHS integration allowed a better hindcast for pelagic functional groups, especially for zooplankton and some planktivorous fish, both for trends and interannual fluctuations, while trends in phytoplankton abundance were poorly reproduced. This poor fit for phytoplankton groups illustrates the difficulty to represent plankton prey-predator interactions and zooplankton grazing pressure in EwE models and the consequence this may have on the prediction of higher trophic levels. Here, high productivity and biomass levels of phytoplankton prevent their control by top-down process and limit their variation in abundance to a very low amplitude. Despite this lack of realism on the dynamic of phytoplankton, the information brought on the bottom-up impact of PP and the development conditions for zooplankton well complemented the insight on top-down impact exerced by higher trophic levels on these secondary producers.
Remote sensing data have been increasingly utilized over the last two decades in marine ecology studies (Chassot et al., 2011). The use of satellite-derived PP estimates is relatively common in marine ecosystem modeling (Christensen et al., 2009;Abdou et al., 2016;Grüss et al., 2016) and has the potential to become a routine procedure in new generations of ecosystem models (Lowerre-Barbieri et al., 2019). The main advantages of satellite-derived PP estimates are that they derived from a direct observation and they are accessible at high spatial resolution, given that regional biogeochemical models are not always locally available. However, satellite-derived PP estimates are only available for recent past (the past two decades) and, so far, they have not been calibrated to provide values for different size classes of phytoplankton (Kramer et al., 2018). In particular, CPR data suggested that the phytoplankton community went through measureable changes in the Celtic Sea over the period 1982-2016: The PCI index, proxy of the total standing biomass, increased, while no significant trend was detected in the diatom and dinoflagellate (both included in large phytoplankton group) counts, which suggested an increase in the biomass of nanophytoplankton between 1985 and 2016.
Despite the relatively short temporal availability of the mesozooplankton niche model of Druon et al. (2019), the endeavor of integrating it in the food-web model turned out to be particularly promising. In particular, estimates from the mesozooplankton niche model allowed us to highlight areas of the Celtic Sea that are favorable to primary production but not to secondary production, and to identify regions with lower food availability for high trophic levels or where the grazing pressure exerted on zooplankton is important. Since the early 2000s, several end-to-end ecosystem models have coupled low and high trophic level models (i.e., have implemented on-line two-way coupling), hence allowing for a better representation of bottom-up and top-down controls within marine ecosystems (Travers et al., 2009(Travers et al., , 2007. However, only a few ecosystem modelers have managed to implement an on-line two-way coupling, and it has been more common to develop alternative ways to link low trophic level and high trophic level dynamics (Libralato and Solidoro, 2009;Piroddi et al., 2017), such as forcing the high trophic level model with outputs from a biogeochemical model (e.g., Halouani et al., 2016). In this context, forcing future ecosystem models by observation-based primary productivity and zooplankton habitat would represent a valuable and conservative alternative approach.

Habitat Modeling
The habitat modeling framework adopted in this study that relied on GAMs was used to assess and predict the response of almost all high trophic level groups to warming and to project the abundance and distribution of each group in the Celtic Sea. The habitat-derived functional response to temperature significantly improved the fit of the Ecosim model. Here again, the improvement of model fit was relatively low, as surface and bottom temperatures appeared to be stable over a substantial fraction of the study period. Still, this improvement was substantial for those functional groups that can cope with extreme temperatures, i.e., which have their temperature optimum outside out of the range of temperature values occuring in the Celtic Sea ecosystem (e.g., plaice, cod and sardine).
Despite their interest in predicting past, current, and future species distributions, habitat modeling approaches have still been rarely coupled to ecosystem models (Grüss et al., 2016;Coll et al., 2019;Moullec et al., 2019b). In this study, we successfully integrated the outputs of a statistical habitat model into a spatial trophic model. The limited complexity of the GAM smoothers employed here, together with the choice of modeling probabilities of presence and not abundances with the GAMs, are in line with the ecological niche theory (Citores et al., 2020). The European IBTS survey data employed in this study encompassed a large range of environmental values, hence providing insights into the potential response of species groups to the warmer conditions that could become the norm in the Celtic Sea by the end of the century. The European IBTS survey dataset recently proved to be highly valuable to study species distributions (Baudron et al., 2020;Moriarty et al., 2020) and characterize ecosystem functioning (Maureaud et al., 2019). The methodology developed by Grüss et al. (2018Grüss et al. ( , 2020 can integrate as many variables as desired (e.g., variables measured at different locations of the water column, percentages of different substrate types) and can be applied to any particular representation of functional groups, including multistanza groups. Thus, the approach of Grüss et al. (2020Grüss et al. ( , 2018 is an interesting alternative to the classical and less flexible automatic approach used in Ecospace ) that imports information from AquaMaps (aquamaps.org), with more limited and much less accurate information (Jones et al., 2012).
Note, however, that habitat models quality remained acceptable but limited for multispecies benthic invertebrate groups, for which presence/absence data was limited here to a few French surveys due to the lack of availability and homogeneity for other surveys. This also stresses the need for improving large scale sampling for these compartments.

Are we Still Missing Something?
While calibrating the Ecosim model, several periods of time with a poorer fit were identified for different functional groups. Early changes in phytoplankton and zooplankton abundances were hardly reproduced by the Ecosim model. Ecosim failed to catch the great decline in small mesozooplankton abundance at the beginning of the time-series, as well as the decline in diatom and dinoflagellate abundances and the increase in nanophytoplankton abundance. Although it could be due to wrong trends in PP estimated by the biogeochemical model in the absence of satellite-derived data, this could also be due to uncaptured environmental effects. This would be expected as our study period starts in 1985, only a few years before the occurrence of abrupt changes in the state of the Celtic Sea ecosystem (e.g., in terms of turbulence, temperature, salinity and plankton abundance) that were identified in a previous study (Hernvann and Gascuel, 2020), and as Ecosim is not capable of tackling such ecosystem regime shift dynamics (Walters et al., 1997;Tomczak et al., 2012). As the stability in the plankton-related time-series did not allow us to explain their dynamics, the oceanographic conditions limiting phytoplankton and zooplankton production may not be the main factors of the regime shift. Several hypotheses have been proposed for identifying the main drivers of the severe drop in mesozooplankton abundance that occurred in the Celtic Sea, including climate drivers or changes in ocean circulation, but the mechanisms responsible for these changes (mediating factor of phytoplankton-zooplankton interactions, direct mortality factor of zooplankton, changes in zooplankton species composition, changes in phytoplankton composition, etc.) remain elusive (Pitois and Fox, 2006). In this context, EwE modeling exercises that estimate artificial "primary production anomalies" may improve the fit of our Ecosim model (Mackinson et al., 2009;Araújo and Bundy, 2012;Bentley et al., 2019). The application of this specific procedure to the Celtic Sea EwE model provided an anomaly time-series significantly correlated to the North Atlantic Oscillation (NAO) index. Nevertheless, these developments were not included in this study as: (i) they did not allow us to better reproduce patterns in plankton abundance in the early years of the period 1985-2016; (ii) they may be confounded with environmental factors that may improve, for instance, fish recruitment; and (iii) they did not provide any indications on the mechanisms at stake.

From Temporal and Spatial to a Fully Spatio-Temporal
Note that a novel spatio-temporal framework, that integrates both temporal and spatial dimensions, has been developed by Steenbeek et al. (2013) and included into the EwE software. The present approach relies on the sequential development of temporal rather than spatial food-web model and does not treat both aspects simultaneously in the same framework. Thus, the changes in functional group abundance predicted temporally do not account for eventual mismatch between prey and predators spatial distribution. The consequences on the overall model predictions may remain limited here due to very limited changes in the distribution of functional groups. However, such novel version of Ecopath spatio-temporal framework should be considered to predict long-term changes in the Celtic Sea in a context of climate change, which would imply substantial shifts in species distributions at various rates according to the functional groups.

Heterogeneity in Biological Production Influenced by Oceanography
The Ecospace model highlighted the heterogeneous spatial patterns of biological production in the Celtic Sea that are mainly driven by contrasting water mixing conditions over the study region. In coastal areas, inner Bristol and English Channel, wellmixed waters allow for early phytoplankton blooms in spring and high primary production during the year depending on nutrient influx/recycling (Pingree et al., 1976;Pingree, 1980;Jordan and Joint, 1984). The seasonally mixed waters spreading over the shelf of the Celtic Sea lead to lower primary production, which occurs later there and at a lower level than in coastal areas (Pingree, 1980;Hickman et al., 2012;Giering et al., 2019). Note however, that PP can be locally enhanced offshore by the combined effects of topography and tidally-induced vertical mixing (through internal waves at the shelf break; Sharples et al., 2007 and water current shear stress in shallow waters of the shelf). The resulting shelfbreak front is well developed in South-West Celtic Sea and finer along the western shelf (Holligan and Groom, 1986;Druon et al., 2001;Joint et al., 2001). In shallower waters, the Ushant front (Pingree et al., 1982) occurs at the interface between the wellmixed and stratified waters of the shelf from the Brittany tip to the Western Channel while the Celtic Sea front (Simpson, 1976) lies between the Bristol Channel and the Irish Sea. According to literature, the most productive areas, and especially productivity fronts, were predicted to concentrate biomass over several trophic levels up to top-predators (Pade et al., 2009;Scales et al., 2014;Druon et al., 2019). Specific functional groups inhabiting the water column displayed aggregated biomass at their location such as mesozooplankton, mackerel and pelagic sharks, as observed at the Ushant front. As about 80% (upwelling) to 90% (oceanic waters) of the PP is transferred to detritus (Raymont, 1980;Joint and Williams, 1985;Parker, 1991), spatial patterns in PP would also affect benthic compartments. The low level of information on benthos in the present study however, prevents any conclusions about these aspects.
The spatial structuration induced by primary production was shown to propagate through the food web. However, except for the above-mentioned functional groups, the spatial contrast in productivity in the Celtic Sea generally dampens with increasing trophic level. The mismatch between higher trophic level species niche and highest primary productivity areas is more important for predators than for prey. In other words, trophic flows progressively operate a biomass transfer from restricted productive areas to lower productive areas in the Celtic Sea. Biomass transfers across the shelf are mainly operated through functional groups characterized by wide ecological niches and large dispersal rates (one input parameter of Ecospace; Romagnoni et al., 2015) such as medium pelagic species mackerel, horse mackerel, pouts and large piscivorous predators in the Celtic Sea. For such functional groups, Ecospace enables the representation of locally enhanced biomass production due to a higher availability of resources when accounting for dispersal toward more suitable environments (Christensen et al., 2014). Thus, through this modeling framework, we highlighted that high fish densities could be maintained in low productive areas of the Celtic Sea through subsidization by prey and predators movements (Polis et al., 1997(Polis et al., , 1996 and transport through residual tidal currents.

Spatial Variation in Functional Groups Assemblages
In agreement with previous in situ studies, the functional groups' distributions predicted by the Ecospace model highlighted the existence of different species assemblages in the Celtic Sea. The species assemblages of the inner Bristol Channel and coastal areas corresponded to those described by Ellis et al. (2000); Kaiser et al. (2004) and Martinez et al. (2013), and these assemblages included herring and sprat as dominant pelagic species, flatfish [plaice, sole, dab (Limanda limanda)], gadoids whiting and pouts, demersal sharks and rays (e.g., catshark Scyliorhinus canicula, rays Raja sp.), gurnards (Eutrigla sp.), small fish (e.g., dragonets Callionymus sp.), large crustaceans (e.g., Maja squinado) and small crabs, pagurids and shrimps. The central Celtic Sea assemblages identified in this study match those reported in Ellis et al. (2013) and Martinez et al. (2013) and are characterized by many gadoid species, including whiting, haddock, pouts and cod, as well as by shrimps. As the present study, Ellis et al. (2013) and Martinez et al. (2013) also characterized the demersal assemblage of the Celtic Sea shelf-edge as being comprised of hake, megrim, anglerfish, pouts, mesopelagic species such as Argentina sp., and benthic cephalopods. The large abundances of pelagic species such as mackerel, horse mackerel, boarfish and blue whiting on the Celtic Sea shelf-edge that were also predicted in this study are consistent with the findings of Reid (2001) and Trenkel et al. (2014).
In agreement with Rees (1999), we found that the Western English Channel displays a particularly diverse benthic assemblage that includes large crustaceans and commercial bivalves [Kaiser et al., 1998;e.g., king scallop (Pecten maximus), queen scallop (Aequipecten opercularis)], as well as suprabenthos (Vallet and Dauvin, 1998). The preponderant community of pelagic species in the Western English Channel has been carefully screened for decades (Southward et al., 1988;Hawkins et al., 2003). The spatial patterns of abundance of most demersal species predicted in this study were consistent with those derived from English Celtic Sea groundfish survey data (Warnes and Jones, 1995;Tidd and Warnes, 2006), even though the latter were not employed for statistical habitat modeling in this study.

Heterogeneous Species Assemblages Result in Different Local Food Webs
The different species assemblages of the Celtic Sea identified in this study provided some new insights to the possible spatial structure of the food web in the Celtic Sea region. While some functional groups are widely distributed, some others are only present in some specific areas of the Celtic Sea and their absence or negligible biomass at particular locations may be interpreted as reflecting local reorganizations of the food web. Thus, the complex Celtic Sea food-web structure could be interpreted as a meta food web made of local, environmentallydriven, food webs connected through space by functional groups with larger dispersion capacities (Kortsch et al., 2019(Kortsch et al., , 2015. Ecospace reveals to be a useful tool to understand the spatial structure of food webs as a complement to more qualitative approaches for characterizing the trophic structure of marine ecosystems (Albouy et al., 2014;Pellissier et al., 2018;Kortsch et al., 2019). Spatial indicators can be produced with Ecospace to determine the spatial structure of the food web in relation to the environment. Our results suggest that shallowest and least saline (affected by rivers) waters of the Celtic Sea, together with the deepest waters of the region, are characterized by a low number of compartments and groups with relatively low variability in terms of TL. The two above-mentioned areas are, however, radically different in terms of functioning. The latter is dominated by a specific trophic pathway where toppredators (i.e., hake, megrim) in high density feed on abundant pelagic species (i.e., blue whiting, boarfish). The former, in deeper environments, is essentially represented by low trophic level benthic and pelagic functional groups with non-diversified diets. The distinction made between the two areas does not exclude the existence of important benthic-pelagic coupling in both areas, but this coupling may be due to the reliance on pelagic primary production rather than to exchanges between upper trophic level species. Intermediate-depth areas of the Celtic Sea were found to be dominated by pelagic functional groups, yet they displayed large concentrations of benthopelagic fish feeding on the benthos and pelagos. Finally, the central Celtic Sea was found to be dominated by predatory and generalist functional groups that feed on both pelagic and benthic compartments, which may render this specific area more stable than the rest of the Celtic Sea. Spatialized indicators highlighted that boundaries between these different areas constituted ecotones, characterized by higher functional group and trophic diversity as they suited to a wider range of species.

The Limited Perception of the Mean Annually-Averaged Ecosystem Functioning
The spatio-temporal dynamic model of this study represents the trophic functioning of the Celtic Sea ecosystem averaged over a year. We proceeded this way because of data availability for model calibration (i.e., survey-derived abundance indices, stock assessment products, and catch statistics). Averaging trophic interactions can be questionable, as trophic functioning over the year is strongly related to seasonal variations in plankton production. These variations are mainly associated with hydrographic features and water mixing processes. In addition, seasonal movements, and subsequent changes in spatial distributions and local abundances have been highlighted for both pelagic [widely distributed pelagic fish (Trenkel et al., 2014) and coastal to a lower degree (Brophy, 2002)] and demersal communities (Dunn and Pawson, 2002;Persohn et al., 2009;Neat et al., 2014). These movements can notably be related to feeding, spawning and overwintering.
Seasonal changes in distributions would explain the variability of fish diets in the Celtic Sea reported in Trenkel et al. (2005) and Chassot et al. (2008). For instance, large Celtic Sea predators consume more blue whiting during the summer and more mackerel during the winter, whereas the consumption of horse mackerel is relatively constant throughout the year. However, modeling changes in fish spatial distributions at a finer temporal scale remains a challenging issue given the data that are available to us.

More Insights Into the Dynamic of a Shelf Ecosystem Under Multiple Environmental and Anthropogenic Pressures
The Celtic Sea Ecosim model best explained the food-web dynamics when it integrated trophic interactions, fishing, environmental conditions controlling plankton production, and warming as a driver of the productivity of higher trophic level groups. The way the Ecosim model fitted to data when different data types were integrated provided a better understanding of the relative impacts of fishing and environmental changes on the Celtic Sea ecosystem during the three last decades.

Fishing on the Frontline
Previous work based on catch and effort reconstruction highlighted the large impacts that fishing has had on the Celtic Sea ecosystem since 1950 (Pinnegar et al., 2002;Guénette and Gascuel, 2012;Hernvann and Gascuel, 2020). In particular, the dramatic increase in fishing pressure that has occurred in the Celtic Sea region over several decades may have drastically reduced the biomass of large predators and may have triggered the deep alteration of the food-web structure through trophic cascades (Hernvann and Gascuel, 2020). The impacts of fishing in the Celtic Sea ecosystem suggested by the present study concur with the conclusions of previous EwE modeling works in the ecosystem by Bentorcha et al. (2017) and Moullec et al. (2017). At the beginning of our study period, in 1985, the biomass of main large bentho-demersal fish and predator species had already been severely impacted by fishing and pelagic and invertebrate fisheries started developing in the Celtic Sea region. Then, between 1985 and 2016, fishing kept on severely impacting the Celtic Sea ecosystem. Despite the integration of environmental variability in the Ecosim model, our findings suggest that fishing had more profound impacts on the Celtic Sea ecosystem than environmental changes.
Compared to previous versions of the Celtic Sea EwE model, the present Ecosim predictions were more accurate due to the use of new survey data (EVHOE and UK-WCGFS) which allowed for a better description of the indirect effects of fishing in the Celtic Sea EwE model. The divergent abundance time-series for smaller and larger fish groups, which increased and decreased, respectively, over the period 1985-2016, were better identified in this study than in the studies that employed previous versions of the Celtic Sea EwE model. The analysis of the mortality estimates suggested that the driving factor of trends in mortality was predation release on smaller species, due to both the continuously increasing pressure exerted on large predatory fish and the reallocation of fishing effort on medium fish in response to the depletion of to large predators. These results concur with the hypotheses formulated in community size-structure studies Shephard et al., 2012Shephard et al., , 2011. Though this phenomenon was only described for demersal species, the large decline in the biomass of horse mackerel in the beginning of the 1990s predicted in this study (a 6-fold decline within only a 10year time frame resulting from a large increase in fishing effort on horse mackerel) may have also contributed to the increasing biomass of smaller species. This may probably be the case (horse mackerel leaving a vacant niche) for other small pelagic species, such as boarfish to boom in the early 2000s, as observed by Blanchard and Vandermeirsch (2005) and Coad (2012) and also correlated to changes in the temperature. Thus, fishing appears as the main driver of changes in species assemblages in the Celtic Sea, especially leading to a decrease in the diversity of diets and a general decline in the mean trophic level in the area.
Finally, the biomass recoveries initiated in the Celtic Sea over the recent period were particularly noticeable, and these biomass recoveries seemed to allow for a progressively increasing dominance of demersal species in the Celtic Sea ecosystem and a rise of the predation pressure exerted on smaller demersal fish, crustaceans, and cephalopods as highlighted by the biomass decline of the latter and the increase in the HTI index.

Climate Change in Ambush
The importance of the fishing driver in the Celtic Sea contributed to mask the impacts of environmental changes on the ecosystem. Indeed, many abundance trends of exploited groups reflected stock depletions and recoveries independently of warming effects. The dominance of the fishing driver in observed time series is consistent with the lack of evidence of any warming effect in the Celtic Sea from studies that analyzed species biomass trends or community size-structure (Genner et al., 2004;. However, through the use of GAM predictions and temperature forcing, we were able to account for the cumulative impacts of fishing and environment changes in EwE and to disentangle the relative contributions of these different drivers of change. The mid-1990s were characterized by a net warming of the Atlantic Ocean that was attributed to a shift of the Atlantic Multidecadal Oscillation (AMO) from a negative to a positive phase, which was likely amplified by anthropogenic climate change (Ting et al., 2009). As a response to this large-scale warming events, both sea surface and bottom temperatures in the Celtic Sea increased. The Celtic Sea EwE with Ecospace model suggests that these temperature changes led to a decrease in the productivity of most of Boreal species, including herring, sprat, cod, whiting and haddock, and to an increase in the productivity of functional groups with warm water affinity, particularly species with extreme thermal preferenda such as sardine, sea bass, mixed medium pelagic fish [mainly anchovy (Engraulis encrasicolus)] and cephalopods, and to a lesser extent, widely distributed pelagic species such as mackerel, horse mackerel and boarfish. In future years, these trends in species productivity may impair the recovery of some current commercially important species (i.e., gadoids) and oblige fishers to target less exploited warm-water species. Additionally, the species distributions predicted by Ecospace indicated a relative concentration of Boreal species' biomasses in a smaller area around the mid-1990s. This modification of Boreal functional groups' distributions persisted until the end of the study period and resulted from a contraction of these groups' niches induced by warming. If the scale of these changes was limited, it confirmed the response of these species to climate change, which was barely detectable in temporal trends. The niche contraction of demersal Boreal species in the south of the Celtic Sea was associated to an increase in functional groups diversity, related to the largest proportion of biomass represented by species with warmer water affinities. The latter, the trophic diversity declined in response to change in assemblage, being mainly benthivorous and of low trophic level.
These responses of functional groups' productivity and distribution to ocean warming predicted in this study were consistent with the findings of multiple fish community studies that were conducted across Europe for a large proportion of both demersal and pelagic species. In such studies, species were shown to respond to warming in the 1990s (Simpson et al., 2011;Montero-Serra et al., 2015). Changes were mostly correlated to species' preferred temperature and central latitude of their range so that the abundance or recruitment of cold water species decreased while that of warm water, widely distributed species increased (Poulard and Blanchard, 2005;Simpson et al., 2011;Auber et al., 2017Auber et al., , 2015Bentley et al., 2020). If no homogeneous poleward shifts of the communities were highlighted over the European shelf (Dulvy et al., 2008;Simpson et al., 2011), the spatial distribution of numerous species was affected by warming (Baudron et al., 2020), some displaying displacements of their center of distribution and others showing migration to deeper waters (Dulvy et al., 2008;Punzón et al., 2016). The potential changes in the distribution of species located south of the Celtic Sea or already present in low abundance in the Celtic Sea are responsible for the increase in species diversity observed by Hiddink and ter Hofstede (2008) and ter Hofstede et al. (2010).

The Relative Stability of the Celtic Sea Ecosystem
When spatializing fish communities' responses to warming over the European shelf, Simpson et al. (2011) pointed out that the Celtic Sea was less impacted than the neighboring areas. This result may be linked to the lower variability of Celtic Sea thermal conditions after the mid-1990s that resulted from the particular hydrographic situation of the Celtic Sea (i.e., Gulf Stream influence, Hughes et al., 2017;Tinker et al., 2020). However, the amplitude of warming in the mid-1990s suggested that local thermal conditions were not the only reason for the weaker response of the Celtic Sea ecosystem to warming compared to neighboring ecosystems. Ecospace results suggest that the greater stability of Celtic Sea communities may be due to the complex spatial structuration and high biodiversity of the Celtic Sea ecosystem. According to various theoretical and modeling studies, the meta-community organization of food webs may act as a stabilizer of trophic dynamics (McCann et al., 2005;Gravel et al., 2011;Mougi and Kondoh, 2016). The spatial arrangement of environmental gradients in the Celtic Sea is complex compared to other ecosystems where gradients are more correlated together (e.g., the Bay of Biscay, the North Sea). Additionally, though they are developed from collated data from multiple origins, harmonized seabed maps at the European scale indicate that substrate patchiness is higher than in neighboring ecosystems (e.g., the North Sea, the Eastern English Channel). The spatial structuration of the Celtic Sea ecosystem may favor a stable meta-community organization and allow for more resistance to the impacts of environmental changes. The complex spatial structure of the Celtic Sea may also reduce the probability of co-occurrence of "autochtonous" and a "allochtonous" species in the ecosystem, hence lowering the risks of competition between these two types of species if they have similar trophic niches. Finally, the substantial contribution of both benthic and pelagic pathways to matter flows in the Celtic Sea ecosystem may also participate in the stability of the ecosystem (Rooney et al., 2006;Blanchard et al., 2011).
We also found that the Celtic Sea is relatively stable in terms of basal productivity. As reflected by biogeochemical models and satellite-derived time-series, the oceanographic conditions favoring phytoplankton and zooplankton abundances over the study period remained overall constant. This may be due to the mechanisms involved in the production of mesoscale structures in the Celtic Sea that favor plankton production. Indeed, a substantial part of total primary production in the Celtic Sea is due to the tidally-induced semi-permanent mixing fronts that occur in the central part of the eastern Celtic shelf, i.e., the Ushant Front and the Celtic Sea Front. Contrary to other ecosystems in which productivity fronts are largely dependent on climatic conditions, tidal mixing fronts in the Celtic Sea can form simply under the mechanical action of wave activity, resulting in a more regular production in the Celtic Sea region, which is consistent with the production of mesozooplankton (Druon et al., 2019).

A Step Further Toward Ecosystem-Based Management?
This new EwE with Ecospace model of the Celtic Sea ecosystem represents a significant step forward in the integration of ecological realism in trophic modeling of the Celtic Sea ecosystem (Guénette and Gascuel, 2009;Bentorcha et al., 2017;Moullec et al., 2017). In the Celtic Sea, highly mixed-species fisheries catch a large diversity of fish assemblages and employ a large diversity of fishing techniques and strategies (Moore et al., 2019) that can result in highly different catch compositions for the same location over time Dolder et al., 2018). In this context, resource mapping and optimization of fishing practices in time and space have been identified as a critical need to reduce the ecosystem impacts of fishing in the Celtic Sea, e.g., by reducing unwanted catch (Dolder et al., 2018;Calderwood et al., 2019;Reid et al., 2019;Robert et al., 2019). The Celtic Sea Ecospace model reproduced the patterns observed during fisheries-independent surveys but also clearly identified the spatial structuration and species assemblages that were identified in some spatial studies of fisheries catch, based on logbooks and vessel monitoring systems (Gerritsen et al., 2012;Mateo et al., 2017).
The EwE with Ecospace model developed in this study represents a valuable tool for the spatial management of fishing activities in the Celtic Sea. In future studies, the Celtic Sea Ecospace model could be employed to produce multiple spatial food-web indicators to describe ecosystem status as part of ecosystem-based management (Tam et al., 2017), or to identify critical areas of the Celtic Sea in terms of diversity and functionalities (Babcock et al., 2005). As other ecosystem models, our EwE with Ecospace model could be also used to inform spatial zoning in the Celtic Sea (e.g., Kaplan et al., 2012;Abdou et al., 2016;Halouani et al., 2016) and other management regulations (e.g., relative to discarding practices, Pennino et al., 2020). Such framework can therefore provide great support to policies such as the Common Fisheries Policy (e.g., spatial fishing capacities) and the Marine Strategy Framework Directive (e.g., descriptor D3 commercially exploited fish, D4 food webs, e.g., Piroddi et al., 2015;Piroddi et al., submitted).
Although fishing appeared to be the main driver of changes in the Celtic Sea ecosystem over the three last decades, this may likely not be the case anymore in the future. While fishing pressure in the Celtic Sea region has been substantially reduced over the recent years and is still currently (slowly) declining, the relative stability of environmental conditions observed over the last two decades may likely not hold; projections of coupled physical and biogeochemical models under various CO 2 emission scenarios predict a global increase in surface and bottom temperatures and a decline in primary production in the Celtic Sea ecosystem (Kay et al., 2018). In this context, similarly to previous temporal or spatio-temporal ecosystem modeling studies Corrales et al., 2018;Moullec et al., 2019a), the calibrated the Celtic Sea Ecospace model would be a powerful tool to project the potential impacts of climate change on the regional living resources and the fisheries that these resources support.

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. The Celtic Sea Ecopath model at the core of this study can be accessed from EcoBase, the open-access database of Ecopath with Ecosim (EwE) models published worldwide in the scientfic literature. EcoBase is an information repository of EwE models freely accessible online at http://ecobase.ecopath.org/ or via the EwE software.

AUTHOR CONTRIBUTIONS
P-YH designed the study, performed the main statistical analyses, and built and ran the ecosystem model. P-YH, J-ND, MR, DG, and IP provided and gathered the data used in the study. AG, J-ND, CP, and DG helped to and participated in designing and applying the habitat and trophic modeling approaches integrated in the present modeling framework. DG, MR, and DK contributed to the design of the integrated modeling framework. DG, MR, DK, AG, J-ND, CP, and IP contributed to results analysis and interpretation. This work was part of a Ph.D. project designed by DG, MR, and DK. P-YH led the drafting of the manuscript with the contributions and revisions from all the authors. All authors contributed to the article and approved the submitted version.

FUNDING
This work was part of a Ph.D. funded by the Région Bretagne and the French Institute for the Exploitation of the Sea (IFREMER). The temporary stay of the principal investigator at the Joint Research Center of the European Commission as part of this research project was funded by an international mobility grant from IFREMER.