Distribution of Megabenthic Communities Under Contrasting Settings in Deep-Sea Cold Seeps Near Northwest Atlantic Canyons

Cold seeps support fragile deep-sea communities of high biodiversity and are often found in areas with high commercial interest. Protecting them from encroaching human impacts (bottom trawling, oil and gas exploitation, climate change) requires an advanced understanding of the drivers shaping their spatial distribution and biodiversity. Based on the analysis of 2,075 high-quality images from six remotely operated vehicle dives, we examined cold seep megabenthic community composition, richness, density, and biodiversity at a relatively shallow (∼400 m water depth) site near Baltimore Canyon (BC) and a much deeper site (∼1,500 m) near Norfolk Canyon (NC), in the northwest Atlantic. We found sharp differences in the megabenthic composition between the sites, which were driven mostly by bathymetric gradients. At both BC and NC there were significant differences in megabenthic composition across habitats. Hard habitats in and around cold seeps had significantly higher values of species richness, density, and biodiversity than soft habitats. Depth and habitat complexity were the leading environmental variables driving megabenthic variability. The presence of microbial mats and gas bubbling sites had a statistically significant contribution to explaining megabenthic variability mainly in the shallower BC and less in the deeper NC areas examined; drivers behind this discrepancy could be related to differences between BC and NC in terms of chemical compound fluxes and megafaunal life history characteristics. Our surveys revealed marine litter, primarily from commercial fisheries. This study highlights the importance of habitat complexity for the proliferation of highly diverse cold-seep ecosystems and underscores the importance of discovery science to inform spatial management of human activities in the deep and open ocean.


INTRODUCTION
Cold seep and hydrothermal vent communities comprise unique deep-sea ecosystems as their survival and proliferation does not depend primarily on photosynthetically produced organic matter; on the contrary these ecosystems rely on reduced chemical compounds (e.g., reduced sulfur, methane) which are released from the subsurface to the seafloor (Levin et al., 2016). Cold seeps have a wide geographical distribution and are found across diverse geological settings such as passive continental margins (Paull et al., 1984), subduction zones (Kulm et al., 1986), and oil/gas seismic wipeout zones (Kennicutt et al., 1985).
Research in cold seep ecosystems has documented the roles of environmental parameters in shaping community structure, mainly for macrofauna (Robinson et al., 2004;Levin, 2005;Cordes et al., 2010a,b;Levin et al., 2010;Bourque et al., 2017;Robertson et al., 2020). For example, macroinfauna have shown depth-related patterns with different communities occupying upper-bathyal (200-1,500 m) and lowerbathyal/abyssal (>1,500 m) depths (Bernardino et al., 2012). The important role of habitat in shaping seep macrofauna has also been examined. Near Baltimore Canyon (∼400 m; BC hereafter) and Norfolk Canyon (∼1,500 m; NC hereafter) in the northwest Atlantic, macrofaunal densities in microbial mats were four times greater than those in mussel beds and slope sediments (Bourque et al., 2017).
Studies on cold seep megafauna have focused on dominant species such as tube worms, bathymodiolin mussels and vesicomyid clams with studies ranging from the whole organism down to molecular levels (Kojima, 2002;Coykendall et al., 2019). The ecology of cold seep megafaunal communities is not very well known yet (Baco et al., 2010;Cordes et al., 2010a;Sen et al., 2019;Zhao et al., 2020;Dueñas et al., 2021). This knowledge gap is especially evident in the United States Atlantic Margin where widespread methane seepage has recently been discovered (Skarke et al., 2014). Habitat diversity within United States Atlantic Margin canyons and seeps was demonstrated to exert an important influence on fish assemblage structure, and some fish species exhibited close associations with cold-water coral and sponge habitats (Ross et al., 2015). However, information on the roles of habitat type (e.g., soft/hard/mixed, biogenic habitats), interspecific relationships, resource availability (e.g., presence of microbial mats) and physical oceanography in shaping megabenthic seep communities is almost absent (Turner et al., 2020). Addressing this gap is an urgent need as areas in the United States Atlantic Margin support intensive fisheries and are under consideration for future wind power and oil and gas exploration (CSA Ocean Sciences Inc. et al., 2017).
Advancing knowledge about the environmental parameters shaping the distribution and composition of cold seeps' megafauna will improve model predictions about the geographic and bathymetric distribution of cold seeps (Cordes et al., 2010b;Georgian et al., 2019) and will assist in the establishment of efficient conservation strategies for their protection .
Considering the major knowledge gaps about seep megabenthos, this study focused on seep communities found in two very different bathymetric regimes, i.e., seeps adjacent to BC and NC. Three major research questions were addressed: (1) Which are the major megabenthic taxa found in the areas?
(2) How do the community composition and structure of the megabenthos change across depth and within each seep system? (3) How do environmental parameters such as the type of habitat and the presence of gas bubbling sites and microbial mats affect the community composition?
This study builds on previous research for BC and NC seep sites on the Mid-Atlantic Bight in the United States Mid-Atlantic Margin (e.g., Ross et al., 2015;Bourque et al., 2017;Demopoulos et al., 2019), advancing our understanding of the structure and functioning of deep-sea seep ecosystems in the northwest Atlantic.

