Quantifying Patterns in Fish Assemblages and Habitat Use Along a Deep Submarine Canyon-Valley Feature Using a Remotely Operated Vehicle

Introduction

The deep sea is the largest environment on earth (Levin et al., 2001) and plays a pivotal role in cycling nutrients and water at global scales. It is also a major source of chemosynthetic primary productivity (Armstrong et al., 2012; Danovaro et al., 2014; Jobstvogt et al., 2014), and it plays a major role in climate regulation by absorbing heat from the atmosphere and sequestering carbon at the seafloor (Armstrong et al., 2012; Rogers, 2015). A number of commercially important industries such as fisheries, oil and gas, minerals and pharmaceuticals operate in deep-sea ecosystems (Armstrong et al., 2012). Deep-sea environments are characterised by a complex suit of geomorphic features such as underwater shoals, banks and canyons (Agapova et al., 1979; Heap and Harris, 2008) which provide diverse and structurally complex environments for marine organisms. The ocean is a major source of undocumented biodiversity, with experts predicting that between one to two-thirds of all marine species are undescribed, many of which reside in the deep sea (Appeltans et al., 2012).

Conducting research in deep-sea environments is particularly challenging due to a lack of natural light, high water pressure, low temperature, low oxygen levels and the need for large vessels and advanced technologies to reach these depths. Despite their remoteness, deep-water habitats are susceptible to numerous anthropogenic impacts on global and local scales, including climate change (Hoegh-Guldberg and Bruno, 2010; Rogers, 2015), overfishing (Clark, 2001), plastic pollution (Woodall et al., 2014), exploration and development activities such as drilling (Jones, 2009) or accidents such as oil spills (White et al., 2012). It is therefore important that quantitative data on these ecosystems is obtained in order to assess their vulnerability and ability to recover from disturbance, and to mitigate anthropogenic impacts on these ecosystems.

Much of the research on deep-water fish populations has utilised destructive sampling methods such as bottom trawling (Williams et al., 2001; Tolimieri and Anderson, 2010; Cruz-Acevedo et al., 2018). Trawl surveys are often completed across large spatial scales and work to aggregate samples, resulting in a decreased ability to provide fine-scale descriptions of organisms and their habitat associations (Cappo et al., 2004). A combination of improved technology and concerns over the use of destructive sampling techniques, has increased the use of non-destructive, camera-based surveys. Methods such as landers and baited remote underwater video (BRUV) are well suited to sampling deep-water fish assemblages given they are not limited by depth and use bait to attract fish to the camera system (Priede and Bagley, 2000; Zintzen et al., 2012; McLean et al., 2015). While remote video is a cost effective and statistically powerful way of sampling fish diversity across a gradient of habitats and depths, it is less suitable for finer scale sampling (<100 s of m) and only provides relative estimates of abundance due to variation in the distance a bait plume travels and the attraction of different fish species to the bait (Cappo et al., 2004; Watson et al., 2010; Galaiduk et al., 2017b). Technologies involved with the use of remote operated vehicles (ROVs) have also developed rapidly, allowing time efficient data collection at a fine-spatial resolution, across a wide range of depths (Jones, 2009; Sward et al., 2019; McLean et al., 2020). They also provide a platform on which multiple scientific instruments can be mounted, to facilitate multi-purpose surveys.

Observations made using ROVs provide contextual information about fishes that relate aspects of their population structure and function to habitat, which traditional trawl methods would otherwise overlook (Adams et al., 1995; Macreadie et al., 2018). ROVs have been used to study the impact of deep-water fishing (Puig et al., 2012; Bo et al., 2014), assess benthic associations of fishes with oil and gas structures (Hudson et al., 2005; Jones, 2009; Bond et al., 2018; McLean et al., 2018b; Schramm et al., 2020a, 2021), provide behavioural observations (Lorance and Trenkel, 2006; Gates et al., 2017), and to collect fragile specimens (Macreadie et al., 2018). Stereo-video, the use of two cameras to facilitate accurate estimates of the length of fish (Harvey et al., 2001) and to standardise a sampling area (Harvey et al., 2004), has also developed rapidly and it is now possible to attach a stereo-video system to an ROV and complete transect based fauna surveys (Schramm et al., 2020b). In addition, GPS position overlay can provide information on the precise location of observations. These complementary data streams are well suited to a spatial analysis framework such as species distribution modelling (SDMs), which allows the fitted models to be extrapolated into un-sampled locations using benthic environmental predictors derived from acoustic surveys. These models can help us to understand the ecology of these understudied species and map their environmental niche associations for future studies and management applications.

