Assessing pinniped bycatch mortality with uncertainty in abundance and post-release mortality: A case study from Chile

The effects of human-caused mortality, such as fisheries bycatch, of endangered, threatened and protected (ETP) species of marine mammals can be evaluated using population model-based stock assessments. The information available to conduct such assessments is often very limited. Available data might include fragmented time-series of abundance estimates, incomplete data on bycatch for the fisheries that interact with ETP species (often few years and low observer coverage), and perhaps some data on scale and trends in fishing effort. Such data are challenging to use as the basis for stock assessments, which generally assume that estimates of removals (bycatch, in our context) through time are available for at least the most recent decade or two. This paper de- scribes a stock assessment method for use with sparse observer data on bycatch mortality, applied within the context of a Bayesian estimation framework. The method produces estimates, with associated measures of precision, of population size and historical time-series of bycatch mortality that are consistent with the observer and abundance data. It provides a rigorous way to account for the uncertainty arising from animals that are caught but released alive and then die subsequent to release, given a post-release mortality rate prior. Observer data from industrial and artisanal purse seine and trawl fisheries and survey data for South American sea lions ( Otaria byronia ) and South American fur seals ( Arctocephalus australis ) off Chile are used to illustrate the method.


Introduction
The main sources of human-caused mortality for many endangered, threatened and protected (ETP) marine species, including marine mammals, are entanglement, entrapment, and hooking in commercial fishing gear (Read, 2005;Reeves et al., 2013), referred to herein collectively as bycatch. Although there are many other important causes of marine mammal mortality, there is a general consensus among scientists and conservationists that bycatch is the most certain and potent driver of human-caused population declines, and the primary barrier to recovery of some depleted populations (Gales et al., 2003;Kovacs et al., 2012;Reeves et al., 2013;Lewison et al., 2014;Avila et al., 2018).
Bycatch management programs require estimates of bycatch mortality and assessment of its effect on the long-term status of ETP populations.
The Fish and Fish Product Import Provisions of the USA Marine Mammal Protection Act (MMPA) (Federal Register, 2020) require that a harvesting nation's regulatory programs address intentional and incidental (bycatch) mortality and serious injury to marine mammals in fisheries that export seafood to the USA. One way for a country to achieve a 'comparability finding' under the associated regulations is to demonstrate that bycatch levels relative to population size are either consistent with maintaining stocks at their optimal sustainable population (OSP) level, or are expected to recover them to that level (OSP has been defined operationally in USA regulations as a population size between the Maximum Net Productivity Level (MNPL), the number of animals at which production is maximised, and the carrying capacity (K) (Wade, 1998)). Methodologies that work with short time-series of abundance and bycatch data will be essential in cases where countries impose new or expanded regulatory programs to address marine mammal bycatch for the purposes of obtaining a comparability finding (Williams et al., 2016).
Several methods have been developed to provide management advice for conserving marine mammal populations subject to humancaused mortality. Assessments based on fitting population dynamics models are conducted by the Scientific Committee of the International Whaling Commission (IWC SC) for several baleen whale populations (e. g., humpback whales: International Whaling Commission (IWC), 2007; gray whales : Punt and Wade, 2012;International Whaling Commission (IWC), 2017; Antarctic minke whales: Punt et al., 2014), by the Scientific Committee of the North Atlantic Marine Mammal Commission (NAMMCO) for regional cetacean and pinniped populations (International Council for the Exploration of the Sea (ICES), 2016,2016,2017,2019), and by some national jurisdictions for a variety of marine mammal populations (e.g., Laake et al., 2018). These assessments provide estimates of current population size, status of the population relative to reference points such as K, MNPL or maximum sustainable yield level, MSYL (Punt, 2017). In addition, assessments can be used to evaluate the performance of management procedures intended to ensure that removals are consistent with management goals (e.g., Wade, 1998;International Whaling Commission (IWC), 2003, (IWC), 2012Punt and Donovan, 2007) The primary data requirements for most population model-based stock assessments are a time-series of removal data, information on biological parameters such as age-specific natural survival rates and ageat-maturity, and one or more indices of absolute or relative abundance. Population models for marine mammals usually assume that humancaused removals are known with little error. This assumption may be realistic for whaling and sealing operations where there is an effective monitoring program (but see Ivashchenko et al., 2011), but not for most marine mammal populations. Because of this, methods have been developed to quantify the implications of uncertainty about levels of removals from marine mammal populations in relation to estimates of current abundance and projections of future abundance (e.g., Givens and Thompson, 1996;Zerbini et al., 2019).
There are several ways to estimate marine mammal bycatch rates in specific fisheries, but best practices involve placing trained observers on fishing vessels to collect information on the species identity and fate of bycaught individuals, and combining this information with estimates of the total effort in the fishery, the latter usually based on logbook and other monitoring programs (Cox et al., 2007;Curtis and Carretta, 2020). Accurate estimation of the bycatch mortality rate for a population also often requires bycatch data from multiple fisheries that affect the same marine mammal population, as well as correction factors to account for spatial or temporal variability in or constraints on monitoring effort. Some fraction of the bycaught animals may be released alive and then die later (as a result of injuries or encumbrance by fishing gear) without being recorded as killed (e.g., Wilson et al., 2014). This form of so-called 'cryptic mortality' (Williams et al., 2011) must also be considered in population assessments. This paper develops a way to include data from observer programs directly into a stock assessment that estimates bycatch mortality and its uncertainty from relatively short data time series (abundance information for 10-20 years and bycatch information for 3-5 years). This approach can use bycatch estimates from an observer program along with associated measures of effort collected from a subset of the fisheries. The approach accounts for the uncertainty associated with both monitoring operations and the impact of post-release mortality within a Bayesian estimation framework. We illustrate our novel method by applying it to South American sea lions (Otaria byronia, SASL) and South American fur seals (Arctocephalus australis, SAFS) in zones off Chile, using short time-series of bycatch data from multiple fisheries and information from an observer program with variable coverage, and we include consideration of hypotheses related to post-release mortality.

