Submerged Carbonate Banks Aggregate Pelagic Megafauna in Offshore Tropical Australia

Bouchet, Phil J.; Letessier, Tom B.; Caley, Michael Julian; Nichol, Scott L.; Hemmi, Jan M.; Meeuwig, Jessica J.

doi:10.3389/fmars.2020.00530

ORIGINAL RESEARCH article

Front. Mar. Sci., 30 June 2020

Sec. Marine Megafauna

Volume 7 - 2020 | https://doi.org/10.3389/fmars.2020.00530

Submerged Carbonate Banks Aggregate Pelagic Megafauna in Offshore Tropical Australia

$\r\nPhil J. Bouchet,,*$ Phil J. Bouchet^1,2,3*

Tom B. Letessier^1,4

Michael Julian Caley^5,6

Scott L. Nichol⁷

Jan M. Hemmi¹

Jessica J. Meeuwig^1,8

¹School of Biological Sciences, The University of Western Australia, Crawley, WA, Australia
²School of Mathematics and Statistics, University of St Andrews, St Andrews, United Kingdom
³Centre for Research into Ecological and Environmental Modelling, University of St Andrews, St Andrews, United Kingdom
⁴Institute of Zoology, Zoological Society of London, London, United Kingdom
⁵School of Mathematical Sciences, Queensland University of Technology, Brisbane, QLD, Australia
⁶Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers, The University of Melbourne, Melbourne, VIC, Australia
⁷Geoscience Australia, Canberra, ACT, Australia
⁸UWA Oceans Institute, University of Australia, Crawley, WA, Australia

The conservation of marine biodiversity is firmly embedded in national and international policy frameworks. However, the difficulties associated with conducting broad-scale surveys of oceanic environments restrict the evidence base available for applied management in pelagic waters. For example, the Oceanic Shoals Australian Marine Park (AMP) was established in 2012 in a part of Australia’s continental shelf where unique topographic features are thought to support significant levels of biodiversity, yet where our understanding of ecological processes remains limited. We deployed mid-water baited remote underwater video systems (mid-water BRUVS) in the Oceanic Shoals AMP to provide the first non-extractive baseline assessment of pelagic wildlife communities in the area. We used these observations and high-resolution multibeam swaths of the seafloor to explore potential relationships between prominent geomorphological features and the (i) composition, (ii) richness, and (iii) relative abundance of pelagic communities. We documented 32 vertebrate species across three sampling areas, ranging from small baitfish to large sharks and rays, and estimated that up to nearly twice as many taxa may occur within the region as a whole. This highlights the Oceanic Shoals AMP as a reservoir of biodiversity comparable to other documented offshore oceanic hotspots. Our results also confirm the AMP as a possible distant foraging destination for IUCN red listed sea turtles, and a potential breeding and/or nursing ground for a number of charismatic cetaceans. Model outputs indicate that both species richness and abundance increase in proximity to raised geomorphic structures such as submerged banks and pinnacles, highlighting the influence of submarine topography on megafauna distribution. By providing a foundational understanding of spatial patterns in pelagic wildlife communities throughout a little studied region, our work demonstrates how a combination of non-destructive sampling techniques and predictive models can provide new opportunities to support decision-making under data shortage.

Introduction

Agencies responsible for the management of biodiversity must decide how to invest limited resources to maximize conservation gains (Meir et al., 2004). Doing so entails strategic trade-offs and value judgments about which aspects of the biosphere should be safeguarded as a priority, and which management actions are most sensible, socially acceptable, financially viable, urgently required, and/or likely to yield the greatest return on investment in light of current capacities and pressures (Bottrill et al., 2008). In order to be effective, such triaging requires adequate knowledge of wildlife biogeography and spatial biodiversity patterns. Without it, the magnitude of what is at risk from anthropogenic threats is difficult to appraise, as is our capacity to avert biodiversity loss (Caley et al., 2014).

Despite efforts to count and document all extant species on Earth (e.g., Census of Marine Life)¹, global estimates of species numbers, including marine ones, have yet to converge (Caley et al., 2014). To date, only a fraction of extant organisms have been discovered (the so-called “Linnean shortfall”), and fewer still have been comprehensively studied (the “Wallacean shortfall”) (Whittaker et al., 2005). Extensive knowledge gaps therefore persist, especially in remote pelagic environments where sampling is often expensive and logistically difficult (Richardson and Poloczanska, 2008; Webb et al., 2010). For instance, the Timor Sea in the tropical Australian northwest is likely significant for marine mammals (Schipper et al., 2008; Double et al., 2014; Whiting, 1999), elasmobranchs (Lucifora et al., 2011; Guisande et al., 2013; Yon et al., 2020), teleost fishes (Edgar et al., 2014; Moore et al., 2017), seabirds (Clarke et al., 2011; Lavers et al., 2014), reptiles (Whiting et al., 2007; Eifes et al., 2013; Lukoschek et al., 2013; Letessier et al., 2015a), and several other taxa (Tittensor et al., 2010; Selig et al., 2014; Gagné et al., 2020), yet remains data-deficient and under-explored from a biological perspective (Moore et al., 2016). With a few exceptions (Lavers et al., 2014; Palmer et al., 2017; Thums et al., 2017), scant information exists on pelagic species diversity, habitat use, or abundance in the area – despite all of these attributes having been identified as Essential Ocean Variables (Miloslavich et al., 2018) and Essential Biodiversity Variables (Jetz et al., 2019) that are key to detecting spatio-temporal changes in marine biodiversity in the face increasing cumulative human impacts (Halpern et al., 2015, 2019). This dearth of information hampers ecological monitoring efforts at local and regional scales, and undermines appropriate evaluations of management performance within the recently established national network of Australian Marine Parks (AMPs)². Novel approaches to data collection and analysis are thus critical to addressing existing knowledge gaps and determining the degree to which pelagic communities are represented within offshore AMPs around Australia and beyond (Letessier et al., 2019).

