QUANTIFYING FISHERIES BALANCE WITH A NOVEL PRODUCTION ESTIMATION METHOD

Balanced harvest is a controversial Ecosystem-Based Fisheries Management (EBFM) concept conceived with intent to minimize ecosystem disruption and maximize human benefits compared to traditional management. However, most marine ecosystems lack comprehensive production estimates necessary for implementation. We developed and tested two new methods for estimating fish production at the species level with minimal data requirements. Application of our techniques to four ecological production units in the Northwest Atlantic (Mid-Atlantic Bight, Georges Bank, Gulf of Maine, and western Scotian Shelf) from 1991-2013 provided a direct estimate of 2032 kt yr-1 of total fish production. The degree of balance between catch and production distributions at the species level, assessed using application of a number of ecological indices, ranged from 0.14 to 0.91 on a scale from 0-1. Increased balance was positively associated with increased yield in the Gulf of Maine (Spearman’s, p < 0.001) but negatively associated in the Mid-Atlantic Bight (Spearman’s, p =0.045). Despite indefinite results about ecological and human impacts, we provide rare empirical exploration of balanced harvest at the species-level and outline new indicators for EBFM.


MANUSCRIPT 1
Quantifying fisheries balance with a novel production estimation method Joseph D. Zottoli 1* and Jeremy S.

