Beyond Post-release Mortality: Inferences on Recovery Periods and Natural Mortality From Electronic Tagging Data for Discarded Lamnid Sharks

Accurately characterizing the biology of a pelagic shark species is critical when assessing its status and resilience to fishing pressure. Natural mortality (M) is well known to be a key parameter determining productivity and resilience, but also one for which estimates are most uncertain. While M can be inferred from life history, validated direct estimates are extremely rare for sharks. Porbeagle (Lamna nasus) and shortfin mako (Isurus oxyrinchus) are presently overfished in the North Atlantic, but there are no directed fisheries and successful live release of bycatch is believed to have increased. Understanding M, post-release mortality (PRM), and variables that affect mortality are necessary for management and effective bycatch mitigation. From 177 deployments of archival satellite tags, we inferred mortality events, characterized physiological recovery periods following release, and applied survival mixture models to assess M and PRM. We also evaluated covariate effects on the duration of any recovery period and PRM to inform mitigation. Although large sample sizes involving extended monitoring periods (>90 days) would be optimal to directly estimate M from survival data, it was possible to constrain estimates and infer probable values for both species. Furthermore, the consistency of M estimates with values derived from longevity information suggests that age determination is relatively accurate for these species. Regarding bycatch mitigation, our analyses suggest that juvenile porbeagle are more susceptible to harm during capture and handling, that keeping lamnid sharks in the water during release is optimal, and that circle hooks are associated with longer recovery periods for shortfin mako.