Background to Chilean fisheries, SASL and SAFS
Chile is one of the major producers of fish and fish products in the world, with around 1.6 million tonnes landed in 2017 (FAO, 2018). In particular, Chile, along with Peru, is the main exporter of anchoveta (Engraulis ringens), the species producing the second-largest catch by weight in the world (FAO, 2018). More than 15,500 vessels participated in all Chilean fisheries during 2016, ranging from artisanal small boats up to large industrial ships (FAO, 2018). The artisanal fleet is heterogeneous, ranging from deckless boats (< 8− 10 m in length) to small and mid-sized vessels with a maximum length of 18 m. The industrial fleet includes vessels larger than 50 gross tonnes and >18 m length (Castilla, 2010).
The largest industrial fisheries in Chile can be divided into a demersal fishery that uses bottom trawls and longlines, and a pelagic fishery that uses mainly purse seines. The demersal fishery operates mainly in central and southern Chile. The primary target species for bottom trawls are South Pacific hake (Merluccius gayi), Patagonian grenadier (Macruronus magellanicus), southern hake (Merluccius australis) and Humboldt squid (Dosidicus gigas). Longlines are used in southern Chile to catch Patagonian toothfish (Dissostichus eleginoides). The pelagic purse seine fishery mainly operates in northern and central Chile, and targets primarily anchoveta, herring (Clupea bentincki) and Chilean jack mackerel (Trachurus murphyi). Bycatch of several species of marine mammals including SASL and SAFS has been documented in these fisheries (e.g. Hückstädt and Antezana, 2003;Reyes et al., 2012;González-But and Sepúlveda, 2016).
The abundance of SASL and SAFS is monitored based on haul-out counts (Venegas et al., 2001;Bartheld et al., 2008;Oliva et al., 2012Oliva et al., , 2020. There are two putative subspecies of SAFS, one in Peru and northern Chile (unnamed subspecies) and the other in southern Chile and the South Atlantic 1 (A. a. australis; Oliveira and Brownell, 2014 Venegas et al., 2002;Oliva et al., 2020). Three macrozones for SASL off Chile are recognized: northern, central and southern (Oliva et al., 2020;Fig. 1), with evidence of population structure within these units [stocks] (Weinberger, 2013;Oliveira et al., 2017) as well as different trends in abundance among units (Oliva et al., 2020). The government of Chile enacted Law No. 20,625 in September 2012, which amended the General Fisheries and Aquaculture Law and introduced the concept of incidental catch (or bycatch) of non-target sea turtles, seabirds and marine mammals, along with control measures and 1 The two subspecies have been proposed but not yet formally described and named; see https://marinemammalscience.org/species-information/list-mari ne-mammal-species-subspecies/list-of-proposed-un-named-marine-mamm al-species-and-subspecies/ sanctions for fisheries with bycatch. These changes were meant to provide legal recognition of bycatch, reduce bycatch levels in Chilean fisheries, and increase awareness of the need to incorporate bycatch into the management of fisheries under an ecosystem approach. Among other mandates, this law required the development of research programs to gather data for establishing fishery-or gear-based bycatch reduction plans. Since 2013, data have been collected by trained scientific observers from the Instituto de Fomento Pesquero (IFOP; Fisheries Development Institute) onboard commercial fishing vessels (Domenech, 2016;Supplementary Table 2). This has included total bycatch by species, determinations of causes of bycatch (i.e., the type of interaction with the fishing gear), and descriptions of fishing practices. The program is voluntary, and the proportion of vessels being monitored by scientific observers varies markedly among fisheries, with <5% coverage for some fisheries (e.g., pelagic purse seining) and >70 % for others (e.g., bottom trawling for South Pacific hake).
A particular challenge for assessing and managing bycatch of SASL and SAFS off Chile is that while there are estimates of abundance  Table 3). As described above, these estimates are based on observer program coverage ranging from <5% to >70 % (Supplementary Table 3), and the voluntary nature of participation by fishers in the observer program could introduce bias. A further challenge for accurately estimating bycatch mortality is that some of the animals that are captured in the gear (i.e., are bycaught and may be recorded as such by observers) escape or are released alive but may eventually die as a result of their injuries (see review in Wilson et al., 2014) or stress-related changes in behaviour (e.g., separation of pups from adults; Edwards, 2002;Noren and Edwards, 2007); the rate of this post-release mortality is unknown.

