Samples were collected along the Atlantic continental margins of Canada and the United States during a research cruise on board the NOAA research vessel Henry B. Bigelow (June 2017). We collected sediment push cores (i.d. = 6.7 cm, L = 35 cm) using the Remotely Operated Vehicle (ROV) ROPOS¹. We sampled 10 different stations (Figure 1 and Table 1) spanning a depth range of 230–996 m and encompassing continental shelf (Western Jordan Basin), continental slope (Fiddler’s Cove and outside Georges Canyon), submarine canyon (Corsair Canyon and Georges Canyon), and inter-canyon (Munson-Nygren Inter-canyon and Nygren-Heezen Inter-canyon) sedimentary environments. We named each station using the location’s followed by the habitat’s abbreviation (Sh for continental shelf; Sl for continental slope; C for canyons; I for inter-canyons), as shown in Table 1. At each station, we collected 4–7 push cores at randomly selected locations 10s of meters away from each other. We dedicated 1–2 cores at every station to analysis of sediment properties and the remaining cores (3–5) to evaluation of benthic nutrient fluxes and macrofaunal diversity. At every station, we also collected bottom water samples using the (2–4) Niskin bottles mounted on the ROV for nutrient analysis and for water exchange during incubations. The ROPOS CTD Seabird 19plus mounted on the ROV recorded depth, bottom water temperature, salinity, and dissolved oxygen, and we also estimated the shortest distance from the shore for each station.

We collected sediment for analysis of organic matter and grain size from the 0–2 cm top layer of the 1–2 dedicated cores at each station, homogenizing the sediment and then placing it in Whirl-Pak bags prior to storage in the dark at −20°C until analyzed. In this study, we use Total Organic Matter (TOM) and Total Organic Carbon (TOC) as a measure of food quantity. We use total nitrogen (TN) and total organic carbon to total nitrogen ratio (C: N) as a measure of organic matter quality over long time scales, with higher TN and lower C: N indicating fresher and higher quality organic matter (Godbold and Solan, 2009; Le Guitton et al., 2015; Campanyà-Llovet et al., 2017). We use concentrations of chlorophyll a (Chl a), phaeopigments (Phaeo), total pigments (Tot Pigm), as well as chlorophyll a to phaeopigments ratio (Chl a: Phaeo), and chlorophyll a to total organic carbon ratio (Chl a: TOC) as measures of phytodetritus input to the seafloor and short-term organic matter quality and freshness. In this sense higher concentrations of pigments and higher Chl a: Phaeo and Chl a: TOC indicate higher inputs of fresh phytodetritus to the seafloor (Pusceddu et al., 2009; Le Guitton et al., 2015). We use % gravel, % sand, % silt, and % clay, as well as mean grain size of the sortable silt fraction (MGS), as a measure of sediment particle size and distribution.

Sediment total organic matter (TOM) was calculated as the difference between dry (desiccated at 60°C for 24 h) and calcinated (muffle furnace at 450°C for 4 h) weight, and expressed as mg⋅g DW^–1 (Danovaro, 2010).

Total organic carbon (TOC) and total nitrogen (TN) were determined by drying a sediment subsample of 1–5 g (wet weight) at 60°C for 24 h, grinding it to a fine powder, and then weighing and acidifying (with pure HCl fumes) for 24 h to eliminate inorganic carbon. Samples were dried again at 60°C for 24 h before starting the analysis. We then weighed an aliquot of dried decarbonated sediments (15 mg) and folded it tightly into a tin capsule. A Carlo Erba NA1500 Series II elemental analyzer (EA) determined the sediment concentration of TOC and TN, expressed as mg⋅g DW^–1.

Sedimentary concentrations of chloroplastic pigments (chlorophyll a and phaeopigments) were determined using a spectrophotometer following Danovaro (2010). Pigments were extracted with 90% acetone (24 h in the dark at 4°C). After centrifugation (800 × g), the supernatant was used to determine the functional chlorophyll-a and acidified with 0.1 N HCl to estimate the amount of phaeopigments. We then dried the sediment at 60°C for 24 h prior to weighing. Sediment concentrations of pigments were expressed as μg⋅g DW^–1. Total phytopigment concentrations were defined as the sum of chlorophyll a and phaeopigment concentrations (Pusceddu et al., 2009).

We digested a subsample of sediment with hydrogen peroxide to eliminate any organic material present and then freeze-dried sediments before analysis with the Beckman Coulter LD13-320 laser diffraction analyzer to determine granulometric properties. Sieving was performed prior to analysis to ensure the elimination of large particles (gravel fraction). For each sample, we determined % of gravel, sand, silt and clay. We also determined mean grain size of the sortable silt fraction (MGS, μm). Sediments were classified following the sediment classification scheme based on the percentages of sand, silt, and clay (Shepard, 1954).

To evaluate fluxes of nitrate, nitrite, phosphate, ammonium and silicate, we incubated sediment cores (sediment volume: 527.6 ± 98.4 cm³; water volume: 705.8 ± 98.4 cm³) and overlying water for approximately 48 h and removed water samples for dissolved inorganic nutrients analysis at regular intervals during the incubation. A total of 32 incubations were run on board during the 2-week cruise.

