Overexploitation, Recovery, and Warming of the Barents Sea Ecosystem During 1950–2013

The Barents Sea (BS) is a high-latitude shelf ecosystem with important fisheries, high and historically variable harvesting pressure, and ongoing high variability in climatic conditions. To quantify carbon flow pathways and assess if changes in harvesting intensity and climate variability have affected the BS ecosystem, we modeled the ecosystem for the period 1950–2013 using a highly trophically resolved mass-balanced food web model (Ecopath with Ecosim). Ecosim models were fitted to time series of biomasses and catches, and were forced by environmental variables and fisheries mortality. The effects on ecosystem dynamics by the drivers fishing mortality, primary production proxies related to open-water area and capelin-larvae mortality proxy, were evaluated. During the period 1970–1990, the ecosystem was in a phase of overexploitation with low top-predators’ biomasses and some trophic cascade effects and increases in prey stocks. Despite heavy exploitation of some groups, the basic ecosystem structure seems to have been preserved. After 1990, when the harvesting pressure was relaxed, most exploited boreal groups recovered with increased biomass, well-captured by the fitted Ecosim model. These biomass increases were likely driven by an increase in primary production resulting from warming and a decrease in ice-coverage. During the warm period that started about 1995, some unexploited Arctic groups decreased whereas krill and jellyfish groups increased. Only the latter trend was successfully predicted by the Ecosim model. The krill flow pathway was identified as especially important as it supplied both medium and high trophic level compartments, and this pathway became even more important after ca. 2000. The modeling results revealed complex interplay between fishery and variability of lower trophic level groups that differs between the boreal and arctic functional groups and has importance for ecosystem management.


