Projecting contributions of marine protected areas to rebuild fish stocks under climate change

Fish and invertebrates stocks in the Northeast Atlantic

We included 739 exploited fish and invertebrates stocks (hereafter referred collectively as fish stocks) in the Northeast Atlantic ocean. Stock is defined here in the fisheries context, which is a fish or invertebrate species that is exploited in a spatial area unit. Following previous studies on assessing the current status of fish stocks⁹ and projecting their future under climate change and fishing scenarios¹⁹, we used marine ecoregions⁵⁷ to delineate the spatial area unit. We focused on species that were exploited in eight marine ecoregions in the Northeast Atlantic (Fig. 5). Specifically, stocks are defined as species occurring in a marine ecoregion with estimated catches. This includes 231 species (202 fishes and 29 invertebrates) with estimated catches between 2000 – 2019 in the Sea Around Us catch database⁵⁸ (www.seaaroundus.org, Table S1). Although uncertainties exist in the spatial allocation and reconstruction of catches in the Sea Around Us data⁵⁹, such uncertainties would only affect the identification of a stock and would not affect the projected relative changes in biomass and catch under climate change.

Red shows protected cells, blue shows surrounding waters. The boundary of marine ecoregions is delineated by white lines.

Climate-fish-fisheries model

We projected changes in biomass and potential catches of exploited species using a linked climate-fish-fisheries model called dynamic bioclimate envelope model (DBEM). The structure of DBEM is described in Cheung et al^19,55. In brief, the model has a horizontal spatial resolution of 0.5^° latitude x 0.5° longitude for the sea surface and bottom, and simulates annual average abundance and catches of each modeled species. DBEM uses spatially explicit outputs from coupled atmospheric-ocean biogeochemical models, including temperature, oxygen level, salinity, surface advection, sea ice extent and net primary production. Sea bottom and surface temperature, oxygen and salinity are used for demersal and pelagic species, respectively. These outputs are then used to calculate an index of habitat suitability for each species in each spatial cell. Other information used to calculate habitat suitability includes bathymetry and specific habitats (coral reef, continental shelf, shelf slope, and seamounts). Changes in carrying capacity in each cell is assumed to be a function of the estimated habitat suitability, and net primary production in each cell.

The model simulates the net changes in abundance in each spatial cell based on logistic population growth, fishing mortality, and movement and dispersal of adults and larvae modeled through advection–diffusion–reaction equations. Specifically, pelagic larval dispersal is dependent on ocean current (simulated by the Earth system models) and pelagic larval duration estimated from an empirical equation. Dispersal of adults is dependent on spatial gradients of density relative to the carrying capacity (dependent on environmental habitat suitability) and species’ mobility (e.g., large pelagic fishes have high movement rate while sessile species have negligible movement rate).

DBEM then calculates biomass from abundance using a characteristic weight representing the average mass of an individual in the cell. The model simulates how changes in temperature and oxygen content would affect the growth of the individual using a submodel derived from a generalized von Bertalanffy growth function. DBEM has a spin-up period of 100 years using the climatological average oceanographic conditions from 1951 to 2000, thereby allowing the species to reach equilibrium before it is perturbed with oceanographic changes. Previous studies have found strong correlation between catch data and DBEM projections from marine regions^55,60. Although DBEM can account for the potential effects of ocean acidification (changes in pH) on growth, reproduction and survivorship^61,62, we did not include ocean acidification in this study because of its large variation of effects on exploited marine species.

In our model, fishing intensity was assumed in the fisheries scenarios, represented by fishing mortality rates (F) relative to the F required to achieve MSY (i.e., F/F_MSY). As the fish model assumes logistic population growth, following the derivation from a simple surplus production model, FMSY is approximately equal to half of the intrinsic growth rate of each species⁶³.

Climate and fishing scenarios and analysis

DBEM was forced with projections from the Geophysical Fluid Dynamics Laboratory (GFDL)-ESM4⁶⁴. The variables that we extracted from the Earth system model simulation include global mean surface atmospheric temperature, sea surface and bottom temperature, dissolved oxygen concentration and salinity, vertically integrated total net primary production, sea ice extent, and surface advection.