Overview
The assessment method presented here is based on fitting a population dynamics model to the abundance estimates (estimates of abundance for animals age 1 and older, and indices of pup numbers), as well as to data on numbers bycaught and bycatch mortality, along with a prior for the post-release mortality by fishery. The removals from the model are predicted by fishery, accounting for (postulated) trends in fishing effort and levels of bycatch mortality by fishery, accounting for post-release mortality. The modelled population is projected from 1990 to 2020 given the calculated numbers dying annually due to bycatch mortality to allow for parameter estimation, and then forward into the future under a scenario of no future human-caused mortality to assess how rapidly the population can rebuild, and to compute a posterior distribution for K. The results of the assessment are used to apply the potential biological removal (PBR, a conservation or sustainability reference point used in the USA system) approach to assess whether current levels of bycatch mortality are consistent with the policy goal of recovering populations to their MNPLs and maintaining them at or above those levels.

Population model
SASL and SAFS off Chile are described using age-and sex-structured models, with separate models for each of three zones (north, central and south for SASL, south for the SAFS) (Fig. 1). The model (see Table 1 for equations and Table 2 for parameter and variable definitions) assumes that the populations in each zone are closed, so no allowance is made for movement among zones. The number of pups produced each year (sex ratio 50:50; Crespo, 1980;Cappozzo and Perrin, 2009) depends on the number of females that have reached the age of first parturition (Eqn T1.2) and a density-dependent birth rate, with the extent of density dependence being a function of the abundance of animals age 1 and older (age 1+), relative to the carrying capacity (Eqn T1.4). Each population is assumed to be in a steady state at the start of the first year considered in the model (1990, a year sufficiently far in the past that the population will be out of steady state once abundance data are available), subject to a human-caused mortality rate on all age 1+ animals. This rate, and hence the number of animals by age and sex at the start of 1990, is computed based on the number of age 1+ animals relative to carrying capacity (International Whaling Commission (IWC), 2017). Human-caused mortality by age and sex from 1990 onwards, assumed to take place at the start of each year, for computational simplicity, is the sum of the removals by age and sex over all fisheries (Eqn T1.6). The human-caused mortality is the sum of the number of animals that are captured by each fishery (Q s,f t,a ) and die immediately and the number that are captured and survive initially (escape or are released) but die post-release (Eqn T1.9). The value of Q s,f t,a depends on the numbers of animals in the population by age and sex, their vulnerability by age and Table 1 Population dynamics equations. Table 2 provides definitions for the symbols.

Eqn No Equation Description
Population dynamics T1.1 Abundance by age and sex Breeding females Estimates of abundance Observed numbers captured by the gear t* denotes the set of years with estimates of total abundance and/or pup counts. t** denotes the set of years with estimated numbers dying and captured by the gear.
sex (ϕ s,f a , assumed to be 1 for all animals age 1 and older and zero for pups; pups are assumed not to be bycaught in the fisheries), the effort for fishery f during year t, and a "catchability" parameter (q f ) (Eqn T1.7). Some proportion, r, of surviving animals die following their capture in the fishing gear, and the value of r varies among fisheries (Eqn T1.9).