After collection, sediment cores were acclimated for about 12 h, allowing sediment particles in suspension to settle back to the sediment surface. Several hours before the beginning of the experiment, the overlying water was carefully (without resuspending the sediment) exchanged with fresh, oxygenated bottom seawater collected in situ, allowing it to overflow in the surrounding water bath. This addition prevented hypoxia in the cores and removed toxic metabolites produced by community metabolism. To start the incubation, all visible bubbles were removed from the surface, chambers were sealed with acrylic caps equipped with magnetic stirrers, and the sediment cores were incubated in a refrigerator kept at in situ temperature (∼4.5°C) and in the dark for about 48 h. Stirrers were working for the duration of the experiment at approximately 3 revolutions per minute to avoid hypoxia near the sediment, without resuspending the sediment. Each 48-h experiment comprised three ∼ 12-h sequential incubation segments. At the beginning of every incubation, we extracted ∼30 ml of water using a 60-ml acid-rinsed plastic syringe. We used ∼ 5 ml of water to rinse the syringe and sample bottle and then removed ∼25 ml of water to store in 30-ml acid-rinsed HDPE plastic bottles at −80°C in upright position for subsequent analysis of dissolved inorganic nutrients. At the end of every 12-h incubation segment, we removed the lids and resampled water for nutrient analysis. During the following 1–6 h, overlying water in the chamber was carefully exchanged with fresh, oxygenated bottom seawater collected at each station to prevent hypoxia in the chambers and to remove toxic metabolites produced by community metabolism.

The concentrations of dissolved inorganic nutrients (nitrate, nitrite, ammonium, phosphate, and silicate) in the water sampled from the incubations as well from the bottom water samples, were determined using a continuous segmented flow analyzer (Seal AutoAnalyzer 3) at the Bedford Institute of Oceanography (Dartmouth, NS, Canada). Analyses were performed following the Industrial Method 186-72W for silicates, the Industrial Method 158-71W [adapted from Armstrong et al. (1967); Grasshoff (1969)] for nitrate and nitrite, the Industrial Method 155-71W [adapted from Murphy and Riley (1962); Aoyama et al. (2012)] for orthophosphate and the fluorometric method developed by Aminot and Kérouel (1997) for ammonium. Nutrient fluxes, expressed as μmol⋅m^–2⋅d^–1, were determined from the measured concentration changes in the overlying water as a function of time, water volume, and sediment area, summed over the three ∼12-h incubation segments.

At the end of each ∼48-h incubation, we removed all cores from the water bath and immediately processed them for subsequent analysis. Following extrusion of the sediment from the cores and sectioning into 0–5 and 5–10 cm sediment layers using inert plastic spatulas, we fixed the unsieved sediment in 4% buffered formaldehyde seawater in 500-ml plastic jars for later faunal identification. In the laboratory, we processed samples over a 300-μm sieve prior to subsequent transfer to 70% ethanol until we could complete microscopic analysis. Before identification, samples were stained with a few drops of Rose Bengal (0.5 g⋅L^–1) to facilitate sorting of the samples and identification of organisms.

For each sample, we sorted macrofaunal organisms and assessed abundances of the major taxa (Classes: Polychaeta, Oligochaeta, Amphipoda, Isopoda, Copepoda, Ostracoda, Bivalvia, Gastropoda, Scaphopoda, Sipunculidea, Ophiuroidea, Asteroidea, Holothuroidea, Echinoidea, and Subclasses: Hexacorallia, Octocorallia). We further identified polychaetes to the family level because they represented the most abundant taxon in the majority of our samples. In addition, the high taxonomic and functional diversity of polychaetes make them good indicators of environmental quality, as well as effective surrogates for total biodiversity in ecological studies (Pocklington and Wells, 1992; Dauvin et al., 2003; Olsgard et al., 2003). The large number of samples and the time-consuming nature of taxonomic analysis precluded species identification, however, we justify our decision to assess polychaete diversity at the family level based on previous studies that demonstrated the efficacy of this approach as a valid alternative to species level analysis when investigating patterns of community structure, functional diversity, species distribution and effects of environmental variables on biological communities (Fauchald and Jumars, 1979; Jumars et al., 2015; Checon and Amaral, 2017). We calculated total macrofaunal densities (ind⋅m^–2), total number of taxa (including classes, subclasses and polychaete families), and Pielou’s evenness (J′) for each sample combining the data from the entire 10-cm cores. Diversity indices were calculated in PRIMER 6 + using the DIVERSE routine.

In order to evaluate biological trait expression, we used Biological Trait Analysis (Bremner et al., 2003), which uses multivariate ordination to describe patterns of biological trait composition over the entire macrofaunal assemblage and quantifies the types of trait present in assemblages and the relative frequency with which they occur, thereby providing a means to explore patterns in assemblage functional structure and functioning (Bremner et al., 2006). We selected 6 biological traits related to different aspects of life histories and ecosystem functioning, and subdivided the 6 traits into 27 categories that characterized behavior/strategies in more detail (Table 2). We used a fuzzy coding approach to assign trait categories to the taxa (classes, subclasses, or families) that allowed each taxon to represent more than one trait, therefore capturing inter- and intraspecific variation in trait expression. We adopted a scoring range of 0 to 5, with 0 reflecting no affinity for the given trait category, 1, 2, 3, or 4 reflecting partial, increasing affinity, and 5 denoting exclusive affinity. We derived information on trait expression for all the taxa from several published sources (Fauchald and Jumars, 1979; Highsmith and Coyle, 1991; Hyne, 2011; Queirós et al., 2013; Jumars et al., 2015; Polytraits Team, 2020) as well as from direct observations on our specimens. When information was unavailable for a given taxon, we obtained information from one taxonomic rank higher (e.g., from orders of polychaetes). In a few cases where no information at all was available for a given taxon, we distributed the 5 scores equally among all the plausible trait categories. To obtain the community trait expression in each sample, we multiplied trait categories for each taxon present in a sample by its density (ind⋅m^–2) in that sample, and then summed over all taxa present in each core to obtain a single value for each trait category in each sample (Bremner et al., 2006). We then explored the resulting matrix using multivariate analysis.

