Multispecies Extensions to a Nonequilibrium Length-Based Mortality Estimator

Abstract

Recent advances in methodology allow the history of the total mortality rate experienced by a population to be estimated from periodic (e.g., annual) observations on the mean length of the population. This approach is generalized to allow data on several species that are caught together to be analyzed simultaneously based on the theory that changes in fishing effort are likely to affect several species; thus, the estimation of times when the mortality rate changes for one species borrows strength from data on other, concurrently caught species. Information theory can be used to select among models describing the degree of synchrony (if any) in mortality changes for a suite of species. This approach is illustrated using data on Puerto Rican handline fishery catches of three snapper species: Silk Snapper Lutjanus vivanus, Blackfin Snapper L. buccanella, and Vermilion Snapper Rhomboplites aurorubens. We identified the best model as the one that provided for simultaneous decreases in mortality rate around the year 1997 and for separate, species-specific magnitudes of change in total mortality. The simultaneous estimation of parameters for multiple species can provide for more credibility in the inferred mortality trends than is possible with independent estimation for each species.

Received May 30, 2016; accepted November 8, 2016

The mean length of animals in a population is a dynamic statistic that reflects the recent and past mortality rates experienced by the population (Baranov 1918, cited by Ricker 1975; Beverton and Holt 1956, 1957). Under stringent equilibrium conditions, the total instantaneous mortality rate (Z; year^–1) can be estimated by the well-known Beverton–Holt equation (Beverton and Holt 1956, 1957),

(1)

where K and L_∞ are parameters of the von Bertalanffy growth curve (K = Brody growth coefficient; L_∞ = asymptotic length), L_c is the smallest length of animals that are fully vulnerable to the fishing gear, and is the mean size of animals that are larger than L_c. Gedamke and Hoenig (2006) generalized this approach by analyzing an annual series of mean length measurements to estimate (1) the years in which the mortality rate changed and (2) the magnitudes of mortality in the separate time periods identified. To do this, Gedamke and Hoenig (2006) derived the transitional behavior of the mean length statistic as the population moved from one equilibrium state toward another equilibrium state after a change in mortality.

The present work was motivated by the consideration that most types of fishing gear catch a variety of fish species. Consequently, if fishing effort changes over time, one might reasonably expect that an assemblage or “complex” of species that are caught together would show synchronous changes in mortality rate. For example, if handline fishing effort doubles from one year to the next, one might suspect that all fish species that are caught by the handline gear would experience an increase in mortality during the same year. Data from all species in the complex could then be used to estimate the common years of change, resulting in greater efficiency in estimation. Note that all species in the complex do not necessarily experience the same time schedule of mortality changes. For example, one species may have a restricted range centered on the main fishing grounds where the increase in fishing occurred, and another species might have a broader distribution occurring on both the main fishing ground and a secondary fishing ground where effort has declined; those two species might not show synchronous changes in mortality. In addition, the magnitude of changes in mortality will not necessarily be the same for all species, as species with low catchability may exhibit a smaller change in mortality than species with high catchability.

In the present work, we sought to model changes in the mean length of several species simultaneously so as to estimate period-specific mortality rates and years of change. We develop four progressively restrictive and nested models for the estimation of mortality when multiple species are considered: (1) trends in total mortality (both the timing and magnitude of the change in mortality) are independent for each species; (2) the year in which the change in total mortality occurs (hereafter, “change point”) is common to all species, but the magnitude of the change is independent; (3) the change point is common to all species, and the magnitude of the change in fishing mortality is the product of a time and species effect; and (4) both the change point and the magnitude of the change in fishing mortality are common to all species. The latter three approaches are attractive because they borrow strength from disparate data with common underlying trends in mortality that otherwise might not be detected by analyzing each species separately. Additionally, these approaches allow us to distinguish the dynamics that are common to all species (i.e., the change points) from those that are not (i.e., the magnitude of the change in mortality). We describe our application of the models to three species that occur together in Puerto Rican handline catches: the Silk Snapper Lutjanus vivanus, Blackfin Snapper L. buccanella, and Vermilion Snapper Rhomboplites aurorubens.

METHODS

