Seafloor biodiversity of Canada's three oceans: Patterns, hotspots and potential drivers

We examined the relationships between bathymetry, latitude and energy and the diversity of marine benthic invertebrates across wide environmental ranges of Canada's three oceans.


| INTRODUC TI ON
Canada has the longest coastline of any country in the world (243,791 km), encompassing more territorial waters (2,687,667 km 2 within 12 nautical miles) than all European countries and the United States of America combined. The vast Canadian waters span the northeast Pacific, Arctic and northwest Atlantic Oceans, and can be further divided into coastal and shelf marine ecoregions ( Figure   1a, modified from Spalding et al., 2007), covering a surface area of approximately 2% of the total area of the world's oceans (DFO, 2008). We hereafter refer to the immense Canadian coastline and the expanse of the oceans as "Canada's Three Oceans" (Archambault et al., 2010). Moreover, as a signatory to the United Nations (UN) Convention on Biological Diversity (CBD) and under its own Oceans Act (1996), Canada is responsible for the "management of estuarine, coastal and marine ecosystems" including implementation of marine protected areas (MPAs) for "the conservation and protection of marine areas of high biodiversity or biological productivity." The rapid advancement in Canada's marine protected areas (MPAs) network in several marine bioregions is partly motivated to fulfil the Aichi 2020 targets to conserve 10% of coastal and marine areas (Aichi Biodiversity Targets, 2010). Despite recent progress in the establishment of marine protected areas (on 1 August 2019, about 794 000 km 2 or 13.9%; DFO, 2019), about half of these areas, as marine refuges, will allow activities such as drilling and fishing and will require good management plans beyond 2020 targets (Gies, 2019) as well as a better understanding of the marine biodiversity for conservation purposes. Existing efforts to monitor macrobenthos within Canadian protected areas vary by ecoregion and specific geographic scope (e.g. within small coastal MPAs vs. across large deep-water MPAs). To our knowledge, there is no systematic, consistently applied methodology.
The lack of synthesis for existing biodiversity datasets across Canada has hindered the assessment of marine biodiversity, a critical step towards identifying ecologically important areas for MPA planning and ecosystem-based conservation and management (Archambault et al., 2010;DFO, 2010;Myers, Mittermeier, Mittermeier, Fonseca, & Kent, 2000). This gap results from limited spatial and temporal coverage in individual sampling efforts that cannot evaluate large-scale patterns or provide comprehensive information on the scale of problems or threats faced within Canada's three oceans. As is true elsewhere, these threats include, among others, overfishing (Costello et al., 2016), invasive species (Goldsmit, Howland, & Archambault, 2014;Molnar, Gamboa, Revenga, & Spalding, 2008), climate change (Cheung et al., 2009), eutrophication and pollution (Diaz & Rosenberg, 2008), habitat destruction (Halpern et al., 2008) and sea-level rise (Rignot, Velicogna, Broeke, Monaghan, & Lenaerts, 2011). These threats punctuate the critical need to bring together existing datasets to establish a comprehensive assessment of biodiversity and form the basis to evaluate any concurrent or future impacts of climate-related changes and human activities on marine ecosystems. Such information can also inform efforts to identify ecologically important areas to aid effective conservation with MPA network design.
Benthic invertebrates are important bioindicators of ocean health and an essential component of monitoring and assessment programmes nationally and internationally (e.g. Borja & Dauer, 2008;Van Hoey et al., 2010). Their relatively low mobility as juveniles and adults generally integrates environmental conditions and perturbations over time (Bilyard, 1987), resulting in prolonged and maximum exposure to disturbances/stressors. Limited mobility also simplifies quantification, thus supporting broad ranges of indicators and diversity measures in developing strategies for ecosystem-based conservation. Moreover, they are trophically diverse and exhibit species-specific tolerances or sensitivity to organic enrichment and eutrophication, typically leading to changes in biodiversity patterns (Borja, Franco, & Pérez, 2000;Pearson & Rosenberg, 1978).
To date, the best available syntheses of benthic macrofaunal (≥0.5 or 1 mm in size) biodiversity for Canada's three oceans reported 2,127 species (Archambault et al., 2010;. However, this number of reported species underestimates biodiversity given: (a) inconsistent taxonomic information across datasets, (b) estimates were not corrected for unequal sampling effort among regions, (c) omission of known published biodiversity datasets, (d) lack of information on key benthic habitats and (e) lack of abundance information. The willingness/capacity for individual project investigators to contribute their biodiversity data further constrained this estimate.
In this study, we examine energy-diversity relationships by modelling the effects of temperature (thermal energy), food supply and the seasonal pulses of food supply (both chemical energy) on diversity across the wide environmental gradients of the Canadian Pacific, Arctic and Atlantic Ocean environments. We extended the previous effort (Archambault et al., 2010) by using additional quantitative, taxonomically standardized and spatially/temporally referenced community data to infer the underlying mechanism(s) driving seafloor biodiversity patterns. This more comprehensive database for Canada's three oceans of benthic biodiversity (TOBB database) provides a more robust baseline from which to evaluate the effectiveness of ecosystem-based management approaches, including the establishment of conservation goals and monitoring programmes for MPA networks and in assessing concomitant impacts of economic development and climate changes on benthic biodiversity within Canada's territorial waters.

