Research Article, J Mar Biol Oceanogr Vol: 12 Issue: 4

A Multi-Stage, Delay Differential Model of Snow Crab Population Dynamics in the Scotian Shelf of Atlantic Canada

Jae S. Choi^*

¹Department of Fisheries and Oceans, Bedford Institute of Oceanography, Dartmouth, Nova Scotia, Canada

*Corresponding Author: Jae S. Choi,
Department of Fisheries and Oceans, Bedford Institute of Oceanography, Dartmouth, Nova Scotia, B2Y 4A2, Canada
E-mail: Jae.Choi@dfo-mpo.gc.ca

Received date: 07 July, 2023, Manuscript No. JMBO-23-105186;

Editor Assigned: 10 July, 2023, PreQC No. JMBO-23-105186 (PQ);

Reviewed Date: 24 July, 2023, QC No. JMBO-23-105186;

Revised Date: 31 July, 2023, Manuscript No. JMBO-23-105186 (R);

Published date: 07 August, 2023, 10.4172/2324-8661.1000276.

Citation: Choi JS (2023) A Multi-Stage, Delay Differential Model of Snow Crab Population Dynamics in the Scotian Shelf of Atlantic Canada J Mar Biol Oceanogr 12:4.

Abstract

Keywords:

Snowcrab dynamics; Ecosystem; Moults; Stenotherms; Zooplankton

Introduction

Snow crab (Chionoecetes opilio) is large Crustaceans, exploited primarily as a food source. They are also used as fertilizer, bait in other fisheries and for the glucosamine polysaccharide derived from chitin, known as chitosan. Chitosan is widely used in medicine as an agent to slow bleeding from wounds [1,2], agriculturally as natural fungicides and bactericides [3], plastic replacement [4] and even as a battery electrolyte [5]. In North America, including in the area of this study, the Scotian Shelf of the northwest Atlantic of Canada, the largest of males (>95 mm carapace width) are preferentially captured due to their higher meat-yield and consumer preferences of larger claws, a trait that occurs only at the terminal molt to maturity (Figure 1). This market-driven protection for the female reproductive crab which never reaches such sizes and the smaller immature crab represents a form of built-in protection for the population which spans 10 years of age or more, from larval release; it helps to offset otherwise heavy exploitation pressures of mature males by fishers with advanced technological and historical knowledge of their environment and the species’ aggregation patterns. Currently, every known population of snow crab is exploited by humans. As such, it is imperative that exploitation occurs in a responsible and sustainable manner.

Unfortunately, characterizing snow crab population dynamics is remarkably difficult because they are long lived, have a complex life history, and large vertical and horizontal shifts in spatial distributions as they grow older (ontogenetic shift in habitat) (Figure 2). Functionally, they participate in the pelagic and subsequently benthic ecosystems and associated nutrient and carbon cycles. Some of the notable life history features of snow crab that make them so interesting but difficult to model, include: sexual dimorphism with mature males being much larger than mature females; pelagic larval stages and benthic pre-adolescent and adult stages; semi-annual and annual molts depending upon size, age and environmental conditions; skipping moults if conditions are poor; terminal molt to maturity; and longevity up to 15 years. They also have a narrow range of temperature preferences [6-8]. They are thought to avoid temperatures above 7°C, as metabolic costs have been shown to exceed metabolic gains above this temperature in laboratory studies. Smaller crab and females also have differences in thermal preferenda[9]. Further, snow crab are generally observed on soft mud bottoms; with smaller-sized and molting crabs showing preference for more complex (boulder, cobble) substrates, presumably as they afford more shelter [10,11]. This life history complexity induces complexity (non-random structure) in their dynamics and spatial distributions.

In the area of study, the continental shelf of the northwest Atlantic Ocean in Figure 2 snow crabs are further exposed to high bottom temperature variability due to the confluence of a number of oceanic currents: the warm Gulf Stream, the cold Labrador current, low salinity outflow from the St. Lawrence river and the cold coastal Nova Scotia current. In this region, snow crabs are generally observed between depths of 50 m to 300 m and between temperatures of -1°C to 11°C [12]. As the focal area is thermally complex, their spatial distributions can fluctuate seasonally and annually. The additional factors of global rapid climate and ecosystem change confounds make this understanding even further.

