Samples and measurements were collected during three cruises with RV Pelagia (May 2014, April 2015, and June/July 2016) around the Rainbow hydrothermal vent field (Figure 1). This vent field, located at the Mid Atlantic Ridge (MAR) at 36°14′N, 33°54′W and at 2,300 m depth was chosen as study site due to its extensive buoyant plume (German and Parson, 1996; Severmann et al., 2004). The vent field is known to produce plumes relatively rich in transition metals resulting in SMS deposits enriched in Cu, Zn, Co, and Ni (Charlou et al., 2002; Douville et al., 2002; Findlay et al., 2015). This plume could still be detected up to 50 km away from the vent site, moving predominantly to the north and east, following the contours of the Rainbow Ridge at ∼2,100 m depth (Severmann et al., 2004). Moreover, the plume could be consistently traced up to 25 km away from the vent source based on turbidity measurements and the chemical analysis of water samples (Haalboom et al., 2019). The study of Haalboom et al. (2019), showed an enrichment of metals compared to the overlying clear water in the neutrally buoyant plume with a high proportion of chalcophile elements closest to the vent site.

Sediment samples were collected using the NIOZ box corer with a 25 cm radius and 55 cm sampling depth (equivalent to a maximum of roughly 100 L of sediment). The box corer is equipped with a tightly closing lid to prevent the sediment from being washed out during retrieval. A total of 15 stations were sampled, of which 13 sites within a 5 km radius of the Rainbow vent site and 2 reference sites at approximately 25 km distance from the vent site. From each box core two samples were collected for molecular identifications. A smaller sediment sample was collected using a sub corer with a 1.25 cm radius and 10 cm sampling depth (equivalent to ±50 mL of sediment) and immediately stored at −80°C. And a larger sample (approximately 500 cm² surface) was collected from the top 15 cm of the remaining sediment, sieved over a 200 μm mesh sieve on board and this was stored at −80°C as well, immediately upon sieving. Next to the fauna samples, a sub core (5 cm radius) was retrieved for sedimentological analysis and stored at 4°C. Collection of sedimentological sample was not possible at sample location 51.

Water samples were collected from the water column 5 m above the seafloor either using 24 Niskin bottles of each 12 L, which were mounted on the CTD-frame (total volume 288l), or using the NIOZ Large Volume Sampler (LVS), with a volume of 1,000 L (Duinveld and Shipboard Scientific Party, 2016). The LVS consists of a cylindrical polyester tank with 0.5 m diameter openings at the top and bottom which can be remotely closed with a lid through the ship’s CTD operation unit. On board the water contained in the tank and Niskin bottles was drained off through a tap at the bottom of the tank and filtered over a 30 μm nylon mesh sieve. The sieve residues were stored at −80°C. In this way, a total of 17 water samples were taken of which 15 within a 5 km radius of the vent site and 2 reference sites (Figure 1).

Major and trace elements were determined for the top 1 cm of all box-cores by Inductively Coupled Plasma Mass Spectrometry (ICP-MS) analyses. A 100 mg of cryodesiccated surface sediments were completely fully digested in 3.5 mL of HF:HNO₃ mixture (2.5:1), to which 1 mL 1M HCl and 1 mL perchloric acid were added. The sample-acid mixtures were heated in closed Teflon bombs for 48 h at 125°C in a hot block, after which the bombs were opened and left to evaporate for 24 h at 125°C. The dried residue was subsequently transferred to in 1 M HNO₃ and an internal standard was added. Before analyses the samples were diluted by a factor ten thousand with milli-Q water. The obtained solution was measured with the Thermo Scientific Element-2 ICP-MS, including two blanks and three standards (PACS-2, JDo-1, and MESS-3). These standards were chosen for their expected element resemblance with the TREASURE cores. ELEMENT 2/XR software was used to convert data into concentrations in ppm. Standard deviations were based on six replications per sample. Detection limit was calculated as three times the standard deviation of six-replicate measurements of a blank, averaged for two blanks.

