# Hierarchical modelling of epibenthic communities on the Scotian Shelf and Gulf of Maine (Atlantic Canada) in support of conservation planning

Francisco Javier Murillo, Benjamin Weigel, Donald Clark, Ellen Kenchington

► To cite this version:

Francisco Javier Murillo, Benjamin Weigel, Donald Clark, Ellen Kenchington. Hierarchical modelling of epibenthic communities on the Scotian Shelf and Gulf of Maine (Atlantic Canada) in support of conservation planning. Canadian Journal of Fisheries and Aquatic Sciences, 2024, 81 (12), 10.1139/cjfas-2023-0326. hal-04780283

HAL Id: hal-04780283

https://hal.science/hal-04780283v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hierarchical modelling of epibenthic communities on the Scotian Shelf and Gulf of Maine (Atlantic Canada) in support of conservation planning

Francisco Javier Murillo<sup>a, ∗</sup>, Benjamin Weigel<sup>b, †</sup>, Donald Clark<sup>c, ‡</sup>, and Ellen Kenchington<sup>b, §</sup>

<sup>a</sup>Bedford Institute of Oceanography, Fisheries and Oceans Canada, Dartmouth, NS; <sup>b</sup>INRAE, EABX, Cestas, France; <sup>c</sup>St Andrews Biological Station, Fisheries and Oceans Canada, St Andrews, NB

Corresponding author: Francisco Javier Murillo (email: Javier.MurilloPerez@dfo-mpo.gc.ca)

## Abstract

Identification of ecologically significant units at different spatial scales is essential for management of biodiversity attributes. This case study illustrates a coupled methodological approach to delineate benthic habitats and associated assemblages at different spatial scales. Two complementary analyses were employed based on the occurrences of 99 epibenthic invertebrate species in Atlantic Canada, sampled during depth-stratified random trawl sets. To identify epibenthic assemblage types, isometric feature mapping and partitioning around medoids (ISOPAM) was used in combination with a joint species distribution model (JSDM), which also produced continuous distribution surfaces across the spatial domain. ISOPAM identified nine significantly spatially coherent assemblages (biotopes) at spatial scales of hundreds of meters to hundreds of kilometers, with associated diagnostic species. Those assemblages were closely matched by eight regions of common profile using the JSDM, indicating strong environmental influences on their distributions. These smaller-scale assemblages were grouped into three biophysical units. The results from the JSDM were validated with independent data showing good congruence. This indicates that the spatial units ascribed to identified assemblages are robust and suitable for marine spatial planning and monitoring.

**Key words:** community structure, spatial management, epibenthos, ISOPAM, JSDM

## Introduction

Accurate, spatially explicit, and ecologically-relevant maps of marine benthic invertebrate communities are essential to support national and international policies aimed at protecting, restoring, and sustainably managing biodiversity (Lecours 2017), under an ecosystem-based approach (Long et al. 2015). In particular, the 2030 targets agreed upon in the Kunming-Montreal Global Biodiversity Framework (KMGBF) under the Convention for Biological Diversity (CBD/COP/DEC/15/4; CBD 2022), and the UN 2030 Sustainable Development Goals (Goal 14 Life Below Water) (A/RES/70/1) exemplify the need for spatial decision-support tools to meet conservation, restoration, and sustainable-use targets.

Target 3 of the KMGBF calls for at least 30% of marine and coastal areas, especially areas of particular importance for biodiversity and ecosystem functions and services, to be effectively conserved and managed through establishing ecologically representative marine protection networks, by 2030. For the North Atlantic, this ambitious international commitment is reflected in national policies such as the EU 2030 Biodiversity Strategy (COM/2020/380 final), the US' 30 × 30 Marine Conservation Goal' (Executive Order 14008), and Canada's 2030 National Biodiversity Strategy