| Data compilation
We compiled thirty-five datasets collected between 1953 and 2012 encompassing 13,172 samples from 6,117 sites among 13 marine ecoregions (Spalding et al., 2007), extending from intertidal zones to 1,072-m depth across Canada's three oceans (see Table S1, metadata information of the TOBB database). Although Fisheries and Oceans Canada defines ecoregions somewhat differently than Spalding et al. (2007), we chose the latter approach because it is well-recognized internationally, thus facilitating comparison with ecoregions elsewhere globally. Here, we define "sites" as the records in a single 2 × 2 km grid based on cylindrical equal-area projection. Various sampling gears were used in different datasets ( Figure 1; Table S1). Approximately 46.5% of records were collected with trawl, sledge or dredge (epifauna) summing to 86.9% of all sites, 24.4% with core or grab (infauna) adding to 11.3% of sites, and 28.6% of records were collected with quadrat (intertidal epifauna) representing 1.5% of all sites (Table   S2). Overall, about 39.1% of records include quantitative information (either count, density or biomass); 60.9% of records report only presence/absence ( Figure S1). The taxonomic information from each dataset was validated to the lowest possible taxonomic levels using the World Register of Marine Species (WoRMS), allowing comparison of taxonomic lists among the datasets. For incorrect original species names (from data providers) that are now assigned (by WoRMS) to more than one accepted species names, we selected the accepted species name with known records in or near Canada's three oceans.
We then used the unique "AphiaID_accepted" from WoRMS as an identifier to define taxonomic units across different datasets.
We downloaded decadal average  bottom temperature from the World Ocean Atlas (WOA) 2013 v2 (from 102 standard depth layers, Table S3; Figure S1) and calculated export POC flux (or food supplies) to the seafloor using Lutz et al.'s algorithm (2007) based on mean and seasonality (standard deviation/mean) of monthly surface net primary production (NPP, see below), as well as the export depth below the euphotic layer. We calculated euphotic layer depth from mean monthly chlorophyll-a (CHLA) concentrations using Morel and Berthon's (1989) case I model and obtained the monthly surface NPP (vertically generalized production model, VGPM) and CHLA concentrations (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) from the Ocean Productivity web portal (Table S3). We calculated export depth by subtracting the euphotic depth from water depth based on the ETOPO1 Global Relief Model (Table S3; Figure 1). All environmental layers were re-gridded by bilinear interpolation to 2 × 2 km grids based on cylindrical equal-area projection and extracted by the geographic coordinates of the sampling "sites." Appendix S1 provides greater detail on biotic and environmental data compilation (see also Figure S2). We did not extract environmental data for intertidal habitats and shallow shelf areas where calculated euphotic depth exceeds water depth.