Environmental DNA (eDNA) is considered to comprise DNA shed from living species and DNA leaking from fragments and dead organisms (Corinaldesi et al., 2011; Taberlet et al., 2012b; Guardiola et al., 2015, 2016). Because of this inclusiveness, for benthic fauna often only small amounts of sediment (up to 10 g) are used for the analysis of the metazoan community (Sinniger et al., 2016; López-Escardó et al., 2018). However, the detection of macrofauna species might be biased as the release and fate of eDNA differs per species (Barnes and Turner, 2016). In this study, DNA was isolated from a small amount of unsieved sediment, to include all smaller specimen and an additional DNA extraction on the sieved sediment fraction was used, to improve the detection of macrofauna species (Brannock and Halanych, 2015; Klunder et al., 2019). Therefore, the sieve residues from the sieved samples were cryodesiccated and ground in an agate mortar using liquid nitrogen. Subsequently, a subsample (10 g) was taken from each sample for further molecular analysis. For the unsieved samples, a subsample (10 g) was taken from each sub core at the following depth intervals: 0–2 cm, 2–3 cm, and 5–6 cm. All subsamples of the sieved (n = 15) and unsieved (n = 3 × 15) fractions served as starting material for further molecular analysis of the benthic fauna. The sieve residue of the water samples served as starting material for further molecular analysis of the pelagic larvae samples. The molecular analysis and subsequent bioinformatics analysis followed the protocols as described in Klunder et al. (2018) using the F04 (5′-GCTTGTCTCAAAGATTAAGCC-3′) and R22mod (5′-CCTGCTGCCTTCCTTRGA-3′) primer pair (Sinniger et al., 2016) except for the DNA extraction of the larvae samples. Here, the DNEASY PowerSoil^TM kit (MoBio Inc.) was used instead of the Powermax Soil^TM kit. Within the bioinformatics analysis, only clusters with at least two reads after de-replication were retained and further clustering of operational taxonomic units (OTUs) was at a 98% similarity cut-off. Taxonomic assignments were performed with a 0.8 minimum confidence against the SILVA 18S rRNA database (release 128, Pruesse et al., 2007) using the RDP Classifier (Wang et al., 2007). The molecular analysis was completed for all samples except for box core location 36 as DNA-extraction failed.

The reproducibility of sequencing output was evaluated using duplicate samples. Each duplicate sample was run in a separate PCR and Illumina sequencing run. The taxonomic composition of metazoan orders was comparable (Supplementary Figure S1). Also, the relative abundance of reads for all metazoan OTUs were highly similar between the two duplicates based on linear regression (R² of 0.88, slope of 0.93; P < 0.001) (Supplementary Figure S2). This similarity means the most abundant OTUs were recovered consistently for each duplicate sample. Therefore, it is assumed potential biases introduced by separate PCR and sequencing runs are small, and reliable estimates about the metazoan diversity can be drawn for all samples within this project.

All data visualizations, calculations and statistical analyses were performed in RStudio. Trace-metal concentrations were calculated in mg/kg (ppm) of the sediment and were normalized to Al, to correct for differences in the aluminosilicate fraction of the sediment. The trace-metal data were prior to ordination with a PCA, normalized using the decostand- function within the R-package ‘vegan’ to distribute the weight of each metal evenly (Oksanen et al., 2018). Bray–Curtis dissimilarities were calculated and used for the PCA. Simper analysis was used to discriminate the effect of each trace-metal.

The read numbers per OTU for all molecular samples were transformed into a relative abundance of reads per OTU per sample. This approach was preferred over rarefaction as it does not omit valid read abundance data (McMurdie and Holmes, 2014; Lanzén et al., 2016). All subsequent analyses were based on these relative data. A threshold for false-positives was calculated from the mock community and was set at a relative read abundance of 0.0005%. The reads were Hellinger transformed (Legendre and Gallagher, 2001) using the vegan function ‘decostand’ and dissimilarity distances were calculated using the Bray–Curtis equation. These dissimilarities were further used for analysis of variance between the locations (permanova) using the adonis2 function and for the simper analysis to discriminate between the effect of each individual OTU.

