The BS is a high latitude LME, covering an area of 2.01 million km² (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 (Fossheim et al., 2015).

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).

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 (Heymans et al., 2016). 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.

The Ecopath model for 1950 was based on an Ecopath model for year 2000 but with biomass, fisheries, P/B and Q/B-values specific for year 1950. Less information regarding many ecological groups was available for the period around 1950 and we chose year 2000 to represent a presumably similar year as 1950 with regard to temperature, and for balancing an annual average year 2000 Ecopath model. The average temperature in the Kola-section in 2000 was similar to 1950 (4.6°C vs. ca. 4.7°C) (Dippner and Ottersen, 2001; Supplementary Figure 1). The water temperature time-series from 1951 to 2013 (average from 0 to 200 m depth, st. 3–7) from the Kola section (70°30′N to 72°30′N along 33°30′E) [source: Knipovich Polar Research Institute of Marine Fisheries and Oceanography (PINRO)]¹ (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 1950-model 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 (Heymans et al., 2016). 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.

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.² To simplify presentation of model results, some output values were presented for aggregated compartments (Table 1).

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.

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 (1950–1996) 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 well-documented 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 time-series were used as input in the time-series fitting, including time-series 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 time-series 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 time-series 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 PPR-proxy 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 (Supplementary Appendix 4 Part A) for the period 1950–2013, and the results from model M11 were compared to model M10 without forcing of the ice-algae. To test if the invasive red king crab and the expanding snow crab may have affected the ecosystem, models with snow-crab and red king crab groups were run for the period 2000–2013 (Supplementary Appendix 4 Part E). The year 2000 model with 26 estimated vulnerabilities from model M10 was run with (model M12) and without (model M13) forcing by observed crab biomass time-series (Supplementary appendix 4 Part A), and the output biomasses from the Ecosim models were compared.

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 pedigree-specified 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).

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 size-selective 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).

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 3). Model M10 forced by PPR-proxy and with 26 estimated vulnerabilities had higher AICc than model M8 forced by constant PPR with 51 estimated vulnerabilities. However, model M10 with 26 vulnerability values had lower prediction SS (SS for the prediction period 1997–2013) than M8 which had far more (n = 51) estimated vulnerabilities. Time-series of model biomasses from exploited fish groups for models forced by the PPR-proxy had a more U-shaped trend during the period 1950–2013 than for models forced by constant PPR.

For M10, all trophic interactions with estimated vulnerabilities included at least one mammal or fish groups with time-series (Supplementary Table 1). Most (n = 18) of the 26 fitted vulnerabilities were low (vulnerability values < 2) indicating bottom-up effects, and the high vulnerabilities indicating top-down effects were estimated for interactions with top-predators; minke whale, harp seals, Northeast Arctic cod (3+), coastal cod (2+), saithe (3+), and long rough dab (Hippoglossoides platessoides) (Supplementary Table 1). For M8 with 51 estimated vulnerabilities, 29 vulnerabilities were >>2 (Supplementary 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 time-series 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 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 time-series 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 (1980–2005) 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 crab-biomass-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).

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 diet matrix of the year 2000 -model including the snow-crab and red king crab groups had a total of 1,029 feeding interactions and a connectance index of 0.095. The FGs with most groups (n) preying on them were; Thysanoessa (n = 43), pelagic amphipods (n = 35), small herring (n = 33), medium sized copepods (n = 31), capelin age 0–2 (n = 30), and northern shrimp (n = 29).

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 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).

The total catch in the BS peaked in the late 1970s mainly driven by the large catches of capelin (Supplementary Figure 5). The catches of mammals decreased during the 1950s and 1960s, stabilized during 1965–1990 and decreased to low levels after 2005 (Supplementary Figure 5). Catches of cod had a decreasing trend from the 1950s and reached a minimum in the period from 1980 to 1990 and then increased to 2013. Other demersal fishes had a peak in the 1970s due to large catches of redfish and Greenland halibut, and had an increase after year 2000 (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, minke-whale, 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).

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.

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 (Gjøsæter et al., 1998; Hop and Gjøsæter, 2013). An alternative explanation for the apparent production-consumption 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 prey-group-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.

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 δ¹⁵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 (125–200 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.

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 (Bogstad et al., 2015), 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, 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/B-values 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 (Bogstad et al., 2015). 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).

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–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øsæter, 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 (Gjøsæter et al., 2012). 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 (Bogstad et al., 2015). 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; Christensen et al., 1992; 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 shrimp-predating 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 open-water 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.

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 PPR-proxy, 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 energy-efficient 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 (Eriksen et al., 2016) 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. Sea-ice 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øsæter, 2013; Bouchard and Fortier, 2020; Supplementary Appendix 2). More knowledge on the dependence of the ice-habitat habitat beyond the effect of ice-algae 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 top-down 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øsæter et al., 2009), but effects were lower during later collapses, likely due to larger abundance of alternative prey (Gjøsæter 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 top-down dominance. Stige et al. (2019) also noted that both bottom-up 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 non-exploited 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.