| Data analysis
We measured benthic invertebrate biodiversity using Hill numbers (Hill, 1973), or the effective numbers of equally abundant species (Chao & Jost, 2010). The Hill number is defined as in which S is the number of species in a sample and p i is the relative abundance of ith species. The order q controls the sensitivity of Hill numbers to species relative abundance, and thus, a larger q gives more weight to common species. For example, q = 0 gives equal weights to all species, and thus, Equation 1 reduces to species richness. When q approximates 1, Equation 1 reduces to Equation 2 is thus the exponential of familiar Shannon entropy, measuring the effective numbers of common species. When q = 2, Equation 1 reduces to This equation is equivalent to the inverse of the Gini Simpson index, measuring the effective numbers of dominant species. We considered the three most widely used Hill numbers, including q = 0 (species richness), q = 1 (exponential of Shannon index) and q = 2 (inverse of Simpson index), allowing us to examine the diversity of total, common and dominant species based on the same unit (hypothetical numbers of species with equal abundance) in a unified framework (Chao & Jost, 2010;Magurran, 2004).
We analysed data separately by sampling gears and whether the data contain occurrence or abundance information. For example, we considered quadrat samples as occurrence data, because some of the quadrat data reported per cent cover and others reported occurrence of each species, whereas the core/grab and trawl/dredge/ sledge samples contain either occurrence or abundance data. Only the data containing abundance information were used to estimate the sample-based Hill numbers. For ecoregion-level analysis, all data were first converted to presence/absence and then to incidence probability within each region before computing the Hill numbers in the same way as the abundance data. In order to standardize sampling efforts, we bootstrap-resampled (1,000 times) one to m numbers of samples (from the occurrence data) and one to m numbers of individuals (from the abundance data) to produce smooth mean sampling curves of the Hill numbers ( Figures S3 and S4). We adapted a rarefaction and extrapolation method developed by Chao et al. (2013) to standardize the Hill numbers among ecoregions ( Figure   S3) or sites ( Figure S4). That is, when m ≤ observed numbers of the sample in an ecoregion/individual in a sample, we extracted the Hill numbers from the rarefaction part of the sampling curve; when m > observed numbers of the sample in an ecoregion/individual in a sample, we extracted the Hill numbers from the extrapolated part of the sampling curves. Given the heterogeneity in data distribution and varying sampling efforts ( Figure 1; Table S2), we computed Hill numbers from m = 1 to 1,000 randomly selected samples to construct sample-based accumulation curves across the three oceans.
For simplicity, the Hill numbers based on m = 200 randomly selected samples from each ecoregion (hundreds to thousands of kilometre scales) and m = 100 randomly selected individuals from each site (kilometre scale) were used for further analysis. We also hereafter refer to the Hill numbers of q = 0 to 2 from standardized sampling effort as "species richness," "exponential Shannon" and "inverse Simpson" indices, respectively. Appendix S1 provides details and rationales for using Hill's numbers and interpolation/extrapolation methods.
Only individual-based Hill numbers calculated from quantitative data (core/grab or trawl/sledge/dredge) were used for statistical modelling. Whereas core/grab data were available for all ecoregions, the quantitative trawl/sledge/dredge data were only available to us from the Arctic Ocean and mostly from Ecoregions 4 & 5 (i.e. the Beaufort Sea) with limited latitudinal variation (Table 1; Figure 1d). Despite dense sampling, the trawl/sledge/dredge in the Gulf of St. Lawrence only contains occurrence information ( Figure 1c). Therefore, we considered latitude and ecoregion in the analyses of core/grab data but not in the trawl/sledge/dredge data.
For the core/grab data, we used linear mixed-effects (LME) model to explore the relationship between Hill numbers (Table 1) and depth or latitudes. We used marine ecoregions as a random factor and set the intercept of each ecoregion to be normally distributed with independent variance structure to account for largescale variations in diversity estimates. Because unimodal diversity patterns typically occur along a large depth gradient (Rex & Etter, 2010), we added a quadratic depth term to the LME as an additional independent variable. Dependent and independent variables were averaged by site (2-km grid) to reduce the influence of sites with numerous replicates or repeated sampling (in the very limited dataset).
To prevent skewness and approximate normal distributions, we first transformed (log 10 ) and then normalized (subtracted by the mean and divided by the standard deviation) water depth and latitude before analysis.
To test for an energy-diversity relationship ( (export POC flux) as an independent variable, acknowledging that interpretation of export POC flux in coastal environments requires great caution. All environmental data were normalized (subtracted by the mean and divided by the standard deviation) before analysis.
We estimated the parameter coefficients of LMEs by model averaging from all possible combinations of the models and selected the best model based on the lowest Akaike information criterion, corrected for small sample size (AICc) (Burnham & Anderson, 2002).
The relative importance of the independent variable was estimated by summing their Akaike weights, which measures the likelihood of the presence of a particular factor in the average model. We assessed the model performance of the best LMEs (with lowest AICc) by 10-fold cross-validation, in which we validated the models re-fitted from every 9/10 of the data with the remaining 1/10 of data to calculate the total variance explained (R 2 ). The mean, standard deviation and 95% confidence interval of the model fit were calculated by 1,000 bootstraps and plotted against the independent variables to evaluate model behaviour.
We used the same sets of independent variables except for latitude to fit multiple linear regressions (without a random factor) to the Hill numbers of the trawl/sledge/dredge data, assessing model performance based on the adjusted R 2 of the linear regression. Even though this analysis used no sledge and dredge data, we nonetheless refer to the terms "trawl/sledge/dredge data" for consistency with the previous analysis.

| Diversity among oceans and ecoregions
Our database contains a total of 3,819 taxa: 2,670 identified to spe-  Figure 1; Table S2).
Pacific species richness accumulation curves in the core/grab analysis were most elevated, followed by the Arctic and Atlantic ( Figure 2a). This trend reverses, however, when weighting the common species. That is, the Pacific diversity curves were lower than the Arctic for exponential Shannon ( Figure 2b) and even lower than the Atlantic for inverse Simpson indices (Figure 2c). In other words, of the three oceans, Arctic samples yielded highest diversity for common and dominant species (Figure 2b,c). For the trawl/sledge/ dredge data, we observed consistently elevated diversity in Arctic accumulation curves (q = 0 to 2) relative to the Atlantic; however, we lack comparable data for the Pacific (Figure 2d-f). For the quadrat data, the Pacific mean diversity curves (q = 0 to 2) were approximately 75%, 70% and 50% more elevated than those for the Atlantic, respectively (Figure 2g-i). Although differences were not dramatic, their 95% confidence intervals were distinctly separate (the 95% CI are too narrow to be visible on the figures).
We also constructed diversity accumulation curves for each gear type and marine ecoregion to compare the three oceans ( Figure S3) and, for simplicity, extracted the results of 200 randomly selected samples (Figure 3). Among the core/grab data in the Pacific, highest species richness occurred in Ecoregion 2 ( Figure 3a); however, differences among ecoregions dissipated with increased weighting on common and dominant species (q = 1 to 2, Figure 3b,c). The highest diversities among the Arctic core/ grab data (q = 0 to 2, Figure 3a-c) were in Ecoregions 4, 7 and 9.
We also observed high diversity in Ecoregion 5, especially for the exponential Shannon and inverse Simpson indices (Figure 3b,c). In the Atlantic, Ecoregion 11 had the highest species richness and exponential Shannon index (Figure 3a,b). For the trawl/dredge/sledge data, highest diversity occurred in Ecoregions 4, 5 and 7 (q = 0 to 2) among the Arctic marine ecoregions (Figure 3d-f). These regional differences also dissipated with increased weighting on common and dominant species (q = 1 to 2, Figure 3e,f). Within the Atlantic (primarily the Gulf of St. Lawrence), we found consistently higher species diversity in Ecoregion 12 than Ecoregion 11 (Figure 3d-f).
For the quadrat data, Ecoregion 2 had the highest Pacific species richness ( Figure 3g); however, differences among ecoregions appeared insignificant for exponential Shannon (Figure 3h) and inverse Simpson indices (Figure 3i). For the Atlantic, the exponential Shannon and inverse Simpson of the quadrat data (Figure 3h,i) were higher in Ecoregion 12 than in Ecoregion 13.