Although many studies have examined the biodiversity patterns of shallow water fishes, comparatively little is known about how these patterns change with increasing depth. It is common to observe a decrease in species richness of fishes with increasing depth (Moranta et al., 1998; Lorance et al., 2002; Tolimieri, 2007; Wellington et al., 2018), though this pattern can be reversed depending on the scale of depth examined (e.g., McClatchie et al., 1997; Mindel et al., 2016). Our study area was a canyon-valley feature located on the continental slope in the northwest shelf region of Western Australia. The study area traversed two Key Ecological Features (KEFs): (1) Canyons linking the Cuvier Abyssal Plain and the Cape Range Peninsula, and (2) Continental Slope Demersal Fish Communities (DSEWPaC, 2012). KEFs are elements of the Australian Commonwealth marine environment that are considered to be of regional importance for either a region’s biodiversity or its ecosystem function and integrity (DSEWPaC, 2012). Previous deep water research in Northwest Australia has identified high fish species and family level richness when compared to assemblages across similar depth ranges elsewhere in the world (Williams et al., 2001), and in a pattern similar to previous studies fish species richness declined with increasing depth from 78 to 825 m (McLean et al., 2018b). Our study set out to assess fish distribution and diversity, to better define these patterns and to add to our understanding of the deep-water ecology of the Northwest system, and the canyon and fish community KEFs.

To increase efficiency in the field, we used a multi-task stereo-ROV platform to survey fish assemblages in the northwest of Australia between the depths of 420–870 m. To the best of our knowledge, there have been no published studies that combine stereo-video technology with ROV transects to assess fish assemblages in deeper continental slope waters. Using a stereo-ROV we collected data on habitat, fish assemblages and water quality simultaneously. Our objective was to describe the abundance, composition and size of deep-water fish assemblages in three areas along a canyon-valley feature. We also fitted SDMs to document environmental niche associations of the four most abundant fish taxa and assessed the applicability of this approach to improve our understanding of the spatial ecology of deep-sea fishes.

Materials and Methods

Site Description

The study was conducted on the continental slope in the Greater Enfield region, located approximately 45 km from the North West Cape of Western Australia (Figure 1) at depths between 420 and 870 m. Three areas (A, B, C) ranging from deepest (A, 870 m maximum) to shallowest (C, 420 m minimum depth) were surveyed along a continuous submarine canyon-valley feature. The feature is not included in the National Submarine Canyon Database, and in comparison to those listed in the database is a relatively small tributary to the north of the Cape Range Canyon KEF. Area A was the deepest site (between 790 and 870 m) and had a wide (250–500 m) and topographically shallow (20–50 m) valley feature. The valley feature at Area B was narrower (600–250 m wide), topographically shallow (50–80 m) and traversed water depths of between 590 and 690 m. Area C was the shallowest survey site, in water depths ranging from 420 to 560 m, and where the canyon was more pronounced, being the narrowest of the three canyons and topographically the deepest (250–300 m wide, 200–250 m deep). The survey encompassed flat seabed and canyon-valley features as identified a priori from multi-beam bathymetry and derivatives. Following the criteria described in Huang et al. (2018) at Area C the feature is a shelf incising canyon, which is narrow, steep walled and deep. However, at areas B and C the feature is wider and topographically shallower, so at these areas it has transitioned into a valley on the continental slope (Huang et al., 2018). All three areas surveyed were predominantly soft sediment bottom, with sparse higher rugosity hard bottom features. For simplicity when describing statistical analysis, the term “canyon” is used to describe both the canyon and valley sections of the feature throughout the methods and results sections.

Survey Method

The study was conducted over three consecutive days from 31 October to 2 November 2015. The survey was conducted using a Centurion QX 312 work class ROV fitted with stereo-video cameras to conduct transects with a standardised width of 5 m. A single downward facing video camera was used to simultaneously record benthic habitat. Water quality data (conductivity/salinity, temperature, and depth; CTD) was collected with a Sea-Bird Electronics (SBE) 19plus V2 SeaCAT Profiler CTD, mounted to the ROV platform. The ROV operated at approximately 1.5 m above the seabed and at a forward speed of approximately 0.5 knots. The total transect length covered by the ROV at Area A was approximately 11,330 m, Area B 3,470 m, and C 5,990 m.

The stereo-video system was made up of two Sony HDR CX550 handycams within custom built underwater housings with a depth rating to 2,000 m. The stereo-cameras were mounted onto a rigid base bar with a separation of 700 mm between the cameras and were inwardly converged at an angle of 8°, following the principles outlined by Harvey and Shortis (1995) and Goetze et al. (2019). This configuration provided an overlapping field of view from approximately 0.5 m in front of the cameras and accurate length measurements out to 8 m (Harvey et al., 2010). Stereo-video footage was recorded at high definition with a 1,920 × 1,080 resolution. A second stereo-video system was also fitted as redundancy in case of failure of the primary system. The stereo-video systems were mounted as low as possible onto the ROV (Figure 2), and camera systems were interfaced with and powered from the ROV systems. The ROV communication channels allowed the systems to be remotely controlled, and a live standard definition video stream from each camera was fed to the surface with live position, depth, date, and time overlays. To maximise illumination of the field of view, a combination of high output Light Emitting Diode (LED) and High Intensity Displacement (HID) lighting was used. These were placed as high as possible above the camera system on the ROV to reduce backscatter from suspended particles in the water column.