Study Areas
Water masses and circulation in the Mid-Atlantic Bight have been studied (Obelcz et al., 2014;CSA Ocean Sciences Inc. et al., 2017 and references therein). The major oceanic source water masses in the Mid-Atlantic Bight are Shelf Water (<200 m water depth), North Atlantic Central Water or West North Atlantic Central Water (<500 m), Labrador Sea Water (>500-1500 m), Western Atlantic Subarctic Intermediate Water (>500-1500 m), Gulf Stream (>500-1500 m), and North Atlantic Deep Water (1500 m-bottom) (see Table 5.1 in CSA Ocean Sciences Inc. et al., 2017 for water mass characteristics).
Two large methane seep environments located near similarly sized shelf-sourced canyons in the Mid-Atlantic Bight, along the United States Atlantic Margin (Obelcz et al., 2014) were explored in 2012 and 2013 (Figure 1). Hecker et al. (1983) first suggested a potential seep on the continental slope near BC, and this was confirmed at water depths between 366 and 450 m (CSA Ocean Sciences Inc. et al., 2017). The second seep located about 20 km south of the mouth of NC, was first identified by sonar by the NOAA vessel Okeanos Explorer at depths between 1,457 and 1,602 m (Skarke et al., 2014). Uranium-thorium dating on authigenic carbonates suggest seepage at Baltimore Canyon between 14.7 ± 0.6 ka and 15.7 ± 1.6 ka and at the Norfolk field between 1.0 ± 0.7 ka and 3.3 ± 1.3 ka (Prouty et al., 2016).
Both sites contained extensive living mussel communities within the Bathymodiolus childressi-complex (Olu-Le Roy et al., 2007;Coykendall et al., 2019). Mussel patch sizes were more variable at NC sites, ranging from a few individuals to several hundred per km-squared . The lineage B. childressi forms a monophyletic group with species in the genus Gigantidas, as reclassified by Thubaut et al. (2013). Since recent studies for this area used B. childressi (CSA Ocean Sciences Inc. et al., 2017;Coykendall et al., 2019;Demopoulos et al., 2019), we will continue with that name for consistency. It should be clarified though that the species name Bathymodiolus childressi is not accepted anymore at the World Register of Marine Species 1 ; the valid species name is Gigantidas childressi 2 . Both BC and NC sites contained microbial mats, authigenic carbonate rocks, and methane bubble activity (Bourque et al., 2017). NC showed high seep activity along two well-defined ridge features separated by about 1 km and contained stable methane hydrates. In contrast, the BC sites were more concentrated in a single area of relatively flat sandy bottom (Skarke et al., 2014;CSA Ocean Sciences Inc. et al., 2017).  Figure 1). The dives began at the deepest location and moved upslope (navigation and video facilitated by multibeam sonar bathymetry obtained in 2011), and the ROV position was recorded continuously with an ultra-short baseline tracking system, with position data synchronized to video data (Ross et al., 2015). A SeaBird SBE 911+ conductivity-temperature-depth (CTD) logger attached to the ROV recorded density (σθ, kg m −3 ), depth (m), dissolved oxygen (DO, ml/l), pH, salinity, turbidity (formazin turbidity units), and temperature ( • C), once per second throughout each 1 http://www.marinespecies.org/aphia.php?p=taxdetails&id=420692 2 http://www.marinespecies.org/aphia.php?p=taxdetails&id=1346725 dive ( Table 1). To standardize operations as much as possible, ROVs ran near the bottom along video transects at speeds of <0.93 km/h, with color video cameras tilted slightly downward and set to wide angle. Two laser pointers spaced 10 cm apart were positioned directly in line with the video camera and were switched on during most dives (Ross et al., 2015). The video recorded continuously and across varied transect lengths within each ROV dive, covering all habitat types. Transects (i.e., anytime the ROV was moving) were interspersed with the ROV stopping to collect regular grab/suction samples of megafauna (to confirm species identifications). Erroneous bottom position track data (e.g., cases where the actual distance traveled by the ROV was greater than the theoretical maximum distance) were removed, following Partyka et al. (2007); Quattrini et al. (2012), and Ross et al. (2015).

Classification of Habitats
Seven cold seep habitat types were defined following Ross et al. (2015) and CSA Ocean Sciences Inc. et al. (2017). The habitats were: (i) sand-mud (S), (ii) sand mixed with dead mussels (SDM), (iii) sand with dead and live mussels (SDM + LM), (iv) mixed hard-soft (Mix), (v) mixed hard-soft with dead mussels (MixDM), (vi) mixed hard-soft with live mussels (MixLM), (vii) mixed hard-soft with dead and live mussels (MixDM + LM) (Figure 2). The habitats sand-mud and sand mixed with dead mussels were regarded as soft habitats. The latter was regarded as a soft habitat as it was mainly composed from soft sediments and dead mussels were a minority. The other five habitats, i.e., sand with dead and live mussels, mixed hard-soft, mixed hard-soft with dead mussels, mixed hard-soft with live mussels, and mixed hard-soft with dead and live mussels were regarded as hard habitats.