| Bathymetric and latitudinal variations
For core/grab diversity calculated from 100 randomly selected individuals (by an average of 1,000 bootstrap resampling), the best TA B L E 1 Top 3 best regression models (LME) for quantitative core/grab data and simple linear regression models (LM) for quantitative trawl/sledge/dredge data Note: Hill numbers of order q = 0 (species richness), q = 1 (exponential Shannon index) and q = 2 (inverse Simpson index) were based on m = 100 randomly selected individuals. Bold font indicates significant parameter coefficient at p < .05.

F I G U R E 2
Species richness (q = 0), exponential Shannon index (q = 1) and inverse Simpson index (q = 2) accumulation curves across Canada's Pacific, Arctic and Atlantic Oceans for benthic samples collected with core/ grab (a-c), trawl/sledge/dredge (d-f) and quadrat (g-i  Table 1a). The best LME predicts that all three species diversity indices increased from intertidal to a maximum diversity at ~ 100 m depth near the shelf break and then declined towards bathyal depths; however, the effects of depth on diversity decline with increasing weight on common and dominant species (i.e. from richness to inverse Simpson; Figure 4a-c). Among the ecoregions, the unimodal diversity-depth relationships were most apparent in Ecoregions 2, 4, 7, 9 and 11 ( Figure S5), probably due to low data coverage in other ecoregions. The LME that explained the higher total variance when considering the random factor (or ecoregion) explained 38%, 23% and 14% of variance for species richness, exponential Shannon and inverse Simpson indices, respectively (R 2 ecoregion , Table 1a). The model averaging all possible LMEs also confirmed hump-shaped depth-diversity relationships with 100% likelihood of including depth in the average models (RI = 1, Table 2).
We found no statistical relationship between latitude and any of the three species diversity indices for either of the top three best LME models (Table 1a) or the average model (Table 2).
For the limited trawl/dredge/sledge data from the Arctic, we found a significant unimodal relationship with depth only for species richness based on 100 randomly selected individuals (Table 1b, Figure 4d). We also found no statistical relationship between depth and exponential Shannon or inverse Simpson indices (Figure 4e,f).