Model Development and Model Fitting

We consider a complex of N species that tend to be caught together and for which we suspect the patterns of fishing effort are similar over time. We hold the following assumptions for each species considered:

1. Individual growth follows the von Bertalanffy growth function, with the parameters K and L_∞ known and constant over time.

2. There is no individual variability in growth.

3. Recruitment is constant and continuous over time; or if recruitment fluctuates, it does so randomly (i.e., with no time trend).

4. Total mortality rate Z is constant with age after the age of recruitment (t_c) that corresponds to L_c.

These assumptions are discussed by Gedamke and Hoenig (2006) and Then et al. (2015a).

Consider first the simplest case in which a population starts at equilibrium with a total instantaneous mortality rate of and then experiences a second mortality rate , where the first subscript indexes the time period of the respective mortality rate and the second subscript indexes species n = 1, 2, . . . , N. Generalizing the results of Gedamke and Hoenig (2006) by adding a subscript for species, the mean length d years after the change point is computed as

(2)

where is the predicted mean length of animals of species n as a function of time d years after the change point; the remaining symbols are as in equation (1), with subscripts added for time period and species as needed.

Gedamke and Hoenig (2006) computed the mean length at any time when there was an arbitrary number of change points over time. Cardinale et al. (2009) analyzed length frequency data from an annual trawl survey conducted for over 100 years, and they estimated nine separate total mortality rates (and eight change points). The general equation for computing the predicted mean length in each year given a history of mortality rates is provided by Gedamke and Hoenig (2006; their Appendix 2).

From a time series of mean length observations, the mortality rates and change points can be estimated as described by Gedamke and Hoenig (2006). By virtue of the central limit theorem, we model the observed mean length of species n in year t = 1, 2, . . . T as a normally distributed random variable with mean and variance , where is the sample size of observed lengths above L_c in year t for species n. A common variance over time is assumed so that the precision of the mean length in year t depends on the sample size in that year for species n. Annual means are modeled as in equation (2) by rewriting the equation in terms of calendar years instead of the number of years since a change in mortality occurred or—more generally—as in Gedamke and Hoenig (2006; their Appendix 2), assuming more than one change point in the time series.

The log-likelihood () of observing mean length for species n over T years is proportional to

(3)

Maximum likelihood estimates of the mortality rates and change points for species n are found numerically by estimating the values of the parameters that maximize the log-likelihood. The maximum likelihood estimate for the residual variance is solved analytically, where

(4)

Modifications for Multispecies Estimation

The total log-likelihood (log_e[Λ]) for all N species is simply the sum of the individual species’ log-likelihoods,

(5)

Assuming that there are I change points, let each species n = 1, 2, . . . , N have its own vector of parameters , where is the vector of period-specific total mortality rates through time (with i = 1, 2, . . . , I + 1); and is the vector of change points in calendar years from mortality rate to (with i = 1, 2, . . . , I). We can then envision a suite of progressively more restrictive models for the patterns of mortality across species: a single-species model (SSM) and three multiple-species models (MSM1, MSM2, and MSM3), as described below.

Single-species model

In the single-species scenario, what happens to one species is not reflected in what happens to other species. For example, fishers may target certain species at certain times so that changes in total mortality for one species are independent of the changes experienced by another species. In this case, the mortality rates () and change points () are all estimated parameters. Maximizing log_e(Λ) (equation 5) is equivalent to applying the mean length estimator to each species independently because there are no parameters in common among species.

Multispecies model 1: common years of change

In MSM1, changes in fishing effort simultaneously affect all species in the complex being considered; however, the magnitudes of the changes in Z are independent. Thus, all species have a common set of change points when the mortality rate changed (i.e., for all periods i and for all pairs of species n and n’). This feature is included in the following two models, and we drop the second subscript for the change points in the subsequent text.

Multispecies model 2: common years of change with species-specific proportional changes in fishing mortality

In MSM2, all species experience synchronous changes in mortality (i.e., common change points as in MSM1), but in addition we employ separability to model changes in fishing mortality with temporal and species components. The total instantaneous mortality rate is broken down into its components,

(6)