(https://www.canada.ca/en/environment-climate-change/services/biodiversity/national-biodiversity-strategy.html; DFO 2023).

The National Framework for Canada's Network of Marine Protected Areas (Government of Canada 2011) provides strategic direction for the design of a national network of marine protected areas (MPAs) that will be composed of a number of bioregional networks. In support of this commitment, Fisheries and Oceans Canada (DFO) is leading a process to create a MPA network for the Scotian Shelf-Bay of Fundy Bioregion (DFO 2012). One of the overarching objectives of that network is to protect representative examples of the different ecosystems and habitat types (DFO 2018).

While for many areas of the Scotian Shelf-Bay of Fundy Bioregion, faunal assemblages of benthic invertebrates have been characterized (e.g., from portions of Georges Bank (Thouzeau et al. 1991), Sambro Bank (Hawkes et al. 2019), the deep waters of the Gully MPA (Kenchington et al. 2014), the eastern Scotian Shelf (Rincón and Kenchington 2016), and the Bay of Fundy (Staniforth et al. 2023) amongst others), there has been no comprehensive analyses and mapping of such communities across the entirety of the region. Classification of seabed features using a benthic habitat template based on "Scope for Growth" and "Natural Disturbance" (Kostylev
and Hannah 2007) has been used as a proxy for benthic habitats (DFO 2012) in the network design, in the absence of detailed analyses of benthic community composition. This classification system conveys information on habitat suitability and vulnerability, but not on the benthic assemblages per se. More recently, O'Brien et al. (2022) described the demersal fish communities in eastern Canada, using agglomerative hierarchical UPGMA clustering of occurrence data derived from trawl sets, and then predicted the distribution of each assemblage using a random forest modelling approach. They included a limited set of invertebrates in their analyses but incomplete and poorly resolved identification of the epifauna meant that the analyses were not representative of the benthic megafauna. However, use of data from the trawl surveys did allow for mapping over large spatial scales of thousands of km and more, while the depth-stratified random sampling facilitated statistical analyses.

Different approaches have been used to identified marine bioregions worldwide (e.g., Bloomfield et al. 2018; Hill et al. 2020; Stephenson et al. 2023). Murillo et al. (2016) identified benthic assemblages from trawl surveys in the Flemish Cap region (NW Atlantic) using a recursive partitioning method, whereby the clusters obtained at one clustering level were subsequently divided to form clusters at the next level of the hierarchy. ISOPAM (isometric feature mapping and partitioning around medioids; Schmidtlein et al. 2010) was used to identify clusters of species with high fidelity to groups of sites. This approach optimizes clusters and cluster numbers for concentration of indicator species in groups, which is a useful property for designing monitoring programmes. Each trawl set was associated with one of 12 spatially coherent epibenthic megafaunal assemblages and mapped as discrete points representing individual records of each assemblage. However, ISOPAM does not allow for extrapolation to unsampled areas and does not directly incorporate environmental covariates in its solutions.

Recent advances in statistical community-level models (Harris 2015; Warton et al. 2015; Ovaskainen et al. 2017) have allowed modelling the distribution of species-rich ecological communities through the use of joint species distribution modelling (JSDM) approaches. The advantage of JSDMs is that species are not modelled independently of each other but with an underlying joint covariance structure that relates to abiotic (environment) and biotic (species co-occurrences) filtering when predicting species responses. Thus, JSDMs can be used to better identify the drivers that control biodiversity across environmental gradients, space, and time, in a similar way to ordination-based methods, with the advantage that they have a predictive nature, allowing for inference to unsampled spaces.

The aims of this study were to describe the composition and map the distribution of epibenthic communities in the Bay of Fundy and Gulf of Maine, and on the Scotian Shelf and adjacent continental slope. Comprehensive identification of the benthic invertebrates from bottom trawl catches, that followed a depth-stratified random sampling design, was undertaken. Spatially coherent ecological assemblages were identified through the use of different, yet complementary, statistical methods, namely ISOPAM and JSDM, and compared. The

resultant maps of benthic communities can be used for the delineation of representative benthic habitats in the nascent MPA network (DFO 2018), while the methodology has broad application to other regions where maps of benthic communities are needed for spatial planning and other purposes.

## Materials and methods

### Study area

The study area (Fig. 1) encompasses the Bay of Fundy, parts of the Gulf of Maine, Georges Bank and the Laurentian Channel, plus the entire Scotian Shelf and a portion of the adjacent continental slope, and is delimited to the southeast by the boundary of Canada's Exclusive Economic Zone, and seaward to the 1500 m isobath. Within that area, a 5 km buffer drawn around all coastline was excluded from the analyses to avoid the inclusion of land points.

The Scotian Shelf and Scotian Slope comprise the continental margin south of Nova Scotia in Atlantic Canada. To the east, the Scotian Shelf is bounded by the Laurentian Channel that separates the Scotian Shelf from the continental shelf off Newfoundland, whereas the Fundian Channel separates it from the Gulf of Maine in the west. The Scotian Shelf extends more than 200 km offshore and is characterized by deep mid-shelf basins and shallow banks in the east and outer portions, with channels and submarine canyons along the slope (Fig. 1). Water mass characteristics, such as temperature, oxygen concentrations, or organic matter supply, in combination with circulation patterns can drive biodiversity and community patterns of benthic fauna (e.g., Neumann et al. 2013; Henry et al. 2014; Yasuhara and Danovaro 2014). The inner shelf is dominated by the outflow from the Gulf of St. Lawrence that travels through Cabot Strait, along the Nova Scotia coastline and around the Gulf of Maine. Mixing with offshore waters from the continental slope modifies the properties of shelf waters. These offshore waters are generally of two types, Warm Slope Water and the colder Labrador Slope Water (e.g., Loder et al. 1998; Hannah et al. 2001; Brickman et al. 2016).

Seafloor topography and geomorphologic features heavily influence benthic faunal composition and diversity (Levin and Sibuet 2012). The surficial geology of the Scotian Shelf is comprised of five formations that occur as an approximately 50 m-thick surficial succession over the bedrock (King 1980; King and Fader 1986). The Scotian Shelf Drift, a poorly sorted sediment of glacial origin, is the lowermost unit occurring sporadically on the Shelf. Between this layer and in some areas overlying it, is the Emerald Silt that is comprised of fine-grained, finely bedded muddy sediment. In the central parts of the basins, surface sediments are most often represented by LaHave Clay, a loosely compacted and homogeneous fine-grained silty clay, whereas on the banks, Sambro Sand and Sable Island Sand and Gravel dominate. These last units are well sorted sand with Sambro Sand containing around 10% silt and 5% clay, which differentiates it from the Sable Island Sand and Gravel. The upper Scotian Slope is overlain by fine sand in the east and sand and gravel in the west. Muddy-fine sand is found on the heads and floors of canyons to 1000 m depth, although outcrops of mudstone can be found on some
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Map showing the locations on the Scotian Shelf and in the Gulf of Maine/Bay of Fundy of the 2017 research vessel survey (NED2017020) fishing stations (solid black circles) used for analyses. The study area (dark black line) includes the DFO Maritimes Region limited by the 1500 m isobath and 5 km from the land. Named banks and basins are labelled. Base map created in ArcGIS (ESRI 2011) with map projection NAD83 UTM 20N. Administrative areas were obtained from GADM database (Global Administrative Areas 2011) and bathymetric data from the Canadian Hydrographic Service (CHS) Atlantic Bathymetric Compilation (ABC) dataset (Wang et al. 2022).</figcaption></figure>

canyon walls. Inter-canyon ridge crests shallower than 800 m depth and ridge crests up to 1500 m are covered by surface-winnowed sand (Piper and Campbell 2002). An updated surficial geology map, based on all available studies that have interpreted and mapped the surficial geology of the seafloor within the Scotian Shelf Bioregion, has been provided by Philibert et al. (2022). This compilation produced a descriptive surficial geology map with consistent map units across the Bioregion which was used to facilitate interpretation and characterization of the benthic assemblages.

## Data sources

### Biological data

Data used in this study were obtained from the 2017 Research Vessel (RV) Summer Ecosystem Survey (NED2017020),

carried out by the DFO on the Scotian Shelf, Scotian Slope, and Bay of Fundy between 26 and 1384 m depth (Fig. 1). The RV survey was conducted on the CCGS *Alfred Needler* between 28 June and 14 August. Fishing stations were allocated using a stratified random sampling design (Halliday and Kohler 1971) and conducted with standardized 30 min bottom tows at vessel speed of approximately 3.5 knots using a Western II-A bottom trawl gear, which has a 19 mm mesh lining in the codend (Tremblay et al. 2007) and swept area of 40 400 m². Fieldwork was carried out under DFO Section 52 fm(G)R License and research involving collection of specimens was conducted in accordance with all applicable laws, guidelines, and regulations.

Thirteen tows of less than 20 min in duration, or where the gear did not perform correctly, were excluded. Additionally, four stations fell inside the 5 km buffer drawn around all land; therefore, 240 stations were used in the analyses (Fig. 1).
<figure><img src="image_3.png" /></figure>

All epibenthic invertebrate fauna retained by the net were sorted on board to the lowest possible taxonomic level. A photographic catalogue of all species caught was created and voucher specimens for subsequent definitive identification in the laboratory were fixed in 70% ethanol or in 4% buffered formaldehyde depending on the taxon (e.g., Williams and van Syoc 2007). Preliminary identifications of all the invertebrate catch and biomass distribution of each phylum were provided in Murillo et al. (2018a). In total, 317 benthic taxa, belonging to 11 phyla, were identified in the study area (Table S1). Of these, 278 (88%) were identified to species level or were recognized as distinct species at a higher taxonomic rank and considered putative species (independent operational taxonomic units).

## Environmental data

Seven covariates, i.e., depth, slope, bottom temperature, bottom salinity, bottom current speed, sediment grain size, and fishing effort, were used in the analyses. These variables were chosen based on quality and availability of data, and assumed relevance to the community assembly processes (e.g., Murillo et al. 2020). Depth was obtained from the Canadian Hydrographic Service Atlantic Bathymetry Compilation and was available for the entire spatial extent with a native resolution cell size of 0.25 arc min or 463 m by approximately 350 m (≈0.16 km²) at the latitudes of our study area. Slope was derived from the depth raster projected in NAD 1983 UTM Zone 20°N coordinate system using the “Slope” tool in ArcMap’s Spatial Analyst toolbox calculated over a 3 × 3 grid cell window. The slope is defined as the maximum change in elevation over the distance between the cell and its eight neighbors and identifies the steepest downhill descent from the cell. Temperature, salinity, and current speed were extracted from a North Atlantic model termed BNAM, which has been shown to accurately reflect regional hydrography (Wang et al. 2018) and captures processes on the shelf and slope which are not discerned with global-scale models (Wang et al. 2021, 2022). BNAM has a horizontal resolution of 1/12°, equating to ~6 km on the Scotian Shelf. Data from the period 1990 to 2015 were extracted and the mean values between all years obtained. These values were then spatially interpolated across the study area using ordinary kriging in ArcMap 10.2.2 software (ESRI 2011) to create continuous surfaces with an ~1 km grid size (see Beazley et al. 2018 for more details). Grain size distribution for surficial sediments was based on Natural Resources Canada database of 10 000 sediment samples taken on the Scotian Shelf (interpolation details in Furlong and Kostylev 2008). Fishing effort from mobile bottom-tending gear for the period 2005–2014 was obtained from the sum of the fishing vessel monitoring system ping time (in hours) within each 1 × 1 km grid cell (see Koen-Alonso et al. 2018 for more details). All layers were posterily up-scaled to a 3 km resolution using the “Aggregate” tool in ArcMap with “mean” used as the aggregation technique. The point value for each station (using start position) was obtained from each layer and Spearman’s rank correlation coefficients (ρ) between predictor variables were calculated (Fig.

S1). Additionally, the surficial geology units from Philibert et al. (2022) were extracted for each station to provide a post hoc characterization of the bottom substrate in each faunal assemblage.

## Data analyses

### Community analyses

Only individual species or putative species were used for the numerical classification to avoid the inclusion of different taxa under the same variable and to aid in evaluation of the model outputs. The resultant species matrix was further reduced prior to analysis, considering only those species that occurred in greater than 3% of all stations, i.e., omitting very rare species. Consequently, the analyses were based on 99 species. The rest of the taxa (179) included species not well sampled by the gear (e.g., Calliotropis ottoi), rare taxa (e.g., Acanthocarpus alexandri), and others that although may be characteristic to their communities (e.g., Anthoptilum grandiflorum, Phormosoma placenta), were not common in our dataset due to limited sampling effort in their habitats (e.g., rocky areas that are avoided in the surveys or from deep stations >200 m that have a lower priority for sampling). After reducing the data, 9 phyla were represented in the analyses.

A cluster analysis of the sampling stations based on species composition was done using Sørensen's Index of Similarity and the isometric feature mapping (Isomap) and partitioning around medioids (ISOPAM) algorithm (Schmidtlein et al. 2010; Murillo et al. 2016). Isomap is an ordination technique that does not require any a priori assumptions with respect to the data structure. The ISOPAM algorithm searches for hierarchical classification of sites considering the overall fidelity of species to groups of sites. The optimality criterion derives from the G statistic (Sokal and Rohlf 1995), where G-values were computed for all species, adjusted using Williams' correction (Williams 1976), and standardized to correct for the influence of the number of clusters. The default standardized G-value (7) and two indicator species per cluster were defined as the stopping criteria for further partitioning of the data in the present study. Pairwise analyses of similarity (ANOSIM) of the resultant ISOPAM cluster groupings were implemented in the “vegan” package to test their statistical significance.

Ordered synoptic tables summarizing the frequency of species in groups along with their significance (assessed at α = 0.01) based on the *phi* coefficient (Φ) of their association to clusters (Tichý and Chytrý 2006) were produced. Φ measures the degree of association of one species to one cluster (fidelity of the species to the cluster) and ranges from 1 to -1 with an associated *P*-value. A threshold value of Φ = 0.34, based on the numbers of clusters, stations, and species (Schmidtlein et al. 2022), was used to distinguish species with high (≥0.34) and low (<0.34) fidelity to the cluster group. The analyses were done using the “isopam” (Schmidtlein et al. 2022) package in R. Each of the cluster groups generated using the ISOPAM algorithm was evaluated using sharpness and uniqueness indices (Chytrý and Tichý 2003). All species exceeding Φ > 0.05
<figure><img src="image_4.png" /></figure>

were considered to impart increased robustness to the calcu-
lations (Chytrý and Tichý 2003). Sharpness is defined as the
number or quality of diagnostic species in a cluster group rel-
ative to its average number of species. A group is considered
“sharp” if most of its species are confined to it, and are mostly
absent or rare in other groups. A cluster group is “unique”
if none of its diagnostic species are considered diagnostic in
other cluster group.

The relationship of the groups obtained with ISOPAM to
environmental gradients was examined using both uncon-
strained and constrained ordination. The unconstrained or-
dination was performed to show an analysis that can be
compared with ISOPAM results as ISOPAM is not constrained
by the environmental variables. The constrained ordina-
tion was used to determine the environmental variables re-
tained in the JSDM. First, detrended correspondence anal-
ysis (DCA) was used to ordinate the data to reduce a se-
vere arch effect detected in preliminary analyses (Hill and
Gauch 1980; Legendre and Legendre 1998). In DCA, Axis 1
was rescaled into units of the average standard deviation of
species turnover (SD) and the length of this axis is thought
to reflect beta-diversity, or species turnover along the gradi-
ent (Legendre and Legendre 1998), where species appear and
disappear over a distance of about 4 SD units. The environ-
mental variables were then fitted to the DCA using the vector
fitting envfit function in the “vegan” package. P-values for the
correlation of each variable with overall benthic community
composition were calculated. Secondly, canonical correspon-
dence analysis (CCA; ter Braak 1986) was performed. In CCA,
the effects of environmental variables on species composition
were tested with permutation tests under the reduced model
with 999 permutations to identify the most important vari-
ables amongst the covariates. Variables retained after apply-
ing the permutation test were used in the JSDM. Both analy-
ses were calculated using the “vegan” package in R (Oksanen
et al. 2018).

Bioregionalization

Hierarchical Modelling of Species Communities (HMSC;
Ovaskainen et al. 2017; Tikhonov et al. 2019b) implemented
in the “Hmsc” R 3.0 package (Tikhonov et al. 2019a) was used
to fit a JSDM to the benthic data to produce continuous dis-
tribution surfaces of the communities.

The workflow of the HMSC analysis followed five steps (Tikhonov et al. 2019b). The response variable was the species presence-absence matrix including the 99 epibenthic taxa. As fixed effects, all variables retained in the CCA analysis were implemented: depth, bottom temperature, bottom salinity, bottom current speed, sediment grain size, and fishing effort. The correlation between predictors was lower than 0.6, except between depth and salinity that presented a correlation of 0.8 (Fig. S1). Prior to the analysis, depth, fishing effort, and grain size were logarithmically transformed to gain a more homogeneous distribution of the variables; all other environmental variables were used in their raw state. A second order term for each environmental covariate was included to allow species to have their optimum at interme-

diate values within their multidimensional niche space in ac-
cordance with niche theory (Austin 2002). The data were then
modelled using the Bernoulli distribution and the probit link
function. A spatially explicit random effect (that also models
co-occurrence among species) at the level of sampling station
was included, using a latent factor approach (Ovaskainen et
al. 2016, 2017). The random effect is calculated as the sum
over the latent factors with which we estimate the site load-
ings and species loadings, where site loadings can be inter-
preted as non-measured environmental covariates or “hidden
covariates” not included in the fixed effects, at the level of
sites, and the species loadings as how species respond to these
hidden covariates (for details see chapter 7 of Ovaskainen and
Abrego 2020). This random effect represents residual species
co-occurrence patterns and/or patterns explained by environ-
mental variables not included in the model, where for the
former, same species loading to the site loadings would re-
sult in positive co-occurrence patterns. The model was fitted,
sampling the posterior distribution of 4 Markov Chain Monte
Carlo (MCMC) chains, of which each was run for 37 500 iter-
ations. Of these, the first 12 500 iterations were discarded as
burn-in. The chains were thinned by 100 to yield 250 samples
per chain, resulting in a total of 1000 posterior samples. Af-
ter model fitting, the MCMC convergence was examined us-
ing the potential scale reduction factor (Gelman and Rubin
1992) for all model parameters. In the case of non-full con-
vergence of parameters across MCMCs, a sensitivity analysis
was performed to evaluate if the RCP grouping was sensitive
to this aspect of the model, i.e., if the posterior samples of
each individual chain produce similar or dissimilar RCPs. If
RCPs remained similar, it was assumed that non-full conver-
gence of chains did not affect the RCPs and hence did not
impact model results on this aspect. The explanatory power
of the model was evaluated by computing the species-specific
area under the curve (Pearce and Ferrier 2000) and the coef-
ficient of discrimination Tjur’s R² (Tjur 2009) for presence-
absence data. This coefficient is defined as the difference be-
tween the average model prediction for successes and fail-
ures (Tjur 2009). The overall explanatory power of the model
was summarized as the mean Tjur R² across species. Addi-
tionally, a two-fold cross-validation assigning the stations in
each fold randomly, was performed to assess the predictive
power of the model. Following Ovaskainen et al. (2017) the
explained variation was partitioned among the fixed and ran-
dom effects for each species. Finally, predictions were made
for each of the 99 taxa in the full study area which included
18 653 grid cells, and regions of common profile (RCPs) were
created using the posterior mean of the predicted occurrence
probabilities for each species.

The main focus of this modelling effort was to compare
community composition and distribution with the results
obtained from ISOPAM, and if consistent, to map the ben-
thic communities as a continuous surface across the spatial
domain using the joint distribution modelling framework
HMSC. HMSC community outputs were run with a range of
cluster options. Further, the clustering in ISOPAM was based
on 240 stations, while for HMSC 18 653 grid cells were used,
invoking time and other constraints associated with cluster-
ing large data sets (Béjar Alonso 2013).
<figure><img src="image_5.png" /></figure>

The number of clusters to generate from the HMSC-derived occurrence probabilities was set at 3, to match the first level of partitioning in ISOPAM, and in a separate analysis, the number of clusters was set to 9, denoting the third level of partitioning from ISOPAM. These two spatial scales appear to respond to different environmental drivers and match different levels of a hierarchical marine ecological classification system (HMECS) used by DFO to classify the structure and distribution of Canada's marine biota and habitats at multiple spatial scales (DFO 2016). Although ISOPAM and HMSC are two different approaches which ideally are expected to converge in their results, ISOPAM uses only the actual observations to create the cluster groupings, and was therefore considered a better representation of the observed faunal associations than the HMSC RCPs, which in turn are model based and also consider environmental conditions as well as the study design when predicting species occurrence probabilities (here explaining 31% of the variation in the data), hence, adding a level of uncertainty to account for, for example, bias in sampling. However, the final grouping from HMSC was constrained to the ISOPAM results, rather than vice versa, assuming that the grouping on empirically collected data reflects the observed community at the sampled sites best. To find the best classification from HMSC matching the one from the ISOPAM algorithm at the third level of partitioning, a contingency table between the groups obtained from the ISOPAM algorithm and from different numbers of groups (7–10) generated with clustering from HMSC results was created. Classification groups obtained by HMSC can differ depending on the clustering algorithm method used. Because our goal was to achieve a spatial map from HMSC that best matches the ISOPAM assemblages, both k-means and partitioning around medoids (k-medoids) clustering algorithms were examined. The partition with the highest percent of similarity (measured as % of correct stations classified) between the ISOPAM algorithm and HMSC classification was selected as the optimal classification. Once the best partitioning from the HMSC was selected, the RCPs were calculated and the average occurrence of each species in each RCP was calculated to identify the species that best characterize each group.

To evaluate the robustness of the RCPs defined, a suite of diagnostic species from each ISOPAM assemblage was selected and corresponding occurrence data extracted from the 2019–2020 RV summer surveys. Those data formed an independent data set to test the accuracy of the community predictions by evaluating whether the station fell into its predicted map unit. This selection included the species with highest Φ in each assemblage that could be identified at sea. New Φ values, for each species were calculated based on the predicted RCPs from HMSC using the R package “indicspecies” (De Cáceres et al. 2023). Figure S2 illustrates a schematic summary of the analysis done in this study.

It is expected that species with limited geographic range and environmental tolerances present better model fit than those with greater ranges (e.g., Guisan et al. 2013; Morán-Ordóñez et al. 2017). These species are likely more diagnostic of the communities and therefore will present higher fidelity values. To study this relationship, a bivariate plot of Tjur R²

**Fig. 2.** Cluster analysis of the fishing stations using the ISOPAM (isometric feature mapping and partitioning around medoids) algorithm based on species composition. Sørensen's Index of Similarity applied to presence-absence data was used to construct the dendrogram. The ISOPAM algorithm generated nine groups and three levels of partitioning.

<figure><img src="image_6.png" /></figure>

for each species obtained from HMSC and maximum fidelity from ISOPAM was constructed and the Pearson correlation coefficient (r) calculated.

# Results

## Community analyses

### Characteristics of faunal groups

#### Major faunal groups

From the reduced species matrix (N = 99), the ISOPAM analysis presented three major classification levels (Fig. 2) with five partitions, resulting in nine cluster groups—each with a strong degree of spatial cohesion (Fig. 3). All groups were significantly different from each other in the ANOSIM of pairwise combinations (*P* < 0.001). The first level of partitioning (Fig. 3A) identified three clearly differentiated major groups of sample stations: (I) Bay of Fundy-Western Scotian Shelf; (II) Central Scotian Shelf and Slope; and (III) Eastern Scotian Shelf.

The Bay of Fundy-Western Scotian Shelf major group (I) was composed of stations from the Bay of Fundy, Georges Bank, and Browns Bank on sand and gravel bottoms between 32 and 269 m depth located in warm average bottom water temperatures and included the stations with higher fishing effort (Table 1, Fig. S3). The American lobster (*Homarus americanus*) was the most diagnostic species with fidelity value of 0.60 (Fig. S4). This group had a high degree of sharpness and moderate uniqueness (Table 1), indicating that there are many diagnostic species with a high *phi* value in this cluster (Fig. S4). The Central Scotian Shelf and Slope major group (II) included stations from Emerald Basin, Western and Sable Island Banks, Northeast Channel, Jordan Basin, and the deepest locations from the slope including the Laurentian Chanel.
**Fig. 3.** (A, C) Maps showing the spatial distribution of the colour-coded ISOPAM (isometric feature mapping and partitioning around medoids)-generated cluster groups (see Fig. 2) at (A) the first, and (C) third level of partitioning. (B, D) Regions of common profile (RCPs) applying clustering (pam algorithm) to the species distributions predicted by Hierarchical Modelling of Species Communities (HMSC). In panel B, the number of clusters has been set to 3, and in panel D, to 9. Each colour corresponds to a different RCP. Differences in colours between the stations and the RCP indicate mismatches between the ISOPAM clusters and the HMSC RCPs. Base map created in ArcGIS (ESRI 2011) with map projection NAD83 UTM 20N. Administrative areas were obtained from GADM database (Global Administrative Areas 2011) and depth contours from the Canadian Hydrographic Service (CHS) Atlantic Bathymetric Compilation (ABC) dataset (Wang et al. 2022).

<figure><img src="image_7.png" /></figure>

This group included stations spanning the largest bathymetric range, down to greater depths (from 33 to 1375 m depth) and therefore with mixed environments, but most of them were located on muddy bottoms (Emerald Silt and LaHave Clay with some over Scotian Shelf Drift). Due to the mix of environmental conditions included in this major group, the shrimp *Pontophilus norvegicus* was the only species with statistically significant association to this group, with a fidelity value of $\Phi = 0.34$ (Fig. S4). This heterogeneity is also reflected in the low sharpness and uniqueness values for this group (Table 1). In contrast, the Eastern Scotian Shelf major group (III) was composed of stations from the St. Anns, Banquereau, and Misaine Banks and was characterized by high sharpness and uniqueness (Table 1), higher than for major group I. This group included the colder and fresher stations between 36

and 312 m depth predominantly in sandy and gravelly bottoms (Table 1). Several species were highly associated with this group, with the snow crab (*Chionoecetes opilio*) and green sea urchin (*Strongylocentrotus droebachiensis*) the most characteristic species with fidelity values of 0.73 and 0.61, respectively (Fig. S4).

### Faunal assemblages within major groups

The second level of partitioning (Fig. 2) identified seven groups, two of which were further split at the third level of partitioning, resulting in nine final assemblages mostly associated with discrete areas (Fig. 3C). Diagnostic species for each assemblage can be consulted in Table S2.
<figure><img src="image_8.png" /></figure>

<table><caption>**Table 1.** Summary statistics (mean ± standard deviation) of variables for each cluster obtained at the first and third level of the ISOPAM (isometric feature mapping and partitioning around medoids) classification.</caption><thead><tr><th></th><th>N</th><th>S</th><th>U</th><th>Depth</th><th>Bottom temperature</th><th>Bottom salinity</th><th>Grain size</th><th>Bottom current speed</th><th>Fishing effort</th><th>Slope</th></tr></thead><tbody><tr><td>I</td><td>65</td><td>84.58</td><td>17.97</td><td>101 ± 42</td><td>5.96 ± 1.00</td><td>33.19 ± 0.97</td><td>1.75 ± 1.96</td><td>0.014 ± 0.007</td><td>11.22 ± 20.907</td><td>0.33 ± 0.41</td></tr><tr><td>II</td><td>111</td><td>58.29</td><td>13.20</td><td>256 ± 273</td><td>6.16 ± 1.19</td><td>34.21 ± 0.72</td><td>0.78 ± 2.64</td><td>0.015 ± 0.011</td><td>4.13 ± 10.50</td><td>0.89 ± 1.32</td></tr><tr><td>III</td><td>64</td><td>100.10</td><td>33.28</td><td>111 ± 57</td><td>3.56 ± 1.01</td><td>33.19 ± 0.68</td><td>1.70 ± 3.12</td><td>0.014 ± 0.007</td><td>0.09 ± 0.29</td><td>0.76 ± 0.78</td></tr><tr><td>I.a</td><td>41</td><td>65.80</td><td>1.98</td><td>94 ± 47</td><td>6.08 ± 1.02</td><td>33.05 ± 1.15</td><td>1.71 ± 2.21</td><td>0.013 ± 0.008</td><td>15.22 ± 25.17</td><td>0.22 ± 0.22</td></tr><tr><td>I.b</td><td>24</td><td>62.80</td><td>1.64</td><td>113 ± 30</td><td>5.75 ± 0.96</td><td>33.41 ± 0.52</td><td>1.83 ± 1.48</td><td>0.014 ± 0.005</td><td>4.38 ± 5.98</td><td>0.50 ± 0.58</td></tr><tr><td>II.a</td><td>69</td><td>51.22</td><td>0.93</td><td>196 ± 91</td><td>6.63 ± 1.00</td><td>34.34 ± 0.53</td><td>0.59 ± 1.00</td><td>0.014 ± 0.007</td><td>6.58 ± 12.73</td><td>0.82 ± 1.35</td></tr><tr><td>II.b</td><td>24</td><td>30.65</td><td>0.46</td><td>91 ± 46</td><td>5.88 ± 1.04</td><td>33.47 ± 0.84</td><td>1.83 ± 5.36</td><td>0.011 ± 0.009</td><td>0.03 ± 0.06</td><td>0.37 ± 0.54</td></tr><tr><td>II.c.1</td><td>8</td><td>51.69</td><td>1.08</td><td>344 ± 113</td><td>5.03 ± 0.71</td><td>34.47 ± 0.47</td><td>0.20 ± 0.27</td><td>0.013 ± 0.004</td><td>0.38 ± 0.57</td><td>0.60 ± 0.56</td></tr><tr><td>II.c.2</td><td>10</td><td>70.75</td><td>4.92</td><td>991 ± 357</td><td>4.47 ± 0.50</td><td>34.93 ± 0.02</td><td>0.03 ± 0.02</td><td>0.036 ± 0.016</td><td>0.13 ± 0.27</td><td>2.83 ± 1.25</td></tr><tr><td>III.a.1</td><td>18</td><td>61.11</td><td>0.73</td><td>154 ± 40</td><td>3.79 ± 1.03</td><td>33.56 ± 0.45</td><td>1.40 ± 3.10</td><td>0.013 ± 0.007</td><td>0.12 ± 0.38</td><td>1.24 ± 0.81</td></tr><tr><td>III.a.2</td><td>15</td><td>58.67</td><td>1.07</td><td>108 ± 41</td><td>2.61 ± 0.63</td><td>32.65 ± 0.81</td><td>3.59 ± 4.61</td><td>0.011 ± 0.005</td><td>0.02 ± 0.04</td><td>0.67 ± 0.78</td></tr><tr><td>III.b</td><td>31</td><td>59.79</td><td>0.98</td><td>89 ± 60</td><td>3.88 ± 0.87</td><td>33.24 ± 0.57</td><td>0.97 ± 1.62</td><td>0.015 ± 0.008</td><td>0.11 ± 0.30</td><td>0.52 ± 0.65</td></tr></tbody></table>

Note: Sharpness (S) and uniqueness (U) are also provided. N, number of fishing stations.

(1) *Georges Bank and inner Bay of Fundy* (I.a). Sampled between 32 and 269 m depth on sandy and gravelly bottoms and associated with the highest fishing effort (Table 1, Fig. S3). It is characterized by hydroids of the family Sertulariidae and erect bryozoans (Fig. S4). The horse mussel is also a characteristic species of this group ($\Phi$ = 0.43).

(2) *Browns Bank and outer Bay of Fundy* (I.b). Found between 74 and 192 m depth and characterized by demosponges of the orders Poecilosclerida and Polymastiida (Fig. S4).

(3) *Emerald and Jordan Basins* (II.a). Occurring in a mix of substrates between 86 and 564 m depth with the warmer bottom temperature and variable fishing effort (Table 1, Fig. S3). Characterized by the glass sponge *Vazella pourtalesii* ($\Phi$ = 0.41).

(4) *Emerald, Western, and Sable Island Banks* (II.b). Between 30 and 234 m depth on sandy bottoms (Table 1, Fig. S3) and does not present any typical species.

(5) *Laurentian Channel* (II.c.1). Found between 164 and 466 m depth on muddy bottoms and characterized by the sea star *Psilaster andromeda* ($\Phi$ = 0.70) and the sea urchin *Brisaster fragilis* ($\Phi$ = 0.61).

(6) *Scotian Slope* (II.c.2). Deepest assemblage located in muddy bottoms between 298 and 1348 m depth associated with higher bottom current speed and slope. It is the most unique group and with the highest association of species. It is characterized by the blind lobster *Stereomastis sculpta* ($\Phi$ = 0.89), the cup coral *Flabellum alabastrum* ($\Phi$ = 0.81), and the red crab *Chaceon quinquedens* ($\Phi$ = 0.77), amongst others.

(7) *Canso and deep Misaine Banks* (III.a.1). Occurring between 100 and 235 m depth mostly on sandy bottoms. It is characterized by the sea anemone *Hormathia digitata* ($\Phi$ = 0.56), the sea star *Ctenodiscus crispatus* ($\Phi$ = 0.50) which is also characteristic of the next group, and the soft coral *Duva florida* ($\Phi$ = 0.47).

(8) *St. Anns and shallow Misaine Banks* (III.a.2). Overlapping with the previous assemblage and occurring between 57 and

230 m depth. However, this group was associated with colder waters and most of the stations were on sand and gravel (Table 1, Fig. S3). It is typified by the basket star *Gorgonocephalus eucnemis* ($\Phi$ = 0.69), the stalked tunicate *Boltenia ovifera* ($\Phi$ = 0.63), and the soft coral *Gersemia rubiformis* ($\Phi$ = 0.63).

(9) *Banquereau Bank* (III.b). Similarly to the previous assemblage, this group was found on sand and gravel between 40 and 313 m depth. It is characterized by the hermit crab *Pagurus arcuatus* ($\Phi$ = 0.59), the sand dollar *Echinarachnius parma* ($\Phi$ = 0.56), and the sea cucumber *Cucumaria frondosa* ($\Phi$ = 0.53).

## Faunal patterns along environmental gradients

The nine assemblages created with the ISOPAM algorithm were separated in the DCA and CCA ordination of faunal stations (Fig. 4, Fig. S5). All environmental variables considered were significantly correlated with the DCA ordination axes (Table 2). Depth, bottom temperature, and bottom salinity showed the highest correlations, followed by fishing effort, sediment grain size, and slope, whereas bottom current speed explained the least amount of the variation ($R^2$ = 0.06). The main gradient (first axis) of the DCA (eigenvalue = 0.486) separated the Central Scotian Shelf and Slope major cluster group (II) from the Bay of Fundy-Western (I) and Eastern (III) Scotian Shelf, and it was strongly associated with depth, salinity, and sediment, the last one negatively. Whereas the second axis of the DCA (eigenvalue = 0.291) separated the Bay of Fundy-Western and Eastern groups and it was positively associated with temperature and fishing effort. Axis 1 length was nearly 7 SD units, indicating a high beta-diversity. Axis 2 also had a high beta-diversity with an axis length over 4 SD units.

Three of the groups (II.a, II.c.1, and II.c.2) of the Central Scotian Shelf were associated with positive values on the first axis, whereas the groups of the Bay of Fundy-Western and
**Fig. 4.** Detrended correspondence analysis ordination of the sampling stations based on species composition classified at the third level of the ISOPAM (isometric feature mapping and partitioning around medoids) analysis (Fig. 2). The arrows indicate the direction and magnitude of the relationship between the environmental variables and the ordination axis fit. Note: Temperature, salinity, and current are all from the bottom (seabed).

<figure><img src="image_9.png" /><figcaption>**Table 2.** Coefficient of determination (R²) and P-value of the environmental variables fitted to the detrended correspondence analysis (DCA).</figcaption></figure>

<table><thead><tr><th colspan="3">DCA</th><th colspan="2">CCA</th></tr><tr><th>Variable</th><th>R<sup>2</sup></th><th>P-value</th><th>AIC</th><th>Variable</th></tr></thead><tbody><tr><td>Depth</td><td>0.73</td><td>&lt; 0.001</td><td>1068.3</td><td>None</td></tr><tr><td>Temperature</td><td>0.65</td><td>&lt; 0.001</td><td>1055.6</td><td>Depth</td></tr><tr><td>Salinity</td><td>0.43</td><td>&lt; 0.001</td><td>1044.2</td><td>Depth + temperature</td></tr><tr><td>Fishing effort</td><td>0.28</td><td>&lt; 0.001</td><td>1041.3</td><td>Depth + temperature + salinity</td></tr><tr><td>Grain size</td><td>0.21</td><td>&lt; 0.001</td><td>1040.5</td><td>Depth + temperature + salinity + grain size</td></tr><tr><td>Slope</td><td>0.17</td><td>&lt; 0.001</td><td>1039.9</td><td>Depth + temperature + salinity + grain size + current speed</td></tr><tr><td>Current</td><td>0.06</td><td>0.003</td><td>1039.5</td><td>Depth + temperature + salinity + grain size + current speed + fishing effort</td></tr></tbody></table>

Note: Akaike information criteria (AIC) and environmental variables included in the canonical correspondence analysis (CCA) listed by decreasing AIC value. In bold model selected in the CCA. Temperature, salinity, and current speed are data from the bottom near the seabed.

Eastern Scotian Shelf and the Emerald, Western, and Sable Island Banks assemblage (II.b) were associated with negative values (Fig. 4). The former were well differentiated from each other, and from stations in groups I and III. Among the stations of the Bay of Fundy-Western and Eastern groups there was a clear differentiation, being the second ones associated with negative values of the second axis. However, the stations of the Browns Bank and outer Bay of Fundy assemblage (I.b) presented some overlap with stations mainly of the Emerald and Jordan Basins assemblage (II.a). Additionally, the DCA was not able to clearly separate stations from the Emerald, Western, and Sable Island Banks assemblage (II.b), which presented some overlap with stations mainly of the Bay of Fundy-Western Scotian Shelf (I) major group and with the Emerald and Jordan Basins assemblage (II.a).

The CCA only explained 16% of the variance and the best model included depth, bottom temperature, bottom salinity,

grain size, bottom current speed, and fishing effort, but adding the last three environmental variables only decreased AIC in less than one unit for each (Table 2). Therefore, the distribution of the epibenthic faunal assemblages on the Scotian Shelf and Slope is related mainly to depth, bottom temperature, and bottom salinity (Fig. 4).

## Bioregionalization

The HMSC model had an average explanatory power (Tjur R²) of 0.31, while the two fold cross-validation-based average predictive power was 0.18. These values are high relative to similar studies elsewhere (Murillo et al. 2020; Montanyès et al. 2023; Weigel et al. 2023). The MCMC convergence was satisfactory for the beta parameters (describing species niches), indicated by potential scale reduction factors (psrf) being on average 1.02 with a mean upper (95%) credible interval = 1.06. However, the alpha parameters (spatial scale of site
**Fig. 5.** Map of the first two spatial latent factors from the Hierarchical Modelling of Species Communities at the sampling station random level. The latent factors comprise spatially-structured site scores of model residuals, representing residual species co-occurrence patterns and/or patterns explained by environmental variables not included in the model. Base map created in ArcGIS (ESRI 2011) with map projection NAD83 UTM 20N. Administrative areas were obtained from GADM database (Global Administrative Areas 2011) and depth contours from the Canadian Hydrographic Service (CHS) Atlantic Bathymetric Compilation (ABC) dataset (Wang et al. 2022).

<figure><img src="image_10.png" /></figure>

loadings) did not fully converge indicated by a psrf of 1.66, well above the 1.1 suggested by Gelman et al. (2013) or 1.01 by Vehtari et al. (2021). The impact of non-converged alpha parameters on RCPs was assessed by a sensitivity analysis comparing the results of chain-specific RCP groupings to evaluate if different chains produced different groupings (Figs. S6 and S7). Since this was not the case, non-convergent alpha parameters did not impact site-specific species predictions and hence did not affect the RCP groupings. Differences between the RCPs generated with different grouping algorithms were higher than the differences observed between MCMC chains using the same algorithm. Bottom temperature explained most of that variability (25.8%), followed by depth (20.4%), and bottom salinity (13.9%). Sediment grain size, fishing effort, and bottom current speed explained 5.7%, 4.1%, and 2.9%, respectively. The spatial random effect at the sample level accounted for 27.1% of the variance, indicating that additional environmental covariates may have further resolved the RCPs. The first latent factor showed higher values around the Sable Island Bank and lower values in the outer Bay of Fundy, Jordan Basin, Fundian Channel, and Browns Bank (Fig. 5A). The second latent factor showed a clear pattern with negative values in the Western and Central Scotian Shelf and higher values as we move away from this area (Fig. 5B). The species with the best model fit, based on Tjur $R^2$ (Table S3), were the blind lobster *Stereomastis sculpta* (Tjur $R^2$ = 0.75), the American lobster *Homarus americanus* (Tjur $R^2$ = 0.70), the shrimp *Metacrangon agassizii* (Tjur $R^2$ = 0.67), the cup coral *Flabellum alabastrum* (Tjur $R^2$ = 0.63), and the snow crab *Chionoecetes opilio* (Tjur $R^2$ = 0.60). Three of these (*Stereomastis sculpta*, *Metacrangon agassizii*, and *Flabellum alabastrum*) were associated with the Scotian Slope assemblage (II.c.2), while the other two (the American lobster and snow crab) were

characteristic of the Bay of Fundy-Western Scotian Shelf and Eastern Scotian Shelf major groups, respectively.

## Group selection

The RCPs, resulting from a cluster with three groups (Fig. 3B), showed different results from ISOPAM at the first level of partitioning (Fig. 3A). The RCP groups were split based on depth and a longitudinal gradient that is related to temperature differences. One group was distributed along the Scotian Slope and Laurentian Channel grouping the deeper stations, whereas the other two groups were confined to the Eastern Scotian Shelf and to the Bay of Fundy-Western/Central Scotian Shelf.

The RCPs generated with a cluster of nine groups (Fig. 3D) were well aligned with the groups obtained with the ISOPAM algorithm at the third level of partitioning (Fig. 3C), except for the stations of the Emerald, Western, and Sable Island Banks assemblage (II.b) that were included in the Banquereau Bank assemblage (III.b). Additionally, a new group splitting the assemblage II.a into two groups was created (Fig. 3D).

The Emerald, Western, and Sable Island Banks assemblage (II.b) obtained with the ISOPAM algorithm did not find an equivalent RCP using the modeling approach. In addition, this group was not typified by any species and presented an overlap with other assemblages in the DCA and CCA ordination. When different grouping options were considered using both k-means and partitioning around medoids (k-medoids), the best solution aligned with the ISOPAM algorithm included eight groups (Table 3). Therefore, we concluded that eight is the optimal number of RCPs to map the benthic communities of the study area (Fig. 6), which corresponds
<table><caption>**Table 3.** Percentage of stations from each assemblage obtained at the third level of ISOPAM (isometric feature mapping and partitioning around medoids) classification correctly classified by Hierarchical Modelling of Species Communities (HMSC) for each of two clustering methods.</caption><thead><tr><th></th><th></th><th>I.a</th><th>I.b</th><th>II.a</th><th>II.b</th><th>II.c.1</th><th>II.c.2</th><th>III.a.1</th><th>III.a.2</th><th>III.b</th><th>Total</th><th>Comments</th></tr></thead><tbody><tr><td></td><td>k\N</td><td>41</td><td>24</td><td>69</td><td>24</td><td>8</td><td>10</td><td>18</td><td>15</td><td>31</td><td>240</td><td></td></tr><tr><td rowspan="4">PAM</td><td>7</td><td>73</td><td>92</td><td>77</td><td>0</td><td>88</td><td>90</td><td>0</td><td>93</td><td>84</td><td>67</td><td>No community II.b. III.a.1, and III.a.2 grouped together in HMSC</td></tr><tr><td>8</td><td>73</td><td>92</td><td>77</td><td>0</td><td>88</td><td>90</td><td>67</td><td>93</td><td>68</td><td>70</td><td>The new cluster aligns with III.a.1 and III.a.2</td></tr><tr><td>9</td><td>68</td><td>88</td><td>54</td><td>0</td><td>88</td><td>90</td><td>67</td><td>93</td><td>68</td><td>62</td><td>The new cluster is created by splitting the II.a into 2 communities</td></tr><tr><td>10</td><td>71</td><td>88</td><td>45</td><td>0</td><td>75</td><td>90</td><td>67</td><td>93</td><td>65</td><td>59</td><td>The new cluster group has no correspondence in ISOPAM</td></tr><tr><td rowspan="4">K means</td><td>7</td><td>71</td><td>79</td><td>58</td><td>0</td><td>88</td><td>90</td><td>0</td><td>93</td><td>84</td><td>60</td><td>Similar to PAM k7 but part of II.a included in I.b</td></tr><tr><td>8</td><td>71</td><td>79</td><td>58</td><td>0</td><td>88</td><td>90</td><td>61</td><td>93</td><td>71</td><td>63</td><td>The new cluster aligns with III.a.1 and III.a.2</td></tr><tr><td>9</td><td>71</td><td>88</td><td>52</td><td>0</td><td>88</td><td>90</td><td>61</td><td>93</td><td>71</td><td>62</td><td>The new cluster is created by splitting the I.b into 2 communities</td></tr><tr><td>10</td><td>71</td><td>83</td><td>54</td><td>54</td><td>88</td><td>90</td><td>67</td><td>93</td><td>65</td><td>67</td><td>The new cluster aligns with II.b</td></tr></tbody></table>

Note: N, Number of stations in the ISOPAM assemblage; k, number of groups defined in the groupings in HMSC; PAM, partition around medoids. The cluster mapping of each grouping can be consulted in Fig. S8.

strongly to the climatic variation observed in the ordination analysis (Fig. 4 and Fig. S5).

Table 4 provides a summary of the 10 most common species in each RCP based on the probability of occurrence and a full list of the probability of occurrences of all the species in each assemblage can be consulted in Table S3. The RCPs identified along the Bay of Fundy, Western and Central Scotian Shelf, Fundian and Laurentian Channels, and the slope (RCPs 1–5) were better characterized by the probability of occurrence of the characteristic species than those from the eastern Scotian Shelf (RCPs 6–8) (Fig. 7). Species characteristics of RCP 6 and 7 presented high probability of occurrence in both groups, although RCP 6 included higher probability of occurrence in the deeper groups (RCPs 4 and 5).

## Validation

When the RCPs obtained from HMSC were validated with data from the 2019 and 2020 summer surveys, i.e., with data not used in the model, a good congruence was found between the maximum fidelity values obtained by the ISOPAM algorithm and the fidelity values obtained using the groups defined by the RCPs (Table 5). All the selected species (except the mud star *Ctenodiscus crispatus*) presented maximum fidelity values in the same groups, indicating that the benthic communities defined are very robust.

## Relation between model performance and species fidelity

When the maximum fidelity of each species was matched with the model fit, a significant positive relationship ($\rho$ = 0.71, $P$ < 0.001) was found (Fig. 8). Species characteristics to the deeper assemblage (II.c.2) that presented fidelity values over 0.5 presented the best fit to the model. Interestingly,

some species with wider distribution that showed highest fidelity associated to high levels of the classification (e.g., *Homarus americanus*, $\Phi$ = 0.60, *Chionoecetes opilio*, $\Phi$ = 0.73 to the first level of partitioning, Fig. S4), presented also high model fit ($R^2$ = 0.70 and 0.60, respectively, Table S3).

## Discussion

In common with many nations, Fisheries and Oceans Canada (DFO) is committed to leading and facilitating integrated and ecosystem approaches to the management of Canada's three oceans (Oceans Act 1996). One important step in defining units for spatial management, and for informing ecosystem approaches to management, is the identification of ecologically significant units (Levin et al. 2009; Livingston et al. 2011). An HMECS, which can ensure that all types of habitat, communities, and ecosystems are effectively represented in marine spatial planning initiatives, and that biodiversity is considered at different spatial scales, has been proposed by DFO (2016). This study illustrates an approach to integrate biological and physical data to delineate benthic habitats at two different spatial scales, and describe their structure using multivariate analyses of biological community composition, and joint species distribution modeling incorporating physical habitat.

At the largest spatial scale, we identified “Biophysical units”, which are areas of distinct physiographic and oceanographic conditions and processes that shape species composition (DFO 2016). This level of classification corresponds to the Level 4 of the HMECS (DFO 2016). We then identified species assemblages at a smaller spatial scale corresponding to Level 6, referred to as “Biotopes” (Habitats and Communities). Biotopes are defined as “discrete taxonomic assemblages characterized by associated substrate and environmental factors” at the scale of hundreds of meters to
**Fig. 6.** Benthic communities predicted with the Hierarchical Modelling of Species Communities on the Scotian Shelf-Bay of Fundy and selected representative species. Base map created in ArcGIS (ESRI 2011) with map projection NAD83 UTM 20N. Administrative areas were obtained from GADM database (Global Administrative Areas 2011) and depth contours from the Canadian Hydrographic Service (CHS) Atlantic Bathymetric Compilation (ABC) dataset (Wang et al. 2022). Epifaunal taxa indicated in each benthic community photo mosaic, starting at the top left photo and following clockwise motion. (I.a) *Homarus americanus*, *Flustra foliacea*, *Placopecten magellanicus*, and *Modiolus modiolus*; (I.b) *Weberella bursa*, *Chirona hameri*, *Polymastia andrica*, and *Aldisa zetlandica*; (II.a) *Vazella pourtalesii*, *Actiniaria* sp.1, *Sclerasterias tanneri*, and *Epizoanthus* spp.; (II.c.1) *Brisaster fragilis*, *Psi-laster andromeda*, *Pennatula aculeata*, and *Gersemia fruticosa*; (II.c.2) *Metacrangon agassizii*, *Stereomastis sculpta*, *Flabellum* (*Ulocyathus*) *alabastrum* and *Chaceon quinquedens*; (III.a.1) *Ctenodiscus crispatus*, *Hormathia digitata*, *Terebratulina septentrionalis*, and *Solaster endeca/Solaster syrtensis*; (III.a.2) *Chionoecetes opilio*, *Gersemia rubiformis*, *Gorgonocephalus eucnemis*, and *Boltenia ovifera*; (III.b) *Cucumaria frondosa*, *Pagurus arcuatus*, *Echinarachnius parma*, and *Chlamys islandica*.

<figure><img src="image_11.png" /></figure>

hundreds of kilometers, and with a spatial resolution of hundreds of m² to 1 km².

Nine benthic invertebrate species assemblages or biotypes were identified using ISOPAM, each statistically differentiated from each other through ANOSIM. The optimal number of RCP from HMSC that best matched the ISOPAM results was eight. Depth, bottom temperature, and bottom salinity were the most influential explanatory variables, although sediment type, fishing effort, slope, and bottom current speed had some influence on the observed patterns, as has been observed in other studies at similar scales (e.g., Murillo et al. 2016).

Those nine assemblages, or biotopes, were organized into three large contiguous spatial biophysical units (Fig. 3A) separating the Bay of Fundy-Western Scotian Shelf from the Central Scotian Shelf and Slope and the Eastern Scotian Shelf. Depth and properties of the dominant water masses (bottom temperature and bottom salinity) were the main drivers of the observed faunal structure both in the ordination analyses (DCA, and CCA) and in HMSC. Two major groups of species with different temperature affinities were defined

around the 44° N or north and east of Sable Island. The western group (ISOPAM major groups I and II excluding the slope assemblage, and RCP 1 from HMSC, Fig. 3A and 3B), included many temperate species like the Jonah crab (*Can-cer borealis*), the American lobster (*Homarus americanus*) or the glass sponge (*Vazella pourtalesii*) that are common across the USA to the south, whereas the eastern group (ISOPAM major group III and RCP 2 from HMSC, Figs. 3A and 3B) included some cold-water species that may find their southern distribution around this biogeographic break (e.g., the bryozoans *Cystisella sacatta* and *Securiflustra securifrons*) or other boreal species that may be retracting in their southern distributions, and are only present in this area on the Scotian Shelf (e.g., *Argis dentata*, *Chionoecetes opilio*, *Drifa glomerata*, and *Hormathia digitata*). Stanley et al. (2018) also identified a genetic break on the Scotian Shelf along the 45° N aligned with a steep climatic gradient and driven by seasonal temperature minima. This biogeographic separation has been previously identified in other genetic (e.g., Van Wyngaarden et al. 2017) and bioregional studies (e.g., O'Brien et al. 2022). The surficial geology within each biotopes was heterogenous (Table 1, Fig. S3A).
<table><caption>**Table 4.** Ten most common species in each Hierarchical Modelling of Species Communities region of common profile (RCP).</caption><thead><tr><th>RCP 1 (~I.a)</th><th>RCP 2 (~I.b)</th><th>RCP 3 (~II.a)</th><th>RCP 4 (~II.c.1)</th></tr></thead><tbody><tr><td><strong>Homarus americanus</strong></td><td>Homarus americanus</td><td><strong>Cancer borealis</strong></td><td><strong>Pontophilus norvegicus</strong></td></tr><tr><td><strong>Asterias rubens</strong></td><td><strong>Polymastia uberrima</strong></td><td>Asterias rubens</td><td><strong>Pennatula aculeata</strong></td></tr><tr><td><strong>Cancer irroratus</strong></td><td>Asterias rubens</td><td>Homarus americanus</td><td>Ctenodiscus crispatus</td></tr><tr><td><strong>Hydrallmania falcata</strong></td><td>Cancer borealis</td><td>Pontophilus norvegicus</td><td><strong>Laetmonice filicornis</strong></td></tr><tr><td><strong>Placopecten magellanicus</strong></td><td>Ophiopholis aculeata</td><td><strong>Vazella pourtalesii</strong></td><td><strong>Psilaster andromeda</strong></td></tr><tr><td><strong>Symplectoscyphus tricuspidatus</strong></td><td>Weberella bursa</td><td>Laetmonice filicornis</td><td><strong>Lithodes maja</strong></td></tr><tr><td><strong>Eucratea loricata</strong></td><td>Mycale lingua</td><td>Pennatula aculeata</td><td><strong>Brisaster fragilis</strong></td></tr><tr><td><strong>Flustra foliacea</strong></td><td><strong>Hippasteria phrygiana</strong></td><td>Ophiopholis aculeata</td><td>Flabellum alabastrum</td></tr><tr><td>Cancer borealis</td><td>Placopecten magellanicus</td><td>Hippasteria phrygiana</td><td>Munidopsis curvirostra</td></tr><tr><td>Ophiopholis aculeata</td><td><strong>Tedania suctoria</strong></td><td>Bolocera tuediae</td><td><strong>Actinauge cristata</strong></td></tr><tr><td><strong>RCP 5 (~II.c.2)</strong></td><td><strong>RCP 6 (~III.a.1)</strong></td><td><strong>RCP 7 (~III.a.2)</strong></td><td><strong>RCP 8 (~III.b)</strong></td></tr><tr><td><strong>Flabellum alabastrum</strong></td><td>Chionoecetes opilio</td><td>Chionoecetes opilio</td><td><strong>Cucumaria frondosa</strong></td></tr><tr><td><strong>Stereomastis sculpta</strong></td><td><strong>Ctenodiscus crispatus</strong></td><td><strong>Strongylocentrotus droebachiensis</strong></td><td><strong>Echinarachnius parma</strong></td></tr><tr><td><strong>Metacrangon agassizii</strong></td><td>Argis dentata</td><td><strong>Crossaster papposus</strong></td><td>Asterias rubens</td></tr><tr><td><strong>Chaceon quinquedens</strong></td><td>Ophiopholis aculeata</td><td>Argis dentata</td><td>Cancer irroratus</td></tr><tr><td><strong>Munidopsis curvirostra</strong></td><td><strong>Mycale lingua</strong></td><td>Boltenia ovifera</td><td>Placopecten magellanicus</td></tr><tr><td>Pontophilus norvegicus</td><td><strong>Gersemia fruticosa</strong></td><td><strong>Solaster endeca/S. syrtensis</strong></td><td>Eucratea loricata</td></tr><tr><td><strong>Atlanlopandalus propinqvus</strong></td><td>Iophon cf. koltuni</td><td><strong>Ophiopholis aculeata</strong></td><td>Strongylocentrotus droebachiensis</td></tr><tr><td>Laetmonice filicornis</td><td>Solaster endeca/S. syrtensis</td><td><strong>Hyas alutaceus/H. coarctatus</strong></td><td>Cancer borealis</td></tr><tr><td>Pseudarchaster parelii</td><td>Crossaster papposus</td><td><strong>Gorgonocephalus eucnemis</strong></td><td>Ophiopholis aculeata</td></tr><tr><td>Chionoecetes opilio</td><td>Strongylocentrotus droebachiensis</td><td>Cucumaria frondosa</td><td>Pagurus arcuatus</td></tr></tbody></table>

Note: Boldfaced numbers represent highest probability of occurrence in RCP. The equivalent assemblage obtained at the third level of the ISOPAM classification is indicated in brackets.

This mix is likely associated to the length of the sampling station (~3 km) which could cross several geological types explaining the low importance of sediment grain size suggested by the HMSC and ordination analyses (DCA, and CCA), besides the importance that seafloor topography and geomorphologic play in the benthic faunal composition and diversity on the seafloor (e.g., Levin and Sibuet 2012).

At lower levels (biotopes), the ISOPAM analysis produced one assemblage in Sable Island Bank not typified by any diagnostic species (Fig. S4). This assemblage was included in the Banquereau Bank RCP 8 defined by HMSC results. O'Brien et al. (2022) included samples from that area in their "Banks/inner Bay of Fundy" assemblage. This area was also found to be similar to the Fundy region based on historical patterns of groundfish distribution (Mahon et al. 1984). Kostylev and Hannah (2007) identified Sable Island Bank as an area with high disturbance and low Scope for Growth due to the large bottom stress caused by waves, suggesting that this area would have relatively few species and all characterized as extreme adversity tolerators or *r* strategists. A closer examination of Table S2 showed that this area was characterized by species with wide distributions diagnostic of other assemblages (e.g., *Cucumaria frondosa*, *Placopecten magellanicus*) or with low fidelities and therefore not diagnostic (e.g., *Asterias rubens*, *Cancer irroratus*, and *Hippasteria phrygiana*). Sable Island Bank was also identified as one single spatial unit in the preliminary analysis of the biogeographic units on the Scotian Shelf conducted by Zwanenburg et al. (2010) where only environmental variables were used and the study area was divided into 12 groups. This area has one of the largest

seasonal changes in bottom temperature on the eastern Scotian Shelf (Han et al. 1999), which make this area of particular interest.

Further, the lower level ISOPAM communities from the slope may be under-represented, as only eight fishing stations deeper than 600 m were sampled during the survey and some of their characteristic species, such as the Peters's sea urchin *Hygrosoma petersii*, the long-nosed armoured shrimp *Glyphocrangon longirostris*, or the sea spider *Colossendeis colossea*, were excluded from the analysis due to their low occurrence in the overall study area (Table S1). Using a threshold depending of the depth range in which the species is caught (e.g., Thompson et al. 2022) could have retained some of these species. However, due to the low number of stations sampled in the deeper slope overall, these selections could not be representative of the whole slope area, and a greater number of stations are needed to adequately characterize this community. This reduction in sampling of the slope also explains the lack of coral records in our samples, despite known habitat and presence of coral species along the slope and in the canyons (e.g., Kenchington et al. 2014; Metaxas et al. 2019; Wang et al. 2022). Increased data collection along the slope and in the deeper bathyal and abyssal waters, ideally with non-destructive sampling tools such as ROVs, or drop cameras, are needed to characterize those regions.

The spatial patterns observed in this study are aligned with the benthic communities and complexes described by Nesis (1965) on the eastern Scotian Shelf and Scotian Slope more than 50 years ago, including the presence of cold waters species off Cape Breton still recorded today (e.g.,
<figure><img src="image_12.png" /><figcaption>**Fig. 7.** Probability of occurrence from Hierarchical Modelling of Species Communities (HMSC) for each species grouped by ISOPAM (isometric feature mapping and partitioning around medoids) assemblages (A) and regions of common profile (B). Species ID can be consulted in Table S2 (A) and Table S3 (B). Note that the 99 species differ in order between (A) and (B).</figcaption></figure>

*Ciliatocardium ciliatum*, *Solaster syrtensis*, and *Grogonocephalus eucnemis*. In addition, the observed spatial patterns at both scales or classification levels are similar to those from the demersal fish assemblages observed on the Scotian Shelf more than 30 years ago (Mahon et al. 1984; Mahon and Smith 1989). A study of the groundfish assemblages at four different decadal periods (Zwanenburg and Jaureguizar, unpublished, in NAFO 2010) showed four persistent areas that are aligned

with the benthic groups at the first level of the ISOPAM classification and with the two major groups observed on the shelf by the HMSC (Fig. 3B).

The groups observed were consistent with the ordination (DCA, and CCA) results, where a high species turnover (beta diversity) was observed. In our analysis, Axis 1 length was nearly 7 SD units and was highly correlated with depth, indicating a full change in species composition moving from
<table><caption>**Table 5.** Fidelity value (Φ) for each of the species from an independent data set (obtained from the 2019–2020 research vessel summer surveys) that were characteristic for each regions of common profile (RCPs) (top row) and for each ISOPAM (isometric feature mapping and partitioning around medoids) assemblage (lower row, from Table S2).</caption><thead><tr><th></th><th></th><th>RCP 1 I.a</th><th>RCP 2 I.b</th><th>RCP 3 II.a</th><th>RCP 4 II.c.1</th><th>RCP 5 II.c.2</th><th>RCP 6 III.a.1</th><th>RCP 7 III.a.2</th><th>RCP 8 III.b</th></tr></thead><tbody><tr><td rowspan="3">I.a</td><td>Homarus americanus</td><td>0.44</td><td>0.33</td><td>0.04</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Modiolus modiolus</td><td>0.24</td><td>0.34</td><td>0.23</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Placopecten magellanicus</td><td>0.39</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>I.b</td><td>Chirona hameri</td><td>0.06</td><td>0.13</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>II.a</td><td>Vazella pourtalesii</td><td></td><td></td><td>0.28</td><td></td><td>0.09</td><td></td><td></td><td></td></tr><tr><td rowspan="2">II.c.1</td><td>Psilaster andromeda</td><td></td><td></td><td>0.66</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Brisaster fragilis</td><td></td><td></td><td>0.70</td><td></td><td>0.19</td><td></td><td></td><td></td></tr><tr><td rowspan="2">II.c.2</td><td>Flabellum alabastrum</td><td></td><td></td><td>0.58</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Chaceon quinquedens</td><td></td><td></td><td>0.61</td><td></td><td>0.04</td><td></td><td></td><td></td></tr><tr><td rowspan="2">III.a.1</td><td>Clenodiscus crispatus</td><td></td><td></td><td>0.25</td><td>0.45</td><td></td><td></td><td></td><td></td></tr><tr><td>Gorgonocephalus eucnemis</td><td></td><td></td><td>0.17</td><td>0.81</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">III.a.2</td><td>Boltenia ovifera</td><td></td><td></td><td>0.04</td><td></td><td>0.24</td><td>0.51</td><td></td><td></td></tr><tr><td>Echinarchnius parma</td><td>0.00</td><td></td><td></td><td></td><td>0.21</td><td>0.69</td><td></td><td></td></tr><tr><td rowspan="2">III.b</td><td>Cucumaria frondosa</td><td>0.06</td><td></td><td></td><td></td><td>0.23</td><td>0.46</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.13</td><td>0.63</td><td>0.15</td><td></td></tr></tbody></table>

Note: Gray fill represents highest association of the species to the assemblage, while boldfaced numbers represent highest fidelity value. Only positive values are shown.

the shelf to the slope, which can also be observed in Table S2. Another complete turnover in species composition could be observed in Axis 2 where the length was also over 4 SD, likely associated with differences in bottom temperature. These results are very similar to the beta-diversity found in the Grand Bank and Flemish Cap area (Murillo et al. 2016), where the length of the two first DCA axes were 7.86 and 3.44 SD units respectively, indicating a complete turnover in species composition over a similar spatial scale of ~300–400 km at similar depths. Interestingly, in the Flemish Cap study and here, fishing effort was correlated with the second axis of the ordinations, indicating that fishing is either a direct (physical abrasion pressure) or indirect (correlation with temperature and hence potentially with target species in this study) determinant of species composition.

A consistency in spatial distribution between the benthic invertebrate and fish assemblages has been observed in other regions of the northwest Atlantic (Murillo et al. 2016; O'Brien et al. 2022) and elsewhere (Neumann et al. 2013), suggesting

a distribution pattern driven by similar responses to common environmental gradients rather than stochastic process or particular interactions between the species. Our HMSC results confirm these previous observations, as most of the explained variance was due to the environmental variables included in the model. Only 27% of the explained variance was associated with the spatial random effect, which could be explained by other assembly processes, such as biotic filtering resulting from interspecific and intraspecific competitive and facilitative interactions; filtering due to small-scale dispersal limitation; and/or filtering caused by environmental variables not included in the model (Montanyès et al. 2023). The first two spatial latent factors (Fig. 5) seem to be reflecting large scale patterns (over hundreds of km), so it is unlikely that they are related to species interactions.

Part of the variability observed in the first latent variable can be explained by the annual range in bottom temperature, where larger values are observed around Sable Island and Banquereau Banks, the areas with the highest positive
<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Bivariate plot relating the maximum *phi* coefficient (Φ) from ISOPAM (isometric feature mapping and partitioning around medoids) and the discriminatory power Tjur R² from Hierarchical Modelling of Species Communities. Species ID can be consulted in Table S2 (A).</figcaption></figure>

residual correlation. However, this is also an area of high frequency sediment mobilization by waves (Li et al. 2021) and was identified by Kostylev and Hannah (2007) as an area of high disturbance. The second latent variable may be indicating the uniqueness of this area in relation to incurrence of Gulf Stream eddies and changes to the strength and position of the Warm Slope Water (e.g., Seidov et al. 2021) that may have influenced the colonization processes.

Defining the number of ecologically significant units in biogeographical studies is a well-known problem in ecology as it depends on the classification method and criteria used (e.g., Stephenson et al. 2018). It is common to use different stopping rules (e.g., Rubidge et al. 2016) or to get an optimal number from different methods based on a majority rule (e.g., Weigel et al. 2023). Our approach of using the results of the ISOPAM algorithm to find the optimal number of clusters for applying HMSC seems to be robust based on the validation results. The ISOPAM algorithm (Schmidtlein et al. 2010) selects a classification which maximizes the overall fidelity of species to groups of sites providing a set of diagnostic or indicator species that can be used for future monitoring (e.g., Murillo et al. 2018b). It is able to cope with a large amount of noise, groups irregularly shaped in attribute space, and species turnover within groups, all characteristics commonly found in data originating from trawl surveys. It was applied for the first time in the marine realm with data from the Flemish Cap and a portion of the Grand Bank of Newfoundland (Murillo et al. 2016), and since then it has been more widely used (e.g., Beckers et al. 2018; Cabido et al. 2018; Gafna et al. 2021), including in other marine studies (e.g., Yu et al. 2022). This algorithm allows the identification of diagnostic species during the classification process at different scale patterns and the validation of the groups derived from HMSC, using independent data, has shown consistency of these

assemblages and diagnostic species. Others have used clustering techniques to identify species assemblages and then used random forest or similar models to predict their distributions (e.g., O'Brien et al. 2022), following a cluster then predict approach. We chose to use JSDM to allow each species to respond independently to the environment and to cross-compare the resultant RCPs with the ISOPAM results. This is only meaningful when the data are recorded to the species level, as has been done in this study.

The strong positive relationship found between individual species model performance from HMSC and the fidelity values from the ISOPAM algorithm allows us to confirm what other authors have found, namely, species with limited geographic range and environmental tolerances present better model fit than widespread species (Guisan et al. 2013; Morán-Ordóñez et al. 2017). Some common species, such as Atlantic lobster or snow crab, that reach maximum fidelity values at the first lower of the classification for the Bay of Fundy-Western Scotian Shelf and Eastern Scotian Shelf major groups, showed very good model fit (Tjur R² of 0.70 and 0.60, respectively). Other less common species, such as the crustaceans *Metacrangon agassizii* and *Stereomastis sculpta*, besides their low occurrence in the study area (both present in only 8 stations) had high performance associated with maximum fidelities at the third level of the classification. These species likely have defined environmental thresholds that can be identified at different scales. For example, the Atlantic lobster and snow crab have defined temperature thresholds (Zisserson and Cook 2017; Greenan et al. 2019) that related them to the different water masses present in the area, allowing the model to delineate their most suitable environmental conditions. Whereas for the slope species, depth seems the main variable delineating their most suitable environmental conditions. Other species, such as the brittle star *Ophiopholis aculeata* (Tjur R² = 0.24, Φ = 0.17) or the sea star *Hippasteria phrygiana* (Tjur R² = 0.2, Φ = 0.27) can be found on a wide range of substrates (Birkeland 1974; Mah 2015; Volage et al. 2021) and are widely distributed along the shelf and slope (Table S2). This results in poor model performance of these species and limited capacity to characterize communities at this scale. It is important to consider that benthic habitats are constituted by a mix of species, some of them with defined environmental thresholds that coexist with others presenting wider distributions. Therefore, approaches that consider species-specific responses offer some distinct advantages over aggregate community responses as they can account for differing levels of environmental affinity (e.g., Ovaskainen and Abrego 2020).

## Taxonomic resolution and sampling tools required for HMECS

Data from the trawl surveys have the advantage of being drawn from a stratified random design using standardized collection methods. The surveys are capable of covering large areas and capture both fish and invertebrates, and so are the best means of collecting data for Level 4 of the HMECS (DFO 2016). However, coverage is not complete, as some areas are avoided due to rough terrain, and there are catchability issues
associated with fishing performance and the morphology of
the species (Staniforth et al. 2023).

Benthic data of the taxonomic resolution obtained for this
study is generally not available over broad temporal and spa-
tial scales, although identification of the invertebrate catch
from DFO research vessel surveys has greatly improved in the
last decade (Nozères et al. 2015). Due to this limitation, a clas-
sification using oceanographic variables alone was proposed
assuming that environmental drivers, weighted by their rel-
evance to biological structure, could be correlated to biotic
composition (DFO 2016). That analysis identified ten physical
units on the Scotian Shelf Bioregion at the scale of of hun-
dreds to thousands of km (Level 4 of the HMECS) that were as-
sumed to be biophysical units (DFO 2016). The results of that
classification had some similarities to the present study, but
failed to capture the similarity in terms of benthic commu-
nities between the inner Bay of Fundy and Georges Bank, or
between Emerald Basin and the Fundian Channel, for exam-
ple. Nevertheless, these physical drivers are monitored annu-
ally and their general correspondence to the major ISOPAM
groups can be used to assess when changes to the benthic
communities may be occurring.

More recently, O'Brien et al. (2022) identified six assem-
blages of demersal fish and benthic invertebrates in the re-
gion using 10 years of data (2007–2016). Benthic invertebrates
from the analyses of O'Brien et al. (2022) included 49 species
from 6 phyla, with 20 of those being arthropods and 17 echin-
oderms. Nomenclature was based on at-sea identifications
from various technicians over the time series. Their analy-
ses did not include species that required taxonomic identifi-
cation in the laboratory such as sponges, bryozoans, or hy-
droids, and therefore did not provide a full characterization
of the benthic invertebrate communities. The 99 taxa from 9
phyla analysed in this study were drawn from 317 taxa col-
lected, which highlights the diversity of the epibenthos in
this region. Further, O'Brien et al. (2022) used random for-
est models based on the cluster grouping of each station to
produce spatial maps, while in this study we applied JSDMs
to jointly model the species in relation to their environment
and the clusters were produced posteriorly.

Despite these very different data sources and analytical ap-
proaches, there was general similarity between the composi-
tion of the assemblages identified by O'Brien et al. (2022) and
those identified here at HMECS Level 4, in part because the
diagnostic species for the three major cluster groups iden-
tified by the ISOPAM analyses were included in the O'Brien
data set (Homarus americanus (I. Western Scotian Shelf); Pon-
tophilus norvegicus (II. Central Scotian Shelf and Slope); Chio-
noecetes opilio (III. Eastern Scotian Shelf)) and similar environ-
mental variables (depth, average maximum and minimum
temperature and average maximum salinity) were the most
important characterizing their groups. At lower levels in the
ISOPAM classification, the congruence between the studies
breaks down, especially where assemblages were character-
ized by sponges and taxa not included in the previous study.
The Browns Bank and outer Bay of Fundy assemblage (I.b)
characterized by sponges is clearly a discrete assemblage, and
the St. Anns and shallow Misaine Banks assemblage (III.a.2)
was grouped with Banquereau Bank by O'Brien et al. (2022).

We conclude that using data with the higher taxonomic
resolution is a prerequisite for identification of robust ben-
thic assemblages at the Level 6 of the HMECS (DFO 2016),
characterized by species with high fidelity to their groups.
However, the use of data from the trawl surveys precludes
delineation of Biological Facies (HMECS Level 7; DFO 2016)
which are groups of biogenic habitat or foundation species
identified by one or more indicator species resolved at spatial
scales of hundreds of meters (DFO 2016), due to the integra-
tion of the catch over a 3 km distance. Previous studies within
the study area conducted at smaller scales using other sam-
pling gears such as a video-monitored sled-dredge on Georges
Bank (Thouzeau et al. 1991), drop cameras on Browns Bank
(Kostylev et al. 2001) and St. Anns Bank (Lacharité and Brown
2019), or remotely operated vehicles in the Gully canyon
(Kenchington et al. 2014). Each identified several benthic as-
semblages where this analysis only found one community
type. These other studies were able to characterize the dif-
ferent benthic habitats at smaller scales and associate them
to one unique substrate type, whereas our sampling stations,
due to the length of the tow can cover different substrate
types and therefore combine species with different substrate
affinities. A nested sampling design within each of the eight
RCPs using different gear types is required to further delin-
eate species assemblages to HMECS Level 7 where biological
facies are identified.

Conclusions

This study produced an authoritative map of benthic
biotopes for the DFO Maritimes region to 1500 m. Nine assem-
blages of benthic invertebrates were identified with ISOPAM
and mapped onto eight RCPs using HMSC. Diagnostic species
for each assemblage were identified and included a wide
range of taxa not considered in previous analyses (e.g., hy-
drozoans, bryozoans, and sponges). These species could be
good candidates for monitoring changes to benthic ecosys-
tems, which have been shown to be undergoing decadal-scale
changes in the Bay of Fundy in response to fishing pressures
and climate change (Staniforth et al. 2023). The identification
of the diagnostic species for each assemblage and its affilia-
tion with its RCP were tested using an independent data set
and found to be robust, producing repeatable results in both
instances. The map of the benthic biotopes can be used by
ocean managers in their selection of representative areas for
MPA network design and marine spatial planning.

Comparison of the results from this study with those of
previous studies suggests that for the delineation of biophys-
ical units (Level 4 of HMECS) the degree of taxonomic reso-
lution in our data is not necessary; however, for the identifi-
cation of biotopes (Level 6 of HMECS) it is essential. The ana-
lytical approach employed allowed for a more rigorous iden-
tification of the species assemblages through ISOPAM than
is obtained from traditional clustering algorithms, while the
use of a JSDM through HMSC enabled spatial mapping of the
ISOPAM assemblages through consideration of the probabil-
ity of occurrence of each of the 99 species and selection of
the number of RCPs through optimizing the matching of the
designations from both approaches.
Data from trawl surveys are useful for defining biophysical units and biotopes of epibenthic species. However, they integrate data at scales of 3 km and so are not able to resolve biological facies (Level 7) which would require different sampling tools such as ROVs, drop cameras, or grabs to sample at smaller spatial scales within the biotopes.

## Acknowledgements

The authors are especially grateful to the DFO science staff and the heads of the surveys for their help at sea and facilitating the data collection. They also acknowledge the Captains and crew on-board the CCGS *Alfred Needler* for their cooperation and professionalism throughout the 2017 summer RV survey. C. Goodwin (Huntsman Marine Science Centre) and E. Rodriguez (American Museum of Natural History) are thanked for their collaboration on the identifications of sponge and sea anemone specimens, respectively. They are also indebted to V. Kostylev (Natural Resources Canada) for providing the sediment grain size data. They thank D. Keith, R. Stanley, and C. Nozères at the Bedford Institute of Oceanography for helpful comments on an earlier version of the manuscript.

## Article information

### History dates

Received: 13 November 2023

Accepted: 25 July 2024

Version of record online: 14 November 2024

### Copyright

© 2024 Author Weigel and © His Majesty the King in Right of Canada, Department of Fisheries and Oceans, 2024. This work is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

### Data availability

All species distribution records environmental variables, and raster layers resulting from HMSC are available at Mendeley Data 10.17632/n8yk8rds9y.1.

## Author information

### Author ORCIDs

Francisco Javier Murillo https://orcid.org/0000-0001-7574-5558

Benjamin Weigel https://orcid.org/0000-0003-2302-5529

Ellen Kenchington https://orcid.org/0000-0003-3784-4533

### Author contributions

Conceptualization: FJM, BW, EK

Data curation: FJM, EK

Formal analysis: FJM, BW

Funding acquisition: FJM, EK

Investigation: FJM, DC

Methodology: FJM, BW, EK

Project administration: FJM, EK

Writing – original draft: FJM, BW, DC, EK

## Competing interests

The authors declare that there are no competing interests.

## Funding information

Funding for this project was through Fisheries and Oceans Canada, Maritimes Region projects "Marine Conservation Target project Mapping Biodiversity and Ecosystem Services of Benthic Communities on the Scotian Shelf in Support of Marine Spatial Planning and Area-Based Management" to EK, FJM (Project ID 385), and a Strategic Program for Ecosystem-Based Research and Advice (SPERA) project "Benthic Community Mapping on the Scotian Shelf" (2017–2019) to EK.

## Supplementary material

Supplementary data are available with the article at https://doi.org/10.1139/cjfas-2023-0326.

## References

Austin, M.P. 2002. Spatial prediction of species distribution: an interface between ecological theory and statistical modelling. Ecol. Modell. 157: 101–118. doi:10.1016/S0304-3800(02)00205-3.

Beazley, L., Wang, Z., Kenchington, E., Yashayaev, I., Rapp, H.T., Xavier, J.R., et al. 2018. Predicted distribution of the glass sponge *Vazella pourtalesi* on the Scotian Shelf and its persistence in the face of climatic variability. PLoS One, 13: e0205505. doi:10.1371/journal.pone.0205505. PMID: 30356324.

Beckers, N., Hein, N., Vanselow, K.A., and Löffler, J. 2018. Effects of microclimatic thresholds on the activity abundance and distribution patterns of alpine Carabidae species. Ann. Zool. Fenn. 55: 25–44. doi:10.5735/086.055.0104.

Béjar Alonso, J. 2013. Strategies and algorithms for clustering large datasets: a review. Computer Science.

Birkeland, C. 1974. Interactions between a sea pen and seven of its predators. Ecol. Monogr. 44: 211–232. doi:10.2307/1942312.

Bloomfield, N.J., Knerr, N., and Encinas-Viso, F. 2018. A comparison of network and clustering methods to detect biogeographical regions. Ecography, 41: 1–10. doi:10.1111/ecog.02596.

Brickman, D., Wang, Z., and DeTracey, B. 2016. Variability of current streams in Atlantic Canadian waters: a model study. Atmos. Ocean, 54: 218–229. doi:10.1080/07055900.2015.1094026.

Cabido, M., Zeballos, S.R., Zak, M., Carranza, M.L., Giorgis, M.A., Cantero, J.J., and Acosta, A.T.R. 2018. Native woody vegetation in central Argentina: classification of Chaco and Espinal forests. Appl. Veg. Sci. 21: 298–311. doi:10.1111/avsc.12369.

CBD. 2022. Kunming-Montreal Global biodiversity framework. CBD/COP/15/L.25, 14p.

Chytrý, M., and Tichý, L. 2003 Diagnostic, constant and dominant species of vegetation classes and alliances of the Czech Republic: a statistical revision. Folia Fac. Sci. Nat. Univ. Masaryk. Brun. Biol. 108: 1–231.

De Cáceres, M., Jansen, F., and Dell, N. 2023. Relationship between species and groups of sites. R Package ‘indicspecies’, version 1.7.13.

DFO. 2012. Marine protected area network planning in the Scotian Shelf bioregion: objectives, data, and methods. DFO Can. Sci. Advis. Sec. Sci. Advis. Rep. 2012/064.

DFO. 2016. Evaluation of hierarchical marine ecological classification systems for pacific and maritimes regions. DFO Can. Sci. Advis. Sec. Sci. Advis. Rep. 2016/003.

DFO. 2018. Design strategies for a network of marine protected areas in the Scotian Shelf bioregion. DFO Can. Sci. Advis. Sec. Sci. Advis. Rep. 2018/006.
<figure><img src="image_14.png" /></figure>

DFO. 2023. Reaching Canada's marine conservation targets. Fisheries and Oceans Canada. Available from https://www.dfo-mpo.gc.ca/oceans/conservation/plan/index-eng.html [accessed 5 August 2023].

ESRI. 2011. ArcGIS Desktop: Release 10. Environmental Systems Research Institute, Redlands, CA.

Furlong, M.H., and Kostylev, V.E. 2008. Methodology and implementation of a Sediment Texture Interpolation (MISTI) Tool: a Scotian Shelf case study. Geological Survey of Canada, Open File 5820, 69p., 1 CD-ROM.

Gafna, D.J., Obando, J.A., Reichelt, M., Schmidtlein, S., and Dolos, K. 2021. Medicinal service supply by wild plants in Samburu, Kenya: comparisons among medicinal plant assemblages. Global Ecology and Conservation 30: e01749. doi:10.1016/j.gecco.2021.e01749.

Gelman, A., and Rubin, D.B. 1992. Inference from iterative simulation using multiple sequences. Stat. Sci. 2: 45-52.

Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., and Rubin, D.B. 2013. Bayesian Data Analysis, 3rd ed. Chapman and Hall/CRC. 675pp.

Global Administrative Areas. 2011. GADM database of Global Administrative Areas, version 2.0 [online]. Available from www.gadm.org.

Government of Canada. 2011. National Framework for Canada's Network of Marine Protected Areas. Fisheries and Oceans Canada, Ottawa. 31pp.

Greenan, B.J.W., Shackell, N.L., Ferguson, K., Greyson, P., Cogswell, A., Brickman, D., et al. 2019. Climate change vulnerability of American lobster fishing communities in Atlantic Canada. Front. Mar. Sci. 6: 579. doi:10.3389/fmars.2019.00579.

Guisan, A., Tingley, R., Baumgartner, J.B., Naujokaitis-Lewis, I., Sutcliffe, P.R., Tulloch, A.I.T., et al. 2013. Predicting species distributions for conservation decisions. Ecol. Lett. 16: 1424-1435. doi:10.1111/ele.12189.

Halliday, R.G., and Kohler, A.C. 1971. Groundfish survey programmes of the St. Andrews Biological Station, Fisheries Research Board of Canada—objectives and characteristics. ICNAF Res. Doc. 71/35.

Han, G., Loder, J.W., and Smith, P.C. 1999. Seasonal-mean hydrography and circulation in the Gulf of St. Lawrence and eastern Scotian and southern Newfoundland Shelves. J. Phys. Oceanogr. 29: 1279-1301. doi:10.1175/1520-0485(1999)029%3c1279:SMHACI%3e2.0.CO;2.

Hannah, C.G., Shore, J.A., and Loder, J.W. 2001. Seasonal circulation on the western and central Scotian Shelf. J. Phys. Oceanogr. 31: 591-615. doi:10.1175/1520-0485(2001)031%3c0591:SCOTWA%3e2.0.CO;2.

Harris, D.J. 2015. Generating realistic assemblages with a joint species distribution model. Methods Ecol. Evol. 6: 465-473. doi:10.1111/2041-210X.12332.

Hawkes, N., Korabik, M., Beazley, L., Rapp, H.T., Xavier, J.R., and Kenchington, E. 2019. Glass sponge grounds on the Scotian Shelf and their associated biodiversity. Mar. Ecol. Prog. Ser. 614: 91-109. doi:10.3354/meps12903.

Henry, L.-A., Frank, N., Hebbeln, D., Wienberg, C., Robinson, E., de Flierdt, T., et al. 2014. Global ocean conveyor lowers extinction risk in the deep sea. Deep Sea Res. Part I. 88: 8-16. doi:10.1016/j.dsr.2014.03.004.

Hill, M.O., and Gauch, H.G., Jr. 1980. Detrended correspondence analysis: an improved ordination technique. Vegetatio, 42: 47-58. doi:10.1007/BF00048870.

Hill, N., Woolley, S.N.C., Foster, S., Dunstan, P.K., McKinlay, J., Ovaskainen, O., and Johnson, C. 2020. Determining marine bioregions: a comparison of quantitative approaches. Methods Ecol. Evol. 11: 1258-1272. doi:10.1111/2041-210X.13447.

Kenchington, E.L., Cogswell, A.T., Maclsaac, K.G., Beazley, L., Law, B.A., and Kenchington, T.K. 2014. Limited depth zonation among bathyal epibenthic megafauna of the Gully submarine canyon, northwest Atlantic. Deep Sea Res. Part II, 104: 67-82. doi:10.1016/j.dsr2.2013.08.016.

King, L.H. 1980. Aspects of regional surficial geology related to site investigation requirements. In Offshore site investigation. Edited by D.A. Ardus. The Society for Underwater Technology. Graham and Trotman Ltd, London, pp. 37-60.

King, L.H., and Fader, G.B. 1986. Wisconsinan glaciation of the Atlantic continental shelf of southeast Canada. Geol. Surv. Can. Bull. 363: 72pp.

Koen-Alonso, M., Favaro, C., Ollerhead, N., Benoît, H., Bourdages, H., Sainte-Marie, B., et al. 2018. Analysis of the overlap between fishing effort and Significant Benthic Areas in Canada's Atlantic and Eastern

Arctic marine waters. DFO Can. Sci. Advis. Sec. Res. Doc. 2018/015. XVII + 270p.

Kostylev, V.E., and Hannah, C.G. 2007. Process-driven characterization and mapping of seabed habitats. In Mapping the Seafloor for Habitat Characterization. Edited by B.J. Todd and H.G. Greene. Geological Association of Canada, Special Paper 47. pp. 171-184.

Kostylev, V.E., Todd, B.J., Fader, G.B.J., Courtney, R.C., Cameron, G.D.M., and Pickrill, R.A. 2001. Benthic habitat mapping on the Scotian Shelf based on multibeam bathymetry, surficial geology and sea floor photographs. Mar. Ecol. Prog. Ser. 219: 121-137. doi:10.3354/meps219121.

Lacharité, M., and Brown, C.J. 2019. Utilizing benthic habitat maps to inform biodiversity monitoring in marine protected areas. Aquat. Conserv. Mar. Freshwater Ecosyst. 29: 938-951. doi:10.1002/aqc.3074.

Lecours, V. 2017. On the use of maps and models in conservation and resource management (warning: results may vary). Front. Mar. Sci. 4: 288. doi:10.3389/fmars.2017.00288.

Legendre, P., and Legendre, L. 1998. Numerical ecology, 2nd ed. Elsevier, Amsterdam.

Levin, L.A., and Sibuet, M. 2012. Understanding continental margin biodiversity: a new imperative. Ann. Rev. Mar. Sci. 4: 79-112. doi:10.1146/annurev-marine-120709-142714.

Levin, P.S., Fogarty, M.J., Murawski, S.A., and Fluharty, D. 2009. Integrated ecosystem assessments: developing the scientific basis for ecosystem-based management of the ocean. PLoS Biol. 7: 23-28. doi:10.1371/journal.pbio.1000014.

Li, M.Z., Wu, Y., Hannah, C.G., Perrie, W.A., Shen, H., and King, E.L. 2021. Modelling seabed disturbance and sediment mobility on the Canadian Atlantic Shelf; Geological Survey of Canada, Open File 8805, 50p.

Livingston, P.A., Aydin, K., Boldt, J.L., Hollowed, A.B., and Napp, J.M. 2011. Alaska marine fisheries management: advancements and linkages to ecosystem research. In Ecosystem Based Management: An Evolving Perspective. Edited by A. Belgrano and C. Fowler. Cambridge University Press, UK. pp. 113-152.

Loder, J.W., Petrie, B., and Gawarkiewicz, G. 1998. The coastal ocean off northeastern North America: a large-scale view. In The Sea, Vol. 11, The Global Coastal Ocean: Regional Studies and Syntheses. Edited by A.R. Robinson and K.H. Brink. Harvard University Press, Cambridge, MA. pp. 105-133.

Long, D.A., Charles, A., and Stephenson, R.L. 2015. Key principles of marine ecosystem-based management. Mar. Policy, 57: 53-60. doi:10.1016/j.marpol.2015.01.013.

Mah, C.L. 2015. A new Atlantic species of *Evoplosoma* with taxonomic summary and in situ observations of Atlantic deep-sea corallivorous Goniasteridae (Valvatida; Asteroidea). Mar. Biodivers. Rec. 8: e5. doi:10.1017/S1755267214001407.

Mahon, R., and Smith, R.W. 1989. Demersal fish assemblages on the Scotian Shelf, northwest Atlantic: spatial distribution and persistence. Can. J. Fish. Aquat. Sci. 46: 134-152. doi:10.1139/f89-285.

Mahon, R., Smith, R.W, Bernstein, B.B., and Scott, J.S. 1984. Spatial and temporal patterns of groundfish distribution on the Scotian Shelf and in the Bay of Fundy, 1970-1981. Can. Tech. Rep. Fish. Aquat. Sci. 130: xi+164.

Metaxas, A., Lacharité, M., and de Mendonça, S.N. 2019. Hydrodynamic connectivity of habitats of deep-water corals in Corsair Canyon, Northwest Atlantic: a case for cross-boundary conservation. Front. Mar. Sci. 6: 159. doi:10.3389/fmars.2019.00159.

Montanyès, M., Weigel, B., and Lindegren, M. 2023. Community assembly processes and drivers shaping marine fish community structure in the North Sea. Ecography, 2023: e06642. doi:10.1111/ecog.06642.

Morán-Ordóñez, A., Lahoz-Monfort, J.J., Elith, J., and Wintle, B.A. 2017. Evaluating 318 continental-scale species distribution models over a 60-year prediction horizon: What factors influence the reliability of predictions? Global Ecol. Biogeogr. 26: 371-384. doi:10.1111/geb.12545.

Murillo, F.J., Kenchington, E., Clark, D., Emberley, J., Regnier-McKellar, C., Guijarro, J., et al. 2018a. Cruise Report for the CCGS Alfred Needler Maritimes Region Research Vessel Summer Multispecies Survey, June 28 to August 14, 2017: Benthic Invertebrates. Can. Tech. Rep. Fish. Aquat. Sci. 3262: v+41.
<figure><img src="image_15.png" /></figure>

Murillo, F.J., Kenchington, E., Tompkins, G., Beazley, L., Baker, E., Knudby, A., and Walkusz, W. 2018b. Sponge assemblages and predicted archetypes in the eastern Canadian Arctic. Mar. Ecol. Prog. Ser. 597: 115-135. doi:10.3354/meps12589.

Murillo, F.J., Serrano, A., Kenchington, E., and Mora, J. 2016. Epibenthic assemblages of the Tail of the Grand Bank and Flemish Cap (north-west Atlantic) in relation to environmental parameters and trawling intensity. Deep Sea Res. Part I, 109: 99-122. doi:10.1016/j.dsr.2015.08.006.

Murillo, F.J., Weigel, B., Bouchard Marmen, M., and Kenchington, E. 2020. Marine epibenthic functional diversity on Flemish Cap (north-west Atlantic)—identifying trait responses to the environment and mapping ecosystem functions. Divers. Distrib. 26: 460–478. doi:10.1111/ddi.13026.

Nesis, K.N. 1965. Biocoenoses and biomass of benthos of the Newfoundland-Labrador region. Trudy VNIRO, 57: 453-489. Fish. Res. Bd. Can. Translation Series No. 1375 (Transl. from Russian by Transl. Bur. Foreign Lang. Div., Dep. Secr. State Can.), 75 p.

Neumann, H., Reiss, H., Ehrich, S., Sell, A., Panten, K., Kloppmann, M., et al. 2013. Benthos and demersal fish habitats in the German Exclusive Economic Zone (EEZ) of the North Sea. Helgol. Mar. Res. 67: 445-459. doi:10.1007/s10152-012-0334-z.

Nozères, C., Bourassa, M., Gendron, M., Plourde, S., Savenkoff, C., Bourdages, H., et al. 2015. Using annual ecosystemic multispecies surveys to assess biodiversity in the Gulf of St. Lawrence. Can. Tech. Rep. Fish. Aquat. Sci. 3149: 126.

O'Brien, J.M., Stanley, R.R.E., Jeffery, N.W., Heaslip, S.G., DiBacco, C., and Wang, Z. 2022. Modeling demersal fish and benthic invertebrate assemblages in support of marine conservation planning. Ecol. Appl. 32. doi:10.1002/eap.2546.

Oceans Act, SC. 1996. c 31, Available from https://canlii.ca/t/5439k [accessed 27 July 2023].

Oksanen, J., Blanchet, F.G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. 2018. Community Ecology Package. R package 'vegan' version 2.5-2. Available from https://cran.r-project.org/web/packages/vegan/vegan.pdf.

Ovaskainen, O., and Abrego, N. 2020. Joint species distribution modelling. Cambridge University Press. 372pp.

Ovaskainen, O., Abrego, N., Halme, P., and Dunson, D. 2016. Using latent variable models to identify large networks of species-to-species associations at different spatial scales. Methods Ecol. Evol. 7: 549-555. doi:10.1111/2041-210X.12501.

Ovaskainen, O., Tikhonov, G., Norberg, A., Blanchet, F.G., Duan, L., Dunson, D., et al. 2017. How to make more out of community data? A conceptual framework and its implementation as models and software. Ecol. Lett. 20: 561-576. doi:10.1111/ele.12757.

Pearce, J., and Ferrier, S. 2000. Evaluating the predictive performance of habitat models developed using logistic regression. Ecol. Modell. 133: 225-245. doi:10.1016/S0304-3800(00)00322-7.

Philibert, G., Todd, B.J., Campbell, D.C., King, E.L., Normandeau, A., Hayward, S.E., et al. 2022. Updated surficial geology compilation of the Scottish Shelf bioregion, offshore Nova Scotia and New Brunswick, Canada; Geological Survey of Canada, Open File 8911, 1. doi:10.4095/330474.

Piper, D.J.W., and Campbell, D.C. 2002. Surficial geology of the Scottish Slope, eastern Canada; Geological Survey of Canada, Current Research 2002-E15, 10p.

Rincón, B., and Kenchington, E. 2016. Influence of Benthic Macrofauna as a Spatial Structuring Agent for Juvenile Haddock (*Melanogrammus aeglefinus*) on the Eastern Scotian Shelf, Atlantic Canada. PLoS One, 20; e0163374. doi:10.1371/journal.pone.0163374.

Rubidge, E.M., Gale, K.S.P., and Curtis, J.M.R. 2016. Community ecological modelling as an alternative to physiographic classifications for marine conservation planning. Biodivers. Conserv. 25: 1899-1920. doi:10.1007/s10531-016-1167-x.

Schmidtlein, S., Tichý, L., Feilhauer, H., and Faude, U. 2010. A brute-force approach to vegetation classification. J. Veg. Sci. 21: 1162-1171. doi:10.1111/j.1654-1103.2010.01221.x.

Schmidtlein, S., Collison, J., Pfannendoerfer, R., and Tichy, L. 2022. Clustering of Sites with Species Data. R package 'isopam' version 1.1.0.

Seidov, D., Mishonov, A., and Parsons, R. 2021. Recent warming and decadal variability of Gulf of Maine and Slope Water. Limnol. Oceanogr. 66: 3472-3488. doi:10.1002/lno.11892.

Sokal, R.R., and Rohlf, F.J. 1995. Biometry: the principles and practice of statistics in biological research. 3rd ed. Freeman, New York, NY, US.

Staniforth, C., MacDonald, B.W., Lirette, C., Murillo, F.J., Kenchington, E., and Kenchington, T. 2023. Drivers of decadal-scale temporal change in epibenthic megafaunal assemblages on scallop fishing grounds in the Bay of Fundy. J. Sea Res. 195: 102419. doi:10.1016/j.seares.2023.102419.

Stanley, R.R.E., DiBacco, C., Lowen, B., Beiko, R.G., Jeffery, N.W., Wyngaarden, M.V., et al. 2018. A climate-associated multispecies cryptic cline in the northwest Atlantic. Sci. Adv. 4: eaaq092. doi:10.1126/sciadv.aaq0929.

Stephenson, F., Leathwick, J.R., Geange, S.W., Bulmer, R.H., Hewitt, J.E., Anderson, O.F., et al. 2018. Using gradient forests to summarize patterns in species turnover across large spatial scales and inform conservation planning. Divers. Distrib. 24: 1641-1656. doi:10.1111/ddi.12787.

Stephenson, F., Rowden, A.A., Tablada, J., Tunley, K., Brough, T., Lundquist, C.J., et al. 2023. A seafloor bioregionalisation for New Zealand. Ocean Coast. Manag. 242: 106688. doi:10.1016/j.ocecoaman.2023.106688.

Ter Braak, C.J.F. 1986. Canonical Correspondence Analysis: a new eigenvector technique for multivariate direct gradient analysis. Ecology, 67: 1167-1179. doi:10.2307/1938672.

Thompson, P.L., Anderson, S.C., Nedin, J., Haggarty, D.R., Peña, M.A., English, P.A., et al. 2022. Disentangling the impacts of environmental change and commercial fishing on demersal fish biodiversity in a northeast Pacific ecosystem. Mar. Ecol. Prog. Ser. 689: 137-154. doi:10.3354/meps14034.

Thouzeau, G., Robert, G., and Ugarte, R. 1991. Faunal assemblages of benthic megainvertebrates inhabiting sea scallop grounds from eastern Georges Bank, in relation to environmental factors. Mar. Ecol. Prog. Ser. 74: 61-82. doi:10.3354/meps074061.

Tichý, L., and Chytrý, M. 2006. Statistical determination of diagnostic species for site groups of unequal size. J. Veg. Sci. 17: 809-818. doi:10.1111/j.1654-1103.2006.tb02504.x.

Tikhonov, G., Ovaskainen, O., Oksanen, J., de Jonge, M., Opedal, O., and Dallas, T. 2019a. Hierarchical Model of Species Communities. R package version 3.0-2. Available from https://cran.r-project.org/web/packages/Hmsc/index.html.

Tikhonov, G., Øystein, O., Abrego, N., Lehikoinen, A., and Ovaskainen, O. 2019b. Joint species distribution modelling with HMSC-R. doi:10.1101/603217. bioRxiv.

Tjur, T. 2009. Coefficients of determination in logistic regression models—a new proposal: the coefficient of discrimination. Am. Stat. 63: 366-372. doi:10.1198/tast.2009.08210.

Tremblay, M.J., Black, G.A.P., and Branton, R.M. 2007. The distribution of common decapod crustaceans and other invertebrates recorded in annual ecosystem surveys of the Scotian Shelf 1999-2006. Can. Tech. Rep. Fish. Aquat. Sci. 2762: iii+74.

Van Wyngaarden, M., Snelgrove, P.V.R., DiBacco, C., Hamilton, L.C., Rodríguez-Ezpeleta, N., Jeffery, N.W., et al. 2017. Identifying patterns of dispersal, connectivity and selection in the sea scallop, *Placopecten magellanicus*, using RADseq-derived SNPs. Evol. Appl. 10: 102-117. doi:10.1111/eva.12432.

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. 2021. Rank-Normalization, folding, and localization: an improved R for assessing convergence of MCMC (with Discussion). Bayesian Analysis, 16: 667-7181. doi:10.1214/20-BA1221.

Volage, F., Hamel, J.-F., and Mercier, A. 2021. Population structure, habitat preferences, feeding strategies, and diet of the brittle star *Ophiopholis aculeata* in nearshore and offshore habitats of the northwest Atlantic. Invertebr. Biol. 140: e12346. doi:10.1111/ivb.12346.

Wang, Z., Lu, Y., Greenan, B., and Brickman, D. 2018. BNAM: an Eddy-Resolving North Atlantic Ocean Model to Support Ocean Monitoring in Can. Tech. Rep. Hydrogr. Ocean Sci, 327: vii+18.

Wang, S., Kenchington, E., Wang, Z., and Davies, A.J. 2021. Life in the Fast Lane: Modeling the Fate of Glass Sponge Larvae in the Gulf Stream. Front. Mar. Sci. 8. doi:10.3389/fmars.2021.701218.

Wang, S., Murillo, F.J., and Kenchington, E. 2022. Climate-Change Refugia for the Bubblegum Coral *Paragorgia arborea* in the Northwest Atlantic. Front. Mar. Sci. 9. doi:10.3389/fmars.2022.863693.
Warton, D.I., Blanchet, F.G., O'Hara, R.B., Ovaskainen, O., Taskinen, S., Walker, S.C., and Hui, F.K.C. 2015. So many variables: joint modeling in community ecology. *Trends Ecol. Evol.* **30**: 766–779. doi:10.1016/j.tree.2015.09.007.

Weigel, B., Kotamäki, N., Malve, O., Vuorio, K., and Ovaskainen, O. 2023. Macrosystem community change in lake phytoplankton and its implications for diversity and function. *Global Ecol. Biogeogr.* **32**: 295–309. doi:10.1111/geb.13626.

Williams, D.A. 1976. Improved likelihood ratio tests for complete contingency tables. *Biometrika*, **63**: 33–37. doi:10.1093/biomet/63.1.33.

Williams, G.C., and Van Syoc, R.J. 2007. The Light and Smith Manual: intertidal invertebrates from central California to Oregon. *In* Methods of preservation and anesthetization of marine invertebrates. *Edited by* J.T. Carleton. University of California Press, Berkeley, CA.

Yasuhara, M., and Danovaro, R. 2014. Temperature impacts on deep-sea biodiversity. *Biol. Rev.* **91**: 275–287. doi:10.1111/brv.12169.

Yu, H., Fang, G., Rose, K.A., Tang, Y., and Song, X. 2022. Examining epibenthic assemblages associated with artificial reefs using a species archetype approach. *Mar. Coast. Fish.* **14**: e10206. doi:10.1002/mcf2.10206.

Zisserson, B., and Cook, A. 2017. Impact of bottom water temperature change on the southernmost snow crab fishery in the Atlantic Ocean. *Fish. Res.* **195**: 12–18. doi:10.1016/j.fishres.2017.06.009.

Zwanenburg, K.C.T., and Jaureguizar, A.J., (Unpublished MS). 2010. Long-term structure and dynamics of demersal fish assemblages on the Scotian Shelf—Implications for ecosystem based management. *In* NAFO. Report of the 3rd Meeting of the NAFO Scientific Council Working Group on Ecosystem Approaches to Fisheries Management (WGEAFM). NAFO Scientific Council Summary Document 10/24, Series No. N5868, 75.

Zwanenburg, K., Horsman, T., and Kenchington, E. 2010. Preliminary analysis of biogeographic units on the Scotian Shelf. NAFO SCR Doc. 10/06.