Regional‐scale patterns of deep seafloor biodiversity for conservation assessment

Mining and petroleum industries are exploring for resources in deep seafloor environments. Lease areas are often spatially aggregated and continuous over hundreds to thousands of kilometres. Sustainable development of these resources requires an understanding of the patterns of biodiversity at similar scales, yet these data are rarely available for the deep sea. Here, we compare biodiversity metrics and assemblage composition of epibenthic megafaunal samples from deep‐sea benthic habitats from the Great Australian Bight (GAB), a petroleum exploration zone off southern Australia, to similar environments off eastern Australia.


| INTRODUC TI ON
Although deep seafloors cover over half the planet (Ramirez-Llodra et al., 2010), few studies have examined patterns of biodiversity for the deep sea over large spatial scales. This is an important knowledge gap as anthropogenic activities are increasingly targeting deep-sea habitats (Glover & Smith, 2003;Ramirez-Llodra et al., 2011), where inadequate data limit effective environmental impact assessment (Clark, Durden, & Christiansen, 2019). There are plans to mine the deep sea for polymetallic nodules on oceanic abyssal plains, cobalt-rich ferromanganese crusts on seamounts and polymetallic sulphide deposits on volcanically active continental margins and mid-ocean ridges (Miller, Thompson, Johnston, & Santillo, 2018).
Deep-sea sedimentary basins on continental margins are also of interest to the oil and gas industry (Zou et al., 2015), while renewable energy and aquaculture increasingly look to offshore areas to dilute environmental and social concern.
Extractive lease or licence areas can be continuous over large areas as they cluster around known resources. For example, exploration areas for polymetallic nodules in the abyssal plains of the Clarion-Clipperton Zone in the NE Pacific cover millions of square kilometres and petroleum protraction areas cover the entire Exclusive Economic Zone of the United States in the Gulf of Mexico. Yet, adequate environmental baselines have been lacking for many of these areas, particularly in deeper water (Cordes et al., 2016;Kaiser, Smith, & Arbizu, 2017). We frequently lack an understanding of how biodiversity is structured at regional scales, such as how much biodiversity could be lost due to catastrophic or accumulated impacts. Our understanding of deep-sea biogeography has been "characterised more by inference than data" with the deep-sea fauna having been characterized as both relatively uniform over oceanic scales and characterized by high species turnover at smaller scales (McClain & Hardy, 2010). The paucity in biological data has often led to the use of untested physical surrogates for biogeographical mapping (UNESCO, 2009;Wedding et al., 2013;Williams et al., 2009).
The Great Australian Bight (GAB) has been recognized as one of the world's most prospective, under-explored oil and gas provinces (Begg, 2018). Lying in the centre of the long southern Australian continental margin, this sedimentary basin was predominantly formed by two progradational river deltas in the Late Cretaceous after Australia rifted away from Antarctica (Krassay & Totterdell, 2003).
Oil and gas lease blocks now cover large parts of the basin along the continental slope .
The GAB has recognized conservation values including a newly proclaimed Australian marine park . The "Benthic invertebrate communities of the eastern Great Australian Bight" also have been designated a Key Ecological Feature (KEF).
Although KEFs are not listed as matters of national environmental significance under the Australia's Environmental Protection and Biodiversity Conservation Act 1999, they are listed as conservation values in Australia's Marine Bioregional Plans (https ://www. envir onment.gov.au/marin e/marine-biore gional-plans ), and so are considered as important components of the marine environment, and are frequently referred to in environment assessments and strategic assessments. The GAB invertebrate communities' KEF designation was based primarily on a study by Ward, Sorokin, Currie, Rogers, and McLeay (2006), who reported relatively high species richness of epibenthic megafauna on the continental shelf of the GAB (0-200 m) which they attributed to the unusual carbonate sediments and overlapping SW and SE Australian faunas. Currie and Sorokin (2014) also suggested that samples of megabenthos from two slope canyons (100-2,000 m) may have elevated richness. Uncertainty about the biodiversity of the bathyal (off-shelf) benthic communities resulted in the spatial and bathymetric boundaries of this KEF being left undefined, yet left the impression that developments in the GAB would need to account for threats to a highly endemic fauna.
Few studies have examined patterns of benthic faunal composition in the GAB from seafloors deeper than the continental shelf or upper slope. The exception is Williams et al. (2018) who analysed diversity and abundance in epibenthic megafaunal samples from five transects (200-3,000 m) across the GAB. However, this study lacked comparable samples from outside the area and could not draw conclusions about comparative richness or abundance. Tanner, Althaus, Sorokin, and Williams (2018) confirmed that taxa found along these transects were typical of museum collections of epibenthic megafauna from the same region and that their affinities were with other regions along the southern coast compared to those off the eastern coast at similar depths, although the museum samples were largely derived from upper bathyal habitats (200-1,000 m).
Over a 3-year period (2015-2017), a number of voyages of the RV "Investigator" collected comparative beam trawl samples from lower bathyal (1,900-3,500 m) and abyssal (3,500-5,000 m) depths from both the GAB (IN2015_C01, IN2015_C02, IN2017_ C01) and off eastern Australia (IN2017_V03) (Figure 1). This provided the opportunity to assess whether (a) the deep-sea benthic diversity of the GAB differed from that of comparable habitats off the east coast, and/or (b) whether the assemblages differed in composition at these spatial scales. In particular, we assessed comparative abundance, richness and evenness using a recently developed community modelling method that models and predicts rank abundance distributions (RADs) from environmental covariates   Dunstan, 2010). These covariates included oceanographic variables such as seafloor water temperature, salinity and dissolved oxygen, carbon flux to the seafloor, mean annual and seasonal variation of net primary productivity at the sea surface, as well as geographical variables latitude, longitude and depth.