INTRODUCTION
Fisheries and climate have been emphasized as major drivers of energy flows in large marine ecosystems (LMEs) (Araujo and Bundy, 2012;Link et al., 2012), and understanding how these drivers interacts and shapes ecosystems is a major challenge and essential to manage these ecosystems. It is important to investigate how these drivers interplay in high-latitude ecosystems, such as the Barents Sea. The Barents Sea (BS) is intensely exploited and profoundly impacted by climate variability. It is characterized by a strong temperature gradient with boreal and sub-arctic conditions in the southwest and high-arctic conditions in the northeast (Loeng, 1991;Smedsrud et al., 2010), and advected heat, nutrients, and biota, along with seasonally migrating fish, seabirds, and mammals strongly impact ecosystem structure and production (Hunt et al., 2013). The trophic structure of the BS ecosystem is similar to other northern high latitude shelf ecosystems (Gaichas et al., 2009;Eriksen et al., 2017).
Climate variability may affect marine ecosystems through effects on primary and secondary production, fish recruitment variability, growth and shifts of populations distribution range (Nilssen et al., 1994;Ottersen and Loeng, 2000;Fossheim et al., 2015). In the BS, a period of warm climate between 1920 and 1960 was followed by a cold period from ca. 1960 to 1980 and then by a period of warming after the 1980s (Loeng and Drinkwater, 2007). Temperature variability has affected recruitment to the commercial fish stocks Norwegian spring spawning herring, Northeast Arctic cod and haddock with larger year-classes produced in warmer years (Ottersen and Loeng, 2000;Sundby, 2000). The warming of the BS ecosystem since the early 1980s (Johannesen et al., 2012) has resulted in northwards shift in distribution and increasing abundance for boreal fish species and a decrease for arctic fish species Kortsch et al., 2015) and an increased importance of benthic invertebrate species with affinity for warmer waters (Jørgensen et al., 2019). After 1980, temporal fluctuations in population sizes have been observed at several trophic levels, e.g., krill, northern shrimp, capelin, and seabirds (Johannesen et al., 2012;Fauchald et al., 2015;Gjøsaeter et al., 2015). Relationships and energetics of major stocks of top-predators, such as Northeast Arctic cod, minke whales, and harp seals have also changed Fauchald et al., 2015).
The change in climatic conditions call for an effort to use and integrate available information to understand the underlying drivers for the ecosystem changes, and ecosystem modeling is a common tool to synthesize quantitative information into a coherent system. This is particularly important in speciesrich systems with complex (e.g., with considerable advection and migration) pathways for impact where statistical modeling may struggle. Previously published food web models of the BS and Norwegian Sea have been fish-centered with relatively few lower trophic level groups and benthic invertebrate groups (Blanchard et al., 2002;Dommasnes et al., 2002;Hansen et al., 2016;Skaret and Pitcher, 2016;Bentley et al., 2017). Arctic and sub-Arctic ecosystems, however, are well-known for the strong role of seafloor communities in regulating carbon cycling pathways (Kędra and Grebmeier, 2021). Based on a dynamic mass-balance model (Ecopath with Ecosim-EwE) and future warming scenarios for the Norwegian and the BS, Bentley et al. (2017) suggested that the biomasses of widely migrating pelagic species, such as mackerel and blue whiting, are expected to increase with future rising ocean temperature. There is some evidence that effects of climate variability on the ecosystem in the BS is largely through bottom-up effects on lower trophic level groups that propagate to higher trophic level groups (Johannesen et al., 2012;Dalpadado et al., 2014Dalpadado et al., , 2020. Fisheries and climate change may act synergistically with each other and/or with other anthropogenic disturbances (Fogarty et al., 2008;Hsieh et al., 2008). Exploited populations may be less resilient to climate variability than unexploited populations due to more truncated age structure and diversity in life history traits (Fogarty et al., 2008;Hsieh et al., 2008). A better understanding of how exploitation and climate variability influence the ecosystem dynamics will support management of marine resources.
Aggregating species into trophic groups may mask complex species interactions and influence calculations associated with food webs and interspecific competition (Thompson and Townsend, 2000). Therefore, to analyze trophic interactions and impacts of harvesting and climate variability in this study, we parametrized an Ecopath food web model for the BS with both Atlantic boreal and Arctic groups, and with a high resolution of lower trophic level groups. This model was evaluated and fitted to time-series of biomasses and fisheries data (catches, fishing mortalities). The main objectives of this study were to evaluate how changes in exploitation and climate have affected ecosystem structure, metrics, and properties of the BS ecosystem during the period 1950-2013. The specific aims were to; (i) quantify carbon flow pathways and production by ecological compartments, (ii) investigate whether past exploitation have reduced biomass and productivity of functional groups and led to trophic-cascaderelated effects, and (iii) assess if climate variability affected the ecosystem productivity and if boreal and arctic groups were affected differently.
We will use available updated information on trophic linkages to parametrize a highly resolved Ecopath with Ecosim food web model in the time-period 1950-2013. The model was fitted and calibrated to group-specific time-series of biomasses and catches and forced by environmental drivers, such as fishing mortality, primary production, and capelin larvae mortality proxies. The effects of exploitation and climate variability was evaluated by the ecosystem metrics and properties produced by this model. We discuss how our findings may support an ecosystem based management of the BS.

Study Area, Environment, and Data Collection
The BS is a high latitude LME, covering an area of 2.01 million km 2 (Skjoldal and Mundy, 2013), extending from the Norwegian Sea and eastwards to Novaya Zemlja and northwards from the coast of Norway and Russia to about 80 • N (Drinkwater, 2011; Figure 1).
Water circulation and currents in the BS are strongly influenced by the bottom topography. The shallowest areas are found around Spitsbergen Bank and in the southeastern part with depths <50 m (Loeng, 1991). The deepest area (deeper than 400 m) is found in the western part where the main influx of relatively warm Atlantic (T > 2 • C) and Coastal (T > 3 • C) waters enters the BS (Loeng and Drinkwater, 2007). Cold Arctic (T < 0 • C) water penetrates the system from the east and north (Hunt et al., 2013). The Polar Front is a transition zone between the warmer boreal southern part and the colder Arctic northern part (Figure 1; Fossheim et al., 2015). During winter, the edge of the seasonal ice cover was normally found just north of the Polar Front (Smedsrud et al., 2013). The ice cover varies both seasonally and inter-annually (Wassmann et al., 2006a;Smedsrud et al., 2013), with maximum coverage typically in March-April and the minimum coverage in August-September (Drinkwater, 2011). The climatic gradient within the BS is reflected in the distribution of organisms (Andriyashev and Chernova, 1995;Jørgensen et al., 2015;Renaud et al., 2018). Boreal fish species have generally expanded northwards at the expense of arctic species during the recent warm period .
Data to parametrize, drive and evaluate the Ecopath and Ecosim models were collected from literature and published data sources from the BS (Supplementary Appendices 2, 4). In cases were data from BS were not available, data from other similar areas were used (Supplementary Appendix 2).
FIGURE 1 | Map of Barents Sea large marine ecosystem. Borders of the ecosystem are shown by red lines Based on (https://www.pame.is/projects/ecosystem-approach/arctic-large-marine-ecosystems-lme-s). The Kola transect stations 3-7 for hydrographic monitoring are shown as black dots. Location of the polar front is shown by blue line.

Model Description
The Ecopath model tracks Carbon as the mass unit to reflect the varying organic carbon content of functional groups. Organic carbon has a much stronger relationship to energy than wet mass (Salonen et al., 1976) and carbon is commonly used as unit in Ecopath models with emphasis on lower trophic levels (Tomczak et al., 2009).
The Ecopath model consist of two master equations (Christensen et al., 2005). The first equation describes how production for a FG i is split into various components Equation 1 can be written as; where P i is the production of group i, Y i is the catch of group i, M2 i is the predation mortality rate on groups i, B i is the biomass (g C m −2 ), (P/B) i is the production/biomass ratio of group i, (Q/B) j is the consumption/biomass ratio of predator j, DC ji is the proportion of prey group i in the diet of predator j, Y i is the catch, E i is net emigration, BA i is the biomass accumulation and B i (P/B) i (1 -EE i ) is other mortality of group i. EE i is the ecotrophic efficiency describing the proportion of production of a group that is consumed within the model. Within each FG i, energy balance is ensured using the equation The BS Ecopath model for year 2000 comprises 108 functional groups (FG) of which 19 groups were multi-stanza groups ( Table 1 and Supplementary Appendices 1-3). Multi-stanza groups contain a set of biomass groups representing life history stages or stanzas for species that have complex trophic ontogeny . Species were grouped in FGs based on their similarities in diet composition, production/biomass ratio (P/B) and consumption/biomass ratio (Q/B), predators and predatory mortalities.  (Figure 1) has been considered as a good indicator for the temperature variability in the BS (Stige et al., 2010).
The annual average Ecopath model representing year 2000 was parameterized by the best available literature data (Supplementary Appendix 2). Data for biomass, catches, diet composition, production/biomass, consumption/biomass, and assimilation efficiency, were mainly derived from the BS. For data-rich groups, such as commercially exploited fish, northern shrimp and top-predators, where abundance or biomass are regularly monitored, biomass and catch time-series could be calculated (Supplementary Appendix 4 Part A,B). For less surveyed groups, diet composition and other parameters were averaged on longer time-periods. Catch data were retrieved from official statistics (ICES) and time-series from publications. After parametrizing and balancing the year 2000 model, it was modified to a year 1950-model by entering year-specific values for biomasses and catches for groups with known data and the 1950model was balanced (Supplementary Appendix 2). The modeled and observed biomass and catch trajectories of Ecosim based on the 1950 Ecopath model as "starting point" were compared to evaluate the parametrization. P/B was often calculated from data on total annual mortality rate (Z, yr −1 ) as P/B = Z . Diet composition data for upper trophic levels (TL ≥ 3) FGs (fish, mammals, and birds), were mostly based on stomach analysis. When several data sets of diet composition were available for a functional group, diet proportions were averaged to provide initial input values (Supplementary Appendix 4 Part C). Diet compositions were converted from wet mass units in the original data sets to carbon mass using group specific carbon/wet mass factors (C/WW) (Supplementary Appendix 2).
To assess uncertainty in the input values for biomass, P/B, Q/B, diet and catch, pedigree scores were allocated to each input value using the system for pedigree indices integrated into EwE (Supplementary Appendix 1 Tables 1-3, 1-4). Pedigree indices are uncertainty scores assessed by the modeler for each input value in Ecopath and were based on either measured uncertainty or assessed from the type of data and the source of the input value. A 95% confidence interval is associated with each index value. 1 http://www.pinro.ru/labs/hid/kolsec22.php?lang

Balancing Ecopath Models
The BS ecosystem is a well-studied and ecological knowledge has accumulated. The Ecopath modeling allows us to evaluate the compatibility of the input data and identify uncertainty in the input data and in the model output because if the Ecopath model is unbalanced, the input data are not compatible. Before and after balancing the Ecopath models, the pre-balance procedure (PREBAL) was used to check if input values were within accepted ecological constraints (Link, 2010; Supplementary Appendix 4 Part D). Biomass, Q/B and P/B decreased with increasing TL, except for some mammals and bird groups that have high Q/B-values.
Balancing the Ecopath models was done manually by checking that EE ≤1 for the mammal, bird and fish groups where biomass input data were available. For groups where the biomass was estimated by the Ecopath model from consumption by its predators and catches, it was checked that biomass values were within the range reported in the literature (Supplementary Appendix 2). Gross efficiencies (P/Q) were checked to be in accordance to the literature, and it was checked that respiration/assimilation were <1. The models were constructed using version 6.6.3.16996 of EwE. 2 To simplify presentation of model results, some output values were presented for aggregated compartments ( Table 1).

Dynamic Simulations Using Ecosim
A dynamic Ecosim model was constructed based on the 1950 Ecopath model. In Ecosim, the biomass growth rate of functional group i is expressed as Where (P/Q) i is the gross efficiency, M i is the natural mortality not caused by predation, F i is the fishing mortality rate, I i is the immigration rate, e i is the emigration rate and B i is the biomass (Coll et al., 2009). The consumption rates in Ecosim (Q ji ) are based on the "foraging arena theory" where the biomass of prey i is divided into a non-vulnerable and a fraction that is vulnerable to predation (Walters and Korman, 1999) and vulnerabilities (v ij ) express the maximum increase in predation mortality when predator abundance is high. Vulnerabilities (v ij ) and a number of other parameters affect consumption rate (Q ij ) of a group i preyed by a predator j (Christensen et al., 2005).
Where a ij is the effective search rate on a prey j, T i is the relative feeding time of prey, T j is the relative feeding time of the predator, S ij is a seasonal or long-term forcing function, M ij is a median function and D j represent effects of handling time.
Low vulnerabilities close to 1.0 are associated with a low increase in predation mortality when predator biomass increase and vice versa for high vulnerabilities. In Ecosim, the additional parameters that limit the consumption rates (Eq. 5) (Christensen et al., 2005) were set to default values during the simulations.

Ecosim Time-Series Fitting, Model Calibration, and Cross-Validation
In fitting and calibrating ecological models, there is a potential risk of overfitting models (Wenger and Olden, 2012). To evaluate Ecosim models, the model fit to time-series data, model behavior and the model's ability to predict a part of the time-series not used in the fitting (cross-validation) were considered. We performed a cross-validation where the available time-series were split into a part (ca. 75%) used for model fitting  and a part used for testing the predictability of the models (1997-2013) (Bergmeir and Benítez, 2012;Wenger and Olden, 2012). A model's ability to predict for a time-period not used in fitting informs about its transferability. Transferability has been emphasized as an important aspect of model evaluation and cross-validation has been suggested as a method to assess transferability and reduce risk of overfitting (Wenger and Olden, 2012).
Ecosim models were fitted to time-series for the period 1950-1996 and calibrated by estimating predator-prey vulnerabilities (v ij ). The number of vulnerabilities that can be estimated is equal to the number of independent time-series minus one (Scott et al., 2016). A common approach in Ecosim calibration/fitting to time-series has been to fit a spline-function which is considered a proxy for primary production, and then relate this spline function to environmental proxies, such as NAO-indices and water temperature (Skaret and Pitcher, 2016;Bentley et al., 2017). However, we found it more appropriate to include a welldocumented relationship for PPR as an environmental driver. Before the Ecosim model was fitted to time-series, the time-series were categorized into "forcing" (forcing the model to time-series values), "absolute" (absolute values were used), "relative" (relative values, such as catch per unit effort were used). A total of 84 timeseries were used as input in the time-series fitting, including timeseries on absolute biomass (n = 32), relative biomass (n = 15), forced biomass (n = 3), forced catch (n = 9), catches (n = 9), fishing mortality (n = 14), harvesting effort (n = 2). This amounts to a total of 56 time-series on absolute and relative biomasses and catches that were used to calculate sum of squares and a potential maximum of 55 vulnerabilities could be fitted. A lognormal error distribution was assumed minimizing the sum of squares of log observed values from log modeled values (Christensen et al., 2005). There were equal weights of each time-series, thus the absolute scale of time-series values did not influence the sum of squares. The mesozooplankton biomass time-serie mainly comprising the FGs medium sized copepods, large calanoids, and small copepods, was not used in the fitting but was compared to the model output of the sum of biomass of its FGs. In addition, various environmental time-series were used as driving forces for the model (Supplementary Appendix 4 Part A,B).
An automatic step-wise fitting procedure was used to calibrate the Ecosim-model to observed time-series for biomass fishing mortalities and catches (Scott et al., 2016). This procedure statistically estimates how much fishery time-series, trophic interactions (predator-prey vulnerabilities) and environmental time-series contribute to model fit. The stepwise fitting procedure constructs a series of model permutations with increasing number of estimated vulnerabilities (v) and determines which combination of vulnerabilities gives the best statistical fit using sums of squares (SS) and Akaike's Information Criterion modified for small samples (AICc) as criteria to select the most parsimonious model (Akaike, 1974;Kletting and Glatting, 2009). The fitting and model evaluation procedure include several steps. (i) The sensitivities of SS for the chosen number of model vulnerabilities for each predator-prey interaction or predator were calculated and FGs with the highest sensitivities were selected. (ii) It was searched iteratively for values of vulnerabilities among the selected vulnerabilities to minimize the SS for the period 1950-1996. (iii) Plots of model-fitted and observed timeseries of biomasses and catches and SS for separate FG were visually inspected to evaluate model fit to observation data and assess if the model behavior was credible (Mackinson, 2014;Heymans et al., 2016). To assess the effects of fisheries, environmental drivers and trophic factors (vulnerabilities), alternative models were tested. Alternative models were fitted without fishery data (baseline models with no catch or fishing mortality), with fisheries data and environmental forcing timeseries and with and without estimating vulnerabilities resulting in a SS and an AICc-value per model. (iv) In the model prediction runs for the period 1997-2013, the predictability of models was assessed by calculating the SS for model output biomass and catches and the corresponding observed data.
To assess the effect of climatically forced phytoplankton primary production (PPR), in Ecosim, two alternative forcing time-series were tested; a constant PPR-proxy and a PPRproxy based on the relationship between phytoplankton primary production and open-water area (Supplementary Appendix 4 Part A). A capelin larvae mortality proxy was calculated based on the relationship between biomass of small herring and capelin larvae mortality rate (Supplementary Appendix 4 Part A) and the proxy was used to force mortality rates of capelin (0-2) in model fitting.
The possible effects of change in ice-coverage on ice-algae primary production were tested by modifying model M10 by forcing ice-algae biomass directly by the ice-cover in a model M11 run (

Monte Carlo Simulations and Model Evaluation
For model(s) considered to have most support assessed by the stepwise fitting for the period 1950-1996, inspection of model behavior and test of predictability for the period 1997-2013, Monte Carlos simulations (MCS) were run to assess uncertainty in output values from Ecopath and Ecosim. In the MCS, input values were randomly sampled from uniform distributions with the width of the distributions corresponding to pedigreespecified input uncertainty level for biomasses, P/B and Q/B values (Supplementary Appendices 1, 2), and the MCS routine included 200 successful trials with balanced models. Each trial had up to 10,000 runs where Ecopath input parameter values were drawn and it was tested if the resulting Ecopath model was balanced. To evaluate uncertainty and compare model outputs with observed data, the 0.025 and 0.975 percentiles of the MCS outputs were calculated.
To evaluate the model fit, Taylor diagrams were used to simultaneously visualize the correlation (Pearson) between observed and modeled time-series, the root-mean-square difference (RMS) and the ratio of the standard deviations of the simulated and the observed time-series (RSD) (Taylor, 2001).

Ecosystem Indicators
It has been advised to use a variety of indicators at the community level to detect ecosystem impacts of fishing (Fulton et al., 2005). The indicators should include groups directly impacted by the fishery, charismatic groups with slow dynamics and response (e.g., mammals) and groups with fast dynamics and response (e.g., zooplankton). We calculated several indicators to assess effects of harvesting and the ecosystem states and changes during the time period 1950-2013 ( Table 2).
The trophic levels of catches and ecosystem biomass may be affected by fisheries and have been used as indicator for ecosystem changes, with both expected to decrease in response to sizeselective exploitation (Branch et al., 2010). Ecopath calculates trophic level (TL) of the FGs, catches and various indices based on TL. The TL j of each predator group j was calculated using the equation: Where DC ij is the proportion of prey i in the diet of predator j and TL i is the trophic level of group i. In Ecopath it is assumed that all the detritus groups have trophic level 1. Trophic transfer efficiency calculated for a given trophic level is the ratio between the sum of exports plus the flow that is transferred from one trophic level to the next and the throughput on the trophic level (Christensen and Walters, 2004).
The sum of all direct and indirect effects of a FG on the food web were quantified applying mixed trophic impacts (Heymans et al., 2014). The mixed trophic impact m ij of a group is the product of all net impacts for all possible pathways that link groups i and j (Christensen and Pauly, 1992;Libralato et al., 2006). The total impact of each ecological group e i is calculated as A modified variant of the Kempton diversity index (Kempton Q) has been developed and implemented in EwE to measure the effects of fishing or climate on species in whole ecosystem models. Kempton Q express biomass diversity of groups with TL ≥3 and was expected to decrease with ecosystem degradation (Ainsworth and Pitcher, 2006;Steenbeek et al., 2018).
Production of each functional group over the modeled time period was calculated and total ecosystem production/biomass ratio and P/B-ratio of non-primary producers FGs and of harvested FGs were calculated. These P/B-ratios were expected to increase if high exploitation intensity decreases the proportion of long-lived exploited FGs biomass to total ecosystem biomass. The fishing mortality rate (F, year −1 ) for each exploited FG on a biomass basis (F = Y/B) was calculated as the ratio of annual catch yield (Y, g C m −2 year −1 ) and biomass (g C m −2 ). Fishing mortality is strongly positively related to fishing effort. The ratio of annual catch yield to annual production (Y/P = catch yield/production) will be used an indicator for intensity of exploitation of exploited groups (Mertz and Myers, 1998). The optimal FG specific exploitation rate (Y/P) that correspond to maximum sustainable yield from single stock considerations have been assessed to be about or slightly below 0.5, i.e., approximately equal fishing and natural mortality rate (Patterson, 1992;Zhou et al., 2012).

Model Parametrization, Evaluation, and Fitting of Ecosim Models
Initially in the balancing procedure, production of pelagic fish prey was less than consumption (i.e., EE > 1.0) in the year 2000 and 1950 models, and the biomass values for capelin and polar cod had to be increased relative to initial values (Supplementary Appendices 2, 4 Part F) to balance the models. In the balancing of the 1950-model, biomass values and total mortality rates for the small herring, capelin and polar cod were increased relative to the initial values to balance the need for prey (Supplementary Appendix 4 Part F). Most FGs except for the mammal and bird groups in the balanced year 2000 and 1950-models had relatively high EE's indicating that most of the production from most groups were consumed by groups within the model.
All models fitted to time-series for the 1950-1996 time-period without estimated vulnerabilities (M1, M3, M5, M7, and M9) had higher AICc and poorer fit than the corresponding models with estimated vulnerabilities (M2, M4, M6, M8, and M10) ( Table 3). Among the former models, the baseline model M2 fitted without fisheries data had much higher AICc than the models (M4, M6, M8 and M10) fitted to fisheries data and with estimated vulnerabilities ( Table 3). The two models (M8 and M10) with estimated vulnerabilities and forced by the capelin (0-2) mortality proxy had lower AICc-values than models (M4 and M6) fitted without the mortality proxy (  Table 2) and many estimated vulnerabilities were from trophic interactions with lower trophic level groups without time-series for biomass or catch. Further results presented were based on the M10 model since it had the lowest prediction SS (Table 3).
For the model M10 forced by PPR-proxy and small-herring induced mortality on capelin larvae, modeled biomass timeseries for most high trophic level (TL > 3) groups [minke whales, harp seals, Northeast Arctic cod (3+) and saithe (3+)] corresponded well with the observed time-series (Figure 2). For haddock (3+), the modeled biomass was lower than the observed biomass after ca. 2005 (Figure 2). Modeled (M10) and observed time-series for the boreal fisheries-exploited FGs Northeast Arctic cod (3+), minke whale, large redfish, large Greenland halibut and saithe (3+), had relatively high (Spearman r s > 0.46) positive correlations with modeled values for the period 1950-2013. In contrast, harp seal and pelagic amphipods, had negative correlations (Figure 3). Capelin groups and polar cod (2+) had moderate (r s from 0.40 to 0.65) positive correlations with modeled values, while polar cod (0-1) and haddock (3+) had no or very low (r s from 0.0 to 0.06) correlation to modeled values. The ratio of standard deviations showed that the observed time-series of haddock (3+), Scyphozoa, pelagic amphipods and Thysanoessa had high temporal variability and were not highly correlated with the modeled values. Observed time-series of long rough dab and northern shrimp (not shown in Figure 3) also had relatively high temporal variability.
In the models (M1-M6) lacking mortality forcing on capelin (0-2), the capelin groups in the Ecosim model did not follow the large observed changes from the early 1980s onwards with a large decrease in biomass 1985 and later periodical ups and downs (Figure 2). For the polar cod groups, the modeled biomasses were  (1950-1996) and prediction (1997-2013) for alternative Ecosim models.

Model
Name Fitting 1950Prediction 1997  Models were fitted to time series for the period 1950-1996 for biomasses, fishing mortalities, fishing effort, and catches and estimating vulnerabilities (V) by the step-wise fitting routine. Time-series for environmental drivers include a series with constant phytoplankton primary production (constant PPR), a series with phytoplankton primary production driven by the ice-cover and open water area (PPR-proxy), and a capelin (0-2) mortality proxy (CapM-proxy). Sum of squares (SS) for the period 1997-2013 was calculated to test the prediction ability of the models fitted for the earlier period. Sum of squares shown in bold for baseline models are calculated for only non-fisheries data and are not comparable to SS for model M3-M10 calculated for all data. SS, sum of squares calculated for model output and time-series observations; AICc, Akaike Information criterion.
Frontiers in Marine Science | www.frontiersin.org larger than the observed and the peaks in observed polar cod group biomasses around 2006 was not reproduced by the model which predicted increases in biomass (Supplementary Figure 2).
For northern shrimp, model M10 did not reproduce the peak in observed biomass around 1980, but both model and observed biomass had similar increasing trends after 1990 (Supplementary  Figure 2). Biomasses of other pelagic lower trophic level groups, such as Thysanoessa and medium-sized copepods, had more complex temporal variability. Both the modeled Thysanoessa and the observed krill time-series showed an increasing trend during the period 1990-2013, but the simulated biomass timeseries did not track the relative large year-to-year changes in the observed krill biomass indices (Supplementary Figure 2). The Russian and Ecosystem survey time-series for krill biomass, were moderately positively correlated for the time period  of overlapping measurements (Spearman r s = 0.39, P = 0.05).
The modeled biomass for medium-sized copepods and large calanoids had increasing trends after ca. 1995 contrasting the relative stable biomass in the observed mesozooplankton biomass time-series (Supplementary Figure 2). Modeled biomass trends during 1950-2015 for many lower trophic level (TL < 3) groups, i.e., detritivorous polychaetes and large bivalves, showed a similar U-shaped trend as the PPR-proxy with a slight dip in the cold period from 1960 to 1980 (Supplementary Figure 2). The long-lived groups, such as large bivalves and large epibenthic suspension feeders had smoother biomass trajectories and showed a more pronounced U-shape than groups with higher P/B and shorter lifespan, such as detritivorous polychaetes.
The comparison of models with (M11) and without (M10) forcing of ice-algae biomass and production showed that sympagic amphipods were strongly negatively affected by the reduction in ice-coverage after year 2000 (Supplementary Appendix 4 Part G). There were much smaller effects on other groups that fed on ice-algae or ice-algae detritus, but noticeable positive effects of high ice-algae production in the cold 1960-1980 period were found for biomasses of ringed and bearded seals, little auk, and Brünnich's guillemot. Effects of variable ice-algae production on polar-cod, pelagic amphipods and harp seals were small.
The increase in red king crab and snow crab in the crabbiomass-forced model M12 affected relatively few groups in the comparison to model M13 without crab-biomass forcing (Supplementary Appendix 4 Part E). The magnitude and direction of the effects were closely related to the importance of snow and red king crab in the diet of predators, and the importance of prey groups in the diet of the crab groups. Increasing snow crab and red king crab biomass in model M13 led to a positive effect on predator biomass (e.g., Northeast Arctic cod 3+) and negative effects for crab prey (Supplementary Appendix 4 Part E).

Food Web Structure and Major Flow Pathways
Trophic levels in the BS ecosystem ranged from 1 for primary producers to 5.1 for Polar bear in the year 2000 and 1950 models (Supplementary Appendix 4 Part F). Total biomass, production, consumption and total system throughput were slightly (0.1-11%) lower in the 2000 than in the 1950 model ( Table 4).
In the year 2000 model, the total ecosystem biomass (13.7 g C m −2 ) was mainly comprised of biomass from detritivorous benthic invertebrates (5.3 g C m −2 ), phytoplankton (2.0 g C m −2 ), other herbivorous zooplankton (1.5 g C m −2 ) and krill (1.1 g C m −2 ) (Supplementary Table 4). Atlantic cod, the main fishery target, had a biomass of 0.10 g C m −2 .
The total ecosystem production was 167 g C m −2 yr −1 with major contributions from the primary producers: phytoplankton (110 g C m −2 yr −1 ) and ice algae (5.3 g C m −2 yr −1 ). At trophic levels between 2 and 3, the main producers were the aggregated compartments microzooplankton and HNAN (19.1 g C m −2 yr −1 ), bacteria (16.1 g C m −2 yr −1 ), other herbivorous zooplankton (8.2 g C m −2 yr −1 ), krill (2.7 g C m −2 yr −1 ) and detritivorous benthic invertebrates (3.1 g C m −2 yr −1 ). At higher trophic levels (TL > 3), major producers were capelin (0.45 g C m −2 yr −1 ), other planktivorous fishes (0.44 g C m −2 yr −1 ) and shrimps (0.17 g C m −2 yr −1 ). Other demersal and benthic fish had a production of 0.17 g C m −2 yr −1 ) and cod had a production of 0.08 g C m −2 yr −1 . Polar bear, whales, seals, cod, other demersal and benthic fishes, capelin and the zooplankton groups had somewhat higher biomass, production and consumption in the 1950 that the 2000 model (Supplementary Figure 3). Four major pathways for carbon flow from lower to higher trophic levels were evident for the year 2000 model (Figure 4); the microbial food web pathway, the copepod pathway, the krill pathway and the benthic invertebrate pathway. With regard to the importance as prey, the krill compartment comprised of the FGs Thysanossa and large krill had the most (n = 8) major prey flows (i.e., among the three largest flows to a predator compartment from prey compartments) (Figure 4). Pelagic planktivorous fishes, herbivore zooplankton, and detritus had the second most connections with predator compartments with five major flows. Whereas, krill had major flows to five top-predator compartments, the herbivorous zooplankton compartments had only one major flows to a top-predator compartment (Figure 4).
The five FGs with highest total trophic impact (see Eq. 7) in the year 2000 model were; (1) diatoms, (2) polar cod (2+), (3) Thysanoessa, (4) medium sized copepods and (5) small benthic molluscs (Figure 5). These FGs had contrasting trophic impacts depending on their trophic position (Figure 5). Diatoms had a positive impact as food source for many lower trophic level FGs and had a much larger impact than autotrophic flagellates, the other phytoplankton group (Figure 5). Medium sized copepods and large calanoids had a positive impact as prey for planktivorous fishes and pelagic predatory groups (chaetognaths, cephalopods, Ctenophora, scyphomedusa, and northern shrimp). The krill groups also had a strong positive impact as prey for demersal fishes, some bird groups and several whale and seal FGs ( Figure 5). However, the krill groups also had negative impact on some other planktonic invertebrate groups. Capelin and polar cod had positive impacts as prey for demersal fish FGs, seals and some whale groups, and negative impacts as predators on krill, the copepod groups and other pelagic zooplankton FGs (Figure 5). Northeast Arctic cod (3+) had negative impacts as FIGURE 4 | Carbon flows between aggregated major compartments based on the Ecopath model for year 2000 with four major pathways for carbon flows from lower to higher trophic levels. (i) the microbial food-web pathway (violet lines), (ii) the copepod pathway (yellow lines), (iii) the krill pathway (red lines), and (iv) the benthic invertebrate pathway (brown lines). Functional group outputs are aggregated according to Table 1. Thick lines shows major flows, i.e., among the three largest flows to each aggregated compartment. Thin lines shows smaller flows. "H" in circles indicate harvested compartments. predator on several other fish FGs and positive impact as prey for Greenland shark and dolphins (Figure 5). Northeast Arctic cod also had negative effects on seal and fish FGs that shared prey with cod ( Figure 5).
Trophic levels estimated by our Ecopath model for year 2000 differed considerably from previously published TL-values from mass-balance models for the BS (Supplementary Table 3). TL-values from our model were on average higher (0.14 and 0.25) than the trophic levels from Dommasnes et al. (2002) and Berdnikov et al. (2019) and lower (-0.21 and -0.25) than the TL's from Blanchard et al. (2002) and Bentley et al. (2017).

Temporal Variation in Ecosystem Properties and Effects of Exploitation and Climate Variability
The total catch in the BS peaked in the late 1970s mainly driven by the large catches of capelin (Supplementary Figure 5). The  (Supplementary Figures 5, 6).
The trend in the PPR-proxy showed an U-shaped trend with low values in the 1960-1980s (Supplementary Figure 1). Similar U-shaped biomass trends were evident for mammal and birds, total fish biomass, demersal fishes, pelagic invertebrate, and benthic invertebrates group (Supplementary Figure 4). Total biomass of the ecosystem decreased from 1950 to the lowest values around 1970 and thereafter increased toward 2013 (Supplementary Figure 4). In contrast to biomass trends of most other groups, pelagic fish biomass was highest in the period 1970-1980 largely driven by the high capelin biomass at this time.
The Kempton's diversity index Q changed moderately during 1950-2015 but had lower values at the end of the time-period than in the beginning (Figure 6). Ecosystem P/B for non-primary producing FGs showed a modest increase in the 1970-1980s, but there was a clear peak in P/B for harvested FGs in the 1970-1980s. Trophic level of the catch (TLc) decreased from 1950 to 1985 followed by three periods of ups-and downs corresponding to periods of opening and closure of the capelin fishery (Figure 6).
Patterns over time for catches and fishing mortalities (F = Y/B) varied substantially between FGs and for Polar bear, minkewhale, harp seals and small herring, catches and catch as proportion of production (Y/P) showed decreasing trends after 1950 (Figure 7 and Supplementary Figures 6, 7). For large Greenland halibut, large redfish, polar cod (2+), northern shrimp and capelin (3+), catches and (Y/P) rose rapidly during the 1970s followed by decreasing trends. Catches of the large gadoids; the Northeast Arctic and coastal cod groups, haddock and saithe, showed temporal variability with low catches in the 1980's. Y/P for these FGs peaked in the period 1970-1980 and were relatively low after 1980 (Figure 7). Y/P were above 0.5 in periods for all exploited FGs except for Northern shrimp and Polar cod (2+) that had low fishing mortalities and low Y/P's (Figure 7 and  Supplementary Figure 7).

Main Findings
Intense harvesting had clear negative effects on the biomasses of exploited and long-lived high trophic level FGs, such as mammals, large gadoids, Greenland halibut and redfish. Decrease in fishing pressure the later years contributed to increases of biomasses for many higher trophic level FGs, indicating a recovery period. An increase in primary production after ca. 1990 due to reduced ice-coverage and larger open water area had a positive impact on production and biomass of boreal FGs at all trophic levels.

Model Evaluation and Fitting to Time-Series
The Ecopath models initial values for biomass and production of capelin and polar cod in year 2000 and 1950 did not match the consumption demand from predators and biomasses had to be increased to match the consumption demand. This could be due to underestimation of capelin and polar cod stocks as shown in previous studies Hop and Gjøsaeter, 2013). An alternative explanation for the apparent productionconsumption mismatch may be that the consumption on planktivorous fish has been overestimated due to bias in the diet composition of predators. Diet data for most predators except for Northeast Arctic cod were not adjusted for possible preygroup-specific differences in digestion rates. Since large fish prey are more slowly digested than smaller prey (Salvanes et al., 1995), predator diet proportions of pelagic fishes may have been overestimated relative to smaller invertebrate prey.
The cross-validation procedure revealed that the model (M8) with fisheries and capelin mortality proxy and constant primary production had a better fit and lower AICc for the fitting period but higher sum of squares for the prediction period than for the model M10 with ice-coverage forced primary production. This led to the conclusion that M10 was the model with most support and this model showed an effect of increasing PPR especially in the last part of the period 1950-2013.
The moderate effects of including snow crab and red king crab FGs in the year 2000 model run with forced biomasses of red king crab and snow crab from year 2000 to 2013 suggest that these FGs were unlikely to have a major effect at the whole ecosystem level during the studied time-period. For red king crab which has a coastal distribution in the southern part of the BS, strong effects on local and regional scale on the bottom fauna have been described (Oug et al., 2011), and similar effects may be expected following ongoing the snow crab expansion.

Comparison With Other Studies
There was a good correspondence between trophic levels estimated for 83 FGs in the Ecopath model for year 2000 and for independent data on trophic levels and δ 15 N from stable isotope data from the BS (Pedersen in prep.). The differences between trophic levels estimated by our and other mass-balance models for the BS may be due to between-model differences in group structure and diets of dominant FGs.
Values for biomass, production and consumption of seals, krill, mesozooplankton, and bacteria differ substantially between our model and those published earlier for the BS (Supplementary Table 4). The twice as high food consumption from seals in our model compared to values given by Sakshaug et al. (1994), was mainly due to lower Q/B-values used by Sakshaug et al. (1994). Further, production from krill in our year 2000 Ecopath model was about twice that given by Sakshaug et al. (1994). For both krill groups, we used a higher P/B than Sakshaug et al. (1994) (2.5 vs. 1.5 yr −1 ). Production by bacteria in our model was less than 25% of the value given by Sakshaug et al. (1994) and both Sakshaug et al. (1994) and Berdnikov et al. (2019) used much higher P/B-values for bacteria  year −1 ) than in our model (21.2 year −1 ). The P/B values for the mesozooplankton, benthic invertebrate and meiofauna FGs in our model were much lower than the values from Blanchard et al. (2002; Supplementary Table 4). Our Ecopath model values for P/B for other herbivore zooplankton including copepods were somewhat higher than estimates derived from other models for Calanus. That our year 2000 model had a lower ecosystem P/B than the 1997 model by Blanchard et al. (2002) (P/B = 12.2 vs. 15.9 year −1 ) is likely caused by the lower P/B-values for several low trophic level FGs in our model (Supplementary Table 4). Differences in ecosystem P/B-values are likely to affect how fast the Ecosim models react to perturbations and such effects should be further investigated.