Identification of Taxa, Three Dimensional Positioning and Length Measurements

Video footage was analysed using EventMeasure Stereo software (SeaGIS, 2014). Identifications of fish were made by one experienced researcher based upon morphological features. In some cases, identification to species level was not possible as identifying features could not be distinguished on the video footage. In such cases, a precautionary approach was taken and taxa were identified to the lowest taxonomic resolution possible. Distinct taxa were given a unique number, which facilitated the assessment of diversity.

The calibrated stereo-video system (see Boutros et al., 2015) allowed accurate and precise measurements of the fork length (tip of nose to the middle caudal fin rays) of fishes, and the identification of the three dimensional position of each fish in relation to the ROV. To ensure fish were within the transect area (a 5 m belt), length measurements further than 2.5 m to either side or greater than 7 m in front of the ROV (based on the minimum visibility), were automatically rejected by the EventMeasure software.

Data Presentation and Statistical Analysis

The collection of data over transects of known position and dimensions allowed the use of two statistical approaches: Firstly, a multivariate analysis of variance approach was used to investigate differences in the fish assemblage structure between the three areas, and between canyon and non-canyon habitats. Secondly, a species distribution modelling approach was used to investigate environmental drivers of the distribution of taxa that were characteristic of each area, and to identify one taxa that was common and widely distributed over all three areas. These models were used to generate continuous predictive maps of the occurrence of these taxa over the survey areas.

Approach 1, Assemblage Patterns

Definition of canyon and non-canyon sections

Using bathymetry and slope data and the known depth profile of the ROV track log, transects were separated into two habitat types: canyon and non-canyon. Canyon habitats were defined as benthic habitat with increased depth, slope and structure, while the flat seabed surrounding the feature was classified as non-canyon habitat. Where a transect crossed over a canyon-valley habitat, the depth was noted outside of the canyon on both sides and the maximum depth inside the canyon-valley feature was also recorded (depths labelled A, B, and C, respectively, on Supplementary Figure 1.1). Taking a conservative approach and to avoid areas of transition between the two habitat types, a buffer section between canyon-valley and non-canyon was excluded from analyses. This buffer was calculated as 10% of the difference between the maximum canyon-valley depth and the depths outside the canyon-valley feature (C - A and C - B on Supplementary Figure 1.1). Video footage from the remaining depth of the canyon-valley feature (orange line in Supplementary Figure 1.1) was used to determine fish assemblage structure inside canyon-valley habitats.

Within each surveyed area (A, B, C), every time the transect crossed a canyon-valley feature, the video footage within the canyon-valley habitat was treated as a single replicate. As such, Area A had six replicate transect sections within canyon-valley habitat (Supplementary Figure 1.2), while Areas B and C both had five replicate transect sections within canyon-valley habitat (Supplementary Figure 1.2). To ensure comparability between habitats, similar replication occurred for non-canyon habitats such that Area A had six replicate transect sections and Area B had five. However, much of the area surveyed at Area C was steeply sloping and only two sections could be classified as non-canyon (Supplementary Figure 1.2; Supplementary Table 1.1).

Standardisation of data for assemblage analysis

To ensure that transect sections were statistically comparable, transect length and densities of fishes underwent three levels of standardisation prior to statistical analyses: Firstly, sections of transect were removed when the seabed was not within the field of view (Supplementary Figure 1.2). Secondly, the length of transects between canyon-valley and non-canyon habitats were standardised to one another. To achieve this, the total length of transect within each habitat was determined within each area, and then averaged. A note was taken of the habitat type with the shortest average transect length at each area. This shorter average length was then used to randomly select replicate segments of transects of that length in the other habitat in that area. If a transect section was shorter than the average length, its original length was used (Supplementary Table 1.1). Lastly, measures of fish assemblages from each transect section (as defined above) were standardised to a 250 m length by 5 m wide transect, or 1,250 m². This transect length was used to increase the number of observations per transect for deep-sea fishes that are relatively sparse compared to shallower water fishes, but is similar to the length used other studies (Watson and Ormond, 1994; Newman et al., 1997; Westera et al., 2003). This allowed for fish assemblages within canyon-valley features of different sizes to be compared to each other and to non-canyon habitats. These standardised transects formed the replicates of the “canyon” and “non-canyon” levels within the factor “Habitats” in the statistical analysis design.

