# Responses of ecological indicators to fishing pressure under environmental change: exploring non-linearity and thresholds

Caihong Fu, Yi Xu, Arnaud Gruss, Alida Bundy, Lynne Shannon, Johanna J.

Heymans, Ghassen Halouani, Ekin Akoglu, Christopher P. Lynam, Marta

Coll, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Caihong Fu, Yi Xu, Arnaud Gruss, Alida Bundy, Lynne Shannon, et al.. Responses of ecological indicators to fishing pressure under environmental change: exploring non-linearity and thresholds. ICES Journal of Marine Science, 2019, 77 (4), pp.1516-1531. 10.1093/icesjms/fsz182. hal-02434138

HAL Id: hal-02434138

https://hal.science/hal-02434138

Submitted on 8 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
ICES Journal of Marine Science (2020), 77(4), 1516–1531. doi:10.1093/icesjms/fsz182

**Contribution to the Themed Section: 'Science in support of a nonlinear non-equilibrium world'**

**Original Article**

**Responses of ecological indicators to fishing pressure under environmental change: exploring non-linearity and thresholds**

Caihong Fu ID <fnref n="1" />*, Yi Xu<fnref n="2" />, Arnaud Grüss ID <fnref n="3" />, Alida Bundy<fnref n="4" />, Lynne Shannon<fnref n="5" />,<fnref n="6" />,
Johanna J. Heymans<fnref n="7" />,<fnref n="8" />, Ghassen Halouani<fnref n="9" />, Ekin Akoglu<fnref n="10" />, Christopher P. Lynam<fnref n="11" />,
Marta Coll<fnref n="12" />,<fnref n="13" />, Elizabeth A. Fulton<fnref n="14" />,<fnref n="15" />, Laure Velez<fnref n="16" />, and Yunne-Jai Shin<fnref n="16" />

Fu, C., Xu, Y., Grüss, A., Bundy, A., Shannon, L., Heymans, J. J., Halouani, G., Akoglu, E., Lynam, C. P., Coll, M., Fulton, E. A., Velez, L., and
Shin, Y.-J. Responses of ecological indicators to fishing pressure under environmental change: exploring non-linearity and thresholds. –
ICES Journal of Marine Science, 77: 1516–1531.

Received 14 May 2019; revised 22 July 2019; accepted 4 September 2019; advance access publication 27 September 2019.

Marine ecosystems are influenced by multiple stressors in both linear and non-linear ways. Using generalized additive models (GAMs) fitted to outputs from a multi-ecosystem, multi-model simulation experiment, we investigated 14 major ecological indicators across ten marine eco-systems about their responses to fishing pressure under: (i) three different fishing strategies (focusing on low-, high-, or all-trophic-level taxa); and (ii) four different scenarios of directional or random primary productivity change, a proxy for environmental change. From this work, we draw four major conclusions: (i) responses of indicators to fishing mortality in shapes, directions, and thresholds depend on the fishing strategies considered; (ii) most of the indicators demonstrate decreasing trends with increasing fishing mortality, with a few exceptions depending on the type of fishing strategy; (iii) most of the indicators respond to fishing mortality in a linear way, particularly for community and biomass-based indicators; and (iv) occurrence of threshold for non-linear-mixed type (i.e. non-linear with inflection points) is not prevalent within the fishing mortality rates explored. The conclusions drawn from the present study provide a knowledge base in indicators' dynamics

<footnote n="1">¹Pacific Biological Station, Fisheries and Oceans Canada, 3190 Hammond Bay Road, Nanaimo, BC V9T 6N7, Canada</footnote>

<footnote n="2">²Fraser River and Interior Area Stock Assessment, Fisheries and Oceans Canada, 3-100 Annacis Parkway, Delta, BC V3M 6A2, Canada</footnote>

<footnote n="3">³School of Aquatic and Fishery Sciences, University of Washington, Box 355020, Seattle, WA 98105-5020, USA</footnote>

<footnote n="4">⁴Bedford Institute of Oceanography, Fisheries and Oceans Canada, PO Box 1006, Dartmouth, NS B2Y 4A2, Canada</footnote>

<footnote n="5">⁵Department of Biological Sciences, University of Cape Town, Private Bag X3, Rondebosch 7701, South Africa</footnote>

<footnote n="6">⁶Marine Research Institute, University of Cape Town, Private Bag X3, Rondebosch 7701, South Africa</footnote>

<footnote n="7">⁷European Marine Board, Wandelaarkaai 7, Oostende 8400, Belgium</footnote>

<footnote n="8">⁸Scottish Association for Marine Science, Scottish Marine Institute, Oban, Argyll PA371QA, UK</footnote>

<footnote n="9">⁹Marine and Freshwater Research Centre (MFRC), Galway-Mayo Institute of Technology (GMIT), Dublin Road, Galway, Ireland</footnote>

<footnote n="10">¹⁰Institute of Marine Sciences, Middle East Technical University, PO Box 28, Erdemli 33731, Turkey</footnote>

<footnote n="11">¹¹Lowestoft Laboratory, Centre for Environment, Fisheries and Aquaculture Science (Cefas), Pakefield Road, Lowestoft, Suffolk NR33 0HT, UK</footnote>

<footnote n="12">¹²Institute of Marine Science (ICM-CSIC), Passeig Marítim de la Barceloneta, n° 37-49, 08003 Barcelona, Spain</footnote>

<footnote n="13">¹³Ecopath International Initiative Research Association, Mestre Nicolau 8, Bellaterra, 08193 Barcelona, Spain</footnote>

<footnote n="14">¹⁴CSIRO Oceans & Atmosphere, GPO Box 1538, Hobart, TAS 7001, Australia</footnote>

<footnote n="15">¹⁵Centre for Marine Socioecology, University of Tasmania, Hobart, TAS 7000, Australia</footnote>

<footnote n="16">¹⁶Institut de Recherche pour le Développement, MARBEC (IRD, University Montpellier, Ifremer, CNRS), Montpellier, France</footnote>

<footnote>*Corresponding author: tel: + 1 250 729 8373; fax: + 1 250 756 7053; e-mail: Caihong.Fu@dfo-mpo.gc.ca.</footnote>
under different fishing and primary productivity levels, thereby facilitating the application of ecosystem-based fisheries management worldwide.

**Keywords:** ecosystem-based fisheries management, generalized additive model, marine ecosystem model, non-linear response, primary productivity

## Introduction

Marine ecosystems are influenced by multiple stressors including climatic, oceanographic, ecological, and anthropogenic (e.g. fishing) changes. At the species level, these stressors may affect survival, growth, and reproduction, which can influence the structure and functioning of the entire ecosystem through complex species interactions. Ecosystem responses to stressors can be either linear or non-linear, depending, *inter alia*, on the properties of individual stressors and the way stressors interact. Particular attention needs to be paid to non-linear responses to stressors, in which thresholds (i.e. tipping points), may exist, beyond which large changes in ecosystem state or properties can occur abruptly (also described as regime shifts; Scheffer *et al.*, 2001; Scheffer and Carpenter, 2003; Groffman *et al.*, 2006; Andersen *et al.*, 2009; Kelly *et al.*, 2015; Moore, 2018). Identifying thresholds in non-linear responses to stressors can facilitate the implementation of risk-based management targets, the avoidance of detrimental ecosystem shifts and the optimization of socio-economic returns from ecosystem exploitation (Foley *et al.*, 2015; Kelly *et al.*, 2015).

In the present study, two ubiquitous stressors are considered: fishing pressure and primary productivity change, where the latter serves as a proxy for environmental change. Fishing has been a major stressor in exploited marine ecosystems around the globe for the last several decades (Boldt *et al.*, 2014), leading to substantial changes in the state of ecosystems (e.g. Jackson *et al.*, 2001). In parallel, changes in ocean temperature, circulation patterns, vertical mixing, and availability of nutrients have resulted in changes in the primary productivity of marine ecosystems (e.g. Brown *et al.*, 2011; Hunt *et al.*, 2011; Poloczanska *et al.*, 2016; Capuzzo *et al.*, 2018). Fishing pressure and primary productivity may work synergistically, affecting the structure and function of marine ecosystems (Kirby *et al.*, 2009; Möllmann *et al.*, 2009) and how they should be managed under future climate change conditions (Serpetti *et al.*, 2017; Baudron *et al.*, 2019).

To explore how an ecosystem responds (e.g. linearly vs. non-linearly) to fishing pressure and primary productivity change, various ecological indicators that represent ecosystem state and properties can be employed. Ecological indicators have been extensively studied using empirical data about their trends (e.g. Blanchard *et al.*, 2010), relationships with fishing pressure and environmental changes (e.g. Link *et al.*, 2010; Shannon *et al.*, 2010; Shannon *et al.*, 2014; Fu *et al.*, 2015), and potential thresholds (e.g. Large *et al.*, 2013, 2015; Tam *et al.*, 2017). Empirical studies have suggested that past and present exploitation strategies, mechanisms of productivity dynamics, and dominant ecological and environmental traits are all likely to influence how ecological indicators respond to various stressors (Shannon *et al.*, 2014). Consequently, it has long been a major challenge to draw general conclusions about the response patterns of indicators to fishing pressure and environmental changes.