| Environmental correlates with diversity
Pearson correlation coefficients among water depth, latitude and the extracted environmental variables varied from −0.5 to +0.5.
For limited core/grab data (below the euphotic depth and with complementary environmental values), export POC flux was the only environmental variable (excluding depth and latitude) shared in the three best LME models in explaining species diversity indices from 100 randomly selected individuals (Table 3a). Predictive power was lower than for the bathymetric-latitudinal model (3%-14%, Table 3). Likewise, among all environmental variables, average LME models were most likely (87%-92%) to include export POC flux (Table 4a), which had a significantly negative effect on all three diversity indices ( Figure 5). In addition, species richness in the average LME model decreased with increasing seasonal variability in primary production (SVI); however, with low relative importance (RI) compared with export POC flux, whereas the exponential Shannon index varied negatively but slightly concaveupward with export POC flux (Table 4a). The non-significant temperature effect in the average model (Table 4a) nonetheless F I G U R E 3 Mean species richness (q = 0), exponential Shannon index (q = 1) and inverse Simpson index (q = 2) from 200 randomly selected samples (based on 1,000 permutations) across 13 marine ecoregions for benthic samples collected with core/grab (a-c), trawl/sledge/dredge (d-f) and quadrat (g-i). Error bar shows the 95% confidence intervals. The numbers of samples, sites and taxa in each category can be found in Table S2. The fill colours indicate temperate Pacific, Arctic and temperate Atlantic Oceans  (Table 3a). Interestingly, latitude played a weak role in the analysis including the environmental parameters.
For the limited trawl/sledge/dredge data (below euphotic depth and with complementary environmental data), only the quadratic term of export POC flux significantly affected the inverse Simpson index (Tables 3b and 4b).