The concentration of 23 different trace elements in the sediment was measured for all sampling stations (Supplementary Table S1). PCA and cluster analysis based on the element/Al ratios showed a division into three main groups (Figure 2). The reference locations (38 and 43), which were supposed to be the least affected by the vent plume because of their large distance from the vent are clustered together with sample locations 58, 56, 54, and 36 which are also situated further away from the vent site, either upstream or downstream of the plume. The two other sample groups that could be distinguished based on their trace elemental composition are also spatially clustered: either north-east (72, 63, 33, and 88) or south-east (70, 30, 69, and 4) of the vent-site. SIMPER analysis showed that the difference between the reference sites and the NE sites was mainly due to relatively higher element/Al ratios of Zn, Cu, and Co in the latter, elements that are typically present as sulfides near hydrothermal vents. The ratios of these elements were also high at the SE sites compared to the reference sites, however, not as high as at the NE sites (Table 1). The relative concentrations of Al, Ca, and Sr were on average higher at the reference locations. All three elements are associated with background sedimentation of pelagic carbonate settling from overlying water masses and/or aluminosilicates transported from large distances as clay particles. Although the sample locations 58, 56, 54, and 36 clustered together with the reference sites in the ordinations, their levels of trace-metals typical for hydrothermal vents such as Zn, Cu, and Co are still higher compared to the reference locations (Table 1, ‘Other’). Therefore, these samples will be labeled as ‘Other’ for the rest of the analysis and will be treated as plume-affected.

A total of 12.3 million reads passed quality filtering, of which 4.8 million could be confidently assigned to metazoans. Clustering at a threshold at >98% and the removal of singletons resulted in 5,367 metazoan OTUs. MDS ordinations of stations were made based on the relative abundance of reads for all metazoan OTUs separately (Figure 3). We subsequently explored whether distances found in Figure 3 were linked to the chemical composition (Figure 2) by grouping fauna accordingly. The ordination showed a statistically significant difference in composition between the plume-affected sites as compared to the reference sites (Permanova; OTU: F₁,₁₃ = 1.285, p = 0.046). Only few of the separate OTUs could be taxonomically assigned beyond the order level. Therefore, Simper analysis of the contribution by OTUs to the dissimilarity between the plume-affected and reference sites mostly gave information at the order level. Simper analysis showed that the difference between the plume-affected sites and the reference sites was mainly due to a lower abundance of an OTU taxonomically assigned to Enoplida (Nematoda) and a higher abundance of Harpacticoida (Arthropoda) at the plume-affected sites. OTUs which could be identified down to genus level and which were important for the dissimilarity between groups were the annelid Myriochele (Sabellida) and Anobothrus (Terebellida). The relative abundance of Myriochele was highest at the distant reference sites. Also at the Other sites Myriochele abundances made up a comparatively large proportion of the OTUs. The relative abundance of Anobothrus was highest at the plume-affected sites, especially at the sites close to the vent (groups NE and SE).

In total, 18 metazoan phyla were identified (Figure 4). The phyla Annelida, Arthropoda, Nematoda, Platyhelminthes, and Kinorhyncha were recovered from all sampling stations. The arthropods were most proportionally abundant at the two groups of plume-affected sites (i.e., 41% of the metazoan reads for NE and 32% for SE) whereas the nematodes (40%) and annelids (36%) were most abundant at the reference sites. Although there are small differences in taxonomic composition between the stations within a clustered group (i.e., NE, SE, Reference, or Other), the patterns found within each group are distinctly similar (Figure 4). The Shannon-Wiener index was calculated for the relative read abundance of all OTUs for each location. Diversity was on average lowest in the SE locations (H = 2.77) whereas the other locations were comparable (NE: H = 3.29, Other: H = 3.20, Ref: H = 3.29).