The fish assemblage abundance data and the total number of fish for each transect were standardised in this way. However, the number of taxa for each transect was not standardised to a 1,250 m² transect area. This was to avoid artificially inflating the number of taxa observed in any transect. The influence of varying transect length on the number of taxa was assessed by including transect length as a covariate in the initial analysis. The interactions between transect length and the factors of Area and Habitat were not significant (all p > 0.20). Therefore we considered that transect length was not a confounding variable for tests of the number of taxa, and was not included in further statistical analysis.

Assemblage analysis approach

To investigate whether there were differences in the fish assemblage composition between the areas and habitats, a two-factor statistical design was used. The two factors tested were Areas (fixed factor, orthogonal; Area A vs. Area B vs. Area C) and Habitats (fixed factor, orthogonal; canyon vs. non-canyon). The standardised fish assemblage composition data was analysed using PERMANOVA (non-parametric analysis of variance; Anderson et al., 2008) add-in to the PRIMER v6 statistical software (Clarke and Gorley, 2006). The data were investigated for homogeneity of variance for the factors Habitat and Area using PERMDISP (permutational analysis of multivariate dispersions, Anderson, 2006). No taxa were numerically dominant and data was not transformed as the data met the assumption of homogeneity of variance. The multivariate fish assemblage composition was tested using PERMANOVA based upon a Bray Curtis similarity matrix. Number of taxa and the total number of fish were analysed separately using a univariate PERMANOVA based on a Euclidean distance resemblance matrix (Anderson, 2006; Clarke et al., 2006a,b). Where a significant difference (p < 0.05) was detected for the factor Area, post hoc tests were conducted to determine which of the Areas were significantly different.

To illustrate the multivariate patterns in the fish assemblage, a Principle Coordinate Ordination (PCO) was presented with vectors overlaid (Anderson et al., 2008). The vectors illustrate the strength and direction of the spearman rank correlation between the density of the taxa and the PCO axes. Taxa with a correlation stronger than 0.4 to either PCO axes were plotted. Results of the univariate analyses (number of taxa and total number of fish) were presented using bar graphs of means and standard error (SE).

Approach 2, Species Distribution Models

We developed individual Species Distribution Models (SDMs) for four fish taxa that were observed most frequently in the ROV videos. These models can help understand the ecology of these understudied taxa and map their environmental niche associations for future studies. The SDMs were developed using boosted regression trees (BRT) and “gbm” package applied in R software (R Core Team, 2014). BRT is a machine learning algorithm for additive numerical optimisation of the loss function to iteratively increase the predictive performance of the final model while gradually emphasising poorly modelled observations in the existing collection of trees (Elith et al., 2006). In the last decade, it has gained popularity within the marine spatial community because it is insensitive to outliers, missing data, or monotone transformations, and can be easily used with any type of predictors such as numeric, binary, categorical (Pittman et al., 2007; Oyafuso et al., 2017; Stamoulis et al., 2018). Fitted values in the final model are computed as the sum of all trees multiplied by the learning rate and are much more stable and accurate than those from a single decision tree model (Elith et al., 2006).

The occurrence data for each of the four taxa was extracted from direct observations of each occurrence from the underwater video recording along the ROV tracks and were using in BRT fitting as records of taxa presence. As BRT models also require a sample of observations to characterise the available environment to discriminate used from available habitat (pseudo-absences, Boyce et al., 2002; Phillips et al., 2009; Franklin, 2010) we derived pseudo-absences for all taxa by randomly sampling all the available data points along the ROV tracks where the study taxa were not recorded. Because presence records for all study taxa were low, we created a final ratio of 1:1 of the observed presences and pseudo-absences of each study taxa along the tracks to effectively estimate unbiased parameters for rare populations (Fielding and Bell, 1997; Franklin, 2010; Galaiduk et al., 2017b). The final datasets were partitioned into training (75%) and testing (25%) data for individual modelled taxa and tested for spatial autocorrelation between observations. The explanatory variables were a set of 9 functionally relevant environmental predictors with Spearman’s rank correlation between them <0.7 which is considered to be acceptable for spatial models (Moore et al., 2011; Galaiduk et al., 2017a; Table 1) that describe the structure, complexity and type of benthic habitat derived from the bathymetry data using Spatial Analyst toolkit in ArcGIS 10.3. The 10th variable, “Habitat type,” was categorical, and described occurrence of benthic biota and signs of bioturbation. This was derived through unsupervised classification procedures in ERDAS using direct observation along the ROV transects and post processed with the Spatial Analyst toolkit in ArcGIS (Table 1).

