Modelling ocean acidification effects with life stage-specific responses alters spatiotemporal patterns of catch and revenues of American lobster, Homarus americanus
Model speciesThe 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 US30,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 fisheries30,49. Canada’s Atlantic marine fisheries provide employment to over 40,000 people in primary harvesting and processing sectors, and supports many rural populations30.Dynamic bioclimate envelope modelWe 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 species51, and simulates future changes in biomass and maximum catch potential. Notably, the DBEM integrates growth models52, ecophysiological models53, advection–diffusion models54, and surplus production population dynamics models55 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 abundanceThe DBEM uses an initial species distribution determined using a rule-based algorithm56,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)59 was used to estimate the global abundance and distributed accordingly with relative abundance values60. 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 growthGrowth is modelled using a derived von Bertalanffy growth model to incorporate how environmental stressors affects body size of lobsters31,32,61. We model the following important life history parameters as a function of relative changes in temperature, oxygen content [O2], and pH [H+]. Growth rate (dB/dt) is dependent on weight-specific anabolic and catabolic rates:$$frac{dB}{{dt}} = H_{i,t} W^{d} – k_{i,t} W$$
(1)
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 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 Icei 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 salinities62.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 Rz 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 suitability50. 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−1 and settlement rate was 0.15 day−1—these values were chosen based on the sensitivity testing of these parameters50.Larval dispersal is modelled using an advection–diffusion54 and a larval duration model based on temperature66, such that abundance Ai,t in each cell is numerically solved for using the equation:$$frac{{partial A_{i,t} }}{partial t} = frac{partial }{partial x}left( {D_{i,t} frac{{partial A_{i,t} }}{partial x}} right) + frac{partial }{partial y}left( {D_{i,t} frac{{partial A_{i,t} }}{partial y}} right) – frac{partial }{partial y}left( {u cdot A_{i,t} } right) – frac{partial }{partial y}left( {v cdot A_{i,t} } right) – lambda cdot A_{i,t}$$
(16)
while adult dispersal is similarly modelled,$$frac{{partial A_{i,t} }}{partial t} = frac{partial }{partial x}left( {D_{i,t} frac{{partial A_{i,t} }}{partial x}} right) + frac{partial }{partial y}left( {D_{i,t} frac{{partial A_{i,t} }}{partial y}} right)$$
(17)
Advection was modelled for larval dispersal using parameters u and v for horizontal (east–west) and vertical (north–south) directions for surface current velocity (m2 s−1), respectively, between neighbouring cells x and y in the east–west and north–south direction, respectively. Instantaneous rate of larval mortality, ML, and settlement, SL was integrated into Eq. (16), where (lambda = 1 – e^{{ – left( {M_{L} + S_{L} } right)}}). The coefficient Di,t is the diffusion parameter:$$D_{i,t} = frac{{D_{i,0} cdot m}}{{1 + e^{{(tau cdot P_{i,t} cdot rho_{i,t} )}} }}$$
(18)
and$$rho_{i,t} = 1 – frac{{emptyset_{i,t} }}{{left( {C_{i,t} /overline{W}_{i,t} } right)}}$$
(19)
where Di,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)50. 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 effectsThe 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 pH32. We used parameters derived from previous experimental studies, where they observed a ~ 15% increase in mortality in larval and juvenile stages37 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 adults65. 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, lmat, 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 survivalOA 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 done32:$$Surv_{t} = Surv_{init} *left[ {1 + left( {p*left( {frac{{left[ {H^{ + } } right]_{t} }}{{left[ {H^{ + } } right]_{init} }} – 1} right)^{w} } right)} right]$$
(20)
Table 2 Scenario settings explored with model projections.Full size tableSurv 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 relationship32.Fishing pressureFishing 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 model67:$$MSY = frac{{B_{infty } cdot r}}{4}$$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, Fi,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:$$F_{MSY,i,t} = frac{r}{2}$$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 limitsFishing 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 maturity44. 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 scenariosWe 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)68. These models are included in the Coupled Model Intercomparisons Projection Phase 5 (CMIP5). We used two Representative Concentration Pathways (RCPs)69 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 scenarios2.All statistical analyses and figures were generated using the programming software R v4.0.370. More
