The vaquita is a critically endangered species of porpoise. It produces echolocation clicks, making it a good candidate for passive acoustic monitoring. A systematic grid of sensors has been deployed for 3 months annually since 2011; results from 2016 are reported here. Statistical models (to compensate for non-uniform data loss) show an overall decline in the acoustic detection rate between 2015 and 2016 of 49% (95% credible interval 82% decline to 8% increase), and total decline between 2011 and 2016 of over 90%. Assuming the acoustic detection rate is proportional to population size, approximately 30 vaquita (95% credible interval 8–96) remained in November 2016.
1. Introduction
The vaquita (Phocoena sinus) is a critically endangered species of porpoise endemic to the northern Gulf of California, Mexico. Vaquitas have been subject to a long history of unsustainable bycatch in gillnets set by small-boat fishers targeting shrimp and finfish (Rojas-Bracho and Reeves, 2013). Recently, there has been a resurgence of an illegal gillnet fishery for an endangered fish, the totoaba (Totoaba macdonaldi), fueled by a lucrative illegal trade with China for totoaba swim bladders. This has raised concerns about increased levels of vaquita bycatch. Vaquita are difficult to monitor using standard visual survey methods (line transects or mark-recapture based on photo-identification) because they are small and visually cryptic, and now are very rare. However, they produce echolocation clicks almost continually, making passive acoustic monitoring of population trends possible. An acoustic monitoring program at a grid of locations throughout the core range of the vaquita was started in 2011; analysis of data to 2015 showed an estimated annual decline of 34% [95% credible interval (CRI) 21%–48%] (Jaramillo-Legorreta et al., 2017). Based on preliminary results through 2014, the government of Mexico enacted an emergency ban on gillnets, which began in 2015. However, despite extensive enforcement efforts, there is evidence that illegal fishing and bycatch continue (CIRVA, 2017). Here we report results from the 2016 acoustic monitoring.
2. Methods
Hardware deployment, acoustic processing, and trend analysis methods followed those described in detail by Jaramillo-Legorreta et al. (2017), and are outlined below.
2.1 Acoustic data collection and processing
Forty-six autonomous acoustic loggers were deployed between June and August each year from 2011 to 2016 in a systematic grid within the vaquita refuge, an area of core habitat designated in 2005 [Fig. 1(A)]. Because vaquita now number far fewer than when acoustic monitoring began, 47 additional loggers were added mid July 2016 within the same area to improve precision in future trend monitoring [Fig. 1(B)]. The loggers were C-PODs (Tregenza et al., 2016): autonomous passive acoustic monitoring instruments designed to detect echo-location clicks of toothed whales and store salient features for offline classification. Upon retrieval, proprietary software (C-POD software version 2.044 with KERNO algorithm) was used to detect coherent sequences (“trains”) of approximately regularly spaced clicks and classify them as possible vaquita click trains; all possible vaquita click trains were then manually validated. Data were summarized as the number of vaquita clicks detected per sampling location per day—the “click rate.”
(Color online) Summary of raw data for 2016. Mean click rate (clicks/day), indicated by shading (note log scale), and days of sampling, indicated by circle size, for (A) the 46 sampling sites grid during the core sampling period (June 19 to August 19), and (B) the augmented 93 sampling sites grid during the part of the core sampling period when additional sensors were deployed (July 11 to August 19).
(Color online) Summary of raw data for 2016. Mean click rate (clicks/day), indicated by shading (note log scale), and days of sampling, indicated by circle size, for (A) the 46 sampling sites grid during the core sampling period (June 19 to August 19), and (B) the augmented 93 sampling sites grid during the part of the core sampling period when additional sensors were deployed (July 11 to August 19).
2.2 Trend analysis
Only data from the original 46 sampling locations were used for trend analysis. Although this grid of locations was designed to give equal coverage in space and time, in practice sampling was uneven due to shifts in annual deployment dates, equipment failure, and loss. The data set was further truncated to a core sampling period during which at least 50% of the detectors were operating across all 6 years: June 19 to August 19, inclusive (62 days). Nevertheless, the sampling effort was still somewhat uneven, and so trend analysis was based on two Bayesian statistical models developed by Jaramillo-Legorreta et al. (2017): a geostatistical model and a non-spatial mixture model. They are described briefly here; full details are given in Jaramillo-Legorreta et al. (2017). Both models estimate average click rate per sampling location and year; the mean of these over locations is the average annual click rate. Between-year change in acoustic activity, denoted λ, is estimated as the ratio of the average annual click rates in successive pairs of years.
The geostatistical model compensates for locations with missing data by “borrowing strength” from those around it: the model assumes the average click rate varies smoothly over space, with a separate smooth surface fit to each year of data but with the amount of smoothness (the spatial autocorrelation) the same across years (i.e., autocorrelation parameters estimated from all years' data). It further accounts for variation in sampling by assuming locations with more sampling days give more precise estimates of the average click rate than those with fewer days.
The post-stratification mixture model probabilistically assigns individual sampling locations to one of three strata. A sampling location is permanently assigned to the same stratum for all years (which is justified based on spatial stability of the data), but each stratum mean click rate is estimated independently for each year. The purpose of stratification is to statistically account for much of the inter-site variance in the number of clicks recorded; the number of strata was chosen subjectively, after exploratory data analysis.
Uninformative prior distributions were used for all model parameters. For each model, inference was performed using Markov chain Monte Carlo methods: a single chain was run for 1 010 000 iterations using the WinBUGS (geostatistical model) and OpenBUGS (mixture model) software packages (Lunn et al., 2012), the first 10 000 samples were discarded and thereafter every 100th sample was retained for posterior distribution summaries, yielding 10 000 samples for each model. Results are reported separately for each model, but also combining samples from both to form a model-averaged estimate. Each model has equal weight in the model averaging.
2.3 Checking acoustic metric
For annual change in acoustic activity to represent change in abundance, there must be no systematic change in animal vocal behavior or range-specific click detection. Data were not available to measure these directly (e.g., through animal-borne recording tags); however, we searched for changes in the number of detected vaquita clicks in minutes when animals were known to be present (a measure of acoustic behavior), for trends in temperature (known to affect propagation) and a proxy for background noise. Details are given in the supplementary material.1
2.4 Abundance estimate for 2016
Vaquita population abundance was estimated from Bayesian analysis of a combined visual and acoustic survey, conducted in the fall of 2015, by Taylor et al. (2016). The posterior distribution for population size was well approximated by a lognormal distribution with mean 66 and standard deviation 33. To project the population forward from 2015 to 2016, 20 000 random samples were drawn from this lognormal distribution, and each sample was multiplied by a sample from the distribution of annual change in acoustic activity, from the trend models. Using November 2, 2015 (the midpoint of the visual survey) as the survey abundance date, the projected estimate represents population size on November 2, 2016.
3. Results
The sampling effort in 2016 was high for the 46 original sampling locations, with most C-PODs being operational for the entire core period [Fig. 1(A)]. As in earlier years, vaquita detections were restricted to only some portion of the refuge, with the highest click rates close to the southwest boundary; the additional 47 locations monitored for the first time in 2016 [Fig. 1(B)] make it clear that detections decline to almost zero on this boundary (and most other boundaries).
The recorded number of vaquita clicks per day in the 46 original locations decreased by 44% from 2015 to 2016 (). However, this statistic does not account for unequal sampling effort across sites. Results from the statistical models, which do account for an unequal effort, and give estimates of statistical uncertainty, are summarized in Fig. 2 (see also the supplementary material1). The values for 2011–2015 are very similar to those reported by Jaramillo-Legorreta et al. (2017), as would be expected since the addition of a sixth year of data should not change substantially results from the previous five. The estimated change in acoustic activity between 2015 and 2016 differs somewhat between the two models and has wide posterior CRIs; combining results from the two produces a model averaged posterior mean estimate of with 95% posterior CRI 0.18 to 1.08. Although this interval includes the “no change” value of 1.0, the posterior probability that acoustic activity decreased between 2015 and 2016 [i.e., ] is 0.96 (Fig. 3).
(Color online) Estimated annual rate of change (λ) and percentage change (%) in mean click rate from the statistical trend models. Solid symbols denote posterior means and vertical bars 95% posterior CRIs. “Geom. mean” is the geometric mean of the annual estimates.
(Color online) Estimated annual rate of change (λ) and percentage change (%) in mean click rate from the statistical trend models. Solid symbols denote posterior means and vertical bars 95% posterior CRIs. “Geom. mean” is the geometric mean of the annual estimates.
Model-averaged posterior probability distribution for the annual rate of change in mean clicks-per-day. The darker gray distribution describes geometric mean annual rate of change from 2011 to 2016. The lighter distribution describes the change between 2015 and 2016. Values less than 1.0 indicate a decline, for example, a value of 0.5 indicates a halving each year.
Model-averaged posterior probability distribution for the annual rate of change in mean clicks-per-day. The darker gray distribution describes geometric mean annual rate of change from 2011 to 2016. The lighter distribution describes the change between 2015 and 2016. Values less than 1.0 indicate a decline, for example, a value of 0.5 indicates a halving each year.
Another way to express the estimated change in acoustic activity is as a percentage rate of decline or increase (Fig. 2). The estimated model averaged percentage rate of decline between 2015 and 2016 was extremely high: posterior mean 49% decline (95% CRI 82% decline to 8% increase).
Over the entire monitoring period 2011–2016, the estimated average annual change in acoustic activity is 0.61 (95% CRI 0.48–0.74), i.e., an average decline of 39% per year (95% CRI 26% to 52%, Fig. 2). This corresponds to an estimated total decline of 90% over this 6-year period. The posterior probability of a decline is 1.00, and there is a >99% chance that the decline in acoustic activity has averaged >20% per year (Fig. 3).
The spatio-temporal pattern of estimated click rates from the geostatistical model is shown in Fig. 4. The spatial pattern of click rates is quite consistent between years, with the southwest area being a relative hotspot in all years. However, the strongest pattern is the overall steep decline in click rates over time.
(Color online) Estimated mean number of clicks per day predicted by the geostatistical model for the 46 numbered sampling sites with data for at least 1 year. Values in the legend are posterior medians (note log scale). Some sites, ⊗, were missing in the indicated year. The size of the circles indicate the number of sampling days on each year (see the legend).
(Color online) Estimated mean number of clicks per day predicted by the geostatistical model for the 46 numbered sampling sites with data for at least 1 year. Values in the legend are posterior medians (note log scale). Some sites, ⊗, were missing in the indicated year. The size of the circles indicate the number of sampling days on each year (see the legend).
No yearly change was found in the metric of acoustic behavior or background noise; water temperature generally increased over time, but the resultant effect on propagation loss is likely to be minimal (see the supplementary material1).
The projected population size estimate for November 2, 2016, using the model averaged acoustic results, is approximately 30 (posterior mean 33; posterior median 27; 95% CRI 8 to 96).
4. Discussion
Despite the gillnet ban, there is a very high probability that acoustic activity within the vaquita refuge has continued to decline between 2015 and 2016. As expected, the estimate for annual change in acoustic activity for a single year (2015–2016) is too imprecise to say whether it differs either from the previous year (2014–2015) or from the series of years (2011–2015). Given the low number of acoustic detections, it will take several years of monitoring to pick up such changes in trend, although the augmented monitoring design will help improve precision. There is little doubt, however, that the decline continues and is rapid.
Relating this decline in acoustic activity to a decline in population size requires the assumption that the acoustic behavior of vaquita and detection performance has not changed; our investigations (see the supplementary material1) provide partial validation of this. Another important assumption in extrapolating the trend observed in the vaquita refuge from acoustic data to total population trend is that population trend outside the refuge is the same as that inside. The combined visual-acoustic survey of Taylor et al. (2016) estimated that, in 2015, approximately 20% (12 of the estimated 59) of the vaquita population were outside the refuge. Hence, if the population change outside the refuge is different from that inside, it will have only a small effect on the overall trend, because it is only a small proportion of the total population. It is possible that animal distribution has changed radically between 2015 and 2016, with animals moving out of the refuge causing the observed decline. However, this seems very unlikely given low detections on the periphery of the range from the augmented grid [Fig. 1(B)] and the relative paucity of detections outside the refuge on previous synoptic surveys (e.g., Fig. 1 of Rojas-Bracho and Reeves, 2013).
If the acoustic change between 2015 and 2016 represents a population change, then we estimate that approximately 30 vaquita remained as of November 2, 2016. Only a fraction of these will be reproductive-age females. Given this estimate, and the ongoing negative population trend, it seems clear that we are facing the imminent extinction of this species unless radical and effective conservation measures are immediately implemented.
Acknowledgments
We thank the following bodies for support and funding: the Mexican Government (through the Mexican Secretaría de Medio Ambiente y Recursos Naturales), especially Minister R. Pacchiano and A. Michel; U.S. Marine Mammal Commission, in particular T. Ragen, R. Lent, and P. Thomas; the World Wildlife Fund (WWF) Mexico, in particular O. Vidal and E. Sanjurjo; Le Equipe Cousteau; The Ocean Foundation; Fonds de Dotation pour la Biodiversité; MAAF Assurances (Save Your Logo); WWF-US; Opel Project Earth; Fideicomiso Fondo para la Biodiversidad; Instituto Nacional de Ecología y Cambio Climático; Comisión Nacional de Áreas Naturales Protegidas; and Directorate of the Reserva de la Biósfera Alto Golfo de California y Delta del Río Colorado. Many thanks to our field staff J. XXIII Osuna, P. Valverde, and R. Arozamena, all the fishers who deployed and recovered the equipment, and to two anonymous reviewers.
See supplementary material at https://doi.org/10.1121/1.5011673E-JASMAN-142-514711 for an investigation of the acoustic metric used and results summary table.