Functional group based marine ecosystem assessment for the Bay of Biscay via elasticity analysis

Observation model

Conditional on the biomass x_i(t) of the ith functional group in year t, the observation model for the French commercial fishery and fishery independent survey CPUE data is given by, (1)${\displaystyle p\left(C_{ik}\left(t\right)|x_i\left(t\right),E_k\left(t\right),q_{ik},{\sigma }_k\right)=LN\left(C_{ik}\left(t\right)|log\left(q_{ik}{E}_k\left(t\right)x_i\left(t\right)\right),{\sigma }_k^2\right),}$

where LN(⋅|m, s²) is a log-normal distribution with mean m and variance s² on the log-scale. C_ik(t) is the catch in metric tonnes for functional group $i\in \left\{K,B,D,G\right\}$ in year t by gear type $k\in \left\{Dredge,Hook,\right.\left.Mixedtrawl,Net,Pelagictrawl,Pot,Purseseine,Survey\right\}$with catchability q_ik and effort E_k(t). The observation error is given by ${\omega }_{ik}\left(t\right)\stackrel{i.i.d.}{\sim }N\left(0,{\sigma }_k^2\right)$ for $i\in \left\{K,B,D,G\right\}$ and the kth gear type. Equation (1) applies to the targeted combinations of gear type described above and as shown in Fig. 1. The catches from incidental combinations of gear type and functional group are assumed ignorable and not included in the likelihood.

The impact of catch from French landings on the unobserved latent state of biomass for the ith functional group in year t is captured by the international landings h_i(t). The catch component C_ik(t) of the observed CPUE from French fisheries thus explicitly enters only the observation model of Eq. (1). With no error in international reported landings, the fished fraction for the ith functional group in year t is given by, (2)${\displaystyle F_i\left(t\right)=\frac{h_i\left(t\right)}{x_i\left(t\right)}=1-e^{-f_i\left(t\right)},}$

where f_i(t) is the instantaneous rate of fishing mortality on the ith functional group in year t. However, it is possible for errors or misreporting of landings to occur. Uncertainty was allowed for by specifying a lognormal distribution for annual international landings such that h_i(t) ∼ LN(log F_i(t) + log x_i(t), 0.001), which approximately provides probability 0.9 of less than 5% reporting error.

The scientific survey data for the relevant functional groups (i.e., excluding pelagic piscivores) is similarly addressed as with the commercial CPUE data, but with different units: average metric tonnes of biomass per square kilometre for each functional group by survey year. The observation error has independent unknown variance, which reflects the fact that the observations from the scientific surveys align to an unknown degree with the fishable biomass of a given functional group.

Process model

Fox (1970) introduced the continuous-time Gompertz-Fox model for a fish population x(t) given constant effort E and catchability q with the ordinary differential equation (3)${\displaystyle \frac{dx}{dt}=bxlog\frac{\beta }{x}-fx,}$

where b is related to the density independent growth rate at a particular population level (see Thompson, 1996), β is the carrying capacity of the environment and f = qE is the instantaneous fishing mortality rate. The first term of Eq. (3) is the Gompertz form of density dependence (Thieme, 2003). Extension of Eq. (3) to multiple variables changes the interpretation of the Gompertz-Fox model parameters because the carrying capacity for the ith variable, which is simply β in Eq. (3), may now depend on other variables. Consider the multivariate extension of the Gompertz-Fox model to the n-dimensional system (4)${\displaystyle \frac{dx}{dt}=\text{diag}\left[x\right]\left(r-f+Alogx\right),}$

where f is the n-dimensional vector of instantaneous fishing mortality rates and A is the n × n matrix of density-dependent parameters. The n-dimensional vector of density independent growth rates r gives the per capita growth rates in the absence of fishing or density dependence. If the variables of the above equation are independent such that A is a diagonal matrix, then the substitutions r_i = b_ilog β_i and a_ii =  − b_i recover the parameterisation of Eq. (3) for variable x_i.

Different choices of discretisation methods applied to Eq. (3) will lead to different discrete-time models (e.g., Turchin, 2003), and so also for multivariate models. A deterministic multivariate discrete-time version of Eq. (4) is given by (5)${\displaystyle x\left(t+1\right)=\text{diag}\left[x\left(t\right)\right]exp\left[r-f+Alogx\left(t\right)\right],}$

where, as for Eq. (4), r is the vector of density independent growth parameters, f is the vector of instantaneous fishing mortalities and the matrix A determines the strength of density-dependence within and among variables.

The discrete-time model of Eq. (5) shares a key property associated with the long term behaviour of the continuous-time model given by Eq. (4). For small time steps in Eq. (5), the biomass of the ith variable at time t + τ is ${\displaystyle x_i\left(t+\tau \right)=x_i\left(t\right)exp\left[\tau \left(r_i-f_i+\sum _j{A}_{ij}log{x}_j\left(t\right)\right)\right].}$