We present the first detailed description of pelagic wildlife occurring within the Oceanic Shoals AMP (ca. 127.5°E, 11.5°S), based on mid-water baited remote underwater video systems (mid-water BRUVS). BRUVS are now a prominent element of Australia’s strategy for standardized ocean monitoring (Bouchet et al., 2018; Przeslawski et al., 2019), and offer a robust, cost-effective, and non-invasive approach to sampling that overcomes the traditional shortcomings of diver-reliant surveys (Langlois et al., 2010; Murphy and Jenkins, 2010). Extensively used in studies of benthic assemblages inside and outside protected areas (Mallet and Pelletier, 2014), they have recently been engineered for mid-water applications (Santana-Garcon et al., 2014a; Bouchet and Meeuwig, 2015; Letessier et al., 2019), making them a suitable choice for surveying the neritic (10–250 m) depths of the Oceanic Shoals AMP. Our objectives were to: (1) document the pelagic wildlife inhabiting the western part of the park, (2) estimate total species richness in the area, and (3) quantify spatial variation in community composition, species richness, and relative abundance, in relation to seafloor geomorphology. These metrics are often broadly used to guide the establishment of marine parks and monitor their effects (Soykan and Lewison, 2015), and can be meaningful for managing zoning within them (Espinoza et al., 2014). Relative to other AMPs of comparable size, the Oceanic Shoals is one of the most geodiverse parks in the Australian marine estate (Heap and Harris, 2008). Its complex topography and extensive fields of submerged limestone banks/pinnacles, terraces, and shoals are believed to act as focal points for marine megafauna in an otherwise largely oligotrophic seascape, although little empirical evidence of such habitat preferences currently exists but see Thums et al. (2017). Here, we build statistical models to test this hypothesis and map predicted patterns in species richness and abundance. In doing so, we offer the first blueprint for understanding pelagic biogeography within the AMP. This information represents an important basis for future hypothesis-driven studies, and may aid the management of human activities in a remote and relatively understudied region subject to mounting pressures.

Materials and Methods

Study Region

The Oceanic Shoals AMP is a large and shallow park (71,740 km², mean depth = 85 m, range = 10 to 250 m) that straddles the administrative boundary between the state of Western Australia and the Northern Territory, ca. 160 km northwest of the city of Darwin. Lying almost exclusively (99.7%) on the continental shelf, it is the sixth most geodiverse park in the Australian network (Harris et al., 2005), and connects two prominent and physically heterogeneous seafloor features: the Sahul Shelf in the west and the Van Diemen Rise to the east (Figure 1). These are characterized by flat-topped carbonate banks/pinnacles and terraces flanked by steep slopes and dissected by narrow, sinuous channels and valleys (van Andel and Veevers, 1967). Banks and pinnacles (hereafter “banks”) are knoll-shaped, hard-substrate structures, sometimes exceeding tens of kilometers in circumference, which rise to depths of 20–25 m (Hovland et al., 1994; O’Brien et al., 2002) and are composed of alternating layers of unconsolidated gravel and sandy sediments with high reactive organic matter (Przeslawski et al., 2011). Banks are thought to have developed from hydrocarbon seepage, and lay embedded in a matrix of uniform, featureless plains covered by silty clays and fine calcareous muds (van Andel and Veevers, 1967). Globally, the region’s terrain reflects a drowned estuarine lowland submerged ca. 18,000 years ago (in the Late Pleistocene to Holocene) following post-glacial rises in sea level (West, 1992). Strong seasonal fluvial runoff and intense tidal motions contribute to high levels of bed shear stress, particle transport and turbidity year-round (Porter-Smith et al., 2004; Nichol et al., 2013). Sporadic upwelling from the Timor Trough to the north sustains a large diversity of habitat-forming, filter-feeding heterotrophs such as corals, sponges, and gorgonians (Przeslawski et al., 2015).

FIGURE 1

Figure 1. Map of the study region showing the survey areas (numbered 1 to 3) in relation to seabed geomorphology and the boundaries of the Oceanic Shoals Australian Marine Park (AMP). Bathymetry appears in gray, and black circles represent previously known carbonate banks. The AMP comprises four zones, namely a Multiple Use zone (IUCN VI, dashed line), a Special Purpose zone (Trawl) (IUCN VI, dotted line), a Habitat Protection zone (IUCN IV, solid line), and a National Park zone (IUCN II, shaded area). Human activities such as recreational and commercial fishing, mining, aquaculture, or shipping are permissible in the former two zones, whereas mining is prohibited in the Habitat Protection zone, and the National Park zone receives full protection.

Data Collection

Field Survey