Model fitting
The parameters of the population dynamics model are: (a) MSYR 1+ , the bycatch rate at which maximum sustainable yield (MSY) is achieved given all age 1+ animals are equally vulnerable to capture (equivalent to the rate of increase when the population is at its MNPL); (b) the age-atmaturity; (c) the maximum pregnancy rate; (d) the survival rate of animals age 1+, S 1+ ; (e) the relative population size (in terms of 1+ individuals) at which maximum productivity is achieved, i.e., the ratio of the MNPL to carrying capacity (MNPL/K); (f) catchability by fishery, q f ; (g) the proportion of animals that are captured in the gear for fishery f that die immediately, p f ; (h) the post-release mortality for fishery f, r f ; (i) the relative vulnerability of animals by age and sex, ϕ s,f a ; (j) the number of age 1+ animals in a reference year relative to carrying capacity in 1990 (first year considered in the model); and (k) a measure of current abundance (either in terms of age 1+ animals or numbers of pups born). It is possible to compute pup survival S 0 and the two parameters of the density-dependence function (Eqn T1.4) given MSYR 1+ , the age-atmaturity, the maximum pregnancy rate, S 1+ , and MNPL/K (Punt, 1999). Parameter combinations for which the calculated S 0 exceeds S 1+ are rejected as biologically implausible.
The model-fitting process involves the application of a Bayesian assessment approach based on the sampling-importance-resampling (SIR) algorithm (Rubin, 1987;Van Dijk et al., 1987;Punt and Hilborn, 1997). The SIR algorithm is implemented separately for each zone (north, central, south) by drawing parameter vectors from priors and computing posteriors with increased numbers of draws from the priors until the posterior is represented by 100 unique parameter vectors (1000 for the base-case model). The results of a Bayesian assessment provide best estimates of, for example, current population size, but also a posterior distribution for current population size from which posterior percentiles can be computed (Punt and Hilborn, 1997).
The data used to fit the model are estimates of total (age 1+) abundance, pup counts (Supplementary Table 1), numbers observed dead in the gear, and numbers observed captured by the gear (Supplementary Table 3). The total abundance estimates are assumed to be absolute indices of abundance while the pup counts are assumed (for the base-case model and most of the sensitivity tests) to be indices of the number of pups immediately after birth and before the effects of density dependence come into play (Eqn T1.10 and T1.11). The estimates of total (1+) abundance and the pup counts are assumed to be log-normally distributed.
The likelihood function used in the SIR algorithm (Eqn T1.10-T1.13) contains contributions for the number of animals observed dead in the gear (Eqn T1.12) and the number captured by the gear alive and dead (Eqn T1.13), with both data sources assumed to be negative-binomially distributed. The expected number of bycaught animals observed dead in the gear is the product of the modelled number of animals that die immediately in the gear (Eqn T1.8) and the proportion of the total catch that is observed, p f (assumed to be known without error; Supplementary Table 3). The expected total number of animals captured in the gear that is observed is the product of the modelled number of animals captured in the gear (Eqn T1.7) and p f .
The likelihood includes parameters for the variance of the estimates of age 1+ abundance and the pup counts, the constant of proportionality for the pup counts (λ), as well as overdispersion parameters for the are "nuisance" parameters and hence are integrated out to provide a marginal likelihood (and to make the application of the SIR algorithm more efficient; Eqn T1.12 and T1.13). The variances of the estimates of age 1+ abundance and the pup counts are integrated out as part of the model-fitting process. Table 3 lists the base-case values for the pre-specified parameters of the population dynamics model and the likelihood function. There is no prior information for most of the parameters so uniform priors that are sufficiently broad to cover the likely range are adopted for several of the parameters. Effort by fishery is assumed to be constant for the base-case model.
The alternative models explored here (Table 4) relate to the prespecified aspects of the model. The factors in Table 4 examine sensitivity to: a) the extent of post-release mortality and whether it differs among fisheries (Alternative models 1-3); b) pre-specifying instead of estimating MSYR 1+ (the value of MSYR 1+ = 0.06 considered in Alternative model 4 is essentially half the maximum rate of increase assumed in the calculations by Wade, 1998); c) the value of MNPL/K (Alternative model 5); d) treating pup counts as measures of absolute (rather than relative) pup numbers (Alternative model 6); and e) the historical trend in effort (and hence bycatch mortality rate), assuming that the trends in effort are the same for all fisheries (Alternative models 7-8).

Removals in relation to policy goals
PBR is the product of three parameters: (1) a minimum estimate of abundance that "provides reasonable assurance that the stock size is equal to or greater than the estimate" 2 (N MIN ); (2) one-half of the maximum intrinsic rate of population growth (0.50 R MAX ); and (3) a recovery factor (F R ) between 0.1 and 1.0 (Wade, 1998): Mortality coefficient for fishery f (i.e., probability of dying due the fishery and being observed as dead) q f Catchability coefficient for fishery f (numbers captured by the gear) r f Proportion of the catch by fishery f that dies after being released.
x Plus-group age (values for the population dynamics parameters, including human-caused mortality rates, are the same from age x onwards)  Table 3 The pre-specified parameters of the base-case population dynamics model and the prior distributions for the estimated parameters for the base-case model.