Ecosystem Structure and Major Carbon Flow Pathways and Compartments
The four dominant carbon flows pathways (microbial food web, copepod, krill, and benthic invertebrate pathways) differed regarding P/B of their contributing FGs and their functions as prey sources. The carbon flows within the microbial food web with high P/B and high turnover merges with the krill and copepod pathway and "hitch-hikes" to higher trophic levels. Microzooplankton made up large proportions (10-32%) of the diets of medium zooplankton, large calanoids, small copepods and Thysanoessa in the model, which is consistent with previous studies suggesting a substantial flow from the microbial food web to higher trophic level pelagic FGs in the BS (Hansen et al., 1996;De Laender et al., 2010).
The importance of the copepod pathway was indicated by the high total impact ranks of medium sized copepods and large calanoid copepods, mainly resulting from impacts as prey for planktivorous fish FGs and pelagic carnivorous invertebrate FGs, but also impacts as consumers of microzooplankton and phytoplankton. Copepods are also major prey for early fish stages. In contrast to the krill pathway, the copepod pathway had less direct impact as prey for birds and mammals except as prey for little auks and bowhead whales (Supplementary Appendix 4 Part C).
The krill pathway contributed to an energy-efficient carbon flow to higher trophic levels and top-predators as seen in other ecosystems (Murphy et al., 2007;Ruzicka et al., 2012). The fact that the krill groups were important as prey, but also as a predator and potential competitor for other zooplankton FGs, suggest that krill may have a wasp-waist function in the BS food web. The increase in biomass of krill and corresponding increase in proportion of krill in the diet of age 1-2 year cod after 1984 , suggest increased importance of krill during the last part of the time period. Advection of krill from the Norwegian Sea to the BS was found to be more prominent in warm years (Orlova et al., 2015), and is likely to have contributed to the increased importance of krill in the BS during the warming period.
The detritus-based benthic invertebrate pathway transports carbon to predatory benthos FGs, demersal and benthic fish and also some birds and mammals FGs. Field-based production estimates for macrobenthos for the BS are scarce and uncertain, but local estimates of production range from 0.1 to 20 g C m −2 year −1 (Kędra et al., 2013(Kędra et al., , 2017. Despite the uncertainty, P/B's (<1.0 year −1 ) of the macrobenthos FGs were much lower than for krill and copepods (P/B of c. 2.5-6 year −1 ) at similar trophic level. Thus, compared to the pelagic pathways, the benthic invertebrate pathway is a "slow" energy channel with low turnover but high biomass (c. 5 g C m −2 ). The presence of pathways or "channels" with different turnover rates and P/Bvalues may enhance stability in ecosystems (Rooney et al., 2006). Our study emphasize the multi-pathway structure in the BS with a fast partly detritus based microbial food web, a slower pelagic krill pathway in addition to the copepod and benthic invertebrate pathway.
Our model confirms the importance of pelagic planktivorous fish, such as capelin as prey and cod as top-predator emphasized in previous studies . However, in our model, the aggregated compartments "other demersal and benthic fishes" comprising 18 FGs had a total biomass, production and food consumption that was about the twice that of cod in the year 2000 model (Supplementary Table 4). That the consumption of fish by this aggregated compartment was similar to cod indicates a potential for top-down effects from FGs in this compartment and potential for significant competition with cod (Supplementary Table 4).

Effects of Harvesting
Fishing had a major impact on the BS ecosystem during the study period 1950-2013 as indicated by model M10 with predicted biomasses that corresponded well to the observed biomasses of most of the boreal and historically most exploited FGs, including during the latest time period (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Optimal exploitation rates (Y/P) depend on life history characteristics, but for fish stocks, optimal Y/P are suggested to be equal to or slightly below 0.5, i.e., fishing and natural mortality being similar (Patterson, 1992;Zhou et al., 2012). The observation that most FGs targeted by the fishery in the BS showed periods when Y/P were larger than 0.5 indicates overexploitation. This is in general accordance with results from single stock assessments (Gjøsaeter, 1998;Nakken, 1998;Toresen and Østvedt, 2000). For several fish stocks, especially the long-lived demersal stocks of Greenland halibut, redfish and Northeast Arctic cod, the increases in fishing mortality in the BS from the 1950s to 1970-80s were evident both in the Ecosim-modeled and observed biomass time-series derived from single-stock assessments (Bowering and Nedreaas, 2000;Johannesen et al., 2012). The increase in ecosystem P/B of harvested FGs in the period of highest fishing mortality support that exploitation had a notable effect on ecosystem structure.
Fishing effort and catches increased rapidly after 1945 (Nakken, 1998), and the Norwegian spring-spawning herring was the first fish stock to collapse due to overfishing in the 1960's (Toresen and Østvedt, 2000). The adult part of this stock is mainly distributed and harvested in the Norwegian Sea, and after the collapse, fishery on juvenile herring which has its major nursery area within the BS was closed in 1971. Following the stock collapse of herring, the capelin fishery expanded in the late 1960s and the 1970s. From the mid-1980s, the BS capelin experienced a series of stock collapses and the variability in total fish catches in the BS and our results show that the trophic level of catches was mainly driven by the state of the capelin fishery. After 1991, the capelin fishery was only open in years when the expected spawning stock was higher than a level estimated by probabilistic assessment accounting for predation from cod on maturing capelin in the pre-spawning period . Northeast Arctic cod has been the major fishery target in the BS in terms of fishing effort and commercial value, but is also a major predator on capelin and other smaller fish and invertebrates . After a period of increasing fishing mortality after 1950, the decrease in fishing mortality of this large stock after ca. 1990 contributed strongly to the recovery and increase in stock biomass, production and catches (Nakken et al., 1996).
The populations of the large baleen whales (bowhead, blue, and fin whales) and walruses were heavily exploited and reduced to levels far below pristine levels prior to 1950 (Reeves, 1980;Weslawski et al., 2000) and Greenland shark was also exploited before and to a low extent after 1950 (Supplementary Appendix 2). For minke whales and harp seals, catches and harvesting mortalities were reduced to lower levels during the period 1950-2013. Biomass trends for these exploited mammal FGs were U-shaped and showed trends of recovery as a result of reduced exploitation and increased ecosystem productivity.
The sequence in which fisheries evolved in the BS resembles a pattern consistent with the "Fishing through the food web" concept (Branch et al., 2010). However, the sequence in which the fisheries targeted various FGs may have been directed by availability rather than trophic level. Fishing pressure on pelagic intermediate trophic level fishes, such as herring and capelin were increasing relatively early in the study period, simultaneously with pressure on Northeast Arctic cod. This may have contributed to a kind of trophic balance when both typical prey and predator FGs were exploited simultaneously. Nilsen et al. (2020) explored balanced harvesting strategies for the Norwegian and BS using the Atlantis model and found that a balanced harvesting regime would only produce marginal increases in total yield of currently exploited FGs compared to the historical exploitation regime after 1980.
Most harvested fish stocks responded with an increase in biomass when fishing mortality (Y/B) decreased in the 1990's (Nakken, 1998) indicating a recovery period. The patterns of catches, biomasses and fishing mortality indicate that the intense exploitation was relaxed around 1990-2000 for many of the targeted fish and exploited mammal FGs and contributed to increases in their biomasses. The Golden Redfish (S. norvegicus) stock and the Coastal cod stock have not fully recovered (ICES, 2019), and may be exceptions to the recovery pattern described above.
Recent studies have shown that the dynamics of cod is tightly linked to the harvesting rate and capelin abundance (Lindstrøm et al., 2009;Koen-Alonso et al., 2021). If their top-down control was dominating in the BS ecosystem, prey groups would expected to be "released" from predation and increase during periods of heavy exploitation i.e., from ca. 1970 to 1990 when some predator stocks, e.g., Northeast Arctic cod and Greenland halibut, were reduced to low levels. One would also expect that the biomass of competitors would increase simply due to reduced competition. The question is if there is any evidence of trophic-cascade and competitive-release effects in the model? The predicted increase in biomass of capelin and long rough dab during this period by the M10 Ecosim model suggest cascade effects due to reduced predation pressure. The period of high biomass of capelin in the period 1970-1985 was likely a result of reduced predation from small herring predation on capelin larvae during the period of collapse of the Norwegian spring spawning herring. In addition, low levels of predation from demersal fishes and cod on older capelin may have contributed to the high capelin biomass level. The increase in modeled biomass of long rough dab during 1970-1990 cannot be verified by survey observations since the stock-assessment time-series started in 1989 (Supplementary Appendix 4 Part A,B).
Among the other FGs, northern shrimp had the longest observational time-series, starting in 1971, and the very low fishing mortality compared to the predation mortality suggests that fishery exploitation had not been a dominant driver for this stock. Decreases in cod stocks in the Northwestern Atlantic have been accompanied by sharp increases in biomass of invertebrate stocks, and this has been interpreted as results of a trophic cascade (Frank et al., 2005). In the BS, however, several shrimppredating FGs may have increased when the cod stock decreased contributing to a relative stable predation pressure and modest biomass changes for northern shrimp. Thus, we suggest that except for capelin, there is little evidence for a strong cascading effect of medium and low trophic level FGs during the period of low top predator abundance in the BS.
The Ecosim-simulated biomass of detritivorous and predatory macrobenthos showed a U-shape during 1960-1990 that could be attributed to the relative low water temperature, small openwater area (=extensive ice-cover) and low primary production. There was no available time-series of macrobenthos biomass at the ecosystem level, but estimates based on grab sampling over a large part of the BS have changed over time from high biomass in 1924-32, to low biomass in 1968-72 and then high biomass in 2003 again (Denisenko, 2001;Jørgensen et al., 2017). The changes in biomass among time-periods have been attributed to changes in both climate and primary production and to changes in bottom trawling effort hypothesized to affect macrobenthic biomass negatively (Denisenko, 2001). To simulate and evaluate the possible direct effects of bottom trawling in Ecosim, a relationship between bottom trawling effort and mortality for various FGs of benthic invertebrates have to be included in future simulations.
The moderate changes in the Kempton diversity index, which is expected to react to intense exploitation during the period 1950-2013, suggest modest changes in ecosystem structure. So why and how has the ecosystem resisted the heavy exploitation of some FGs? The high fishing effort during the period 1965-1990 in the cod fishery in the BS have likely also increased fishing mortality of other fish FGs which were caught as bycatch (Denisenko, 2001;Rusyaev and Orlov, 2013). In the BS, there may have been few species that could take over for cod as major piscivore during the period of intense exploitation, and candidates, such as Greenland halibut, redfish and harp seals had been extensively targeted and reduced by harvesting prior and during the cold period.

Effects of Climate Variability
The U-shaped trend in biomass of many FGs is interpreted as mainly a result of trends in primary production. Models forced by the PPR-proxy performed better during the warm period (after mid 1990s) suggesting that more open water is linked to increased PPR and food web productivity. Earlier studies have emphasized the positive effects of relatively warm Atlantic water in the southern part of the BS on fish stock recruitment and individual fish growth (Sundby, 2000;Ottersen et al., 2002), especially for Northeast Arctic cod, haddock and Norwegian spring spawning herring partly spawning in the Norwegian Sea (Bogstad et al., 2013). High variability in haddock recruitment lead to relatively low predictability by Ecosim for haddock (3+) biomass. The relatively low temporal model fit to the indices (relative biomasses) for the youngest stanza of the cod and haddock groups suggest that recruitment may be driven partly by other mechanisms than represented by the PPRproxy, and further testing of other drivers for recruitment may improve model fits and predictability. Variability in advection of Atlantic water and copepods, krill and young fish stages into the southwestern part of the BS affect the ecosystem (Drinkwater, 2011), and is likely to contribute to variability that was not captured by the Ecosim models.
The correlations between modeled and observed data for the unharvested lower trophic-level FGs were lower than for the most heavily exploited FGs and may be due to a lower signal to noise ratio for observed data for these FGs. The lower trophic level FGs also had higher temporal variability than for the more long-lived higher trophic level FGs. This indicates that the Ecosim model did not fully reproduce short-term variability in lower trophic level FGs but could still capture the main long-term trends.
The proxy for PPR that was applied in the calibration and fitting to time-series in Ecosim was based on a well-founded relationship between primary production and open-water area (Dalpadado et al., 2020). This relationship is supported by model studies showing lower primary production at lower temperature, and temperature and open water area were strongly positively correlated (Wassmann et al., 2006b;Slagstad et al., 2011). The improvement in model prediction by including our PPR-proxy as an environmental driver in the model suggests that changes in PPR driven by changes in open-water area and indirectly by water temperature had an effect on the development of the BS ecosystem during the 1950-2013 period.
The lower trophic level FGs for which we had observed biomass time-series (krill and mezozooplankton, pelagic amphipods, and scyphomedusae) showed contrasting trends after year 1995. The observed biomasses of the two krill groups and scyphomedusae showed increasing trends from year 2000 to 2013 while there was a stable biomass of mesozooplankton (Dalpadado et al., 2020) and a decreasing trend for pelagic amphipod biomass. The Ecosim model, however, predicted an increase in biomass of medium sized copepod and large calanoid biomasses after 1995. Mesozooplankton biomass is dominated by the mainly boreal medium sized copepods and the arctic large calanoids, and these FGs may have responded differently to warming in the western part of the BS after 1995 (Aarflot et al., 2017). Biophysical modeling with warming scenarios show expectations of increased production the boreal C. finmarchicus and decreased production in the arctic C. glacialis (Slagstad et al., 2011). Predation from capelin has been emphasized to have a major top-down effect on mesozooplankton in the BS, but after a peak in observed biomass mesozooplankton around 1994 when capelin biomass was low, mesozooplankton biomass has been stable despite ups and downs in capelin biomass (Dalpadado et al., 2020). Stige et al. (2019) suggested that less sea-ice coverage may have a negative effect on the arctic large calanoid C. glacialis.
The krill biomass in the BS was dominated by Thysanoessa, and krill had the longest observed time-series among lower trophic level FGs. Increases in both observed and modeled krill biomass in the period after ca. 1995 indicates that the energyefficient krill pathway may have strengthened during the period 2000-2013. Krill as prey may have contributed to shorten the food chains and enhance production at high trophic levels. The temporal year-to-year variability in the observed krill time-series was not well-reproduced by the model. That the observational time-series for the two krill groups were moderately positively correlated may indicate that they both represent temporal trend in the krill biomass but with a relatively low signal to noise ratio. It is challenging to estimate biomass of krill precisely due to very patchy spatial distribution  and varying advection of krill into the BS may also contribute to variability (Orlova et al., 2015).
The modeled effects of a decreasing trend in sea-ice coverage and reduced ice-algae production (model M11) after ca. 1980 notably affected biomasses of ringed and bearded seals, little auks, and Brünnich's guillemots. Polar cod and pelagic amphipods were less affected and variable ice-algae production could not explain the observed decreasing time-trends in these predator FGs after ca. year 2000. This may indicate that sea-ice may be a limiting habitat for these FGs beyond the production of ice-algae. Seaice coverage is important for polar cod during reproduction and recruitment and both large calanoid copepods and pelagic amphipods feed on ice algae and these FGs are important prey for polar cod (Hop and Gjøsaeter, 2013;Bouchard and Fortier, 2020; Supplementary Appendix 2). More knowledge on the dependence of the ice-habitat habitat beyond the effect of icealgae production and other effects of warming, may be needed to improve model input and performance.
The Ecopath model output suggests that scyphomedusae did not have a major predatory effect in the BS ecosystem despite its increase in the warm period after year 2000 (Eriksen, 2016). For Ctenophora, there was no time-series or precise biomass estimate, but recordings of frequency of occurrence of Ctenophora in Northeast Arctic cod stomachs shows a clear increase after 1996 in the southwestern BS (Eriksen et al., 2018), and may suggest an increase in biomass of Ctenophora in the area during this period.
The inclusion of mortality from small herring on capelin larvae in the Ecosim model increased the model fit to observed data by primarily improving the fit for the capelin groups but not for the other FGs. This may suggest that topdown and bottom-up effects of capelin in the Ecosim model were moderate during the modeled time-period. An apparent top-down effect from capelin as predator on krill has been observed (Eriksen and Dalpadado, 2011), and field measurements revealed that biomasses of capelin and total mesozooplankton varied inversely during 1989-1997 but not in the period after 1997 (Dalpadado et al., 2020). Strong negative effects of low capelin biomass on predators, such as Northeast Arctic cod and harp seals were observed during the first capelin collapse in 1985-1988(Gjøsaeter et al., 2009), but effects were lower during later collapses, likely due to larger abundance of alternative prey (Gjøsaeter et al., 2009). This inconsistency in correlations suggests complex trophic interactions and potential indirect effects that are difficult to identify from modeled or observed time-series.
The patterns of mixed trophic impacts for various Ecopath model FGs showed that most FGs had both bottom-up and top-down impacts, suggesting that both types of trophic control have been important in the BS and other studies also point in this direction (Johannesen et al., 2012;Lindstrøm et al., 2017;Stige et al., 2019). By examining predator-prey correlations from the BS, Johannesen et al. (2012) found shifts between negative and positive correlations during the time period 1977-2002, indicating shifts in trophic control between bottom-up and topdown dominance. Stige et al. (2019) also noted that both bottomup and top-down effects were present when considering pelagic fish and zooplankton interactions.
Lower fishing mortalities coinciding with warming and increasing primary production during the recovery period after around 1990 may have strengthened the role of cod and other demersal fishes as top predators. The coincidence of the period of overexploitation of fish stocks with the cold low-productive climatic period during 1960-1980 may have prevented other species to take over when the stocks of large gadoids and the long-lived redfish and Greenland halibut had been intensively exploited and reduced. The relatively low diversity of the nonexploited fish FGs in the BS may also have contributed to the lack of success of other species to replace exploited stocks. How ecosystem management can be used to preserve structure and mitigate negative climatic effects should be investigated in future studies.

CONCLUSION
Four major carbon pathways were identified in the BS and the modeling results indicated increased productivity at lower trophic levels during warm years with large ice-free open-water area after ca. 2000. This contributed to higher productivity for most high trophic level FGs. The krill pathway was important for both medium and high trophic level compartments, and krill biomass and production increased during the warm period. There were signs of decrease in observed biomasses of some high arctic FGs that were not reproduced by the models even after forcing the model ice-algae with ice coverage time-series.
In the low-productive period from 1960 to 1985, fishery exploitation reduced biomasses of FGs in a sequential pattern causing reductions of biomasses for mammals, large gadoids and other long-lived demersal fishes. The increased biomass for capelin during this period was interpreted as a trophic cascade effect of relaxed predation. When exploitation was relaxed, biomasses of many exploited FGs increased during the recovery period after about 1990. Despite heavy exploitation, the basic ecosystem structure seems to have been preserved in the BS during the periods of overexploitation and recovery.

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.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because this was a modeling study using historical data and as such we did not need ethical review and approval.

AUTHOR CONTRIBUTIONS
TP developed the initial model concept and structure and implemented the models. NM contributed to the data structuring. TP, NM, UL, and PR contributed to the ideas in this manuscript and drafted the manuscript. All authors contributed to the data input, provided edits to the manuscript, contributed to the article, and approved the submitted version.

FUNDING
This Norway-Canada collaborative project "A Transatlantic innovation area for sustainable development in the Arctic" (CoArc) provided partial support for this work (Project Number 8048 at Akvaplan-Niva and QZ-15/0457 at Norwegian Ministry of Foreign Affairs). The Research Council of Norway also provided partial support through the project "The Nansen Legacy" (RCN#276730).

ACKNOWLEDGMENTS
We would like to thank Bjarte Bogstad and Martin Biuw at Institute of Marine Research, Norway for providing cod stomach summary data and harp seal data, respectively. We thank Leif C. Stige at University of Oslo for providing time-series on pelagic amphipods.