Snow crab eggs are brooded by their mothers for up to two years, depending upon ambient temperatures as well as food availability (which also depends upon temperature as it influences overall system primary and secondary production) and the health condition and maturity status of the mother (up to 27 months in primiparous females (first breeding event) event; and up to 24 months in multiparous females (second and subsequent breeding events) [13]. More rapid development of eggs, from 12 to 18 months, has been observed in other areas [14,15]. Over 80% of the female snow crab on the Scotian Shelf is estimated to follow an annual cycle, possibly due to extended periods of exposure to warmer temperatures. A primiparous female of approximately 58 mm carapace width produces between 35,000 to 46,000 eggs, which are extruded between February and April. Multiparous females are thought to be more fecund, with more than 100,000 eggs being produced by each female. Eggs are hatched from April to June when larvae are released. This pelagic larval stage lasts for three to five months (Zoea stages I and II) during which snow crab feed upon zooplankton. Pelagic stages seem to have highest abundance in October and so may begin settling in January. Thereafter, they settle closer to the ocean bottom in their Megalopea stage. Very little is known of survival rates at these early life stages, but they are thought to depend highly upon temperature.

There have only been a small number of attempts at modeling the dynamics of snow crab populations. Being long-lived, sexually dimorphic, pelagic and benthic organisms, usually living far offshore, they are not easily surveyed. Most population assessments approach abundance estimation with a fishery or survey-based relative index of a narrow segment of the population (the fishable component). Further, aging of snow crab is not possible due to a lack of retention of calcified body parts through molts. This has encouraged adoption of size-based models [16-18], often with strong and questionable assumptions. For example, assume the male fishable component is the spawning stock biomass, though of course it is the females that are the reproductive component and completely unexploited with a shorter life span. Similarly, implicitly assumes that the fishable biomass regenerates itself by application of a biomass dynamics model [19]. Importantly, sex and size-selectivity of sampling gear are almost always assumed to be a constant or some smooth monotonic function of body size and so static across time and space or known without bias. These are significant and biologically problematic assumptions. The problems encountered are fundamental issues shared with all other attempts at assessment. Specifically, they are the interplay between snow crab life history and behavior; our inability to observe them without bias due to sampling design being non-random with respects to the environmental factors controlling their distribution and abundance; and the sampling machinery (mesh size, trawls that cannot access rugose areas) that can only observe a small fraction of the population.

This paper identifies a novel approach that overcomes some of these difficulties in modelling snow crab, in particular, by incorporating observation error across sex and stage-structure (often called a latent or state-space model) and non-stationary habitat variability to simultaneously model population dynamics and infer population parameters. We examine its utility and demonstrate the mechanism by which we can make these and more complex models operationally viable through the use of efficient and modular computing that simultaneously models and infers parameters by leveraging the Julia programming environment and the supporting Turing library for Bayesian parameter inference and the Differential Equations library for dynamical modeling [20-22]. This examination is conducted in the Scotian Shelf Ecosystem where the variability of ocean climate is known to be high and where we also have a consistently sampled population since the late 1990s [23,24]. The principal conclusion of this study is that the novel proposed model is informative and demonstrates utility in understanding snow crab dynamics, integrating all available information and versatility in incorporating/aggregating ecosystem effects through its influence upon habitat viability. It is sufficiently flexible in approach to permit other ecosystem factors such as classical inter-specific interactions and movement, which are planned for future studies.

Methodology

Data collection and subsequent index estimation are described in [25,26]. In short, sampling in an unbiased manner under complex hydrodynamic conditions is particularly difficult as the random (spatially) stratified sampling that is usually adopted to account for such variability, fail to do so, due to the large dynamic (i.e.,spatiotemporal) structures operating on scales comparable to the domain. A model-based approach Conditional Auto Regressive Spatiotemporal Models, (CARSTM), a simple extension of Generalized Linear Models (GLM) that accounts for spatiotemporal auto correlated random effects in a Bayesian context and computed with Integrated Nested Laplace Approximations (INLA), is used to address the estimation problem [27]. More specifically, biases induced by differential spatiotemporal variability in habitat for differing life stages which, upon aggregation across space, permits an informative index of abundance comparable across time (years). This study focuses upon the indices of abundance derived from this process for the period from 1999 to 2022. Importantly, no survey was conducted in 2020 due to COVID-19 related uncertainties. The index generation method being Bayesian, imputes these estimates, but variability attached to this year is large. These aggregate time series data and supporting Julia models used in this paper. The well understood phenomenological single component model of Verhulst [is the usual starting point for most models and has been used to model snow crab for many years [28-30].

It is parameterized by an intrinsic rate of increase (r) representing the maximum exponential rate of increase and carrying capacity (K) the upper threshold. When numerical abundance N → 0, the loss term also When N → K, then N → 1 and so the loss rate approaches rN and the overall rate of change approaches zero dN → 0. Parenthetically, this can be simplified further by dividing both sides by K to give By focusing upon n=N as a non-dimensional number ranging from (0, 1), this becomes For parameter estimation/inference, this form is useful as it is simpler and the magnitudes of variables are mostly similar, with the exception of K. This renders beneficial properties to numerical optimization procedures. The sigmoid nature of the equation has rendered it a frequently encountered model that shows versatility in phenomenologically describing population dynamics. Its discrete form is particularly common in fisheries settings and has also been used to describe snow crab dynamics in the area of study [31,32]. The focus is usually upon biomass B and its normalized value, b=BK-1, rather than numerical abundance: b_t+1=b_t+Δ b_t, where Δ b_t=rb_t(1-b_t)- Fishingt K-1.

Here, Fishingt is scaled by K to make it non-dimensional and also in the interval (0,1), it is usually treated as an external factor (perturbation), measured without error. Of course, there is usually observation error and/or illegal removals that will get erroneously absorbed as a biological process; in this case by elevating the intrinsic rate of increase to compensate for the additional losses. This will have an effect of creating bias and uncertainty in related parameters (see below). Note also that the discretization to an annual basis is a significant assumption as the time span is large relative to the processes of most biota. This has the effect of averaging out sub-annual dynamics. This means temporal aliasing or discretization errors and censoring are introduced which ultimately increases process and observation errors (see sub annual dynamics, below). From this phenomenological view, the minimal parameters required to estimate biological reference points, and the relative distance a system is from such reference points help delimit the status of a population. Values such as Maximum sustainable yield and the fishing mortality associated with such a yield are commonly used to help define some consistent landmarks of scale for use as management reference points to guide the implementation of a consistent precautionary approach for fishery exploitation.

Many approaches exist to estimate these model parameters. Currently, Maximum Likelihood approaches dominate due to their computational speed. However, it has been the author’s experience with this data that they do not navigate and optimize very high- dimensional parameter spaces reliably, especially in the delay differential equation models that we explore, below (Figure 3). This renders their utility in an operational setting, minimally useful.

The related Maximum A-Posteriori solutions, where using the same optimization techniques but with the addition of “prior-like” constraints can result in marginally more stable results, but they still have tremendous difficulty with high dimensional parameter spaces and the associated numerous local optima/multiple equilibria.

Here we use Bayesian inference, as informative priors for these parameters are explicit and variance propagation of index variables can be accomplished by using them as priors to respective error terms. They facilitate parameter estimates that are stable and credible, given some prior knowledge of latent processes and observation models. Previously, Just Another Gibbs Sampler (JAGS) and System Trace Audit Number (STAN) were used to compute the posteriors via Markov chain Monte Carlo (MCMC) [33,34]. The latter uses the more efficient No-U-Turn Sampler (NUTS) which significantly speeds up the estimation process of these discrete difference equation models. Presently, we use Julia programming environment and the supporting Turing library for parameter inference and the Differential Equations library for modeling delay differential equation models; they demonstrate efficiency and flexibility to simultaneously model and infer parameter values in such problems, relative to basic MCMC procedures due to use of automatic differentiation and heavily optimized solution engines.

Specifically, the latent (“real” but unobserved) biomass was assumed to have a process model error that is a recursive Gaussian or Normal distribution (N) with a standard deviation σp (Bolker 2008) such that:

Here, “~” indicates “distributed as”. Marginally informative priors were assumed: r~N(1.0, 0.1), and K~N(K, 0.1. K). The prior mean of the carrying capacity (K=[5.0, 60, 1.25], in kt, for the North-Eastern Nova Scotia (NENS), South-Eastern Nova Scotia (SENS), and Crab Fishing Area (CFA) 4x, respectively; were based upon historical analyses with other analytical procedures, namely Geostatistical Kriging with External Drift and Generalized Additive Models [35]. It is further assumed that the real unobservable (latent) process b is sampled with error (observation standard deviation σ0). The observation error model in most fishery applications is a simple multiplicative factor, often referred to as a “catchability” coefficient q. In this case we assume that: q~N (1.0, 0.1), and so the observation error model becomes:

This means that the observed biomass indices Y are some constant fraction of the true abundance bt, across all sampling time (years, season) and locations. The recursive logistic process model (eq.2) in tandem with the observation model (eq.3), we will refer to it as, “Model 1”.

From the graphical representation of the processes of Model 1 (Figure 4), one can immediately see that there is no input term from outside of the system: it is a simple phenomenological self- loop model with two outputs (death and fishing). The term rb can be referred to as the (net) “birth process” [36-38], and represents a first order processes that increases b (via birth, growth, enhanced survival, improved food availability, movement, strong year class strength, etc.), through the single (static) parameter r. The mechanism is not identified and is implied to be some internal recycling of b. In reality, there are biological mechanisms involved and associated parameters that are not-static and non-stationary (in first and second order; across time and space). Growth in biomass or increase in numbers is seasonal and pulsed, due to the molting process and body mass increases, that progress at variable rates depending upon time of year and environmental variability. Across years, there are strong and weak year classes due to match-mismatch type processes and so also pulsed [39,40]. These simple models cannot address these more realistic dynamics and so we obtain instead, an average state, with poor resolution of peaks and troughs and associated parameter estimates that are likely to be in error.

The strength of Model 1 is simplicity; but this simplicity is also a weakness. As the biomass dynamics model is phenomenological; it is heuristic, without mechanism. For example, fishable biomass “creates” fishable biomass at a maximum exponential rate (r) to a limit of K. In reality, fishable biomass being mature male snow crab, does not produce more mature male snow crab; instead, it is the female mature population that creates the eggs and broods them to larvae, some fraction of which eventually mature into the fishable component. The dynamics and longevity of these females are very different from those of the males, often showing unbalanced sex ratios and differential utilization of habitat and of course complete unfished. The rate of increase in biomass is not a simple exponential rate; instead, a certain number of fishable crab die from various causes (predation, disease, competitive fighting), a certain number of recruits terminally molt to a size in the fishable range to increase their numbers, they grow in weight as they eat, some fraction move in and out of a given area, and some are fished; these rates are not constant in time nor space and there are many lags in time. Fundamentally it is a numerical process. So Model 1 from the perspective of “biomass dynamics” is less than satisfactory in terms of biological realism. However, in the heuristic perspective, it implies that the biomass in one year is related to the biomass in the previous year within some constraints, constraints that are only diffusely/indirectly related to the “real” mechanistic processes represented by individual-individual interactions. As such, they can be seen as a temporal autoregressive model.

These constraints unfortunately result in a model that is sometimes not responsive enough to large fluctuations in dynamics caused by intrinsic and/or extrinsic factors. For example, due to the pulse-like dynamics observed since the late 1990s, no recruits entered fishable components for a number of years, even though their abundance was low. Model 1 expected recruitment simply because biomass was low. This was an erroneous expectation due to the low numbers of pre-recruits that lasted for an extended period of time; the extreme warm bottom temperatures during this period and associated shifts in their spatial distributions; the increase in abundance of predators that followed which potentially reduced the strength of that recruitment. Model 1 is too simple to express the expectations of such highly nonlinear pulse-like dynamics and interactions with environment and predators. In another example, there was an extreme warming event in 2021 that significantly altered and constricted the spatial distribution of snow crab. Again a simple model such as Model 1, with static expectations of environmental conditions (stationarity) is not able to account for such effects.

In an attempt to begin addressing these issues, we entertain Model 2 which resolves some size groups, sex and maturity and permits time-varying viable habitat area (non-stationarity) (Figure 4). Model 2 is, therefore, an intermediately complex, marginally more structured and mechanistic (ontological) model relative to Model 1.

There is approximately an 8+ year period that is required for females F_t−8 to produce the next generation of instar 9 females and males, denoted by F_t and M_5,t respectively. They represent crabs in the range of approximately 40 to 60 mm carapace width. The “birth” rates represent a combination of egg production and larval survival to instar 9+. If the population is stable, “birth” rates are expected to be similar in magnitude to overall death rates. Note that “stage” is defined based upon size-sex-morphometric traits and so misclassification is likely. The error in such knife-edge stage determinations can be substantial if growth and maturity schedules vary significantly. It is assumed, here, that in the aggregate, this has been stable in the survey record (1999 to 2022). This may of course be incorrect if persistent size/stage/age related mortality occurs due to exploitation/predation patterns and/or directional environmental change such as climate warming or food availability shifts.

Some fraction of instar 9 males transition to instar 10 (M₄), and then 11 (M₃) to 12 (M₂) and 13+(M₁); each transition rate parameterized by v, with numeric indices identifying the instar pairs. Due to the knife edge-cut of stages, there will be misclassification issues which will be absorbed by these transition rates. For example, a small fraction of instar 11 crab (M3) will be large enough to be considered a fishable size (>95 mm carapace width); but most will molt into instar 12 (M₂); this would result in a reduced molt transition rate. M2 will represent a composite group of most crab that have entered fishable size but are still immature. Some fraction of this group will mature into the fishable component in the next year; the fractional nature will result in a smaller molt transition rate. Others will continue to molt to instar 13 and higher; for our purposes, all morphologically mature males of fishable size will be considered M1 and so it represents a range of different ages. This is especially the case as terminally molted crab can live for up to another 5 years; this aggregation will effectively decrease mortality rate of the group.

Survey sampling and estimation of instar 9 to 13+ (40 mm to 130+mm carapace width) is reasonably informative. Earlier instars tend to have divergent habitat preferences relative to those of the fishable component and are known to be poorly/erratically sampled due to size approaching mesh size of sampling nets. This is because the snow crab surveys are primarily optimized for sampling the fishable component (>95 mm carapace width), both in terms of mesh size and choices of sampling location. The utility of these additional compartments is that they also permit more stable parameter estimation and forward projections that are biologically more reasonable (mechanistic). The challenge is to keep track of much more information due to the increased realism and computational limits. Each state variable u=(M₁,M₂,M₃,M₄,M₅,F), is scaled by their respective carrying capacity K., such that they are non-dimensional numbers: . Thus, for example, in the molt transition process from i=2 to j=1, the instantaneous rate is:

The multiplier is required to convert normalized density of category 2 to that of category 1 via the ratio of the respective normalizing constants, Ki. As with Model 1, first order birth/survival rates bi are assumed: βi=biui.

So far, most of the core model has been the same as Model 1, only applied to more components with some appropriate time lags. Where Model 2 diverges from Model 1 is by decomposing mortality into two components: the first component is a simple background mortality parameterized as a first order decay (ωiui,t) the second component is mortality associated with fluctuations of viable habitat, parameterized as in the logistic equation as in Model 1 as a second order process ui,t. As such, it incorporates all density dependent processes such as competition and interaction between stages. Note the congruence of this term with that of the loss term of Model 1 (see Fig. 4, eq. 1). With normalization to u, we obtain:

Here, H=V/max(V), represents the viable habitat surface area V normalized to a maximum value of 1, where the maximum is defined in the reference/focal time period. This modifies the carrying capacity proportionately and so KH can be seen as the “effective carrying capacity”, adjusting for viable habitat fluctuations. The assumption here is that when viable habitat area declines by some fraction H, so too does the effective carrying capacity, proportionately.

The inclusion of predator-prey and other ecological processes are therefore straightforward extensions as additional terms in this type of model formalism are well understood [41-44]. Such extensions are planned for future research and evaluation.

The full set of delay differential equations are, therefore:

Model 2 is, therefore, a more mechanistic (structured) extension of Model 1 that brings to bear additional information available about the abundance of other size and sex classes and transition rates between them, and mortality rates that depend upon variations of numerical density in viable habitat surface area. Fishery landings were converted to number based upon annual average weight of the fishable component estimated from surveys and discretized to a two week time step within each year. The latter was to avoid strong discontinuities which facilitates differential equation modeling and parameter inference. Survey-based indices of numerical abundance were estimated using a Hurdle process, where probability of observation and the numerical abundance of positive-valued observations were modeled as an extension of small-area analysis Conditional Auto Regressive (CAR) Models in space and time; see methods in Modeled average weights, also analyzed using the same small-area space-time analyses were used to convert numbers back to biomass to make comparisons with Model 1.

In Model 1, the meaning of q is clear, a survey catches only some (fixed) fraction of the true abundance. This is a reasonable first approximation. However, what is implied is: first, this fraction is unchanging, applying equally in high abundance and low abundance years and areas, whether they are in the mating or molting season or not. Ignored are issues such as differential net saturation, aggregation behave or, even when bottom temperatures can force aggregation and dispersal away from or into other areas. The second issue implied by q is that when a survey fails to capture snow crab, that there truly is no snow crab in the area: an observation of zero is a true zero. This is clearly false. Visual observation of trawl operation through video monitoring has shown that crabs are missed: when they are in a slight depression, when they are close to protruding rocks or on bedrock, when bottom topography is complex or simply by the gear jumping off of the bottom due to tidal currents. That is, sampling gear is biased. Furthermore, it is not just sampling gear that is problematic. Survey design is biased: in design, each areal unit is expected to be homogeneous in space and time. External factors being factored out such as bottom temperatures, food availability and predation, aggregation behavior are not homogeneous, no matter how well designed a sampling design may be. Further, there is the notion of trawlable bottom: some areas cannot be sampled without tearing or losing nets. These areas are of course not sampled and only habitats that are easier to sample will be used and consequently, over represented. Though allocation may be random in design and theory, in application, it is at best, almost random and depending upon bottom type, usually biased.

As such, the observation model is also more complex in Model 2. The assumption of zero values in survey abundance indicating a true zero is a very strong assumption. This is especially tenuous when it is known that surveys sample many size ranges very poorly. The observation model of Model 2 is, therefore, a simple linear regression. Of course, higher order terms and other covariates can and should enter into the observation model to account for potential survey and behavior induced bias. In this paper,this is mostly accomplished via the statistical abundance index model. Some additional freedom is given in the observation model in that the prior for the intercept term ci was informed by μi, the fraction of the minimum observed value relative to the maximum value.

Here, y_ti=Y_ti/max(Y_i), constrains abundance to the interval (0, 1) and helps to confer better numerical properties (faster and more stable optimization, integration) in that variables are of similar magnitudes. Additionally, the variability of the observations was propagated into Model 2 by assuming that the coefficient of variation of observations were modal estimates of the observation model using a Log Normal (LN) prior with a mode of log Coefficient Variation (CV) specific to each component i and survey time t and a Standard Deviation (SD) of 0.25: σ_ti ~ LN(log(CV_ti), 0.25).

This works as y_ti are already normalized to (0, 1). The other priors used for the Model 2 are also informative based upon expected biological constraints and some very wide distributions for the variance components to provide well mixed posterior samples as determined by Effective Sample Size (given autocorrelation in samples within chains) and Rhat (convergence criterion across chains):

Fixed reference points analogous to the concepts of Maximum Sustainable Yield (MSY) and Maximum Rate of Fishing Mortality (FMSY) are not readily identifiable in Model 2 as “production” is now dependent upon multiple system components each of which have externally imposed effective habitat viability trends and related natural and fishing mortalities (perturbations). That is, they are nonstationary, in the dynamical sense. It is, however, possible to compute a related concept, which we will call the fisheries footprint; it is computed as the difference in trajectories between the system state estimated for the fishable component with and without fishing. The biological parameters are estimated under conditions of fishing activity. In the absence of fishing, these parameters may be expected to be different, especially when fishing activity is extreme. In some cases, maturity and growth schedules can shift [45,46]. The assumption here is that extreme fishing activity has not occurred with sufficient pressure nor time to generate phenotypic or genotypic change. The difference in predicted trajectories between the fished and unfished conditions, therefore, provides a crude, first order estimate of the impact of fishing, without assumptions beyond that the model reasonably approximates reality. This fisheries footprint is placed on a relative scale from 0 to 1, by normalizing with the expected unfished abundance. As such, the fisheries footprint identifies the fraction of potential biomass that was reduced by fishery exploitation.

Results and Discussions

Estimates of abundance (biomass) of the exploited component are shown for each region for each model in Figure 5. In NENS, the discrete Model 1 showed two troughs in (pre-fishery) abundance in 2005 and 2017. The surveyed index after adjustment for the observation model tends to be lower relative to the pre-fishery abundance, which is consistent with the removal of biomass during the fishing season. The exception to this pattern was in late 2013 and 2019 when strong recruitment was also observed. Model 2 shows a similar pattern of fall fishable biomass, though with a reduced magnitude. The post 2019 period showed significant divergence relative to the survey index and especially in 2020 when no surveys were completed due to Covid-19 related disruptions. SENS also demonstrated a similar periodicity to NENS in both models In Figure 5, model 2 suggests a dampened time series, attributable to the dynamics of the recruitment being estimated to be flat. In CFA 4X, Model 1 suggests an important decline following a peak in 2010. Model 2 suggests a similar trajectory, though one that is much more variable. This area is subject to crab movement from the adjoining area as well as spatial aggregation due to extreme temperature conditions, elevated mortality from other predators and disease. Model 2 also suggests overall abundance that is much lower than the very optimistic expectations of Model 1.

It is important to note that Model 1 will project an increase in abundance, even if recruitment does not exist. It will only project a decline if abundance is above carrying capacity. This structural expectation is of course unrealistic. Model 2, which accounts for recruitment, navigates this problem a little more sensibly as an increase can only occur if recruitment exceed mortality. For example, it projects for CFA 4X, a decline, even in the absence of fishing. In SENS, abundance is projected to increase slightly in the absence of fishing. Only NENS was projected to have a rapid increase in abundance.

Another observation is that abundance may be overestimated by Model 1. For example, in CFA 4X where we have supporting information of very high natural mortality rates associated with predation and environmental variability, through the loss of strong year classes, abundance estimates are extremely optimistic. The same issues also exist for NENS, where strong adolescent crab year classes seem to disappear at a rate faster than might be expected, possibly due to predation. By extension, SENS also likely exhibits overly optimistic time trends, given the very strong year classes diminishing rapidly before entry into the fishable component. As a result of this potential overestimation of abundance, fishing mortality estimates by Model 1 may in fact be overly low. Model 2’s solutions seem more reasonable given the supporting contextual information known of the different areas. The overall relative shapes of abundance and fishing mortality are, however, similar (Figure 6).

The continuous form of model 2 also permits us to infer the sub annual dynamics of the snow crab. The saw tooth pattern of abundance suggests that during the fishing season, the rate of exploitation far surpasses the slower rate of growth of the fishable component. It also identifies how important temporal aliasing might be if observations of landings or survey sample times are not properly accounted. The fishery footprint across the different areas indicates that they are very similar to fishing mortality in overall form (as they should be, as the numerator, catch are the same). The former is on a scale that is the same as abundance and so simpler to understand than the exponential scale of an instantaneous rate. Fishery footprint also has a reasonably direct relationship to a question that is often asked in applied situations: How much of an effect is fishing activity having? In the case of the areas studied, the fishery footprint has declined to conservative and manageable levels (that is, a small fraction of the total available abundance) since the mid-2000’s. SENS, in particular, has been consistently sustainable in the historical record. In the presence of strong environmental variability, it will be necessary to continue to maintain a small fisheries footprint to reduce the susceptibility of the species such forces as has been the case for Atlantic cod and other collapsed fisheries. This is especially the case as heavy exploitation in the long-term can result in reduced reproductive success of females; loss of dominance by invasion of habitat by competitors and even genetic/phenotypic selection for reduced size at maturity [47] (Figure 7).

The first ever estimate of potential recruitment (M2) and mature female abundance (F) for the area in focus is shown in Figure 8. It shows that Model 2 has difficulty tracking these components. However, peaks and troughs are identified and they seem to be coherent across areas. The lack of concordance with the observations suggests that other, as yet un modeled processes likely need to be better accounted. The most likely missing processes include movement (inshoreoffshore, shallow-deep) as well as predation, as the areas of study are not well mixed and show spatiotemporal structure and increases in the relative abundance of predators such as Atlantic Halibut [48]. Smaller areal units would also better parameterize the spatiotemporal heterogeneity observed in the latent biological processes known to occur in the area [49] (Figure 8).

Conclusions

In studies of exploited marine animal populations, most models focus upon age structure. Most also tend to be temporally discrete (annual) difference equation approximations. This is largely due to the annual cycle in temperate regions of exploitation, and data assimilation, survey and assessment cycles and historically due to limitations to computational capacity. With such discretization’s, come various approximations and assumptions and ultimately potential error or bias from temporal (and spatial) aliasing. Space is treated as an externality or more usually, completely ignored, though of course, advection-diffusion differential equation models are readily formulated.

The continuous approach that we developed in this study addresses mechanisms in a structured manner. The additional model complexity succeeds in providing a similar and perhaps more reasonable understanding of snow crab population dynamics than the classical biomass dynamics model. It also permits a promising way forward towards describing the fishery footprint measured relative to potential, non-stationary abundance in the absence of fishing. It also identifies the importance of accounting for temporal dynamics with minimal aliasing in survey and landings, and the difficulties that may still need to be addressed: most notably the non-well-mixed nature of the areas studied show spatial and temporal structures that require further parameterization and development.

References

1. Zhang YJ, Gao B, Liu XW (2015) Topical and effective hemostatic medicines in the battlefield. Int J Clin Exp Med 8(1):10–19.

   [GoogleScholar][PubMed]

2. Moghadas B, Dashtimoghadam E, Mirzadeh H, Seidi F, Hasani-Sadrabadi MM (2016) Novel chitosan-based nanobiohybrid membranes for wound dressing applications. Rsc Advances. 6(10):7701-7711.

   [Crossref][Google Scholar]

3. Stoner RJ (2007) Pre-harvest application of proprietary elicitor delays fruit senescence. InAdvances in Plant Ethylene Research: Proceedings of the 7th International Symposium on the Plant Hormone Ethylene 301-302. Dordrecht: Springer Netherlands.

   [Crossref][Google Scholar]

4. Tampieri A, Celotti G, Landi E, Sandri M, Roveri N, et al. (2003) Biologically inspired synthesis of bone‐like composite: Self‐assembled collagen fibers/hydroxyapatite nanocrystals. Biomed Mater Res A 67(2):618-625.

   [Crossref][Google Scholar][PubMed]

5. Poosapati A, Vadnala S, Negrete K, Lan Y, Hutchison J, et al. (2021) Rechargeable zinc-electrolytic manganese dioxide (EMD) battery with a flexible chitosan-alkaline electrolyte. ACS Appl Mater 4(4):4248-4258.

   [Crossref][Google Scholar]

6. Foyle TP, O’Dor RK, Elner RW (1989) Energetically defining the thermal limits of the snow crab. J. Exp. Biol 145(1):371-393.

   [Crossref][Google Scholar]

7. Sainte-Marie B, Lafrance M (2002) Growth and survival of recently settled snow crab Chionoecetes opilio in relation to intra-and intercohort competition and cannibalism: a laboratory study. Mar Ecol Prog Ser 244:191-203.

   [Crossref][Google Scholar]

8. Kuhn PS, Choi JS (2011) Influence of temperature on embryo developmental cycles and mortality of female Chionoecetes opilio (snow crab) on the Scotian Shelf, Canada. Fish Res 107(1-3):245-252.

   [Crossref][Google Scholar]

9. Choi JS, Cameron B, Christie K, Glass A, MacEachern E (2022) Temperature and depth dependence of the spatial distribution of snow crab. en-CA. In: In prep

   [Crossref][Google Scholar]

10. Sainte-Marie B, Hazel F (1992) Moulting and mating of snow crabs, Chionoecetes opilio (O. Fabricius), in shallow waters of the northwestern Gulf of Saint Lawrence. Can J Fish Aquat Sci 49(7):1282-1293.

    [Crossref][GoogleScholar]

11. Comeau M, Conan GY, Maynou F, Robichaud G, Therriault JC, et al. (1998) Growth, spatial distribution, and abundance of benthic stages of the snow crab (Chionoecetes opilio) in Bonne Bay, Newfoundland, Canada. Can J Fish Aquat Sci 55(1):262-279.

    [Crossref][Google Scholar]

12. Choi JS (2009) Habitat preferences of the snow crab, Chionoecetes opilio: where stock assessment and ecology intersect. Biol Manag of Exploit 10:361.

    [Crossref][Google Scholar]

13. Sainte-Marie B (1993) Reproductive cycle and fecundity of primiparous and multiparous female snow crab, Chionoecetes opilio, in the northwest Gulf of Saint Lawrence. Can J Fish Aquat Sci. 50(10):2147-2156.

    [Crossref][Google Scholar]

14. Elner RW, Beninger PG (1995) Multiple reproductive strategies in snow crab, Chionoecetes opilio: physiological pathways and behavioral plasticity. J ExpMarBiolEcol 193(1-2):93-112.

    [Crossref][Google Scholar]

15. Webb JB, Eckert GL, Shirley TC, Tamone SL (2007) Changes in embryonic development and hatching in Chionoecetes opilio (snow crab) with variation in incubation temperature. Biol Bull  213(1):67-75.

    [Crossref][Google Scholar]

16. Siddeek MS, Wade E, Moriyasu M (2009) Proposed harvest control rule for managing the snow crab stock in the southwestern Gulf of St. Lawrence, Canada. Fish res 95(2-3):268-279.

    [Crossref][Google Scholar]

17. Cadigan NG, Wade E, Nielsen A (2017) A spatiotemporal model for snow crab (Chionoecetes opilio) stock size in the southern Gulf of St. Lawrence. Can J Fish Aquat Sci 74(11):1808-1820.

    [Crossref][Google Scholar]

18. Ianelli J, Punt A, Webber D, Martell S, Szuwalski C (2022) General Model for Assessing Crustacean Stocks (GMACs).
19. Cook AM, Zisserson BM, Cameron BJ, Choi JS (2015) Assessment of Scotian shelf snow crab in 2014. CSAS

    [Google Scholar]

20. Bezanson J, Edelman A, Karpinski S, Shah VB (2017) Julia: A fresh approach to numerical computing. SIAM Rev. 59(1):65-98.

    [Crossref][Google Scholar]

21. Ge H, Xu K, Ghahramani Z (2018) Turing: a language for flexible probabilistic inference. Int. Conf Artif intel stat 1682-1690.

    [Google Scholar]

22. Rackauckas C, Nie Q (2017) Differential equations. jl–a performant and feature-rich ecosystem for solving differential equations in julia. J Open Res Softw 5(1).

    [Crossref][Google Scholar]

23. Choi JS, Frank KT, Leggett WC, Drinkwater K (2004) Transition to an alternate state in a continental shelf ecosystem. Can J Fish Aquat Sci 61(4):505-10.

    [Crossref][Google Scholar]

24. Drinkwater KF (2005) The response of Atlantic cod (Gadus morhua) to future climate change. Ices journal of marine science 62(7):1327-1337.

    [Crossref][Google Scholar]

25. J.S. Choi (2020) “A framework for the assessment of snow crab (Chionoecetes opilio) in Maritimes Region (NAFO div 4VWX)”. Can Sci Advis Sec Res Doc 68
26. Choi JS (2022) Reconstructing the decline of Atlantic Cod with the help of environmental variability in the Scotian Shelf of Canada.

    [Crossref][Google Scholar][

27. Rue H, Riebler A, Sørbye SH, Illian JB, Simpson DP, et al. (2017). Bayesian computing with INLA: a review. Annu Rev Stat Appl. 4:395-421.

    [Google Scholar]

28. Verhulst PF (1845) Resherches mathematiques sur la loi d'accroissement de la population. Nouv MemAcadRSci 18:1-41.

    [Google Scholar]

29. Bacaër N (1838) Verhulst and the logistic equation. A short history of mathematical population dynamics. 35-39.

    [Crossref][Google Scholar]

30. Smith SJ, Bourdages H, Choi J, Dawe E, Dunham JS, et al. (2012) Technical guidelines for the provision of scientific advice on the precautionary approach for Canadian fish stocks: Section 7–invertebrate species. DFO Can Sci Advis Sec Res Doc 117.

    [Google Scholar]

31. Schaefer MB (1954) Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. IATTC 1

    [Crossref][Google Scholar]

32. Meyer R, Millar RB (1999) BUGS in Bayesian stock assessments. Can J Fish Aquat Sci 56(6):1078-1087.

    [Crossref][Google Scholar]

33. Plummer M (2003) JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. 124(10):1-10.

    [Google Scholar]

34. Stan Development Team (2022) Stan Modeling Language Users Guide and Reference Manual, v 2.30.1.
35. Choi JS, Zisserson BM, Reeves AR (2005) An assessment of the 2004 snow crab populations resident on the Scotian Shelf (CFAs 20 to 24). Fish Ocean Can Sci

    [Google Scholar]

36.  Gillespie DT (1992) A rigorous derivation of the chemical master equation. Phys. A: Stat. Mech 188(3):404-425.

    [Crossref][Google Scholar]

37. Qian H, Bishop LM (2010) The chemical master equation approach to nonequilibrium steady-state of open biochemical systems: Linear single-molecule enzyme kinetics and nonlinear biochemical reaction networks. Int J Mol Sci 11(9):3472-3500.

    [Crossref][Google Scholar][PubMed]

38. Jaroslav Albert (2016) “A Hybrid of the Chemical Master Equation and the Gillespie Algorithm for Albert J. A hybrid of the chemical master equation and the Gillespie algorithm for efficient stochastic simulations of sub-networks. PloS one 11(3):e0149909.

    [Crossref][Google Scholar][PubMed]

39. Cushing DH (1969) The regularity of the spawning season of some fishes. ICES J Mar Sci 33(1):81-92.

    [Crossref][Google Scholar]

40. Durant JM, Hjermann DØ, Ottersen G, Stenseth NC (2007) Climate and the match or mismatch between predator requirements and resource availability. Clim Res 33(3):271-283.

    [Crossref][Google Scholar]

41. Lotka AJ (1920) Analytical note on certain rhythmic relations in organic systems. Roc Natl Acad Sci 6(7):410-415.

    [Crossref][Google Scholar][PubMed]

42. Von Bertalanffy L (1950) The theory of open systems in physics and biology. Science. 111(2872):23-29.

    [Crossref][Google Scholar][PubMed]

43. Pianka ER (1970) On r-and K-selection. Am Nat. 104(940):592-597.

    [Google Scholar]

44. Wangersky PJ (1978) Lotka-Volterra population models. Annu Rev Ecol Syst 9(1):189-218.

    [Crossref][Google Scholar]

45. Zwanenburg KC, Bowen D, Bundy A, Frank K, Drinkwater K, et al. (2002) 4 Decadal changes in the Scotian Shelf large marine ecosystem. LMEs 10:105-150.

    [Crossref][Google Scholar]

46. Choi JS, Frank KT, Petrie BD Leggett WC (2005) “Integrated assessment of a large marine ecosystem: a case study of the devolution of the eastern Scotian Shelf, Canada”. Oceanogr. Mar. Biol Ann Rev 43:47–67.

    [Crossref][Google Scholar]

47. Comeau M, Conan GY (1992) Morphometry and gonad maturity of male snow crab, Chionoecetes opilio. Can J Fish Aquat Sci. 49(12):2460-2468.

    [Crossref][Google Scholar]

48. DFO “Stock Status Update of Atlantic Halibut (Hippoglossus hippoglossus) on the Scotian Shelf and Southern Grand Banks in NAFO Divisions 3NOPs4VWX5Zc for 2020”. In: DFO Can Sci.Advis Sec Sci Resp.
49. Quinn TJ (2003) Ruminations on the development and future of population dynamics models in fisheries. Nat Resour 16(4):341-392.

    [Crossref][Google Scholar]