INTRODUCTION
In contrast to traditional management, where single-species yields are considered individually, a number of ecosystem-based harvest strategies have been proposed with intent to reduce the negative ecosystem impacts associated with fishing and to increase total yield. These approaches include harvesting equal proportions of everything above a certain size (Larkin 1977), in proportion to each species' rate of natural predation (Fowler 1999), all harvestable species from all trophic levels in equal proportions to their production (Bundy et al. 2005), and in proportion to the productivity of all species, stocks, sexes, and sizes (Zhou et al. 2010).
The final two approaches listed above illustrate differing definitions of what has become known as balanced exploitation or balanced harvest. Balanced exploitation has garnered attention and controversy. A number of modelling studies advocate the utility of balanced harvest (Jacobsen et al. 2014, Law et al. 2014 Smith 2017) while others caution faults in feasibility, economic tradeoffs, and ecological impacts . Empirical data for assessment of balanced harvest is limited to a few freshwater case studies (Kolding & van Zweiten 2014) largely due to the requirement of extensive production (or productivity depending on how balanced harvest is defined) estimates and other uncommonly collected metrics. Bundy et al. (2005) did empirically assess the balance of the eastern Scotian Shelf with their index method but limited their assessment to balance at the trophic level and relied on static productivity estimates from an Ecopath model. Its apparent that reliable, dynamic, species-specific production estimates are essential for implementing balanced harvest at the species level, with subsequent partitions required for implementation at the level of size and sex.
Although productivity has been an important determinant of catch limits in single-species management (Ricker 1975), acquiring estimates for all species in an ecosystem, both harvested and unexploited, is daunting. Ecosystem energetics modelling has produced estimates of total fish production in ecosystems for decades beginning with simple models that utilized total primary production and trophic transfer efficiencies (Ryther 1969). These models advanced with computing power into modern end-to-end models that calculate fish production indirectly by balancing energy budgets with combinations of top-down and bottom-up processes and various assumptions (Fulton 2010). However, these modelling approaches do not allocate production by species and direct estimates of whole ecosystem fish production are rare due to the extensive information required for estimating production for each species individually. Time-varying production estimates are currently restricted to the widely accepted but data-expensive increment-summation technique or application of production-to-biomass ratios to biomass estimates (Cusson & Bourget 2005, Dolbeth et al. 2005. These approaches are limited by data availability and accuracy of production-to-biomass ratios respectively, which prevent wide-scale application across whole ecosystems.
Surplus-production models require minimal data inputs enabling ecosystemwide application. Estimates of surplus production can be obtained with just a fisheries dependent or independent index of abundance and an index of removals or effort.
However, surplus production is fundamentally different from production because it excludes biomass that dies naturally between sampling events (Figure 1). Estimation and reincorporation of this dead biomass should produce incremental estimates of production for any species given accurate estimates of biomass and mortality.
We demonstrate the effectiveness of reincorporating dead biomass to produce estimates of production using two new techniques with comparisons to real data and age-structured simulations. We also apply these techniques to dominant exploited fish and invertebrate species in four ecological production units (EPUs) to produce direct whole-ecosystem estimates of fish production. Comparison of the proportions of catch and production across EPUs with adapted ecological indices provides rare empirical assessment for one balanced harvest approach. Specifically, we relate balanced harvest, defined as harvesting all species in equal proportions of their production, to total landings and production evenness in order to assess potential benefits for yield and ecosystem health respectively.

Study Area and Data
Our study area comprised four Ecological Production Units (EPUs) as designated by the Northeast Fishery Science Center (NEFSC) located on the continental shelf along the northeastern United States (Figure 2 invertebrates. Estimates of commercial discards from the NEFSC were corrected for discard mortality and combined with commercial landings to obtain catch estimates.
We selected the species that constituted the top 95% of biomass or catch in each EPU for further analysis.

Production Estimation
Annual surplus production and biomass were estimated for all species in each EPU then again for all EPUs combined by fitting stochastic surplus-production models in continuous time (SPiCT) as described in Pedersen and Berg (2017). Informative prior estimates of the catchability parameter were applied when initial estimates exceeded one order of magnitude (<0.1 or >10) under the assumption that the prescaled survey data were reasonably accurate. Robust fits were utilized when catch or biomass indices contained outliers that greatly altered model fits. A number of species were excluded due to poor model fits, often caused by poor survey representation.
To obtain production estimates from surplus-production model fits we exploited the interrelatedness of the two terms. According to Hilborn & Walters (1992) "the term surplus production is generally used to represent the difference between production and natural mortality." It can therefore be represented as: Where SP is surplus production, P is production, t is year, and D is dead biomass from natural mortality. D was estimated by applying the conditional rate of natural mortality: Where M is the instantaneous rate of natural mortality and ) is estimated biomass at the beginning of the year, which can be acquired from fitted a SPiCT model. Using additional fitted SPiCT parameters m (maximum sustainable yield), K (carrying capacity), and n (controls shape of production curve) we can let: Eq.3 and estimate annual surplus production with: Eq. 4 Combining equations 1, 2 and 4 gives the following equation that can be used to calculate annual production and is further referred to as the SP conversion method: To estimate production without surplus-production estimates we applied the definition of production when considering solely the fate of biomass during a time period (Holme & McIntyre 1984): Equation 6 illustrates that production is equal to the difference in biomass B between sampling events plus dead biomass lost from all sources of mortality (A), which is subdivided into natural losses (eq. 2) and catch. This approach, subsequently referred to as the direct method, enables the use of a time series of biomass estimates such as those outputted from a SPiCT model fit or stock assessment, given that accurate catch and natural mortality data are also available. Constant estimates of instantaneous natural mortality were derived from the most recent stock assessment or calculated with Jensen's (1996) estimate of the second Beverton and Holt invariant using published von Bertalanffy growth parameters.

Methods Comparison
To assess the accuracy of our production estimation techniques we compared our estimates to those from the increment-summation method as described by Gillespie & Benke (1979) for eight stocks within the study area that had the necessary age-structured information readily available.
Additional assessment of our techniques was performed using simulation studies. We developed an age-structured model that simulated the population biomass, annual catch, and an annual survey index for two representative species at two levels of stochasticity (see Appendix for full set of governing equations). The large, longlived species represented a predator species while the small, short-lived species represented a forage fish or small predator. Results from the SPiCT model fitted to the simulated catch and survey indices were compared to the actual production values calculated with the increment-summation method.

Ecological Analysis
To investigate regional and temporal patterns of fishing, a number of ecological indices were applied to the species comprising the top 90% of catch or production each year from 1989 to 2015. Production was estimated using the SP conversion method except for striped bass, ocean quahog, and Atlantic surfclam, which were estimated with the direct method applied to stock assessment biomass.
The first and last two years were excluded due to missing data and unrealistic SPiCT fits at boundary years.
The evenness of catch and production were calculated with Simpson's (1949) reciprocal index corrected by the number of species N: where c is the proportion of species i in the catch and p is the proportion of total production of species i. Granger's bivariate causality test (Granger 1969) with a one- year lag was applied to catch and production evenness values to look for a predictive relationship. Preference for each species by the fishery in relation to its production was assessed with an adaption of Manly's (1972) selection index (a): In which a species with a = 1/N is neutrally selected relative to the other species. The whole-ecosystem catch-production balance was assessed by treating catch composition as a predator diet selected from the available production in the ecosystem and applying Levins' (1968) measure of niche breadth normalized for comparison between time periods and ecosystems: Niche breadth ranges between 0 and 1, with 1 representing an EPU harvested exactly in proportion to its production.
Associations among ecological indices were explored with Spearman's rank correlation to avoid assumptions about relationship structure and to minimize the influence of outliers.

Methods Comparison
Mean annual production estimates from both the SP conversion and direct methods applied to survey data were similar to those calculated with the incrementsummation method for eight stocks with age-structured assessments. There was no consistent pattern of bias for either new technique and both produced less or equal temporal variability relative to the increment-summation estimates. In some instances, the SP conversion and direct methods produced similar estimates, as exhibited by GOM haddock, summer founder, and scup ( Figure 3). This minimal discrepancy between methods resulted from relatively small and stable annual changes in biomass estimates and catch that was consistently near equal to surplus production.
The direct method applied to survey data slightly outperformed the SP conversion method but application of the direct method to stock assessment biomass was consistently most accurate overall. The Normalized Root Mean Squared Deviation (NRMSD) was lower for the direct method applied to survey data than for the SP conversion method for 5 of the 8 species but differences were typically minimal (Table 1). Noticeably lower NRMSD resulted from application of the direct method to stock assessment biomass except for butterfish. This minimal difference for butterfish may have resulted from use of spawning stock biomass rather than total biomass when applying the direct method to stock assessment output or from difficulties in estimating the high M for this species.
Production estimates from the simulation study further validated the general accuracy of our techniques in cases when the SPiCT model accurately estimated biomass (NRMSD < 0.2). The direct method performed better than the SP conversion method based on NRMSD in scenarios with moderate variability but both methods performed similarly with low variability (Table 2, Figure 4). The SP conversion method failed and produced negative estimates of production when biomass was estimated to exceed carrying capacity (Figure 4-B1).
The simulation results also demonstrated the importance of accurate natural mortality estimates and highlighted the difference between production and surplus production. Production estimates were most accurate when the applied mean adult instantaneous natural mortality (M) matched the value utilized in the simulation (0.25,0.325) in 3 of 4 situations for each method according to NRMSD compared to other M inputs (0.1,0.175,0.25,0.325,0.4). In the special case where M = 0 in the SP conversion method (eq. 5), results are equivalent to surplus production. These instances produced median estimates that were 53-64% less than actual median production and demonstrated that surplus production is not an acceptable stand-alone proxy for fish production. The severity of the difference was influenced by population size and would be exaggerated with higher real M.
In addition to varying adult natural mortality, varying juvenile mortality also determined the degrees of bias. The simulation utilized natural mortality that decreased exponentially with age. When initial juvenile natural mortality was not sufficiently high, production estimates employing the correct mean adult natural mortality were positively biased for reasons covered in the discussion section.

Ecological Analysis
Our estimates of production and exploitation revealed a number of latitudinal patterns. Estimated total mean annual production from 1991 to 2013 ranged from 201 thousand metric tons (kt) in the most northerly SS EPU to 1008 kt in the most southerly MAB EPU (Table 3) with no clear temporal trends except for an apparent decline and rebound of production in the GOM during the 2000s and decline in the MAB during the late 2000s ( Figure 5). When scaling for area, GB was the most productive EPU while the SS returned the highest yield. Production and catch were both dominated by fewer species in the more northerly EPUs (Figure 6). Relative exploitation of production was also greater in northerly EPUs with 15-40% of production harvested in most years (Figure 7).
Production-to-biomass (P:B) ratios calculated from production estimates and biomass from model outputs were not static. Although most estimates were sensible, ranging from 0.02 for GB ocean quahog to 1.87 for MAB northern shortfin squid, many species exhibited inconsistent ratios across EPUs and, to a lesser extent, through time (Table 4).
Species with the highest biomass in each EPU according to the trawl survey were not the most productive species in 76% of years across all EPUs. This confirms that comparisons of catch to production are more appropriate than comparisons to biomass. Species-level analysis showed that a number of species including butterfish, Atlantic croaker, and ocean quahog are under-selected in proportion to their relative production across EPUs while others like sea scallops and Acadian redfish are over-selected (Figure 8). The selectivity of most species varied considerably over time with some species reversing the direction of their selectivity. A few species, including spiny dogfish and pollock, were over-selected in at least one EPU and under-selected in another during the same year (Figure 8).
The balance between production and catch of all dominant species in each EPU varied regionally and temporally. Levins' Niche breadth values were highest in the SS and GOM (Table 3). Although there were no dramatic long-term changes, there was a steady decrease in niche breadth in the last 3 to 5 years of the study window in all EPUs individually but not when all EPUs were modeled together ( Figure 9). These decreases were caused by disproportionate increases in herring catch in the MAB and GB, lobster catch in the SS and both herring and lobster catch in the GOM. The SS decline in particular demonstrated the niche breadth calculation's sensitivity to changes or errors in a single species.
The relationship between catch and production evenness varied temporally and showed that catch composition is not independent of ecosystem structure in some EPUs. Spearman's rank correlation indicated that the GOM (p = 0.011) and MAB (p = 0.002) EPUs both had significant positive relationships between Simpson's evenness values for catch and production. These EPU's exhibited relatively stable evenness for catch and production but unlike the GOM and GB, catch evenness in the MAB is greater than production evenness ( Figure 10). Catch evenness was also initially higher in the SS but the relationship reversed over the study period.
There was some evidence for ecosystem response to exploitation patterns.
Granger's causality test indicated that catch evenness is a significant predictor of production evenness during the following year in 3 of the 4 EPUs at α = 0.1 (Table 5).
This suggests that ecosystem structure is sensitive to catch composition in these ecosystems. In the GOM, increased balance between catch and production produced greater total yield as demonstrated by the significant positive relationship (p < 0.001) between niche breadth (balance) and landings (Figure 11). Production was also positively associated (p = 0.004) to landings ( Figure 5). Conversely, MAB niche breadth was negatively associated with landings (p = 0.045) despite also exhibiting a positive relationship between production and landings (p < 0.001). The MAB also had a negative relationship between niche breadth and production evenness (p = 0.013) ( Figure 12). No significant balance-yield relationships were found in other EPUs.

Method Comparison and Limitations
Results demonstrate that the direct and SP conversion methods can produce relatively accurate estimates of annual production when accurate model fits are achieved, but the direct method appears superior. In addition to outperforming the SP conversion method in both simulations and comparisons to increment-summation, the direct method has greater flexibility. It can be applied to a time series of biomass and catch from a stock assessment or other source without fitting a surplus-production model. This enables utilization of biomass estimates from more complex and informed models that should, in theory, produce more accurate results. The direct method also better accounts for temporal changes in productivity. Unlike the SP conversion method that is restricted to the shape of the surplus-production curve, the direct method can produce different estimates of production for the same level of biomass, depending on the change in biomass and catch.
Both methods have room for improvement. For the purposes of this work, both approaches made the unlikely assumption that natural mortality is constant over time.
In the case of the SP conversion method, density-dependent natural mortality was accounted for by the surplus-production curve, but the reincorporation of dead biomass was based on constant M. Perhaps this M value could follow a relationship based on the shape of the surplus-production curve for future implementation.
However, the production curve is fixed through time by default in the SPiCT model, such that temporal variability in productivity and natural mortality from regime shifts, predator-prey relationships, etc. are not accurately reflected without a time-varying curve. In the case of the direct method, density dependence is not directly considered but is indirectly exhibited through differences in the incremental change of biomass, so temporal changes in productivity are incorporated. Temporal variability in natural mortality could easily be accounted for by varying the M value but as is also true with the SP method, producing reliable estimates of M is extremely difficult.
Utilizing conditional M is also problematic because it can lead to overestimation of natural deaths as noted by Ricker (1975). However, M estimates are generally restricted to represent harvested and observed ages, excluding younger ages with high associated mortality. We argue that inclusion of the high-mortality earlystages of the age spectrum in our analysis adequately negates the effects of utilizing conditional M under the assumptions that juvenile fish have higher M and contribute a sufficiently large proportion of production. This is supported by the comparable estimates to the increment-summation method using real data and by the simulations because, as mentioned in the results section, positive bias prevailed unless juvenile mortality was sufficiently high.
These methods rely heavily on accurate estimates of biomass. Careful model fitting is imperative for valid results. The SPiCT model is quite flexible and affords modelers a number of tools to obtain sensible fits. Validating biomass trends with other sources is highly recommended and can be aided with the use of confidence intervals for biomass estimates provided by SPiCT model fits. Consideration of the estimated confidence intervals for biomass also has potential for incorporation of uncertainty into estimates of production and the applied ecological indices.

Ecological Findings
Our estimates of total production of all harvestable species appear comparable but consistently lower than estimates from ecosystem energetics modelling. Au (1973) estimated with a bottom-up approach that, depending on the number of trophic steps, between 1500 and 3800 KT of strictly fish biomass was produced annually within shelf waters of the International Commission for the Northwest Atlantic Fisheries subareas 5 and 6. This region roughly covers our study EPUs except for the SS; excluding the SS we estimated annual mean production of 1830 KT. However, our estimates exclude species that did not comprise the top 95% of survey biomass or top 90% of catch and species that are poorly represented by trawl surveys such as pelagic and net-avoiding species. Conversely, we included some exploitable invertebrate species, which likely negates some of the impact of excluding fish species so, by our estimation, real fish production is likely at the lower end of Au's range. Sissenwine et al. (1983) utilized the relationship between consumption-tobiomass and P:B ratios in the GB fish community to balance an energy budget and produce estimates of 2210 and 3650 KT of annual fish production for the mid-1960s and mid-1970s respectively. These estimates far exceed our 1991-2013 mean estimate of 486 KT. Sissenwine et al.'s estimate may be elevated due to inclusion of young fish (< age 1) not represented in our study or because of overestimation of primary production. Collie et al. (2009) compiled an end-to-end energy budget that appraised GB fish production (including young fish) at a lower rate of 3.562g carbon/m 2 (1445 KT total assuming 1g carbon = 8 g wet weight) using 15% less primary production.
Other than the impacts of excluding species as mentioned previously, our estimates may be lower due to the effects of migration or misrepresented catch. It is possible that seasonal migrations of species reflected in the fisheries were not properly represented by the biannual survey, leading to over-or under-estimation of production.
GB and the SS may be particularly vulnerable to effects of migration due to their small size and in the case of the SS, largely politically defined boundary. This could also explain why relative exploitation in the SS was noticeably higher than in other EPUs.
Low estimates of relative exploitation and production can also be caused by underestimated catch (Omori et al., 2016). We may have compounded the influence of any underreported catch by excluding recreational catch in our analysis. However, most species included in the analysis experience negligible recreational fishing pressure.
Early estimates of P:B calculated using the increment-summation method applied to results of virtual population analyses by Grosslein et al. (1980) provides a unique comparison to our results for some GB species. As was true with our estimates, Grosslein et al.'s reported time-variant P:B ratios with geometric means for GB Atlantic cod (0.60), haddock (0.41), yellowtail flounder (0.63), silver hake (0.59), and Atlantic herring (0.29). Despite the three-decade time difference, our P:B ratios were comparable for the first three species (0.39, 0.63, 0.55) but much higher for silver hake (1.11) and Atlantic herring (0.88). It is surprising that herring had the lowest P:B ratio because herring are relatively small and productivity tends to increase allometrically with body size (Banse & Moser 1980) so it appears that Grosslein et al. underestimated herring P:B. The same argument can be made for silver hake but the GB ratio is 30% higher than the next largest EPU estimate and falls at the high end of our range of P:B ratios so the true ratio likely lies between 0.63 and 1.11. Regardless of the mean P:B ratios, persistent variability raises concerns about the use of static P:B ratios in other studies.
Assuming that increased yield is beneficial to humans, some empirical evidence suggests that there is potential benefit from harvesting species in proportion to relative production. The GOM provides direct evidence for human benefit because highest production and yield occurred when the catch-production balance was highest.
However, the inverse relationship was true in the MAB providing a counterexample.
Looking further, both the MAB and GOM experienced substantial declines in production and landings. In the MAB, niche breadth remained stable during the decline in landings and production resulting in the significant negative relationship, while GOM niche breadth declined in tandem with landings and production resulting in a positive relationship. In the MAB, the majority of the decline in production was explained by two species (butterfish and spiny dogfish) that were previously underutilized according to Manly's preference index whereas the decline in the GOM was mostly caused by three species, one underutilized (spiny dogfish), one overutilized (silver hake), and one harvested nearly in proportion to its production (Atlantic herring). When the production of the underutilized species in the MAB declined, this reduced the imbalance between catch and production causing niche breadth values to improve while the opposite effect occurred in the GOM as the balance increased. The decline in production of overutilized species in the GOM is logical, but the decline of underutilized species in the MAB is surprising at first glance. Other than traditional explanations such as recruitment failure or climatic cycles, the declines in underutilized species could still be caused by overharvest because of the nature of our proportion-based approach. If a few or even just one species is heavily harvested relative to its production this causes other species to appear to be under-harvested relative to their production even if they are being harvested at a sustainable level for that species. This emphasizes the importance of accurately determining the production of all major species in order to implement an all-ecosystem production management approach.
Expansion of the investigation of potential human benefits to include analyses of relationships among EPUs found further positive results. EPUs with greater niche breadth values had greater landings except for GB. Georges Bank is known to be one of the most productive marine ecosystems in the world because of its unique geographic and physical characteristics (Cohen et al. 1982) so the high landings despite low niche breadth observed there may be due to differences in available energy or habitat quality among EPUs.
Results from our study show that harvesting in proportion to production appears to have weak, negative, if any, effects on the ecosystem contrary to the minimum disturbance in trophic structure predicted by some modelling studies (Jacobsen et al. 2014, Zhou & Smith 2017. Many traditional indices used for assessing ecosystem health are not applicable to our EPU production estimates because of truncation from the selection of the dominant-species subset. It can be argued that evenness is still a viable representation of EPU health regardless of truncation. If harvesting in proportion to the production of species increased ecosystem health, one would expect a positive relationship between Levins' niche breadth values and evenness of EPU production assuming that greater evenness is representative of a healthier ecosystem. Contrary to expectations, most EPUs showed no significant trends except for the MAB, which demonstrated a negative relationship.
When expanding analysis among EPUs expectations were defied again as the two EPUs with the highest niche breadth values had the lowest production per area.
However, measuring ecosystem health is a topic of much debate and there is skepticism about the practicality of health indices (Suter 1993).
Despite some evidence of human benefit and ecosystem harm, our study region and period may not be adequate to observe the full effects of ecosystem-based harvest approaches. In most EPUs, the variability in catch-production balance was greater than any long-term trends. Intentional EBFM, rather than coincidental patterns with minimal temporal discrepancies in balance, may be required to overcome any lasting effects of the long history of heavy exploitation of the northwest Atlantic. Although there was some evidence for rapid response of ecosystem evenness to catch evenness, multi-year periods of high catch-production balance are likely required to overcome processes like variability in recruitment and oscillatory predator-prey interactions, competition etc. Additionally, our indices assess partitioning of landings and catch in this study but exclude the magnitude of total fishing pressure. Unsustainable total removals from the ecosystem could easily negate any benefits of ideal harvest patterns and appropriate community harvest must be carefully determined for implementation of EBFM.
In summary, we demonstrated that the direct method for estimating fish production is an effective approach for providing production estimates for individual species based on simulations and comparison to other works using real data.
Application of ecological indices to catch and production estimates provides a relatively simple framework for empirically assessing the degree of balance between catch and production at the species level. A case study in the Northwest Atlantic showed that increased catch-production balance was associated with increased landings in the GOM but the highly productive GB and the MAB provide regional and temporal counterexamples respectively, leading to uncertainty about potential human benefits of balancing harvest with production. The relationship between balance and production-evenness was unclear and other metrics for ecosystem health and intentional balanced harvest may be necessary for observation of ecological benefits from balanced harvest.  Table 1. Differences in normalized root mean squared deviation of production estimates between the increment-summation method and SP conversion method, direct method applied to fall survey data, and direct method applied to stock assessment biomass.   --> MAB production 0.02** MAB production --> MAB catch 0.72 GOM catch --> GOM production 0.08* GOM production --> GOM catch 0.34 SS catch --> SS production 0.08* SS production --> SS catch 0.12 GB catch --> GB production 0.95 GB production --> GB catch 0.78 Figure 1. Hypothetical annual production and annual surplus production as functions of population biomass from 0 to the carrying capacity (K). The difference between curves is due to biomass that dies naturally between sampling events (represented by the right-hand side of equation 5).      Annual relative exploitation (total catch divided by total production) in each EPU during the study period.

Figure 8.
Log-scaled Manly's preference index scores for each EPU adjusted so that 0 represents a species harvested perfectly in proportion to its relative production (selected neither for or against by the fishery).

APPENDIX
The following section walks through the simulation model structure and functions followed by tables listing symbol definitions and values utilized for this study.
Simulations were performed with R and Rstudio. The first year of all simulations began with the same starting population abundance determined by multiplying the unfished equilibrium recruitment number by the equilibrium population structure vector (l).
($,WX6) = Z * The equilibrium population structure vector comprised proportions, beginning with 1 and declining to 0 following the function for instantaneous natural mortality (M). The natural mortality function began high and declined exponentially toward an asymptote equal to adult natural mortality with a normally distributed random term representing the quality of growth and mortality in a year (Q) and an additional normally distributed random term for each year ( . ).
($,W) = _`a b -c$ + dea f W + . (9) . (9) = N(0, The equation for instantaneous fishing mortality followed a logistic form with a recursive total fishing pressure term ( n ) and an additional random term for each age class, each year ( n 5 ).
(W,$) = s 6 + m − 6 1 + -tW u n (o) + n 5 n (o) = n (o45) + n E (o) n (ovw) = 1 n 5 ~ N(0, n 5 m ) n E ~ N(0, n E m ) In order to model recruitment, we calculated a maturity ogive vector m. For simplicity, we used static maturity at age throughout each simulation. Weight at age followed the von Bertalanffy growth equation for body weight with the same quality coefficient for natural mortality that was explained previously, making growth and survival not independent.
($) = Ž (1 − -l o •$ ) • The weight at each age for the first year was randomly determined using Q = 1 and inputted into the first row of matrix .
WXZ,$ ~ N( (lX6) , ℎ WXZ (lX6) ) The weight at age for each recruitment class was randomly determined using the previous year's quality of growth assuming it would influence recruitment through parental condition or some other mechanism.
W,$XZ ~ N( (l (o45) ,$XZ) , ℎ (l (o45) ,$XZ) ) The weight at age for each cohort increased recursively following the von Bertalanffy growth equation with another random term in addition to the randomly determined year quality.
W,$ = (W,$) + (1 − -(n (o,9) p. (o,9) ) ) SPiCT models were then fit to the survey and catch indices and converted to production using the SP conversion and direct methods for comparison to real production calculated using the increment-summation method. SPiCT models that accurately predicted biomass were selected for full analysis because the goal of this study was to evaluate our new production estimation methods, not to validate the effectiveness of the SPiCT model.