Parameter
Value or prior Source Pre-specified parameters MNPL/K 0.6 Common assumption underlying marine mammal assessments (Taylor and DeMaster, 1993 Uniform over the allowable range ℓnλ

U[∞,-∞] Non-informative for a scale parameter
Calculations of PBR herein use: N MIN = the lower 20th percentile of the log-normal distribution of the most recent abundance estimate; R MAX = 0.12, the default value for pinnipeds, as applied in the USA; and F R = 1, which is appropriate for stocks at OSP (Wade and Angliss, 1997;Wade, 1998). In defining N MIN , the lower 20th percentile was chosen by simulation to meet the policy goals of: 1) remaining at or above MNPL (i. e., at the lower bound of OSP) for 20 years given the population started at MNPL, and 2) recovering to at least MNPL within 100 years, with a 0.95 probability, given the population started at 0.3 K (Wade, 1998).
The results of the model lead to three ways to compute PBR: (a) the product of 0.06 (half of the default value of R MAX for pinnipeds and F R = 1) and the lower 20th percentile of the posterior for 1+ abundance from the model; (b) the product of the posterior median for MSYR 1+ (which for an age-aggregated model with logistic growth is 0.5 R MAX ), the lower 20th percentile of the posterior of age 1+ abundance from the model, and F R = 1; and (c) the product of 0.06 and the lower 20th percentile of the sampling distribution of the most recent estimate of total abundance (where the CV is set as the median of the posterior for the sampling CV obtained by fitting the model). Option (c) is most consistent with the way PBR is typically computed because it uses only the most recent estimate of abundance, while options (a) and (b) represent using the model results as the primary basis for computing PBR. Fig. 2 shows the model-predicted time-trajectories of age 1+ abundance and pups (posterior medians and 50 % and 90 % intervals), along with the abundance estimates, the model-predicted observed number of bycaught animals, and the model-predicted number of animals observed dead in the gear for the base-case model. The model fits the abundance data well, although the pup counts for SASL in the north and south zones are highly variable (Fig. 2, 1st and 3rd rows; Supplementary Fig. 1). The model is able to replicate most of the data on observed bycaught animals and the proportion of those animals that die (data points often well within the 90 % probability intervals from the model), and the posterior medians of numbers dying and bycaught for each fishery almost always lie within the range of observations for that fishery. The poorest fits are when bycatch estimates are very low, in which case the model predicted higher values than the observations (e.g., numbers of SAFS and SASL observed dying in and captured by the gear in the trawl fishery in the south zone; Fig. 2, 3rd and 4th rows). The fits to total bycatch for the artisanal purse seine fishery for the SAFS (Fig. 2, 2nd row), and to total and dead bycatch for industrial trawl in the South for both species (Fig. 2, 4th row) are also poor, which is due in part to the high interannual variation in the estimates.

Fit diagnostics
The alternative models also fit the data adequately ( Supplementary  Figs S2− 9). Visually, the fit of the model is not noticeably poorer when productivity is defined by MSYR 1+ = 0.06 rather than being estimated or when the pup counts are assumed to be absolute rather than indices of abundance (Fig. 2 vs. Supplementary Figs S5 and S7), but this is because the estimated sampling variance is higher for these two alternative models.
Comparison of the alternative models with the base-case model using Bayes factors (Table 5) does not provide strong evidence that the alternative models are superior to the base-case model (values in Table 5 substantially larger than 1). The models with MSYL/K = 0.5 are best supported by the data, though the difference between even the best fitting of these models (SASL in the north zone) and the base-case model is still negligible ("not worth more than bare mention" according to Kass and Raftery [1988]). In contrast, alternative models 4 (MSYR 1+ = 0.06) and 6 (pup counts as absolute measures of abundance) consistently lead to Bayes factors less than ~0.3 among species / zones (implying "positive" evidence against these two alternative models; Kass and Raftery, 1988). Table 5 shows the posterior distributions for a subset of the parameters of the models for the four species/zones for the base-case model. The variance parameters in Eqn T1.11 and T1.12 are integrated out prior to the calculation of the likelihood and are therefore not shown in Fig. 3. The posteriors for the parameters determining the rates of post-release mortality by fishery, species and zone match the priors very closely (and therefore are not shown in Fig. 3), largely because post-release mortality is unobserved.

Fig. 3 and Supplementary
The prior for MSYR 1+ [0.02, 0.15] is updated substantially by the data for all species / zone pairs. The posteriors for MSYR 1+ for SASL by zone are centered below 0.06 (well below for the north and central zones) and that for SAFS in the south zone is centered well above 0.06 (half of the default value for the maximum rate of increase for pinnipeds when applying the PBR approach -Wade [1998]).
The posteriors for the sampling CVs for the abundance estimates are also updated substantially. The sampling CV for the estimates of age 1+ abundance is poorly defined for SASL in the central zone and SAFS in the south zone (wide posteriors encompassing most possible values). In contrast, the posteriors for the remaining sampling CVs exclude higher values.
The posterior for initial (1990) population size relative to carrying capacity is quite flat for SASL in the central zone, while it favours values larger than 0.5 for SASL in the north and south zones, and values lower than 0.1 for SAFS.