A 21-day interdisciplinary field expedition (GA0339/SOL5650) to the Oceanic Shoals AMP was undertaken on the RV Solander in September/October 2012 by the Australian Institute of Marine Science, Geoscience Australia, the University of Western Australia, and the Museum and Art Gallery of the Northern Territory. Details of daily scientific activities were captured in a blog available at https://www.nespmarine.edu.au/rv-solander-blog. Three approximately equidistant and equal-sized areas (mean ± SD = 164.9 ± 8.3 km²) were targeted for surveying (Figure 1). Each was selected to capture a range of depths, seabed facies, and geomorphic units (Figure 2; Heap and Harris, 2008), and was sampled with a collection of geophysical, geochemical, and biological techniques, including both seabed and mid-water BRUVS (Nichol et al., 2013). Bathymetric and seabed acoustic reflectance (backscatter) data were also acquired using a Kongsberg EM3002D (300 kHz) multibeam sonar system to provide 100% coverage of the seafloor at high spatial resolution (1 m, i.e., a resolution 250× finer than the best available bathymetric maps) (Caley et al., 2015). The multibeam swath maps were interpreted in tandem with in situ sedimentological samples to generate polygon shapefiles of the region’s underlying seabed geomorphology. Geomorphological classes were as follows: bank/pinnacle, depression, mound, plain, scarp, and terrace (Figure 2). Full details of field protocols and site characteristics are given in Nichol et al. (2013).

FIGURE 2

Figure 2. Bathymetry and geomorphology of the three survey areas (1–3, from left to right). Three-dimensional images of the seafloor can be seen in the upper panel. The lower panel shows the distribution of geomorphic units, interpreted from multibeam mapping. Carbonate banks cover 40.8, 19.3, and 28.2% of each area, respectively. Black circles mark the positions of mid-water BRUVS sites. Arrows indicate the direction of compass north.

Incidental visual point count surveys of air-breathing vertebrates (e.g., seabirds, cetaceans) were also carried out using compass-fitted Bushnell 7 × 50 reticle binoculars, both within and in transit between sampling areas. Whenever possible, radial distances and absolute bearings to sightings were measured, and the animals’ locations inferred from the known height of the viewing platform (top deck, 7 m above water line), time-stamped computer records of the ship’s GPS coordinates, and established trigonometry theory (Lerczak and Hobbs, 1998). Group size, group composition (if identifiable), and dominant behaviors were logged for each species (Supplementary Table S1). Due to the opportunistic nature of the counts, the resulting data were excluded from the analysis, but are reported in Supplementary Table S1.

Pelagic Videography

Mid-water BRUVS were deployed at multiple stations inside each sampling area (n₁ = 38 in Area 1, including 45% on banks, n₂ = 38 in Area 2, including 34% on banks, and n₃ = 40 in Area 3, including 50% on banks) during daylight hours (07:00–17:00 GMT+9.30). As per Letessier et al. (2013), the instruments were suspended at a fixed depth of 10 m and moored with a plow anchor (Supplementary Figure S1). Two high-definition Hero2 GoPro cameras (set to wide field-of-view mode) mounted in a stereo-pair configuration and fitted with battery packs and 32 GB memory cards allowed the recording of wildlife activity. These are low-cost alternatives to traditional handheld video technologies (e.g., Sony HDR-CX12) that have proven effective for the ecological monitoring of natural fish communities (Letessier et al., 2015b). Deployments lasted for an average of 190 min (SD = 28 min, range = 46 to 217 min). Although soak times are conventionally shorter in baited video research, no studies had ever investigated optimal deployment durations for mid-water BRUVS at the time of the survey. Preliminary data from a pilot study conducted off Dirk Hartog Island earlier that year indicated that deployments of 120–135 min were insufficient for capturing all species (i.e., no flattening of the associated species accumulation curve) (Letessier et al., 2013). Longer soak times were also deemed necessary to allow cautious animals that first remain in the distance to more closely enter the field of view, which is critical to improve the identification of species in turbid environments such as those of the Timor Sea. To minimize pre-soak disturbance (e.g., from engine noise), mid-water BRUV units were launched while the vessel was underway. Bait arms were attached to the main frames, and made of stainless steel with a perforated PVC-coated bait canister fastened to one end (ca. 1.5 m from the cameras). Canisters were filled with pilchards (Sardinops sagax, 2.5 kg, minced using an electric meat processor), in line with conventional protocols employed throughout Australasia (Wraith et al., 2013; Bouchet et al., 2018).

Station locations were determined according to a generalized random tessellation stratified design (GRTS; Stevens and Olsen, 2004) which generates statistically robust and geographically-balanced samples that can accommodate heterogeneous inclusion probabilities (Jiménez-Valencia et al., 2014). Although its implementation is still relatively novel in Australia, GRTS has been successfully harnessed as a monitoring tool in other parts of the marine park network (Hill et al., 2014; Foster et al., 2017). The design was weighted toward shallow habitats (<50 m) to maximize data collection opportunities in the vicinity of carbonate banks whilst retaining adequate spatial coverage overall (Figure 2). As no detailed charts of regional bathymetry were available prior to the survey (Caley et al., 2015), this was done based on an interpolated nine arc-second (ca. 250 m at the equator) depth grid of the country’s exclusive economic zone³. However, high turbidity and minimal ambient light made concurrent demersal surveys (i.e., benthic BRUV deployments conducted by partners from the Australian Institute of Marine Science) logistically difficult and little effective beyond 60 m (Nichol et al., 2013). Sampling plans were therefore revised in Areas 2 and 3 to limit surveys within this depth range. Stations were separated by a minimum distance of 500 m to curtail spatial auto-correlation and reduce the potential overlap of bait plumes emanating from adjacent units (Santana-Garcon et al., 2014b).