To determine the effect of environmental predictors and their importance on the probability of occurrence of four fish taxa, we fitted BRT models on training datasets for these taxa following the procedures outlined in Elith et al. (2008). Optimal model settings were chosen using 10-fold cross-validation by optimising learning rate, bagging rate and tree complexity (Leathwick et al., 2006). The optimal model was considered a model that produced the lowest cross-validated residual deviance with at least 1,000 fitted trees. Selected models were then simplified to remove less informative predictors which achieved more parsimonious models without degradation of model fits (Elith et al., 2008). The importance of predictor variables in the simplified BRT models was determined using the variable importance score based on the improvements of all splits associated with a given variable across all trees when this variable was added in the model (Leathwick et al., 2006).

To assess the predictive performance, and discrimination and accuracy of fitted models, a set of common evaluation metrics of predictive performance was calculated on the test datasets. We used Receiver Operating Characteristic (ROC) and the area under the curve (AUC) as graphical means to test the sensitivity (true positive rate) and specificity (false positive rate) of a model output (Fielding and Bell, 1997). The AUC is a measure of the ability of a model to discriminate between a presence or absence observation (Elith et al., 2006). In addition, we calculated a threshold dependent Kappa statistic which is commonly used in ecological studies with presence-absence data and provides an index between 0 and 1 of how much a model predicted actual classes versus a guess (Cohen, 1960). A probability threshold that balances sensitivity and specificity was chosen as it provides a measure of how well the model predicts both presences and absences (Liu et al., 2005). After evaluation, the final models for individual taxa were predicted on 4 × 4 m grid using all available observations.

Results

Oceanographic Situation

The CTD profiles revealed a steady decline in water temperature with depth across all Areas, from ∼25°C at the surface to ∼6°C at the deepest point (∼800 m) (Figure 3A). Salinity increased from ∼35 ppt at the surface to a maximum of ∼35.5 ppt at ∼150–200 m, and then decreased to a minimum of ∼34.6 ppt at depths greater than ∼500 m across all three surveyed Areas (Figure 3A).

The relationship between temperature and salinity (T-S plot; Figure 3B) at the deepest site (Area A) was used to evaluate the localised water masses, with reference to recent evaluations of deep-water hydrography off the Gascoyne region (Woo and Pattiaratchi, 2008). Based on the analyses of Woo and Pattiaratchi (2008), three water masses were identified. Firstly, a lower salinity tropical surface water (TSW) with a temperature range of ∼22–25°C was found in the upper ∼100 m of water column. Secondly, a higher salinity South Indian Central Water (SICW) with a temperature range of 12–22°C was centred on ∼200 m depth. Lastly, a lower salinity Antarctic Intermediate Water (AAIW) with a temperature range of ∼6–9°C was identified ∼550 m depth.

Near-bottom time series profile from Area A (depths ∼800–875 m) showed a temperature range of ∼5.7–6.3°C and relatively uniform salinity (∼34.6; Figure 4). At Area B, near-bottom (∼560–680 m) temperatures ranged ∼6.6–7.8°C and salinity was relatively uniform (∼34.6; Figure 4). The shallowest transect (∼350–560 m) in Area C showed a temperature range of ∼6.7–13°C and a salinity range of ∼34.6–35.1 (Figure 4). In general, the coldest, most dense water was found as the ROV traversed the deepest portions of the canyon-valley features.

General Description of Fish Assemblages

Across the three areas surveyed, the total transect length was ∼20,790 m giving a total transect area of 103,950 m² or 10.4 hectares. Along the entire transect, a total of 610 individual fish were recorded belonging to 80 unique taxa and 41 families. A full list of the taxa recorded by area, the total number of each taxa and their mean lengths are shown in Supplementary Table 1.2. Many of the taxa and individual fish recorded were small bodied (<30 cm fork length; Supplementary Table 1.2), however, a large 63 cm morid cod (Moridae sp3), was recorded, along with a number of larger bodied elasmobranchs. A 1.3 m whaler shark (Carcharhinus sp1) was the largest elasmobranch identified. The western gulper shark (Centrophorus westraliensis) had the largest mean fork length (86 ± 2.5 cm SE). The indigo legskate (Sinobatis caerulea) was the second largest taxa recorded on average at 74 ± 11.3 cm disc length (excluding the tail). One large 71 cm Chimaera Hydrolagus lemures (blackfin ghost shark) was also measured (Supplementary Table 1.2).

Approach 1, Assemblage Analysis