The opportunistic nature of our sampling resulted in an unbalanced design with different habitats represented differently at different depths (e.g., most shelf and slope stations were at shallower depths and most inter-canyon stations were at deeper depths). This confounding complicated the comparison between habitats because of well-established effects of depth (and its correlates) on biological communities and functional processes (Rex et al., 2006). To address this concern, we ran a preliminary analysis to understand the effect of depth on macrofaunal community structure (chosen because community structure showed the most obvious differences among stations) among all the habitats. To do so, we assigned depth classes to each station as follows: WJBSh 200 m; OGCSl and GC 600 m; FCoSl, MNI-1, and NHI-1 800 m; MNI-2 and NHI-2 900 m; CC and MNI-3 1000 m. We then assessed variation in macrofaunal community composition among depth classes using multivariate analysis of variance (PERMANOVA) on Bray-Curtis similarity matrices of square-root transformed density data and we ran the pair-wise comparisons as post hoc analyses to identify which depth classes differed in terms of macrofaunal composition. Based on these results, we defined four a posteriori depth classes (Table 3) that we used to investigate variation in benthic nutrient fluxes and macrofauna across habitats. We also note the unbalanced distribution of our stations, with clear geographic separation of WJBSh from all other stations, and the clustering together of all inter-canyon stations apart from other stations (see Figure 1). We acknowledge that this distribution represents a limitation of this study and might have affected our findings.

To investigate variation in each benthic nutrient flux (nitrate, nitrite, phosphate, silicate, and ammonium), taxonomic diversity index (total macrofaunal density, total number of taxa, and Pielou’s evenness), and vertical distribution of macrofauna (% organisms in the 0–5 cm layer) among habitats we used two-way type III univariate analysis of variance (ANOVA). We ensured that our data met homogeneity of variance (with Levene’s tests) and normality (with Q-Q plots of residuals) assumptions prior to analysis. We then performed post hoc pair-wise comparisons of significant effects (p < 0.05) using standard Tukey’s tests. We also investigated variation in multivariate benthic nutrient fluxes, macrofaunal community composition, and biological trait expression (separately) among habitats using multivariate analysis of variance (PERMANOVA) performed with 9999 random permutations of appropriate units. We ran the pair-wise comparisons as post hoc analysis whenever we found significant differences (p < 0.05). We verified homogeneity of multivariate dispersions using the PERMDISP routine. For both ANOVA and PERMANOVA analyses we used a nested design with the fixed factors “depth” (3 levels), and “habitat (depth)” (2 levels within the 600 m depth class: slope, canyon; 3 levels within the 800–1000 m depth class: slope, inter-canyon, canyon).

We visualized separation of multivariate benthic nutrient fluxes in ordination space using Principal Component Analysis (PCA), and separation of macrofaunal community composition and biological trait expression (separately) in ordination space using non-metric multidimensional scaling (nMDS) ordinations of similarity matrices. We also identified the taxa and biological trait categories that distinguished assemblages among habitats using a percent similarity procedure (SIMPER) analysis (Clarke and Gorley, 2006).

We next ran a set of analyses to identify the environmental drivers of benthic nutrient fluxes and macrofauna. We first explored correlations between total macrofaunal density, total number of taxa, and Pielou’s evenness (separately) and all available environmental variables using Draftman’s plots and correlation analysis. We next used a stepwise distance-based linear model permutation test (DistLM; McArdle and Anderson, 2001) to identify which set of environmental variables predicted variation of multivariate benthic nutrient fluxes, macrofaunal community composition, and biological trait expression (in 3 separate analyses). We used resemblance matrixes of multivariate benthic nutrient flux data (based on Euclidean distances), and macrofaunal community composition and biological trait expression (based on Bray-Curtis similarity, calculated from square root transformed data) as a measure of between-samples similarities. The predictive environmental variables allowed to enter the models were: latitude, distance from shore, depth, bottom water temperature, salinity, oxygen concentration, bottom water concentrations of nitrate, silicate, phosphate ad ammonium, sedimentary concentrations of TOM, TOC, TN, Chl a, Phaeo, and Tot Pigm, C: N, Chl a: Phaeo, Chl a: TOC ratios, % gravel, % sand, % silt, % clay, and MGS. Variables were standardized to mean 0 and standard deviation 1 prior to analysis. We assessed normality and collinearity of predictor variables using Draftsman’s plots, ensuring that highly correlated variables did not appear simultaneously in the final models.