Data Treatment

Video Processing

Mid-water BRUV videos were post-processed and annotated using the SeaGIS software EventMeasure⁴. All animals were identified to the lowest possible taxonomic level, and counted to yield standard assemblage metrics including species richness and MaxN, a conservative index of relative abundance that avoids double counting and corresponds to the maximum number of individuals of each species visible at any one time (Priede et al., 1994). MaxN values were then summed to generate a measure of total relative abundance for each deployment (Total MaxN) (Rees et al., 2014). Because sympatric Australian (Carcharhinus tilstoni) and common (Carcharhinus limbatus) blacktip sharks are known to hybridize and co-occur in equal frequencies (Ovenden et al., 2010; Morgan et al., 2012) but cannot be confidently distinguished other than by genetic sampling, we labeled all blacktip sharks (and potential hybrids) as C. tilstoni (Espinoza et al., 2014). Furthermore, Australian blacktip, spot-tail (Carcharhinus sorrah) and graceful (Carcharhinus amblyrhynchoides) sharks can be difficult to tell apart on video and were thus lumped together as a tilstoni/sorrah/amblyrhynchoides species complex. Narrow-barred (Scomberomorus commerson) and broad-barred (Scomberomorus semifasciatus) mackerels were grouped together for the same reasons (Supplementary Figure S2).

Trophic Levels

We also retrieved the trophic level (TL) of all species from FishBase (Froese and Pauly, 2018) using the R package rfishbase (Boettiger et al., 2012; see Supplementary Table S2). However, TLs are not documented or known for all species. Where TL values were missing, species were assigned the mean trophic level of all their congeners (or all species within the same family, as appropriate) whose known distributions encompassed the Eastern Indian Ocean (FAO Area 57) and/or the Western Central Pacific (FAO Area 71). FishBase contains two trophic position metrics (FoodTroph and DietTroph), and we relied on FoodTroph because of its higher abundance of records compared to DietTroph. Cetaceans and reptiles are not included in FishBase and were assigned TL values based on expert knowledge and a review of the published literature. We assigned killer whales (Orcinus orca) a TL of 4.5 (Pauly et al., 1998), sea snakes (Hydrophiidae spp.) a TL of 4.0 (based on Hydrophis ornatus as a best approximation of taxonomic identity), and olive ridley sea turtles (Lepidochelys olivacea) a TL of 3.1 (Peavey et al., 2017).

Geomorphometrics

A suite of nine geomorphometrics was derived from the multi-beam bathymetry grids (Supplementary Figure S3), and chosen to capture a range of meaningful terrain attributes (Bouchet et al., 2015). Among them were indices of seabed aspect, slope, and curvature, as well as several measures of terrain variability retained to encapsulate independent facets of topographic complexity, and calculated on both small (100 m) and large (1000 m) scales (Lecours et al., 2017; Supplementary Table S3). The topographic position index (TPI) was preferred over the relative difference to mean value (Lecours et al., 2017) as it could be more easily calculated at different spatial scales, given available software. Sea surface temperature, salinity, and fluorescence (a proxy for primary productivity) were also initially considered but not retained due to collinearity, inadequate resolution, insufficient spatial coverage, or a combination of the above. Each raster grid was upscaled to a 50 m resolution via bilinear interpolation before analysis. This was accomplished both to speed up computing and because the focal algorithm used to calculate seabed curvatures was mathematically constrained to a fixed 3 × 3 neighborhood around a center cell (i.e., giving a 3 m window on a 1 m resolution raster initially). This scale was deemed too fine to match the processes likely affecting the distribution of mobile pelagic species. The choice of 50 m was arbitrary, but seen as an acceptable compromise between minimizing information loss and maintaining a resolution that is appreciably greater than the best available national bathymetry dataset for the region (Caley et al., 2015).

Data Analysis

All analyses were undertaken in R v3.6.0.

Species Rarefaction

We constructed a sample-based rarefaction curve (Gotelli and Colwell, 2001) from the mid-water BRUVS data using the R package iNEXT (Hsieh et al., 2016). The curve gives the statistical expectation of the corresponding species accumulation, and was extrapolated to N = 5,000 individuals using the numerical derivations presented in Colwell et al. (2012). We assumed observations originated from a wider statistical universe (i.e., the entire theoretical assemblage), and generated robust unconditional confidence intervals that do not converge to zero at the full-sample end of the curve (Colwell et al., 2012).

PERMANOVA

