Estimating IUCN Red List population reduction: JARA—A decision‐support tool applied to pelagic sharks

The International Union for Conservation of Nature's (IUCN) Red List is the global standard for quantifying extinction risk but assessing population reduction (criterion A) of wide‐ranging, long‐lived marine taxa remains difficult and controversial. We show how Bayesian state–space models (BSSM), coupled with expert knowledge at IUCN Red List workshops, can combine regional abundance data into indices of global population change. To illustrate our approach, we provide examples of the process to assess four circumglobal sharks with differing temporal and spatial data‐deficiency: Blue Shark (Prionace glauca), Shortfin Mako (Isurus oxyrinchus), Dusky Shark (Carcharhinus obscurus), and Great Hammerhead (Sphyrna mokarran). For each species, the BSSM provided global population change estimates over three generation lengths bounded by uncertainty levels in intuitive outputs, enabling informed decisions on the status of each species. Integrating similar analyses into future workshops would help conservation practitioners ensure robust, consistent, and transparent Red List assessments for other long‐lived, wide‐ranging species.


INTRODUCTION
Identifying species at risk is essential to prioritize conservation efforts against accelerating biodiversity loss (Butchart et al., 2010;Hoffmann et al., 2008). The International Union for Conservation of Nature's (IUCN) Red List of Threatened Species is the global standard for quantifying extinction risk (Hoffmann et al., 2008;Mace et al., 2008). In the ocean, overexploitation is the overwhelming driver of biodiversity loss (McClenachan, Cooper, Carpenter, & Dulvy, 2012), so a population reduction estimate (Red List Criterion A) is usually applicable when assessing marine vertebrates (Dulvy et al., 2014). However, few marine timeseries are long enough to directly calculate the percent reduction over three generations lengths (3GL) required under Criterion A (IUCN Standards and Petitions Subcommittee, 2019). Fewer still are broad enough in spatial coverage to easily undertake a global Red List Assessment, especially for wide-ranging species (Boyd, 2010;Dulvy, Sadovy, & Reynolds, 2003;Godfrey & Godley, 2008).
Bayesian state-space models (BSSM) offer a powerful and flexible framework to model variable population trends (e.g., Boyd, DeMaster, Waples, Ward, & Taylor, 2017;Meyer & Millar, 1999;Winker, Carvalho, & Kapur, 2018). BSSMs also have three key properties that could help improve the objectivity of IUCN Red List assessments. First, the posterior probabilities provide an intuitive way to express uncertainty around rates of population change to conservation practitioners (Sherley et al., 2018). Second, missing values can be estimated automatically, providing a robust and transparent way of dealing with timeseries of differing lengths and quality, and of forecasting future population trajectories (Kindsvater et al., 2018). Third, it is simple to combine the posterior probabilities for regional population trends (and their uncertainty) into a global reduction rate. BSSMs have been used to assess extinction risk under the U.S. Endangered Species Act (Boyd et al., 2017), but rarely to assign IUCN categories (Regehr et al., 2016;Rueda-Cediel, Anderson, Regan, & Regan, 2018).
Here, we outline how a BSSM tool, Just Another Red List Assessment (JARA, https://github.com/henningwinker/JARA; Winker & Sherley, 2019;Sherley et al., in press), facilitated the assessment of 13 wide-ranging pelagic and coastal-pelagic sharks at an IUCN Species Survival Commission Shark Specialist Group (IUCN SSC SSG) Red Listing workshop (Dallas, USA, 5-9 November, 2018). We illustrate our approach using four circumglobal shark species: Blue Shark (Prionace glauca), Shortfin Mako (Isurus oxyrinchus), Dusky Shark (Carcharhinus obscurus), and Great Hammerhead (Sphyrna mokarran). These species differed in data quality and availability, generation length, and previously published Red List category (Dulvy et al., 2014). Finally, we discuss the wider applicability of tools like JARA for the Red List assessment of long-lived, widely distributed taxa.

Workshop application
Regional relative abundance datasets for each species were analyzed using JARA, a generalized BSSM tool for global extinction risk estimates under IUCN Red List Criterion A (Winker & Sherley, 2019;Sherley et al., in press). The input timeseries were either formal stock assessment outputs (trends in biomass), or standardized or nominal catch per unit effort (CPUE) from scientific surveys, fisheries data or bather protection nets, depending on the data available for each species and region (Table 1, see Supporting Information for details). Initial results (e.g., Figure 1) were presented to workshop participants (see Acknowledgements) based on the available timeseries ( ), generation lengths (GL ; Table 1), and the proportional geographic area that each region comprised of a species' global range (Figures 2  and 3). Participants either approved the data choices or suggested alternative datasets with better temporal or spatial coverage for additional model runs. The final choices were made by consensus. Participants were then presented with easyto-interpret outputs for each species showing (a) fits to the observed regional data, (b) observed annual rates of change, (c) any projections (where necessary), (d) how the posterior distribution for the percentage change over 3GL aligned against the thresholds for Red List criteria A2 (see Supporting Information for definition; e.g., Figure 1), and (e) the most likely IUCN Red List category (Table 1): Critically Endangered (CR), Endangered (EN), Vulnerable (VU) (the threatened categories), Near Threatened (NT) or Least Concern (LC) (IUCN, 2012). Here we highlight the decisionsupport nature of JARA and accordingly outline the choices made and decisions taken by consensus of the workshop participants. In particular, because other forms of information, such as geographic range, habitat and ecology, threats, use and trade, and conservation actions, are considered in assessments, in addition to the modeled population trend, the final Red List category proposed may may differ from the one suggested by JARA. Moreover, after the workshop results (which we report here) were finalized, it was occasionally necessary to conduct further JARA runs following comments from reviewers and consultation with 166 IUCN SSC SSG members (https://www.iucnssg.org/who-we-are.html). Thus, the published IUCN Red List assessments (e.g., www.iucnredlist.org/species/39341/2903170) may ultimately differ from those presented here.

State-space model formulation
Each timeseries ( ) was assumed to follow an exponential growth model: +1 = + where is the logarithm of the expected abundance in year t, and the normally distributed annual rate of change with mean̄, the estimable T A B L E 1 Population change (%) and posterior probabilities (%) for changes falling within the IUCN Red List categories Least Concern (LC), Near Threatened (NT), Vulnerable (VU), Endangered (EN), and Critically Endangered (CR) for Blue Shark (Prionace glauca), Shortfin Mako (Isurus oxyrinchus), Dusky Shark (Carcharhinus obscurus), and Great Hammerhead (Sphyrna mokarran). The most likely status based on criteria A2 is indicated according to the category that contained the highest proportion of the posterior probability for the rate of change over three generation lengths (3GL), with the exception that VU is also indicated as a precaution in cases where LC obtained the highest probability but with <50% of the total posterior probability. The global change is based on weighting the regional posterior probabilities by the proportional area (PA) weighting, an area-based proxy of the percent of the global population in each region based on each species current geographic range (see Supporting Information) Notes: 1 Stock assessment output (ICCAT, 2016); 2 stock assessment output (Carvalho & Winker, 2015); 3 stock assessment output (ISC Shark Working Group, 2017); 4 stock assessment output (Takeuchi, Tremblay-Boyer, Pilling, & Hampton, 2016); 5 stock assessment output (Rice, 2017); 6 stock assessment output (ICCAT, 2017); 7 stock assessment output (ISC Shark Working Group, 2018); 8 standardized catch per unit effort (Francis, Clarke, Griggs, & Hoyle, 2014); 9 stock assessment output (Brunel et al., 2018); 10 stock assessment output (SEDAR, 2016); 11 standardized catch per unit effort (Braccini & O'Malley, 2018); 12 catch per unit effort (Dudley & Simpfendorfer, 2006); 13 stock assessment output (Jiao, Cortes, Andrews, & Guo, 2011); 14 catch per unit effort (Carlson J.K. & Driggers W.B. unpubl. Data); 15 catch per unit effort (Dudley & Simpfendorfer, 2006). N. = north; S. = south; NW. = northwest; NE. = northeast; E. = eastern; W. = western. Species and regionally specific GL were calculated from female age at maturity ( mat ) and maximum age ( max ) as = (( max − mat ) × 0.5) + mat (see Supporting Information). For Great Hammerhead, Global 1 used North Atlantic 1 data to generate the weighted change and Global 2 used North Atlantic 2 data. mean rate of change for a population, and process variance 2 (see Supporting Information for details). We linked the logarithm of the observed relative abundance , for index expected abundance trend, using the observation equation and observation variance assumptions and priors presented in the Supporting Information. We used a noninformative normal prior for̄∼ Normal(0, 1000). Priors for the process error variance were 2 ∼ 1∕gamma(0.001, 0.001), or approximately uniform in log space (Winker et al., 2018).

Regional change
The percentage change in abundance in each regional index was calculated from the posteriors of the estimated population timeserieŝ= exp( ). If the span of was longer than 3GL, the percentage change was automatically calculated as the difference between a 3-year average around the final observed data point , and a 3-year average around year − (3GL) (e.g., Figure S1). The year + 1 is always projected to obtain a 3-year average around (to reduce the influence of short-term fluctuations; Froese, Demirel, Coro, Kleisner, & Winker, 2017). When was < 3GL, JARA projected forward, by passing the number of desired future years to the BSSM until̂> (3GL) + 2 (e.g., Figure 1C). Projections were based on the posterior of medians across all T observed years ( Figure 1B). The projection gives similar results to extrapolating backward to attain a 3GL period, as recommended when estimating population reduction from annual rates of change (IUCN Standards and Petitions Subcommittee, 2019). As these are not the "moving window" reductions required for criterion A4 (IUCN Standards and Petitions Subcommittee, 2019), we assessed results against criterion A2. JARA visualizes the extent and uncertainty in these projections, to aid final category choice.  Pac.) and Indian Ocean. The South Pacific Ocean timeseries is observed data from 1994 to 2014 and projected from 2015 onward (see Figure  S4). All others are fits to observed timeseries from stock assessment outputs (see Table 1

Weighted global change
To assess each species globally, we subsampled from the posterior probability distribution for each regional percentage change according to the proportional area (PA) weighting of that species' geographic range occurring in each ocean region (Table 1, Figures 2 and 3). Parts of a species range with no trend data were treated as missing regions. To meet the IUCN guidelines for uncertainty associated with missing regions, we assumed that these regional populations had declined by between 0 and 100% (IUCN Standards and Petitions Subcommittee, 2019) by subsampling from a uniform distribution Uniform(−100, 0) according to the PA weighting for the missing regions. All of these sub-samples were combined to produce the weighted global rate of change (Figures 2 and 3).

Model implementation
We ran JARA (v. 1.33) in JAGS (v. 4.3.0) via the "R2jags" library (v. 0.5-7) for R (v. 3.5.0), using three Monte Carlo Markov chains of 220,000 iterations, burn-in of 20,000 and thinning to every second observation. Each chain was initiated by assuming a prior distribution on the initial state centered around the first data point of each abundance timeseries ( = 1 ), 1 ∼ Norma(log( 1 ), 1000). Convergence was diagnosed using the "coda" R library (v. 0.19-1) and minimal thresholds for Geweke's diagnostic thresholds of p = 0.05 (Geweke, 1992). All models unambiguously converged.

Timeseries availability and use of projections
For Blue Shark (Table 1, Figure 2) all timeseries used were stock assessment outputs and were > 3GL except in the South Pacific where projections were ∼33% (∼1GL) of the timeseries ( Figure S1). For Shortfin Mako, stock assessment outputs were available for four regions, with standardized catchper-unit-effort (CPUE) available for the South Pacific; all timeseries required some projection to span 3GL (Table 1). However, South Atlantic biomass estimates were reported to be highly uncertain and implausible (ICCAT, 2017), so the workshop participants decided to use the North Atlantic stock assessment to represent the Atlantic Ocean. For the Shortfin Mako, Dusky Shark, and Great Hammerhead all timeseries were < 3GL (Table 1), so projections were required to reach 3GL ( Figures S7-S18). Data were missing for ∼50% of the ranges of Dusky Shark and Great Hammerhead and input data were a mix of CPUE and stock assessments (Table 1), which the workshop participants needed to take into consideration when finalizing their assessments.   The "Missing" polygons (in c and e) show the samples drawn from a uniform distribution, Uniform(−100, 0) according to the summed proportional area (PA) of the species' geographic range contained within those regions for which no regional trend data were available (see Table 1). The Global polygon shows the weighted global rate of change, with each regional posterior weighted according to the species-specific PA calculation (right-hand panels and Table 1). For Great Hammerhead (e), only the outputs using the more recent North Atlantic dataset (North Atlantic 2, Table 1)

The globally distributed, unmanaged, data-rich Blue Shark
In the Pacific and Indian Oceans, the medians and >80% of the percentage change posterior probabilities fell into LC (increases or reduction <20%; Figure 2, Table 1). In the South and North Atlantic, the medians and posterior probabilities suggested VU and EN respectively ( Figure 2, Table 1), but the populations were stable over at least the last GL (Figure 2, S1 and S2). Consequently, the combined regional posteriors produced a bimodal distribution for the global percentage change ( Figure 2); the median was −7.3%, and ∼66% of the posterior distribution was consistent with LC (Table 1). Despite the slight population reduction, catch and trade of this species is largely unregulated and the workshop participants cited a suspected global population reduction of ∼20-30% over the last 3GL to suggest retaining a NT assessment.

The globally distributed, partially managed, Shortfin Mako
In the Atlantic, which has the longest observed timeseries, the median and 90% of the posterior probability indicated a status of EN (Figures 1 and 3A). The last 2GL also had a faster annual rate of decline than the whole timeseries ( Figure 1B). In both the North Pacific and Indian Oceans, the median and the largest proportion of the posterior distributions (∼44%) suggested VU (Table 1, Figure 3A), although ∼44% of the posterior also indicated EN in the Indian Ocean (Table 1). For the South Pacific, uncertainty was very high ( Figure 3A; trend predominately based on projections), but a slightly positive annual rate of change (+0.48, Figure  S9) meant ∼70% of the posterior distribution suggested LC ( Figure 3A). The combined regional posteriors produced a global median of −46.6%, but 43.3% of the posterior distribution indicated a 50-80% reduction over 3GL, meeting the threshold for EN (Table 1). These findings, combined with steep historic declines in the Mediterranean Sea, global underreporting of catches, limited control of fishing mortality, highly valued meat, and international trade in moderate-value fins led workshop participants to recommend a global status of EN.

Steep but uncertain decline in the Dusky Shark
Both Indian Ocean CPUE datasets yielded broad posterior distributions spanning LC to CR ( Figure 3C), but annual rates of change (−2.8% and −0.9%, Figures S13 and S14) meant that >70% of the posterior overlapped the threatened categories (VU-CR; Table 1). Similarly, the 56-year stock assessment output for the North Atlantic produced an annual change of −2.6%, which had worsened in the most recent GL to −3.9% ( Figure S12). This generated a median change of −89.9% over 3GL, with 99.1% of the posterior in CR (Figure 3, Table 1). When weighted, the regional posteriors suggested a global listing of CR, though the median (−71.6%) did not exceed that categories' threshold (≥80%). The workshop participants inferred a global population reduction of 50-80% over 3GL and recommended an EN assessment for Dusky Shark, based on a balance of the estimated declines and the possibility for stabilization and slow recovery in parts of its range following management action. This was confirmed after the analysis of an additional dataset during review of the final Red List assessment.

Steep declines worldwide with signs of Atlantic recovery in the Great Hammerhead
An initial model using stock assessment output for the North Atlantic showed a change of −29% over 3GL (North Atlantic 1, Table 1, Figure S16). However, some workshop participants noted that important U.S. management interventions had been implemented since that dataset ended in 2005 (NMFS, 2006). This led to an additional analysis replacing the North Atlantic stock assessment output with one nominal and one standardized CPUE timeseries underpinning that stock assessment (updated to 2017). This model run indicated population recovery at 6.6% per annum ( Figure S17); despite high uncertainty, 100% of the change posterior fell into LC (North Atlantic 2, Table 1, Figure 3C). By contrast, the Indian Ocean (the only other dataset) showed an annual change of −6.6% ( Figure S18) and a 100% probability of exceeding the CR threshold (Table 1, Figure 3C). Using the newer North Atlantic dataset (North Atlantic 2), and accounting for the three regions with missing data using uniform distributions (see section 2.4), the posterior of the weighted global rate of change spanned a threat status from LC to CR, with ∼36% at each extreme (Table 1). While recognizing that management can be successful, even for species with low biological productivity, the workshop participants considered the North Atlantic to be the exception globally. Given the trajectory in the Indian Ocean, scarcity, quality and spatial coverage of the datasets relative to the regions they represented, and ongoing fishing pressure and trade, the workshop participants inferred a global reduction of >80% over the last 3GL and listed Great Hammerhead as CR.

DISCUSSION
We present a straightforward solution to the persistent challenge of rigorously undertaking Red List assessments for marine species with sparse data and widespread distributions. We use a classic statistical approach to estimate population growth rate from timeseries (Dennis, Munholland, & Scott, 1991), while incorporating the power and flexibility of BSSM. Although updated frequently, the IUCN guidelines (IUCN Standards and Petitions Subcommittee, 2019) were laid out long before the widespread availability of Bayesian techniques, and provide varying levels of guidance on a range of key issues (e.g., model specification, combining multiple indices, coping with uncertainty). Below, we discuss our approach to these issues and outline the opportunities presented by using JARA to guide Red List decision-making. Below, we highlight five useful features of our approach: (1) avoids the need for assessors to choose between different models of population change (e.g., linear vs. complex), (2) reduces the influence of outliers in individual timeseries, (3) outputs are intuitive and easy to communicate, (4) calculation of global status from the weighted sum of regional statuses, accounting for regions without data using uniform distributions, (5) the application of the tool to policy decisions. The IUCN guidelines provide extensive specification on how to estimate population reduction based on four different models of decline (linear, exponential, accelerating, and complex;IUCN Standards and Petitions Subcommittee, 2019). With JARA, assessors need never agonize over model choice again; they are provided with posterior distributions of growth rates decomposed into GL subsets to diagnose whether population trajectories are accelerating or complex ( Figure 1B). For example, in the North Atlantic, Blue Shark population change has stabilized over the last 2GL, while the population reduction of Shortfin Mako has accelerated slightly (cf. Figures 2B and S1B).
Moreover, timeseries are noisy with outliers (due to process and observation error) and the classical regression approaches often used to undertake Red List assessments in the past (Fox et al., 2019, IUCN Standards andPetitions Subcommittee, 2019) can result in biologically implausible rates of change. Consequently, the debate on status often vacillated on the inclusion (or not) of single data points (NKD, HW, pers. obs.), necessitating sensitivity and counterfactual analyses. A more objective solution is needed and although JARA is somewhat constrained by the values at the start and end of a timeseries (e.g., Figure S16), it is far less sensitive to outliers and more accurately captures the rate of population change than traditional regression approaches (Fox et al., 2019;Winker & Sherley, 2019). At the workshop, we found JARA especially useful in preventing the need to subjectively exclude putative outliers from individual timeseries.
During an assessment, the ramifications of the differing degrees of process and observation error in the data type used (stock assessment vs. CPUE) and the use (or not) of projections must be clearly articulated (Carvalho, Lee, Piner, Kapur, & Clarke, 2018;d'Eon-Eggertson, Dulvy, & Peterman, 2015). Here, the visualizations provided by JARA (e.g., Figure 1) allowed participants to rapidly examine the various input timeseries and provide more appropriate regional data if needed. This allowed the workshop to pause, re-run analyses with new data, then reopen discussions and evaluate the sensitivity of conclusions to differing inputs (e.g., as for Great Hammerhead; cf. Figures S19 and S20). The JARA outputs also helped to contextualize the uncertainty within each regional and global assessment (Figures 2  and 3). Participants with little or no experience with Bayesian statistics were quickly able to interpret posterior distributions of percentage change, based on the median, category with the greatest probability, and how the balance of probability (the posterior distribution) was divided between LC and the threatened categories. For example, the workshop participants could easily discriminate the reliability of the modeled trends for the North Atlantic Dusky Shark (56 years of stock assessment output) versus the Eastern Indian Ocean (10 years of CPUE data, cf. Figures S12 and S13). They took this uncertainty on board and spent considerable time debating a recommendation, with discussions on life-history characteristics, geographic range, habitat preferences, catchability, degree of overlap with fisheries, and positive impacts of management evident in slow recovery used to reach consensus in a process similar to previous elasmobranch assessments (Dulvy et al., 2014).
A key challenge of undertaking Red List assessments of widely distributed marine taxa has been objectively combining different regional or subpopulation assessments into a global reduction estimate (Godfrey & Godley, 2008;Reynolds, Dulvy, Goodwin, & Hutchings, 2005). The guidelines recommend calculating a global average, weighted by the number of mature individuals in each subpopulation 3GL ago (IUCN Standards and Petitions Subcommittee, 2019). Detailed abundance data are occasionally available for marine mammals and seabirds. However, for fishes, it is more common to have population reconstructions from stock assessments or indices of abundance, like CPUE (e.g., Carvalho et al., 2018). A key feature of JARA is the ability to store and combine posterior distributions by proportional weighting to generate a representative global posterior (see Supplementary Information). After much consideration, workshop participants choose to use IUCN distribution maps to calculate proportional areas as an indirect approximation of the number of mature individuals 3GL ago, but this approach could be further refined using species distribution models (Ready et al., 2010).
The IUCN Criterion A population reduction-the rate of decline scaled by generation length-is a robust method of assessing extinction risk of exploited marine fishes (e.g., d' Eon-Eggertson et al., 2015), which tends to align well with stock assessments statuses (e.g., Dulvy, Jennings, Goodwin, Grant, & Reynolds, 2005;Davies & Baum, 2012). The proposed listing of Shortfin Mako Shark on Appendix II of the Convention on the International Trade in Endangered Species (CITES) was not supported by the United Nations Food and Agriculture Expert Panel, who "found no evidence that populations meet the CITES criteria, whether based on historical extent of decline or recent rates of decline" (FAO, 2019). Here, we have shown, using all available data (while incorporating uncertainty), that the recent rate of decline (over the last 3GL) is steep, particularly in the North Atlantic (Figures 1  and 3A). Consequently, this species is globally Endangered with "a very high risk of extinction" due to insufficiently regulated fisheries and absence of international trade regulation (www.iucnredlist.org/species/39341/2903170). These results contributed to the IUCN/TRAFFIC analyses of the proposal and were presented to the Parties intersessionally (IUCN and TRAFFIC, 2019). As such our findings have immediate application to policy decisions, for example at future CITES Convention of Parties.
The need for robust decision-support tools to facilitate management and evidence-based policy making in the face of uncertainty has never been greater (Polasky, Carpenter, Folke, & Keeler, 2011). Misclassification on the IUCN Red List can have real consequences for species conservation; coupling BSSM with expert oversight has great capacity to ensure more robust, consistent, and transparent assessments. Although still being refined, JARA is widely applicable and using it to support expert decision making would benefit the assessment of other long-lived, wide-ranging marine taxa, including sea turtles, seabirds, fishes, and marine mammals (Boyd, 2010;Godfrey & Godley, 2008).

ACKNOWLEDGMENTS
We thank Caroline Pollock (IUCN Red List Unit), members of the IUCN SSC SSG who contributed to and reviewed species assessments and the workshop participants; Rodrigo Barreto, Daniel Fernando, Sonja Fordham, Malcom Francis, Rima Jabado, Kwang-Ming Liu, Andrea Marshall, and Evgeny Romanov. William Driggers kindly provided unpublished Great Hammerhead data from NOAA Fisheries Service's Southeast Bottom Longline Survey. The scientific results and conclusions, as well as any views or opinions expressed herein, are those of the authors and do not necessarily reflect those of their institutions. Funding was provided by Natural Science and Engineering Research Council Discovery and Accelerator Awards, Canada Research Chairs Program, US National Science Foundation (DEB-1556779), and the Shark Conservation Fund. This is a Global Shark Trends Project output.