Annotations of Megafauna
In total, ∼ 57 h of videos were recorded, which were viewed in QuickTime Player (version 10.5) on Apple ProRes 422 HD in resolution 1920 × 1080 at 29.7 frames per second (FPS). The videos were paused in "frame grabs" to enumerate and identify TABLE 1 | Environmental parameters for water depth (m), water temperature ( • C), salinity (psu), dissolved oxygen (DO; ml/l), density (σθ, kg m −3 ), and pH for each of the ROV dives at the Baltimore Canyon (BC) and the Norfolk Canyon (NC).
benthic invertebrate megafauna to the lowest possible taxonomic level. Each frame was given an Eastern Daylight Time (EDT) stamp following Ross et al. (2015). The initial number of extracted frames was 2858. From these we excluded those frames with either low image quality or frames having no fauna (n = 514). We also excluded those frames where the ROV was stationary (n = 269). This resulted in 2075 high-quality frames that we analyzed in the present study. Taxonomic identifications were based on the best available literature (CSA Ocean Sciences Inc. et al., 2017;Turner et al., 2020) and guidance from world-class taxonomic specialists.
The number of megabenthic individuals was counted in each image. Counting was conservative to avoid misidentifications and was carried out for almost all megafauna seen in the images. Presence/absence data were recorded for colonial species (e.g., encrusting sponges) and for four solitary species (the polychaete Hyalinoecia sp., the planktonic tunicates Salpidae sp., the red shrimps, and the comb jellies) because their massive numbers made counting unreliable. When it was too challenging to identify individuals to species-level, they were combined and reported as a higher taxonomic rank, following Quattrini et al. (2012). Those individuals which could not be identified as a group due to the lack of taxonomic identification, were grouped as morphotypes, e.g., sponges were grouped into two morphotype categories (1, encrusting, cushion; 2, massive, papillate), following Kazanidis et al. (2019).
Depth measurements from the SeaBird SBE 911+ data logger, latitude/longitude (universal transverse mercator; UTM), temperature, density, salinity, pH, oxygen concentration, and turbidity were time synchronized to each annotated frame. This synchronization was based on the time stamps (Hours: Minutes: Seconds) of each image frame and collected environmental data. The presence of gas bubbling sites, microbial mats, and frozen hydrates were also noted (Figure 1 and Supplementary Tables 1-3) and the distances between gas bubbling sites, microbial mats, frozen hydrates and the points that high quality images were extracted, were calculated in QGIS v.3.16.2 using the Analysis Tool "Distance Matrix" (QGIS Development Team, 2019). Fish trap lines, potential trawl marks and litter were also noted, including their depth, latitude and longitude values.