For total number of fish, a significant difference was found (at α = 0.05) between Areas and between Habitats, but there was no interaction between the two factors (Table 2). On average, the greatest number of fish per transect (1,250 m²) was recorded at Area B (Figure 5; Area B vs. Area A t₍₁₈₎ = 2.56, p = 0.003, Area B vs. Area C t₍₁₂₎ = 2.20, p = 0.04). The total number of fish recorded per transect were similar at Areas A and C (Figure 5, t₍₁₄₎ = 1.59, p = 0.134). On average, a significantly greater number of fish were recorded in canyon-valley Habitats compared to non-canyon Habitats (Table 2; Figure 5). Canyon-valley Habitats were also found to support a significantly greater number of taxa on average than non-canyon Habitats (Table 2; Figure 5). No significant difference was found in number of taxa between the survey Areas (Table 2), although the composition of the assemblages did vary.

Significant differences in the multivariate fish assemblage composition were found between both Areas and Habitats, but no interaction between these factors (Table 2). Post hoc tests on the factor Area revealed that all Areas differed to one another (all t > 1.14, all p < 0.03). The multivariate PCO showed three distinct groupings (Figure 6). The PCO shows transects at Area A clustered to the left of the plot were correlated to unidentified anguilliform (eel-like) fishes, and Paraliparis sp1 (a snail fish), as both fish were abundant in this area. In comparison, Paraliparis sp1 made up only 2% of the individuals recorded at Area B and were absent at Area C (Supplementary Table 1.2). Samples to the right side of the PCO plot are arranged into two clusters (Figure 6). The bottom grouping is characterised by four Macrouridae (rat tail) taxa which were characteristic of Areas B and C. The grouping toward the upper right side of the plot is characterised by a mixed grouping of taxa. Chaunax sp1 (a coffin fish) were more abundant at both Areas B and C than at Area A. Ariosoma sp1 (a conger eel), Coelorinchus sp1 (a rat tail), and Moridae sp2 (a morid cod) were recorded only at Area B (Supplementary Table 1.2). One sample from Area C is grouped to the left of the plot with Area A samples. This transect contained a single fish, an unidentified anguilliform fish.

Approach 2, Species Distribution Models

Model Performance

Boosted regression trees models provided “good” model predictions (AUC = 0.80–0.89) for two taxa and “poor” predictions (AUC = 0.60–0.69; Table 3) for the other two taxa according to the criteria of Hosmer et al. (2013). Similar trends in the performance of fitted models were evident for all the other evaluation metrics with sensitivity, specificity and the total proportion of correct predictions being highest for the Macrouridae sp5 and Paraliparis sp1. These performance measures were further corroborated by Kappa statistics with models for Macrouridae sp2 and Macrouridae sp4 performing “fair” (Kappa = 0.21–0.40) and models developed for Macrouridae sp5 and Paraliparis sp1 providing “substantial” (Kappa = 0.61–0.80; Table 3) predictions (Landis and Koch, 1977).

Habitat Associations

Depth had the greatest influence on the probability of occurrence of all taxa (Figure 7). For Macrouridae sp2, Macrouridae sp5 and Paraliparis sp1, higher probability of occurrence was associated with water depth over 600 m. In contrast, Macrouridae sp4 was predicted to be most likely to occur at depths between 500 and 600 m (Figure 7). Various indices of topographic complexity of the relief (plan, profile, curvature and range 50) were also important in influencing the probability of occurrence of Macrouridae sp2, Macrouridae sp4 and Macrouridae sp5. Higher probability of occurrence of Macrouridae sp2 and Macrouridae sp4 was predicted on the seabed surrounding the canyon-valley slopes, whereas Macrouridae sp5 was predicted to be deeper and associated with both flat areas of low structural complexity, and within the feature on canyon-valley slopes and the canyon floor (Figures 7, 8). A higher probability of occurrence of Paraliparis sp1 was predicted on deep non-canyon habitats (Figures 7, 8).

Spatial Patterns of Distributions

A high probability of occurrence of Macrouridae sp2, Paraliparis sp1 and Macrouridae sp5 was predicted in the deepest area A (Figure 8; Supplementary Material 2, 3, 5). A high probability of occurrence of Macrouridae sp5 was predicted in Area A, and also over the deeper areas (>500 m) and within the canyon-valley of Area B. In contrast, higher probability of occurrence of Macrouridae sp2 and Macrouridae sp4 in area B was predicted for regions outside the canyon-valley and on its slopes, respectively (Figure 8). Macrouridae sp4 was also predicted to occur within the deeper canyon-valley of Area C and the shallower regions of Area A, reflecting its predicted relationship with depth (Figure 8; Supplementary Material 4).

Discussion

Deep-water environments and subsea canyons support diverse fish assemblages which are different to nearby shallow water assemblages (Williams et al., 2001; Fernandez-Arcaya et al., 2017). Here, we combined novel survey methods and advanced analytical approaches to document patterns of deep-water fish assemblage composition and spatial distributions. We also showed that subsea canyon environments support fish assemblages characterised by a higher number of taxa and overall number of individuals when compared to the areas outside the canyon, which likely reflects the heterogeneity of habitat and oceanographic conditions within the canyon (Klinck, 1996; McClain and Barry, 2010). A combination of non-extractive survey methods and SDMs provided a useful tool for improving our understanding of fine-scale spatial distribution patterns of deep-water fishes and the associated environmental drivers.