As the time step τ approaches zero, ${\displaystyle \frac{d{x}_i}{dt}=\underset{\tau \to 0}{lim}\frac{x_i\left(t+\tau \right)-x_i\left(t\right)}{\tau }=x_i\left(t\right)\underset{\tau \to 0}{lim}\frac{exp{\left(r_i-f_i+\sum _j{A}_{ij}log{x}_j\left(t\right)\right)}^{\tau }-1}{\tau }}$ (6)${\displaystyle =x_i\left(t\right)\left[r_i-f_i+\sum _j{A}_{ij}log{x}_j\left(t\right)\right].}$ Equation (6) is equivalent to the continuous-time model given by Eq. (4). Note that if r_i − f_i + ∑_jA_ijlog x_j(t) = 0 then both dx_i∕dt = 0 in Eq. (4) and x_i(t + τ)∕x_i(t) = 1 in Eq. (5). The continuous-time model of Eq. (4) and the discrete-time model of Eq. (5) therefore share the same interior equilibrium given by $x^{\ast }=exp\left[-A^{-1}\left(r-f\right)\right]$. In this sense, Eq. (5) is an interpretable discretisation of the Fox model extended to a multidimensional system because its parameters have the same functions as those of Eq. (4).

Elasticities for one year look-ahead predictions

In economics, elasticity is defined as the proportional response of a variable of interest to a proportional increase of another variable (Hicks, 1946, Ch. 16). Elasticities are also important in population biology (Caswell, 2001). Elasticities are useful because the focus on proportional change results in a non-dimensional analysis, which provides a common scale to support comparison of how the variable of interest responds to changes in parameters having different units. A fundamental quantity of interest is how the biomass of the ith functional group x_i(t + 1) proportionally responds to a proportional change in a variable v(t) during the preceding year t. This quantity is the elasticity for a one-year look-ahead prediction given by ${\displaystyle \frac{\partial x_i\left(t+1\right)∕x_i\left(t+1\right)}{\partial v\left(t\right)∕v\left(t\right)}=\frac{\partial log{x}_i\left(t+1\right)}{\partial logv\left(t\right)},}$

where x_i(t + 1) is the growth equation defined by Eq. (9). The above equation is an example of a transient elasticity (Caswell, 2007), here applied to a stochastic dynamic model. For a small perturbation of v(t), an elasticity of +1 predicts that the next year’s biomass will increase by the same percentage as the percentage increase in v(t). If the elasticity is greater than +1 then the next year’s biomass will increase by a greater percentage, and if the elasticity is less than +1 (but greater than zero) then the next year’s biomass will increase by a lesser percentage. Negative elasticities are similarly interpreted but instead predict decreases in the next year’s biomass given a small increase to v(t).

From Eqs. (8) and (9), the logarithm of the growth equation for functional group i from year t to t + 1 is given by ${\displaystyle log{x}_i\left(t+1\right)=\left(1+A_{ii}\right)log{x}_i\left(t\right)+\sum _{j\ne i}{A}_{ij}log{x}_j\left(t\right)+r_i-f_i\left(t\right)+{ϵ}_i\left(t\right),}$