where is the fishing mortality rate for species n during period i; and M_n is the time-invariant and age-invariant natural mortality rate for species n. A change in fishing effort in the next time period would cause the fishing mortality for a reference species to change by a factor δ. Not all species can be expected to experience the same proportional change in fishing mortality (e.g., due to different catchability coefficients), so we incorporate species-specific effects (). Subsequent mortality rates for species n are

(7)

where is the proportional change in fishing mortality for a reference species n = 1 at time D_i; and is the multiplicative species effect relative to the reference species for all other species n = 2, . . . , N.

In MSM2, the following parameters are estimated: the first total mortality rate for all N species; the residual variances for each species; the common change points D_i; the proportional changes in fishing mortality for the reference species; and the species-specific effects for all other species in the complex. Successive total mortality rates (i > 1) are derived by propagating equation (7). The values of M_n are obtained externally to the analysis of the mean length data—for example, via the methods described by Then et al. (2015b) or Hamel (2015). If there is only one change point, MSM2 is equivalent to MSM1.

Multispecies model 3: common years of change and common proportional changes in fishing mortality

The MSM3 assumes that changes in mortality are concurrent among species in the complex, with the same proportional changes in fishing mortality for all species. For all species in the complex, we modify equation (7) such that

(8)

where δ_i is the proportional change in fishing mortality at time D_i. By dropping the species subscript, the proportional change δ_i is common to all species.

In MSM3, estimated parameters include the first total mortality rate for each species, the residual variances for each species, the common change points D_i, and the corresponding δ_i; the values of M_n are again obtained externally. Successive total mortality rates (i > 1) are derived by propagating equation (8).

Model Complexity and Model Selection

Knowledge of fishing practices for species in the complex (e.g., spatial extent; depth strata fished) can be used to guide the choice of model for estimating total mortality rates. If regulations on fishing effort are implemented in different years or if there are changes in species targeting, this information may be considered when selecting the appropriate model. In the absence of such knowledge or if an empirical approach is desirable, information theory (Akaike’s information criterion corrected for small sample sizes [AIC_c]) can be used to select the best model (Burnham and Anderson 2002). Accounting for both parsimony and goodness of fit, the model with the lowest AIC_c value (AIC_c difference [ΔAIC_c] = 0) can be chosen as the most plausible model, whereas the support for other candidate models decreases with higher AIC_c values. Model selection is conditional on the values of natural mortality (M) that are specified in MSM2 and MSM3. Uncertainty in M may be addressed with a sensitivity analysis using alternative values.

By expanding the mean length analysis to incorporate multiple species and common parameters between species, the number of estimated parameters is reduced. Table 1 presents the formulas for the number of estimated parameters given I changes in mortality and N species. For example, with three species and one change in mortality, a total of 12, 10, 10, and 8 parameters would be estimated by the SSM, MSM1, MSM2, and MSM3, respectively. Incorporating an additional change point would increase the number of estimated parameters by six in the SSM but only by four in MSM1. Adding another species in the analysis would increase the number of estimated parameters by four in the SSM but only by three in MSM1.

TABLE 1. Number of estimated parameters for the single-species model (SSM) and multispecies models (MSM1, MSM2, and MSM3), where N is the number of species, I is the number of change points (years during which a change in total mortality occurred), and I + 1 is the number of estimated mortality rates. Values include the estimated residual variance for each species.

Model attribute	SSM	MSM1	MSM2	MSM3
General formula	2N(I + 1)	N + I + N(I + 1)	3N + 2I – 1	2(N + I)
One change point	4N	3N + 1	3N + 1	2N + 2
Additional parameters with
Additional species	4	3	3	2
Additional change point	2N	N + 1	2	2

Application to Deepwater Snappers in the Puerto Rican Handline Fishery

We demonstrate the application of the multispecies mean length mortality estimators by using data for the Silk Snapper, Blackfin Snapper, and Vermilion Snapper in the Puerto Rican handline fishery. The distributions of these three species are overlapping in terms of habitat and depth strata (Sylvester 1974; Boardman and Weiler 1980; Claro et al. 2001). Catches in the deepwater snapper complex have historically been dominated by Silk Snapper, followed by Vermilion Snapper and then Blackfin Snapper (Claro et al. 2001; SEDAR 2011). Currently, the three species are managed together as species complex (Snapper Unit 1) by the Caribbean Fishery Management Council (USOFR 2005). The present analysis is intended to be a demonstration of our methods and not an assessment of the three stocks.