Population status and removals
The posterior for abundance in 2020 relative to carrying capacity for the base-case model suggests that SASL and SAFS in the south zone and SASL in the north zone were above MNPL with high probability at the start of 2020. In contrast, this posterior is quite uninformative for SASL in the central zone, with a probability that the number of age 1+ animals exceeds MNPL (half of carrying capacity) close to 0.5 (Fig. 3, right column, second row). Carrying capacity is poorly estimated for all species/zones, reflected in a wide distribution for the population size at equilibrium, in the absence of future removals (e.g., Fig. 2 first two columns).
All of the populations, except SASL in the central zone, are predicted by the base-case model to have increased somewhat between 1990 and 2020, with age 1+ abundance estimates ranging from 50 to 96% (posterior median) of carrying capacity in 2020. The 2020 relative population size varies, however, among the alternative models. The most optimistic results for SASL in the north, central and south zones are when MSYR 1+ = 0.06 (Model 4), and the best case for SAFS in the south zones is decreasing effort (Model 8) ( Fig. 4; Supplementary Table 6). The more optimistic results for SASL when MSYR 1+ = 0.06 are not No post-release mortality 4 MSYR 1+ is set to 0.06 5 MNPL/K = 0.5 6 The pup counts are assumed to be absolute indices of pup numbers 7 Effort doubles from 1990 to 2019 8 Effort drops by 50 % from 1990 to 2019 A.E. Punt et al. Fisheries Research 235 (2021) 105816 unexpected given that the posteriors for MSYR 1+ for the base-case analysis support values less than 0.06 for MSYR 1+ (Fig. 3). For SASL in the north and central zones, the least optimistic case is when MNPL is reduced from 0.6 to 0.5 (Model 5), while decreasing effort (Model 8) and MSYR 1+ = 0.06 (Model 4) are the least optimistic cases for SASL in the central zone and SAFS, respectively ( Fig. 4; Supplementary Table 6). Again, the least optimistic results for SAFS in the south zone when MSYR 1+ = 0.06 are due to the posterior in the base-case model supporting values for MSYR 1+ greater than 0.06. It is important to note, however, that the models with MSYR 1+ = 0.06 are not well supported by the data (Table 5). Fig. 5 shows the posteriors for the time-trajectories of bycatch mortality by species, zone and fishery for the base-case model. The distributions are generally very wide, reflecting that catchability (parameter q in Equation T1.7) is not well determined by the data. The   Fig. 2. Time-series of total population (age 1+ abundance) and age-0 (pup) numbers by species/zone, along with the data used to fit the model (left two columns). The dark lines in the left two columns are posterior medians, while the dark and light shading covers the 50 % and 90 % probability intervals, respectively. The projections beyond 2019 are based on the assumption of no future removals. The sampling intervals for the abundance indices account for estimated sampling variance. The two right columns show the fit to the bycatch data (data open symbols; closed symbols posterior medians and lines posterior 90 % intervals). The results are for the base-case model for which effort is constant over time.

Table 5
Bayes factors relative to the base-case model (marginal likelihood of the alternative models divided by that of the base-case model) by species / population. A Bayes factor larger than 1 implies a better fit than the base-case model and vice versa.  Table 3). Uncertainty is lower for the fisheries with higher levels of observer coverage (e.g., the factory trawl fishery in the south zone). Fig. 6 shows the posteriors for the time-trajectories of numbers captured by the gear by species, zone and fishery in the base-case model. As expected, the posterior time-trajectories of numbers captured and released (Fig. 6) qualitatively match those of numbers dying (Fig. 5), except for the purse seine fishery where most of the caught animals are released alive and the scale of the number captured can be substantially larger than the number dying.
The time-trajectories of bycatch and bycatch mortality for the alternative models are qualitatively identical to those for the base-case model, although estimates of removals in absolute terms depend on the assumptions of the model (Supplementary Table 6). There is no obvious pattern in which models lead to the highest or lowest levels of bycatch mortality between species and zones, owing to the interaction among time-trends in abundance, and values for catchability. Table 6 lists the posterior mean estimates of bycatch mortality during 2010-2019 for the base-case model and the eight alternative models. The posterior mean estimates of bycatch mortality vary by model (118-157 for SASL in the north zone; 493-561 for SASL in the central zone; 240-282 for SASL in the south zone; 95-121 for SAFS). The bycatch mortality for SASL in the north, SASL in the south and SAFS is well below the PBR value in each case, irrespective of how PBR is computed. In contrast, the bycatch mortality for SASL in the central zone is less than PBR when PBR is based on R MAX = 0.12 but above PBR when PBR is computed based on the estimated value of MSYR 1+ .

