Each BRUVS was considered an individual deployment and analyzed as a single site, although the BRUVS on the line were probably not independent from each other. Our survey effort corresponded to an east to west (72° E−134.5° W) and north to south (7.6° N − 33° S) transect, straddling the Coral Triangle ( Fig 1A ), and spanned a range of human pressures [ 56 ]. All sampling was conducted between the 17th of April 2012 and the 25th of January 2015, during daylight hours, between 09:30 and 17:00 local time.

Marine predators were surveyed using standardized midwater stereo-BRUVS [ 49 – 51 ] ( S1 Fig ) across nine regions in the Indo-Pacific (n = 1,041; Fig 1A and Fig 1F ). The field survey was undertaken under ethics approval and permit RA/3/100/1166 from the Animal Ethics Committee of the University of Western Australia, following guidelines under the Animal Welfare Act 2002 (WA) and the Australian Code for the Care and Use of Animals for Scientific Purposes. The data collected involved passive observation of animals using baited video systems, and no animals were manipulated directly. BRUVS were typically deployed as longlines of five rigs, where each rig was suspended at 10 m, 200 m from its nearest neighbor and within 300 km of the nearest coastline. The entire line was left to drift freely for 2 hours. The rigs are made up of a vertical pole and a horizontal crossbar that supported two GoPro underwater action cameras [ 52 ]. The two cameras converged with an inward angle of 8 degrees on a bait canister, suspended at ca. 1.5 m from the cameras at the end of an adjustable arm. The bait canister contained 1 kg of crushed sardines (Sardinops spp.). The recorded video footage allows taxonomic identification of individuals and estimates of relative abundances as the maximum number of individuals of a given species in a single frame (MaxN) [ 53 ]. As with any sampling methodology, BRUVS are unlikely to fully capture the species pool [ 54 ]. However, BRUVS have been widely used to generate reliable and consistent estimates of richness, size, and abundance. Moreover, despite interspecific differences in bait response and animal mobility, and variation in bait plumes [ 55 ], BRUVS remains one of the most reliable ways of standardized sampling of large predators, such as sharks, for testing spatial and temporal variation. The BRUVS sampled both the midwater assemblages over a range of seabed depths and conditions, including near coastal habitats (<100 m depth, <50 m from the coast; n = 201), raised banks and shallow shoals remote from the coast (<300 m depth, >10 km, n = 199), deep and shallow seamounts (summit 1,100 m and 70 m depth, respectively; n = 156), and abyssal plains (>2,000 m depth, n = 80).

We modeled the spatial variation in three attributes for each individual BRUVS, which reflects different aspects of the predator community as recorded on video ( S2 Fig ). Vertebrate species richness ( Fig 1B ) was recorded as the total number of vertebrate species observed per 2-hour deployment. We calculated mean maximum body size for each deployment, weighting for abundance [ 57 ] (hereafter simply “body size,” Lmax in cm; Fig 1C ). This attribute is commonly used for assessing the state of coral reefs, as an indicator of the overall fish and shark community [ 57 , 58 ], and the degree to which the trophic pyramid is dominated by large individuals and species. Ecosystems with abundant and large individuals tend to exert greater top-down control and require high nutrient input [ 8 ]. In addition, body size is a highly sensitive indicator of fishing pressure [ 59 ], as larger individuals and species are preferentially targeted and removed. Body size for each deployment was computed using the following relationship: in which Lmax i is the maximum length recorded for species i according to FishBase [ 60 ] or in the literature, and MaxN i is the MaxN abundance, the maximum amount of individual observed during the 2-hour recording for species i, and n is number of species. Finally, we estimated the total relative abundance of sharks, across all species. While the ecological roles of reef sharks as apex predators remain a topic of debate [ 9 ], sharks are particularly vulnerable to exploitation and emblematic symbols of conservation [ 16 , 61 ]. Moreover, they are recognized as important indicators of marine health, potentially controlling lower trophic levels [ 9 , 62 , 63 ] through trophic cascades in both reef [ 64 ] and pelagic systems [ 65 ]. Shark abundance was calculated as the sum of MaxN across all shark species for each deployment ( Fig 1D ) and modeled as log[sumMaxN + 1].

Drivers of predator diversity and abundance