We tested for differences in the structure of pelagic communities detected on, and away from, banks (coded as a fixed, two-level factor, geom) using a one-way non-parametric multivariate analysis of variance (PERMANOVA) implemented in the vegan package (Oksanen et al., 2018). A large body of research illustrates the benefits of this technique in teasing out ecological patterns from BRUVS-based experiments conducted across a variety of habitat conditions and water column positions (e.g., Claudet et al., 2010; Gladstone et al., 2012; Zintzen et al., 2012; Santana-Garcon et al., 2014a). Computations were performed on fourth root-transformed data and based on the zero-adjusted Bray Curtis dissimilarity measure, which ignores (uninformative) joint absences (Clarke et al., 2006). Variance partitioning was achieved after conditioning upon deployment length as an input covariable. Robust tests of individual terms were obtained by permutation under a reduced model (n = 9,999; type III conditional sums of squares). If any significant p-values (α = 0.05) were returned, a test for homogeneity of dispersions (PERMDISP) was done, and a similarity percentage routine (SIMPER) run to gauge the contribution of each species to discriminating assemblages as a function of geomorphology (Wakefield et al., 2013).

Spatial Modeling

We developed generalized additive mixed models (GAMMs) of both species richness (henceforth S) and Total MaxN (henceforth N) (Pedersen et al., 2019) in the mgcv R package (Wood, 2004) assuming Poisson distributions with log link functions, as is appropriate for count data. Our focus here was only on higher-order megafauna, which we defined as those species with TLs ≥ 3.5, in keeping with Tremblay-Boyer et al. (2011) and Lassalle et al. (2012). We chose GAMMs due to their success in modeling non-linear species-environment relationships (Brodie et al., 2020) and in trading off good predictive performance in unsampled conditions with reasonable ecologically intelligibility (Derville et al., 2018). Our underlying expectation was that both S and N would increase toward shallower banks (i.e., lower depth, higher backscatter), where animals supposedly aggregate. To test for this effect, we coded depth and backscatter as a tensor product term (Wood, 2006). We adopted a full-subsets approach, whereby all combinations of candidate geomorphometrics were considered up to a maximum of four uncorrelated terms per model (Pearson correlation cut-off <| 0.7|) (Lambert et al., 2014) and ranked according to their second-order Akaike’s Information Criterion (AICc) scores. We limited basis sizes to k = 4 (Best et al., 2012) in order to capture dominant ecological relationships (Mannocci et al., 2017) and we relied on restricted maximum likelihood as the criterion for model optimization as it penalizes overfitting and leads to more pronounced optima (Wood, 2011). Each model was built using thin plate splines and included an offset term for effort (i.e., the log of the duration of the video clips, in minutes), as well as a random effect term identifying each sampling area. All outputs were inspected for overdispersion, autocorrelation, and violations of residual assumptions. Finally, we created combined prediction maps from the GAMM models for each sampling area. We did this by adapting the methods of Dunstan et al. (2012) and splitting predictions of S and N into deciles, before assigning to each grid cell an additive score equal to the normalized sum of its decile ranks. This approach provides a simple approximation of the bivariate density distribution of richness and abundance values. While the decision to use decile splits was arbitrary, it was made to allow a reasonable trade-off between interpretability and the ability to capture patterns in the bivariate density (Dunstan et al., 2012). The resulting index goes from zero (low S and low N) to one (high S and high N), and can be seen as an index of pelagic “hotspots”. To mitigate extrapolation, predictions were constrained within the univariate range of each covariate (Sequeira et al., 2018a) and within an alpha-hull encompassing depth and backscatter values at sampling sites.

Uncertainty Propagation

We assessed uncertainty due to both model selection and parameter estimation using Monte Carlo techniques, as described in Wenger et al. (2013). We started by resampling the BRUVS data (non-parametric Bootstrap with replacement, n = 100 times) and building GAMMs for each replicate dataset, as per the methods above. With every run, we recorded which model minimized the AICc, allowing us to build a frequency table of the best model formulations (Supplementary Table S4). We then randomly selected a model out of this best set, with a probability equal to its frequency value. This strategy is consistent with current practices in Bayesian model averaging (Wenger et al., 2013). Next, we randomly drew values from the multivariate distribution of all fixed effects included in the chosen model. We used these to build replicate GAMMs, keeping track of all bootstrap outputs and predictions, and quantifying uncertainty as the robust coefficient of variation (rCV) across Bootstrap replicates. rCV is a modified form of the traditional CV which is less sensitive to skewness and outliers (Arachchige et al., 2019) and is calculated as rCV = 100 × normalized IQR/median, where IQR is the interquartile range (25–75 percentile) and the normalized IQR is 0.7413 × IQR (Temnerud et al., 2013).

Results

Species Richness

In a total of 368.5 h of sampling effort, we recorded 1,693 individuals from 30 identified species (corresponding to 13 taxonomic families), nearly a third of which were sharks and rays (Supplementary Table S2 and Supplementary Figure S4). The most abundant species were mackerel scads (Decapterus macarellus), unicorn leatherjackets (Aluterus monoceros), and sharks from the Carcharhinus complex, together accounting for 84% of N. Among sharks, species from the Carcharhinus complex were also the most prevalent, being observed at 44% of all sampled sites. One mobulid (the oceanic manta ray, Manta birostris), one cetacean (the killer whale, Orcinus orca) and one hydrophiid species (likely the ornate reef sea snake, Hydrophis ornatus, although identification was uncertain) were also filmed, as well as three olive ridley sea turtles (Lepidochelys olivacea), as detailed elsewhere (Letessier et al., 2015a). Both species richness and total MaxN were highest in Area 2 (S = 19, N = 589), followed by Area 3 (S = 18, N = 480) and Area 1 (S = 17, N = 570). S varied from zero to eight across sites, with a mean of 2.9 ± 1.5 SD. Abundance ranged from zero to 84 across sites, with a mean of 14.4 ± 15.8 SD. Three species of cetaceans were additionally documented in 17.5 h of opportunistic visual surveys: the bottlenose dolphin (Tursiops truncatus), the killer whale, and the false killer whale (Pseudorca crassidens) (Supplementary Figure S5). Their pods ranged in size from three to 25 individuals across species, and frequently comprised calves or juveniles seen resting, foraging and socializing (Supplementary Table S1).