Relative abundance of reads at the order level for the three most abundant phyla (Annelida, Arthropoda, and Nematoda) in each cluster are shown in Figure 5. The most pronounced difference within the annelids, is the high abundance of the order Terebellida at the SE sites whereas the order Spionida was most abundant at the other three groups of sites. The arthropod and nematode orders did show a clear distinction between the plume-affected sites (NE and SE) and the reference sites. The arthropod order Harpacticoida was the dominant order at the plume-affected sites whereas the Malacostraca was the most abundant arthropod group at the reference sites. Similarly, the nematode order Monhysterida was most abundant at the plume-affected sites whereas the Enoplida was more dominant at the reference sites.

A total of 2.3 million metazoan reads were obtained for the pelagic larvae samples. The majority of the reads (99.6%) were assigned to pelagic groups such as the planktonic copepod taxon Calanoida; the free swimming tunicate taxon Doliolida and the krill taxon Euphausia. The remaining reads (0.4% ≈ 9,000 reads) were assigned to genuinely benthic taxa. As the number of reads for the benthic taxa were low, only the presence-absence of these taxa were scored (Table 2). Five benthic annelid taxa were present at the NE sites whereas no annelid taxa were found at the SE sites. The only benthic order found at the SE sites was the Actiniaria. Benthic taxa within the Mollusca phylum were only found at the reference sites.

The combination of high-throughput sequencing and metabarcoding has been shown before as a useful methodology for community analysis (i.e., Chariton et al., 2015; Guardiola et al., 2016; Lanzén et al., 2017) and this study further exemplifies its applicability for exploring benthic deep-sea communities which are otherwise hard to analyze. Still, the method is not without caveats. A major shortcoming of metabarcoding are the considerable gaps in the reference sequence databases (Wangensteen et al., 2018). Despite growing barcoding efforts in the past decades, reference sequences for deep sea benthic species are rare (Dell’Anno et al., 2015). Very low numbers of OTUs could be identified down to the family or genus level, especially for OTUs within groups such as Bryozoans, Platyhelminthes, and Nemerteans. Although more specific taxonomic identifications were possible for the Annelida, Arthropoda, and Nematoda, most OTUs within these phyla could still not be identified beyond the family level.

The 18S rRNA gene was used to both infer the benthic community as well as the larval community. The V1-V2 region of the 18S rRNA locus targeted in this study allows identification of species across a wide taxonomic range for benthic communities (Fonseca et al., 2010; Sinniger et al., 2016). The cytochrome c oxidase subunit I gene (COI) has been used in barcoding studies for single larvae within a sanger-sequencing approach (i.e., Brandão et al., 2016) or for the analysis of larvae from secluded groups (i.e., Gollner et al., 2016). However, the COI gene has less-conserved primer binding sites as compared to the 18S rRNA gene, which might lead to false-negatives in a metabarcoding study due to primer mismatches (Deagle et al., 2014). As we were interested in larvae from all benthic species present, we chose the more conserved 18S rRNA gene as our marker. Also, the use same marker for both samples, sediment and larvae, facilitated the direct comparison of OTUs between the samples.

Within this project, we have chosen the SILVA reference database to perform the taxonomic assignment instead of GenBank. The SILVA database is a curated database and few mis-annotations exist in this database as opposed to the GenBank database (Pruesse et al., 2007; Bik et al., 2012). Although the SILVA database is argued as being too strict (Balvočiūtė and Huson, 2017; López-Escardó et al., 2018), we favored a small loss of diversity over the taxonomic annotations of false positives as these false reads might imply the incorrect presence of rare taxa.

Within this study, the relative abundance of reads for all metazoan OTUs was used as proxy for abundance. The use of amplicon data for abundancy measurements has been disputed (Porazinska et al., 2010; Elbrecht and Leese, 2015; Barnes and Turner, 2016; Pawlowski et al., 2018). However, several studies have found positive correlations between the relative abundance of reads and biomass measurements for metazoan groups including marine invertebrates (Lejzerowicz et al., 2015; Aylagas and Rodríguez-Ezpeleta, 2016). The match between actual (relative) biomass and the relative abundance of reads will not be perfect, but we believe that the relative abundance does reflect abundancy differences between major taxonomic groups. Also, the reproducibility of the PCR and Illumina sequencing as shown in this paper support our approach for abundancy measurements.

