Placentia Bay is a deep-water embayment located on the southern side of the Island of Newfoundland, in the province of Newfoundland and Labrador (Atlantic Canada). Following an overview of its ecosystem status, present fish distribution and spawning locations, and a species inventory, Placentia Bay was designated by the Government of Canada Department of Fisheries and Oceans as an EBSA (Templeman, 2007; DFO, 2016). Any severe disturbance (i.e. oil spill) would be of greater ecological consequence within the EBSA boundary than if the same perturbation occurred beyond its extent (DFO, 2004). Placentia Bay hosts important capelin spawning beaches, seagrass meadows, seabird colonies, and herring aggregations (Sjare et al., 2003). It is also a location where charismatic megafauna such as whales (DFO, 2018) and leatherback turtles (DFO, 2012b) may be observed. Paired with its ecological relevance, Placentia Bay has local socio-economic importance, including traditional and commercial fisheries of groundfish and shellfish (Robichaud and Rose, 2006; DFO, 2012a; DFO, 2019), a growing finfish aquaculture industry (LGL, 2018), and frequent marine traffic, including oil tankers. Placentia Bay is considered to be amongst marine areas at greatest risk of oil spill across Canada (DFO, 2007).

Four survey areas along the west coast of Placentia Bay were selected for this study: Rushoon, D’Argent Bay, Burin and St-Lawrence (Figure 1). These sites are notable for their local ecological importance. Rushoon was the northernmost site surveyed, and is likely one of the more pristine of the four due to its remoteness. This may be subject to change, as development of salmon aquaculture progresses in the area, which may impact the local environment (LGL, 2018). Further southwest, D’Argent Bay was selected for its capelin spawning sites, herring aggregations, and high occurrence of whales (Sjare et al., 2003). Burin is near a relatively large cluster of towns, including the towns of Burin and Marystown. Finally, closest to the mouth of Placentia Bay, St. Lawrence was selected for its capelin spawning sites, as well as its role as an important seabird habitat (White, 2018).

Multibeam echosounder (MBES) surveying was carried out between winter 2018 and summer 2020 (Table 1) on board the 27’ Fisheries and Marine Institute research vessel D. Cartwright. MBES uses sound waves that travel from the sensor to the seafloor and back again, forming a swath of depth soundings across the path of the vessel. As sound waves echo back from the seafloor, the strength of their return is also measured, and this acoustic backscatter can be indicative of substrate types. Harder substrates like bedrock or boulders return a higher backscatter signal, while softer sediments like mud or sand will produce a weaker signal (Lamarche et al., 2011; Lurton and Lamarche, 2015; Calvert et al., 2015). Bathymetry and backscatter were collected using a Kongsberg EM 2040P, operating at 400 kHz. Positioning data were obtained using a Fugro 3,610 differential GPS with a Seastar subscription that provided spatial accuracy of up to 8 cm. Sound Velocity Profiles (SVP) between the transducer and seafloor were obtained using an AML BaseX sound velocity profiler, while sound velocity measurements at the MBES transducer head were obtained by an AML Micro SV sound speed sensor.