Length data for the three deepwater snapper species from 1983 to 2013 were obtained from the commercial handline fishery through the Trip Interview Program (TIP) of the National Marine Fisheries Service’s Southeast Fisheries Science Center. Fishing occurred at depths of up to 519 m (280 fathoms), with most animals caught above 370 m (200 fathoms). The K and L_∞ parameters of the von Bertalanffy growth function for each species were obtained from the literature (Table 2). For Silk Snapper and Vermilion Snapper, there was considerable variability in the reported growth parameters. Ultimately, we used estimates from the 2011 assessment for Silk Snapper (SEDAR 2011) and estimates from Caribbean stocks for Blackfin Snapper and Vermilion Snapper (Table 2). Analyses of sensitivity to the misspecification of von Bertalanffy parameters have shown that the mean length mortality estimator is most sensitive to the overestimation or underestimation of L_∞ (Gedamke and Hoenig 2006).

TABLE 2. Von Bertalanffy growth parameters (K = Brody growth coefficient; L_∞ = asymptotic length) for the three deepwater snapper species.

Species	K	L∞ (mm)	Source
Silk Snapper	0.10	794	SEDAR 2011
Blackfin Snapper	0.10	635	Espinosa and Pozo 1982
Vermilion Snapper	0.13	532	Manickchand-Heileman and Phillip 1999

The L_c parameter was determined for each species by examining the respective length frequency data spanning the entire time period (Figure 1). In general, the length near the peak of the histogram was chosen. For Silk Snapper, the data showed an increase in the peak over time from 260 mm to 310 mm; in this case, we selected the peak of the most recent time period (i.e., 310 mm) as the L_c (Table 3).

FIGURE 1. Length frequency histograms for three deepwater snapper species captured in the Puerto Rican handline fishery from 1983 to 2013. Dashed vertical lines indicate the length at full fishery selectivity (L_c), above which the annual mean lengths were calculated for the multispecies models.

TABLE 3. Estimates of the length at full fishery selectivity (L_c), which was used to calculate mean lengths and natural mortality rates (M) for the three deepwater snapper species.

Species	Lc (mm)	M
Silk Snapper	310	0.18
Blackfin Snapper	250	0.19
Vermilion Snapper	220	0.25

The annual mean length of animals above L_c (i.e., ) and the respective sample sizes were then calculated; 7,497 records for Silk Snapper, 1,902 records for Blackfin Snapper, and 3,836 records for Vermilion Snapper were available after we removed records of lengths smaller than L_c. Sample sizes generally increased through time, although there were several intermittent years without length data for a given species (Figure 2).

FIGURE 2. Observed mean lengths (points and thin lines) and predicted mean lengths (bold lines) from the single-species model (SSM), multispecies model 1 (MSM1), and multispecies model 3 (MSM3) for Silk Snapper, Blackfin Snapper, and Vermilion Snapper. The gray shaded region indicates the 95% confidence interval of the predicted mean length from MSM1 using the derived asymptotic SEs. Concentric circles indicate the annual sample size of observed lengths (small circles = 100–249; medium circles = 250–499; large circles = 500 or more). No circles were drawn for sample sizes less than 100. The observed mean length in 1988 for Silk Snapper (514 mm from 29 samples) is not shown but was used in the analysis.

Using the SSM, MSM1, and MSM3, mortality rates were estimated by assuming that there was one change point in the time series. For Blackfin Snapper, the total mortality rate in the SSM was also estimated by assuming zero change points (equilibrium conditions), as there was equal support for zero and one change point. Because the data only suggested one change in mortality over time (Figure 2), we did not implement MSM2 due to redundancy (see Modifications for Multispecies Estimation). For MSM3, estimates of M were obtained via the method of Then et al. (2015b) by using von Bertalanffy parameters rather than maximum ages because the latter were not available (Table 3).