Finally, we identified the biological drivers of benthic nutrient fluxes in order to assess underlying biodiversity and ecosystem functioning (BEF) relationships. To do so, we used a stepwise distance-based linear model permutation test (DistLM; McArdle and Anderson, 2001), with macrofaunal taxonomic diversity indices, as well as community composition and biological trait expression matrices as predictive variables. The use of matrices of community composition and biological trait expression in this analysis allowed testing of the relevance of each taxon density and trait category expression, respectively, on multivariate flux variation. We used the resemblance matrix of multivariate benthic nutrient flux data based on Euclidean distances as a measure of between-samples similarities. Noting that the number of predictor variables in this case greatly exceeded the number of samples, we did preliminary testing of the influence of each group of biological variables (taxonomic diversity indices, macrofaunal community composition, and macrofaunal biological trait expression) on benthic nutrient fluxes separately. Taxa that appeared in fewer than three samples (echinoids, amphinomids, heterosporous, phyllodocids, serpulids) were removed from the community composition matrix to further reduce the number of predictor variables and the number of zeros. Taxonomic diversity indices were standardized to mean 0 and standard deviation 1, whereas community composition and biological trait expression data were square-root transformed prior to analysis. For each group, the variable(s) that correlated best with benthic nutrient fluxes were selected and combined in a final analysis to determine the best biological model explaining variation in benthic nutrient fluxes.

For all DistLM analyses, we ran stepwise routines with 9999 permutations and used AICc (Akaike’s information criterion corrected) selection criterion, which is recommended for analyses with a small number of samples relative to the number of predictor variables (Anderson et al., 2008). We examined R² to identify the best model and determine the proportion of the variation explained by that model, visualizing results with distance-based redundancy analysis (dbRDA; Anderson et al., 2008). All multivariate analyses, including PERMANOVA, PERMDISP, SIMPER, PCA, nMDS, DistLM, and dbRDA analyses were performed in PRIMER v6 (Anderson et al., 2008).

In terms of sedimentary organic matter quantity (Table 4) GC, CC, and WJBSh were characterized by the highest concentrations of TOM and TOC. In terms of sedimentary organic matter quality (Table 4), WJBSh was characterized by the highest long-term quality (e.g., high TN and low C: N), whereas GC had the highest input of phytodetritus (high Chl a, Chl a: TOC) of all stations. Other stations presented intermediate levels of organic matter quantity and quality, with the MNI stations showing the lowest input of phytodetritus.

In terms of granulometric properties (Supplementary Table 1), most sediments were classified as silty sand according to the Shepard (1954) classification scheme, except for sandy sediments at GC and MNI-2, sandy silt at MNI-1 and clay-silt at WJBSh. MGS ranged from 20 μm to 39 μm at WJBSh and NHI-1, respectively.

Nitrate flux (Figure 2A) ranged from −725.4 μmol⋅m^–2⋅d^–1 to 1288.4 μmol⋅m^–2⋅d^–1 across our replicate cores and the average values were directed toward the water column (release of nitrate) at all stations. Nitrite flux (Figure 2B) ranged from −245.7 μmol⋅m^–2⋅d^–1 to 62.0 μmol ⋅m^–2⋅d^–1 across our replicate cores and the average values were directed toward the sediment (uptake of nitrite) at all stations except MNI-2 and CC. Ammonium flux (Figure 2C) ranged from −2705.9 μmol⋅m^–2⋅d^–1 to 404.0 μmol⋅m^–2⋅d^–1 across our replicate cores and the average values were directed toward the sediment (or uptake of ammonium) at all stations except GC and FCoSl. Phosphate flux (Figure 2D) ranged from −108.9 μmol⋅m^–2⋅d^–1 to 778.3 μmol⋅m^–2⋅d^–1 across our replicate cores and the average values were directed toward the water column (release of phosphate) at all stations except FCoSl. Silicate flux (Figure 2E) ranged from −7225.4 μmol⋅m^–2 d^–1 to 6473.5 μmol⋅m^–2 d^–1 across our replicate cores and the average values were directed toward the water column (release of silicate) at all stations except FCoSl. We detected significant differences in silicate fluxes between habitats nested in depth classes (ANOVA, p < 0.05, Supplementary Table 2) and post hoc pair-wise analysis revealed significantly higher silicate effluxes in inter-canyon compared to slope habitats within the 800–1000 m depth class (Figure 2E). Other fluxes did not differ significantly among habitats nested in depth classes (ANOVA, p > 0.05, Supplementary Table 2).

Multivariate analysis combining all the fluxes provided comparison of the overall remineralization function of the sediments, for which we detected significant differences (PERMANOVA, p < 0.05, Supplementary Table 3) among habitats nested in depth classes. Post hoc tests revealed significantly different benthic nutrient fluxes between inter-canyon and slope habitats within the 800–1000 m depth class. In the PCA plot (Figure 3), the first two PC axes explained 68% of the variation of benthic nutrient fluxes and analysis of the eigenvectors showed that no single flux dominated the multivariate similarity pattern among samples, with nitrite flux correlating most strongly with the first PCA axis (positively), and silicate flux correlating most strongly with the second PCA axis (negatively).

We sorted 1,257 macrofaunal organisms in total, representing 40 different taxa. Densities varied from 2,554 to 27,526 ind⋅m^–2 at MNI-2 and WJBSh, respectively. Polychaetes were the dominant taxon in most samples (average 49% of total macrofauna), followed by amphipods (average 15% of total macrofauna). Among the polychaetes, the families Cirratulidae, Paraonidae, and Polynoidae dominated across stations, representing average 22, 12, and 9% of total polychaetes, respectively.