INTRODUCTION
Quantifying fishing-related (F) and natural (M) mortality continues to be one of the main challenges in understanding and managing marine fauna. Representative starting values and priors for M are needed for demographic analyses (e.g., Cortés, 2016), evaluating resilience to population decline (e.g., Gedamke et al., 2007;Au et al., 2015), estimating extinction risk (e.g., García et al., 2008), and stock assessment (e.g., Cortés, 1998Cortés, , 2002. For elasmobranchs in particular, M is typically approximated from life history information, using previously derived functional relationships with longevity, growth or size (Kenchington, 2014;Cortés, 2016;Pardo et al., 2016). To varying extents, common methods rely on age determination, and are calculated from theoretical longevity, length-at-age and weight-at-age relationships, and/or von Bertalanffy growth function parameters (reviewed in Kenchington, 2014). This means all methods are sensitive to the level of uncertainty in age determination for elasmobranchs, where longevity specifically may be systematically underestimated Harry, 2018;Natanson et al., 2018). Underestimation of maximum age results in an overestimation of M from life-history based methods. There is a pressing need to move away from lifehistory based estimates of M to more direct estimates derived from species-specific data. Electronic tagging is an important source of information on movement, habitat associations and survival of large pelagic fishes (Hammerschlag and Sulikowski, 2011;Hazen et al., 2012), and provides an opportunity to directly estimate natural mortality from survival data (e.g., Benoît et al., 2015Benoît et al., , 2020a. Nonetheless the substantial cost associated with archival tags still constrains sample sizes (Hazen et al., 2012) and poses a particular challenge for reliable estimation of M for long-lived species.
Pelagic sharks tend to have high interaction rates with highseas fisheries targeting swordfish and tunas, and the majority of global shark catches represent bycatch (Lewison et al., 2004;Oliver et al., 2015). The magnitude of shark bycatch and the need for mitigation to reduce population declines (Dulvy et al., 2014) have driven recent research on shark survivorship following release (Ellis et al., 2017;Miller et al., 2020). In the North Atlantic, shortfin mako (Isurus oxyrinchus) and porbeagle (Lamna nasus) are two species for which landings have decreased in recent years and discard rates are increasing as a result of national and international management measures. A large proportion of discards have the potential to be released alive, given that estimated at-vessel mortality rates range from 35-56% for shortfin mako and 21-44% for porbeagle (reviewed in Ellis et al., 2017). Although quantifying rates of post-release mortality (PRM) remains a priority for future stock assessments to improve estimates of total removals, additional consideration of variables that affect survivorship is critical to develop effective bycatch mitigation measures (Davis, 2002;Ellis et al., 2017).
Capture and handling are two separate processes that can influence survivorship of bycatch (Benoît et al., 2012). For the majority of species, different handling protocols in addition to tagging effects are very rarely evaluated because they are assumed to be negligible in relation to capture effects (Musyl et al., 2009;Molina and Cooke, 2012;Jepsen et al., 2015). For sharks, research on survivorship tends to consider only covariates with at-vessel and/or post-release mortality. In general, lamnid sharks appear to be quite resilient to various types of capture and handling . However, sublethal effects on behavior and/or physiology are likely even though individuals survive (Skomal, 2007). Several studies report changes in swimming and dive behavior upon release, indicative of a recovery period (e.g., Skomal and Chase, 2002;Sippel et al., 2011;Wilson et al., 2014;Whitney et al., 2016). Any behavioral changes associated with recovery from physiological stress may ultimately contribute to mortality by making animals more susceptible to disease or predation, less able to forage, and/or more susceptible to recapture (Davis, 2002;Jepsen et al., 2015). Thus, mitigation measures designed to reduce the duration of any recovery period following release or to minimize capture and handling effects could be relevant when developing best practices to reduce shark bycatch mortality.
For this study, we compiled data from satellite tagging on porbeagle and shortfin mako sharks in the North Atlantic. Deployments were conducted by Canada, the United States, Portugal, and by the International Commission for the Conservation of Atlantic Tunas (ICCAT) through the Shark Research and Data Collection Program. Of the shark species whose status is regularly assessed at ICCAT, shortfin mako and porbeagle are currently considered overfished with a very high probability (ANON, 2019(ANON, , 2020. Recovery planning for both species would benefit from improved mortality estimates for stock assessment, as well as from the development of bestpractices for mitigation of bycatch mortality. For these purposes, our objectives were to infer M from survivorship data in light of relatively small sample size, to characterize any recovery period following tagging from changes in dive depths and periodicity, and to evaluate covariate effects on PRM and/or the duration of any recovery period.

MATERIALS AND METHODS
This study combined data from 177 archival satellite tag deployments during 2001-2019 in the North Atlantic (Figure 1), 73 on porbeagle and 104 on shortfin mako (Supplement 1). Both species were captured during regular commercial fishing activities by pelagic longline fleets (N = 134), scientific cruises using pelagic longline (N = 38) or commercial trawl trips (N = 5) and tagged by fisheries observers, science personnel, or fishermen trained by science personnel. Tags were attached to the sharks by tethering a dart anchor into the dorsal musculature, immediately beside the posterior end of the first dorsal fin (Campana et al., 2016;Musyl et al., 2011). Anchors consisted of either nylon umbrella darts (Domeier Anchor) or titanium darts, excluding the single deployment with an experimental fin clamp. Stainless steel wire or 400 lb test monofilament line (∼15 cm) was used to tether the tags to the anchor and the wire/line was sheathed in high temperature heat-shrink tubing to prevent chaffing at the point of attachment and to protect the leader. The PSAT tags were programmed to release from the sharks with the anchor and wire assembly remaining attached to them. Computers), and X-tags (N = 11; Microwave Telemetry). All tags recorded depth, either directly or through pressure, which was used to evaluate behavior and survival following tagging (PSATLIFE tags 0.05% resolution for pressure; survivalPAT, miniPAT, PAT 4, and PAT Mk 10: ± 0.5 m depth; and X-Tags: 0.34 m depth resolution). Each tag type reported the archived depth data at a different temporal resolution, ranging from a single daily maximum and minimum from survivalPAT tags to values at 5-min intervals from the PSATLIFE tags. Deployments were a maximum of 28 days for PSATLIFE tags and 30 days for survivalPAT tags. The PAT tag was deployed for 19 days with the experimental fin clamp. Longer-term deployments were possible from the other tag types and maximum deployment durations were 255 days for miniPATs, 204 days for X-tags, and 356 days for PAT MK10s (Supplement 1).

Post-release Behavior and Inferring Mortalities
Behavioral changes following tagging were assessed from recorded depth (pressure) profiles. SurvivalPAT tags provided no information on daily dive variability and were not included in the behavioral analyses. Daily dive variance (σ 2 ) was calculated from dive amplitudes and initially used to characterize behavior following tagging. For example, dive depth was calculated as the maximum minus the minimum depth recorded for each summary interval (4, 6, 8, or 12 h summaries) for miniPAT and PAT Mk 10 tags, and then variance was calculated from these depths for each day. No attempt was made to impute missing values for days without transmitted data. Porbeagle have been shown to exhibit limited vertical movement (i.e., low variability in dive depths) and residency at the surface following the physiological stresses associated with capture and release, indicative of a recovery period (Hoolihan et al., 2011). In our data, low variability in dive depths upon release was always associated with residency in the top 60 m of the water column. Thus, we identified the animals that exhibited a recovery period following capture and handling as those with low variability in dive depth coupled with residency in the top of the water column at the start of the deployment (Supplement 2: Supplementary Figures 1C, 2B). To quantify the duration of the recovery period, we identified the day on which dive variability markedly increased. Variance increased substantially once an animal started to dive more regularly and more deeply (i.e., maximum dive depths and periodicity increased). We identified the end of the recovery period as the day with the maximum difference between dive variance at the start of the deployment vs the remainder of the deployment. This involved sequentially calculating the difference in variance among time periods throughout the track, i.e., comparing day 1 vs day 2 onward, days 1-2 vs 3 onward, days 1-3 vs 4 onward, and so on (Supplement 2: Supplementary Figure 2C). Compared to analyses that use eigenfunctions and orthogonal axes to determine irregular post-release behavior (e.g., Hoolihan et al., 2011), using variance was computationally simpler and had a direct ecological interpretation in terms of how behavior was changing over time. Inconsistent sampling frequencies among the tag types and programmed settings prevented analyzing dive behavior using more sophisticated statistical methods such as Wavelet analyses (e.g., Thorburn et al., 2019) or the fast Fourier transform (e.g., Shepard et al., 2006). It is important to note that some animals exhibited similarly restricted diving behavior at other times during monitoring, which may have been related to geographical position. However, if restricted diving behavior upon release was solely a function of geographical position, it would not be expected to be functionally related to tagging covariates. Although we report the estimated duration of recovery periods for each individual (Supplement 1), our analyses of recovery time is focused on comparisons of mean recovery time between two groups.
Mortality events were inferred from continual records at a constant depth for multiple days (indicative of a dead animal on the bottom) or pop-ups following progressively increasing depth records up to the tag crush depth (indicative of an animal that is sinking; e.g., Musyl et al., 2011). Thus, the tag data tracked survival in continuous time (days until death) with right-censored observations from individuals that lived until the end of the observation period. The observations were censored because the ultimate time of death of the individual is unknown, yet the animal was known to be alive until the end of the observation period (Cox and Oakes, 1984). To separate postrelease mortality events (i.e., mortality associated with capture and handling) from natural mortality events (i.e., independent from the capture process), we evaluated patterns in dive behavior for animals that ultimately died. Similar to the evaluation of dive tracks from individuals that lived, we identified animals that were negatively affected by capture and handling as those with near-zero variability in dive depth coupled with residency in the top of the water column upon release. A mortality event that followed such a period of restricted diving behavior, with minimal evidence of re-establishment of cyclical movement, was considered post-release mortality and directly related to capture and handling. There was a single instance where an individual abruptly died yet had not exhibited any prior behavior that could be attributed to capture and handling. This mortality event was sudden and preceded by dive depths and periodicity consistent with those observed from animals that lived until the end of the observation period (comparison in Supplement 2: Supplementary Figure 1). This mortality was suspected to represent a natural mortality event.

Factors Influencing Recovery and Survival
There were several characteristics of the capture and handling process that could be evaluated from these tag deployments. The covariates that were considered included fork length, stage (juvenile, adult), sex (male, female), gear type (longline, trawl), hook type (circle, J), hooking location (mouth, gut), and handling location (in water, on-board) (Supplement 1). Note that the gut category for hooking location included gut-hooked (5 shortfin mako, 10 porbeagle) and foul-hooked individuals (0 shortfin mako, 6 porbeagle; Supplement 1). When categorizing life stage, we used sex-specific length at 50% maturity to separate juveniles from adults, with values of 182 cm and 280 cm fork length (FL) for male and female shortfin mako  and 174 cm and 218 cm FL for male and female porbeagle . There were only 58 deployments on shortfin mako and 57 on porbeagle that had information for the entire suite of covariates (Supplement 1).
The properties of the data on recovery times (e.g., sparse, zeroinflated) made typical parametric regression analyses unsuitable, so we used a randomization test to evaluate relationships with covariates. The main assumption underlying this approach is that the observed sample is representative of the larger population. We ran 10,000 samples to characterize the distribution for the mean difference in recovery times between factor levels of each covariate, implemented in the "simpleboot" package in R (Peng, 2019). The distribution of differences would be centered on zero if there was no effect of the covariate on recovery time and the proportion of samples with means that fell below zero represented the p-value for the comparison. To evaluate any association between recovery time and the continuous covariate FL, we used a Spearman Rank Correlation test. Relationships with hooking location and gear type could not be examined because there were insufficient data in one of the categories.
The influence of covariates on survivorship for each species was assessed using Cox proportional hazards models (CPHM; Cox, 1972;Therneau and Grambsch, 2000). CPHM are a well-established semi-parametric approach that estimates the multiplicative effect of covariates on a common hazard function, which describes the time-specific instantaneous probability of dying at a given time t, conditional on having survived to t. For each CPHM, the proportional hazards assumption was tested based on trends in the Schoenfeld residuals and was assessed visually by plotting the log of the negative log survivor function vs the log of event time. To provide the best inferences possible in light of missing covariate data, we undertook two series of analyses of the influence of covariates using CPHM. In the first, each covariate was modeled individually using all available observations, with no attempt to impute missing values. In the second series, we limited the data to observations for which values were available for all covariates (N = 58 for shortfin mako and N = 57 for porbeagle). This second series of analyses was intended to identify the suite of covariates associated with survivorship. A forward-selection scheme based on Akaike's Information Criterion (AIC) was employed, For shortfin mako, a model with hook-related injury resulted in a decrease in AIC of 10.5 compared to an intercept only model, but no other single or multiple covariate models were found to be comparable or superior based on AIC. For porbeagle, there were no models including a single covariate that resulted in a reduction in AIC of at least two units compared to an intercept-only model. Therefore, models incorporating multiple covariates were not pursued further and we report the results for individual covariates only, using all available observations. As in the behavioral analyses, statistical significance was accepted at p < 0.05.

Estimating Post-release and Natural Mortality
A CPHM does not distinguish between components of mortality, so we used the parametric mixture model of Benoît et al. (2015) to estimate separate rates of catch-related post-release mortality (PRM) and natural mortality (M). Specifically, the survivorship to time t, S(t), was modeled as: where α and γ are parameters of a Weibull survival function that describes the attrition of fish that will die after release due to the capture and release event, π is the post-release mortality rate, and M is the instantaneous annual rate of natural mortality (for a derivation see Benoît et al., 2015). Model parameters were estimated using maximum likelihood (details in Benoît et al., 2015Benoît et al., , 2020b. The non-parametric Kaplan-Meier (KM) estimator (Cox and Oakes, 1984) was used to visualize the survivorship of the two species, providing a basis for visually assessing model fit.
The model in Equation 1 effectively parses out mortality into capture-related (PRM) and natural (M) components based on their assumed time course. Post-release mortality is considered to asymptote over a finite timespan, typically within hours or days (reviewed in Musyl and Gilman, 2019). Meanwhile, released individuals are continuously at risk of dying from natural causes such as disease or predation, and an exponential function is commonly assumed in population modeling. The model can freely and reliably estimate the two mortality components provided that sufficient observations are available for both early rapid mortality and the later time periods (Benoît et al., 2015). Alternatively, the estimation can be aided by specifying the cause of mortality for some or all observations (Benoît et al., 2020a). In this study, patterns in dive depths and periodicity suggested that 33 of the mortalities of shortfin mako were catch-related, and only one had the potential to be natural. All mortalities of porbeagle appeared to be catch-related. Therefore, we fit the parametric model above with three variations: (1) fixing M = 0, which attributes all observed mortality events to PRM, (2) estimating M using the full model above, and (3) using the full model with cause-specific classifications of mortality (shortfin mako only). The cause-specific estimation for shortfin mako was accomplished by specifying different likelihood equations for the different classes of event observations (Benoît et al., 2020a;Kneebone et al., 2020). Specifically, observations for which the cause of death was inferred to be catch-related employed a likelihood in which M was fixed at 0, those for which the cause was assumed to be natural employed a likelihood in which π was fixed at zero, and those of uncertain cause employed the full likelihood for Equation 1. Similarly all censored observations employed the full likelihood as these individuals were at risk of dying from both catch-related and natural causes.

Simulation Modeling to Further Infer Natural Mortality Rates
Life history-based estimates of M for pelagic sharks are very low relative to other fish species (Cortés, 2002), suggesting natural mortality events are rare. The probability of observing natural deaths during the course of a tagging experiment should be correspondingly low, particularly when the median deployment duration from all tag types was 28 days for both porbeagle and shortfin mako. This likely explained why a natural mortality event was only observed once in these data. We used a simulation model to allow for inferences on the probable magnitude of M for each species given the observations made in this study. Our approach determined the probability of observing no natural deaths during the experiments for porbeagle and the probability of one or fewer for shortfin mako, as a function of the natural mortality rate.
Following the method by Bender et al. (2005), each iteration of the simulation proceeded as follows. Vectors of mortality probabilities, Z(t), with lengths corresponding to the total number of mortality event observations for each species were generated by randomly selecting values from a uniform distribution over the interval [0,1]. Assuming exponential natural mortality, mortality event times from each individual, t M,i , in days, associated with each value of Z(t) i for a given simulated annual natural mortality rate M s were calculated as: A censoring time, t C,i , was simulated for each individual by sampling with replacement from among the mortality event times for each species in the tagging experiments. Instances in which t M,i ≤ t C,i reflect a simulated instance in which an individual died from natural causes while or before dying from catch-related causes or having its tag detach. The proportion of iterations for which no individuals (porbeagle) or one or no individuals (shortfin mako) died from natural causes for a given value of M s is the estimated probability that the observed number of natural deaths occurred at that rate of natural mortality. It becomes less likely to observe no natural mortality events over the duration of the study as the magnitude of M increases. We also simulated the probability of observing no natural deaths for shortfin mako to illustrate the extent to which a single observation can change the probabilities associated with different natural mortality rates. Ten thousand iterations were undertaken for each M s value, which ranged from 0.02 to 0.70, with increments of 0.02.

RESULTS
The opportunistic tagging resulted in a range of sizes of both species and sexes, with slight oversampling of shortfin mako < 100 cm FL and fewer than expected porbeagle between 150 and 170 cm FL (Figure 2) relative to typical length-frequency distributions from landings data (Coelho et al., 2018;Santos et al., 2020). The vast majority of tagging occurred on juvenile animals, consistent with the selectivity patterns in longline fisheries (ANON, 2019(ANON, , 2020. Animals ranged in size from 78 to 249 cm FL (mean = 163 cm FL) for porbeagle and 66-240 cm FL (mean = 144 cm FL) for shortfin mako. The sex ratio of tagged animals was skewed in both species, with more females for porbeagle and more males for shortfin mako. Sample sizes varied substantially across the different tagging covariates (Table 1), as was expected from opportunistic tag deployments. Given the tag types used, we could determine the recovery period following tagging for 59 shortfin mako and 53 porbeagle. We distinguished pre-and post-recovery periods using a sharp change in variance (Supplement 2: Supplementary Figure 2C),  and a comparison of mean daily dive variance during and following the recovery period showcases the substantial increase in dive depths and periodicity following recovery. During the inferred recovery period, the median variance was near-zero for both species, while it increased to ∼5,000 following the recovery period (Figure 3). For individuals that remained in the top of the water column following tagging, there were no instances where dive variance was greater immediately following tagging than in the remainder of the deployment. The majority of sharks that died after tagging did so relatively quickly, many within hours. All mortalities of porbeagle occurred within 45 days of release and there were many long-term survivors, some with monitoring for up to a year ( Figure 4A and Supplement 1). Similarly, all mortalities of shortfin mako shark occurred within 50 days of release, and there were many long-term survivors, including some with monitoring times in excess of 200 days ( Figure 4B and Supplement 1).

Factors Influencing Recovery and Survival
The estimated durations of recovery for porbeagle (mean = 9.1 days) were similar to previous evaluations of recovery periods based on dive behavior for multiple species of pelagic teleosts (mean = 7.1 days) and pelagic sharks (mean = 10.8 days) (Hoolihan et al., 2011;Musyl et al., 2015). The estimated durations of recovery for shortfin mako (mean = 3.8 days) tended to be lower. For the animals that survived, there were no differences in mean recovery time between the sexes of either species, between different hooking injury types for both species, or between juvenile and adult shortfin mako (Table 1). There was a negative relationship between recovery time and fork length for both species, which was significant only for shortfin mako ( Table 1). For porbeagle, mean recovery times were 4.35 days longer (95% CI: 1.22-7.60) for juveniles and 3.63 days longer (95% CI: 0.87-6.77) when tagged onboard a vessel (Table 1 and Supplement 2: Supplementary Figures 3, 4). For shortfin mako, recovery duration was 1.97 days longer (95% CI: 0.66-3.39) when captured on circle hooks as compared to J hooks, and 3.31 days longer (95% CI: 1.61-5.22) when tagged onboard a vessel as compared to in the water ( Table 1 and Supplement 2: Supplementary  Figures 4, 5). The effect of hook type (circle, J) on porbeagle, or hook injury (foul, gut, and mouth) on shortfin mako or porbeagle could not be assessed due to extremely low or zero sample sizes in one of the categories ( Table 1). Survivorship of shortfin mako was significantly lower (at the 5% level) for individuals hooked in the gut rather than the mouth (Table 1), with a hazard ratio of 3.47 (95% CI: 1.04-11.61) (Supplement 2: Supplementary Figure 6A). Survivorship was also significantly lower for individuals tagged onboard rather than in-water, with a hazard ratio of 2.96 (95% CI: 1.01-8.67) (Supplement 2: Supplementary Figure 6B). Survivorship of porbeagle was significantly affected only by the manner in which fish were hooked on the fishing gear (Table 1). Compared to individuals hooked in the mouth, survivorship was reduced for individuals hooked in the gut or more generally foul-hooked. Combining the latter two categories to increase sample size, the hazard ratio for fish hooked elsewhere than in the mouth was 8.49 (95% CI: 2.21-32.46), constituting an important reduction in survival (Supplement 2: Supplementary  Figure 7). For both species, risk of mortality was negatively associated with increased fork length, though the effect was not statistically significant.

Estimating Post-release and Natural Mortality
Visual evaluation suggested fits from the parametric survival model were comparable to those from the non-parametric Kaplan-Meier (KM) estimator (c.f. Figures 4A-D). For porbeagle shark, the parametric model with only post-release mortality (i.e., M = 0) fit the trends in survivorship very well, producing an estimate of PRM of 0.171 (95% confidence interval: 0.099-0.277) (Table 2 and Figure 4D). An identical estimate of PRM was obtained when M was estimated in the model because the estimate of the M parameter was essentially zero, with an exceedingly wide confidence interval (results not shown). It is therefore not possible to directly estimate natural mortality for porbeagle using data from these tagging experiments. For shortfin mako, the parametric model with only post-release mortality fit the trends in survivorship very well, producing an estimate of PRM of 0.358 (95% confidence interval: 0.259-0.479) (Table 2 and Figure 4C). As with porbeagle, the model that attempted to freely estimate natural mortality produced an identical PRM estimate and an estimate of the M parameter that was essentially zero, with an exceedingly wide confidence interval (results not shown). In contrast, the cause-specific estimation (i.e., when one natural death event was identified in the data, see Figure 4E) produced an estimate of post-release mortality of 0.339 (0.246-0.453) and an estimate of natural mortality of 0.101 (0.016-0.659) ( Table 2). This model also provided a good fit to the survivorship trend, although the uncertainty around survivorship at later times was greater and increasing in time compared to the model excluding natural mortality. This pattern reflected the uncertainty associated with the additional and ongoing natural mortality component (Figure 4E).

Simulation Modeling
The absence of observations of natural deaths for porbeagle during the tagging experiments is consistent with the species' having low natural mortality. The simulation model suggests probabilities of ≤0.10 associated with each natural mortality rate above 0.15 (Figure 5). In contrast, the observation of a single natural death for shortfin mako resulted in substantially higher probabilities for the same natural mortality rates. The probability of observing one or no natural deaths only dropped below 0.10 when M was greater than 0.3. We also ran the simulations for a scenario assuming no natural deaths had been observed for shortfin mako, to evaluate how assigning cause to mortality events affects predicted PRM rates as well as how overall monitoring duration affects the simulations. In that scenario, the  panels (A,B) the circles represent right-censoring times, where the size of the circle indicates the number of censored observations, and the filled squares and crosses, respectively, indicate an inferred natural death event and death events with uncertain cause. All other mortality events were inferred to be related to capture and handling. In panels (C-E) the dotted line is the KM estimate, plotted as a reference. simulated probability for shortfin mako was essentially double that of porbeagle at 0.2 when M = 0.15 (Figure 5). There were 27 porbeagle that were monitored for > 90 days (∼3 months) as compared to only 15 shortfin mako (Supplement 1). Extended monitoring periods using archival tags increases the chances of observing mortality from natural causes. If such mortality events are not observed despite longer monitoring, there is greater certainty that the rate of natural mortality is low, as was the case for porbeagle.

DISCUSSION
Large sample sizes involving extended monitoring periods would be optimal to directly estimate M from satellite tagging data for lamnid sharks. Yet it remains possible to infer probable values or to constrain estimates, even in the absence of direct observations. For porbeagle, there was less than a 10% probability associated with values of M higher than 0.15 based on the simulation modeling. From the survival mixture model and a  single suspected natural mortality event, the maximum likelihood estimate of M was 0.101 for shortfin mako. To put these rates in perspective, approximately 1.5% of a population is expected to live to maximum age (Hewitt and Hoenig, 2005). Under a simple exponential model for mortality, 1.5% of the population would live to be ∼28 years (M = 0.150) for porbeagle and ∼41 years (M = 0.101) for shortfin mako. This would be on the lower end of longevity estimates for porbeagle in the Northwest Atlantic (24-43 years; Natanson et al., 2002), as might be expected from an upper limit of M. Our estimate of longevity for shortfin mako in the North Atlantic also falls within the expected range of longevity of 20-52 years (Natanson et al., 2006;Rosa et al., 2017). In all, there was fair correspondence with rates derived from life history for both species, even though our data came primarily from juvenile animals. Estimates of M based on age and growth parameters, maturity, and longevity most often yield a single value, and variability is generated by applying different types of estimators (e.g., Cortés, 2002Cortés, , 2016 or by allowing for variability in longevity when using a single estimator (e.g., Bowlby and Gibson, 2020). Using the Then et al. (2015) suite of estimators based on longevity and growth data, M ranged from 0.081 to 0.267 for porbeagle and from 0.068 to 0.318 for shortfin mako for males and females combined. Our estimates of M from survival data fell within these ranges, which lends credence to the natural mortality values currently being used in stock assessment (ANON, 2019(ANON, , 2020 and gives independent support that our current understanding of these species' biology is largely representative. Our results highlighted the types of information on survivorship that can be gained from long-term vs shortterm tag deployments. Estimating post-release mortality and evaluating the influence of covariates with survivorship was the original goal for the majority of the tagging contributing to this study. This explains the predominance of satellite tag types optimized for ≤30 days; chosen to reduce cost and increase sample size, given that shorter monitoring periods are generally sufficient for estimating post-release mortality in large pelagic fish (e.g., Musyl and Gilman, 2019;Benoît et al., 2020b). Although our PRM rates were similar to those of previous studies on these species, there were still limited and unbalanced data relative to covariates, which reduced statistical power and thus detectability of relationships (Sippel et al., 2015). Mortality related to capture and handling apparently extended beyond a 30-day monitoring period, which suggests PRM rates for porbeagle and shortfin mako could have been slightly underestimated if derived exclusively from short-term deployments (e.g., Marçalo et al., 2018). Finally, short-term deployments reduced the potential for natural mortality events to be observed over the duration of the experiment, making them suboptimal for species characterized by low M. Our results support several of the discussion points from a recent FIGURE 5 | Simulated probability of observing n or fewer mortality events resulting from natural causes during the tagging experiments, as a function of the annual natural mortality rate for porbeagle (n = 0, black solid line) and shortfin mako (n = 0, dashed gray line; n = 1, solid gray line). The inset histogram summarizes the observed mortality event times for the tagged porbeagle (black) and shortfin mako (gray). meta-analysis of post-release mortality in pelagic sharks , demonstrating that a substantial number of tag deployments is required to tease apart fishing-related mortality from M, that long-term deployments are necessary to increase precision of M estimates, and that a minimum 3-month pop-up period (>90 days) would be useful when trying to separate post-release from natural mortality as opposed to relying on short-term archival tags.
There is some debate on whether delayed mortality for pelagic sharks can be linked to capture and handling or whether PRM would be expected to asymptote relatively quickly. Several survivorship studies in addition to ours have reported delayed mortality, up to 50 days following release (summarized in Musyl and Gilman, 2019). Although mortalities that occur within hours of the tagging event are readily ascribed as PRM (Sulikowski et al., 2020), it is less clear if longer-term mortalities should be attributed to capture and handling (Hutchinson et al., 2015). Sublethal effects, such as reduced activity levels upon release (e.g., Raoult et al., 2019), reflex impairment and physiological damage (e.g., Jerome et al., 2018), or measureable changes in distribution (e.g., Bullock et al., 2015) could ultimately result in delayed mortality due to increased susceptibility to disease and predation, or cessation of feeding (Davis, 2002;Campana et al., 2016). Also, M is continuous and can occur at any time, irrespective of the length of time since the capture event. Instead of categorizing mortality as PRM and M based on a subjective timeframe, we used dive behavior to indicate whether mortality was likely related to capture and handling. We felt this was appropriate given the definitive contrast in dive variability that characterized recovery, and the correspondence between our estimates of the duration of recovery and previous evaluations of dive behavior from archival satellite tagging data (e.g., Campana et al., 2009;Hoolihan et al., 2011). Interestingly, the only likely natural death was recorded from a mouth-hooked male shortfin mako (157 cm FL) within 17 days of tagging. This individual's dive track exhibited variability equivalent to recovered individuals until the mortality event. Given that this was one observation, it would be beneficial to explore the utility and robustness of behavior-based classifications of mortality in future research.
Evaluating capture and handling covariates relative to recovery time as well as survivorship allowed for a more fulsome use of the tagging data and strengthened the inferences that could be made. As in other PRM studies, the majority of tagged animals survived, giving relatively few observations of mortality events from which to infer the effect of covariates (Sippel et al., 2015). Incorporating behavioral analyses of dive patterns from surviving individuals was an inexpensive and straightforward way to increase the amount of information gained, with some of the differences in recovery period being significant even when differences in survivorship were not. For example, the estimated coefficients from the CPHM suggested a non-significant increase in survivorship of porbeagle with fork length. There was a corresponding significant increase in mean recovery time following release for juveniles as compared to adults (i.e., juveniles took longer to recover from capture and handling) as well as a significant decrease in recovery time with fork length for shortfin mako. Taken together, we conclude that capture and handling was more detrimental to smaller juveniles of both species, although it is important to recognize that decreased swimming performance caused by carrying the tag would also be affecting these smaller animals (Todd Jones et al., 2013). Although enhanced international cooperation and additional tagging would be optimal to bolster sample sizes (Ellis et al., 2017;Harcourt et al., 2019), we suggest that quantifying recovery periods from surviving individuals is an additional avenue to explore the effect of covariates with capture and handling. Ideally, better standardization among tag types and programmed settings would support more complex statistical analyses of behavior (e.g., Shepard et al., 2006;Thorburn et al., 2019), and would be useful to give more precise estimates of the duration of recovery periods. In the absence of this, our comparisons provide meaningful information on sublethal effects that arise from specific characteristics of the capture and handling process for porbeagle and shortfin mako.
Revealing covariates with injury and mortality is important for developing mitigation options for non-retained bycatch (Molina and Cooke, 2012;Ellis et al., 2017). One of the most consistent relationships in our study was related to handling, specifically onboard vs in-water tagging. The > 3 day difference in recovery time for both species in addition to the significant reduction in shortfin mako survivorship suggests that physiological stresses associated with removal from the water significantly outweigh any benefit of gear removal following capture. Although trailing gear is commonly thought to contribute to PRM (Gilman et al., 2016), all animals that were tagged in the water for this study were released by cutting the gangion, thus retaining the hook plus an unquantified amount of monofilament leader (no weights or steel leaders). When tagged onboard, the shark remained under duress for a longer period and may have been subject to physiological damage when lifted out of the water and/or from the animal's inability to support its own weight while onboard (Musyl et al., 2009). Studies that directly evaluate handling effects (in isolation from capture effects) are rare, but longer handling times and increased exposure to air have also shown a significant negative effect on activity levels upon release for Squaliform and Carcharhiniform species (Raoult et al., 2019). Scientific work benefits from deep and precise insertion of the tag anchor to reduce the probability of pre-mature tag loss, which is easier to accomplish when the animal is onboard (Biais et al., 2017). In our study, in-water tagging of porbeagle used PSATLIFE tags only (28 day maximum deployment) and it was difficult to determine if pre-mature pop-ups were related to anchor placement. This batch of tags had a 40% non-transmission rate (Bowlby et al., 2019) indicating other tag construction and/or software issues. Shortfin mako that were tagged in the water had a longer median monitoring duration as compared to those tagged onboard (c.f. 60 days vs 57.5 days). Although it was not possible to determine the specific characteristics of boarding that increased recovery time for both species and decreased shortfin mako survival (e.g., the method of lifting the animal out of the water, the duration the animal was onboard, the method of gear removal), in-water release from commercial captures and in-water tagging for scientific work appears optimal for lamnid sharks.
In terms of best-practices for the release of bycatch from commercial interactions, our results support the recommendation to release sharks immediately upon capture, leaving embedded hooks and as little trailing line as possible . Contrary to earlier suggestions that handling practices have little influence on the condition of sharks upon release (Campana et al., 2009;Musyl and Gilman, 2019), handling in and of itself was associated with substantial sublethal effects. Our results also support the general recommendation to increase protection of the juvenile life stages of bycaught species (Ellis et al., 2017), optimally by minimizing the potential for interaction through spatial management. However, they are less clear relative to optimal hook type. On one hand, increased gut hooking is expected from capture on J hooks (Epperly et al., 2012;Gilman et al., 2016), where gut and foul hooking were associated with significantly higher post-release mortality for both species. However, shortfin mako exhibited longer recovery times following release when caught on circle as opposed to J hooks, possibly because circle hooks are harder to remove (Cooke and Suski, 2004) or may not be expelled from the jaw as quickly (Poisson et al., 2019). Such apparently contradictory results underscore the multi-facetted nature of bycatch mitigation, where it is often unclear if benefits relative to one component of the capture process are outweighed by detriments to another (Reinhardt et al., 2017). Ultimately, taking a holistic approach to bycatch mitigation is necessary, particularly to make any trade-offs explicit in the overall management approach .

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The animal study was reviewed and approved by Fisheries and Oceans regional Animal Care Committee.

AUTHOR CONTRIBUTIONS
HDB conceived research, undertook data collection and fieldwork, developed analyses, wrote manuscript, and contributed in funding. HPB developed models, undertook analyses, and reviewed manuscript. WJ, CS, and BA undertook data collection and fieldwork and reviewed manuscript. JS, RC, AD, EC, FH, DM, and GB undertook data collection and fieldwork, reviewed manuscript, and contributed in funding. All authors contributed to the article and approved the submitted version.