Modelling ocean acidification effects with life stage-specific responses alters spatiotemporal patterns of catch and revenues of American lobster, Homarus americanus

Model species

The American lobster Homarus americanus, is a crustacean found in temperate regions across the Northwest Atlantic Ocean. It is a highly valuable fishery species, caught off the coast of Eastern Canada and Northeast USA. For the past decade, they have been the most valuable single-species fishery in all of Canada and the US^30,48. In Canada, catch was estimated at 104,000 tonnes and 14% of all of Canada’s total marine fisheries catch in 2019. However, their landed value was almost $1.6 billion and 44% of the Canada’s total marine fisheries landed value, and over 49% of the landed value from Canada’s Atlantic coast fisheries^30,49. Canada’s Atlantic marine fisheries provide employment to over 40,000 people in primary harvesting and processing sectors, and supports many rural populations³⁰.

Dynamic bioclimate envelope model

We used a dynamic bioclimate envelope model (DBEM)^31,50 to assess spatial and temporal changes in the abundance and fisheries catch under different scenarios of climate change and fishing pressure. The model infers the environmental preference of the modelled species⁵¹, and simulates future changes in biomass and maximum catch potential. Notably, the DBEM integrates growth models⁵², ecophysiological models⁵³, advection–diffusion models⁵⁴, and surplus production population dynamics models⁵⁵ to determine ocean change effects on species distribution, abundance, and catch. Since H. americanus is primarily a benthic invertebrate, we used values at sea floor for environmental variables in our model. However, they also have a pelagic larval phase and we used surface environmental variables to model this life stage.

Initial distribution and abundance