Total macrofaunal density (Figure 4A) varied significantly between habitats nested in depth classes (ANOVA, p < 0.05, Supplementary Table 4) and post hoc tests reveled higher density in canyon compared to slope habitat within the 600 m depth class. Total number of taxa (Figure 4B) and Pielou’s evenness (Figure 4C) did not vary significantly among habitats (ANOVA, p > 0.05, Supplementary Table 4).

The percentage of macrofaunal organisms in the upper 5-cm layer ranged from 83 to 100% and varied significantly (ANOVA, p < 0.05, Supplementary Table 4) between habitats nested in depth classes. Post hoc tests revealed significantly lower percentage of organisms in the upper sediment layer in slope compared to inter-canyon habitats (Figure 4D).

Macrofaunal community composition differed significantly (PERMANOVA, p < 0.05, Supplementary Table 5) among habitats nested in depth classes and post hoc tests revealed different macrofaunal composition between slope and canyon habitats within the 600 m depth class, and between inter-canyon and other habitats within the 800–1000 m depth class. Stations and habitats also separated in ordination space based on our nMDS analysis (Figure 5A).

According to our SIMPER analysis, ophiuroids, amphipods, cossuridae, polynoids, ostracods, cirratulids, and bivalves were the taxa that mostly contributed to differences between slope and canyon habitats within the 600 m depth class, and all were more abundant in the canyon, whereas opheliids were more abundant on the slope. The average dissimilarity between Georges Canyon and adjacent slope was 43% and these taxa contributed to 52% of the variation. Bivalves, cirratulids, eunicids, maldanids, ostracods, lumbrinerids, amphipods, and hesionids contributed most to differences between inter-canyons and other habitats within the 800–1000 m depth class were more abundant in slope and canyon habitats, whereas ophiuroids, paraonids, and nereids were more abundant in inter-canyons. The average dissimilarity between inter-canyon and other habitats was 62% and these taxa contributed 52% of the variation.

In terms of biological traits, motile, deposit-feeders, surface modifier organisms with indirect, pelagic larval development, and medium lifespan dominated sediments overall. Macrofaunal community biological trait expression differed significantly (PERMANOVA, p < 0.05, Supplementary Table 6) among habitats nested in depth classes and post hoc tests revealed different trait expression between slope and canyon habitats within the 600 m depth class, and between inter-canyon and slope habitats within the 800–1000 m depth class. Stations and habitats also separated in ordination space based on our nMDS analysis (Figure 5B).

According to our SIMPER analysis, greater expression of motile, pelagic, indirect and direct larval development, high and very high fecundity, surface modifiers and lifespan 3–5 years in canyon habitats contributed most to trait modality differences between Georges Canyon and the adjacent slope. The average dissimilarity between slope and canyon was 18% and these trait modalities contributed 50% of the variation. Greater expression of indirect, pelagic larval development, sessile, discretely motile and motile, suspension and filter feeder, upward and downward conveyors, lifespan 3–5 years, and medium fecundity in slope habitats contributed most to trait modality differences between inter-canyon and slope habitats within the 800–1000 m depth class. The average dissimilarity between inter-canyon and slope was 20% and these trait modalities contributed 51% of the variation.

The best environmental distance-based linear model (DistLM), explaining 39% of the overall variation in benthic nutrient fluxes, included the variables bottom water oxygen concentration, distance from shore, and percentage of gravel (Table 5 and Supplementary Figure 1a). The best biological distance-based linear model (DistLM) based on all biological variables, explained 49% of the overall variation in benthic nutrient fluxes, and included the relative densities of orbinid, onuphid, and maldanid polychaetes and sipunculids (Table 5 and Supplementary Figure 1b).

Draftsman plots and correlation analysis identified sedimentary concentration of chlorophyll a as the single variable that best correlated with total macrofaunal density (R = 0.85, Figure 6), total number of taxa (R = 0.47), and Pielou’s evenness (R = −0.56).

The best environmental distance-based linear model (DistLM), which explained 35% of the overall variation in macrofaunal community composition, included the variables latitude, Chl a: TOC ratio, sedimentary concentration of TOM, and depth (Table 6 and Supplementary Figure 2a). The best environmental distance-based linear model (DistLM) explained 72% of the overall variation in macrofaunal biological trait expression and included the variables sedimentary concentration of Chl a and TN, latitude, and sedimentary concentration of Phaeo (Table 6 and Supplementary Figure 2b).

Benthic nutrient fluxes did not show clear variation patterns among the different habitats sampled. Fiddler’s Cove, on the continental slope, was the station that differentiated most from the others, mostly because of the higher intake and/or lower release of phosphate and silicate. The source of such differences remains open to investigation as this station did not differentiate from others in terms of physico-chemical variables, granulometric properties, sedimentary organic matter, or macrofaunal diversity. At the same time, sediments in Georges Canyon, which were characterized by distinct sedimentary organic matter composition and macrofaunal communities, did not differ from others in terms of benthic nutrient fluxes, pointing to different mechanisms shaping biological communities and benthic remineralization processes. The complexity of biogeochemical processes occurring in marine sediments also complicates understanding them and our ability to identify their patterns clearly (Hall et al., 1996). For example, up to 75% of the ammonium formed during mineralization of sedimentary organic matter is nitrified and subsequently denitrified in the sediment, consistent with the absence of significant ammonium in- or effluxes across the sediment-water interface during intense organic matter remineralization (Devol and Christensen, 1993; Enoksson, 1993). A similar decoupling of remineralization and nutrient fluxes has been reported for other nutrients (Hall et al., 1996).

