Addressing the dichotomy of fishing and climate in fishery management with the FishClim model

Data

Sea Surface temperature (1850–2019)

Sea Surface Temperature (SST, °C) from 1850 to 2019 originated from the COBE SST2 1° × 1° gridded dataset⁷⁴, https://psl.noaa.gov/data/gridded/data.cobe2.html. SST data were interpolated on a 0.25° latitude × 0.25° longitude grid on a monthly scale from 1850 to 2019.

Bathymetry

Bathymetry (m) came from GEBCO Bathymetric Compilation Group 2019 (The GEBCO_2019 Grid—a continuous terrain model of the global oceans and land). Data are provided by the British Oceanographic Data Centre, National Oceanography Centre, NERC, UK. doi:10/c33m. (https://www.bodc.ac.uk/data/published_data_library/catalogue/10.5285/836f016a-33be-6ddc-e053-6c86abc0788e/). These data were interpolated on a 0.25° latitude × 0.25° longitude grid.

Biological data

Daily mass concentration of chlorophyll-a in seawater (mg/m³) originated from the Glob Colour project (http://www.globcolour.info/). The product merges together all the daily data from satellites (MODIS, SeaWIFS, VIIRS) available from September 1997 to December 2019, on a 4 km resolution spatial grid. These data were interpolated on a daily scale on a 0.25° latitude × 0.25° longitude grid. These data were only used to map the average maximum standardised SSB (mdSSB) around the North Sea (Fig. 1a). When long-term changes in mdSSB were examined, we used modelled chlorophyll data (see section “Climate projections” below).

Cod recrutment at age 1, Spawning Stock Biomass (SSB) and fishing effort F for 1963–2019 originated from ICES³⁵.

We used a plankton index of larval cod survival, which was an update of the index proposed by Beaugrand and colleagues³³. Based on data from the Continuous Plankton Recorder (CPR)⁷⁵, the index is based on the simultaneous consideration of six key biological parameters important for the diet and growth of cod larvae and juveniles in the North Sea:^76,77 (i) Total calanoid copepod biomass as a quantitative indicator of food for larval cod, (ii) mean size of calanoid copepods as a qualitative indicator of food, (iii-iv) the abundance of the two dominant congeneric species Calanus finmarchicus and C. helgolandicus, (v) the genus Pseudocalanus and (vi) the taxonomic group euphausiids. A standardised Principal Component Analysis (PCA) is performed on the six plankton indicators for each month from March to September for the period 1958–2017 (table 60 years × 7 months-6 indicators). The plankton index is simply the first principal component of the PCA³³.

Climate projections

Climate projections for SST and mass concentration of chlorophyll in seawater (kg m⁻³) originated from the Coupled Model Intercomparison Project Phase 6 (CMIP6)⁵ and were provided by the Earth System Grid Federation (ESGF). We used the projections known as Shared Socioeconomic Pathways (SSP) 245 and 585 corresponding respectively to a medium and a high radiative forcing by 2100 (2.5 W m⁻² and 8.5 W m⁻²)⁷⁸. The daily simulations of four different models (i.e. CNRM-ESM2-1, GFDL-ESM4, IPSL-CM6A-LR, and UKESM1-0-LL) covering the time period 1850–2014 (historical simulation) and 2015–2100 (future projections for the two SSPs scenarios) were used. All the data were interpolated on a 0.25° by 0.25° regular grid. Key references (i.e. DOI and dataset version) are provided in Supplementary Text 1. Long-term changes in modelled SSB were based on these data (including modelled daily chlorophyll data).

The FishClim model

Let K_t be the maximum standardised Spawning Stock Biomass (mdSSB hereafter) that can be reached by a fish stock at time t for a given environmental regime φ_t. X_t+1, standardised SSB (dSSB hereafter) at time t+1 was calculated from dSSB at time t as follows:

α is the fishing intensity that varies between 0 (i.e. no fishing) and 1 (i.e. 100% of SSB fished in a year). It is important to note that α (see Eq. (10)) should not be mistaken with ICES fishing effort F⁷⁹ (calculated from SSB). The second term of Eq. (1) is the intrinsic growth rate of the fish stock that is a function of both K_t and the population growth rate r (r was fixed to 0.5 in most analyses, but see Fig. 3d however where r varied from 0.25 to 0.75). The population growth rate r is highly influenced by the life history traits of a species⁸⁰ but also by environmental variability^54,55,81. Here, the population growth rate was assumed to be constant in space and time and the influence of environmental variability occurred exclusively through its effects on K_t. We made this choice to not multiply the sources of complexity and errors (i.e. population growth rate is very difficult to assess and varies with age⁸⁰). The third term reflects the part of dSSB that is lost by fishing. Note that natural mortality is not explicitly integrated in Eq. (1) because this process is difficult to assess with confidence at the scale of our study. Here, we assumed that the second term of Eq. (1) implicitly considered this process; when K increases, it is likely that natural mortality diminishes, especially at age 1³⁴. We tested this assumption below. Most of the time when fishing occurs, X_t<K. But in case of a strong negative environmental forcing at a time of small fishing intensity, X_t can be transitory above K.

Maximum dSSB K_t at time t was assessed using a niche model based on the MacroEcological Theory on the Arrangement of Life (METAL)⁸² using SST, an index of food availability based on daily mass concentration of chlorophyll in seawater and bathymetry. The model was therefore based on a three-dimensional niche: thermal, bathymetric and trophic niches.

The thermal niche was asymmetrical. Asymmetric niches can be modelled by using a Gaussian function⁸³ with the same ecological optimum y_opt but two different standard deviations t₁ and t₂, i.e. two different ecological amplitudes:

Here y_opt= 5.4 °C and t₁ and t₂ were fixed to 5.7 °C and 4 °C, respectively, so that the thermal niche was close to that assessed by Beaugrand and colleagues³¹ (Supplementary Fig. 2). This Supplementary Figure compares the thermal response curve we chose in the present study with the data analysed in Beaugrand and colleagues³¹. The figure shows that the response curve (red curve) is close to the histogram showing the number of geographical cells with a cod occurrence as a function of temperature varying between −2 °C (frozen seawater) and 20 °C.

Because t₁ > t₂, the niche was slightly negative asymmetrical (Supplementary Fig. 1). U₁(y) was the first component of mdSSB along the thermal gradient y. c was the maximum value of mdSSB; c was fixed to 1 so that mdSSB varied between 0 and 1^84,85. y was the value of SST. Slight variations in the different parameters of the niche did not alter either the spatial patterns in the distribution of mdSSB nor the correlations with recruitment.

To model the bathymetric niche of cod, we used a trapezoidal function. Changes in mdSSB, U₂, along bathymetry, were assessed using four points (θ₁, θ₂, θ_3, θ₄):

With θ₂ ≥ θ₁, θ₃ ≥ θ₂ and θ₄≥ θ₃ and y the bathymetry; θ₁ = 0, θ₂ = 10⁻⁴, θ₃ = 200 and θ₄ = 600 m (Supplementary Fig. 1). These parameters were retrieved from the litterature^86,87. Here also c, the maximum abundance reached by the target species was fixed to 1 and U₂ varied between 0 and 1. Trapezoidal niches have been used frequently to model the spatial distribution of fish and marine mammals^88,89.

The trophic niche was modelled by a rectangular function on a daily basis. To the best of our knowledge, no information on the trophic niche is available. We modelled the trophic niche by fixing U₃ to 1 when chlorophyll-a concentration was higher than 0.05 mg m⁻³ during a minimum period of 15 days and 0 otherwise (Supplementary Fig. 1). This minimum of chlorophyll was implemented as a proxy for suitable food, which has been shown to be important in the North Atlantic for cod recruitment and distribution^6,33.

There exists two ways to combine the different ecological dimensions of a niche: (i) use an additive or (ii) a multiplicative model^82,90. We used a multiplicative model because when one dimension is associated to a nil abundance, the resulting abundance combining all dimensions is also nil in contrast to an additive model; therefore only one unsuitable environmental value may explain a nil abundance. All dimensions were associated to abundance values that varied between 0 and 1⁹⁰.

Therefore, maximum dSSB, K, for a given environmental regime E was given by multiplying the three niches (thermal, bathymetric and trophic):

where p = 3, the three dimensions of the niche.

Analyses

Mapping of maximum standardised SSB

mdSSB is close to the “dynamic B₀” approach; B₀ is the SSB in the absence of fishing (generally expressed in tonnes)⁵¹ whereas mdSSB is the SSB in the absence of fishing standardised between 0 and 1 and assessed from the knowledge of the niche of the species. We first assessed mdSSB in the North-east Atlantic (around UK) at a spatial resolution of 0.25° latitude × 0.25° longitude on a daily basis from 1850 to 2019. For this analysis, FishClim was run on monthly COBE SST (1850–2019), mean bathymetry and a climatology of daily mass concentration of chlorophyll-a in seawater from the Glob Colour project (see Data section). We then calculated an annual average based on the main seasonal productive period around UK, i.e. from March to October⁹⁰. Finally, we averaged all years to examine spatial patterns in mean mdSSB (Fig. 1a).

Temporal changes in maximum standardised SSB

We assessed average long-term changes in mdSSB in the North Sea (51°N–62°N and 3°W–9.5°E); the annual average was calculated from March to October because this is a period of high production⁹⁰ . We compared long-term changes in mdSSB with cod recruitment at age 1, a plankton index of larval cod survival based on the period March to October³³, and ICES-based SSB³⁵ for 1963-2019 (Fig. 1b–d).

Correlation analyses with modelled maximum standardised SSB

Pearson correlations between long-term changes in mdSSB (average for the North Sea, 51°N–62°N and 3°W–9.5°E) and cod recruitment at age 1 in decimal logarithm³⁵, a plankton index of larval cod survival in the North Sea³³, and observed ICES SSB in decimal logarithm³⁵ for the period 1963–2019 were calculated (Fig. 1b–d). The same analysis was performed between assessed fishing intensity α from our FishClim model and fishing effort F³⁵ in the North Sea (Fig. 1e). The probability of significance of the coefficients of correlation was adjusted to correct for temporal autocorrelation⁹¹.

Assessment of fishing intensity from ICES spawning stock biomass

Using North Sea ICES SSB, we applied Eq. (1) to assess fishing intensity α:

With X_t+1 and X_t the ICES dSSB (in decimal logarithm). Standardisation of ICES SSB, necessary for this analysis, was complicated because many different kinds of standardisation were achievable so long as X remained strictly above 0 (i.e. full cod extirpation, not observed so far³⁵) and strictly below min(K) (i.e. all black curves always below all points of the blue curve were possible, Supplementary Fig. 3). Indeed, ICES SSB includes exploitation and environmental fluctuations whereas K (i.e. mdSSB) integrates only environmental forcing; the difference is mainly caused by the negative influence of fishing. We chose the black curve (ICES SSB) that maximised the correlation between α (fishing intensity in the FishClim model) and F (ICES fishing effort)³⁵.

Reconstruction of long-term changes in ICES spawning stock biomass

The estimation of α allowed us to reconstruct long-term changes in cod (ICES) dSSB and to examine the respective influence of fishing and CIEC by means of Eq. (1) (“Methods”) using four hypothetical scenarios (Fig. 1f). First, we fixed fishing intensity and considered exclusively environmental variations through its influence on dSSB. (i–ii) We assessed long-term changes in dSSB from long-term variation in observed mdSSB (called K_t in Eq. (1)) with a constant level of exploitation fixed to (i) minimum (upper blue curve, i.e. the lowest fishing intensity observed in 1963–2019) or (ii) maximum (lower blue curve, i.e. the highest fishing intensity observed in 1963–2019).

Second, we fixed the environmental influence on dSSB and considered variations in fishing intensity. We estimated long-term changes in dSSB from long-term variation in estimated α with a constant mdSSB fixed to (iii) minimum (lower red curve, i.e. the lowest mdSSB observed in 1963–2019) or (iv) maximum (upper red curve, i.e. the highest mdSSB observed in 1963–2019). It was possible to compare long-term changes in reconstructed (ICES) dSSB (thick black curve in Fig. 1f) with these four hypothetical scenarios (Fig. 1f); note that these comparisons were not affected by the choice we made earlier on the standardisation of (ICES) SSB.

Quantification of the respective influence of fishing and climate/environment on spawning stock biomass

Using the previous curves, we examined the respective influence of fishing and CIEC on reconstructed (ICES) dSSB (Fig. 2). First, the influence of fishing was investigated by estimating the residuals between reconstructed (ICES) dSSB based on long-term changes in mdSSB (i.e. K_t in Eq. (1)) and α (thick black curves in Fig. 1f) and modelled dSSB based on fluctuating fishing intensity α and invariant K (average of the two red curves in Fig. 1f). This calculation led to the red curve in Fig. 2b. Next, we performed the opposite procedure to examine the influence of CIEC on dSSB (i.e. invariant fishing intensity α based on the two blue curves in Fig. 1f). This calculation led to the blue curve in Fig. 2b.

A cluster analysis, based on a matrix years × three time series with (i) long-term changes in reconstructed standardised (ICES) SSBs, (ii) fishing and (iii) CIEC, was performed to identify key periods (vertical dashed lines in Fig. 2). We standardised each variable between 0 and 1 and used an Euclidean distance to assess the year (1963–2019) × year (1963–2019) square matrix so that each variable contributed equally to each association coefficient. We used an agglomerative hierarchical clustering technique using average linkage, which was a good compromise between the two extreme single and complete clustering techniques⁹². In this paper, we were only interested in the timing between the different time periods (i.e. the groups of years) revealed by the cluster analysis (Fig. 2).

We also calculated an index of fishing influence (ε, expressed in percentage) by means of two indicators γ and δ, which were slightly different to the ones we used above. The first one, γ, was modelled dSSB with fluctuating fishing intensity and a constant mdSSB based on the best suitable environment observed during 1963–2019 (only the upper red curve in Fig. 1f; fishing influence). The second one, δ, was modelled dSSB based on fluctuating environment and fishing intensity (black curve in Fig. 1f) on modelled dSSB based on a fluctuating environment but a constant fishing intensity fixed to the lowest value of the time series (only the upper blue curve in Fig. 1f; environmental influence). The index of fishing influence (ε, expressed in percentage) was calculated as follows:

For each period of 1963–2019 identified by the cluster analysis, we quantified the influence of fishing (and therefore the environment) using a Jackknife procedure^93,94. The resampling procedure recalculated ε by removing each time 1 year of the time period, which allowed us to provide a range of values (i.e. minimum and maximum) in addition to the average value (bar{varepsilon }) calculated for each interval, including the whole period (Fig. 2c).

Long-term changes in modelled spawning stock biomass (1850–2019, 2020–2100 and 2020-2300)

We modelled mdSSB (K_t in Eq. (1)) using outputs from four Earth System models (ESMs) based on two scenarios of SST/Chlorophyll changes (i.e. SSP245 and SSP585) for the period 1850–2100 (and for one scenario and one ESM until 2300; Fig. 3).

For the period 1850–2019, we used daily SST/Chlorophyll changes from the four ESMs to estimate potential changes in mdSSB (thin dashed black curves in Fig. 3a). An average of mdSSB was also calculated (thick green curve in Fig. 3a).

For the period 2020–2100, we showed all potential changes in mdSSB based on the four ESMs and both scenarios SSP245 (thin dashed blue curves in Fig. 3a) and SSP585 (thin dashed red curves). An average of mdSSB was also calculated for scenarios SSP245 (thick continuous blue curve) and SSP585 (thick continuous red curve). In addition, we assessed dSSB based on a constant standardised catch fixed to the average of 2008–2019, the last period identified by the cluster analysis (G5, i.e. (alpha X) = 0.03 in Eq. (1)), and the average values of all ESMs for SSP245 (thick dashed blue curve in Fig. 3a) and SSP585 (thick dashed red curve). This analysis was performed to show how a constant catch might alter long-term changes in mdSSB. When X_t (Eq. (1)) reached 0.1, the stock was considered as fully extirpated.

Understanding how fishing and climate/environment interact now and in the future

We modelled dSSB as a function of fishing intensity α and CIEC to show how fishing and the environment interact (Fig. 3b, c). We calculated dSSB for fishing intensity between α = 0 and α = 0.5 every step Ɵ = 0.001 and for mdSSB between K = 0 and K = 1 every step Ɵ = 0.001 to represent values of dSSB as a function of fishing and CIEC. We then superimposed reconstructed ICES dSSB (1963–2019) on the diagram for three periods: 1963–1985 (high SSB), 1986–1999 (pronounced reduction in SSB), and 2000–2019 (low SSB). Maximum standardised SSB for 2020–2100 (or 2300 exclusively for Scenario SSP 585 of IPSL ESM) assessed from four ESMs and scenarios SSP245 and SSP585 were also superimposed. Fishing intensity is unpredictable for 2020–2100 and so we arbitrarily fixed it constant between 0.08 and 0.17 in increments of 0.1 for display purposes, starting by ESMs based on scenario SSP 245 followed by scenario SSP 585 (Fig. 3b). When X_t (Eq. (1)) reached 0.1, the stock was considered as fully extirpated.

We calculated an index of sensitivity of dSSB as a function of fishing intensity and CIEC. To do so, we first calculated sensitivity of dSSB to fishing intensity α. Index ζ_i was calculated at point i from dSSB X and fishing intensity α at i−1 and i+1 (see also Eq. (1)):

With min(α) = 0, max(α) = 0.5 and Ɵ = 0.001.

Similarly, we calculated sensitivity of dSSB to K. Index η_j was calculated at point j from dSSB X and mdSSB K at j−1 and j+1 (see also Eq. (1)):

With min(K) = 0, max(K) = 1 and Ɵ = 0.001.

Then, we summed the two indices to assess the joint sensitivity of dSSB to fishing intensity Z and mdSSB H:

Matrix I was subsequently standardised between 0 and 1:

With I* the matrix of sensitivity of dSSB to fishing intensity and mdSSB standardised between 0 and 1 (Fig. 3c).

Number of years needed for recovery after stock collapse

We investigated how the number of years needed for a stock to recover after stock collapse (i.e. dSSB=0.05 in Eq. (1); i.e. 10% of mdSSB) varied as a function of mdSSB (between 0 and 1 by increment of 0.001); this was only influenced by the environmental regime φ_t and population growth rate r. For this analysis, we fixed a target dSSB of 0.4 (vertical dashed green vertical line in Fig. 3d) and three different values of r: 0.25, 0.5 and 0.75. We simulated a hypothetical moratorium with a fishing intensity α = 0 in Eq. (1).

Here, stock collapse was defined as dSSB ≤ 0.1 × mdSSB, i.e. when the dSSB reached less than 10% of the unfished biomass mdSSB. This threshold corresponds to values usually defined in the literature; e.g. Pinsky and colleagues⁹⁵ defined a collapse when landings are below 10% the average of the five highest landings recorded for more than 2 years, Worm and colleagues⁶⁹ defined stock collapse when the biomass becomes lower than 10% of the unfished biomass, Andersen⁹⁶ when it is lower than 20% and Thorpe and De Oliveira⁶⁷ when it is lower than 10–20%.

Potential consequences of fisheries management and climate-induced environmental changes

We examined how fishing and CIEC may affect cod stocks and their exploitation around UK with a focus in the North Sea (Figs. 4, 5). For these analyses, we averaged long-term changes in modelled dSSB corresponding to each scenario (all thin dashed blue and thin red curves in Fig. 3a for SSP245 and 585, respectively). In these analyses, the stock was considered fully extirpated when X_t (Eq. (1)) reached 0.1.

Year of cod extirpation for 2020–2100

We estimated year of cod extirpation from 2020 to 2100 in each geographical cell based on (i) a constant fishing intensity (α = 0.04) in time and space, and (ii) an adjusted fishing intensity using the concept of Mean Sustainable Yield (MSY). The choice of α = 0.04 did not alter our conclusions; a lower or a higher value delayed or speed cod extirpation in a predictable way, respectively.

In fisheries, MSY is defined as the maximum catch (abundance or biomass) that can be removed from a population over an indefinite period with dX/dt = 0, with X for dSSB and t for time. Despite some criticisms about MSY⁶⁶, the concept remains a key paradigm in fisheries management^35,63. We used this concept to show that controlling fishing intensity delayed cod extirpation. From Eq. (1), we calculated fishing intensity, called α_MSYt, so that X remained above X_MSYt at all time t:

In this analysis, we fixed X_{MSY t} = K_t/2.

We assessed ({alpha }_{{{{{rm{MSY}}}}t}}) from Eq. (16) and then estimated dSSB from ({alpha }_{{{{{rm{MSY}}}}t}}) and K_t (based on averaged SSP245 and SSP585) by means of Eq. (1).

Although results were displayed at the scale of the north-east Atlantic (around UK), we calculated the difference in year of cod extirpation between scenarios of warming (SSP245 and SSP585) and between scenarios of cod management (constant versus adjusted—MSY— fishing intensity). Differences were presented by means of histograms (Fig. 4). From each histogram, we calculated the median of the differences in year of cod extirpation E⁹⁷.

Pooled standardised catch by 2100 (2020–2100)

In term of fishing exploitation, we assessed pooled standardised catch (i.e. pooled dSSB) in 2100 (2020–2100), again for two scenarios of CIEC (SSP245 and 585) and two scenarios of cod management (constant versus adjusted—MSY—fishing intensity; Fig. 5). We then calculated the percentage of reduction in pooled standardised catch caused by fishing or the intensity of warming. Finally, we assessed the median of the percentage of reduction in pooled standardised catch for the North Sea area (51°N–62°N and 3°W–9.5°E). The goal of this analysis was to demonstrate that controlling fishing intensity optimises cod exploitation.

Statistics and reproducibility

All statistical analyses can be reproduced from the equations provided in the text, the cited references or the data available in Supplementary Data.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com