Analytical rarefaction of the mid-water BRUVS samples yielded estimates of species richness per standardized sample of 40 and 50 individuals of S₄₀ = 6.7 (95% CI = 6.3–7.1) and S₅₀ = 7.4 (95% CI = 6.9–7.9), respectively (Figure 3), in line with previous studies conducted off Australia on the same array of taxonomic groups within a similar latitudinal range (Worm et al., 2003, 2005; Morato et al., 2010). Extrapolation suggested an asymptote of 38 species (95% CI = 28–48). A sample-based rarefaction based on species incidence data (see R code) indicates that an additional ca. 411 h of video footage would be required to document this number of species.

FIGURE 3

Figure 3. Individual-based rarefaction (solid line) and extrapolation (dashed line) of species richness (S) from mid-water BRUVS samples, with 95% unconditional confidence interval (shaded area). The reference value (i.e., observed species richness) is shown as a filled black circle. Diamonds mark the estimated value of S corresponding to samples of 40 (gray) and 50 (white) individuals, for comparison with previous studies (see section “Species Richness”).

Community Structure

Pelagic communities associated with banks were significantly distinct from those elsewhere (geom, p = 0.0109, Table 1). Twelve species accounted for >90% of this dissimilarity, with mackerel scads, unicorn leatherjackets, and sharks from the Carcharhinus complex again best discriminating assemblages between banks and surrounding plains (cumulative contribution >50%, Supplementary Table S5). C. tilstoni and C. sorrah are both commercially exploited, the former being a major component of the northern Australian gillnet fishery (Ovenden et al., 2010; Harry et al., 2013), and occurred in greater numbers on banks (Av. Ab = 0.60) than away from them (Av. Ab = 0.38). By contrast, more mobile and wider-ranging species like the tiger shark (Galeocerdo cuvier) and the scalloped hammerhead shark (Sphyrna lewini) were observed in deeper waters, consistent with their known movement and diving behaviors (Duncan and Holland, 2006). However, individual contributions remained low (<10%) for most species, indicating that changes in their individual abundance, rather than occurrence, were likely responsible for the observed patterns.

TABLE 1

Table 1. Permutational multivariate analysis of variance (PERMANOVA) of the influence of seabed geomorphology on the structure of pelagic species assemblages in the Oceanic Shoals Australian Marine Park, based on zero-adjusted Bray-Curtis distance matrices derived from fourth-root transformed data. geom is a two-level fixed factor (on bank, off bank). Duration represents the duration of each video clip, in minutes.

Spatial Models

Generalized additive mixed models revealed a significant interaction between seabed depth and acoustic backscatter (p-values <0.05), with an increase in both S and N generally predicted in the vicinity of shallower, harder-substrate banks (Figures 4, 5). Some variation was apparent at the scale of individual banks, with higher values of S and N often fringing the banks’ margins (Supplementary Figures S6–S8). The best model for S explained 12.3% of the deviance, and only retained the tensor product term of depth × backscatter. By contrast, the best model for N explained 52.1% of the deviance, and included not only the tensor product term but also significant smooths for seabed curvature and the large-scale topographic position index, suggesting higher abundance over more convex and complex terrain (Supplementary Figure S9). Uncertainty was greater in model predictions of N, particularly along bank margins.

FIGURE 4

Figure 4. Tensor product smooths for the non-linear interaction between seabed depth (in meters) and acoustic backscatter (in decibels) in the models of (A) species richness, S, and (B) species relative abundance, N, minimizing the second-order Akaike’s Information Criterion (AICc). Models explained 12.3 and 52.1% of the deviance, respectively.

FIGURE 5

Figure 5. Maps of pelagic “hotspots” within sub-areas of the Oceanic Shoals Australian Marine Park (AMP), as predicted from the best GAMM models. Top: Median additive decile ranks (normalized to 0–1) calculated from predictions of species richness (S) and abundance (N) in each sampling site (1 to 3 from left to right). Yellow tones represent areas where both S and N are high (i.e., hotspots), whereas blue tones denote areas where both S and N are low (cold spots). Areas outside the covariate range of environmental predictors are left white (see main text for details). Bottom: Bivariate partitioning of S and N used in deriving the hotspot index. Gray cells denote combinations that were not observed in the data. The outlines of carbonate banks are superimposed as solid black lines. Arrows indicate the direction of compass north.

Discussion