The survey area had a high diversity of taxa for a deep-sea environment with 80 taxa of fish from 41 families recorded, which is consistent with previous observations of fish diversity in this region (Williams et al., 1996, 2001). The family Macrouridae (rat tails) dominated the fish assemblage across the survey area, which is also supported by earlier research (Williams et al., 1996, 2001). We recorded 11 Macrouridae taxa whereas Williams et al. (1996) reported 16 species of Macrouridae using trawls within the same latitudinal and depth range (21–22°S, 300–900 m depth). There was a consistent pattern between areas in that the canyon-valley features supported a greater number of fish and a higher diversity of taxa than non-canyon habitat. This is likely due to variation in the habitat structure of the canyon-valley feature (McClain and Barry, 2010). At Areas A and B, the canyon-valley habitat was wide (250–500 m) and topographically shallow (20–50 m) in comparison to Area C which was a more typical canyon feature. Submarine canyon habitats can support higher abundances of benthic feeding fishes than nearby continental shelf slopes (De Leo et al., 2010). Greater prey availability in sediment infauna and epifauna may support the increased abundance of fishes, particularly Macrouridae, which were characteristic of canyon-valley habitats (De Leo et al., 2010).

A different fauna characterised the deeper Area A compared to Areas B and C. This was likely to be due to the influence of the AAIW deep water current at area A. Submarine canyons are important in linking continental slope and shelf areas and water masses. Canyon features change nutrient dynamics through increased upwelling and interactions with alongshore currents (Klinck, 1996), and are important for supplying nutrients to deep ocean through down-canyon transport (Fernandez-Arcaya et al., 2017). Area A was differentiated from Areas B and C by a higher abundance of the species group of unidentified anguilliform fishes and also a higher abundance of Paraliparis sp1 (a snail fish). Species from the genus Paraliparis generally inhabit the deep sea, and have been recorded in this region at depths of 1,030 m (Williams et al., 1996), and 2,821 m in New Zealand (Roberts et al., 2015). A similar pattern of a faunal break between fish assemblages on the upper and the mid-slope was reported by Last et al. (2011), which can be attributed to the influence of the AAIW (Williams et al., 2001). This deep-water current has a significant influence on the distribution of fish communities in the south-west region between the Great Australian Bight and north of the Ningaloo Reef where its core is known to fluctuate with depth from approximately 875 m at 27.50°S to 520 m around 21.50°S (Williams et al., 1996, 2001; Woo and Pattiaratchi, 2008). A full assessment of the influence of this current on fish assemblages within our study area would require further specifically targetted research.

Our study contributes additional information to understanding fish communities and their relationship to benthic habitats in two KEFs for the Australian northwest marine region; Canyons linking the Cuvier Abyssal Plain and the Cape Range Peninsula and Continental Slope Demersal Fish Communities (DSEWPaC, 2012). Deep-sea KEFs are often data-poor and their definition is often based on limited data (Falkner et al., 2009). The Canyons linking the Cuvier Abyssal Plain and the Cape Range Peninsula are representative features of the region but are not unique in a wider Australian context (Falkner et al., 2009). In terms of the Continental Slope Demersal Fish Communities KEF, Last et al. (2005) recorded 500 species making it the most biodiverse slope region in Australia. In comparison, we recorded 80 taxa, but our study focus was on a tributary to a canyon system, and not the entire continental slope assessed by Last et al. (2005).

We elucidated the distribution patterns of the four most abundant taxa within the assemblages using SDMs. SDMs are a useful tool for exploring habitat associations of the deep-sea fishes because they provide insights into the ecology of these rarely observed taxa and data can be extrapolated into unsurveyed areas. Often depth and habitat complexity drive environmental niche associations of demersal fishes (Monk et al., 2010; Moore et al., 2013; Galaiduk et al., 2018). This was observed in our study, where depth was the most important environmental predictor for all modelled taxa. Our models identified Paraliparis sp1 to be primarily associated with water depth over 800 m. In contrast, the three other taxa were predicted to occur at water depth greater than 500 m with Macrouridae sp4 predicted to mainly occupy depths between 500 and 600 m. Habitat complexity calculated at different scales was also an important predictor of habitat associations of the modelled taxa, suggesting that these taxa use various sized nooks and crevices (as estimated from the bathymetry derivatives) for shelter and protection from predation (Kelley et al., 2006). Macrouridae sp2 and Macrouridae sp4 were predicted to occupy the canyon slopes, whereas Macrouridae sp5 was predicted to be more associated with areas within the canyon and the canyon floor, suggesting within-family partitioning of the available environmental niche space to avoid resource overlap (Ross, 1986).

