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

Responses of ecological indicators to ﬁshing pressure under environmental change: exploring non-linearity and Marine ecosystems are inﬂuenced by multiple stressors in both linear and non-linear ways. Using generalized additive models (GAMs) ﬁtted to outputs from a multi-ecosystem, multi-model simulation experiment, we investigated 14 major ecological indicators across ten marine ecosystems about their responses to ﬁshing pressure under: (i) three different ﬁshing 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 ﬁshing mortality in shapes, directions, and thresholds depend on the ﬁshing strategies considered; (ii) most of the indicators demonstrate decreasing trends with increasing ﬁshing mortality, with a few exceptions depending on the type of ﬁshing strategy; (iii) most of the indicators respond to ﬁshing 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 inﬂection points) is not prevalent within the ﬁshing mortality rates explored. The conclusions drawn from the present study provide a knowledge base in indicators’ dynamics under different ﬁshing and primary productivity levels, thereby facilitating the application of ecosystem-based ﬁsheries management worldwide.


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 socioeconomic 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. nonlinearly) 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., 2013Large et al., , 2015Tam 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., 2018Fu et al., , 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.

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(Shannon et al., , 2008, the southern Catalan Sea (Coll et al., 2006(Coll et al., , 2013, the western Scotian Shelf (Araú jo and Bundy, 2011, 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(Grüss et al., , 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 k. In other words, fishing mortality rates ) were applied to each species/taxon of interest i, where k 2 {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 lowtrophic 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, timeseries 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 c 2 {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 l of 1 and a range of standard deviations r ¼ {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 r, 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 r ¼ 0.1, 0.2, or 0.3) were simulated to explore interactions of these two factors.  (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 community-level biomass indicators, the indicators Pred, Lifesp, TLco were also derived from survey biomass data (i.e. biomassbased); 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): where Y is the ecological indicator; a is the intercept held constant; X is the predictor, representing either fishing mortality or primary productivity level; sðÞ is the smoothing function; and e 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 r ¼ 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 deviance 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 timeseries 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 Responses of ecological indicators to fishing pressure under environmental change 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: 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 ; Á Á Á ; a 5 are coefficients of the polynomial model. The analytical root(s) of the second derivatives, where ecological thresholds reside, are derived from: 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.

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    Table 1. 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 r 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    Table 1.
Responses of ecological indicators to fishing pressure under environmental change (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      Table 1.
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-linearmixed 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 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 r ¼ 0.1, r ¼ 0.2, and r ¼ 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.

Responses of ecological indicators to fishing pressure under environmental change
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 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 (F_htl) fishing strategy using fishing mortality as the predictor variable under 4 different scenarios of primary productivity (directional, random changes with r ¼ 0.1, r ¼ 0.2, and r ¼ 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.

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) 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 r ¼ 0.1, r ¼ 0.2, and r ¼ 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.
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  Table 1. 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. 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.
( 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 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 r ¼ 0.1, 0.2, or 0.3), which suggests that the indicators considered in the present study are robust to environmental variability for the Black Sea ecosystem.
The   Table 1.
Responses of ecological indicators to fishing pressure under environmental change indicators. These results concur with the findings of Grüss et al. (2016aGrüss et al. ( , 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 P 2 values compared with R F 2 , 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 . All indicators except for Lifesp, TLc, TLcVar, and MTI have very low R F 2 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 P 2 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 P 2 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 P 2 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 demonstrates 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 nonlinearity 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., 2013Large et al., , 2015Samhouri 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., 2017Lockerbie et al., , 2018Briton et al., 2019).

Concluding remarks
Using GAMs combined with linear regression and polynomial models to analyse the outputs from the multi-ecosystem, multimodel 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.