Discussion
The primary aim of this paper is to develop a method of marine mammal stock assessment that makes use of bycatch data in several fisheries when monitoring rates are uneven across fisheries. This is a common situation for marine mammal species, unlike many fish species for which there is usually a landed component of the catch that is well monitored. The method can be applied to situations in which some abundance data (ideally at least two or three estimates of absolute abundance) are available, along with some data from an observer program. Our approach was applied in a complex but common situation, where multiple fisheries using different types of fishing gear interact with multiple marine mammal species over a broad geographical range. This approach is intended to support the development and implementation of bycatch reduction programs.
Our method is similar to the common approach of dividing observed dead animals per unit of observed effort by the proportion of observer coverage. However, we integrate this calculation into the model-fitting process and hence capture the associated uncertainty more fully. Our method can also use information on the proportion of captured animals that are released alive but subsequently die from the effects of being bycaught. The Bayesian approach employed here allows for the use of priors in estimating the parameters of the model (although the only strongly informative prior for the Chilean case-study was that for recent abundance) and the assessment results can be summarized using probability distributions for model parameters, time-trajectories of population numbers (in principle by sex and age), and numbers dying by fishery. This presentation allows not only an appraisal of the precision with which quantities of management interest are estimated but also estimation of the probability that management and conservation reference points will be exceeded. In the context of the USA's Fish and Fish Product Import Provisions of the MMPA (50 CFR § 216.24), the relevant thresholds are the status of the population relative to PBR, and whether current removals meet the conservation objective (population recovery to, or maintenance at, OSP).
The methods presented in this paper are similar to those of Givens and Thompson (1995) and Zerbini et al. (2019) for cetaceans, with the exception that, unlike those approaches, which account for uncertainty in human-induced mortality by placing priors on the annual removals, this paper treats the number of observed marine mammals dead and alive in fishing gear as data to be fitted as part of the parameter estimation process. Unlike these alternatives, our method does not require estimates of removals for all years, but rather uses the model to make inferences about them, given an assumption regarding trends in effort, and fits the model to the available data on abundance. Cook (1997) developed a method for conducting assessments of fish stocks that does not require estimates of removals but relies on having information from surveys for population numbers by age for several yearsdata that are not available for most marine mammals. Porch et al. (2006) outlined a method of stock assessment that does not require catch data but instead relies on prior information on rates of change in abundance, and indices of relative abundance. The Porch et al. method could, in principle, be applied to species such as SASL and SAFS, but there is no obvious way to obtain prior information on the change in population size except by using the data included in the present analysis and expert opinion. In addition, unlike in Porch et al. (2006), estimates of absolute abundance are available for several stocks of the species considered here; such data are generally amongst the most informative when conducting population assessments.

Population status and trends
The estimates of population size relative to quantities such as K and MNPL are sensitive to assumptions regarding some of the features of the model. Moreover, population size relative to carrying capacity varies substantially among the species/zones considered here, with SASL in the north and south zones and SAFS in the south zone close to carrying capacity, and SASL in the central zone close to MNPL. However, it is possible to conclude from the evaluations here that abundance is increasing for all species/zones. Thus, pup production must exceed current rates of removal. The estimates of current depletion from the model are consistent with the classification of the SASL as a 'least concern' species according to the Chilean regulation for the classification of wild species in conservation categories 3 , and also according to the IUCN Red List (Cárdenas-Alayza et al., 2016a). SAFS in the south zone is classified as a 'least concern' species by the IUCN Red List (Cárdenas-Alayza et al., 2016b) but as 'near threatened' according to the Chilean regulation for the classification of wild species in conservation The USA National Oceanic and Atmospheric Administration (NOAA), which administers the portions of the MMPA pertaining to most marine mammals, evaluates status of stocks (populations) and issues regulations for reducing incidental take in commercial fisheries (bycatch). The Seafood Import Provisions, promulgated in 2016, require that imported fish and fish products be evaluated with respect to USA standards. Those regulations require countries that export fish and fish products to the USA, and that are identified by NOAA as having fisheries that involve or may involve marine mammal bycatch, to have a program for monitoring and mitigating bycatch (including prevention) that is "comparable in effectiveness" to that applied to USA fisheries. The implementation of the Seafood Import Provisions is likely to improve programs for estimating bycatch (e.g., increasing the observation effort). However, in many countries the time-series of bycatch estimates will remain short for several years after initiating improved monitoring programs, and therefore approaches to assessment that require timeseries of estimates of removals of ~10+ years in length cannot currently be applied.
It was possible to calculate PBR for SASL and SAFS using only the most recent estimate of abundance given a CV based on the model (Table 6). However, the assessment method of this paper allows a more integrative way to determine whether bycatch mortality exceeds PBR because the estimates of abundance used to calculate PBR are based on multiple years of data on abundance and bycatch as well as estimates of biological parameters.