Biological data collection in deep-sea environments is expensive and labour-intensive. Large vessels capable of deploying multiple data collection platforms (e.g., ROVs, CTD probes, sediment corers, and camera systems) through the water column to the seabed, often to depths of 1,000’s of metres, are required (McLean et al., 2020). In this study, we streamlined the collection of multiple data types simultaneously (stereo-video imagery of fishes, benthic habitat imagery, and CTD) by deploying multiple instruments on a single ROV platform. This efficient configuration is recommended for future environmental studies of both shallow and deep seas as it minimises the time required to collect data. The methods are particularly applicable to baseline studies that support environmental impact assessments for offshore facilities (e.g., oil and gas) and inform the planning and management of marine reserves.

Remotely Operated Vehicle operated surveys are an industry-supported method with potential for large data collection in remote marine environments rarely accessible to researchers (McLean et al., 2017, 2018a; Macreadie et al., 2018). ROVs present a platform which can be fitted with multiple tools to collect a variety of environmental data from challenging remote environments (McLean et al., 2017, 2018a; Macreadie et al., 2018). Previously, sampling has relied on extractive trawls, or more recently deep-water baited cameras (Priede and Bagley, 2000; Jamieson et al., 2009; Marouchos et al., 2011; Wellington et al., 2021). These techniques are better suited to describing broader-scale spatial patterns (100 s of m) than the fine-scale sampling (<10 s of m) using the transect and SDM approach described in this study. The ROV transect approach allows for the in situ observation of fishes (Macreadie et al., 2018), and the investigation of fine-scale fish-habitat relationships through SDMs. It also eliminates the bias associated with the use of baited video systems (i.e., the bait plume effect; see Galaiduk et al. (2017b) and Monk et al. (2012) for a comparison between methods and applications in SDMs).

The ROV, however, has its own biases associated with noise, and the bright lighting required to collect quality video footage in dark environments which may affect the behaviour of mobile fauna (Ryer et al., 2009). We did not observe adverse behavioural reactions or fleeing, however, it is possible that certain taxa may simply avoid the ROV entirely and were thus not recorded (Linley et al., 2013). It is likely that avoidance by mobile fauna is greatest over soft sediment habitats such as those sampled during this study (Schramm et al., 2020a). These artefacts and avoidance may have implications for SDMs, however, they were consistent throughout the survey and should not affect statistical comparisons between Habitats or Areas. Traditional extractive methods for sampling deep-sea environments have biases and limitations relating to the gear and the spatial resolution of samples obtained. Previous studies from the region used a large netted paired-warp trawl rig (Williams et al., 1996, 2001) that may have under-sampled small species (Williams et al., 2001). Smaller meshed nets used in this region have caught a different suite of fishes, including more Ophidiidae (cusk eels) and Congridae (conger eels) (Williams et al., 2001). Many of the fishes measured during our ROV survey were between 10 and 30 cm in length, and so may have been under represented in previous trawl surveys. Similarly, in our stereo-video ROV survey the positive identification of many small bodied (<10 cm) individuals was difficult because morphological features were indistinct. Therefore, both stereo-video ROV and trawl methods may under sample the diversity of small bodied fishes.

Conclusion

The submarine canyon-valley feature sampled here supports a characteristically rich fish assemblage with a greater number of taxa and individuals than areas outside the feature. This canyon-valley feature is a tributary to the larger Cape Range canyon system which may play an ecological role in linking the Cuvier Abyssal Plain and the Cape Range Peninsula (Huang et al., 2018). Here, we present a novel application of ROV stereo-video transects in deep water, demonstrating the utility of this technique for facilitating fine-scale fish species distribution modelling that has the potential to feed into spatial management of deeper offshore waters.

Data Availability Statement

The datasets presented in this article are not readily available because the data are commercial in confidence. Requests to access the datasets should be directed to environment@bmtglobal.com.

Ethics Statement

Ethical review and approval was not required for the animal study because at the time that the study was conducted ethical review and approval was not required for observational studies.

Author Contributions

BS, EH, LT, MW, KI, and DM conceived and designed the study. BS, KI, and LT collected the data. BS, RG, and KI analysed the data. BS, RG, EM, JG, and KI developed the initial version of the manuscript. All authors contributed to interpreting the data and developing the full manuscript.

Funding

This research was supported by funding and data from the Woodside-operated Greater Enfield Project, a Joint Venture between Woodside Energy Ltd. and Mitsui E&P Australia Pty Ltd.

Conflict of Interest

KI, MW, and LT was employed by the company BMT Commercial Australia. DM was employed by the company Woodside Energy Ltd.

The remaining 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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

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

References