Raw sonar files for each survey area were imported into the Quality Positioning Services (QPS) Qimera v2.0.3 software. Bathymetric data were adjusted for Rushoon and D’Argent Bay using observations from the Argentia tide station (Station #835, https://www.tides.gc.ca/en/stations/835) and for Burin and St. Lawrence using the St. Lawrence tide station (Station #755, https://www.tides.gc.ca/en/stations/00755).). SVPs were imported into each respective Qimera project: 25 for D’Argent Bay, 19 for Rushoon, 19 for Burin, and 17 for St. Lawrence. Spline filters were employed to automatically remove outliers, and once processing was completed, a bathymetric surface was exported as a Floating Point GeoTIFF Grid at 5 × 5 m spatial resolution for each area. Processed files were additionally exported to. GSF format, and were imported into QPS Fledermaus Geocoder Toolbox (FMGT) v.7.8.4. For backscatter processing. Backscatter mosaics were exported at 5 × 5 m spatial resolution. All environmental rasters were projected to UTM Zone 21.

Bathymetric surfaces were used to derive terrain features that have been found to influence benthic biota distribution. These features can act as surrogates for variables that are difficult to measure directly (e.g. slope can act as a proxy for finer-scale currents). Six terrain features identified by Lecours et al., 2016a that capture a large amount of topographic information were calculated in ESRI ArcGIS using the Terrain Attribute Selection for Spatial Ecology toolbox (TASSE) (Lecours, 2015). These included slope (change in elevation), eastness and northness (orientation, calculated as the sine and cosine of slope), relative difference to the mean value (RDMV; relative position), and standard deviation (SD) of bathymetry (a measure of rugosity). Additionally, bathymetric position index (BPI) and vector ruggedness measure (VRM) were also selected to potentially provide information on seafloor structure. BPI is an adaptation of the topographic position index (Weiss, 2001) that measures the relative position of an area to the surrounding seabed (Lecours et al., 2016b). Positive and negative values represent peaks and troughs, respectively. VRM incorporates both slope and aspect in a single measure of surface roughness that is independent of slope (Hobson, 1972; Sappingtom et al., 2007; Martín-García et al., 2013), unlike SD. BPI, and VRM were derived from the bathymetric surfaces using the Benthic Terrain Modeller (BTM) Version 3.0 (Wright et al., 2012). Additionally, a layer measuring the distance from the coast was calculated within each survey area using the ‘Euclidean Distance’ tool in ESRI ArcGIS, which has previously been used to inform on benthic assemblage distributions (Degraer et al., 2008; Richmond and Stevens, 2014; Vassallo et al., 2018).

Features were derived from backscatter mosaics based on the spatial distribution of the varying shades of grey denoting backscatter intensity (Haralick et al., 1973). Three features were calculated using a grey-level co-occurrence matrix via the R package ‘GLCM’ (Zvoleff, 2020): contrast (local variation), homogeneity (closeness of distribution), and entropy (randomness) (Haralick et al., 1973). These features were selected for their common application in previous studies of similar systems (Blondel and Sichi, 2009; Samsudin and Hasan, 2017; Shang et al., 2021).

Terrain features were calculated across a range of spatial scales (i.e. using moving-windows of varying sizes), as no single scale can account for all ecological processes of a benthic ecosystem (Dolan, 2012; Lecours et al., 2015). All features (except for BPI) were derived using a 3x3-cell window, with additional scales obtained using the “calculate-average” approach (Dolan and Lucieer, 2014; Misiuk et al., 2021), by focal averaging using increasing window sizes (13 × 13-and 35 × 35-cell windows; Table 2). BPI layers were calculated by setting the outer and inner radii (in cells) of an annulus (ring-shape) and calculating the mean elevation value of all cells within the analysis neighborhood (Walbridge et al., 2018). BPI was calculated using the Benthic Terrain Modeler Toolbox, with outer and inner radii of 60 and 3, 150 and 45, 260 and 13, 630 and 195, 700 and 35, and 1750 and 525 cells.

Ground-truthing sites (n = 192) across all four survey areas were selected using a Generalized Random Tessellation Stratified (GRTS) survey design (Supplementary Table S1) (Stevens and Olsen, 2004), stratified by bathymetry and backscatter. For the St. Lawrence sites, MBES data were unavailable prior to ground-truthing, and bathymetry was interpolated from depth soundings on a nautical chart (e.g. Great St. Lawrence Harbour, Marine Chart CA4642_2). Ground-truthing consisted of underwater video with two systems: a Sony FDR-X3,000 Action Cam (1920 × 1,080 pixels, 60 frames/sec), contained in a Deep Blue Abysso waterproof housing, paired with two 3,500 lumen neutral white light Cree LED bulbs and two green lasers spaced 10 cm apart; and a Deep Trekker DTPod (1920 × 1,080 pixels, 30 frames/sec) with an integrated light and two red lasers spaced 2.5 cm apart. For both systems, a live video feed was used to adjust the height of camera above the seabed. Boat positioning at the start and end of each transect was obtained using the onboard Garmin GPS 16x; the midpoint of each transect was used to extract point-wise data from terrain features.

Two-minutes of continuous video were annotated using the Monterey Bay Aquarium Research Institute’s (MBARI) Video Annotation and Reference System (VARS) (Schlining and Stout, 2006). The presence or absence of Agarum clathratum and non-geniculate coralline algae (CCA) were recorded. Because of varying coverage across different ground-truthing sites, CCA were divided into three classes: full coverage, partial coverage, and absent. The full-coverage class was assigned when most of the seafloor was covered by CCA. Partial coverage was given when at least half of the substrate was covered. If little to no CCA was present, the site was marked as absent. Because abundance of S. droebachiensis may influence macroalgae distributions, counts at each site were included as a possible predictor variable.

All epifauna larger than 2 cm were counted and identified to the lowest possible taxonomic level, using expert knowledge and published species guides (Gosner, 1979; Harvey-Clark, 1997; Martinez and Martinez, 2003; Fox et al., 2014; Salvo et al., 2018). When species identification was not feasible, a morphotype approach was employed (Howell et al., 2019). The total number of individual species/morphotypes per transect was calculated and converted to densities by estimating the total area recorded for each transect. Using the Blender v.2.8.2, 12 frames at 10 s intervals were extracted and the distance between lasers was measured and used to calculate transect width with ImageJ. This was averaged and multiplied by the total length of the transect. Species/morphotypes that had an abundance of <5 individuals across all ground-truthing sites were removed to reduce the influence of low-abundance species. The dominant assemblages for ground-truthing sites with <5 faunal counts could not be characterized (De la Torriente et al., 2018); these sites were removed prior to analysis.

The species matrix was Hellinger transformed to reduce the importance of larger epifaunal abundances (Legendre and Gallagher, 2001; Borcard et al., 2011). Faunal assemblages were clustered using the average hierarchical clustering method ‘Unweighted Pair-Group Method using arithmetic Averages’ (UPGMA) (Sokal and Michener, 1958), with hierarchical relationships plotted as a dendrogram. The UPGMA approach is fast, simple, and may even outperform other clustering algorithms (Kreft and Jetz, 2010). Dissimilarity values between two branches of the derived dendrogram (i.e., “fusion levels”) were plotted to determine the optimal number of faunal assemblages. Silhouette widths were additionally calculated and plotted for each fusion level, providing a metric for distinguishing assemblages. Silhouette widths range from -1 to 1, with one representing assemblages that are clearly distinguished (Borcard et al., 2011). Once faunal assemblages were identified, the characteristic species or morphotypes of each assemblage were identified using an IndVal procedure (Legendre, 2013). A species accumulation curve was developed for each assemblage to identify if species richness was effectively captured (Ugland et al., 2003; Bevilacqua et al., 2018).

For both the assemblages and for each marine algae taxa model, a Boruta Feature Selection (Kurse and Rudnicki, 2010) algorithm was run separately to include terrain features grouped by scale (window of analysis: 3 × 3, 13 × 13, 35 × 35 cells) in order to reduce the number of candidate terrain features and promote model parsimony and support performance (Nemani et al., 2021). Important variables are identified by the Boruta wrapper as it compares the importance of a variable with a randomly shuffled version containing the same distribution of values (i.e., “shadow features”). A variable is deemed important if it consistently contributes more to the model than its shadow variable. Degenhardt et al., 2019 found that the Boruta algorithm generally outperformed other selection methods, and previous successful applications can be found in Li et al., 2016, Diesing and Thorsnes, (2017) and Nemani et al., 2021. Variables identified as “important” or “tentative” were selected for model training here. Variables were further dropped that had absolute correlation values exceeding 0.7, as determined using the R function ‘corrplot’ (Wei, 2013). In these cases, the variable with the lower impact on model accuracy was dropped.

Separate multiclass classification Random Forest (RF) and Light Gradient Boosting Machine (LightGBM) models were developed for the epifaunal dataset, the A. clathratum dataset, and the CCA dataset. Each model was trained with 2/3 of the samples, which included a proportional representation of classes (i.e. presence and absence of A. clathratum or epifaunal assemblage). The remaining 1/3 was reserved to test model performance. Model accuracy was assessed using the test data by computing a confusion matrix of predicted and observed classes (Congalton, 1991), from which overall accuracy and Cohen’s kappa were derived. Overall accuracy is the number of accurately predicted classifications divided by the total number of observations. The kappa statistic incorporates the chance of random agreement between observed and predicted classes based on their prevalence (Cohen, 1960; Allouche et al., 2006). Both the RF and LightGBM models were run with and without fine-scale (<1 m) substrate % coverage to understand the importance of fine-scale features in explaining spatial patterns. Using both the RF and LightGBM models and the environmental data layers, full-coverage predictive maps were developed for the distribution of epifaunal assemblages, presence/absence of A. clathratum, and for the absence/partial/full coverage of CCA. Since substrate % coverage features were not spatially continuous, full-coverage predictions were based on the terrain feature-only models.

RF is a well-established ensemble machine-learning algorithm that builds a ‘forest’ of classification trees from which predictions are obtained through majority voting. It employs ‘bagging’, where the data are repeatedly bootstrapped to train different classification trees (Quinlan, 1986), which are uncorrelated. Additionally, each tree split uses a random subset of variables, instead of the entire set. RF models were trained using the R package ‘randomForest’ (Liaw and Weiner, 2002), with the default number of variables included in each split (‘mtry’) retained for each model. Random Forest has been used successfully for predictive models with limited sample sizes similar to the number of observations in this study (Stephens and Diesing, 2014; Robert et al., 2015; McLaren et al., 2019; Misiuk et al., 2019; Ilich et al., 2021; Janowski et al., 2021).

LightGBM is a relatively recent Gradient Boosting Decision Tree algorithm (Ke et al., 2017). Unlike RF, LightGBM uses ‘boosting’ rather than ‘bagging’ to enhance performance. Boosting is another ensemble method in which classification trees are grown on the residuals from previous trees (Vaghela et al., 2009; Schapire and Freund, 2013). Tree growth and learning rate are controlled by a suite of tunable parameters (Table 3). The parameters recommended in the LightGBM release 3.2.1.99 guide (Zhang et al., 2012) were chosen to tune the model.

A variable importance plot was derived for each model, ranking variables based on their predictive importance. The RF model used Gini Importance, which measures the total decrease in node impurity by calculating the sum over the number of splits including a variables, across all trees (Friedman et al., 2001). LightGBM implements a gain-based method that is similar to the Gini Importance used by RF (Lundberg et al., 2020), but excludes unimportant features. Univariate partial dependence plots were derived from the best-performing model to visualize the relationship between an individual variable and a faunal assemblage or marine algae (Marini et al., 2015; Vassallo et al., 2018).

By comparing model accuracy and ranked variable importance, preliminary modelling using the derived terrain features suggested that important variables may be missing. Substrate heterogeneity information extracted from videos of the seafloor was tested to improve model performance. Substrate observations were obtained from the underwater video collected for biological analysis (described below). An image was extracted every 10 s from each video using Blender v.2.8.2 and images were imported to ImageJ (Image processing and analysis in Java). A 50-square grid was superimposed on each image and each square was labelled as one of seven classes: bedrock, boulder, gravel-mix, fine sediment (including both mud and sand), red algae (excluding coralligenous algae), Agarum clathratum, and Saccharina latissima (Figure 2). Following Connell et al., 2014, red algae can be qualitatively categorized as turfs (i.e. loosely to densely aggregated filamentous algae <15 cm tall and covering an area greater than 1 m²). Biogenic substrates were included because observing the sediment underneath was often not possible and algae presence is likely to influence the presence of other epifauna. The number of occurrences of each class was divided by the total number of grids for a ground-truth site, yielding fine-scale (<1 m) percent cover for each substrate class.

Of the 192 ground-truthing sites collected, 117 fulfilled the criteria for assemblage identification. One site (Rushoon-53) exhibited two distinct epifaunal assemblages and substrate types, and was split into two, producing 118 total samples. From these, a total of 12,096 individuals were counted and 55 taxa were identified (Supplementary Table S2), 14 of which (25.5%) were identified to species level. The brittle star morphotype ‘Ophiuroidea spp.’ had the highest overall abundance (6,693 individuals; 55.5% of all counts), followed by ‘Hormathia sp.1’ (1,631 individuals; 13.5%), Strongylocentrotus droebachiensis (766 individuals; 6.3%), and Echinarachnius parma (627 individuals; 5.2%). Porifera had the highest taxa diversity with 14 morphotypes. Echinodermata was second in taxa diversity (7 morphotypes; six species), followed by Cnidaria (8 morphotypes; two species).

Four epifaunal assemblages were identified (Figure 3; Table 4). Sixty-four of 119 sites (53.8%) were identified as OPH, typified by ‘Ophiuroidea spp.‘, ‘Porifera sp.5’, and ‘Hydrozoa sp.1’. OPH contained the highest taxa richness (51 taxa: 39 morphotypes, 12 species), as well as the most taxa unique to a single assemblage (8 taxa). SDR represented 21 sites (17.6%) and contained 13 morphotypes and seven species, with S. droebachiensis and ‘Stauromedusae sp.2’ as the typifying taxa. MIX was found at 26 sites (21.8%) and had the second highest taxa richness, with 30 morphotypes and 12 species. It was typified by ‘Hormathia sp.1’, ‘Cerianthidae sp.3’, ‘Cnidaria sp.1’, and ‘Sagittidae sp.1’. EPA was the rarest assemblage, with only eight representative sites (6.7%). EPA had the lowest taxa richness, with only eight morphotypes and seven species, and was typified by E. parma and Pseudopleuronectes americanus, with the former unique to this assemblage. While OPH exhibited the highest taxa richness, species accumulation curves indicate that the species/morphotypes richness of SDR, MIX, and EPA are underrepresented (Figure 4). Greater surveying effort could reveal additional species/morphotypes for each assemblage. However, this lack of plateau may also be the result of removing species with <5 occurrences across all ground-truthing sites.

Following the Boruta algorithm and assessment of collinearity, 18 variables were included in the epifaunal models, 18 in the A. clathratum models, and 24 in the CCA models. Without fine-scale substrate % features, the RF model accuracy was 61.0% (kappa = 0.31) and the LightGBM model accuracy was 68.3% (kappa = 0.47). With the addition of fine-scale substrate features, the epifaunal RF test accuracy increased to 78.1% (kappa = 0.62), and the LightGBM model accuracy increased to of 78.1% (kappa = 0.62). Both the RF and LightGBM models overestimated the occurrence of OPH, and as a result, a large number of SDR and MIX observations were incorrectly classified as OPH. This is evidenced in the RF and LightGBM predictive maps, where OPH was predicted to cover the majority of each survey area, excluding Rushoon (Figure 5). RF, however, included a greater coverage of SDR along the coastal side of D’Argent Bay and of MIX along the interior of St. Lawrence. LightGBM underestimated the coverage of both SDR and MIX in D’Argent Bay and St. Lawrence, respectively, even though ground-truthing sites assigned as those assemblages were observed. RF correctly identified all EPA observations while LightGBM was unable to correctly predict any of the EPA sites.

Only ground-truthing sites within the observed depth range for each marine algae were included. For A. clathratum model development, 115 sites at a depth range of 15–66 m were included. A. clathratum was present at 60 of 115 ground-truth sites, the majority of which occurred in Burin (25 of 60 sites; 41.7%). A. clathratum was present at over 50% of sites in all survey areas except for St. Lawrence. Model performance without fine-scale substrate % features was 82.1% (kappa = 0.64) for the RF model and 87.2% (kappa = 0.74) for the LightGBM model. Maps derived from these models predicted a majority coverage of A. clathratum across Burin, though a patchwork of absences is predicted by LightGBM in the interior of Burin’s extent. Differences in predicted A. clathratum extent is apparent in Rushoon and St. Lawrence. The greatest discrepancies were observed in D’Argent Bay though, with very little coverage predicted by LightGBM compared to RF (Figure 6).

132 sites were included in the CCA model at a depth range of 15–78 m. For CCA, 42 sites had full coverage (31.8%) and 31 sites had partial coverage (23.5%). Rushoon had the highest number of ground-truthing sites with full coverage off CCA (15 of 27). Aside from one instance of full coverage, the northeastern coast of St. Lawrence was characterized by partial coverage of CCA. Few absences of CCA were observed in D’Argent Bay (7 out of 36 sites) and Burin (11 out of 45 sites). When models were developed without fine-scale substrate % coverage, the CCA RF model had an accuracy of 68.9% (kappa = 0.49) and the LightGBM model had an accuracy of 68.9% (kappa = 0.51). Predictive maps derived from these differed in the concentration of absences, with LightGBM predicting an underestimation across all four survey areas of both partial and full coverage (Figure 7).

The accuracy of all marine algae models was increased with the introduction of fine-scale substrate % coverage features. The A. clathratum RF and LightGBM model accuracies increased to 89.7% (kappa = 0.80) and 92.3% (kappa = 0.85), respectively and both the CCA RF and LightGBM accuracies increased to 82.2% (kappa = 0.72).

Depth was the most important variable for explaining the distribution of epifaunal assemblages, according to the RF model (Figure 8). OPH and MIX were associated with deeper portions of the surveyed areas, while SDR and EPA were generally observed at depths shallower than 50 m. Relationships between EPA and the five most important variables appeared opposite to those observed for the other assemblages (according to the best-preforming RF model) (Figure 9). EPA was the only assemblage that was more prevalent with increased coverage of fine sediment and flat terrain. Assemblages SDR, MIX, and OPH appeared more prevalent at areas of rugged and boulder terrain, according to models. EPA was most prevalent within 2 km of the coast, with decreasing occurrence with greater distance. Unlike SDR and MIX, OPH was most prominent at sites closer to shore, but also occurred at the furthest offshore sites that were sampled.

Fine-scale substrate % coverage was identified as being even more important than depth to explain the spatial patterns of both A. clathratum and CCA models (Figure 8). Specifically, the presence of even low % coverage of fine sediment influences the presence of both A. clathratum and CCA. The predicted presence of either marine algae was higher with increasing coverage of harder substrates, including both boulder (%) and gravel (%) coverage (Figure 10).

According to both RF and LightGBM models, bathymetry was the second most important variable influencing the distribution of A. clathratum, with the majority of occurrences constrained to depths of 15–48 m. Bathymetry appeared less important for predicting CCA distribution in the RF model, and was not included at all in the CCA LightGBM model. The model suggested that softer sediments were more influential than harder substrate in determining the presence of A. clathratum (Figure 10). For CCA models, partial and full coverage were associated with higher backscatter values, while sites with absences were characterized by lower backscatter (Figure 11).