Caveats and future work
The results presented here are sensitive to some of the assumptions of the assessment, in particular to the values assumed for MNPL/K (current relative stock size is lower of lower values for MNPL/K). This sensitivity occurs because there is no evidence for a reduction in growth rate from the (admittedly few) data on abundance so a population cannot be much larger than its MNPL. The results are predictably sensitive to MSYR 1+ , although the data provide some information on this parameter ( Table 5). The results depend on assumptions regarding post-release mortality, but the relationship between key model outputs (e.g., 2020 relative population size and estimated bycatch mortality) is not straightforward (e.g., Alternative Model 3 with no post-release mortality does not always lead to the lowest estimated bycatch mortality and highest 2020 depletion).
One key input for the assessment, to which the results are sensitive, is the proportion of captured animals that escape or are released but subsequently die (i.e., the post-release mortality rate). There is no information about this proportion for SASL and SAFS off Chile so sensitivity to the prior assigned to this parameter was explored (Table 4; Supplementary Figures S2-S4) 5 . We would expect the post-release (cryptic) mortality in the purse seine fisheries to be very low, because most of the animals escape from the fishing gear during net hauling, i.e. just a few are released alive on deck. In contrast, trawling may lead to higher post-release mortality at least for nets without escape panels. Post-release mortality could be assessed using a tagging program or observers could attempt to assess visually whether the animals released 'alive' would survive; such information could be used to develop an informative prior for this quantity. The USA has developed a policy for doing this in its fisheries (NMFS, 2012, and see Appendices in Carretta et al., 2018). Post-release mortality is not the only source of unobserved mortality, and in principle a prior could be placed on other such sources.
Another caveat is that the Chilean observer program is voluntary (some vessels do not have space for scientific observers and some companies do not allow observers on board) so reported bycatch rates likely do not reflect those in the rest of the fisheries, as is assumed in the present approach. This is of greatest concern for the fisheries that lead to highest mortality (likely trawls), although the observer coverage for the trawl fisheries is relatively high (29-100 %) when compared to most observer programs, and the trawl fishery, unlike the purse-seine fishery, involves a small number of vessels. This is also of considerable concern when coverage is very low because bias in bycatch estimations can be significant (Hall, 1999). In principle, the bycatch data included in the analyses could be restricted to those vessels considered to be operating "normally" (if such an appraisal can be made) because the method only requires a data set for which the observer coverage and bycatch observations are representative. A more serious concern is that the observed vessels are not necessarily spatially or temporally representative of the entire fishery. In addition, the method of this paper assumes that fishing practices were effectively the same between 1990-2014 and 2015-2017, i.e., the catchability parameter is assumed to be constant within a fishery over time. Whether this is the case is unknown, but the trawl fishery had higher effort in the past and consequently fishing practices may have changed.
There are no historical data on effort for the fisheries considered in this paper. It would have been straightforward to include data on effort for all or some fisheries, had such information been available. However, at least for SASL and SAFS off Chile, the posteriors for the number of animals at the start of 2020 relative to carrying capacity were not very sensitive to alternative assumptions regarding effort trends. This robustness was likely a consequence of the availability of abundance estimates that constrain the time-trajectory of population size.
The model was implemented within a Bayesian framework; thus, priors can be placed on the parameters of the model. Priors were placed on many parameters but some of the parameters (such as adult survival rate and age-at-maturity) were assumed to be known. In principle, metaanalyses could be used to develop priors for these parameters. However, previous applications of age-and sex-structured models such as that employed in this paper suggests that estimates of population size and population size relative to reference points will be robust to the prior for these quantities if MSYR 1+ is estimated and there are no data on population age composition or on the age composition of the bycatch.
We conducted assessments for four species-zone combinations; however, the parameters were estimated separately by species and zone. In principle (and perhaps at substantial computational cost), parameters such as the MSYR 1+ could be shared among species and zones or the MSYR 1+ for each species and zone could be assumed to be drawn from a hyper-prior. Either of these approaches would allow data from relatively data-rich species and zones to inform similar, data-poor situations.
The number of estimated parameters ranges between 12 (SASL in the north and south zone and SAFS in the south zone) and 18 (South American sea lions in the central zone). The parameters related to the variance of the numbers captured by the gear, and the constant of proportionality for the pup counts are "nuisance" parameters that are marginalized out, reducing the number of parameters that are included when applying the SIR algorithm. However, the number of parameters in the models is close to the limit of how many parameters can be estimated using SIR. If additional parameters were to be estimated, it could become necessary to use an alternative estimation framework such as Markov chain Monte Carlo (MCMC), perhaps implemented using Stan (Carpenter et al., 2017). Implementation of the model using MCMC algorithms would allow the application of other model selection methods such as the Deviance Information Criterion (Spieghalter et al., 2002) and Leave-one-out cross-validation information criterion (Vehtari et al., 2017) as well as other approaches for determining the adequacy of model fit. The model currently involves the solution of non-linear equations (e.g., the relationship between MNPL/K and MSYR 1+ , and the parameters of the density-dependence function) each time the likelihood is calculated, which may necessitate the development of a case-specific MCMC algorithm.

Final considerations
This paper has developed an approach that can be used for population model-based assessments of marine mammal stocks when some data on abundance are available but only limited information is available on bycatch. The approach we developed should be particularly helpful for fisheries that have some bycatch data including post-release mortality (and ideally information on the rate of post-release mortality), and could provide a way to assist countries needing to evaluate whether bycatch mortality is comparable with conservation policy goals. The method is applied to data for South American sea lions and South American fur seals off Chile and confirms that the populations are increasing, but that the 2020 population size relative to carrying capacity varies among zones and species.
The results of assessments such as those in this paper could be used as the basis for management strategy evaluations, such as determining whether a given alternative to the PBR approach for computing limits on bycatch levels will similarly ensure that marine mammal populations are maintained at or above, or if necessary allowed to increase towards, a desired benchmark level (e.g., MNPL).

Author credit statement
AP and JM developed the concept of the paper; SM and DO created the data sets; M Sepulveda developed the framework for the graphics; all authors contributed to writing and editing, and the experimental design of the simulations.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table 6
Posterior mean estimates of seals and sea lions dying due to bycatch mortality