Effective ocean conservation requires knowledge of the spatial ecology of marine species and the identification of essential habitats supporting their populations (Ward-Paige et al., 2012). We used non-destructive baited videography to begin the task of understanding patterns in the structure, richness, and abundance of pelagic communities across a vast, geodiverse seascape: the Oceanic Shoals AMP. We identified 32 taxa (30 on mid-water BRUVS, another two during visual point counts) spanning a broad array of life histories, phylogenies, and trophic guilds, from small schooling planktivores to large-bodied carnivores (Supplementary Table S2). Our estimate of asymptotic richness suggests that up to 48 species may be present in the area (Figure 3), consistent with other previously described open-ocean hotspots at similar latitudes (despite methodological discrepancies; Worm et al., 2003, 2005; Morato et al., 2010). Importantly, our work aligns with previous observations that pelagic hotspots are associated with prominent topographic features at local scales (Worm et al., 2003; Fontes et al., 2014; Vassallo et al., 2018). It also complements a growing body of research documenting the Timor Sea as a conspicuous and globally significant epicenter of biodiversity (Tittensor et al., 2010; Moore et al., 2017).

Marine organisms commonly form predictable aggregations at sites of abrupt topography such as submarine canyons or seamounts (Nur et al., 2011; Aïssi et al., 2012; Garrigue et al., 2015; Bouchet et al., 2017; Sutton et al., 2019). The combined importance of seabed backscatter and water depth in our models suggests that the region’s carbonate banks may play a similar role in attracting pelagic megafauna. This pattern is consistent with recent evidence that flatback (Natator depressus) and olive ridley turtles migrate to the complex seafloor of the Oceanic Shoals AMP to forage during the post-nesting season (Whiting et al., 2007; Thums et al., 2017). While the animals’ preferred prey are unknown, studies of sessile epifauna across the Van Diemen Rise and the eastern Joseph Bonaparte Gulf point to unique and dense mixed sponge and coral gardens on banks, mounds, shoals, and other raised geomorphic structures (Przeslawski et al., 2014). By promoting structural heterogeneity and providing substrate and shelter for numerous benthic fish and invertebrates, banks may concentrate prey and act as critical links in regional benthic-pelagic coupling pathways through enhanced predation opportunities. Similarly to Moore et al. (2017), our models also predicted elevated wildlife abundance around bank slopes and edges [although note that models of N performed far better than models of S, contrary to recent studies in other, perhaps more biologically complex, systems; Sequeira et al. (2018b)]. This could be explained by enhanced productivity caused by turbulent mixing and cryptic upwelling events likely taking place as tidal currents are directed around the banks, as evidenced by tidal scouring of the seabed pockmarks in their periphery (Picard et al., 2014). However, the region’s highly turbid waters will typically block exposure to ambient light beyond 60 m (Nichol et al., 2013), such that nutrient enrichment from physical forcing might be insufficient to fully sustain pelagic food webs where banks lie too deep to allow photosynthesis to occur. Additional research examining fine-scale interactions between ocean currents, complex seafloor geomorphology, and benthic/pelagic productivity is needed to disentangle the mechanisms underlying patterns of species occurrence inside this AMP (Moore et al., 2017).

Overall, carcharhinids were among the most diverse and dominant taxa, and we documented consistent patterns of co-occurrence between blacktip (C. tilstoni) and spot-tail (C. sorrah) sharks. These two species are generally site-faithful with small home ranges (Stevens et al., 2000). It is plausible, therefore, that resident populations may occupy the AMP year-round, with implications for the sustainability of the commercial fishing activities of which they are targets (Field et al., 2012). By contrast, tiger sharks (G. cuvier) display characteristically variable movements. While some animals undertake nomadic journeys that cover thousands of kilometers annually, others concentrate on discrete reefs for feeding and reproduction (Speed et al., 2010). In the absence of spatially and temporally explicit tracking data, the drivers of their presence in AMP are unclear. Likewise, how cetaceans use the park also remains unresolved, as no estimates of abundance, demographic trends or geographic dispersal presently exist for offshore odontocetes in the Timor Sea (Allen et al., 2012; but see Palmer et al., 2017). In this context, several observations of young calves and juveniles in multiple survey locations offer tentative preliminary evidence that several cetacean species may occupy the area at some point during their early life stages.

Due to its simplicity and intuitiveness, species richness has long remained a favored metric in conservation priority setting. A critical yet risky assumption underlying its use is that richness values can be readily compared, both spatially and temporally. In practice, sampling error, intensity, coverage, and detectability can all affect richness estimates to some degree (Fleishman et al., 2006). The usefulness of richness for monitoring has therefore been met with mounting controversy (Hillebrand et al., 2018) prompting calls to develop composite indicators of biodiversity across scales from genes through to ecosystems (Lyashevska and Farnsworth, 2012). However, in the face of uncertain budgets and unpredictable opportunities to invest in conservation, simple decision rules such as protecting richness hotspots may actually prove effective, especially when management actions are implemented incrementally over time (Meir et al., 2004). This is arguably the case for the Australian Marine Park network, which was proclaimed in 2012 but whose management zonation took more than 6 years to come into effect. Fleishman et al. (2006) argue that species richness retains value only if considered in consort with metrics of species endemism, threat, or abundance. With this in mind, we combined analyses of community structure with spatial models of both S and N, and reported on a list of species attributes (e.g., exploitation status, group composition when known, vulnerability; Supplementary Table S2) to facilitate improvements in our ecological understanding of the Oceanic Shoals AMP for monitoring and decision-making. This is an important initial step toward supporting park management, yet we caution that our data only represent a snapshot of species occurrence in a single season, meaning that any temporal trends cannot currently be resolved from this initial survey. The dynamics of species immigration, extinction, and taxonomic turnover are complex, and apparent stability in diversity patterns can mask potentially substantial changes in community composition (Hillebrand et al., 2018). Repeat surveys will thus be required before temporal trends in pelagic biodiversity within the AMP can be estimated with any confidence.