We examined relationships between the predator attributes and spatial drivers classified broadly under three categories: environment, geomorphology, and human pressures. Hypothesized environmental drivers extracted for each deployment were i) median SST (22−29.19°C, NOAA's Multiscale Ultra-high Resolution [MUR] SST http://coastwatch.pfeg.noaa.gov/erddap/wms/jplMURSST/index.html), a proxy for latitudinal patterns in species diversity universally observed across taxa [17]; ii) SST standard deviation (0.53−2.36°C), an indicator of frontal dynamics generating nutrient mixing and multilevel productivity [34,66]; iii) median chlorophyll-a concentration (0.03−1.15 mg m−3, 8-day AQUA MODIS http://coastwatch.pfeg.noaa.gov/erddap/wms/erdMHchla8day/index.html), an indicator of primary productivity and available trophic energy [25]; and iv) distance to the center of the Coral Triangle (211−6,667 km), the epicenter of fish diversity [67]. Geomorphological drivers were i) seabed depth (6−3,638 m), a dimension that fundamentally structures and constrains marine habitats vertically [68]; ii) distance to the nearest coast (0−326 km), a measure of terrestrial energy availability [24] and a physical barrier restricting the horizontal extent of the marine habitat [69]; and iii) distance to the nearest seamount [70] (summit depth <1,500 m, 1.5−505 km), the presence of which is known to attract predators [26].

For each deployment, we quantified human pressure using a range of metrics. Total industrial fishing effort for all gear were estimated for each deployment (https://globalfishingwatch.org/). We used published records of fishing hours [37]. In these records, fishing hours were estimated from AIS vessel position fixes and algorithms that determine fishing behavior based on movements. AIS records is limited to vessels greater than 15 m and can be unreliable, both in terms of people turning it on/off and falsifying records when they do not want to be monitored, as well as in the frequency of transmissions. We extracted averaged values over 0.5° (hours fished, 0–0.8 hr km−2), corresponding with the sampling period (2012–2016). We also computed the minimum distance to the nearest human population using the LandScan 2016 database (0.1−829 km)[71], and minimum distance to the nearest human density center (hereafter “market,” 11–1,450 km), using the World Cities spatial layer (ESRI). This layer defines human density centers as provincial capital cities, major population centers, landmark cities, national capitals, and shipping ports. These two distance metrics are derived to indirectly capture the many cumulative effects which humans have on ecosystem predators [32,72] including noise pollution [73], nonreported fishing [1], vessel strikes [74], infrastructure development [72], and direct exploitation [75]. Moreover, these metrics encompass some aspects of the historical impacts that have occurred before the onset of modern record keeping [1,76]. We explored both distance to population and distance to market because while small populations can have notable impacts on regional predators [77], pressures scale substantially when supported by an industrialized market [78]. For each deployment, we also estimated the human population in a 50 km and 500 km buffer region (0−111,295 and 0−3,018,935 humans, respectively), and the human development index of the nearest country (HDI, 0.61–0.93, http://hdr.undp.org/en/statistics/hdi/), which takes into account health and education status. Finally, we tested the impact of management by taking advantage of the different protection level implemented for each deployment (Fig 1A). Using MPA coverage from the World Database of Protected Areas [79], we assigned each a protection category corresponding to whether it was 1) unprotected and open to fisheries (n = 340), 2) inside a small no-take MPA (IUCN class I–IV, <1,000 km2) or inside an MPA that allowed some extractive pressure (IUCN class III–IV, n = 311), or 3) inside a large no-take MPA (IUCN class I and II, >1,000 km2, n = 390). These three broad categories offer incrementally more effective and strict protection on predators. Large no-take MPAs were assessed specifically since previous studies have documented that they meet some of the unique conditions (both large and no-take) necessary for protecting large species [6], choosing 1,000 km2 as a conservative threshold [80]. The protection categories were unbalanced and did not vary independently with distance to human market or distance to human populations. We note that the Commonwealth Marine Reserve network inside the Australian EEZ have recently undergone a review (2014−2015) of the reserves implemented in 2012 (https://www.environment.gov.au/marinereservesreview/about). These new changes have been implemented and are locked in for the next 10 years.

We used BRTs [22] to estimate the relative strengths of the effects of environmental conditions, geomorphology, and human pressures on the three predator attributes. BRTs can detect nonlinear relationships between response variables and their drivers, e.g., SST and species richness [17]. Further, BRTs are robust to codependencies amongst drivers, which are common in ecology. Codependencies can arise when the effect of a driver is conditional on another driver meeting a certain value. For example, the net effect of seamounts on predators aggregation is highly conditional upon regional frontal features and eddies [81]. Finally, BRTs are considered reasonably robust to collinearity, arising from correlated drivers. For example, the absence of human populations in the middle of the ocean renders distance to coast and distance to human population correlated in our data (r = 0.84).

To select the best BRT model, we chose the optimal combination of tree complexity, learning rate, and bag fraction as the one minimizing the out-of-bag (OOB) estimates of error rate [82]. The bag fraction term introduces stochasticity into the models and controls overfitting [83]. Model features were chosen depending upon goodness-of-fit via a cross-validation (CV) procedure. The contribution of each driver (%) was estimated as the proportion of times each driver was selected to split the data among all the trees, weighted by the squared improvement to the model as a result of each split, and averaged over all trees. Standard deviation around the contribution of each driver was generated by a hundred random seed iterations of the model selection computation. We retained drivers with more than 10% contributions to the fully saturated model in order to generate simplified and parsimonious final model [22]. Model parameters for the best models are reported in S2 Table.

Due to complex interactions between bait diffusion rates, current speed, and fish attraction, it is difficult to determine the sampling range of the individual BRUVS. Here, a separation of 200 m between each BRUVS on the same string was a trade-off between practicalities in the field and maximizing the distance between each rig. This distance may be insufficient to guarantee independence [84]. In order to account for potential spatial autocorrelation, we introduced a spatial autocovariate term [85], calculated from the residuals of our simplified BRTs. The residual autocovariate was calculated by arranging the residuals of the simplified BRT models on 0.001 degree grid (111 m) and applying a focal mean using a first-order neighborhood [86]. This approach has been shown to significantly reduce the spatial autocorrelation in the model residuals [87]. Model residuals of the final autologistic models were checked for spatial autocovariation, using a Moran’s I test. For out-of-sample predictions, the autocovariate was set at the median value. For purposes of reporting the percent contribution of each driver, the percentage contribution was rescaled without the contribution of the autocovariate term.

All BRTs were built in R (R Development Core Team 2011 version R version 2.15.2) using the gbm package version 1.6–3.1 and custom code available online (http://cran.r-project.org/web/packages/gbm). To detect potential thresholds within the relationship between the three attributes and spatial drivers, we tested for the presence of nonlinear relationships. When the null hypothesis of no change of slope was disproven (Davies’s test) [88], we performed breaking point regressions [89] to identify the threshold values. BRT model predictions were rendered on a 10-km by 10-km grid between 7° N and 32° S and between 30° E and 75° W for each predator attribute, showing the world’s land masses and countries’ EEZs [90]. Out-of-sample predictions were restricted to sites where conditions were similar to the conditions of the deployment sites (S2 Fig). Potential hotspots were defined as the highest 5% of predicted values.