Percentage of gravel, bottom water oxygen concentration, and distance from shore were the best environmental predictors of multivariate benthic nutrient fluxes. The effect of presence of gravel in the sediments (detected in some of the samples from WJBSh and NHI-1) on flux rates might be related to the higher permeability of gravel that affects solute transport (e.g., dominance of porewater advection in gravel and sand and dominance of molecular diffusion in muddy sediments; Janssen et al., 2005). Bottom water oxygen concentration was found to correlate positively with fluxes of nitrate, nitrite, and ammonium, and negatively with fluxes of silicate and phosphate. The higher release of phosphate into the water column in lower oxygen environments (e.g., OGCSl) can be explained by sedimentary redox conditions and depletion of ferric iron that diminish the capacity of sediments to retain phosphate (Ingall and Jahnke, 1994; Paytan and McLaughlin, 2007). Finally, we observed higher nutrient fluxes offshore, which contrasts with other studies that detected higher fluxes nearshore on the continental shelf (e.g., Friedrich et al., 2002) and could relate to lower bottom water oxygen concentrations at our nearshore stations.

Surprisingly, variables related to sedimentary organic matter quantity and quality did not contribute substantially to variation of benthic nutrient fluxes, in contrast to other studies (e.g., Link et al., 2013a, b; Belley et al., 2016; Miatta and Snelgrove, 2021), indicating the overriding importance of other environmental factors in determining benthic flux variation among our stations. Importantly, topographic features, such as canyons and inter-canyons, influence sedimentary organic matter deposition, hydrodynamic forces, and disturbance events (e.g., sediment deposition and resuspension) and might affect benthic nutrient fluxes, decoupling organic matter input and its remineralization. For instance, previous studies reported that sediment resuspension influence oxygen and nutrient effluxes in marine sediments (Tengberg et al., 2003; Niemistö et al., 2018), and intense sedimentation rates in slope sediments were also found to alter benthic nutrient fluxes (Hensen et al., 2000). The absence of other studies measuring benthic nutrient fluxes in canyon and inter-canyon habitats limits our capacity to contrast our findings with other regions, and points to the need for more studies investigating benthic processes in these understudied environments. Other environmental factors not considered in our study might also contribute to the variation of benthic fluxes. These include, for instance, sediment concentrations of manganese and iron, linked to variation in benthic nutrient fluxes (Link et al., 2013b) through the sequestration of phosphate by ferric iron (Paytan and McLaughlin, 2007), as well as density and diversity of sedimentary bacteria (Belley and Snelgrove, 2016), as they are directly responsible for organic matter remineralization and nutrient regeneration in marine sediments (Jorgensen, 2006).

Overall, biological variables explained more variation in benthic nutrient fluxes than environmental variables, with relative densities of orbiniid, onuphid, and maldanid polychaetes and sipunculids among the best biological predictors. These taxa likely affect fluxes through their activities, including bioturbation. Orbiniid polychaetes and sipunculids contribute to bioturbation mostly by biodiffusion and up- and downward conveyor (Queirós et al., 2013), similarly to maldanid polychaetes, which are also considered funnel feeders and were reported as important contributors to benthic fluxes in a previous study (Belley and Snelgrove, 2016). Onuphids are tube-forming polychaetes that build thin, flimsy burrows made from sand grains embedded in a polysaccharide matrix (Waldbusser and Marinelli, 2009) that are permeable to diffusive exchange of solutes (Aller, 1983; Hannides et al., 2005). Despite a lack of available specific information, the characteristics of onuphid tubes and burrows suggest potential bioirrigation activities, which markedly increase nutrient exchange between pore water and the overlying water column (Kristensen and Andersen, 1987; Aller, 1988; Heilskov et al., 2006; Meysman et al., 2006). We observed heterogeneity in the effect of different organisms on each nutrient flux, underscoring the complexity of biodiversity and ecosystem functioning relationships, which is further complicated by interactions between organisms and environmental characteristics. For example, the effect of bioturbation activities on benthic nutrient fluxes is usually more pronounced in muddy than sandy sediments, where hydrological processes tend to determine porewater advection (Mermillod-Blondin, 2011). However, Waldbusser and Marinelli (2009) reported that, whereas environmental variables such as sediment granulometry and physical forces may drive large-scale variability in porewater advection in permeable sediments, type and abundance of bioturbating infauna significantly affect smaller-scale variation. The effect of bioturbation on benthic fluxes can also be solute-specific, further complicating its understanding. For instance, whereas bioirrigation generally increases nutrient fluxes, it can negatively affect silicate and phosphate effluxes (Waldbusser and Marinelli, 2009), potentially by increasing oxygenation of the sediment that increases the capacity of sediments to absorb and retain inorganic phosphorus (Ingall and Jahnke, 1994; Paytan and McLaughlin, 2007).