| Benthic biodiversity predictions
We conducted predictive biodiversity mapping only with core/ grab data, because other types of samples did not have sufficient coverage across the three oceans. We input ETOP1 global relief model (i.e. water depth >0) and 13 ecoregions into the best LME models (see parameter coefficients in Table 1a; i.e. only depthrelated parameters) to estimate relative species richness, exponential Shannon and inverse Simpson indices ( Figure 6). More complex models based on full environmental variables (Table 3a) performed no better than the simple model (based on R 2 , cf. Tables 1a and 3a). In addition, our analyses did not predict intertidal soft-bottom benthic diversity because our input depth started below sea level. Despite the difficulty in mapping variation, average predictions (i.e. not considering the random factor ecoregion) for relative diversity increase away from shore to the continental shelf break (Figure 6a-c). The huge expanses of continental shelves in the Canadian Arctic and Atlantic support large areas of high predictive diversity. Beyond the shelf breaks (~100-m depth), predictive diversity decreases rapidly along the continental slope.
Considering the random factor (i.e. ecoregion) increased model prediction accuracy (higher R 2 , Table 1a) but with discontinuous predictive diversity around ecoregion boundaries (Figure 6d-f) Note: Hill numbers of order q = 0 (species richness), q = 1 (exponential Shannon index) and q = 2 (inverse Simpson index) were calculated from 100 randomly selected individuals. Bold font indicates significant coefficient at p < .05.

F I G U R E 4
Species richness (q = 0), exponential Shannon index (q = 1) and inverse Simpson index (q = 2) as functions of depth from the quantitative core/grab (a-c) and trawl/sledge/dredge samples (df). The Hill numbers (q = 0, 1 and 2) were calculated from 100 randomly selected individuals with 1,000 permutations. The solid circles indicate the rarefied Hill numbers, whereas the open circles indicate extrapolated Hill numbers. The solid line and shaded area indicate 1,000 bootstrap mean and 95% confidence intervals of the best linear mixed-effect (LME) model fit, respectively Note: Hill numbers of order q = 0 (species richness), q = 1 (exponential Shannon index) and q = 2 (inverse Simpson index) were based on m = 100 randomly selected individuals. Bold font indicates significant parameter coefficient at p < .05. Significance codes: *p < .05, **p < .01, ***p < .001.
when weighing more on the common and dominant species and thus considering species evenness (Figure 6e,f).

| D ISCUSS I ON
The TOBB database compilation increased the previous report of Canadian marine benthic infauna taxon lists from 2,127 (Archambault et al., 2010) to 2,829 taxa (core/grab data in this study). Adding the benthic epifauna (from trawl/sledge/dredge, quadrat and camera data) increases the total taxon list (3,337 taxa) more than 50% above the previous estimate. However, this number nonetheless clearly underestimates the total numbers of taxa present even within the existing database, because (a) only ~70% of taxa were identified to species; (b) for consistency, we only considered species names reported by WORMS (with accepted AphiaID); and (c) huge depth ranges and geographic areas of the three oceans remain un-sampled. Therefore, our compilation represents a minimum and very conservative estimate. Further, ongoing efforts to sample and identify marine invertebrates in Canadian waters have already increased the total number of taxa reported here (e.g. Arctic sponges: Murillo et al., 2018).