BRUVS are prone to a number of biases, including imperfect detection in turbid environments, and preferential recording of animals that are more active or defensive of the bait (Dunlop et al., 2014). Note, for example, that it is very difficult to determine the rate at which bait degrades and releases into the water column, making bait dispersal a key unanswered question in BRUVS research currently (Whitmarsh et al., 2017). The problem is exacerbated by the fact that oceanographic conditions (e.g., swell, current speed, etc.) may vary from site to site, and it is our experience that some of the bait canisters seemed to empty more quickly than others. There is no simple solution to this, and the implications of these shortcomings have been discussed at length elsewhere (Pelletier et al., 2011; Taylor et al., 2013; Mallet and Pelletier, 2014; Bouchet and Meeuwig, 2015). However, in deploying mid-water BRUVS inside the Oceanic Shoals AMP, we seized a critical opportunity to trial and validate baited videography in open ocean conditions away from coasts. In doing so, we addressed an essential methodological need (Santana-Garcon et al., 2014c), as few other non-invasive techniques would have yielded insights into pelagic communities at an equally appropriate level of detail. Large-scale monitoring programmes for offshore marine reserves demand rapidly deployable tools that can achieve high levels of replication (Letessier et al., 2019). Mid-water BRUVS are highly effective in this regard, and now indeed form part of a national strategy for standardized biological sampling in Australian waters (Bouchet et al., 2018; Przeslawski et al., 2019).

The place-based management of the open ocean requires the maintenance of biodiversity values in the face of competing industry interests. At one end of the spectrum, pristine and remote sites are often chosen as logical candidates for protection as they make conservation possible with minimal interruption or displacement of human activities (Singleton and Roberts, 2014), while buying insurance against future change. At the other end, pressing needs for threat reduction and population recovery prevail and focus is directed on chronically perturbed systems affected by high anthropogenic use. The pelagic waters of the Timor Sea are some of the least degraded “wilderness areas” on the planet, owing to historically low levels of resource extraction by international standards (Trebilco et al., 2011). Yet, ongoing offshore petroleum industry activity and recent evidence of widespread illegal fishing and ghost netting indicate rising risks of disturbance off northern Australia (Field et al., 2009; Wilcox et al., 2013; Edyvane and Penny, 2017). The capacity of the region to withstand these pressures is difficult to ascertain given current gaps in our understanding of species distributions and resilience. Furthermore, the submerged shoals and banks of the Oceanic Shoals AMP have been earmarked as key ecological features (KEFs) by the Australian Government but remain drastically under-represented in no-take sanctuaries (Moore et al., 2016). This limited protection is concerning as few marine industries in the north-west operate at depths greater than 200 m, making continental-shelf habitats such as the ones found within the AMP those most exposed to anthropogenic activities (Moore et al., 2016). The development of robust and cost-effective strategies for conducting regional-scale assessments of biodiversity assets that can support spatial planning throughout the Australian EEZ is therefore both timely and necessary. Integrating different observational techniques such as swath mapping and baited videography, as illustrated here, can be useful for disentangling species-habitat relationships and is thus a key step in addressing this shortfall.

Data Availability Statement

Example video footage can be viewed online at https://goo.gl/DGa9hE. R code for running the analyses can be downloaded from GitHub at https://github.com/pjbouchet/oceanicshoals. The pelagic BRUVS data are available through the Australian Ocean Data Network at http://data.imas.utas.edu.au/portal/search.

Ethics Statement

The animal study was reviewed and approved by the University of Western Australia’s Animal Ethics Committee (AE Permit RA/3/100/1090).

Author Contributions

MC, JM, PB, TL, and SN conceived the ideas and designed the methodology. PB, TL, and SN collected the data. JH provided software and technical assistance. PB conducted the analyses and wrote the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Funding

This research was undertaken for the Marine Biodiversity Hub, a collaborative partnership supported through funding from the Australian Government’s National Environmental Research Program (NERP). PB was the recipient of a scholarship for international fees during the course of this work. SN publishes with permission of the Chief Executive Officer, Geoscience Australia.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to all expedition members and importantly to the master and crew of the RV Solander for their expertise, dedication and hospitality during fieldwork. We thank L. Groves for logistical support as well as L. Ferreira, F. Garcia-Solis Marchant, and S. Hodge for assistance with video processing. G. Fry and K. Sanders helped with sea snake identifications.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.00530/full#supplementary-material

Footnotes

References

Aïssi, M., Fiori, C., and Alessi, J. (2012). “Mediterranean submarine canyons as stepping stones for pelagic top predators: the case of the sperm whale,” in Mediterranean Submarine Canyons: Ecology and Governance, ed. M. Würtz (Gland: IUCN), 99–104.

Google Scholar

Allen, S. J., Cagnazzi, D. D., Hodgson, A. J., Loneragan, N. R., and Bejder, L. (2012). Tropical inshore dolphins of north-western Australia: unknown populations in a rapidly changing region. Pacific Conserv. Biol. 18, 56–63. doi: 10.1071/PC120056