The analysis of major- and trace metals from the sediments around the Rainbow hydrothermal vent field showed a spatial zonation separating the sample locations, which is related to a spatial pattern in plume fall out. The pattern found is likely caused by a combination of residual currents and distance from the vent site. The plume emitted by the Rainbow hydrothermal vent field has been shown to move predominantly in the NE direction, following the Rainbow Ridge (Severmann et al., 2004). Sample locations in the direct path of the plume (33, 62, 72, and 88; approximately 3 km NE of the vent site) were found to have highest element/Al molar ratios for the chalcophile elements (Zn, Cu, and Co) compared to the other sample locations. These elements precipitate rapidly from vent fluids as sulfides (Edmonds and German, 2004) and higher concentrations of these metals in the proximity of the Rainbow hydrothermal vent have been found before (Cave et al., 2002). The sample locations SE of the vent site (4, 33, 69, 70; approximately 1 km distance) were not in the trajectory of the plume depicted by Severmann et al. (2004). However, the element/Al ratios of the chalcophile elements were still higher than the ratios measured for the reference locations at 25 km from the vent sites. Plume dispersion in the SE direction was already shown (German et al., 2010; Haalboom et al., 2019), which could explain the enrichment by chalcophiles of the sediments at the SE locations. The element/Al ratios for the sediments collected at locations further away from the vent site (36, 54, 56, and 58) were more similar to the ratios found at the reference locations rather than the NE and SE locations. Only sample locations 56 showed slightly elevated levels of both the chalcophiles as well as Ba, Ni, and Mo compared to the reference locations. This location is probably both influenced by the dispersion of the plume in the NE as well as the SE direction.

The spatial patterns found in the distribution of benthic communities match the patterns found for the major- and trace metals analysis. The difference between the benthic communities at the plume-affected sites and the reference sites show an impact of the plume on the benthic communities surrounding the vent area (up to several kilometers). Spatial patterns within benthic communities surrounding hydrothermal vent areas have been shown before (Podowski et al., 2010; Boschen et al., 2016; Sen et al., 2016) and are thought to be mainly defined by chemistry rather than topography (Podowski et al., 2010; Bowden et al., 2016). The fall out of the plume toward the underlying sediments is partly due to the interaction between the chemistry of the vent fluids and the microbial communities in the plume. The vent fluids emitted differ for each vent location and are influenced by the underlying host rock (Wetzel and Shock, 2000; Cave et al., 2002). Subsequently, microbial activity alters the metal composition in the plume by scavenging and oxidation (Cowen et al., 1990; Dick et al., 2009) whereas the microbial community itself also show a shifting composition and diversity with distance from the vent site (Haalboom et al., 2019). The fall out of major elements and trace metals as well as organic matter synthesized by the chemosynthetic microbial communities could both affect the benthic communities (Sarrazin et al., 2015; Sen et al., 2016).

A clear difference was found between the benthic communities at the locations influenced by the plume compared to the reference sites. The relative abundances of Arthropoda were higher at the plume affected sites, whereas the Annelida and Nematoda were more dominant in the reference areas. Also, the plume-affected sites NE and SE of the vent differed from each other. Especially at the sample locations NE of the vent site (i.e., the locations most influenced by the plume), the relative abundances of Annelida were low whereas the relative abundance of Annelida in the SE were slightly higher due to a higher relative abundance of order Terebellida (Annelida). Species within this order are predominantly deposit-feeders (Zhadan and Tzetlin, 2002; Jumars et al., 2015) and rely on organic matter deposited on the seafloor. As the SE sample locations were the nearest locations to the vent site, the influx of nutrients could possibly be higher compared to the other locations (Sarrazin et al., 2015; Boschen et al., 2016). Another striking difference seen between the plume-affected sites and the reference site is the relative read abundance of Harpacticoida. This high abundance of copepods relative to other species has been shown before at nutrient rich environments such as hydrothermal vents and cold seeps (Gollner et al., 2010; Plum et al., 2015). The faunal assemblages at the other sites within this cluster (54, 56, and 58) did significantly differ from the reference location as well as from the sample locations in the direct path of the plume which implies a spatial zonation away from the vent site. Peripheral benthic communities of an active vent have been shown to be distinct up to 200 m from the vent, forming a so called ‘transition-zone’ between the vent-communities and the regular deep-sea benthic communities (Sen et al., 2014; Boschen et al., 2016). Here, we show that this pattern might be extended up to several kilometers from the vent site. As hydrothermal vents are rare and often spatially isolated within the entire deep sea, the plume-affected communities could be a potentially unique community and thereby valuable in terms of overall benthic diversity.