Macrofaunal density and taxonomic diversity indices correlated poorly with multivariate benthic nutrient fluxes, adding evidence for the greater importance of taxonomic and functional identity over community diversity in regulating ecosystem processes (Hooper et al., 2005). In our study, the lower importance of biological trait expression in determining variation of benthic nutrient fluxes compared to the relative densities of key taxa is likely a consequence of our analytical method. For instance, the low taxonomic resolution used in our study, the lack of data on trait expression for certain taxa, the exclusion of potentially important traits (e.g., bioirrigation, biodeposition) due to lack of information, the use of literature trait data rather than direct observation of individuals traits, and the use of relative densities rather than biomass to weight trait expression might reduce the accuracy of our trait analysis. Interestingly other studies reported a low correlation between bioturbation intensity and nutrient generation in sediments (Solan et al., 2008), suggesting that the relationships between ecosystem functioning and bioturbation might not be so straightforward and that different species might interact with the benthic environment differently, causing variability in how bioturbation influences benthic fluxes (Ieno et al., 2006; Solan et al., 2008).

In contrast to benthic nutrient fluxes, macrofauna displayed clear variation patterns among habitats, with community composition proving more sensitive than diversity indices and biological trait expression in characterizing patterns. Differences in macrofauna appeared to be driven mostly by the quantity and quality of sedimentary organic matter. In particular, the sedimentary concentration of chlorophyll a, indicative of input of fresh marine-derived, highly labile organic matter to the seafloor, was the best predictor of macrofaunal density and taxonomic diversity. These findings support the role of phytodetritus as an important source of nutrition for benthic organisms and the importance of food quality over quantity in sustaining benthic communities (Campanyà-Llovet et al., 2017; Leduc et al., 2020). Input of labile organic matter has been previously reported to sustain high macrofaunal biomass (Pilditch et al., 2015; Leduc et al., 2020), and higher biodiversity through the creation of new niches that enhance the coexistence of different taxa (Levin et al., 2001). Lower Pielou’s evenness associated with higher Chl a suggests a particular advantage for a few opportunistic taxa that respond to the pulses of highly labile food (Levin et al., 2001). Depth was less important overall than food availability in determining macrofaunal diversity, and we attribute the disruption of the typical depth-diversity relationship to the presence of lateral or down-slope transport of organic matter at some of our stations (Flach and Heip, 1996; Levin et al., 2001; Wei et al., 2010).

Input of phytodetritus, together with latitude, organic matter quantity and quality, and depth were the best environmental predictors of macrofaunal community composition and biological trait expression. These findings align with other studies (e.g., Levin and Gage, 1998; Wei et al., 2010; Robertson et al., 2020). For example, Wei et al. (2010) found that patterns of macrofaunal composition correlated strongly with depth and flux of particulate organic carbon from surface production, whereas Käß et al. (2021) reported a great effect of labile, phytodetrital organic matter in determining macrofaunal biological traits. Similarly, sedimentary concentration of Chl a alone explained nearly half of the total variation of biological trait expression among our stations. The low proportion of community composition variability explained by our model points to the important role of other factors in determining community assemblages. For instance, our study did not explicitly consider hydrodynamic patterns and physical disturbance, which are important drivers of benthic community composition along continental margins (Levin et al., 2001; McClain and Barry, 2010; Cunha et al., 2011; Robertson et al., 2020). The presence of sedimentary mega-epifauna (e.g., sea pens) can also influence macroinfaunal communities (Miatta and Snelgrove, submitted) and we observed differences in megaepifaunal density and diversity among stations during our ROV dives, such as the presence of sea pen fields at FCoSl (Figure 8D).

The shallowest sediments on the continental shelf at Western Jordan Basin (Figure 7, continental shelf), where sedimentary organic matter properties point to the accumulation of organic material (Campanyà-Llovet et al., 2018) that indeed typically characterizes sediments on the continental shelf, displayed relatively high macrofaunal density and diversity, with distinct taxa and trait composition (statistical results not reported here). Here, greater presence of predators and a reduced role for deposit feeders probably relate to higher environmental stability and long-term availability of organic matter that favors the occupation of multiple ecological niches and the presence of taxa with contrasting feeding strategies and higher trophic levels (Simboura et al., 2000). Higher presence of biodiffusers living closer to the sediment-water interface, as well as sessile mega-epifauna (author’s personal observation; Figure 8A) likely benefited from greater availability of particulate organic matter, coupled with lower hydrodynamic forces (Harris, 2014; Pierdomenico et al., 2019), as suggested by lower mean grain size at this station (Van Rijn, 1993).