A sensitivity analysis of MSM3 to the prescribed M was performed. A factorial design was implemented to specify M for each of the three species at 60, 80, 100, 120, and 140% of the base values in Table 3, resulting in a total of 125 factorial combinations for all levels in all three species. The percent deviation (%DEV) in the estimated total mortality rate from each factorial combination was calculated as

(9)

where Z_sens is the estimated total mortality rate in the respective sensitivity run; and Z_base is the total mortality rate that was estimated by using the base values of M.

RESULTS

Application to Snappers in the Puerto Rican Handline Fishery

Model fits to the mean length data are shown in Figure 2. Total mortality rates were first estimated independently for each of the three snapper species under the assumption of one change in mortality (i.e., two mortality rates estimated per species) with the SSM. For all three species, the data suggested a decrease in total mortality during the observed time series (Table 4). However, the change points varied widely: the year 1996.8 for Silk Snapper, 1985.9 for Blackfin Snapper, and 1997.7 for Vermilion Snapper. Change points are estimated in continuous time with the decimal representing tenths of a year.

For Blackfin Snapper, the equilibrium Z in the data was estimated as 0.46 when no change point was specified; this value was intermediate to the two mortality rates that were estimated with one change point. Using AIC_c, there was almost equal support for the equilibrium model and the one-change-point model, but the former produced a slight trend in the residual fit. Thus, we proceeded with the analysis for Blackfin Snapper under an assumption of one change point.

TABLE 4. Estimates (SE in parentheses) of the total instantaneous mortality rate (Z) and change points (years during which a change in total mortality occurred; Z₁ = total mortality before the change point; Z₂ = total mortality after the change point) from application of the single-species model (SSM) and multispecies models 1 and 3 (MSM1 and MSM3) for the three deepwater snapper species. The proportional change in fishing mortality (i.e., δ) for MSM3 was estimated as 0.52 (SE = 0.08).

Parameter	SSM	MSM1	MSM3
Silk Snapper
Z1	0.63 (0.08)	0.62 (0.07)	0.76 (0.09)
Z2	0.49 (0.02)	0.49 (0.02)	0.49 (0.02)
Change point	1996.8 (2.93)	1997.5 (0.81)	1994.8 (2.26)
Blackfin Snapper
Z1	0.77 (0.28)	0.50 (0.05)	0.60 (0.06)
Z2	0.43 (0.03)	0.43 (0.03)	0.41 (0.03)
Change point	1985.9 (1.66)	1997.5 (0.81)	1994.8 (2.26)
Vermilion Snapper
Z1	1.89 (0.27)	1.90 (0.27)	1.39 (0.16)
Z2	0.60 (0.06)	0.61 (0.06)	0.85 (0.09)
Change point	1997.7 (0.74)	1997.5 (0.81)	1994.8 (2.26)

Next, total mortality rates were estimated by assuming a common change point using MSM1. The estimated mortality rates for Silk Snapper and Vermilion Snapper were virtually unchanged, and the estimated common change point in year 1997.5 did not noticeably depart from the change points generated by the SSM for Silk Snapper and Vermilion Snapper. For Blackfin Snapper, the first mortality rate Z₁ was estimated to be lower and the change point occurred much later in MSM1 than in the SSM. The estimated second mortality rate Z₂ for Blackfin Snapper was unchanged between the two models.

Multispecies model 3 estimated an earlier common change point (i.e., 1994.8) than was estimated by MSM1. The corresponding common value of δ was estimated to be 0.52. A decrease in mortality over time was still inferred, but the fit to the data was distinctly different from that of the SSM and MSM1 (Figure 2). For Silk Snapper and Blackfin Snapper, values of Z₂ from MSM3 were virtually unchanged compared to those from the SSM and MSM1, whereas values of Z₁ varied (Table 4). For Vermilion Snapper, Z₁ was lower, and a smaller reduction in mortality during 1994 was inferred.

Multispecies model 1 was the best-fitting model among the three, closely followed by the SSM (Table 5). There was little support for estimating common changes in fishing mortality, as the ΔAIC_c value for MSM3 was more than 10 units. Between MSM1 and SSM, there was not much support for estimating additional parameters (i.e., species-specific change points).