A small group of genuinely benthic taxa could be detected from the water samples taken at several locations surrounding the hydrothermal vent area. Most invertebrate groups have both continuous and seasonal reproduction. Whereas groups like isopods, amphipods, cumacea, and tanaids are direct developers (brooding eggs), pelagic larvae are still common in the deep-sea, however, mainly lecithotrophic (non-feeding) (Gage and Tyler, 1991). Sea anemones (Actiniaria, Cnidaria) were detected both upstream and downstream of the vent area. The majority of species within this taxonomic group require hard substrata for their adult stages; a substrate found aplenty in the MAR area. Moreover, several species of sea anemones have been documented living on hydrothermal vent edifices (López-González et al., 2003). Actiniaria were not detected in the reference area although hard substrata were abundant in these areas as well. Their absence here was confirmed in video footage as acquired during the research cruises, and might be explained to a relatively lower abundance of Actiniaria in the reference areas, below the detection threshold, compared to the vent areas.

Most of the pelagic taxonomic groups detected are known to contain species living in the deep-sea. However, also a reasonable part of the OTUs were assigned to groups living in the epipelagic zone such as the Siphonophorae. The Siphonophorae consist of mainly pelagic species living in the upper few hundred meters of the water column, while fewer species within this group live in the deeper parts of the ocean (Mapstone, 2014). The relative abundance of reads for Siphonophorae in the water samples obtained with the LVS was close to zero, whereas the relative abundance for this taxonomic group was considerably higher for the water samples obtained with Nisking bottles, mounted on a CTD-carrousel. Many specimens of a pelagic member of this group, the Physalia physalis (Portuguese man o’war), a species bearing long sticking tentacles, were detected at the water surface during the sampling cruises and one was even found attached to the CTD frame as the sampling gear was recovered on deck. We assume the higher abundance of Siphonophorae in the Niskin-water samples was due to a difference in volume/surface ratio between the two sampling devices as a lower volume/surface ratio is prone to sampling these Physalia physalis as they might stick to the surface during the downcast and probably release vast amounts of eDNA within the sampled water. Therefore, the detection of these free swimming and floating species seems to be mainly a sampling artifact.

This study shows the spatial distinction of benthic communities surrounding the Rainbow hydrothermal vent. It highlights that distinctive and potentially valuable communities for the overall deep sea benthic fauna exist on a small geographic scale. Deep sea mining is postulated to create mining plumes due to resuspension and discharge (Boschen et al., 2013; Van Dover, 2014; Miller et al., 2018). These plumes can affect the soft-sediment background fauna in two ways: (1) via increased sedimentation (Boschen et al., 2013) or, (2) via toxification by (heavy) metals becoming bioavailable upon oxidation (Van Dover, 2011; Levin et al., 2016; Weaver et al., 2018). Biological diversity was lowest at locations closest to the vent site, possibly as a result of the toxicity of the vent effluent. Mining plumes will possibly disperse by ocean currents, and the area of impact on the communities might therefore be similar to the area affected by a natural continuous vent plume, depending on the settling rates of the particles. Due to the spatial variation of benthic communities that we found in the plume’s trajectory, conservation and restoration of all these communities in the event of deep-sea mining will be difficult. The actual effects of these plumes on the soft-sediment communities are yet unknown, but our results show that these communities are different under the plume than at reference sites, a strong indication that they are affected by the plume. Therefore, it cannot be excluded that the mining of SMS deposits and the thereby created plumes will affect these unique benthic communities.