To overcome this major challenge, we used a comparative multi-ecosystem, multi-model simulation experiment to explore the differential responses of ecological indicators to fishing mortality and primary productivity change. We took advantage of the work conducted by the IndiSeas working group (Shin *et al.*, 2012) in which researchers across different continents carried out unified simulation experiments for ten marine ecosystems around the globe. Prior to the simulation experiments, each of the ten ecosystems already had a working ecosystem simulation model developed (Shin *et al.*, 2018; Fu *et al.*, 2018, 2019), thus having overcome the time constraint and resource limitation of developing brand new ecosystem models. Relying on different ecosystem modelling approaches applied to different marine ecosystems allowed us to draw conclusions independent of the specific structure of ecosystems or the specific structure and assumptions of ecosystem simulation models.

The present study primarily addresses the following questions: (i) How do ecological indicators respond to fishing mortality in trends and directions under different types of primary productivity change? (ii) How frequently do non-linear responses to fishing mortality occur in different ecological indicators? (3) How frequently do thresholds occur in non-linear indicators' responses to fishing mortality? Results from the present study will provide insights on specific responses of indicators to fishing mortality; the thresholds of non-linear responses that we identify can inform the development of management strategies integrating environmental considerations, thereby facilitating the move towards ecosystem-based fisheries management worldwide.

## Material and methods

### Ecosystem models

In the present study, ten marine ecosystems are considered. The structure and functioning of these ecosystems were simulated using one of four different ecosystem simulation modelling approaches: Ecopath with Ecosim (EwE; Christensen and Walters, 2004), OSMOSE (Shin and Cury, 2004), Atlantis (Fulton *et al.*, 2004), or a multispecies size spectrum model (Blanchard *et al.*, 2014). Five of the ten ecosystems were modelled with EwE: the Black Sea (Akoglu, 2013; Akoglu *et al.*, 2014), the southern Benguela (Shannon *et al.*, 2004, 2008), the southern Catalan Sea (Coll *et al.*, 2006, 2013), the western Scotian Shelf (Araújo and Bundy, 2011, 2012), and western Scotland (Alexander *et al.*, 2015). OSMOSE was employed to model three of ten ecosystems: the West Coast of Canada (Fu *et al.*, 2013), the West Florida Shelf (Grüss *et al.*, 2016a, 2016b), and the Gulf of Gabes, Tunisia (Halouani *et al.*, 2016). Finally, the Southeastern Australian ecosystem was modelled with Atlantis (Fulton *et al.*, 2014), whereas the North Sea was modelled with the size spectrum model (Blanchard *et al.*, 2014). All the ecosystem models used in the present study have been published and validated against observed or estimated abundance, biomass and/or catch data. Details on the structure, assumptions, and species/taxa represented in the
ten ecosystem models are provided in Supplementary Tables S1 and S2.

## Fishing mortality

For each ecosystem, various levels of fishing mortality rate relative to $F_{MSY}$ (the fishing mortality rate corresponding to the maximum sustainable yield of a specific species/taxon within a specific ecosystem) were implemented to capture the impacts of different fishing levels. The same design was carried out consistently across all ecosystems for valid comparisons. Prior to running the experimental simulation within each ecosystem, $F_{MSY}$ was estimated for each exploited species/taxon, by constructing the yield to fishing mortality curve (fisheries yield as a function of fishing mortality rate) of that species/taxon while holding the fishing mortality rate of all other species/taxa constant at their respective current fishing mortality rates ($F_{curr}$). Then, during the experimental simulation runs, fishing mortality rates were incrementally increased relative to $F_{MSY}$ for each species/taxon of interest through the use of a multiplier $\lambda$. In other words, fishing mortality rates $F_i = \lambda \times F_{MSY}^i$ (year$^{-1}$) were applied to each species/taxon of interest $i$, where $\lambda \in \{0.25, 0.5, 0.75, 1, 1.25, 1.5\}$. This range of fishing mortality rates covers representative values for the yield to fishing mortality curves (Fu et al., 2018; Shin et al., 2018).

For each ecosystem, we consider three different fishing strategies: a “low-trophic-level” (F\_ltl), a “high-trophic-level” (F\_htl), and an “all-trophic-level” (F\_all) strategy. With the F\_ltl fishing strategy, fishing experiments focus on all forage species/taxa retained in commercial or subsistence fisheries. Here, forage species/taxa are defined as species/taxa whose adults mainly feed on plankton (phyto-, zoo-, or ichthyoplankton) or similar low-trophic level organisms (e.g. meiobenthos). With the F\_htl fishing strategy, fishing experiments focus on predatory taxa, which include large demersal and large pelagic taxa. Finally, the F\_all fishing strategy represents broad-scale exploitation, where the focus taxa are all taxa retained in commercial or subsistence fisheries, including both high-trophic-level (HTL) and low-trophic-level (LTL) taxa. Any pre-recruit stages that are represented in the ecosystem models are excluded from the fishing scenarios, as are marine mammals, marine turtles, and seabirds. The HTL and LTL species/taxa represented in each of the ten ecosystem models are detailed in Supplementary Table S1.

## Primary productivity

For all ten ecosystems, changes in phytoplankton biomass, comparable across models and ecosystems, were used to represent changes in primary productivity. For the Atlantis ecosystem model, however, it is not possible to directly scale phytoplankton biomass, because of the role of phytoplankton in the biogeochemical cycles represented in Atlantis. Therefore, in this case, time-series of nutrient inputs were scaled so that the resulting changes in phytoplankton biomass were in line with those implemented in other models.

Two types of changes in primary productivity were simulated: directional and random. For directional changes in primary productivity, a multiplier $\gamma \in \{0.85, 0.9, 0.95, 1, 1.05, 1.1\}$ was directly applied to modelled phytoplankton biomass. This range of variability encompasses the range of changes observed in the ten ecosystems in the last decade (Boyce et al., 2014). With respect to random changes in primary productivity, the modelled phytoplankton biomass was multiplied by a random multiplier drawn

from a lognormal distribution with a mean $\mu$ of 1 and a range of standard deviations $\sigma = \{0.1, 0.2, 0.3\}$. These standard deviation values are consistent with the observed annual satellite-derived chlorophyll-a concentration levels from the MODIS (Moderate Resolution Imaging Spectroradiometer) Aqua spectral data (Shin et al., 2018). A set of 30 random multipliers was generated for each value of $\sigma$, so as to adequately sample the random distribution.

Consequently, within each ecosystem, 12 combinations of fishing strategies (F\_ltl, F\_htl, and F\_all) and primary productivity change (directional change, and random change while assuming $\sigma = 0.1, 0.2, \text{or } 0.3$) were simulated to explore interactions of these two factors.

## Ecological indicators

Following Fu et al. (2019), 14 ecological indicators were explored (Table 1), including biomass to fisheries catch ratio (B/C), proportion of predatory fish (Pred), mean intrinsic vulnerability (IVI), mean lifespan (Lifesp), mean trophic level of the community (TLco), and marine trophic index (MTI), the trophic level (TL) of catch calculated assuming a constant TL for taxa (TLc), the TL of catch calculated assuming a variable TL for taxa (TLcVar; Reed et al., 2017), and six community-level biomass indicators [biomass of all-trophic-level taxa (B_all), the biomass of HTL taxa (B HTL), the biomass of LTL taxa (B LTL), the ratio of B HTL to B ALL (B H2A), the ratio of B LTL to B ALL (B L2A), and the ratio of B LTL to B HTL (B L2H)]. Aside from the six

<table><caption>Table 1. List of the indicators considered in this study.</caption><thead><tr><th>Indicator</th><th>Definition</th><th>Abbreviation</th></tr></thead><tbody><tr><td>Biomass to catch ratio</td><td>B/C</td><td>B/C</td></tr><tr><td>Proportion of predatory fish</td><td>B (predatory fish)/<br>B (surveyed)</td><td>Pred</td></tr><tr><td>Mean intrinsic vulnerability</td><td>Σ<sub>i</sub> IV<sub>i</sub>C<sub>i</sub>/Σ<sub>i</sub>C<sub>i</sub></td><td>IVI</td></tr><tr><td>Mean lifespan</td><td>Σ<sub>i</sub> Age<sub>max</sub>,<sub>i</sub> B<sub>i</sub>/Σ<sub>i</sub>B<sub>i</sub></td><td>Lifesp</td></tr><tr><td>Trophic level of catch</td><td>Σ<sub>i</sub> TL<sub>i</sub>C<sub>i</sub>/Σ<sub>i</sub>C<sub>i</sub></td><td>TLc</td></tr><tr><td>Trophic level of catch with<br>variable TL</td><td>Σ<sub>i</sub> TL'<sub>i</sub>C<sub>i</sub>/Σ<sub>i</sub>C<sub>i</sub></td><td>TLcVar</td></tr><tr><td>Marine trophic index</td><td>Σ<sub>i</sub> s(TL<sub>i</sub>>3.25)/Σ<sub>i</sub>s(TL<sub>i</sub>>3.25) TL<sub>i</sub>C<sub>i</sub>/Σ<sub>i</sub>C<sub>i</sub></td><td>MTI</td></tr><tr><td>Mean trophic level of<br>community</td><td>Σ<sub>i</sub> TL<sub>i</sub>B<sub>i</sub>/Σ<sub>i</sub>B<sub>i</sub></td><td>TLco</td></tr><tr><td>Biomass of all surveyed species</td><td>Σ<sub>i</sub>s(all) B<sub>i</sub></td><td>B_all</td></tr><tr><td>Biomass of high-trophic-level<br>(htl) species</td><td>Σ<sub>i</sub>s(htl) B<sub>i</sub></td><td>B HTL</td></tr><tr><td>Biomass of low-trophic-level<br>(ltl) species</td><td>Σ<sub>i</sub>s(ltl) B<sub>i</sub></td><td>B LTL</td></tr><tr><td>Ratio of htl biomass to total<br>biomass</td><td>B<sub>htl</sub>/B<sub>all</sub></td><td>B_H2A</td></tr><tr><td>Ratio of ltl biomass to total<br>biomass</td><td>B<sub>ltl</sub>/B<sub>all</sub></td><td>B_L2A</td></tr><tr><td>Ratio of ltl biomass to htl<br>biomass</td><td>B<sub>ltl</sub>/B<sub>htl</sub></td><td>B_L2H</td></tr></tbody></table>