Data Analysis
Each high-quality image extracted from the ROV videos was set as the sampling unit. Samples where no faunal records were made were excluded from further analysis. In addition, three taxa (the bubblegum coral Paragorgia arborea, the soft coral Parazoanthus sp. and unidentified hydrozoans) were excluded from data analysis as they were absent in >90% of the samples Ross et al., 2015).
In each high-quality image we measured the number of species/morphotypes (SR), total number of individuals (N) and the biodiversity indices Margalef (d), Pielou's evenness (J ), and Shannon (H ). The calculation of d, J , and H was carried out in the software Primer v.7 (Clarke and Gorley, 2015). Statistically significant differences across habitats in terms of SR, N, d, J , and H were examined using analysis of variance in R, following Kazanidis et al. (2018). In the case of one-way ANOVA, multiple comparisons were carried out through the Tukey's test. In the case of one-way analysis of means the comparisons were carried out through the Games Howell test (Burk, 2018). When the Kruskal-Wallis test was used, the pairwise comparisons were carried out through the Dunn test (Dinno, 2017). Accounting for the multiple comparisons, the p values were adjusted using the Bonferroni correction (Armstrong, 2014).
Analyses on the faunal composition and structure of megabenthic associations were also completed using Primer. Data on the abundance of megabenthos were square-root transformed and were used in the calculation of Bray-Curtis similarities and similarity matrices. Based on these matrices, nonmetric multi-dimensional scaling (nMDS) 2-dimensional (2D) plots were constructed. Analysis of similarities (ANOSIM) was used to identify the existence of statistically significant differences in community composition between BC and NC as well as across habitats for each seep system. In addition, SIMPER analysis was used to identify species responsible for the average dissimilarity between (a) BC and NC, (b) pairs of habitats within BC and NC.
Finally, the roles of environmental parameters (i.e., depth, latitude, longitude, type of habitat, temperature, salinity, density, oxygen concentration, pH, turbidity, presence of microbial mats, presence of frozen hydrates, and presence of gas bubbling sites) in shaping megabenthic associations within each of the six ROV dives, were analyzed. The role of microbial mats, frozen hydrates and gas bubbling sites in shaping megabenthos was examined by measuring their numbers within a 20 m radius for each single image frame. This radius was set noting that bacterial mats and gas bubbling sites affect megafauna communities on a scale of tens of meters (Levin et al., 2016;Sen et al., 2019). The role of environmental parameters in shaping megabenthos was examined considering megafauna and environmental data recorded in each image frame (see section "Annotations of Megafauna" for details). The role of environmental variability in shaping megabenthic communities was examined in Primer using distance-based linear modeling ("Dist-LM" routine) (Anderson et al., 2008). This routine quantifies the relative contribution of environmental variables in shaping the variability of biological communities (Anderson et al., 2008). To avoid Type I error inflation, explanatory variables were checked for correlation (using Draftsman plots in Primer) and variables with significant correlation scores (Pearson's r > 0.7 or <−0.7) were discarded (Wei and Simko, 2017). A stepwise selection was applied to retain the statistically significant explanatory variables and the Akaike information criterion (AIC) was used as a selection criterion (Anderson et al., 2008).

Distribution of Habitat Types
The distribution of the seven types of habitats differed between the six ROV dives (Figure 3). In NF-07, NF-08, NF-14, and J2-689 dives there was a strong presence of sand-mud habitat while sand with dead and live mussels was dominant in dives J2-682 and J2-683.

Community Composition and Structure
In total 40 discrete taxa or morphotypes were recorded across the two seeps, and 37 of these were included in the analyses. The three that were excluded -as they were seen only in 13 images in total-were Paragorgia arborea, Parazoanthus sp., and hydrozoans. The 37 taxa/morphotypes were distributed as follows: 24 (∼65%) were exclusively found at BC, nine (∼24%) were exclusively found in NC and four (∼11%) were found in both areas ( Table 2). The distribution of the total number of taxa was as follows: 11 arthropods (∼30%), ten echinoderms (∼27%), seven anthozoans (∼19%), five mollusks (∼14%), two sponges (∼5%), one polychaete (∼3%), and one tunicate (∼3%) (Figure 4 and Table 2). 2 | Megafaunal categories included in the still image analysis, with cross-references to the example images in Figure 4.
At BC, habitat complexity had a prominent role in shaping megafaunal composition ( Figure 5B). There was a clear clustering of communities found in soft habitats (sand-mud, sand mixed with dead mussels) versus those found in hard habitats (e.g., mixed hard-soft, mixed hard-soft with dead mussels). ANOSIM analysis revealed statistically significant differences for 14 out of the 21 groups examined (Supplementary Table 4). Table 2. See Table 2 for description of each image. Black arrows in panel 8 highlight the specimens of Alvinocaris markensis, in panel 25 the ophiuroid specimens, and in panel 35 highlights encrusting sponge.

FIGURE 4 | Photographs illustrating the megafaunal categories listed in
Overall, the differences between the BC groups were mainly driven by the sponge Polymastia sp., the anthozoans Actinoscyphia sp., Hormathia sp. (white morphotype), and Bolocera tuediae, the arthropod Chaceon quinquedens, "red shrimp, " "white spider crab, " Tomopaguropsis sp., and the polychaete Hyalinoecia sp. (Supplementary Table 6). For example, the dissimilarity in the pair sand-mud vs. mixed hardsoft with dead and live mussels (the pair with the highest level of dissimilarity) was driven by the anthozoans B. tuediae (16.15% contribution in dissimilarity) and Hormathia sp. white morphotype (15.33%), and the sponge Polymastia sp. (11.47%). Hormathia and Polymastia had higher average densities in mixed hard-soft with dead and live mussels than sand-mud, while B. tuediae had higher average density in sand-mud than mixed hard-soft with dead and live mussels (Supplementary Table 6).
The type of habitat also had a major contribution in shaping megabenthic communities in NC ( Figure 5C). ANOSIM analyses showed significant differences for sand-mud vs. sand with dead mussels (R = 0.197, p = 0.019), sand-mud vs. mixed hard-soft with dead and live mussels (R = 0.255, p = 0.001), sand with dead mussels vs. sand with dead and live mussels (R = 0.076, p = 0.001), sand with dead mussels vs. mixed hard-soft with dead and live mussels (R = 0.632, p = 0.001) and sand with dead and live mussels vs. mixed hard-soft with dead and live mussels (R = 0.057, p = 0.049) (Supplementary Table 5). The dissimilarity in the pair sand with dead mussels vs. mixed hard soft with dead and live mussels (the pair with the highest level of dissimilarity) was driven by Ophiuroidea (43.20% contribution in the average dissimilarity) and Echinoidea 1 white morphotype (28.52%). Ophiuroids had higher average density at sand mixed with dead mussels while Echinoidea 1 white morphotype had higher average density at habitat mixed hard-soft with dead and live mussels (Supplementary Table 7).
Values of SR, N, d, J , and H showed statistically significant differences across the types of habitats, especially at BC dives. Overall, higher values of megafaunal density and biodiversity were found in hard (e.g., mixed hard-soft with dead mussels, mixed hard-soft with dead, and live mussels) than soft (sandmud, sand mixed with dead mussels) habitats (Figure 6 and Table 3).

Role of Environmental Parameters in Shaping Seep Community Structure
Dist-LM analyses showed that the amount of biological variability explained by certain key environmental factors ranged from 17.31% (J2-682 dive) to 31.41% (NF-07 dive). In five out of the six ROV dives, habitat was the leading parameter explaining biological variability (Table 4). At NC (J2-682 and J2-683 dives), the type of habitat, microbial mats, pH, and gas bubbling sites were the parameters that had a statistically significant contribution in explaining variability, while at BC (i.e., NF-07, NF-08, NF-14, and J2-689 dives) there was a larger number of variables shaping megabenthic communities ( Table 4).
The presence of microbial mats and gas bubbling sites had a minor but statistically significant contribution toward explaining community structure at all four BC dives. In contrast, the NC microbial mats and gas bubbling sites had a statistically significant contribution only in J2-683 dive ( Table 4). SIMPER analysis revealed the species with the higher contribution to megafaunal dissimilarities between locations with and without microbial mats/gas bubbling sites (Supplementary Tables 8, 9). At BC the species with higher average density in locations with microbial mats were mainly the hermit crab Tomopaguropsis sp. and the anthozoan Hormathia sp. white morphotype (Supplementary Table 8). SIMPER showed the species with higher average density at BC locations with gas bubbling sites were mainly Hormathia sp. white morphotype, the red crab Chaceon quinquedens and the polychaete Hyalinoecia sp. (Supplementary Table 9). At NC, the species with higher average density in locations with microbial mats and gas bubbling sites were the "Echinoid red morph" and "Echinoid white morph" (Supplementary Tables 8, 9).

Records of Marine Litter in the Baltimore and Norfolk Canyons
The presence of marine litter in the two seeps was limited to fifteen images containing litter items: two images in NF-07, two Habitat types: sand-mud (S), sand mixed with dead mussels (SDM), sand with dead and live mussels (SDM + LM), mixed hard-soft (Mix), mixed hard-soft with dead mussels (MixDM), mixed hard-soft with live mussels (MixLM), mixed hard-soft with dead and live mussels (MixDM + LM).

Changes in Megabenthic Associations Across Depth
Differences in depth played a major role in structuring the megafaunal communities between the shallower Baltimore Canyon (340-496 m) and deeper Norfolk Canyon (1,444-1,611 m) seep sites. There were only four common taxa/morphotypes (e.g., Lithodes sp., Actinoscyphia sp., "red shrimps, " Cerianthus sp.) shared between BC and NC seeps and overall, densities of such species were higher at BC than NC. Two of the relatively common faunal components that were found solely at BC were comb jellies and salps. Cacchione et al. (1978) observed salp carcasses close to the seabed in the Hudson Canyon at ∼3,400 m water depth; salp carcasses might supply more than half of the daily energy requirements of benthic microfauna in that area (Wiebe et al., 1979). On top of northwest Atlantic canyons, jellyfish have been shown to play an important role in pelagic carbon cycling (Condon et al., 2011;Fuentes et al., 2018) and can also form a substantial food input for benthic communities (Luo et al., 2020) including those in the deep sea (Billett et al., 2006;Lebrato et al., 2012;Sweetman et al., 2014). Our findings about sharp differences in megafaunal composition between BC and NC agree with previous work from Brooke et al. (2017) showing differences in the densities of hexacorals and octocorals between BC and NC. Our findings are also in line with those from Turner et al. (2020) reporting almost no overlap between BC and NC megabenthos. Studies on BC and NC demersal fishes showed depth was the key factor separating communities in two zones with a boundary around 1,400 m (Ross et al., 2015). Studies at BC and NC macroinfauna (∼180-1200 m water depth) showed a departure from the expected depth-density relationship, driven by higher faunal abundance at 800 and 900 m which was a deposition zone for organic matter (Robertson et al., 2020). Studies by Robertson et al. (2020) have also shown that NC and adjacent slope were more organically enriched than BC and slope. Furthermore, in Robertson et al. (2020) it was shown that upper slopes in both areas contained macrofaunal communities that were more disturbed than those communities in deeper areas; it was speculated that this disturbance may be due to impacts from fishing activities or interactions of the shelf with hydrography (CSA Ocean Sciences Inc. et al., 2017).
Our results agree with findings in the deep northwestern Atlantic and with findings about zonation of deep-sea biota on continental margins (Rowe et al., 1982;Carney, 2005 and references therein). Specifically, evidence has been provided for the occupation of different bathymetric zones between mussels: in the Gulf of Mexico Bathymodiolus childressi was found in shallower waters than Bathymodiolus heckerae (1,005-2,284 and 2,180-2,745 m, respectively), with the two species never found together (Cordes et al., 2010a). In the Gulf of Mexico's seeps little or no overlap was found between upper and lower slope macroor megafaunal communities. A striking difference between the two zones was the dominance of the ophiuroid Ophioctenella acies in the deeper communities; the shift in faunal composition seems to take place between 1,300 and 1,700 m (Cordes et al., 2007). Thus, there is further reinforcement for a widespread, faunal bathymetric transition zone in the region near 1,400-1,500 m as suggested by Ross et al. (2015).
The examination of biodiversity gradients across depth in the North Atlantic have revealed different patterns. Firstly, studies in northwestern Atlantic have shown the presence of a unimodal distribution for macrofauna (gastropods, polychaetes, protobranchs, and cumaceans) as well as for invertebrate and fish megafauna (e.g., Rex, 1981;Pineda and Caswell, 1998;Stuart and Rex, 2009). In this pattern, relatively low values of species diversity are seen in the continental shelf, followed by an increase and a peak at mid to lower bathyal zone and then decreases toward the abyssal plains (Rex, 1981). Comparisons among recorded unimodal distributions have revealed variability regarding the depth of the peak (Rex, 1981;Etter and Grassle, 1992;Allen and Sanders, 1996;Levin et al., 2001). The unimodal distribution of diversity is possibly explained by the dynamic equilibrium model predictions (Huston, 1979) where species interactions have an important role in shaping patterns of diversity (Rex, 1981;Stuart et al., 2003). Furthermore, studies on benthic megafauna on the continental margin south of New England (242-2,394 m depth) defined four megafaunal slope zones: upper, upper-middle, transition on the lower-middle, and a lower slope zone (Hecker, 1990). Faunal density was high in the upper and lower slope and low in the two middle slope zones; density and species composition patterns across depth were similar at three of the locations but they were significantly different at the eastern edge of Georges Bank, an area characterized by steep topography, glacial erratics, and outcrops (Hecker, 1990). Furthermore, studies on megafaunal assemblages off Cape Hatteras (157-1,924 m) showed that populations in the upper and lower slope were similar in terms of density and species composition, to those found at other locations; in contrast, megafaunal abundances were elevated in mid slope, mainly reflecting dense populations of fish, eels and anemones, probably related to enhanced food availability (Hecker, 1994). In the Rockall Trough (northeast Atlantic), Paterson and Lambshead (1995) have shown a parabolic distribution in polychaetes' species richness with a peak located at about 1,800 m water depth. The authors mentioned that disturbance, in terms of the frequency of high energy currents, and intrataxocene predation was important on the Hebridean Slope. Gage et al. (2000) have shown reduced macrobenthic species diversity at ∼400 m depth in the continental slope off Scotland, probably due to hydrodynamic disturbance, but no clear mid-slope peak in species diversity was found. The examination of macrobenthic communities in the Goban Spur region (northeast Atlantic continental slope) showed an increase in species diversity from ∼200 down to 4,500 m water depth; species richness, however, had a parabolic pattern with a peak in the upper slope (1425 m). Differences in lifehistory characteristics and food supply could explain to some extent the diversity patterns found across the depth gradient (Flach and De Bruin, 1999). Finally, in the Porcupine Seabight and adjacent Abyssal Plain, Olabarria (2005) recorded a pattern of increasing bivalve diversity from ∼500 to 1,600 m followed by a decrease to minimum at about 2,700 m, followed by an increase up to maximum at ∼4,100 and then a drop to ∼4,900 m. The low values of diversity between 2,100 and 2,700 m were associated to an extent to the high abundance of the bivalve Kelliella atlantica while high values of diversity in the abyssal depths were associated in part with the source-sink hypothesis (Rex et al., 2005), low predation (Allen and Sanders, 1996), and strong fluxes of organic matter (Billett et al., 2001).

Changes in Megabenthic Associations Across Habitats
Habitat type played an important role in shaping the composition of the seep megabenthic communities. At both BC and NC, the communities found in soft habitats were clearly separated from those communities on hard habitats. At BC the sponge Polymastia sp. was among the species with the highest contribution in the dissimilarities found between soft and hard habitats. The close association of Polymastia sp. with hard habitats in coastal southern Brazil (Bumbeer et al., 2016) and northern Norway (Dunlop et al., 2020) facilitates the successful settlement of larvae and survival of juveniles. Soft sediments may not be ideal for the proliferation of Polymastia sp. due to the negative impacts caused by sedimentation on filtering and oxygenation in sponges (Kutti et al., 2015). Hormathia sp. anthozoans had higher densities in hard habitats than soft habitats which agrees with findings on Hormathia pacifica recorded living upon manganese encrusted basalt and dead sponge tissue in the Taney Seamounts (Sanamyan et al., 2015). Our results fit well with previous studies on Belgica Province carbonate mounds where on-mound habitat had higher biodiversity than off-mound (Henry and Roberts, 2007).
Contrary to Polymastia sp. and Hormathia sp., the red shrimps had stronger presence in soft than hard habitats. Due to low taxonomic resolution it is challenging to define the drivers that shape the recorded patterns. Some Penaeidae living in deep waters (e.g., Parapenaeus longirostris found in the eastern Atlantic from 20 to 750 m water depth) inhabit sand-mud bottoms and feed on a large variety of bathypelagic, benthic and endobenthic prey (Kapiris, 2004). The Hyalinoecia sp. polychaetes also had stronger presence in habitats composed of soft sediments and dead mussels than those containing a mix of hard and soft sediments and live and dead mussels. Our data agree with studies on the ecology of Hyalinoecia artifex from the United States Continental Margin, noting this species is a predator/scavenger found mainly on soft sediments and is negatively associated with authigenic carbonates (Meyer et al., 2016). Bolocera tuediae anemones in our study were found both in hard and soft sediments; in previous studies they have been recorded in high abundance in bioturbated soft sediments in the Franken Carbonate Mound (Wienberg et al., 2008). In NC, SIMPER analysis provided evidence that "Echinoid white morphotype" had higher densities in hard than soft habitats. The drivers behind this pattern are not known, but we propose this is related to feeding. Studies on the carbon stable isotope values of Echinus sp. in deep-sea cold seeps in the Mediterranean provided some preliminary indication of reliance on chemosynthesis originated carbon (Olu-Le Roy et al., 2004). Further work is needed to  KW, Kruskal-Wallis rank sum test; F, oneway test not assuming equal variances. p values are given (*0.01 < p < 0.05, **p < 0.01, ***p < 0.001, "ns" non-significant). Types of habitats: sand-mud (S), sand mixed with dead mussels (SDM), sand with dead and live mussels (SDM + LM), mixed hard-soft (Mix), mixed hard-soft with dead mussels (MixDM), mixed hard-soft with live mussels (MixLM), mixed hard-soft with dead and live mussels (MixDM + LM).
determine whether echinoids in NC contain symbiotic bacteria assimilating chemical compounds, graze on bacteria or feed on mussel tissues. Apart from differences between soft and hard habitats, differences were also recorded between other habitats. The red crab Chaceon quinquedens had higher density in habitats with live Bathymodiolus mussels than habitats without live mussels, suggesting an affiliation between the two. During this project, C. quinquedens was observed in dense live mussel beds at BC, apparently feeding upon mussel tissue (S. W. Ross, pers. obs.). Stable isotope analysis from the United States Continental Margin showed the proportion of methane-derived carbon within crab tissue is highly variable (∼0% to ∼30-50% of nutrition), depending on location (Turner et al., 2020); they also suggested that some seep mussel beds may act as a nursery ground for C. quinquedens.
A relatively small number of species (mainly echinoids and ophiuroids) were responsible for the dissimilarities found between habitats at NC, in contrast to BC. Although direct evidence is lacking, we hypothesize this observation is related to the life-history characteristics of species recorded at NC. Most of them were mobile (e.g., decapods, sea spiders, echinoids, ophiuroids, and holothurians) and apparently do not have equally strong affiliations with a specific type of substrate as seen at BC (e.g., Polymastia and Hormathia). The holothurian Chiridota heheva seen at NC is a cosmopolitan chemosynthetic ecosystem species (Thomas et al., 2020) feeding on various food sources (Carney, 2010). Pantopoda and Lithodes crabs in chemosynthetic environments (e.g., Lithodes longispina) have been characterized as predators/scavengers (Bowden et al., 2013;Wang et al., 2016;Dietz et al., 2018) moving across habitats and searching for prey.
The existence of clear differences among habitats both at BC and NC agrees with previous studies. At BC there were higher mean densities of macroinfauna in microbial mats (∼80,000 ind/m 2 ) than background sediments (∼15,000 ind/m 2 ) (Bourque et al., 2017). At NC, Bourque et al. (2021) recorded similar macroinfaunal densities among habitats sampled at comparable depths, but diversity was higher near the hard substrate environments, likely due to their increased habitat heterogeneity and enhanced food supply. Differences in demersal fish communities among BC habitats were also reported with assemblages in sandy habitats being significantly different from the most complex habitats, e.g., mixed hard substrates with soft sediments (Ross et al., 2015). In the Barents Sea taxonomic richness and abundance of megabenthos were much higher at the Svanefjell seep site compared with a non-seep site; crusts of seep carbonates drove the higher diversity of the seep site (Sen et al., 2019). Ross et al. (2015) reported a decrease in the role of habitat complexity in shaping demersal fish communities from the shallower BC to the deeper NC. They suggested the decreased availability of complex habitat (seamounts excluded) and limited food supply as depth increases argue against strong habitat association in deep environments, where flexibility and opportunism might be more advantageous. Our findings do not support a decrease in the importance of habitat complexity for megabenthos with depth; similarly, in Bourque et al. (2017) the type of habitat had a statistically significant contribution in explaining the structure of macroinfaunal communities both at BC and NC. The apparent explanation about this  discrepancy between benthic and benthopelagic ecosystem components may be that macroinfaunal and megabenthic communities have a higher dependency on the seafloor compared to pelagic organisms, and this is reflected in the type of environmental variables explaining faunal distribution across space (Gallucci et al., 2020).

Role of Microbial Mats and Gas Bubbling Sites in Shaping Megabenthic Associations
The presence of microbial mats and gas bubbling sites had a minor but statistically significant contribution in the explanation of variability in megabenthos. Microbial mats are formed where methane flux and anaerobic oxidation of methane rates are high, and therefore the upper layer (∼1 cm) of sediment is rich with sulfide (Foucher et al., 2009). At BC, the hermit crab Tomopaguropsis sp. was among the species with a higher contribution in explaining dissimilarities between locations with and without microbial mats (Supplementary Table 8). Without access to biological samples it is impossible to draw conclusions on the type (if any) of association between these hermit crabs and microbial mats. However, investigation of carbon sources on Paralomis sp. crabs in cold seeps in Costa Rica revealed microbial mat-derived carbon provided an important contribution to the crab's nutrition. Lipid analyses showed Paralomis were feeding on other 13 C-depleted organic matter, probably symbiontbearing megafauna as well as sedimentary detritus (Niemann et al., 2013). Hermit crabs have also been reported feeding on Beggiatoa mats at the Gullfaks seeps in the North Sea (Hovland, 2007). Apart from the hermit crab Tomopaguropsis sp., the anthozoans Hormathia sp., and Actinoscyphia sp. had higher average density in locations closer to mats and gas bubbling sites, respectively. Hormathia in Concepción Methane Seep Area, Chile was characterized as a non-endemic seep species (Sellanes et al., 2008). Although direct evidence for the BC anthozoans is lacking, it is likely these species can assimilate, through their microbial symbionts, chemical compounds released to the water column (Rodríguez and Daly, 2010;Goffredi et al., 2021). Macroinfaunal densities at BC and NC microbial mats were four times greater than those at mussel beds and slope sediments, dominated by polychaete families; multivariate analysis showed specific sediment properties shaped macrofaunal communities, including larger grain sizes at NC microbial mats (Bourque et al., 2017). In contrast to BC, the presence of microbial mats and gas bubbling sites at NC had a statistically significant contribution in shaping megabenthic associations only in one dive, i.e., J2-683. SIMPER analysis showed that echinoids (white and red morphs) had higher densities in samples with, than those without, microbial mats/gas bubbling sites. Again, without access to biological samples it is impossible to draw firm conclusions on the associations, if any, between these echinoids and the microbial mats/gas bubbling sites. It should be mentioned though that studies analyzing gonad fatty acids have suggested that the urchin Spatangus purpureus feeds on phytodetritus, faunal detritus, microbial mats or a combination of all three sources (Barberá et al., 2011). The possible discrepancy regarding the role of microbial mats and gas bubbling sites between BC and NC is an intriguing finding considering that the numbers of microbial mats recorded at NC were comparable to those at BC (Figure 1 and Supplementary Table 2). The reasons behind this discrepancy are not known. Possibly they are related to spatial variability in the distribution, intensity and/or composition of chemical fluxes and/or the mat microbial species composition. The mapping of 570 seeps in the United States Atlantic Margin showed approximately 440 of them lie on the upper continental slope between ∼180 m (shelf break) and 600 m water depth (Skarke et al., 2014). Furthermore, studies have shown the role of sulfide fluxes in shaping macrofaunal community structure at gas hydrates deposits with flux being highest in Beggiatoa, slightly less in Calyptogena clams, and low in areas with Acharax bivalves (Sahling et al., 2002). Horizontal and vertical distribution of sulfide also shaped the fine-scale distribution, structure, and composition of macrofaunal associations in Pacific methane seeps off California (Levin et al., 2003), while studies in the Gulf of Mexico revealed significant differences in the structure of infaunal communities across the Arcobacter, Thioploca, and Beggiatoa microbial mats (Robinson et al., 2004). The concentric zonation of mussels, clams and Siboglinidae polychaetes around pockmarks in the Blake Ridge Diapir seep field and Congo-Angola Margin suggested the influence of chemical gradients on megafaunal distribution (Olu-Le Roy et al., 2007). The relatively limited role of microbial mats in shaping megabenthic associations at NC could be due to the most abundant species recorded at NC (e.g., ophiuroids, the sea cucumber Chiridota heheva, two morphotypes of sea urchins) not relying heavily on microbial mats as a food source. Studies from Gulf of Mexico cold seeps showed asteroids and ophiuroids had tissue isotope values reflecting the input of free-living microbial detritus in their diet or predation (Carney, 2010); C. heheva can feed on various sources (Carney, 2010) and partially preys on bivalve tissues ; sea urchins on the Oregon and California Margin seem to have photosynthesis-based diets (Levin and Michener, 2002). Stable isotope analysis of BC and NC megafauna will provide an insight to their diets and dependency on microbial mats .

Marine Litter at the BC and NC Sites
Image analysis revealed a relatively small number of marine litter items at BC and NC sites. Some of them were clearly related to commercial fishing activities (e.g., abandoned nets) such as for the red crab C. quinquedens and fishes (Ross et al., 2015;Turner et al., 2020). We note here that both seep sites were outside the main axes of the two major canyons. Since canyons, including BC and NC, are known to accumulate, concentrate and transport debris and litter (Van den Beld et al., 2017;Dominguez-Carrió et al., 2020;S. W. Ross, unpubl. data), we might expect less density of anthropogenic debris on the open slope outside the main canyon axes than within. Despite the small number of items observed, their negative impacts can be long-lasting, particularly in the case where long-lived ecosystem engineers are affected (e.g., see the damaged bubblegum coral, Figure 7G; Williams et al., 2010). In addition to the direct impacts on deep-sea fauna, bottom trawling can further affect the health status of canyon fauna through the alteration of physical properties of surface sediments (Martin et al., 2014), sediment resuspension (Paradis et al., 2017), downslope deposition and smothering (Puig et al., 2012); close monitoring of health status and the establishment of area-based management tools (e.g., marine protected areas) will serve the conservation of rich seep biodiversity and associated ecosystem services for the future.

CONCLUSION
The present study advanced our understanding about the key role that environmental parameters (e.g., depth, type of habitat, and microbial mats) play in shaping the distribution of cold seep megabenthos. This information is crucial for improving model predictions on the spatial distribution of seep fauna and for the successful protection of deep-sea habitats, e.g., through the designation of networks of ecologically coherent marine protected areas. In 2016 the Northeast Canyons and Seamounts Marine National Monument was designated aiming to protect an area of high diversity and ecological connectivity across depths and along the United States Atlantic continental margin (Auster et al., 2020 and references therein; see also links below for the Frank R. Lautenberg Deep-Sea Coral Protection Area which took effect in early 2017) 3 . In addition to its high diversity, the Monument is relatively undisturbed and could serve as a reference point, which is vital for the assessment of deep-sea environmental status . Deep-sea ecosystems are currently threatened by multiple anthropogenic stressors. Development of basin-scale strategic research and establishment of interdisciplinary collaborations are keys for the efficient protection and conservation of fragile deep-sea ecosystems for the future.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SWR designed and co-led the ship and ROV surveys. SWR, JMR, and JC coordinated the data collection. JC analyzed the video material and completed the megafauna annotations. GK and JC carried out the statistical analysis and wrote the manuscript with editorial assistance from SWR and JMR. All authors contributed to the article and approved the submitted version.

FUNDING
This study received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement nos. 678760 (ATLAS) and 818123 (iAtlantic) to JMR. The field data collection was funded largely by the United States Bureau of Ocean Energy Management (BOEM) (contract with CSA Ocean Sciences Inc., who was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.). Ship and ROV time were provided by the NOAA Office of Ocean Exploration. The USGS provided additional field support.