Projections followed two contrasting scenarios—shared socio-economic pathway (SSP) 1—representative concentration pathway (RCP) 2.6 (SSP1-2.6) and SSP5-8.5^65,66. The SSP1-2.6 and SSP5-8.5 represent a “strong mitigation” low-emissions pathway and a “no mitigation” high-emissions pathway, respectively. The simulation time frame is from 1950 to 2100. The Earth system model projections consider 1950 to 2010 as a historical period that diverges into the two climate change scenarios from 2011 to 2100.

We re-expressed the simulated changes in ocean conditions and fish stocks under the SSPs according to the respectively projected global atmospheric warming levels. Different Earth system models and their versions vary in their projected intensity of climate change⁶⁷. One commonly-used option to account for such variabilities in projections, especially for impacts and risk assessments, is to express response variables relative to global mean atmospheric warming level instead of time frame⁶⁷. In this study, the simulated annual response variables from the climate-fish-fisheries model were related to the global mean atmospheric warming levels at the respective year and SSP. Such approaches to analyze projected climate impacts on fish stocks and fisheries have been used in previous studies^19,30,68. We included four fishing scenarios and three no-take marine protected area scenarios. The F/F_MSY scenarios included 0.5, 0.75, 1 and 1.5, with 1 being at MSY level, and 1.5 at over-exploited level (Table 1). We implemented the fishing scenarios across the simulation time frame (1950 – 2100). Thus, these are idealized fishing scenarios intended for theoretical explorations of the effects of no-MPAs under climate change.

No-take marine protected areas scenarios and analysis