| Samples
All samples were collected using an identical 4-m-wide beam trawl with a 25-mm mesh net on soft-sediment substrata in several expeditions to the Great Australian Bight and the eastern continental margin of Australia on the RV Investigator ( For this study, we have restricted samples to those that were collected at seafloor depths between 1,900 and 5,000 m to ensure inter-regional comparability. The depth of 1,900 m was chosen as the lower limit so as to include two samples from the GAB, whose mean depth along the tow was slightly shallower than the target depth of 2,000 m. We restricted the taxonomic scope to the following megafauna groups that were identified by the same experts (see acknowledgements) across voyages: hexactinellid and demospongid sponges, anthozoans, barnacles, decapods, pycnogonids, F I G U R E 1 Map of sample sites (red) and net primary production (NPP) around southern and eastern Australia Environmental variables were interpolated from global datasets to produce datasets for model training and prediction (see below). The environmental data were interpolated to mid-sample locations (latitude, longitude and depth) for the model training dataset, and at a spatial resolution of 0.1° between 129 and 159°E, 42 and 23°S, and water depths between 2,000 and 5,000 m, for a prediction dataset. Prediction depths were derived from the ETOPO1 dataset (Amante & Eakins, 2009). Annual mean and standard deviation of seafloor water temperature (°C), salinity (psu) and dissolved oxygen (ml/L) were interpolated from the CARS2009 dataset (Ridgway, Dunn, & Wilkin, 2002). Mean annual net primary productivity (NPP, g C m −2 year −1 ) and the seasonal variation of net primary productivity (SVI, g C m −2 year −1 ) were generated from a vertically generalized production model (VGPM; Behrenfeld & Falkowski, 1997) using satellite-derived chlorophyll (SeaWiFS) data from the years 2003 to 2010 (see http://www. scien ce.orego nstate.edu/ocean.produ ctivi ty/). Carbon flux to the seafloor (C flux, g C m −2 year −1 ) was estimated using NPP and SVI data and a productivity export model (Lutz, Caldeira, Dunbar, & Behrenfeld, 2007).

| Statistical analyses
Most statistical analyses were performed using the R statistical environment v 3.4.3 (see Data S1 in Appendix S1), and maps were produced by QGIS v 3.4.3-MadeIra. Exploratory analyses of the data were performed using non-metric multivariate statistics. The species-site abundance data (see Table S1 in Appendix S1) were converted into density (m −2 ) measurements by dividing abundance by the sample area and then log-transformed to down-weight the influence of abundant taxa. A triangular dissimilarity matrix was created using the Bray-Curtis coefficient with the vegdist (method="bray") function, clustered using hclust(method="ward.D2") and ordinated using the non-metric multidimensional scaling (nMDS) function met-aMDS() in the R package "vegan" v2.4.5 (Oksanen et al., 2016 random permutations of sample data amongst the factor groups. We compared regional diversity patterns using rank abundance distributions using the R Package "RAD" v0.3 (Dunstan & Foster, 2011;Foster & Dunstan, 2010). This package models three com- Species richness of each sample was also calculated for a set number of individuals (the smallest number of individuals in any sample = 38) using the rarefaction function rarefy() also in the R package "vegan," which uses the algorithm of Hurlbert (1971) and the standard error procedure of Heck, Belle, and Simberloff (1975). These results were stratified into bathomes (lower bathyal vs. abyssal) and analysed using one-way ANOVA with region (SE, NE and GAB) as the categorical factor (Gotelli & Colwell, 2011) using the Excel (v2013) function "ANOVA: Single factor" in the Analysis ToolPak add-in.
The number of species that were collected in one, two or more regions was tabulated by phylum and depth strata (lower bathyal and abyssal, excluding the unmatched shallow-1,900 to 2,300 m-samples from the GAB). As species richness varies with collection effort, the number of species was adjusted to the mean number of species per sample multiplied by a standard number of samples (median = 6) per region/depth layer.

| RE SULTS
The final data matrix was 666 operational taxonomic units (OTUs) across the 49 samples (see Table S1 in Appendix S1).
Only 251 of the OTUs were assigned species-level taxonomic names, and potentially 60% of the fauna is undescribed. A hierarchical cluster analysis primarily divided samples into lower bathyal and abyssal groups (Figure 2a). Within each of these groups, regional (GAB, SE Australia, NE Australia) subclusters were evident. The exception was for the bathyal GAB samples which were divided into (a) mid-bathyal (1,900-2,300 m) and ( the Porifera (a taxonomic group for which there were no species in common between the GAB and east coast samples) but showed little difference in pattern (not shown). The seven clusters outlined above ( Figure 2) were treated as a categorical variable "biome" in subsequent analyses.
A forward selection process in the RAD modelling procedure identified second-order polynomials of carbon flux (C flux) and seafloor salinity linked to abundance (N); carbon flux and latitude (°S) with an interaction term linked to richness (S); and second-order polynomials of seafloor salinity, Seasonal Variation Index (SVI) of NPP and longitude linked to evenness (η) (Table 2, Figure 3). The model residuals are given in Figure S1 in Appendix S1. Predictive maps ( Figure 4, Table S5 in Appendix S1) resulted in elevated abundance and richness at the shallower end of the study bathymetric range (2,000-2,500 m depth), particularly around SE Australia which is characterized by elevated NPP (Figure 1). Evenness also was consistently higher around SE Australia (Figure 4).

Rarefied richness (mean estimated richness per 38 individuals)
was not significantly different between biomes when analysed as a single-factor ANOVA ( Table S6 in Appendix S1), suggesting that abundance drives regional richness patterns. The number of species that were collected from only one region varied according to phylum and depth strata ( Figure 5, Table S7 in Appendix S1), ranging from no sponge species being shared between the GAB and eastern Australia to over half the arthropod species being shared.

| Patterns of regional-scale diversity
Samples with the highest standardized abundance and richness for benthic megafauna within our study area (see Table S1 in Appendix S1) were found at depths of 2,500-3,000 m off SE Australia (SE).
In general, abundance varied inversely with latitude (lowest in the north) across both the lower bathyal (2,300-3,200 m) and abyssal (>3,200 m) depth strata, although abyssal richness was higher for sponges ( Figure 5). Richness at abyssal depths showed a similar regional pattern, but richness at bathyal depths was generally higher off NE Australia (NE) than for the Great Australian Bight (GAB).

TA B L E 2
Summary of final rank abundance distribution (RAD) models for total abundance, species richness and evenness of samples Note: AIC, Akaike's information criterion; C Flux, mean annual carbon flux to the seafloor (g C m −2 year −1 ); NPP, net primary productivity (g C m −2 year −1 ); SVI, seasonal variation of net primary productivity (g C m −2 year −1 ); Temp, seafloor water temperature (°C). a Final model.
b From models with terms higher in the table.

F I G U R E 3
Key variable responses using final RAD models (Table 2)   These patterns do not conform to the typical latitudinal diversity gradient (LDG) reported from shallow waters (Tittensor et al., 2010) and terrestrial environments (Hillebrand, 2004), where richness decreases away from tropical areas. Instead, the patterns conform to that reported by Woolley, Tittensor, et al. (2016) for ophiuroids (brittle stars), who found the peak of regional richness for the lower Carbon flux has been repeatedly recognized as the key driver of spatial patterns of LBA richness, from local to global scales (Rex & Etter, 2010;Woolley, Tittensor, et al., 2016). Carbon flux to the seafloor was the most important environmental variable driving our RAD models for both abundance and richness (Table 2, Figure 3). A band of elevated phytoplankton density occurs at austral temperate latitudes (~40°S) across the Atlantic, SE Indian and SE Pacific Oceans, including off southern Australia and the Tasman Sea (Lutz et al., 2007). Many of our samples off SE Australia were collected from directly under this phytoplankton bloom (Figure 1). Conversely, the GAB samples were located to the north of this band and our samples from the NE were from relatively oligotrophic subtropical waters (Radke et al., 2017). It must be emphasized that our carbon flux data are modelled from surface chlorophyll data and depth Previous studies have shown that the bathymetric diversity gradient (BDG) for seafloor fauna is generally unimodal, with a diversity peak in the mid-bathyal (~2,000 m) in the North Atlantic Ocean (Rex & Etter, 2010)  the mid-bathyal to abyss across oceans (Vinogradova, 1962). Our data are consistent with this pattern within each region (GAB, SE, NE). The key variable in our abundance and richness models, carbon flux, exponentially declines with depth and is again a plausible partial explanation for these patterns. The source-sink hypothesis of Rex et al. (2005) relates low faunal density to low food supply at abyssal depths. The low density causes species extinction through Allee effects which is only partially balanced by dispersal from bathyal sources on continental margins (Rex et al., 2005) or more productive abyssal areas (Hardy, Smith, & Thurnherr, 2015). Under this scenario, SE Australia could be the source of much of the diversity in more oligotrophic regions, although we did find that the GAB and NE regions contained many species not present in the SE ( Figure 5). The decline in abundance and richness at high values of carbon flux (Figure 3a,c) may be a spatial artefact of the lack of 1,900-2,300 m sites of SE and NE Australia and needs to be verified by further sampling.
There are also numerous other environmental variables that vary with depth, and it is complex to partition out their effect on diversity. While pressure doubles between 2,000 and 4,000 m within our study area, water temperature only declines by an average of 1.1°C, salinity by 0.05 psu and oxygen by 0.7 ml/L, and it is unclear to what extent these differences drive biological patterns.
Much of the environmental variation is related to the presence of distinct water masses at different depths. At 1,000 m, low-salinity Antarctic Intermediate Water (AAIW) flows westwards across the GAB in an offshore flow known as the Flinders Current (Davis, 2005;Oke, Griffin, Rykova, & Bastos de Oliveira, 2018 (Lee et al., 2019;Talley, 2013). Water masses in the North Atlantic have been shown to be spatially dynamic at decadal (Yasuhara et al., 2019) to millennial temporal scales (Yasuhara, Hunt, Cronin, & Okahashi, 2009). These flows are consistent with a potential source-sink relationship between the diversity-rich SE and the NE and GAB regions, although this requires further assessment of species range limits.
Richness is known to be driven by macro-evolutionary processes as well as ecology. The peak at upper to mid-bathyal (200-2,000 m) depths is, at least partly, due to the age of these environments. The tropical upper to mid-bathyal has relatively high diversity but a low lineage diversification rate and thus appears to have been a long-term refuge for deep-sea animals . Conversely, the abyss is characterized by low richness for many groups (Rex & Etter, 2010). From an evolutionary perspective, abyssal taxa can amount to a small disparate subset of the bathyal lineages (Christodoulou, O'Hara, Hugall, & Arbizu, 2019;O'Hara et al., 2019), suggesting multiple infrequent range expansions from bathyal to the abyss over time and little subsequent diversification. The isopod family Munnopsidae is often cited as an exception to this trend, having radiated extensively in the deep sea (Lins, Ho, Wilson, & Lo, 2012). Deep-sea faunal turnover of benthic foraminifera has occurred throughout the Cenozoic Era, possibly due to climate-driven modifications to the thermohaline circulation (Thomas, 2007).
Seasonal variation in seafloor water temperature and salinity did not explain patterns of abundance or richness in our data.
However, seasonality in net primary production (SVI, as defined by Lutz et al., 2007) was an informative variable for our model of evenness, along with salinity and longitude. Seasonality in NPP is considered important as organic matter forms aggregates (marine "snow") with increased density, leading to enhanced flux to the seafloor during blooms (Bax, Burford, Clementson, & Davenport, 2001;Rex & Etter, 2010). Carcasses and faecal pellets from seasonal salp blooms can provide two thirds of carbon input to the seafloor in the Tasman Sea (Henschke et al., 2013). In our data, evenness peaked at intermediate to high seasonality (Figure 3f,4c). This is unlike the foraminiferan assemblage North Atlantic, where marine snow is exploited by a few dominant species (Corliss, Brown, Sun, & Showers, 2009). Also, we cannot rule out some inter-annual variation in our data, as our samples were collected in both 2015 and 2017 in the GAB, but in a single 2017 expedition for the SE and NE regions.

| The Great Australian Bight
Great Australian Bight soft-sediment assemblages at LBA depths do not have elevated abundance or richness compared to equivalent habitats off Australia's eastern coast. We did not find quantitative evidence for a "Key Ecological Feature" based on elevated biodiversity metrics that extend to LBA depths.
The evidence for the presence of elevated species richness in the GAB region is largely derived from the studies of Ward et al. (2006) and Currie and Sorokin (2014). Both these studies, however, only compared total richness of their survey samples (gamma diversity) with a few ad hoc surveys in other parts of the world. Such comparisons are known to be confounded by differences in habitat heterogeneity, gear type, scale and number of samples, inter-annual variation and number of individuals (Gray, 2002). For example, claims of elevated richness for deep over shallow seas (Grassle & Maciolek, 1992) have been shown to be incorrect once the number of individuals and sample area were factored out (Gray et al., 1997 Williams et al. (2018) also reported more of their GAB species had been recorded from SE than SW Australia; however, they acknowledge the lack of sampling at lower depths in the west may bias these results. Across the GAB, Williams et al. (2018) found that there was no longitudinal change in species composition, biomass or density at similar depths and considered it a single biogeographical province.
We found that the degree of similarity of the GAB samples to these off SE Australia varies considerably between taxonomic groups ( Figure 5). Most notably, none of the 68 recorded sponge species were collected from both GAB and eastern Australian surveys [although two of the eastern hexactinellid sponge species were found in the GAB using RoV collection devices, Hyalonema sp QM4976 and Lophophysema inflatum, the latter also known from elsewhere around the Indian Ocean (Tabachnick & Levi, 1999)]. At the other extreme, over half of the arthropod species from the LBA of the GAB were recorded from off the east coast ( Figure 5). At least some of the apparent regional endemism may be due to under-sampling (Coddington, Agnarsson, Miller, Kuntner, & Hormiga, 2009) of these inaccessible habitats.

| CON CLUS ION
The lower bathyal and abyssal fauna around southern and eastern Australia shows regional diversity and compositional differentiation at the scale of 1-2,000 km along both latitudinal and longitudinal gradients. This has management implications as mineral and petroleum exploration/extraction activities can occur over similar spatial scales, potentially spanning the entire range of some species. Some groups (such as sponges) show higher turnover than others (e.g., arthropods) and appear to be better indicators of regional endemicity at these depths.
Additional surveys outside areas of prospective resource development are required to better understand regional-scale patterns of biodiversity. This would support more informed evidence-based management of local developments. We were fortunate in this example covering 5,000 km of shoreline to have consistent access to expert biologists and taxonomists. In many instances, data for this kind of analysis will need to come from a multitude of academic and commercial surveys. Welldocumented best practice survey methods, appropriate metadata and open data will be required to characterize regional biogeography in these instances or to add to even larger scale analyses in the future.

ACK N OWLED G EM ENTS
The authors wish to thank the CSIRO Marine National Facility

DATA AVA I L A B I L I T Y S TAT E M E N T
All data and R scripts for the analyses used in this study are available in the Appendix S1.