| Latitudinal variation
Although several hypotheses predict a decline in diversity with increasing latitude ( Note: Hill numbers of order q = 0 (species richness), q = 1 (exponential Shannon index) and q = 2 (inverse Simpson index) were calculated on m = 100 randomly selected individuals. Bold font indicates significant coefficient at p < .05.

TA B L E 4
Model average parameter estimate for selected quantitative core/ grab and trawl/sledge/dredge data F I G U R E 5 Partial residual plots showing the relationship between food supply (i.e. export POC flux) and infaunal diversity after taking into account other environmental factors in the best linear mixed-effects (from Table 3a). Panels show the partial residuals of (a) species richness (q = 0), (b) exponential Shannon (q = 1) and (c) inverse Simpson (q = 2) indices. The solid line and shaded area indicate the linear regression fit and its 95% confidence intervals between partial residuals and export POC flux, respectively F I G U R E 6 (a-c) Average and (d-f) ecoregion-level prediction of relative species richness, exponential Shannon and inverse Simpson indices of quantitative core/grab data. The grey solid lines indicate 100-m isobath. We predicted the infaunal diversity by inserting depth (ETOPO1 global relief model) and 13 marine ecoregions (numbers within dashline polygons) into the best linear mixedeffect model (Table 1). The predictions are then scaled to between zero and one Jablonski et al., 2006;McClain & Schlacher, 2015;Rex et al., 1993;Tittensor et al., 2010), we did not find higher diversity accumulation curves in the Canadian Pacific or Atlantic than in the Arctic  Spalding et al., 2007). Nevertheless, evidence from sample-based and individual-based analyses suggests that we previously underestimated Canadian Arctic benthic diversity, because low sampling intensity may have biased the compilation of previous taxonomic lists (Archambault et al., 2010). The Arctic covers the greatest total area and the largest number of ecoregions among Canada's three oceans, encompassing high complexity in habitat type and physical oceanography (Archambault et al., 2010;. These characteristics are consistent with our assertion that, of Canada's three oceans, the Arctic harbours significantly greater benthic diversity. Nevertheless, we acknowledge such patterns may not apply to all taxa (e.g. cold-water corals, sponges) and can be affected by different delineations of Arctic marine bioregions.