High-trophic-level and low-trophic-level taxa for all the ecosystems are listed in Supplementary Table S1.

B, biomass (in tonnes); C, catch (in tonnes); s, species; TL, trophic level; TL', variable TL; IVI, intrinsic vulnerability index.
community-level biomass indicators, the indicators Pred, Lifesp, TLco were also derived from survey biomass data (i.e. biomass-based); in contrast, the indicators IVI, MTI, TLc, and TLcVar were calculated from catch data (i.e. catch-based). Within each of the 10 ecosystem models, these 14 ecological indicators were averaged over the last 10 years of simulation period under each scenario of varying fishing mortality rate and primary productivity, which implies that these indicator data points corresponding to different fishing mortality rates and primary productivity levels are independent of each other.

## Generalized additive models: linearity/non-linearity and directions

Nonparametric regression methods such as generalized additive models (GAMs; Hastie and Tibshirani, 1990) are powerful tools for exploring linear or non-linear responses of a variable to predictors without being constrained to an underlying parametric model of a specific form, which is particularly useful when ecological thresholds of non-linear responses are of interest (Lintz et al., 2011; Nicolaou and Constandinou, 2016). Therefore, GAMs were employed in this study to explore the responses of ecological indicators to fishing mortality and primary productivity change. Each of the two stressors was modelled individually in order to avoid convergence issues because of the small number of fishing mortality and primary productivity levels considered in this study, as well as to avoid potential correlation between fishing mortality and primary productivity. The GAM model is formulated as (Large et al., 2013):

$$Y = \alpha + s(X) + \varepsilon, \quad (1)$$

where Y is the ecological indicator; α is the intercept held constant; X is the predictor, representing either fishing mortality or primary productivity level; s() is the smoothing function; and ε is the error term of a normal distribution with zero mean and finite variance. The smoothing parameter was estimated using the generalized cross-validation criterion (GCV criterion; Wood et al., 2016) implemented in R package “mgcv” (Version 1.8-24; Wood, 2017). Thin plate regression spline smoothing terms with a ridge penalty were employed, such that the smoothing term could be eliminated and the predicted curve could be reduced to a linear function when the smoothing term did not improve fit. Confidence intervals were calculated for each GAM with a naive bootstrap using random sampling and replacement [see Large et al. (2013) for more details on the naive bootstrap approach].

For the purpose of exploring linear or non-linear responses of ecological indicators to fishing mortality, the present study focuses on the development and exploration of 10 × 14 × 3 × 4 = 1680 GAMs, covering 10 ecosystems, 14 ecological indicators, 3 fishing strategies (F_ltl, F HTL, and F_all), and 4 scenarios of primary productivity change (directional change, and random change while assuming σ = 0.1, 0.2, or 0.3). The performance of a GAM model is often gauged by examining the amount of variance it explains (i.e. the proportion of the deviation in the data explained). Comparisons of GAM predictions across the different scenarios of primary productivity change provide perspectives on how the nature and degree of environmental variability affects the capacity of GAMs to capture the detrimental impacts of fishing, as well as whether indicators'

responses to fishing pressure are robust to varying levels of environmental variability.

Predicted GAM curves based on empirical time-series have often been explored for non-linear responses to drivers using first and second derivatives (e.g. Large et al., 2013; Burthe et al., 2016). However, empirical time-series are typically faced with statistical problems such as temporal autocorrelation, confounding among unknown sources of variability and inadequacy in replication, resulting in difficulty in identifying linearity/non-linearity (Litzow and Hunsicker, 2016). In contrast to empirical time-series data explored in Large et al. (2013) and Burthe et al. (2016), the ecological indicator data in the present study were obtained through individual simulations with different fishing mortality and primary productivity and are thus independent data points and not time-series, as stated above. Thus, the linearity/non-linearity was explored in a more literal sense using specifically designed methods for these datasets, as described below.

For each of the 1680 GAMs with respect to fishing mortality, we first determine whether the indicator's response to fishing mortality is linear by fitting a linear regression to the GAM predicted values. If the correlation coefficient between the GAM predicted values and the linear regression fitted values is >0.95, the indicator's response to fishing mortality (derived from the GAM) is considered linear. Otherwise, the indicator's response to fishing mortality is assumed to be non-linear. If the relationship between the indicator and fishing mortality is linear and if the slope of the linear regression is statistically significant (p < 0.05), the direction of indicator's change in response to fishing mortality is defined as a “linear increase” when the slope is >0 and a “linear decrease” when the slope is <0.

When the relationship between the indicator and fishing mortality is non-linear, the GAM curve is further investigated for identifying monotonicity and determining the direction of indicator's change in response to fishing mortality by examining the first derivative of the curve. A “non-linear-monotonic increase” is associated with all first derivatives being non-negative; whereas a “non-linear-monotonic decrease” is associated with all first derivatives being non-positive. However, if the first derivatives have both positive and negative values, the GAM curve is defined as non-linear-mixed. For this type of non-linearity, the direction of an indicator's change in response to fishing mortality is determined by fitting a linear regression to GAM predicted values. When a significant linear relationship exists (i.e. the correlation coefficient between the GAM predicted values and the linear regression fitted values is >0.95), the non-linear-mixed curve is defined as a “non-linear-mixed increase” when the linear regression fitted slope is >0 and as a “non-linear-mixed decrease” when the linear regression fitted slope is <0.

## Thresholds

For non-linear-mixed GAM curves, there are regions along the fishing mortality gradient where the first derivative of the response function changes sign at an inflection point. Inflection points on GAM curves have been considered as ecological thresholds (e.g. Large et al., 2013; Tam et al., 2017), and such ecological thresholds, also known as tipping points (e.g. Moore, 2018), may have implication for fisheries management (e.g. indicate the occurrence of a regime shift) and, therefore, deserve special attention. Because of the fact that there are only six discrete fishing
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems under the low-trophic-level (F_ltl) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

mortality rates ($F_{MSY}$ multiplier = 0.25, 0.5, 0.75, 1.0, 1.25, 1.5), we do not use the threshold-determination approach developed by Large et al. (2013), rather, we fit a fifth order polynomial model [degree of freedom = 6, Equation (1)] in all the cases where the indicator's response to fishing mortality is identified as non-linear-mixed:

$$f(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + a_4x^4 + a_5x^5, \quad (2)$$

where x is the predictor, representing the levels of fishing mortality; f(x) is the ecological indicator as a function of x; and $a_0, \dots, a_5$ are coefficients of the polynomial model. The analytical root(s) of the second derivatives, where ecological thresholds reside, are derived from:

$$f''(x) = 2a_2 + 6a_3x + 12a_4x^2 + 20a_5x^3 = 0. \quad (3)$$

The analytical root(s) ($x = x_i$) are solved using the “rootSolve” package in R (R Core Team, 2019). Only root(s) within the range of $F_{MSY}$ multiplier (0.25~1.25) are retained for consideration.

In addition, root(s) with their predicted values $f(x_i)$ being outside the range defined by the indicator value on either side of the $x_i$ value are not considered as thresholds, as they are inflection points resulting from local maxima or minima artificially introduced by high-order polynomial fits.

# Results

## Variance explained

The amount of variance explained by fishing mortality (i.e. $R_F^2$), though varying among the ten ecosystems under each of the three fishing strategies (Figures 1–3), has some common themes. First, the $R_F^2$ values for all indicators tend to be the highest under the F_all fishing strategy with all four scenarios of primary productivity change across all ecosystems, whereas those under the F_ltl fishing strategy tend to be the lowest (the higher the $R_F^2$ values, the stronger the indication of fishing impacts). Within a given ecosystem, an individual indicator can have drastically different $R_F^2$ values among the three fishing strategies, implying that fishing patterns are important when interpreting the dynamics of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems under the high-trophic-level (F_htl) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

ecological indicators. Second, the $R_F^2$ value of a specific indicator under any fishing strategy is generally smaller under directional than under random primary productivity change, and the $R_F^2$ value consistently decreases with increasing primary productivity variability (i.e. when σ increases from 0.1 to 0.3), suggesting that an indicator's response to fishing mortality can be increasingly obscured when environmental variability increases. Third, comparing across the 14 indicators, the biomass indicators B_all, B_htl, and B_ltl tend to have the lowest $R_F^2$ values under each of the three fishing strategies, which can be more clearly seen from their values averaged over the ten ecosystems (Figure 4). These ecosystem-averaged $R_F^2$ values indicate the ability of the 14 indicators to reflect fishing impacts. Fourth, the $R_F^2$ value for the indicator B/C under the F_all fishing strategy is close to one under all four scenarios of primary productivity change for all ecosystems except the southern Catalan Sea, indicating that the trends for the indicator B/C are usually well explained by fishing mortality when both HTL and LTL taxa are targeted by fishing.

## Responses of indicators to fishing mortality: linearity/non-linearity

Overall, the shapes of the GAM curves (linear, non-linear-monotonic, or non-linear-mixed) for most indicators are consistent among the four scenarios of primary productivity change (Figures 5–7), suggesting that most indicators’ responses to fishing mortality are largely insensitive to primary productivity variability. In contrast, the shapes of the GAM curves for most indicators across the ten ecosystems tend to vary with the three different fishing strategies, indicating the importance of fishing patterns when considering the responses of ecological indicators to fishing mortality.

Across the ten ecosystems, the response to fishing mortality is largely linear (over 50% for most indicators; Figure 8). Within a given ecosystem, the shapes of GAM curves tend to differ among the 14 indicators. The catch-based indicators (e.g. B/C, IVI, TLc, MTI) tend to have larger proportion of non-linear responses. In contrast, the community and biomass-based indicators
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators in the 10 study marine ecosystems under the all-trophic-level (F_all) fishing strategy under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Bar plots of variances explained from the generalized additive models (GAMs; with fishing mortality as the predictor) for each of the 14 indicators for the 10 study marine ecosystems that were averaged over all the 10 ecosystems under each of the 3 fishing strategies: (top) F_ltl (the low-trophic-level fishing strategy), (middle) F_htl (the high-trophic-level fishing strategy), and (bottom) F_all (the all-trophic-level fishing strategy) under 4 different scenarios of primary productivity (directional, σ = 0.1, σ = 0.2, and σ = 0.3). The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

(e.g. TLco, B_all, B_htl, B_ltl) are largely linear in their response to fishing mortality (Figure 8), which may suggest that the responses of the community and biomass-based indicators to fishing mortality are more predictable.

## Responses of indicators to fishing mortality: direction

Similar to the shapes of GAM curves, one common theme across the ten ecosystems is that the response direction of GAM curves (“decrease” or “increase”) for a specific indicator depends on the type of fishing strategy (Figures 5–7 and 9). Under the F_ltl fishing strategy, the indicators Pred, Lifesp, and B_H2A show increasing GAM trends in more than half of the study ecosystems under both directional and random primary productivity changes. However, these three indicators all show decreasing GAM trends in most ecosystems under the F_htl and F_all fishing strategies. Under the F_htl fishing strategy, there are six indicators (IVI, TLc, MTI, B_ltl, B_L2A, and B_L2H) that show an increasing trend with increasing fishing mortality in over 50% of the ecosystems under both directional and random primary productivity changes. However, these six indicators predominantly
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the low-trophic-level (F_ltl) fishing strategy using fishing mortality as the predictor under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

decrease with increasing fishing mortality under both the F_ltl and F_all fishing strategies. Under the F_all fishing strategy, the indicator B_L2H is the only one that predominantly shows an increasing GAM trend with increasing fishing mortality, demonstrating that the response direction of all other indicators becomes more consistent (i.e. always decreasing with increasing fishing mortality) when fishing targets both HTL and LTL taxa (i.e. under the F_all fishing strategy). It is worth noting that, among all the indicators, B/C is the only one that predominantly exhibits a decreasing trend with increasing fishing mortality under both directional and random primary productivity changes and under all three fishing strategies, which indicates that B/C is the most informative indicator of fishing impacts.

## Thresholds

Overall, under all three fishing strategies, the occurrence of threshold for the non-linear-mixed type is not prevalent for most ecosystems (Figures 5–7; Supplementary Figures S1–S3). Exceptions include the Black Sea (under the F_ltl fishing strategy), Catalan Sea (under the F HTL fishing strategy), and Southeastern Australia and West Coast of Canada (under the F_all fishing strategy), where more than 50% of the non-linear-mixed GAM curves have thresholds (Table 2). The infrequency of threshold occurrences may imply that the root(s) of second derivatives either fall outside the fishing mortality range (0.25–1.5\*F\_MSY) for most non-linear-mixed type (aside from the possibility of imaginary roots), indicating that future
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the high-trophic-level (FHTL) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with $\sigma = 0.1$, $\sigma = 0.2$, and $\sigma = 0.3$, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

simulations may desire a broader range of fishing mortality rate in order to avoid missing true ecological thresholds.

The occurrence of inflection points resulting from local maxima or minima of the fifth order polynomial fits (not signalling thresholds) is rare except in the Southeastern Australia ecosystem, where there are 16 such type of inflection points (consisting 35% of the non-linear-mixed type contrasting with 52% of thresholds; Table 2). This implies that thresholds identified within the fishing mortality range are generally appropriate.

## Discussion

The multi-ecosystem, multi-model simulation experiment conducted within the IndiSeas programme (Shin et al., 2018; Fu et al., 2019) has generated a comprehensive simulation database

of ecological indicators under two ecosystem stressors, fishing pressure, and primary productivity change. The GAMs employed in the present study have allowed us to produce response curves for 14 major ecological indicators evaluated for 10 ecosystems in response to fishing mortality, revealing linearity or non-linearity in the responses, with decreasing or increasing trends and thresholds associated with non-linear responses. This comprehensive and comparative research thus provides a rich knowledge base of commonalities and differences among ecosystems in indicators' dynamics and performance at different levels of fishing mortality under different primary productivity scenarios, as well as thresholds associated with non-linear responses of indicators, thereby providing essential knowledge to facilitate the move towards ecosystem-based fisheries management.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Schematic plots of the shapes, including linear (shown as line), non-linear-monotonic (solid curve), and non-linear-mixed (dashed curve), and directions, including increase (upward arrowhead), decrease (downwards arrowhead), and no significant trend (no arrowhead), of the generalized additive models (GAMs) for each of the 14 indicators for the 10 marine ecosystems under the all-trophic-level (F_all) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with σ = 0.1, σ = 0.2, and σ = 0.3, respectively). The locations of thresholds are shown as dark dots and those for inflection points that do not constitute thresholds are shown as light dots. The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

## Comparison with previous work

Using the gradient forest method, Fu et al. (2019) investigated the performance of the same set of ecological indicators about their sensitivity, specificity, and threshold responses to both fishing pressure and primary productivity change. However, the gradient forest method does not allow quantitative descriptions of the directions (“decrease” or “increase”) and shapes of the indicators’ responses to pressures (e.g. fishing mortality). Samhouri

et al. (2017) employed a multi-model inference framework, a combination of the gradient forest method and GAMs, for exploring directions, shapes, and thresholds in the responses of ecosystem state indicators to pressures of environmental change and human activities. The present study also uses GAMs to expand the previous work with the gradient forest method (Fu et al., 2019). Additionally, the present study has expanded the investigation from one fishing strategy (F_all) in Fu et al. (2019) to all
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Proportion of non-linear GAM (generalized additive model) types (including monotonic and mixed) for each of the 14 indicators averaged over the 10 marine ecosystems under each of the 3 fishing strategies (F_ltl for the low-trophic-level fishing strategy, F HTL for the high-trophic-level fishing strategy, and F_all for the all-trophic-level fishing strategy) under 2 types of primary productivity changes: (top) directional, and (bottom) random (the results for random were averaged over 3 levels of variability with σ = 0.1, 0.2, and 0.3). The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

<table><caption>**Table 2.** Total number of non-linear-mixed GAM (generalized additive model) curves and identified thresholds for each of the ten ecosystems and under each of the three fishing strategies (F_ltl, F HTL, F_all), as well as the percentage of non-linear-mixed GAM curves that have thresholds.</caption><thead><tr><th rowspan="2">Ecosystem model</th><th colspan="3">Number of non-linear-mixed GAM</th><th colspan="3">Number of thresholds</th><th colspan="3">Number of fake thresholds</th><th colspan="3">Percentage of threshold occurrences</th></tr><tr><th>F_ltl</th><th>F HTL</th><th>F_all</th><th>F_ltl</th><th>F HTL</th><th>F_all</th><th>F_ltl</th><th>F HTL</th><th>F_all</th><th>F_ltl</th><th>F HTL</th><th>F_all</th></tr></thead><tbody><tr><td>BlackSea-EwE</td><td>21</td><td>8</td><td>38</td><td>19</td><td>5</td><td>18</td><td>1</td><td>0</td><td>4</td><td><strong>90</strong></td><td>63</td><td>47</td></tr><tr><td>GulfGabes-OSMOSE</td><td>16</td><td>22</td><td>9</td><td>8</td><td>2</td><td>2</td><td>0</td><td>0</td><td>0</td><td>50</td><td>9</td><td>22</td></tr><tr><td>NorthSea-SizeS</td><td>24</td><td>7</td><td>11</td><td>11</td><td>4</td><td>7</td><td>0</td><td>0</td><td>0</td><td>46</td><td>57</td><td>64</td></tr><tr><td>CatalanSea-S-EwE</td><td>22</td><td>23</td><td>26</td><td>9</td><td>15</td><td>5</td><td>0</td><td>0</td><td>0</td><td>41</td><td><strong>65</strong></td><td>19</td></tr><tr><td>Australian.SE-Atlantis</td><td>26</td><td>12</td><td>46</td><td>8</td><td>0</td><td>24</td><td>0</td><td>0</td><td>16</td><td>31</td><td>0</td><td><strong>52</strong></td></tr><tr><td>Benguela-S-EwE</td><td>9</td><td>4</td><td>16</td><td>6</td><td>1</td><td>4</td><td>0</td><td>0</td><td>1</td><td>67</td><td>25</td><td>25</td></tr><tr><td>Canada.W-OSMOSE</td><td>6</td><td>6</td><td>20</td><td>2</td><td>4</td><td>19</td><td>0</td><td>0</td><td>0</td><td>33</td><td>67</td><td><strong>95</strong></td></tr><tr><td>Scotland.W-EwE</td><td>12</td><td>22</td><td>7</td><td>1</td><td>3</td><td>1</td><td>1</td><td>0</td><td>0</td><td>8</td><td>14</td><td>14</td></tr><tr><td>FloridaShelf.W-OSMOSE</td><td>7</td><td>10</td><td>12</td><td>1</td><td>5</td><td>4</td><td>0</td><td>0</td><td>0</td><td>14</td><td>50</td><td>33</td></tr><tr><td>ScotianShelf.W-EwE</td><td>21</td><td>21</td><td>4</td><td>0</td><td>2</td><td>3</td><td>0</td><td>0</td><td>1</td><td>0</td><td>10</td><td>75</td></tr></tbody></table>

The scenarios with more than 50% chance of having thresholds are bolded, excluding those that have <19 (one-third of 56) non-linear-mixed GAM curves within a particular ecosystem.

three fishing strategies (F_ltl, F HTL, F_all). The expansion has been proven to be valuable as it has led to the important discovery that for all ten ecosystems, all aspects of the GAM curves, including model performances in the form of variance explained (Figures 1–3), shapes (Figures 5–7), and the response directions (Figure 9) are dependent on the different fishing strategies considered. Despite the importance of considering fisheries exploitation history, as pointed out by Shannon et al. (2014) and Fu et al. (2018), ecological indicators have hitherto been investigated

largely in the absence of such knowledge (e.g. Blanchard et al., 2010; Coll et al., 2016). Future data collection and analysis related to fisheries exploitation should be more conscious about fishing patterns, i.e. what trophic levels fisheries target.

## Responses of indicators to fishing mortality: variance explained

Our results demonstrate that responses of indicators to fishing mortality differ across fishing strategies in all ten ecosystems
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Proportion of GAM (generalized additive model) directions (increase and decrease) for each of the 14 indicators averaged over the 10 marine ecosystems under each of 3 fishing strategies (F_ltl for the low-trophic-level fishing strategy, F HTL for the high-trophic-level fishing strategy, and F_all for the all-trophic-level fishing strategy) under 2 types of primary productivity changes: (top) directional, and (bottom) random (the results for random were averaged over 3 levels of variability with σ = 0.1, 0.2, and 0.3). The full name and meaning of ecological indicators are provided in Table 1.</figcaption></figure>

(Figures 1–3, 5–7, and 9). In particular, the $R_F^2$ values for all indicators tend to be highest under the F_all fishing strategy whereas lowest under the F_ltl fishing strategy, which may suggest that the impacts of fisheries exploitation on LTL taxa could be obscured because of the relatively closer coupling between the dynamics of LTL taxa and primary productivity. In order to further understand the commonalities and differences among the ten ecosystems, we explored how they would respond to primary productivity change. We conducted an additional 1680 GAMs (10 ecosystems × 14 ecological indicators × 3 fishing strategies × 4 scenarios of primary productivity change) with primary productivity as the predictor. Unlike the GAMs with fishing mortality as the predictor, the GAMs with primary productivity as the predictor generally resulted in low $R_p^2$ under all fishing strategies across all ecosystems with only a few exceptions (Supplementary Figures S10–S12), which add insight to the results for fishing mortality.

One interesting exception is the southern Benguela ecosystem under the F_ltl fishing strategy (Supplementary Figure S10). The southern Benguela is an upwelling ecosystem with abundant lower trophic level taxa. Upwelling indices of this ecosystem have shown increased variability in the 1990s and 2000s, presumably reflecting climate change (Blamey et al., 2012), which have caused ecosystem changes mediated through the lower trophic level fish community of this wasp-waist foodweb (Cury et al., 2000). The exceptionally high $R_p^2$ values (but low $R_F^2$ values) for most indicators under the F_ltl fishing strategy suggests that this ecosystem is extremely sensitive to changes in primary productivity (see Ortega-Cisneros et al., 2018), and any potential ecosystem change evoked through lower trophic level fisheries is primarily driven

by environmental changes. This conclusion is consistent with previous indicator-based studies that have concluded that interpretation of ecological indicator trends in response to fishing in the southern Benguela requires careful unpacking of changes in environmental conditions in this ecosystem (Shannon et al., 2010; Lockerbie et al., 2016). The importance of environmental drivers in upwelling systems has been highlighted in the northern Benguela ecosystem, where environmental conditions and overfishing together have caused apparent regime shifts (Heymans and Tomczak, 2016).

The Black Sea ecosystem also displays unique features, particularly the large discrepancies in the GAM curves between the directional and random primary productivity change scenarios indicating that this ecosystem is very sensitive (or responsive) to changes in primary productivity, more so than to changes in fishing mortality, when the primary productivity changes are directional. Similar to the southern Benguela ecosystem, the changes in the Black Sea ecosystem are predominantly mediated by the wasp-waist control exerted on the foodweb by small pelagic fish with a significant interplay of bottom-up control (Akoglu et al., 2014), which determines the endurance of small pelagic fish under the influence of fisheries overexploitation. However, the large discrepancies disappear when comparing the results of the three different scenarios of random primary productivity changes (where $\sigma = 0.1, 0.2, \text{ or } 0.3$), which suggests that the indicators considered in the present study are robust to environmental variability for the Black Sea ecosystem.

The West Florida Shelf ecosystem contrasts with the southern Benguela ecosystem in that the West Florida Shelf has very high $R_F^2$ values for most indicators and very low $R_p^2$ values for all
indicators. These results concur with the findings of Grüss et al. (2016a, 2016b), which show that the West Florida Shelf ecosystem is under such a strong top-down control that elevated natural mortality because of large environmental changes (e.g. harmful algal blooms) has negligible impacts on ecosystem structure and functioning compared with elevated mortality resulting directly or indirectly (via trophic interactions) from changes in fishing patterns. Other ecosystems such as the western Scotian Shelf and West Coast Scotland have an intermediate response. Although the community-level biomass indicators (B_ltl, B HTL, B_all) tend to have higher Rₚ² values compared with Rₓ, other indicators, such as B/C and IVI, indicate strong responsiveness to fishing mortality. These ecosystems have a long exploitation history; however, their temperate location may be providing some resilience to exploitation (Frank et al., 2007).

In contrast to other ecosystems in this study, the southern Catalan Sea stands out as one that has quite different indicator dynamics in response to fishing mortality and primary productivity change. This ecosystem has a long fishing history, and current fishing mortality rates in the southern Catalan Sea (northwestern Mediterranean Sea) are very high (Fernandes et al., 2017). All indicators except for Lifesp, TLc, TLcVar, and MTI have very low R<sub>F</sub><sup>2</sup> values under both directional and random primary productivity changes. Essentially, the southern Catalan Sea ecosystem has been fished down and the foodweb of this ecosystem is currently under bottom-up control. Consequently, the variations in many of the ecological indicators are no longer responsive to changes in fishing pressure, and are instead primarily driven by directional changes in primary productivity with high R<sub>p</sub><sup>2</sup> for most indicators when there is directional primary productivity change, as already highlighted in previous studies (Shannon et al., 2014; Coll et al., 2016; Lockerbie et al., 2017).

As shown above, the five ecosystems (Black Sea, southern Benguela, southern Catalan Sea, western Scotian Shelf, and western Scotland) demonstrate diverse responses of indicators to fishing mortality and primary productivity change, despite the same ecosystem modelling approach EwE applied to these ecosystems. This suggests that the types of ecosystem modelling approach may not contribute to the differing results among the ten ecosystems. Because of the commonly low R<sub>p</sub><sup>2</sup> of the GAMs with primary productivity as the predictor, there is no need to explore other properties of the GAMs, such as shapes and directions. However, the generally low R<sub>p</sub><sup>2</sup> values could be caused by the rather narrow range of primary productivity (multiplier = 0.85, 0.9, 0.95, 1, 1.05, 1.1). Future research with a broader range of primary productivity scenarios are warranted for investigating the responses of ecological indicators to extreme climate change (primary productivity change in this case).

## Responses of indicators to fishing mortality: direction

The majority of the GAM curves demonstrate decreasing trends with increasing fishing mortality (Figures 5–7 and 9; Supplementary Figures S4–S9). These decreasing trends reflect the expected behaviour of indicators in response to increasing fishing mortality (Rice and Rochet, 2005; Bundy et al., 2010; Shin et al., 2010). However, the trends (“decrease” or “increase”) of the GAM curves of a specific indicator can depend on the type of fishing strategy considered. For instance, the indicators Pred, Lifesp, and B_H2A increase with increasing fishing mortality under the F_ltl fishing strategy, but decrease under the F HTL and

F_all fishing strategies. This implies that caution is required when interpreting the dynamics of these indicators when fishing primarily targets lower trophic level taxa, and even more so in ecosystems such as upwelling systems (e.g. the southern Benguela) where environmental variability strongly influences ecological indicators.

On the other hand, other indicators, including IVI, TLc, MTI, B_ltl, B_L2A, and B_L2H, increase with increasing fishing mortality under the F HTL fishing strategy, but decrease under F_ltl and F_all in most ecosystems. The decrease of indicator TLc in response to fishing has been observed in many ecosystems (Pauly and Watson, 2005; Fu et al., 2012), resulting from “fishing down the foodweb.” It is likely that the ecosystems for which these studies were carried out are ecosystems where fishing activities primarily target HTL taxa. However, the present study suggests that the indicator TLc (along with IVI, TLc, MTI, B_ltl, B_L2A, and B_L2H) may increase in response to increasing fishing mortality if fishing activities in the ecosystem of interest primarily target LTL taxa, providing potentially false assurance of ecosystem health. Therefore, our multi-ecosystem, multi-model simulation experiment supports the assertion of Branch (2015) and Shannon et al. (2014), i.e. multiple working hypotheses are needed to explore how fishing affects marine foodwebs rather than assuming “fishing down the foodweb” applies in all cases.

In contrast to the situations under the F_ltl and F HTL fishing strategies, all indicators except B_L2H decrease with increasing fishing mortality for most ecosystems under the F_all fishing strategy. Therefore, we conclude that when fishing expands over time on all-trophic-level taxa, the direction of change in ecological indicators is more predictable, with a generally decreasing trend in response to increasing fishing levels emerging.

## Responses of indicators to fishing mortality: linearity/non-linearity and thresholds

Identifying non-linearity is a premise for applying early warning signals to ecosystem-based fisheries management that can be used in a similar way to limit reference points in the precautionary approach (Gabriel and Mace, 1999). However, formal tests of non-linearity are often hampered by typically short and autocorrelated time-series in empirical studies (Litzow and Hunsicker, 2016). Most studies on non-linearity have been based on the identification of threshold responses on GAMs curves fitted to empirical relationships between time-series of response (e.g. indicators) and drivers (e.g. environmental variables and fishing pressure; Large et al., 2013; Samhouri et al., 2017; Tam et al., 2017). The present study explores non-linearity in a literal sense because it deals with independent data points from ecosystem simulations for which autocorrelation is not a concern. Our results suggest that most indicators tend to respond to fishing mortality in a linear way, particularly for the community and biomass-based indicators (e.g. TLco, B_all, B HTL, B_ltl; Figure 8). Although the indicator B/C predominantly responds to fishing mortality in a non-linear way, the majority of trends in the B/C indicator is non-linear-monotonic under the F_all fishing strategy with distinct response direction and straightforward interpretation of response to fishing mortality. This is in line with the suggestion by Litzow and Hunsicker (2016) that most observed ecosystem changes may be parsimoniously explained by linear and reversible tracking of perturbation.
For potentially non-linear relationships between response and drivers, GAMs and/or the gradient forest method have often been employed to quantify thresholds (e.g. Large et al., 2013, 2015; Samhoury et al., 2017; Tam et al., 2017; Fu et al., 2019). The previous studies of quantifying thresholds of ecological indicators in their responses to fishing pressure and environmental change have contributed important knowledge to identify operational reference points for moving towards ecosystem-based fisheries management by avoiding regime shifts caused by climate change and/or overfishing (Tam et al., 2017). Although the previous studies identified all inflection points (second derivative of the GAMs curve = 0) or significant shift in the gradient forest analysis as thresholds, the present study deliberately focuses only on the cases where fitted GAMs curves are non-linear-mixed and ignores all other cases of linear tracking and monotonic non-linearity. Fu et al. (2019), using the 14 indicators explored here, concluded that under the F_all fishing strategy, all 14 indicators had threshold responses and over 50% of the threshold responses were around 0.6\*F\_MSY for most ecosystems. Compared with Fu et al. (2019), thresholds found in the present study under the F_all fishing strategy are generally similar, indicating general consistency between these two different approaches. However, by ignoring cases of linear tracking and monotonic non-linearity, the present study indicates that the occurrence of non-linearity is less frequent. In addition, thresholds identified using the fifth order polynomial curves are also infrequent for the non-linear-mixed type in most ecosystems under all fishing strategies (Figures 5–7; Supplementary Figures S1–S3). The elimination of threshold responses resulting from potentially linear and reversible tracking of perturbation is a step further towards ecosystem-based fisheries management by helping managers to identify sudden ecological transitions of real ecological concerns, because over-application of non-linearity and thresholds has caused confusion among managers (Litzow and Hunsicker, 2016).

Overall, it is becoming clear that tailoring decision-support systems to the unique biological and environmental characteristics of an ecosystem is critical, and, once again, we stress the importance of taking cognizance of the fisheries exploitation history of an ecosystem (e.g. Shannon et al., 2014; Coll et al., 2016; Lockerbie et al., 2017, 2018; Briton et al., 2019).

## Concluding remarks

Using GAMs combined with linear regression and polynomial models to analyse the outputs from the multi-ecosystem, multi-model simulation experiment conducted under the IndiSeas programme is useful for uncovering indicators' responses to fishing mortality in direction, non-linearity, and threshold. The present study has led to the following conclusions: (i) responses of indicators to fishing mortality in shape, direction, and threshold depend on the fishing strategies considered; (ii) the majority of the indicator response curves demonstrates negative relationships with fishing mortality, with a few exceptions which depend on fishing strategies, suggesting the importance of considering the history of fisheries exploitation when interpreting indicators for a particular marine ecosystem; (iii) most indicators tend to respond to fishing mortality in a linear way, particularly for community and biomass-based indicators, indicating that responses of these indicators to fishing mortality are more predictable; (iv) the infrequent occurrence of threshold for the non-linear-mixed response in most ecosystems may suggest that our approach has helped

eliminate threshold responses of linear and reversible tracking of perturbation and focus attention only on thresholds of real ecological concerns. Although future simulations should be carried out over a broader range of fishing mortality rate to further explore ecological thresholds, the various responses (linearity, non-linearity, and threshold) demonstrated here should be used to start to define management targets (levels and/or directions) and recommendations for sustainable levels of fishing pressure.

## Supplementary data

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

## Acknowledgements

This is a contribution to the IndiSeas Working Group, co-funded by IOC-UNESCO (www.ioc-unesco.org) and EuroMarine (http://www.euromarinenetwork.eu), to the project EMIBIOS (FRB, contract no. APP-SCEN-2010-II) and to the IOC-UNESCO GOOS Programme and GOOS Biology and Ecosystems Panel. The work on Canada West Coast ecosystem was sponsored by Fisheries & Oceans Canada under the Aquatic Climate Change and Adaptation Services Programme. YJS and LV were supported by the project EMIBIOS (FRB, contract no. APP-SCEN-2010-II). AG was supported by NOAA's Integrated Ecosystem Assessment (IEA) programme (http://www.noaa.gov/iea/). MC was supported by the Marie Curie Career Integration Grant Fellowships—PCIG10-GA-2011-303534—to the BIOWEB project and the Spanish National Project PELWEB (CTM2017-88939-R). All other authors were supported by their respective affiliations. We thank the participants of the IndiSeas Working Group for the useful discussions and insights since the early stages of the study. We are grateful to Huizhu Liu and Philippe Verley for their technical support on OSMOSE, Penny Johnson for her modelling work on Atlantis, Jeroen Steenbeek for his help with EwE, and Jennifer Houle for the size spectrum model. We are also thankful to Yu Bai for reviewing an earlier version of this manuscript and the R codes, as well as providing suggestions on threshold detection. At last but not the least, we would like to thank Editor Jason Link and three anonymous reviewers whose valuable comments have greatly improved the manuscript.

## References

* Akoglu, E. 2013. Nonlinear dynamics of the Black Sea ecosystem and its response to anthropogenic and climate variations. PhD thesis, Middle East Technical University, Ankara, Turkey.

* Akoglu, E., Salihoglu, B., Libralato, S., Oguz, T., and Solidoro, C. 2014. An indicator-based evaluation of Black Sea food web dynamics during 1960–2000. Journal of Marine Systems, 134: 113–125.

* Alexander, K. A., Heymans, J. J., MaGill, S., Tomczak, M., Holmes, S., and Wilding, T. A. 2015. Investigating the recent decline in gaidoid stocks in the west of Scotland shelf ecosystem using a food-web model. ICES Journal of Marine Science, 72: 436–449.

* Andersen, T., Carstensen, J., Hernandez, G. E., and Duarte, C. M. 2009. Ecological regime shifts: approaches to identification. Trends in Ecology and Evolution, 24: 49–57.

* Araújo, J. N., and Bundy, A. 2011. Description of the ecosystem models of the Bay of Fundy, Western Scotian Shelf and NAFO Division 4X. Canadian Technical Report of Fisheries and Aquatic Sciences, 2952: xii + 189 p.

* Araújo, J. N., and Bundy, A. 2012. The relative importance of climate change, exploitation and trophodynamic control in determining
ecosystem dynamics on the western Scotian Shelf. Canada. Marine Ecology Progress Series, 464: 51–67.

Baudron, A. R., Serpetti, N., Fallon, N. G., Heymans, J. J., and Fernandes, P. G. 2019. Can the common fisheries policy achieve good environmental status in exploited ecosystems: the west of Scotland demersal fisheries example. Fisheries Research, 211: 217–230.

Blamey, L. K., Howard, J. A. E., Agenbag, J., and Jarre, A. 2012. Regime-shifts in the southern Benguela shelf and inshore region. Progress in Oceanography, 106: 80–95.

Blanchard, J. L., Andersen, K. H., Scott, F., Hintzen, N. T., Piet, G., and Jennings, S. 2014. Evaluating targets and trade-offs among fisheries and conservation objectives using a multispecies size spectrum model. Journal of Applied Ecology, 51: 612–622.

Blanchard, J. L., Coll, M., Trenkel, V. M., Vergnon, R., Yemane, D., Jouffre, D., Link, J. S., and Shin, Y. J. 2010. Trend analysis of indicators: a comparison of recent changes in the status of marine ecosystems around the world. ICES Journal of Marine Science, 67: 732–744.

Boldt, J. L., Martone, R., Samhouri, J., Perry, R. I., Itoh, S., Chung, I. K., Takahashi, M. et al. 2014. Developing ecosystem indicators for responses to multiple stressors. Oceanography, 27: 116–133.

Boyce, D. G., Dowd, M., Lewis, M. R., and Worm, B. 2014. Estimating global chlorophyll changes over the past century. Progress in Oceanography, 122: 163–173.

Branch, T. A. 2015. Fishing impacts on food webs: multiple working hypotheses. Fisheries, 40: 373–375.

Briton, F., Shannon, L., Barrier, N., Verley, P., and Shin, Y. J. 2019. Reference levels of ecosystem indicators at multispecies maximum sustainable yield. ICES Journal of Marine Science, doi: 10.1093/icesjms/fsz104.

Brown, C. J., Schoeman, D. S., Sydeman, W. J., Brander, K., Buckley, L. B., Burrows, M., Duarte, C. M. et al. 2011. Quantitative approaches in climate change ecology. Global Change Biology, 17: 3697–3713.

Bundy, A., Shannon, L. J., Rochet, M. J., Neira, S., Shin, Y. J., Hill, L., and Aydin, K. 2010. The good (ish), the bad, and the ugly: a tripartite classification of ecosystem trends. ICES Journal of Marine Science, 67: 745–768.

Burthe, S. J., Henrys, P. A., Mackay, E. B., Spears, B. M., Campbell, R., Carvalho, L., Dudley, B. et al. 2016. Do early warning indicators consistently predict nonlinear change in long-term ecological data? Journal of Applied Ecology, 53: 666–676.

Capuzzo, E., Lynam, C. P., Barry, J., Stephens, D., Forster, R. M., Greenwood, N., McQuatters-Gollop, A. et al. 2018. A decline in primary production in the North Sea over twenty-five years, associated with reductions in zooplankton abundance and fish stock recruitment. Global Change Biology, 24: e352–e364.

Christensen, V., and Walters, C. 2004. Ecopath with Ecosim: methods, capabilities and limitations. Ecological Modelling, 72: 109–139.

Coll, M., Navarro, J., and Palomera, I. 2013. Ecological role of the endemic Starry ray Raja asterias in the NW Mediterranean Sea and management options for its conservation. Biological Conservation, 157: 108–120.

Coll, M., Palomera, I., Tudela, S., and Sardà, F. 2006. Trophic flows, ecosystem structure and fishing impacts in the south Catalan Sea, northwestern Mediterranean. Journal of Marine Systems, 59: 63–96.

Coll, M., Shannon, L. J., Kleisner, K. M., Juan Jordà, M. J., Bundy, A., Akoglu, G., Banaru, D. et al. 2016. Ecological indicators to capture the effects of fishing on biodiversity and conservation status of marine ecosystems. Ecological Indicators, 60: 947–962.

Cury, P., Bakun, A., Crawford, R. J., Jarre, A., Quinones, R. A., Shannon, L. J., and Verheye, H. M. 2000. Small pelagics in upwell-ing systems: patterns of interaction and structural changes in

wasp-waist ecosystems. ICES Journal of Marine Science, 57: 603–618.

Fernandes, P. G., Ralph, G. M., Nieto, A., Criado, M. G., Vasilakopoulos, P., Maravelias, C. D., Cook, R. M. et al. 2017. Coherent assessments of Europe's marine fishes show regional divergence and megafauna loss. Nature Ecology & Evolution, 1: 0170.

Foley, M. M., Martone, R. G., Fox, M. D., Kappel, C. V., Mease, L. A., Erickson, A. L., Halpern, B. S. et al. 2015. Using ecological thresholds to inform resource management: current options and future possibilities. Frontiers in Marine Science, 2: 95.

Frank, K. T., Petrie, B., and Shackell, N. L. 2007. The ups and downs of trophic control in continental shelf ecosystems. Trends in Ecology & Evolution, 22: 236–242.

Fu, C., Gaichas, S., Link, J. S., Bundy, A., Boldt, J. L., Cook, A. M., Gamble, R. et al. 2012. Relative importance of fishing, trophodynamic and environmental drivers in a series of marine ecosystems. Marine Ecology Progress Series, 459: 169–184.

Fu, C., Large, S., Knight, B., Richardson, A., Bundy, A., Reygondeau, G., and Boldt, J. 2015. Relationships among fisheries exploitation, environmental conditions, and ecological indicators across a series of marine ecosystems. Journal of Marine Systems, 148: 101–111.

Fu, C., Perry, R. I., Shin, Y.-J., Schweigert, J., and Liu, H. 2013. An ecosystem modelling framework for incorporating climate regime shifts into fisheries management. Progress in Oceanography, 115: 53–64.

Fu, C., Travers-Trolet, M., Velez, L., Grüss, A., Bundy, A., Shannon, L. J., Fulton, E. A. et al. 2018. Risky business: the combined effects of fishing and changes in primary productivity on fish communities. Ecological Modelling, 368: 265–276.

Fu, C., Xu, Y., Bundy, A., Grüss, A., Coll, M., Heymans, J. J., Fulton, E. A. et al. 2019. Making ecological indicators management ready: assessing their ability to detect impacts of fishing and environmental change. Ecological Indicators, 105: 16–28.

Fulton, E. A., Parslow, J. S., Smith, A. D., and Johnson, C. R. 2004. Biogeochemical marine ecosystem models II: the effect of physiological detail on model performance. Ecological Modelling, 173: 371–406.

Fulton, E. A., Smith, A. D. M., Smith, D. C., and Johnson, P. 2014. An integrated approach is needed for ecosystem based fisheries management: insights from ecosystem-level management strategy evaluation. PLoS One, 9: e84242.

Gabriel, W. L., and Mace, P. M. 1999. A review of biological reference points in the context of the precautionary approach. NOAA Technical Memorandum NMFS-F/SPO-40. pp. 34–45.

Groffman, P. M., Baron, J. S., Blett, T., Gold, A. J., Goodman, I., Gunderson, L. H., Levinson, B. M. et al. 2006. Ecological thresholds: the key to successful environmental management or an important concept with no practical application? Ecosystems, 9: 1–13.

Grüss, A., Harford, W. J., Schirripa, M. J., Velez, L., Sagarese, S. R., Shin, Y.-J., and Verley, P. 2016b. Management strategy evaluation using the individual-based, multispecies modeling approach OSMOSE. Ecological Modelling, 340: 86–105.

Grüss, A., Schirripa, M. J., Chagaris, D., Velez, L., Shin, Y.-J., Verley, P., and Oliveros-Ramos, R. 2016a. Evaluating natural mortality rates and simulating fishing scenarios for Gulf of Mexico red grouper (*Epinephelus morio*) using the ecosystem model OSMOSE-WFS. Journal of Marine Systems, 154: 264–279.

Halouani, G., Lasram, F. B. R., Shin, Y.-J., Velez, L., Verley, P., Hattab, T., Oliveros-Ramos, R. et al. 2016. Modelling food web structure using an end-to-end approach in the coastal ecosystem of the Gulf of Gabes (Tunisia). Ecological Modelling, 339: 45–57.

Hastie, T., and Tibshirani, R. 1990. Generalized Additive Models. Chapman and Hall/CRC Press, London, UK.

Heymans, J. J., and Tomczak, M. T. 2016. Regime shifts in the northern Benguela ecosystem: challenges for management, in: Ecopath
30 years—modelling ecosystem dynamics: beyond boundaries with EwE. *Ecological Modelling*, 331: 151–159.

Hunt, G., Coyle, K. O., Eisner, L. B., Farley, E. V., Heintz, R. A., Mueter, F., Napp, J. M. et al. 2011. Climate impacts on eastern Bering Sea foodwebs: a synthesis of new data and an assessment of the oscillating control hypothesis. *ICES Journal of Marine Science*, 68: 1230–1243.

Jackson, J. B. C., Kirby, M. X., Berger, W. H., Bjorndal, K. A., Botsford, L. W., Bourque, B. J., Bradbury, R. H. et al. 2001. Historical overfishing and the recent collapse of coastal ecosystems. *Science*, 293: 629–638.

Kelly, R. P., Erickson, A. L., Mease, L. A., Battista, W., Kittinger, J. N., and Fujita, R. 2015. Embracing thresholds for better environmental management. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 370: 20130276.

Kirby, R., Beaugrand, G., and Lindley, J. 2009. Synergistic effects of climate and fishing in a marine ecosystem. *Ecosystems*, 12: 548–561.

Large, S. I., Fay, G., Friedland, K. D., and Link, J. S. 2013. Defining trends and thresholds in responses of ecological indicators to fishing and environmental pressures. *ICES Journal of Marine Science*, 70: 755–767.

Large, S. I., Fay, G., Friedland, K. D., and Link, J. S. 2015. Critical points in ecosystem responses to fishing and environmental pressures. *Marine Ecology Progress Series*, 521: 1–17.

Link, J. S., Yemane, D., Shannon, L. J., Coll, M., Shin, Y-J., Hill, L., and Borges, M. F. 2010. Relating marine ecosystem indicators to fishing and environmental drivers: an elucidation of contrasting responses. *ICES Journal of Marine Science*, 67: 787–795.

Lintz, H. E., McCune, B., Gray, A. N., and McCulloh, K. A. 2011. Quantifying ecological thresholds from response surfaces. *Ecological Modelling*, 222: 427–436.

Litzow, M. A., and Hunsicker, M. E. 2016. Early warning signals, nonlinearity, and signs of hysteresis in real ecosystems. *Ecosphere*, 7: e01614-10.

Lockerbie, E., Shannon, L. J., and Jarre, A. 2016. The use of ecological, fishing and environmental indicators in support of decision making in southern Benguela fisheries. *Ecological Modelling*, 69: 473–487.

Lockerbie, E. M., Coll, M., Shannon, L. J., and Jarre, A. 2017. The use of indicators for decision support in northwestern Mediterranean Sea fisheries. *Journal of Marine Systems*, 174: 64–77.

Lockerbie, E. M., Lynam, C. P., Shannon, L. J., and Jarre, A. 2018. Applying a decision tree framework in support of an ecosystem approach to fisheries: IndiSeas indicators in the North Sea. *ICES Journal of Marine Science*, 75: 1009–1020.

Möllmann, C., Diekmann, R., Müller-Karulis, B., Kornilovs, G., Plikshs, M., and Axe, P. 2009. Reorganization of a large marine ecosystem due to atmospheric and anthropogenic pressure: a discontinuous regime shift in the central Baltic Sea. *Global Change Biology*, 15: 1377–1393.

Moore, J. C. 2018. Predicting tipping points in complex environmental systems. *Proceedings of the National Academy of Sciences of the United States of America*, 115: 635–636.

Nicolaou, T., and Constandinou, G. 2016. A nonlinear causality estimator based on non-parametric multiplicative regression. *Frontiers in Neuroinformatics*, 10: 19.

Ortega-Cisneros, K., Cochrane, K. L., Fulton, E. A., Gorton, R., and Popova, E. 2018. Evaluating the effects of climate change in the southern Benguela upwelling system using the Atlantis modelling framework. *Fisheries Oceanography*, 27: 489–503.

Pauly, D., and Watson, R. 2005. Background and interpretation of the 'Marine Trophic Index' as a measure of biodiversity. *Philosophical Transactions of the Royal Society B*, 360: 415–423.

Poloczanska, E. S., Burrows, M. T., Brown, C. J., Garcia Molinos, J., Halpern, B. S., Hoegh-Guldberg, O., Kappel, C. V. et al. 2016. Responses of marine organisms to climate change across oceans. *Frontiers in Marine Science*, 3: 62.

R Core Team, 2019. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

Reed, J., Shannon, L., Velez, L., Akoglu, E., Bundy, A., Coll, M., Fu, C. et al. 2017. Ecosystem indicators—accounting for variability in species’ trophic level. *ICES Journal of Marine Science*, 74: 158–169.

Rice, J. C., and Rochet, M-J. 2005. A framework for selecting a suite of indicators for fisheries management. *ICES Journal of Marine Science*, 62: 516–527.

Samhoury, J. F., Andrews, K. S., Fay, G., Harvey, C. J., Hazen, E. L., Hennessey, S. M., Holsman, K. et al. 2017. Defining ecosystem thresholds for human activities and environmental pressures in the California Current. *Ecosphere*, 8: e01860.

Scheffer, M., and Carpenter, S. R. 2003. Catastrophic regime shifts in ecosystems: linking theory to observation. *Trends in Ecology and Evolution*, 18: 648–656.

Scheffer, M., Carpenter, S., Foley, J. A., Folke, C., and Walker, B. 2001. Catastrophic shifts in ecosystems. *Nature*, 413: 591–596.

Serpetti, N., Baudron, A. R., Burrows, M. T., Payne, B. L., Helaouët, P., Fernandes, P. G., and Heymans, J. J. 2017. Impact of ocean warming on sustainable fisheries management informs the ecosystem approach to fisheries. *Scientific Reports*, 7: 13438.

Shannon, L., Coll, M., Bundy, A., Gascuel, D., Heymans, J. J., Kleisner, K., Lynam, C. P. et al. 2014. Trophic level-based indicators to track fishing impacts across marine ecosystems. *Marine Ecology Progress Series*, 512: 115–140.

Shannon, L. J., Christensen, V., and Walters, C. 2004. Modelling stock dynamics in the southern Benguela ecosystem for the period 1978–2002. *African Journal of Marine Science*, 26: 179–196.

Shannon, L. J., Coll, M., Yemane, D., Jouffre, D., Neira, S., Bertrand, A., Diaz, E. et al. 2010. Comparing data-based indicators across upwelling and comparable systems for communicating ecosystem states and trends. *ICES Journal of Marine Science*, 67: 807–832.

Shannon, L. J., Neira, S., and Taylor, M. 2008. Comparing internal and external drivers in the southern Benguela and the southern and northern Humboldt upwelling ecosystems. *African Journal of Marine Science*, 30: 63–84.

Shin, Y-J., Bundy, A., Shannon, L. J., Blanchard, J. L., Chuenpagdee, R., Coll, M., Knight, B. et al.; the IndiSeas Working Group 2012. Global in scope and regionally rich: an IndiSeas workshop helps shape the future of marine ecosystem indicators. *Reviews in Fish Biology and Fisheries*, 22: 835–845.

Shin, Y-J., and Cury, P. 2004. Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing. *Canadian Journal of Fisheries and Aquatic Sciences*, 61: 414–431.

Shin, Y-J., Houle, J. E., Akoglu, E., Blanchard, J., Bundy, A., Coll, M., Demarcq, H. et al. 2018. The specificity of marine ecological indicators to fishing in the face of environmental change: a multi-model evaluation. *Ecological Indicators*, 89: 317–326.

Shin, Y-J., Shannon, L. J., Bundy, A., Coll, M., Aydin, K., Bez, N., Blanchard, J. L. et al. 2010. Using indicators for evaluating, comparing, and communicating the ecological status of exploited marine ecosystems. 2. Setting the scene. *ICES Journal of Marine Science*, 67: 692–716.

Tam, J. C., Link, J. S., Large, S. I., Andrews, K., Friedland, K. D., Gove, J., Hazen, E. et al. 2017. Comparing apples to oranges: common trends and thresholds in anthropogenic and environmental pressures across multiple marine ecosystems. *Frontier in Marine Science*, 4: 282.

Wood, S. N. 2017. Generalized Additive Models: An Introduction with R, 2nd edn. Chapman and Hall/CRC Press, London, UK.

Wood, S. N., Pya, N., and Säfken, B. 2016. Smoothing parameter and model selection for general smooth models (with discussion). *Journal of the American Statistical Association*, 111: 1548–1575.