The DBEM uses an initial species distribution determined using a rule-based algorithm^56,57,58. This algorithm determines species’ distribution range base on a series of geographic constraints, including latitudinal range, depth range, occurring ocean basins, and published or expert provided ‘bounding box’. It assumes that the relative abundance of the species distributes along gradients within these geographic limits with the centroid of the range having the highest relative abundance (Palomares et al. 2016). Species distributions are mapped on a global grid with a resolution of 0.5° longitude by 0.5° latitude cells with values representing relative abundance. Historical reconstructed catch data (http://www.seaaroundus.org)⁵⁹ was used to estimate the global abundance and distributed accordingly with relative abundance values⁶⁰. These catch data are based on various government and non-government reports, primary and grey literature, and also mapped on a 0.5° longitude by 0.5° latitude grid.

The initial species distribution is then overlaid on climatologies of historical environmental conditions (e.g. temperature, salinity, oxygen, pH) simulated from outputs of Earth system models (e.g. Bopp et al. 2013; Dunne et al. 2013). The DBEM assumes that species distribution is in equilibrium with the average historical environmental condition (1971–2000 average) and abundance in cells are assumed to be at carrying capacity.

Modelling individual growth

Growth is modelled using a derived von Bertalanffy growth model to incorporate how environmental stressors affects body size of lobsters^31,32,61. We model the following important life history parameters as a function of relative changes in temperature, oxygen content [O₂], and pH [H⁺]. Growth rate (dB/dt) is dependent on weight-specific anabolic and catabolic rates:

where H and k represent the coefficients for oxygen supply (anabolism) and oxygen demand for maintenance metabolism (catabolism), respectively, for cell i at time t. Body weight is scaled to anabolism with the exponent d < 1 while it is scaled linearly with catabolism. Values of d typically fall between 0.5 and 0.95 across invertebrate species, and we assume d = 0.7 for our analysis on H. americanus. Other values of d have been tested in previous studies; larger values resulted in much higher sensitivity to environmental stressors, while smaller values resulted in a minimal decrease in sensitivity^32,61. We can rearrange Eq. (1) and solve for H when growth rate, dB/dt = 0 and maximum body size W_∞ is reached, such that:

We use parameter values from SeaLifeBase (http://www.sealifebase.ca)⁶² for maximum body length, l_∞, and growth rate, K, from the von Bertalanffy growth equation, (l_{t} = l_{infty } (1 – e^{{ – K(t – t_{0} )}} )), where l is the length and t is the age in years (Table 1)⁶³. Maximum body weight, W_∞, is calculated using the length–weight conversion equation, (W = al^{b}), where a and b are coefficients also taken from SeaLifeBase (Table 1). Growth rate, K, is related to catabolic coefficient k:

Parameters H and k were estimated as a function of environmental stressors:

and

where (e^{ – j/T}) represents the Arrhenius equation to model the change in chemical reactions as a function of temperature T in degrees Kelvin. The parameters j are equal to E_a/R where E_a is the activation energy and R is the Boltzmann constant, respectively; activation energies are estimated to be 0.388 eV and 0.689 eV, based on Cheung et al. (2011), resulting in a j₁ and j₂ of 4500 K and 8000 K, respectively. Coefficients g_i and h_i are fixed parameters throughout the simulation and estimated by rearranging Eqs. (4) and (5), and substituting Eqs. (2) and (3) for H and k:

and

given initial values of maximum body size W_∞,0 and von Bertalanffy growth parameter K₀, and initial environmental conditions of temperature, oxygen concentration, and hydrogen concentration.

We measured impacts of changes in environmental conditions on growth by estimating new H_i,t and k_i,t coefficients using Eqs. (4) and (5). Equation (2) can be rearranged to solve for a new maximum body size W_∞ using new values of H_i,t and k_i,t, while Eq. (3) can be used to calculate a new growth parameter K.

To simulate change in mean body size using DBEM, we used a size transition matrix, X, to model the probabilities of an individual growing from one length class to other size classes in one time-step (year) and each grid cell where the American lobster was predicted to occur^50,64:

and

where l and l′ are the length of a particular size class and the adjacent length size classes, l_∞ is the asymptotic length, y is the age of an individual, and K is the von Bertlanffy growth parameter. Variation in growth, σ, assumed to have a coefficient of variation of 20% and is independent of length and age⁵⁰. Our model applies this general size transition model and makes no assumptions of growth stages specific to lobsters (e.g. moulting) or sex.

Mean body size (g), (overline{W}), is calculated:

where S is a relative distribution length-age frequency matrix from age class (y – 1) at size class l, and initial relative distribution at age 0 (when y = 1) across length classes was assumed to be (S_{i,t,0,l} = 1;;0;;0 ;; ldots ;;0_{{l_{infty ,i,t} }}). Parameter M is the population natural mortality, calculated from maximum body size W_∞, von Bertalanffy growth parameter K, and temperature T_Celsius (in degrees Celsius) using a model developed by⁵²:

Spawning biomass is estimated using the size transition matrix, X, and the mean weight of each size class for size classes greater than the size at maturity, l_mat⁶⁵:

Length at maturity is determined for each cell based on the maximum body size l_∞ as determined by Eq. (2) and the length–weight conversion equation.

Modelling population biomass

Population biomass is modelled using a logistic growth model⁵⁵, where biomass (B) can converted to relative abundance (A) using mean weight ((overline{W})) with the formula (A_{i,t} = B_{i,t} /overline{W}_{i,t}). Thus abundance in cell i between time t and t + 1 was estimated using:

where r is the intrinsic population growth rate of the species (Table 1), C_i,t is the carrying capacity for each cell i, (L_{j,i,t}) and (I_{j,i,t}) are the settled larvae and net migrated adults, respectively, into cell i from surrounding cells j, and F_i,t is the fishing mortality rate. Average mortality, (overline{M}), for each cell was weighted by size class specific mortality rates tested in this study. Grid cells are assumed to be at carrying capacity from the start of the simulation, and carrying capacity changes as a function of habitat suitability, P, and primary production, PP, from initial conditions (t = 0) to the current timestep, t, such that:

Habitat suitability is dependent on five environmental factors in combination with species specific traits, such that:

Habitat suitability is determined by: T is temperature (Kelvin) and TPP is the species’ temperature preference profile; Bathy is the bathymetry and MinD and MaxD is the minimum and maximum depth of the species range; Habitat_i,j is the proportion of total area of a cell with a specific habitat j (e.g. inshore, offshore, coral, estuarine, etc.); Salinity is the salinity class of the cell based on Thalassic series—metahaline (> 40 ppt), mixoeuhaline (> 29 ppt), polyhaline (> 18 ppt), mesophaline (> 5 ppt), oligohaline (> 0 ppt)—and SAssoc is the association of the species with each salinity class; and Ice_i is the sea ice % area coverage in a cell and IceP is the ice-dependency of the species. For H. americanus, they are not specifically associated with any habitat and thus only restricted by depth parameters. However, they are limited to mixoeuhaline and polyhaline salinities⁶².

The TPP was estimated using the initial predicted relative abundance (described above) overlaid with the inputs of earth system models of initial environmental conditions. The relative weight for each temperature class z of the temperature preference profile was calculated as (TPP_{z} = R_{z} /sum R_{z}), where R_z is the relative abundance in each temperature class.

A fuzzy logic model was used to model the movement between neighbouring cells based on differences in habitat suitability⁵⁰. Emigration into a cell is favoured if habitat suitability is higher than surrounding cells, and immigration out of a cell is favoured if habitat suitability is lower than surrounding cells.

We estimated larval production as 30% of spawning population biomass for each cell i, while larval mortality was 0.85 day⁻¹ and settlement rate was 0.15 day⁻¹—these values were chosen based on the sensitivity testing of these parameters⁵⁰.

Larval dispersal is modelled using an advection–diffusion⁵⁴ and a larval duration model based on temperature⁶⁶, such that abundance A_i,t in each cell is numerically solved for using the equation:

while adult dispersal is similarly modelled,

Advection was modelled for larval dispersal using parameters u and v for horizontal (east–west) and vertical (north–south) directions for surface current velocity (m² s⁻¹), respectively, between neighbouring cells x and y in the east–west and north–south direction, respectively. Instantaneous rate of larval mortality, M_L, and settlement, S_L was integrated into Eq. (16), where (lambda = 1 – e^{{ – left( {M_{L} + S_{L} } right)}}). The coefficient D_i,t is the diffusion parameter:

and

where D_i,0 is the initial diffusion coefficient and a function of the spatial grid size (GR): (D_{i,0} = left( {1.1 times 10^{4} } right) cdot GR cdot 1.33). Parameters m and (tau)—both set at 2 in the model—determine the curvature of the functional relationship between D, P, and (rho)⁵⁰. Parameter (rho_{i,t}) represents density-dependent factors and a function of population density (number of individuals), (emptyset_{i,t}), carrying capacity ((C_{i,t})), and mean body weight ((overline{W}_{i,t})) in each cell i.

Models of ocean acidification effects

The DBEM operates to model larval dispersal using advection–diffusion models. Survival is calculated at each time step (biweekly) based on a static annual survival rate. We recently tested a simple linear relationship between survival rate and pH, represented by percent changes in the survival rate given a change in pH³². We used parameters derived from previous experimental studies, where they observed a ~ 15% increase in mortality in larval and juvenile stages³⁷ from a doubling of hydrogen ion concentration.

We explore the OA effects by modelling variations in life stage-specific sensitivities to OA. Larvae, juveniles, and adults are modelled based on size classes, and the weight at maturity determines the size at which juveniles transition into adults⁶⁵. Therefore, impacts of OA on survival can be modelled for various size classes. We model the effects of OA on the three major life history stages—larval, juvenile, and adult—and use a correlative approach to link changes in ocean acidity with changes in survival.

The length transition matrix (Eqs. (8) and (9)) is split up into 40 length size classes, divided evenly from 0 to (l_{infty ,i,t}). We assume larvae transition from the pelagic phase to the growth transition matrix, and enter as the ‘larval’ stage for only the first size class. Next, juvenile size classes comprise all size classes below the length at maturity, l_mat, as determined in Eq. (12), and lobster in any size classes greater are considered adults. While our models do not incorporate lobster-specific life cycle traits (i.e. transitioning between larval stages then to juvenile stages), we use more general models that can be broadly applied to many species.

Modelling effects on survival

OA effects can be modelled as relative changes in survival rate for all life stages in Table 2. In other words, percent changes in acidity (i.e. hydrogen ion concentration) from baseline initial conditions results in a percent change in baseline survival rate. We use a model structure similar to that of previous work we have done³²:

Surv is the survival rate per year and used here as an example but can be applied to other life histories affected by OA (e.g. growth, reproduction). Survival rate in year t is derived from the initial (init) survival rate and the relative change in [H⁺] between year t and initial [H⁺] conditions. Note that in our previous model, p represents the point value of the percent change effect size with a doubling of [H⁺]. This model utilizes single point effect size estimates that have no underlying assumed relationship between acidity and survival. In our model, we used an exponent value, w, equal to 1, which assumes a linear relationship³².

Fishing pressure

Fishing mortality was assumed to be at maximum sustainable yield (MSY), which is the theoretical maximum biomass that can be sustainably removed from the population indefinitely. MSY is calculated using a Gordon Schaefer population growth model⁶⁷:

where B_∞ is the population carrying capacity and r is the intrinsic population growth rate. We use this measure of MSY as a proxy for the maximum catch potential (MCP) into the future, thus we assume that fisheries management are optimized and operate at MSY.

Furthermore, the fishing mortality rate, F_i,t—i.e. the annual proportion at which biomass is taken from the current population biomass—in each cell i at time t at MSY can be calculated as:

The fishing mortality rate was adjusted to explore scenarios of reduced fishing pressure and interactions with climate change and OA on the population dynamics of lobster. Any reductions in fishing pressure began in 2010 to represent how changes in fishing implemented now could change the state of lobster populations with the added stressors of climate change.

Fishing size limits

Fishing size limits were set to represent management scenarios and to observe their effects on lobster populations and size distribution of the population. We set four scenarios of minimum body size restrictions of lobster catch: no limit, canner small (> 0.5 lb, 220 g), canner large (> 0.75 lb, 320), and market (> 1 lb, 430 g). Canner lobsters are smaller lobsters that are often sold at a cheaper price. They range from 0.5 to 1 lb, and are largely caught in Northumberland Strait where size limits are currently set lower due to warmer waters and smaller size at maturity⁴⁴. For these scenarios of fishing size limits, we continue to assume catch is at maximum sustainable yield. Therefore, the same catch biomass (calculated using fishing mortality rate, F) will be the same for the various size limit restrictions, and more biomass will be taken from upper size classes where size limits are implemented. The no size limit results in fishing mortality to all size classes (including undersized lobsters).

Climate change scenarios

We use outputs from three Earth system models for projections of various future climate change outcomes: NOAA’s Geophysical Fluid Dynamics Laboratory (GFDL-ESM), Institute Pierre Simon Laplace Climate Modelling Centre (IPSL-ESM), and Max Planck Institute for Meteorology (MPI-ESM)⁶⁸. These models are included in the Coupled Model Intercomparisons Projection Phase 5 (CMIP5). We used two Representative Concentration Pathways (RCPs)⁶⁹ which are greenhouse gas (GHG) concentration trajectories derived to reflect possible combinations of various socioeconomic assumptions, RCP 2.6 and RCP 8.5. The number associated with RCPs represent the radiative forcing values by the year 2100 based on greenhouse gas concentrations. RCP 2.6 corresponds with a low climate change scenario and assumes immediate mitigation of GHG emissions where annual emissions peak by mid-decade (year 2025) but is reduced considerably. This scenario is more in line with the 2015 Paris Agreement to limit warming to + 1.5 °C relative to other emissions scenarios applied in most CMIP5 Earth system models. Conversely, RCP 8.5 corresponds to a high climate change scenario and our current trajectory where we continue to use fossil fuels, have little to no change to switch to renewable energy sources, and GHG emissions continue to increase with no implementation of any mitigation action. We chose these three Earth system models as they provide sea surface and bottom layers and the full range of environmental variables required by the DBEM (i.e. sea temperature, dissolved oxygen, primary production, pH, current advection, salinity, sea ice extent) for both RCP scenarios².

All statistical analyses and figures were generated using the programming software R v4.0.3⁷⁰.

Source: Ecology - nature.com