| Bathymetric variation
Despite widespread recognition of a hump-shaped diversity-depth relationship in deep-sea communities (Levin et al., 2001;McClain & Schlacher, 2015;Rex & Etter, 2010), amplitudes and depths of peak diversity vary among oceanic basins depending on the productivity of the overlying water (Rex & Etter, 2010). In the productive North Atlantic, benthic diversities peaked between 2,000 and 3,000 m; in contrast, the diversity maximum in the oligotrophic Gulf of Mexico and the Mediterranean Sea appears shallower at around 1,000 to 2,000 m (Rex & Etter, 2010;Wei & Rowe, 2019 Labrador cover a wide bathymetric range ( Figure S5). We, therefore, set the ecoregions as a random factor in our modelling to account for regional variations and examine the average effect; however, it is still possible that our observed bathymetric diversity patterns are influenced by the ecoregions with higher data coverage. Availability of more deep-water data (>1,000 m) in more ecoregions might well shift the shape and depth of peak diversity into deeper water. The vast unknown diversity in the deep-sea (Appeltans et al., 2012;Grassle & Maciolek, 1992;Mora, Tittensor, Adl, Simpson, & Worm, 2011) will also likely raise the amplitude of the diversity peak. Moreover, the shallower maximum depths of the Canadian Arctic relative to the Atlantic and Pacific likely affect the shape of the diversity-depth relationship through the boundary (or mid-domain) effect (Pineda & Caswell, 1998). Besides the global effect of bathymetry, benthic biodiversity prediction will be enhanced by regionalization of models.

| Diversity-energy relationship
Export POC flux, a proxy for food supply to the seafloor, appeared to be the dominant environmental factor controlling all species diversity measures in the core/grab data (Tables 3a and 4a) and the inverse Simpson index of the trawl/sledge/dredge data (Tables 3b   and 4b), albeit without improving diversity predictions. This variable must be interpreted with caution given the challenges of inferring POC flux in nearshore waters (Behrenfeld & Falkowski, 1997;Lutz et al., 2007), ice-covered regions (Arrigo et al., 2012) and locations with high terrestrial POC inputs (e.g. Pacific west coast and Arctic shelf) (Arrigo, Matrai, & Dijken, 2011). In addition, Lutz's POC flux and bottom temperature represent a long-term ocean condition not synced with the specific timing of faunal data collection at various seasons and years. We, therefore, offer our interpretation below as a hypothesis for future testing rather than as a definitive conclusion. In contrast, the influence of temperature and seasonality of primary production were only evident for core/grab data, with reduced relative importance compared to export POC flux. The higher relative importance of chemical (i.e. food supply) versus thermal energy (i.e. temperature) on seafloor biodiversity confirms recent studies on seafloor diversity-energy relationships (McClain et al., 2012;Tittensor et al., 2011;Woolley et al., 2016). However, in contrast to these reports of a unimodal POC flux-diversity relationship, we found a negative relationship between food flux and seafloor biodiversity. The key distinction is that our analysis covers a much shallower depth gradient on the continental shelf and upper slope, and thus likely reflects the negative portion of the unimodal relationship (McClain & Schlacher, 2015). The POC flux algorithm also removed intertidal and subtidal data from our analysis (see Methods).
Several hypotheses could explain the negative relationship between food flux and diversity, including competitive exclusion under high food supply (Huston, 1979;Rex & Etter, 2010), increasing megafaunal bioturbation (associated with increasing productivity) and suppressing infaunal species diversity (McClain & Barry, 2010), or high productivity leading to limitations in other resources such as oxygen or biogenic habitat (Tilman, 1982). Moreover, the negative relationship between seasonality of surface production and species richness in the core/grab data suggests that disturbance (
Moreover, the creation of such a database also identifies key knowledge gaps, and significant challenges for building evidence-based assessments for conservation planning and ongoing monitoring of protected areas, for example, in continental slope and Arctic intertidal environments. The lack of quantitative epifaunal data anywhere in the Arctic except for sparse data in a few of its ecoregions points to a need for better coverage. State-of-the-art statistical approaches can overcome uneven sampling and extrapolate diversity patterns to under-sampled regions; however, future research activities must fill these knowledge gaps in order to improve the accuracy and resolution of biodiversity pattern predictions. Having in mind the very large areas, the forecasted environmental changes across a vast array of latitudes (about 37 degrees of latitude differences in the presented data), providing new comparative biodiversity values from the Canadian seafloor with a mechanistic approach using chemical and thermal energy data, are valuable for the prediction of future diversity distribution. In that regard, our study marks the first steps in constructing the most comprehensive, quantitative seafloor biodiversity database in Canada to date and demonstrates its potential for benthic biodiversity predictions and substantial information for conservation management. More thorough assessment and understanding of Canada's three oceans of seafloor biodiversity will require continued efforts to acquire and synthesize existing or new benthic diversity datasets, including attention to develop integrated, and nationally consistent benthic survey approaches. Given the extent of Canada's territorial waters and potential extent of the national set of conservation networks, development and implementation of a coherent set of benthic monitoring approaches will require involvement from DFO (the principal regulatory authority), as well as from various stakeholders, including other federal departments, provincial/municipal government agencies, Indigenous groups, academia and community groups.

ACK N OWLED G EM ENTS
This study was supported by the Natural Sciences and Engineering

Research Council of Canada (NSERC) Canadian Healthy Oceans
Network. Additional support was provided by the Ministry of Science and Technology of Taiwan (MOST 108-2611-M-002-001) to CLW and NSERC Discovery Grants to MC and PS. Data sources were numerous and many organizations contributed information (see Table S1). We thank Kevin Conley, Mylène Bourque, Julie Lemieux, Mélanie Lévesque, Adeline Piot and Annie Séguin for contributing datasets (see Table S1). We thank the two anonymous referees for fruitful comments. Utility of large regional database for understanding abundance and diversity characteristics of natural marine soft substrate fauna.