TABLE 5. Difference in Akaike’s information criterion corrected for small sample sizes (ΔAIC_c) from application of the single-species model (SSM) and multispecies models 1 and 3 (MSM1 and MSM3) to the three deepwater snapper species.

Model	ΔAICc	Parameters
MSM1	0.0	10
SSM	2.2	12
MSM3	11.9	8

As an indicator of model certainty and the benefit of using multiple species to infer trends in mortality, we examined the asymptotic SE of the change point in the models (Table 4). For MSM1, the asymptotic SE of the change point was 0.84, compared to 2.26 for MSM3 and a mean SE of 1.77 from the three species-specific change points in the SSM. The likelihood profile for the change point was also used to examine the contribution of the data from each species to the goodness of fit in MSM1 (Figure 3). The likelihood profile indicated that the data for Vermilion Snapper contributed the most information for estimating the change point. Together, these diagnostics suggested that the mean length data for the three species were best modeled by a common change point and independent trends in mortality.

FIGURE 3. Likelihood profile for the change point (year during which the change in total mortality occurred) from multispecies model 1 in application to the Puerto Rican deepwater snapper complex.

Sensitivity Analysis of Natural Mortality Specification

The %DEV values of the estimated total mortality rates for the three snapper species were all less than 15%, indicating that the estimates in MSM3 were not considerably affected by the specification of M (Table 6). There was no consistent trend in the estimated Z across the range of M values (Figure 4). For Silk Snapper and Blackfin Snapper, there was little to no trend in the estimation of Z₂ in the sensitivity analysis, whereas more variability and a slight trend were observed in Z₁. For Vermilion Snapper, both Z₁ and Z₂ were generally more variable with some trend given M, although the trends were in opposite directions for the two total mortality rates (Figure 4). In comparison with the SSM and MSM1, the ΔAIC_c values obtained from the sensitivity analysis for MSM3 were all greater than 6.2 and therefore did not affect model selection.

TABLE 6. Range of percent deviation (%DEV; min = minimum; max = maximum) in estimates of total mortality (Z₁ = total mortality before the change point; Z₂ = total mortality after the change point) from the sensitivity analysis of natural mortality specification in multispecies model 3 as applied to the three deepwater snapper species.

	Silk Snapper %DEV		Blackfin Snapper %DEV		Vermilion Snapper %DEV
Parameter	Min	Max	Min	Max	Min	Max
Z1	–10.9	13.8	–6.9	7.3	–6.3	10.7
Z2	–1.9	1.6	–2.2	2.7	–13.3	10.9

FIGURE 4. Estimates of the total mortality rate (Z) for the three deepwater snapper species based on the sensitivity analysis of multispecies model 3 to different specified values of natural mortality (M). The x-axis is jittered to enhance visibility of the Z-values obtained in each “bin” of M.

DISCUSSION

There is interest in developing methods for multispecies assessments, especially for ensembles of data-rich and data-poor stocks or species (Punt et al. 2011). Here, we have demonstrated the application of the mean length mortality estimator in a multispecies context.

Mean length models that incorporate common trends in mortality can use data for well-sampled species to inform trends in species with fewer observations. This can be accomplished because sample sizes are incorporated into the likelihood function to give more weight to the trends in species with more observations. In our example, the length data for Vermilion Snapper were valuable for estimating the change point. The sample sizes for Vermilion Snapper were relatively large and consistent for the duration of the time series, whereas the observations of Blackfin Snapper were historically sparser and the sample sizes for Silk Snapper were low before the year 2000.

When a common change point was estimated for all three snapper species, the trends in mortality for Blackfin Snapper and, to a lesser extent for Silk Snapper, seem to have “borrowed strength” (Punt et al. 2011) from those for Vermilion Snapper. This behavior can be diagnosed with the likelihood profile for the change point in MSM1. As a result, there was higher confidence in the timing of the common change point than in the timing of the individual change points estimated by the SSM. More credibility can be given for the change point in Blackfin Snapper to occur synchronously with those of the other two species in MSM1. The SSM did not considerably improve the goodness of fit relative to MSM1 because the estimated mortality rates for Silk Snapper and Vermilion Snapper and the most recent mortality rate for Blackfin Snapper were very similar between the two models. This was reflected in the ΔAIC_c, which did not indicate more support for the estimation of separate change points.