No-take marine protected area scenarios were expressed as the area of the Northeast Atlantic Ocean with no fishing (5%, 15% and 30% of the modeled ocean area). The analysis was limited to the Food and Agricultural Organization of the United Nations (FAO) major fishing zone 27-Atlantic Northeast (https://www.fao.org/fishery/en/area/27). The region was gridded into 0.5° latitude x 0.5° longitude grid cells. We identified the specific grid cells to be protected from fishing using the statistical software R. The 5% protection scenario includes the proportion of cells within FAO area 27 that are currently occupied by MPAs (by calculating the spatial coverage of protected areas relative to the total area of FAO zone 27 following the world database of protected areas⁶⁹). The computed area (5.5% of the total area of the FAO area) is different from the actual designated area of the MPAs because of the coarser grid resolution of our model relative to the size of the MPA. The larger no-take MPAs coverage scenarios were then built on the smaller ones by progressively adding protected areas through randomly designating locations to be no-take MPAs from the 5, to 15 and 30% protection scenarios. Thus, the 5 and 15% scenarios are subsets of the 30% scenario. The average sizes of no-take MPAs patches (i.e., group of neighboring pixels that were assigned as MPAs) were 1246 ± 1056 km² (standard deviation, n = 614 patches), 1555 ± 1364 km² (n = 1379 patches) and 1646 ± 1510 km² (n = 2571 patches) for the 5, 15 and 30% scenarios, respectively.

Under each MPAs scenario, we reallocated fishing mortality from the protected area grid cell to the immediate surrounding grid cells, simulating the redistribution of fishing effort that often occurs when an area is closed to fishing (Fig. 5). To do this, we assigned which grid cells were protected and classified cells immediately surrounding them as surrounding, all other grid cells were classified as unprotected ²⁷. Fishing mortality was then proportionally redistributed from protected cell to surrounding cells (prop), with prop defined as:

For example, if there are 4 cells surrounding a protected cell, then \({prop}\) \(=\) \(1.25\), if there are 2 surrounding cells then \({prop}\) \(=\) \(1.5.\)

Based on the computed prop, we re-estimated fishing mortality (\({f}_{{mort}}\)) (i.e., reallocate fishing effort) in the surrounding cells as:

where \({\hat{f}}_{{mort}}\) is the fishing mortality adjusted for the ‘spill-over’ of fishing effort from the protected grid cell. If the cell is not protected nor surrounding a protected grid cell then \({\hat{f}}_{{mort}}\,\)= \({f}_{{mort}}\), that is, fishing mortality will not be adjusted for the no-take MPAs effects.

Statistical analysis

We analyzed the relative contributions of no-take MPAs, fishing intensity and climate change on stock biomass and potential catches in the Northeast Atlantic region. Firstly, we analyzed the relationship between the simulated changes in biomass of the studied fish stocks, individually and aggregated across the stocks in the Northeast Atlantic region, relative to the projected global atmospheric warming levels. Secondly, we analyzed the mean responses of the biomass across the studied stocks in the Northeast Atlantic region under the various no-take MPAs, fishing and climate change scenarios using linear mixed effect models (‘lme4’ package of R). Here, the fixed effects were the area protected from fishing (no-take MPAs = 5%, 15% and 30%, as a factor), F/F_MSY (0.5, 0.75, 1, 1.5, as a factor) and global atmospheric warming levels (GWL, as a continuous variable). We set F/F_MSY = 1 (i.e., fishing at maximum sustainable rate) as the base factor in the models. To minimize the confounding effects of the algorithm that reallocated fishing effort of no-take MPAs to its surrounding area, we analyzed the effects of no-take MPAs relative to the 5% no-take MPAs scenario instead of a 0% no-take MPAs scenario. The different fish stocks were considered as random effects (1|stock). We used a multi-model approach to compare the following model structure:

Model 1: Biomass~factor(F/F_MSY) + (1|stock)

Model 2: Biomass~factor(F/F_MSY) + factor(MPA) + (1|stock)

Model 3: Biomass~factor(F/F_MSY)+ factor(MPA) + GWL + (1|stock)

Model 4: Biomass~factor(F/F_MSY)*GWL+factor(MPA) + GWL + (1|stock)

We tested the interactions between fishing and global warming level because preliminary exploration of the data indicated stocks’ sensitivity to climate change may vary at different stock sizes in the Northeast Atlantic region. We computed and compared the Akaike Information Criteria (AIC) of the models using the outputs from the climate-fish-fisheries models across all the scenarios and the ‘AIC’ function in R. We selected the model with the lowest AIC.

We also apply the same set of model structures to examine the effects of the variables and factors on catch potential. For catch, we expected that over-exploited fish stocks have lower biomass production, and no-take MPAs may have different effects on the catch of these fish stocks compared to those that are under- or fully- exploited (i.e., fishing mortality below or at MSY level). Thus, we added a model structure with a term that represented whether the stock is subjected to over-exploitation or not (i.e., fish stocks were considered being over-exploited if F/F_MSY > 1). This over-exploited term was included as a factor that interacted with MPAs.

Model 5*: Catch~factor(F/F_MSY)*GWL+factor(MPA)*factor(over-exploited) + GWL + (1|stock)

We repeated the analysis of the linear mixed effects models by changing the MPA factor from the area protected in the Northeast Atlantic (factor) to the area of the geographic range of each stock that was protected from fishing (continuous variable). In this study, we randomly assigned spatial cells as protected areas. Thus, the no-take MPAs scenarios were not designed for optimizing biomass or catch potential, or other ecological, social or economic objectives. We therefore undertook a post-hoc by calculating the proportions of distribution of the fish stocks that were within the designated no-take MPAs under each scenario and testing the effects of such proportions on the projected changes in biomass and catch. The results of such additional analysis would help understand how alternative no-take MPAs designs that protect different proportions of the stocks may contribute to biomass rebuilding under climate change.

Based on the outcomes of selected models, we can identify the relative contribution of each factor and variable, and their combinations, to changes in biomass and catch potential. Biomass and catch potential changes were expressed relative to the estimated biomass or catch at F = FMSY (i.e., dividing the estimated effects of each variable or factor by the intercept of the model).

Sensitivity analysis

We ran a sensitivity analysis by repeating the selection locations three times for each of the three no-take MPAs scenarios with F/F_MSY set to 1 under SSP RCP8.5. We applied these alternative sets of locations to a subset of species (N = 10) that broadly represent different ecology (pelagic and demersal, coastal and oceanic), life history characteristics (small-bodied/fast-growing and large-bodied/lowered growing) and taxonomy (crustaceans, molluscs, finfish and elasmobranchs) (Table S1). We ran an additional mixed effect model and examined the effects of the no-take MPAs locations (locations) on biomass:

Model (sensitivity analysis): Biomass ~ factor(MPA) + GWL + factor(locations)+(1|stock)

This study does not require ethical approval from the Research Ethic Board.