where A_ij accounts for both the sign and the magnitude a_ij of the density dependent effect of functional group j on functional group i. The elasticity of the one-year look-ahead prediction for the biomass of the ith functional group in year t + 1 given a change to the biomass of the jth functional group in year t is then (10)${\displaystyle \frac{\partial log{x}_i\left(t+1\right)}{\partial log{x}_j\left(t\right)}=\left\{\begin{array}{cc}1+A_{ii}=1-a_{ii}\hfill & ifi=j\hfill \\ A_{ij}=\text{sgn}\left(A_{ij}\right)a_{ij}\hfill & ifi\ne j,\hfill \end{array}\right.}$

where sgn(z) is +1, 0 or −1 depending on whether the argument z is positive, zero or negative. Note that the density dependent elasticities for the one-year look-ahead predictions are invariant over time.

Elasticities are defined for positive quantities. The density independent growth rate r_i may however be negative, in which case it instead describes the rate of decay. For example, the growth of a consumer functional group may solely depend on the biomass of a producer functional group. Without its prey or any alternative sources of growth or mortality aside from density independence, then only a proportion exp[r_i] < 1 of the consumer biomass will survive into the next year. To allow for negative r_i, the elasticity of biomass to a change in the magnitude |r_i| is instead evaluated by substituting sgn(r_i)|r_i| for r_i in the logarithm of the growth equation log x_i(t + 1). The elasticity of the one-year look-ahead prediction for a change in the magnitude of a non-zero density independent growth (decay) rate r_i is then given by (11)${\displaystyle \frac{\partial log{x}_i\left(t+1\right)}{\partial log\left|r_i\right|}=\frac{\partial x_i\left(t+1\right)}{\partial \left|r_i\right|}\frac{\left|r_i\right|}{x_i\left(t+1\right)}=x_i\left(t+1\right)\text{sgn}\left(r_i\right)\frac{\left|r_i\right|}{x_i\left(t+1\right)}=\text{sgn}\left(r_i\right)\left|r_i\right|=r_i,}$

which is invariant over time. In contrast, the elasticities of the one-year look-head predictions for the instantaneous fishing mortality rates vary by year, ${\displaystyle \frac{\partial log{x}_i\left(t+1\right)}{\partial log{f}_i\left(t\right)}=-f_i\left(t\right),}$ where f_i(t) > 0.

Long term dynamics

Consider the long-term implications of a management policy that holds the instantaneous fishing mortality constant for each functional group such that f(t) = f for all t. This policy is equivalent to specifying a constant annual fished fraction of biomass, see Eq. (2). The logarithmically transformed dynamic system of Eq. (9) with time-invariant fishing mortalities is given by, (11)${\displaystyle logx\left(t+1\right)=\left(I+A\right)logx\left(t\right)+r-f+ϵ\left(t\right).}$

A dynamic system is strictly stationary if realisations of the process for a set of times $\tau =\left\{t_1,t_2,\dots ,t_T\right\}$ have the same joint probability distribution as realisations at times τ + s for constant s. The dynamic system of Eq. (11) is strictly stationary if the spectral radius of the matrix I + A is less than one. Given stationarity, the joint distribution for log x(t) is multivariate normal with location vector m =  − A⁻¹(r − f) and covariance matrix P that solves P = (I + A)P(I + A)^⊤ + S (Harvey, 1989, Ch. 3). The distribution for x(t) is therefore multivariate lognormal, x(t) ∼ LN(m, P), such that $\hat{x}=exp\left[m\right]$ is the vector of long-term median biomasses of the functional groups and $\bar{x}=exp\left\{m+\text{diag}\left[P\right]∕2\right\}$ gives the long-term average biomasses.

Cumulative impacts of fishing pressure

Whereas the above one-year look-ahead elasticities address transient dynamics, it is also of interest how the system responds in the long-term to a scenario with constant fishing mortality rates. If the density dependence matrix A permits a stationary system, then the elasticities of certain attributes of the stationary distribution may be investigated. In particular, it will be of interest to investigate how the stationary distribution responds to a proportional change in the unfished fraction of a functional group, U_i = 1 − F_i. From Eq. (2), note that dU_i∕(U_idf_i) = dlogU_i∕df_i =  − 1. Therefore a small perturbation of size Δ in the value of f_i results in a −Δ proportional change in the unfished fraction U_i, which provides a useful interpretation of the fishing mortality rate.

Let $\tilde{x}\left(\theta ,y\right)=exp\left[-A^{-1}\left(r-f\right)+y\left(A,S\right)\right]$ denote attributes of the long-term stationary distribution of x(t), where y(A, S) is a function of the density dependent parameters A and the process uncertainty matrix S. The vector $\tilde{x}\left(\theta ,y\right)$ is equivalent to the long-term averages of the functional groups if y(A, S) = diag[P]∕2, and is equivalent to the long-term medians if y(A, s) = 0. In fact, due to the properties of the lognormal distribution (see, e.g., Aitchison & Brown, 1957, Ch. 2), the marginal p-quantiles of the stationary distribution of x(t) may for any p ∈ (0, 1) be obtained by choosing y(A, S) = diag[P]^1∕2ν_p, where ν_p is the p-quantile of a standard normal distribution. The elasticities of $\tilde{x}\left(\theta ,y\right)$ with respect to the unfished fractions are then given by, (12)${\displaystyle \frac{\partial log\tilde{x}\left(\theta ,y\right)}{\partial logU}=\frac{\partial log\tilde{x}\left(\theta ,y\right)}{\partial f}\frac{\partial f}{\partial logU}=-A^{-1},}$

where the proportional changes in $\tilde{x}\left(\theta ,y\right)$ given perturbed fishing mortality rates are $\partial log\tilde{x}\left(\theta ,y\right)∕\partial f=A^{-1}$ and ∂f∕∂log U =  − I. The elasticity matrix of Eq. (12) applies to the long-term means or the marginal p-quantiles, including the medians, of x(t) presuming that the long-term distribution is stationary. The elasticity of functional group i with respect to the unfished fraction of functional group j is given by the (i, j) entry of −A⁻¹. The elasticities of all functional groups with respect to the unfished biomass of the jth functional group are given by the jth column of the matrix −A⁻¹. This corresponds to a scenario where a fishery or set of fisheries (gear types) decrease aggregate fishing pressure on a particular functional group.

On the other hand, the scenario where the fishing pressure of a particular fishery or gear type is restricted may be of interest. In the Bay of Biscay, a particular gear type often targets more than one functional group (Fig. 1). Since elasticities can be added (Caswell, 2001, Ch. 9), the elasticities to simultaneous proportional increases of equal magnitudes in the unfished fractions of more than one functional groups are obtained by summing the columns of −A⁻¹ that correspond to the functional groups with decreased fishing pressure. Note also that an increase in fishing pressure corresponds to decreased unfished fractions, which simply reverses the sign of Eq. (12). Correspondingly, the response of a functional group to increased fishing pressure in one functional group and decreased fishing pressure in another is given by the difference of the respective columns of −A⁻¹. Therefore both the direct and indirect long-term effects of fishing pressure on functional group biomasses due to changes in either single or multiple fisheries may be predicted from Eq. (12).

Four long-term management scenarios are explored, where the long-term stationary distribution is assessed for four choices of constant fishing mortality rates:

f⁰

This scenario investigates the long-term stationary distribution given cessation of all fishing.

f^Avg

This scenario sets the fishing mortality rates to the estimated average over the years 1999–2015.

f^G

This scenario is as for f^Avg except that the fishing mortality rates for pelagic piscivores G is set to zero. This scenario examines the amount of compensation that might occur in the Bay of Biscay ecosystem from increased top-down control by pelagic piscivores released from fishing pressure.

f^MT

This scenario is as for f^Avg except that the fishing mortality rates are reduced proportional to the reduction in the average landings by French fisheries for ICES divisions 8a and 8b of a functional group from 1999 to 2015 with landings from mixed trawls excluded. The average proportion of benthivore landings attributed to mixed trawls was 62%, and the average proportion of demersal piscivore landings attributed to mixed trawls was 45%. The fishing mortality rates for benthivores B and demersal piscivores D are reduced to 38% and 55% of their f^Avg values, respectively. The fishing mortality rates for pelagic planktivores and piscivores are as for f^Avg because mixed trawls did not substantially contribute to landings for these latter two functional groups.

Priors

Table 1 summarises the prior information for the unknown parameters, concatenated into the vector of unknown parameters θ, which includes the latent states x(t = 1999), to form the prior density p(θ). A description of the information contained by each prior distribution is given below.

Posterior density

The posterior distribution is defined by the process model of Eq. (9) for the latent states x(t) in the years t = 2000, …, 2015, by the observation model of Eq. (1) and by the prior distributions for the parameters, θ, that include the unknown initial latent states, x(t = 1999). The posterior density is given by, ${\displaystyle p\left(x\left(1999:2000\right),\theta |C,E,h\left(1999:2015\right)\right)\propto p\left(\theta \right)p\left(h\left(1999\right)|x\left(1999\right),F\left(1999\right)\right)}$ (13)${\displaystyle \times \prod _{t=2000}^{2015}p\left(x\left(t\right)|x\left(t-1\right),\theta \right)\left(\prod _{i}p\left(h_i\left(t\right)|x_i\left(t\right),F_i\left(t\right)\right)\left[\prod _{k}p\left(C_{ik}\left(t\right)|x_i\left(t\right),E_k\left(t\right),\theta \right)\right]\right),}$

where x(1999:2000) gives the latent biomasses of the functional groups for the years 1999–2000, C is the array of annual catches by French commercial fisheries and scientific surveys that target specific functional groups for years 2000–2015, E is the array of effort recorded by French commercial fisheries and scientific surveys for years 2000–2015, and h(1999:2015) gives the recorded international landings for each functional group for years 1999–2015. The first line of Eq. (13) incorporates the priors for the parameters θ and the likelihood of the observed international landings in 1999 given the initial states. The second line of Eq. (13) incorporates the prior for the latent dynamic process of the unknown annual biomasses of the functional groups for years 2000–2015. The second line also includes the likelihood for the years 2000–2015 for harvest of functional groups $i\in \left\{K,B,D,G\right\}$ and gear type $k\in \left\{Dredge,Hook,Mixedtrawl,Net,Pelagictrawl,Pot,Purseseine,Survey\right\}$. The C_ik(t) are assumed non-negligible only for the targeted combinations of gear type and functional group shown in Fig. 1; catches for incidental combinations of gear type and functional group are assumed ignorable and do not contribute to the likelihood. The posterior distribution of a function g(θ) of the parameters θ, such as the long-term stationary distribution of the latent states or the elasticities, is given by p(g(θ)|C, E, h).

Estimation

Simulations and posterior inferences used R (R Core Team, 2018) with posterior samples obtained by Markov chain Monte Carlo (MCMC; see e.g., Robert & Casella, 2004; Gelman et al., 2014) using JAGS (Plummer, 2003) from R port rjags (Plummer, 2016). Five independent parallel chains were initialised by drawing realisations from the prior distributions except for the log catchabilities and the latent states. The log catchabilities were initialised to the average of the log CPUE indices. The latent state initialisations must not conflict with the censored observations, that is, the initial latent state trajectory must allow for a feasible harvest sequence. The latent states were therefore initialised from a Pareto distribution with cumulative distribution function given by P(x_i(t = 1999) < w) = 1 − (h_i(t)∕w)² if w ≥ h_i(t) and 0 otherwise. A priori there is probability 0.75 that the fished fraction of functional group i in year t was less than 0.5. The Pareto prior is heavy-tailed and, with index parameter 2, has finite mean that is twice the observed value of landings and infinite variance.

The five independent parallel chains were each run for 500,000 adaptive iterations followed by 500,000 non-adaptive iterations with thinning rate of 10 applied to the non-adaptive phases. Trace plots were examined and Gelman and Rubin’s convergence diagnostic (Gelman & Rubin, 1992) was applied using the coda package (Plummer et al., 2006). The upper bounds of the 95% confidence intervals for the estimated potential scale reduction factors were near 1 for every parameter. Results presented use the 250,000 posterior samples retained from the thinned non-adaptive phase. The code and data used for the estimation of the model are included as an R package (see Supplemental Information).

Observation model

The observation model uncertainty provided intuitively reasonable results (Fig. 2). As noted in the prior description, an observation standard deviation of 1 corresponds to a 50/50 chance that the observed CPUE in year t is within one-half or 2 × the true biomass. A standard deviation of 0.3 gives 0.5 probability that the CPUE is within ±20%, and a standard deviation of 0.03 give 0.5 probability of being within a ±2% margin. The observation error standard deviation estimates were highest (had greatest uncertainty) for purse seines and were lowest for mixed trawls (Fig. 2). Thus the trawl fishery yielded greater precision than the purse seiners. The related parameter of catchability was lowest when looking at rare functional group and gear combinations. For example, estimated catchability was higher for both benthivores and demersal piscivores than pelagic piscivores in the netter fishery (Table S2). On the other hand, catchability for pelagic planktivores was higher than demersal piscivores in pelagic trawls. Estimated catchability was higher for planktivores then pelagic piscivores in the purse seine fishery.

Ecosystem dynamics

There are two main conclusions for latent biomass states (Fig. 3). First, there were common qualitative trends evident for benthivores and both demersal and pelagic piscivores. These functional groups appear to have declined from 2000 to the late 2000s before increasing for the remaining duration of the observation period. This same pattern is shown for the total biomass (Fig. S3). The pattern of increase and decrease in the one-year look-ahead elasticities for fishing (Fig. 4) was similar. Figure 4 shows annual fishing elasticities were increasingly negative for all functional groups into the late 2000s, followed by an easing toward zero, except for planktivores, where the estimated fishing elasticities became increasingly negative again after 2011. Unlike the other functional groups, planktivores appeared to plateau and perhaps slightly decrease toward the end of the the observation period (Fig. 3).

The one-year look-ahead elasticities for fishing were closely related to the amount of biomass in the Bay of Biscay ecosystem. The median elasticities of the fishing mortality rates were of the greatest magnitude for all functional groups during the years 2006 or 2007 (Fig. 4). These years corresponded to a period of low biomass (Fig. 3) and high fished fractions (Fig. S4). For planktivores, the magnitude of the fishing mortality elasticity was also high in 2015. This year corresponded to a slight decrease in the estimated biomass of planktivores and an elevated fished fraction relative to the preceding years. In contrast, the fishing mortality elasticity was at a low magnitude for pelagic piscivores in 2015. Pelagic piscivores had at this time an estimated upward trend in biomass and a relatively low fished fraction compared to preceding years.

The density dependence parameters are directly interpretable as elasticities of the one-year look-ahead predictions with respect to functional group biomasses. The posterior mean elasticities for the density dependence relationships are shown in Fig. 5. On average, intra-functional group density dependent elasticities were greater for the pelagic functional groups than for the demersal functional groups. Uncertainty in the estimated intra-functional group density dependent elasticities was however high for both pelagic functional groups (Fig. S5), and comparatively low for demersal functional groups. For pelagic functional groups, it is possible that migratory behaviour increased the variability of intra-functional group density dependence.

Turning to inter-functional group elasticities, Fig. 5 shows that the density-dependent elasticities were higher between pelagic piscivores and benthivores than between pelagic piscivores and planktivores. The positive effect of planktivores on the consumer functional groups on average outweighed the negative predation effect of consumers on planktivores (Fig. 5). The mean estimates shown in Fig. 5 were however accompanied by uncertainty (Fig. S5) with higher magnitudes of inter-functional group density dependence apparently more likely for bottom-up rather than top-down effects. The bottom-up versus top-down impact was therefore compared for each density relationship between the consumer and resource functional groups. The density dependent influence of planktivores on its consumers was somewhat more likely to be stronger than top-down predation pressure for both demersal (P(a_DK > a_KD|C, E, h) = 0.73) and pelagic piscivore (P(a_GK > a_KG|C, E, h) = 0.73) trophic pathways. However, the comparison between top-down and bottom-up magnitudes was more equivocal for benthivores. The posterior probability that the density dependent influence of benthivores on demersal piscivores was greater than the predation pressure was only P(a_DB > a_BD|C, E, h) = 0.58. The posterior probability of a stronger bottom-up benthivore effect on pelagic piscivores was equally likely to a stronger top-down effect (P(a_GB > a_BG|C, E, h) = 0.49).

A further analysis examined the probability that the elasticity magnitudes for inter-functional group density dependence exceeded 0.2 (Table 2). A density dependent elasticity with a magnitude of a_ij = 0.2 predicts that increasing functional group j by a proportional amount would change the biomass of functional group i in the next year by 20% of the same proportional amount. For each pairwise trophic relationship, the bottom-up effect of a resource (prey) functional group on its predator had a greater chance of exceeding 0.2 than the magnitude of the top-down predation effect. This result shows that high magnitudes of density dependence were more likely for positive bottom up effects of prey functional groups (K and B) on their predators (D and G) than for negative predation effects of the consumer functional groups on the prey functional groups. The greater magnitude of bottom-up versus top-down density dependence was however more pronounced for planktivores than for benthivores (Table 2). Whereas the probabilities that the bottom-up elasticities exceeded the 0.2 threshold were an order of magnitude greater than the top-down elasticities for planktivores, for benthivores the exceedance probabilities for bottom-up elasticities were only about double that of the top-down elasticities.

Density independent elasticities were lowest for the pelagic piscivores compared to the other functional groups (Fig. S6). This result may reflect the greater magnitude of process error standard deviation estimated for pelagic piscivores compared to the other functional groups (Fig. S7). High process uncertainty may lead to dynamics driven relatively more by stochasticity rather than the deterministic dynamics of Eq. (5). As noted above, uncertainty was also high for estimates of intra-functional group density dependence of pelagic piscivores, which together with high process uncertainty may reflect an increased importance of migration for this functional group.

Cumulative impacts of fishing pressure

The density dependence captured by the matrix A determines whether or not the ecosystem model is stationary. For stationarity, the maximum of the moduli of the eigenvalues of the transition matrix I + A must be less than one. The Bay of Biscay functional group ecosystem model was stationary with posterior probability of 1. Relative to the unfished long-term distribution of biomass (scenario f⁰ in Fig. 6), the scenario of average fishing mortality f^Avg reduced the stationary distributions of biomass for all functional groups (Fig. 6).

The estimated posterior median elasticities of the long-term averages, medians or other arbitrary quantiles of the stationary distribution (denoted by $\tilde{x}\left(\theta ,y\right)$, see Eq. (12)) with respect to the unfished fractions U were

Table S3 provides the posterior interquartile distances as a measure of the uncertainty associated with the above estimates of the long-term elasticities. The scenario f^G that removes all fishing of pelagic piscivores corresponds to increasing the unfished fraction of pelagic piscivores, and so the last column of the above matrix applies. These elasticities predict that increasing the unfished fraction of pelagic piscivores (equivalently, decreasing the fishing mortality rate f^G) substantially increases pelagic piscivore fishable biomass and also slightly decreases benthivore fishable biomass, in the long-term, relative to scenario f^Avg. Figure 6 shows that the biomass of pelagic piscivores increased relative to the scenario f^Avg of average fishing mortality. Also, the median biomass of benthivores was slightly lowered in scenario f^G because of the greater biomass of predators retained in the ecosystem. The landings of benthivores were also reduced relative to the average scenario f^Avg (Fig. S8). The biomass and landings of planktivores and demersal piscivores were relatively less affected.

From Fig. 1, the mixed trawl fishery targets both benthivores and demersal piscivores. Eliminating effort from the mixed trawl fishery (scenario f^MT) corresponds to a reduction of the fishing mortality rate to 38% of its f^Avg value for benthivores and to 55% of its f^Avg value for demersal piscivores. The unfished fraction is thereby increased for both of these functional groups, and so the cumulative elasticities are obtained by summing the two middle columns of the above matrix −A⁻¹. The resulting predictions are of substantial increases for both benthivores and demersal piscivores, and also a slight increase for pelagic piscivores for the scenario f^MT relative to f^Avg. Figure 6 shows that the biomass of benthivores increased to roughly the same level as in the unfished scenario (f⁰), and the biomass of demersal piscivores also increased. The relative difference between the posterior medians of the long-term biomass of pelagic piscivores also increased by ∼9% when moving from scenario f^Avg to f^MT. The landings of benthivores had the greatest reduction in this last scenario (Fig. S8).

Predictive cross-validation and model adequacy

Predictive cross-validation of the CPUE data was performed across all possible single-year hold-outs from the 2000 to 2015 observation period. Accounting for all targeted gear type and functional group combinations yielded a total of 288 CPUE observations and predictions across the 16 years (Fig. 7). About one-half of the observations were below the corresponding cross-validation predictive medians. About 56% of the observations were within the corresponding 50% central cross-validation predictive intervals, and about 90% of the observations were within the corresponding 90% central cross-validation predictive intervals.

Modelling

The Gompertz model parameters were shown to be directly interpretable in terms of the one-year look-ahead prediction elasticities, which measure the magnitude of the proportional change in next year’s biomass given a proportional change to a parameter or variable in the current year. In particular, these elasticities provided an intuitive interpretation of the time-invariant density dependent parameters within the model, where from Eq. (10) each entry A_ij in the matrix A is the elasticity of functional group i given a small perturbation to the biomass of the jth functional group. The negative of the inverse of the density dependent matrix, −A⁻¹ is on the other hand shown to determine the long-term elasticities of functional group biomasses with respect to their unfished fractions. These long-term elasticities apply to either the long-term means or any choice of marginal quantiles from the stationary distribution.

An advantage of the modelling approach is that it provides coherent estimates of model parameters and functional group biomasses as well as ecosystem process and observation errors. The analysis allowed joint estimation of unknown fishing, process and observation model parameters in addition to the unobserved latent states of functional group biomasses. The joint estimation of these unknown factors is in contrast to most food web models currently used. The price to pay for this coherence is that model complexity needs to be reasonable. The Gompertz model used in this paper, which is a linear process model with multivariate normal process and observation error on the logarithmic scale, is a common choice for the analysis of multivariate ecological time series (e.g., Ives et al., 2003; Spencer & Ianelli, 2005; Lindegren et al., 2009; Bell, Fogarty & Collie, 2014; Torres et al., 2017). Depending on how fishing pressure is modelled, the log-linear normal model may allow for application of analytical updating schemes for the latent states (e.g., Thompson, 1996). For Bayesian state space ecological models with unknown magnitudes of observation error and process noise, alternative choices may for example accommodate non-Gaussian observation error (e.g., Hosack, Peters & Hayes, 2012) or non-linear and non-Gaussian process models (e.g., Hosack, Peters & Ludsin, 2014). On the other hand, the Gompertz model, as a log-linear normal model for the latent states, is a reasonable assumption for many ecological populations (Dennis et al., 2006). The stationary distribution of the Gompertz model is also easily accessible, which accommodates long-term predictions of ecosystem response to management scenarios while accounting for uncertainty.

As an assessment of model adequacy, the Gompertz model was shown to have predictive capability by comparing the cross-validation predictive densities for hold-out years against observed CPUE (Fig. 7). The proportions of observations below the cross-validation predictive medians of the hold-out years was about 50%, as expected. Similarly, the proportions of observations within the 50% and 90% central cross-validation predictive intervals were close to 50% and 90% respectively, even though sometimes conflicting information was contributed by the various commercial and scientific survey CPUE time series. For example, the purse seine CPUE of pelagic piscivores (G) decreased moving into the last year of the time series (2015), whereas the observed CPUE increased for all other targeted gear types of pelagic piscivores in that year. Nevertheless, the cross-validation predictive densities for the hold-out years achieved by the model appeared reasonable (Fig. 7). The demonstrated predictive capability of the multivariate Gompertz model increases confidence in the usefulness of the interpretable one-year look-ahead elasticities and the importance of jointly estimating unknown process and observation uncertainty in the state space model framework.

It should be noted that alternative discretisation schemes are commonly applied to continuous-time theoretical models originally formulated as an ordinary differential equation (ODE) when deriving discrete-time models (e.g., Turchin, 2003), such as that used in this paper. In fisheries, a common form of discretisation has the structure ${\displaystyle x_i\left(t+1\right)=x_i\left(t\right)+G\left(x_i\left(t\right)\right)-C\left(t\right),}$

where C(t) is the harvest or catch in year t. The choice of the function G(⋅) captures different forms of density dependent growth (e.g., Hilborn & Walters, 1992; Haddon, 2011). The function G(⋅) is usually inspired by the form of density dependence as represented by continuous-time models such as the Schaefer model (Schaefer, 1954), Pella-Tomlinson model (Pella & Tomlinson, 1969) or Gompertz-Fox model (Fox, 1970). For the form of discretisation above, high levels of harvest in year t can lead to the prediction that next year’s stock size is negative, which is a biological impossibility. The discretised form of the Gompertz model used in the paper avoids this problem by instead choosing the exponential per capita growth function of Eq. (8), which precludes biologically impossible predictions of negative stock size. The resulting equilibrium of the deterministic version of the discrete multivariate Gompertz model used here is shown to be identical to the multivariate generalisation of the original deterministic continuous-time Gompertz-Fox ODE model introduced by Fox (1970).

Ecosystem dynamics

After a decrease in the early 2000s, model results estimate that total (fishable) biomass steadily increased in the Bay of Biscay ecosytem from around 2008. This pattern was supported by increases in all functional groups, except planktivores, which remained stable during the latter portion of the observation period. The model results reflect the decreasing trends that are evident for several of the CPUE time series early in the study period (Fig. S2). In the Bay of Biscay during the early 1990s, ten out of 20 assessed fish stocks were considered overexploited and 28 fish species were classed as subject to concern by the IUCN (Lorance et al., 2009). Further, empirical analysis indicates that total demersal biomass had remained stable between the late 1980s and the beginning of the study period (Rochet et al., 2005). Thus, the Bay of Biscay ecosystem was clearly impacted by overfishing prior to the study period and seemed to have remained so during the first part of the study period, despite the fact that fishing power has decreased since the early 1990s (Mesnil, 2008) with an acceleration in the early 2000s (Rochet, Daurès & Trenkel, 2012). However, time delays and functional groups responding in different time frames to a reduction in fishing pressure is expected due to differences in life history traits (Collie, Rochet & Bell, 2013). Planktivores are expected to rebound faster, which could explain the stability of biomass found here for this group.

A previous study from the neighbouring Celtic Sea found similar patterns of biomass fluctuations for a short period of overlapping years in the early 2000s (Bell, Fogarty & Collie, 2014). That study analysed the same functional groups in the years 1987–2004 while also using a Gompertz model, although the analysis differed in other respects such as the way in which fishing pressure was included and the inferential framework. Nevertheless, Bell, Fogarty & Collie (2014) estimated decreasing total biomass in the Celtic Sea during the early 2000s. All functional groups were estimated to be decreasing in the Bay of Biscay during this period (Fig. 3 and Fig. S3). The Celtic Sea study did suggest an increase in planktivores during the early 2000s. However, only trawl survey data and no commerical CPUE was included (Bell, Fogarty & Collie, 2014). The survey CPUE from the Bay of Biscay also suggest a small increase in planktivore biomass in the early 2000s, whereas the CPUE from nets and purse seines suggest a decrease (Fig. S2). The estimated latent biomass of planktivores (Fig. 3) appears to track these latter CPUE indices rather than the scientific survey data during this period in the Bay of Biscay. The scientific surveys, although fishery-independent, also have a limited spatio-temporal footprint compared to the commercial fisheries. The scientific survey may also target different age classes, size classes or species from fisheries that target fishable biomass. The commercial fishery and scientific surveys were therefore a priori given equal weightings by placing independent and identical priors on the observation error standard deviations associated with each gear type. The relative weighting of the scientific survey, or any of the commercial fisheries, could however be altered by changing these priors to capture additional information contributed by expert opinion or empirical data.

Significant, though variable, intra-functional group and inter-functional group density dependence was identified in the Bay of Biscay. The estimated intra-functional group density dependence was more uncertain for the pelagic functional groups than the demersal functional groups. This difference between demersal and pelagic groups could be a result of density mediated food and habitat limitations being more important for the demersal groups, whereas environmental conditions are known to be important drivers of population dynamics for pelagic species, in particular small and medium sized pelagics (Trenkel et al., 2014). The predator functional group of pelagic piscivores also exhibited high uncertainty in process model uncertainty (e.g., the process noise standard deviation). For pelagic piscivores, it is possible that migratory behaviour increases the uncertainty of intra-functional group density dependence, where year to year changes in biomass may be more sensitive to allochthonous import and export rather than endogenous production.

Top-down control by predators is expected to dampen ecosystem level fluctuations in marine ecosystems (Cury, Shannon & Shin, 2003). Although some food-web models are pre-determined as either bottom-up or top-down by construction e.g., Steele (2009), the ability to estimate the strength of top-down and bottom-up relationships for each predator–prey pair is an advantage of the state-space modelling approach. There was evidence from the one-year look-ahead elasticities for top-down control being more important for benthivores in the Bay of Biscay than top-down control of planktivores (Fig. 5 and Fig. S5). Overall, however, high magnitudes of elasticities were more probable for bottom-up effects of producers on consumers than top-down effects (Table 2). Using a detailed simulation model, Oken & Essington (2015) found that a top-down predator effect can be difficult to detect in a biomass production model, in particular if recruitment variation and observation error are high, or a second predator also targets the species. All of these factors occur in the Bay of Biscay ecosystem, which may have made detection of top-down control difficult.

Fishing

Seven fishing fleets were defined by broad gear categories that exploited from one functional group (potters) to three functional groups (pelagic trawlers and netters). To make catchability estimates as much as possible comparable among functional groups within fleets, but also among fleets, the fishing effort was defined by the common currency of days-fishing. The fleet with the highest catchability was purse seine fishing on planktivores (posterior median catchability equal to 0.18, Table S2). This result is not surprising given the type of fishing method and the relatively small size range of the species making up the group. The lowest catchability was estimated for netters fishing pelagic piscivores (median catchability equal to 0.001). Interestingly, the catchabilities of demersal piscivores was comparable for mixed trawlers (median catchability 0.006) and the bottom trawl used in the fishery independent survey (median catchability 0.005). A higher catchability might have been expected for mixed trawlers given that this fishery targets large individuals above the minimum landing size. However, the observation error standard deviation was lowest for mixed trawls compared to all other gear types including the scientific survey.

The model results presented are of fishable biomass exploitable by the commercial fisheries that operate in the Bay of Biscay. Estimates of total biomass would require independent estimates of functional group catchabilities for the commercial fisheries. To the best of our knowledge, this is the first time catchability coefficients have been estimated on the functional group level within a state space model. Typically, catchability is estimated at the species level either experimentally (e.g., Myers & Hoenig, 1997; Millar & Fryer, 1999) or at the population level in stock assessment models. Time-varying catchability may also be estimated within a state-space model framework (Wilberg et al., 2009). For example, Hosack, Peters & Ludsin (2014) found evidence that catchability of commercial fisheries may vary with a proxy for habitat quality in a joint two-species state space model with unknown observation error and process uncertainty. However, the difficulties of parameterising complex models of catchability are heightened for ecosystem models that take into account many species with alternative life history strategies and behaviour. For the scientific survey used in this study, catchability has been estimated experimentally revealing large variability among species (Trenkel & Skaug, 2005). Such species level differences are averaged out at the functional group level estimates obtained here. Indeed, functional group and community level survey estimates have been found more robust to variability (season or gear type) compared to the species level (Trenkel et al., 2004).

The time-varying elasticities for fishing mortality showed how small perturbations of fishing mortality rates have temporally varying strong or weak effects on the one-year look-ahead predictions (Fig. 4). When biomass was low and the fished fraction of a functional group was high, then the changes in the fishing mortality rate had stronger effects on the one-year look-ahead predictions compared to periods when biomass was high and the fished fraction low. At first it may seem surprising that pelagic piscivores supported comparably high rates of fishing mortality as other functional groups (Fig. 4) despite having the lowest rates of density independent growth (r_G in Fig. S6). However, pelagic piscivores received positive biomass inputs from the density dependent trophic relationships with planktivores and benthivores (Fig. 5 and Fig. S5) that counteracted the low rates of density independent growth. The implication is that the level of sustainable fishing pressure for a given functional group may depend on other functional groups in the ecosystem.

Assessing the cumulative effects of multiple fisheries on marine ecosystems continues to be a concern for fisheries management (Hobday et al., 2011; Link & Marshak, 2019; Zhou et al., in press). However, the complex and nonlinear dynamics of marine ecosystems limit their predictability (Glaser et al., 2014; Planque, 2016). Therefore the long-term predictions of only four rather extreme fishing scenarios were studied. The results nevertheless provide insight into the potential effects of interactions between fisheries management and food web dynamics caused by direct and indirect effects. The posterior stationary distribution under each scenario conformed with predictions derived from the long-term elasticities of the biomasses with respect to changes in fishing pressure (see Eq. (12)). The scenario with no fishing on pelagic piscivores (f^G), which might be expected to increase top-down control on planktivores, led to a similar stationary biomass distribution of planktivores as status quo fishing (f^Avg) or excluding mixed trawlers (f^MT). The scenario f^G instead decreased the stationary biomass of benthivores: the inter-functional group density dependence relationship with the greatest (median or mean) magnitude was the negative density dependent effect of pelagic piscivores on benthivores (Fig. 5 and Fig. S5). Both of these results agreed with the long-term elasticities of the biomasses with respect to changes in fishing pressure. Not surprisingly, the scenario with no mixed trawlers (f^MT) mainly benefited benthivores and demersal piscivores relative to the status quo fishing scenario (f^Avg) due to the release of fishing pressure on these two functional groups, although the biomass of pelagic piscivores also increased slightly. Again these results agreed with predictions from the long-term elasticities. From a fisheries management point of view the scenario results indicate that indirect effects might be expected in the Bay of Biscay when taking strong management actions for certain fleets or functional groups.