The data did not provide support for modeling synchronous changes in fishing mortality with MSM3, and the most noticeable change in model fit was observed with the Vermilion Snapper data. From the SSM and MSM1, total mortality was inferred to be much higher in Vermilion Snapper than in Silk Snapper and Blackfin Snapper, which implies that the changes in total mortality occurred with different proportional changes in fishing mortality in Vermilion Snapper relative to the other two species since the M values for all three species were assumed to be very similar. Estimating a common proportional change in fishing mortality altered the goodness of fit to the data. There was not enough information (i.e., more than one change point) in the data to explicitly model species-specific effects with MSM2. Thus, MSM1 provided the empirically best fit to the data in our application.

The sensitivity analysis of MSM3 in our snapper example showed that the range in Z-estimates was smaller than that specified for M. This behavior occurs because the model fit is determined by Z, and any overestimation of M is expected to be compensated for by an underestimation of fishing mortality and vice versa. Since values of M are obtained externally to the models and can vary widely depending on the empirical method used to derive those values (Hamel 2015), we recommend sensitivity analyses for applications of MSM2 or MSM3 and any possible effects on model selection.

For Vermilion Snapper, there was some uncertainty in the rather large magnitude of the estimated mortality rates. Examination of the length frequency histogram showed a truncated length structure, which could also arise from the overestimation of L_∞ specified in the model or from size-selective fishing that did not catch large animals. Identifying the cause of the truncation will require more information on the life history of the stock. As long as they are time invariant, these factors should not affect the fact that the mortality rate decreased (Gedamke and Hoenig 2006). However, knowledge of the absolute magnitude of mortality is still desirable for management purposes. This case study highlights the importance of choosing the appropriate life history parameters for mean length mortality estimators and for data-limited assessment methods in general. In our study, the growth function for Vermilion Snapper was obtained from a different stock. A high-quality growth study of the Vermilion Snapper stock in Puerto Rican waters would have conferred greater certainty of the causes underlying the large total mortality estimates in our analysis.

The advantage of using the mean length mortality estimator on a multispecies basis is to identify concomitant trends in mortality. The work we have presented can be further modified to account for particular exploitation patterns. For example, for a given species complex, certain change points can be synchronous while others can remain independent if there is information available to guide that model specification. A switch in the fishery’s target species or a regulation to reduce total effort in the fishery may produce a synchronous change in mortality, whereas the opening of a new fishery may elicit an independent change in mortality. The general framework presented here for expanding the mean length mortality estimator to multiple species by using models of varying complexity and synchrony in mortality trends allows for improved inference for one species within the context of companion species in a complex, especially if the companion species are better sampled and studied.

Selection of the Minimum Length of Vulnerability to the Fishery

The Beverton–Holt nonequilibrium methods presented here assume time-invariant, knife-edge selection in the observed length data. However, selectivity is more likely to be logistic, with partially selected small animals in the length frequency distribution. To conform to the assumption of knife-edge selection in the model, a value of L_c is chosen for use in calculating the of the truncated distribution. In a stock that is fully exploited or that has a high M/K ratio, the length distribution above L_c is monotonically descending (Figure 5.13 of Pauly 1984; Figure 7 of Hordyk et al. 2015). Thus, the choice of modal length for L_c would be appropriate, with lengths smaller than the modal length assumed to be incompletely selected. However, a stock that is lightly exploited and that has a low M/K ratio will have a modal length near L_∞, with fully selected animals comprising a significant portion of the ascending limb to the left of the mode. In this scenario, the modal length would not be suitable if the stock is exploited over a large size range; an L_c value smaller than the modal length would be more appropriate. Based on the modal length relative to L_∞ over the time series, the modal length was used as the L_c, assuming that the length distributions for all three deepwater snapper species reflected fully exploited stocks (Figure 1).