The organic enrichment measured in the sediments from Georges Canyon (Figure 7, shelf-incising canyon), together with signs of deposition (e.g., presence of large detritus in some areas of the canyon, Figure 8C) and strong currents (e.g., presence of ripples, Figure 8C; higher mean grain size) point to potential occurrence of turbidity and tidal currents, sediment gravity flows and other lateral transport mechanisms (de Stigter et al., 2007; Puig et al., 2014). Reports from other canyons indicate similar findings (e.g., Pusceddu et al., 2010; Amaro et al., 2015; Campanyà-Llovet et al., 2018; Pierdomenico et al., 2019; Leduc et al., 2020). The high concentrations of phytopigments and high Chl a: Phaeo ratio in the sediments from GC point to rapid downward transport of marine-derived organic matter along the canyon axis, contributing to the supply of labile organic matter at depths far below the sinking of particulates from productive surface waters (Campanyà-Llovet et al., 2018). The presence of sea pigs (class Holothuroidea, order Elasipodida) on the seafloor in some parts of Georges Canyon (author’s personal observation, as well as on the outer slope (Figure 8B), also confirms recent input of fresh organic matter, because sea pigs are mobile deposit feeders that form transient, dense aggregations on the seafloor in response to high influx of fresh organic matter that they consume in a timespan of hours to days (Miller et al., 2000). The increased organic matter quantity and quality in Georges Canyon likely contributed to high macrofaunal density, low evenness and dominance of deposit feeders that characterized sediments in this canyon compared to the adjacent slope. These findings align with other studies (e.g., De Leo et al., 2010, 2014; Cunha et al., 2011). Organic enrichment usually contributes to high biomass and depressed biodiversity and evenness because it favors high densities of a small number of opportunistic species (Levin et al., 1991, 2001; Rosenzweig and Abramsky, 1993; Levin and Gage, 1998). Moreover, the occurrence of turbidity flows can affect biological communities through periodic physical disturbance (de Stigter et al., 2007; De Leo et al., 2014; Puig et al., 2014; Liao et al., 2017). For instance, we observed scarce presence of sedimentary mega-epifauna (Figure 8C), but occurrences of corals and sponges on the canyon’s walls and boulders, where they can take advantage of abundant flux of particulate organic matter with minimal disturbance through sediment deposition and resuspension (Harris, 2014; Pierdomenico et al., 2019). Moreover, macrofaunal traits such as dominance of mobile deposit feeders with high fecundity and shorter life span and presence of opportunistic taxa in the sediments from Georges Canyon, supports the effect of disturbance events. Periodic disturbance caused by turbidity flows creates repeated opportunities for recolonization and maintains the benthic fauna in an early successional state dominated by opportunists, with high abundance and biomass of organisms that can tolerate disturbance and exploit the increased food availability (Levin et al., 2001; McClain and Barry, 2010; Vetter et al., 2010). Disturbance can also help maintain high biodiversity by supporting small-scale heterogeneity that maintains habitat niches (Snelgrove, 1999; Levin et al., 2001), which likely contributed to the relatively high number of taxa in Georges Canyon. Our findings support the argument that gradients of productivity, sedimentary processes, and physical disturbance play an important role in shaping biological communities in canyon habitats (Levin et al., 2001; McClain and Barry, 2010; Robertson et al., 2020).

In contrast, macrofaunal communities from Corsair Canyon (Figure 7, blind canyon, Figure 8F), did not differentiate from other stations at comparable depth, suggesting that this canyon is unaffected by lateral transport mechanisms such as turbidity currents that influence biological communities through organic enrichment and physical disturbance. Previous studies report no evidence for turbidity currents or down-axis mass movement of sediment in Corsair Canyon (Dillon and Zimmerman, 1970). Granulometric properties with higher presence of mud in our sites in Corsair Canyon contrasting with sandy sediments in Georges Canyon, also suggests that Corsair Canyon is a slope canyon that does not indent the continental shelf edge (Jobe et al., 2011), therefore limiting the lateral transport of organic materials along this canyon (Harris and Whiteway, 2011). However, we cannot rule out that differences in depth, sediment type, or other variables could also contribute to differences between Georges and Corsair canyons [as reported in other systems by Williams et al. (2009); Robertson et al. (2020)]. We also acknowledge evidence of substantive small-scale heterogeneity in resource availability and biodiversity in canyons, not only between canyon axes but also within the same canyon axis at scales <100 m (McClain and Barry, 2010; Cunha et al., 2011; Campanyà-Llovet et al., 2018). The limited sampling opportunities in our study within each canyon constrains our ability to draw clear conclusions on the physical and ecological processes in Georges and Corsair Canyons at larger scales.

Finally, higher percentage of organisms in the upper sediment layers characterized inter-canyons (Figure 7, inter-canyon), along with differences in community composition and biological trait expression compared to slope sediments at comparable depth. Moreover, even though differences were not significant, density and diversity tended to be lower in inter-canyons (especially at Munson-Nygren Inter-canyon) than in slope sediments. These findings suggest the presence of different mechanisms influencing macrofaunal communities in inter-canyons, although the absence of comparative studies limits our capacity to draw strong inferences. Perhaps a combination of limited resource availability and intense bottom currents determined the observed patterns. Previous studies have reported negative effects of intense hydrodynamics on benthic communities (e.g., Levin et al., 2001) and in our study strong currents are suggested by the coarse mean grain size in inter-canyons (especially NHI), as well as the lower quantity and quality of sedimentary organic matter, that can result from decreased rates of organic matter deposition in relation to strong hydrodynamics (Vetter et al., 2010). We also observed very high densities and diversity of mega-epifauna (e.g., corals, sponges, and bivalves) on the inter-canyon walls and boulders, whereas sediments were generally bare (Figure 8E), as also reported by Quattrini et al. (2015), another indication of strong hydrodynamics that can disturb sedimentary mega-epifauna through sediment resuspension (Pierdomenico et al., 2019). Other findings, such as high C: N ratio values at some inter-canyon sites, point to the presence of non-canyon related lateral transport mechanism that conveys organic matter of terrestrial origin (Flach and Heip, 1996) from the continental shelf along inter-canyons. A previous study also documented the presence of unconsolidated sediments, and deposition of larger material, suggesting a highly geologically dynamic nature for inter-canyon areas (Quattrini et al., 2015). In our study, the high contribution of terrestrial material to the organic carbon pool in inter-canyons may also contribute to the depressed density and diversity of macrofauna (Cunha et al., 2011; Leduc et al., 2020), as suggested by Leduc et al. (2020) who recently identified marine organic matter as the main limiting factor shaping macrofaunal communities in New Zealand submarine canyons.