For Silk Snapper, an increase in the modal length was detected from the annual length frequency distributions in the mid-2000s, whereas the modal lengths for Blackfin Snapper and Vermilion Snapper were more stable over time (Figure 5). This shift in selectivity may have resulted from the temporary 16-in (406.4 mm) minimum size limit implemented in 2004–2006 (SEDAR 2011). The large modal lengths for Silk Snapper in 2005 and 2006 reflect this management regulation. A smaller shift in selectivity appeared to have persisted after the regulation was repealed. The mortality analysis presented here parameterized L_c to reflect the most recent selectivity pattern for Silk Snapper. The size limit regulation did not appear to affect the mortality estimation procedure, as no change in mortality was inferred to have occurred when the regulation was in effect.

FIGURE 5. Modal lengths from annual length frequency distributions for the three Puerto Rican deepwater snapper species. The dashed horizontal line in each panel shows the length at full fishery selectivity (L_c), which was used for mortality estimation.

Nonetheless, to examine the implications of a change in modal length for Silk Snapper, we analyzed the sensitivity of the mortality estimates to the chosen value of L_c by using the SSM under the assumption of one change point (Figure 6). The early mortality rate Z₁ fluctuated widely because the high values of L_c compared to the modal lengths early in the time series led to the truncation of large proportions of the length data. On the other hand, the recent mortality rate Z₂ was relatively stable when values near the recent modal length were used for L_c. For a situation in which selectivity has changed for a stock, it may be preferable to configure the mean length mortality estimator to reflect the most recent conditions so as to estimate current exploitation while allowing for greater uncertainty when inferring past trends in mortality.

FIGURE 6. Sensitivity analysis of the estimated total mortality rates (Z₁ = total mortality before the change point; Z₂ = total mortality after the change point) in Silk Snapper when different lengths at full fishery selectivity (L_c) are used.

Other Assumptions and Considerations

Improvements in fishing technology may also increase the spatial extent of exploitation of fish stocks. With serial depletion of coastal stocks, fishing effort generally moves further offshore over time as inshore areas become depleted. Changes in mean length could be a result of changes in targeting rather than changes in mortality. The ability to detect serial depletion requires high-resolution spatiotemporal data (Walters 2003; Cardinale et al. 2011). However, data on the location of catch were generally unavailable in the TIP database (SEDAR 2009). If a fishery has expanded spatially over time, one could restrict the mean length analysis to areas fished within certain depth strata or spatial strata if such data are available. If the stock actually consists of several distinct units with little interchange, then separate analyses would allow for estimation of the mortality experienced within each stratum. However, if individual fish inhabit different depths as they grow, then the observed size range in the strata would also reflect movement with age in addition to mortality. In such a scenario, effective management of the entire resource would have to rely on expert judgment to determine the extent of the spatial expansion and appropriate management measures. Analyses of additional data types, such as CPUE, can provide further insight on serial depletion (Cardinale et al. 2011).

Since the models we have examined assume that recruitment is constant, large pulses of recruitment can confound estimates of mortality. Although changes in mean length alone cannot differentiate between a recruitment pulse and a change in mortality, the length frequency distribution may provide more information. A large recruitment cohort will progress through the length structure of the catch over time; this cohort would cause the mean length to temporarily decrease. On the other hand, if there is a poor recruitment year-class, a gap in the size distribution may be observed and will also progress over time. It is important to ensure that changes in mortality are not inferred from the model when large or small year-classes are concurrently observed. Variability in recruitment may elicit small trends in mean length (Gedamke and Hoenig 2006), but model selection using AIC_c provides the balance between detection of long-term changes in mortality and spurious changes in mortality due to overfitting; the latter may be confounded with recruitment.

ACKNOWLEDGMENTS

We thank Larry Beerkircher and Josh Bennett for assistance with the data. This work was funded in part by the National Oceanic and Atmospheric Administration (Grant Number NA14OAR4170297). Q.C.H. was supported by a National Marine Fisheries Service/Sea Grant Population and Ecosystem Dynamics Fellowship (NA15OAR4170184). This is Contribution Number 3609 of the Virginia Institute of Marine Science, College of William and Mary.