An Assessment of Environmental Metabarcoding Protocols Aiming at Favoring Contemporary Biodiversity in Inventories of Deep-Sea Communities

The abyssal seafloor covers more than 50% of planet Earth and is a large reservoir of still mostly undescribed biodiversity. It is increasingly targeted by resource-extraction industries and yet is drastically understudied. In such remote and hard-to-access ecosystems, environmental DNA (eDNA) metabarcoding is a useful and efficient tool for studying biodiversity and implementing environmental impact assessments. Yet, eDNA analysis outcomes may be biased toward describing past rather than present communities as sediments contain both contemporary and ancient DNA. Using commercially available kits, we investigated the impacts of five molecular processing methods on eDNA metabarcoding biodiversity inventories targeting prokaryotes (16S), unicellular eukaryotes (18S-V4), and metazoans (18S-V1, COI). As the size distribution of ancient DNA is skewed toward small fragments, we evaluated the effect of removing short DNA fragments via size selection and ethanol reconcentration using eDNA extracted from 10 g of sediment at five deep-sea sites. We also compare communities revealed by eDNA and environmental RNA (eRNA) co-extracted from ∼2 g of sediment at the same sites. Results show that removing short DNA fragments does not affect alpha and beta diversity estimates in any of the biological compartments investigated. Results also confirm doubts regarding the possibility to better describe live communities using eRNA. With ribosomal loci, eRNA, while resolving similar spatial patterns than co-extracted eDNA, resulted in significantly higher richness estimates, supporting hypotheses of increased persistence of ribosomal RNA (rRNA) in the environment and unmeasured bias due to overabundance of rRNA and RNA release. With the mitochondrial locus, eRNA detected lower metazoan richness and resolved fewer spatial patterns than co-extracted eDNA, reflecting high messenger RNA lability. Results also highlight the importance of using large amounts of sediment (≥10 g) for accurately surveying eukaryotic diversity. We conclude that eDNA should be favored over eRNA for logistically realistic, repeatable, and reliable surveys and confirm that large sediment samples (≥10 g) deliver more complete and accurate assessments of benthic eukaryotic biodiversity and that increasing the number of biological rather than technical replicates is important to infer robust ecological patterns.

The abyssal seafloor covers more than 50% of planet Earth and is a large reservoir of still mostly undescribed biodiversity. It is increasingly targeted by resource-extraction industries and yet is drastically understudied. In such remote and hard-to-access ecosystems, environmental DNA (eDNA) metabarcoding is a useful and efficient tool for studying biodiversity and implementing environmental impact assessments. Yet, eDNA analysis outcomes may be biased toward describing past rather than present communities as sediments contain both contemporary and ancient DNA. Using commercially available kits, we investigated the impacts of five molecular processing methods on eDNA metabarcoding biodiversity inventories targeting prokaryotes (16S), unicellular eukaryotes (18S-V4), and metazoans (18S-V1, COI). As the size distribution of ancient DNA is skewed toward small fragments, we evaluated the effect of removing short DNA fragments via size selection and ethanol reconcentration using eDNA extracted from 10 g of sediment at five deep-sea sites. We also compare communities revealed by eDNA and environmental RNA (eRNA) co-extracted from ∼2 g of sediment at the same sites. Results show that removing short DNA fragments does not affect alpha and beta diversity estimates in any of the biological compartments investigated. Results also confirm doubts regarding the possibility to better describe live communities using eRNA. With ribosomal loci, eRNA, while resolving similar spatial patterns than co-extracted eDNA, resulted in significantly higher richness estimates, supporting hypotheses of increased persistence of ribosomal RNA (rRNA) in the environment and unmeasured bias due to overabundance of rRNA and RNA release. With the mitochondrial locus, eRNA detected lower metazoan richness and resolved fewer spatial patterns than co-extracted eDNA, reflecting high messenger RNA lability. Results also highlight the importance of using large amounts of sediment (≥10 g) for accurately surveying eukaryotic diversity. We conclude that eDNA should be favored over eRNA INTRODUCTION Environmental DNA (eDNA) metabarcoding is an increasingly used tool for biodiversity inventories and ecological surveys. Using high-throughput sequencing (HTS) and bioinformatic processing, it allows the detection or the inventory of target organisms using their DNA directly extracted from soil, water, or air samples (Taberlet et al., 2012a). As it does not require specimen isolation, it represents a practical and efficient tool in large and hard-to-access ecosystems, such as the marine realm. Besides allowing studying various biological compartments simultaneously, metabarcoding is also very effective for detecting diversity of small organisms (microorganisms, meiofauna) largely disregarded in visual biodiversity inventories due to the difficulty of their identification based on morphological features (Carugati et al., 2015).
The deep sea, covering more than 50% of Planet Earth, remains critically understudied, despite being increasingly impacted by anthropogenic activities and targeted by resourceextraction industries (Ramirez-Llodra et al., 2011). The abyssal seafloor is mostly composed of sedimentary habitats containing high numbers of small (< 1 mm) organisms, and characterized by high local and regional diversity (Grassle and Maciolek, 1992;Smith and Snelgrove, 2002). Given the increased time efficiency offered by eDNA metabarcoding and its wide taxonomic applicability, this tool is a good candidate for large-scale biodiversity surveys and environmental impact assessments (EIAs) in the deep-sea biome.
eDNA is a complex mixture of genomic DNA present in living cells, extra-organismal DNA, and extracellular DNA originating from the degradation of organic material and biological secretions (Torti et al., 2015). Extracellular DNA has been shown to be very abundant in marine sediments, representing 50-90% of the total DNA pool (Dell'Anno and Danovaro, 2005;Corinaldesi et al., 2018). However, this extracellular DNA compartment may not only contain DNA from contemporary communities. Indeed, nucleic acids can persist in marine sediments as their degradation rate decreases due to adsorption onto the sediment matrix (Corinaldesi et al., 2008;Torti et al., 2015). Low temperatures, high salt concentrations, and the absence of UV light are additional factors enhancing long-term archiving of DNA in deep-sea sediments (Torti et al., 2015;Nagler et al., 2018). Decreased rates of abiotic DNA decay can thus allow DNA persistence over millennial timescales. Indeed, up to 125,000-year-old ancient DNA (aDNA) has been reported in oxic and anoxic marine sediments at various depths (Boere et al., 2011;Coolen et al., 2013;Lejzerowicz et al., 2013a). As extracellular DNA fragment size depends on its state of degradation (Nagler et al., 2018 report overall size ranges from 80 to over 20,000 bp), aDNA fragments have generally been reported to be <1,000 bp long (Boere et al., 2011;Coolen et al., 2013;Lejzerowicz et al., 2013a;Lennon et al., 2018). Restricting molecular biodiversity assessments to large DNA fragments may thus allow avoiding the bias of aDNA in biodiversity assessments aiming at describing contemporary communities using eDNA metabarcoding.
Environmental RNA (eRNA) has been viewed as a way to avoid the problem of aDNA in eDNA biodiversity inventories because RNA is only produced by living organisms and quickly degrades when released in the environment due to spontaneous hydrolysis and the abundance of RNases (Torti et al., 2015). Few studies have investigated this in the deep-sea, with contrasting results. Investigating foraminiferal assemblages, Lejzerowicz et al. (2013b) found similar taxonomic compositions with DNA and RNA, although highlighting that RNA is more appropriate for targeting the active community component. Contrastingly, Guardiola et al. (2016) detected marked differences between RNA and DNA inventories for most eukaryotic groups but found that both biomolecules detected similar patterns of ecological differentiation, concluding that "dead" DNA did not blur patterns of community structure. Laroche et al. (2018Laroche et al. ( , 2017 found stronger responses to environmental impact in alpha diversity measured with eRNA, while eDNA was better at detecting effects on community composition. Finally, long-term archived and even fossil RNA were also reported in sediment and soil Cristescu, 2019), casting doubts as to its advantage over DNA to inventory contemporary biodiversity.
The design of a sound environmental metabarcoding protocol to inventory biodiversity on the deep seafloor relies on a better understanding of the potential influence of aDNA on the different taxonomic compartments targeted. Using commercially available kits based on 2 and 10 g of sediment, we studied samples from five deep-sea sites encompassing three different habitats and spanning wide geographic ranges in order to select an optimal protocol to survey contemporary benthic deep-sea communities spanning the tree of life. We analyze eDNA and eRNA extracts via metabarcoding, targeting the V4-V5 regions of the 16S ribosomal RNA (rRNA) barcode (Parada et al., 2016) for prokaryotes, the 18S-V4 rRNA barcode region for micro-eukaryotes (Stoeck et al., 2010), and the 18S-V1V2 rRNA (thereafter 18S-V1) and Cytochrome c Oxidase I (COI) barcode markers for metazoans (Leray et al., 2013;Sinniger et al., 2016).
Our objectives were threefold: (1) Evaluate the effect of removing short DNA fragments from DNA extracts obtained using a 10-g extraction kit, (2) Compare eDNA and eRNA inventories resulting from the same samples via a 2-g joint extraction kit, and (3) Assess the aforementioned kits in terms of repeatability and suitability for different taxonomic compartments.

Collection of Samples
Sediment cores were collected from five deep-sea sites from various habitats (mud volcano, seamounts, and an area close to hydrothermal vents; Supplementary Table S1).
Triplicate tube cores were collected with a multicorer or with a remotely operated vehicle at each sampling site. The sediment cores were sliced into layers, which were transferred into zip-lock bags, homogenized, and frozen at −80 • C onboard before being shipped on dry ice to the laboratory. The first layer (0-1 cm) was used for the present analysis. In each sampling series, an empty bag was kept as a field control processed through DNA extraction and sequencing.

Nucleic Acid Extractions and Molecular Treatments
eDNA With the 10-g PowerMax Kit DNA extractions were performed using ∼10 g of sediment with the PowerMax Soil DNA Isolation Kit (MO BIO Laboratories, Inc., Qiagen, Hilden, Germany). To increase the DNA yield, the elution buffer was left on the spin filter membrane for 10 min at room temperature before centrifugation. For field controls, the first solution of the kit was poured into the control zip lock before following the usual extraction steps. DNA extracts were stored at −80 • C.

Size Selection of eDNA Extracts
Size selection of total eDNA extracted as detailed above from ∼10 g of sediment was carried out to remove small DNA fragments. NucleoMag NGS Clean-up and Size Select beads (Macherey-Nagel, Düren, Germany) were used at a ratio of 0.5 × for removing DNA fragments <1,000 bp from 500 µl of extracted eDNA. The target fragments were eluted from the beads with 100 µl elution buffer, and successful size selection was verified by electrophoresis on an Agilent TapeStation using the Genomic DNA High ScreenTape kit (Agilent Technologies, Santa Clara, CA, United States).

Ethanol Reconcentration of eDNA Extracts
A 3.5 ml aliquot of eDNA extracted from ∼10 g of sediment was reconcentrated with 7 ml of 96% ethanol (EtOH) and 200 µl of 5 M sodium chloride (NaCl), according to the guidelines in the Hints and Troubleshooting Guide of the PowerMax Soil DNA Isolation Kit. As this protocol does not include any incubation time, it favors large DNA fragments. The DNA pellet was washed with 1 ml 70% EtOH, centrifuged again for 15 min at 2,500 × g, and air-dried before being resuspended in 450 µl elution buffer.

PCR Amplification and Sequencing
Nucleic acid extracts were normalized to 0.25 ng/µl, and 10 µl of standardized samples were used in PCR. Four primer pairs were used to amplify one mitochondrial and three rRNA barcode loci targeting metazoans (COI, 18S-V1), micro-eukaryotes (18S-V4), and prokaryotes (16S-V4V5 for homogeneity; Supplementary  Table S3). Two metazoan mock communities (detailed in Brandt et al., 2020) were included for 18S-V1 and COI. For each sample and marker, triplicate amplicon libraries (see Supporting Information for amplification details) were prepared by ligation of Illumina adapters on 100 ng of amplicons following the Kapa Hifi HotStart NGS Library Amplification Kit (Kapa Biosystems, Wilmington, MA, United States). After quantification and quality control, library concentrations were normalized to 10 nM, and 8-9 pM of each library containing a 20% PhiX spike-in were sequenced on a HiSeq2500 (System User Guide Part #15035786) instruments in a 250 bp paired-end mode.

Bioinformatic Analyses
All bioinformatic analyses were performed using a Unix shell script (Brandt et al., 2020), available on Gitlab 1 , on a home-based cluster (DATARMOR, Ifremer). The details of the pipeline, along with specific parameters used for all metabarcoding markers, are given in Supplementary Table S4 and in Brandt et al. (2020). Pairs of Illumina reads were corrected with DADA2 v.1.10 (Callahan et al., 2016) following the online tutorial for paired-end data 2 and delivered inventories of amplicon sequence variants (ASVs). Metazoan data were further clustered into operational taxonomic units (OTUs) with swarm v2, a singlelinkage clustering algorithm (Mahé et al., 2015) that aggregates sequences iteratively and locally around seed sequences based on d, the number of nucleotide differences, to determine coherent groups of sequences, independent of amplicon input order, allowing highly scalable and fine-scale clustering. ASVs were swarm clustered at d-values of 4 for 18S-V1 and 6 for COI, using the FROGS pipeline (Escudié et al., 2018).
We chose to evaluate micro-eukaryote and prokaryote diversity at the ASV level due to its increasing use in the literature (Callahan et al., 2017). Although the use of OTUs may also be justified for microbial diversity depending on study objectives (Brandt et al., 2020), we did not expect an alteration of alpha and beta diversity patterns between ASV and OTU levels for the different molecular treatments investigated. ASVs and OTUs were taxonomically assigned via BLAST + (v2.6.0) based on minimum similarity and minimum coverage (-perc_identity 70 and -qcov_hsp 80). For ASVs, sequences obtained with DADA2 were subsequently assigned with blastn. For OTUs, BLAST assignment in FROGS was performed using the affiliation_OTU.py command. It is not uncommon for deep-sea taxa to have closest relatives in databases (even congenerics) exhibiting nucleotide divergence exceeding 20% (Shank et al., 1999;Herrera et al., 2015). Considering our interest in diverse and poorly characterized communities, more stringent BLAST thresholds were thus not implemented at this stage. However, additional filters were performed during downstream bioinformatic processing described below, and taxonomic information was used at phylum level, only when the assignment was deemed reliable at this taxonomic level. The Silva132 reference database was used for taxonomic assignment of rRNA marker genes (Quast et al., 2012), and MIDORI-UNIQUE (Machida et al., 2017) was used for COI.
Molecular inventories were refined in R v.3.5.1 (R Core Team, 2018). A blank correction was made using the decontam package v.1.2.1 (Davis et al., 2018), removing all clusters that were more prevalent in negative control samples than in true or mock samples. Unassigned and non-target clusters were removed. Additionally, for metazoan loci, all clusters with a terrestrial assignment (groups known to be terrestrial-only) were removed. Samples with fewer than 10,000 target reads were discarded. We performed an abundance renormalization to remove spurious ASVs/OTUs due to random tag switching . The COI OTU table was further curated with LULU v.0.1 (Frøslev et al., 2017) to limit the bias due to pseudogenes, using a minimum co-occurrence of 0.93 and a minimum similarity threshold of 84%.

Statistical Analyses
Sequence tables were analyzed using R with the packages phyloseq v1.22.3 (McMurdie and Holmes, 2013), following guidelines in online tutorials 3 , and vegan v2.5.2 (Oksanen et al., 2018). Alpha diversity between molecular processing methods was estimated with the number of observed target clusters in rarefied datasets. Cluster abundances were compared via analyses of deviances (ANODEV) on generalized linear mixed models using negative binomial distributions, as the data were overdispersed. Pairwise post hoc comparisons were performed via Tukey honestly significant difference (HSD) tests using the emmeans package.
Differences in community compositions resulting from molecular processing were evaluated with Mantel tests (Jaccard and Bray-Curtis dissimilarities for metazoan and microbial taxa, respectively; Pearson's product-moment correlation; 1,000 permutations). Permutational multivariate analysis of variance (PERMANOVA) was performed on normalized datasets to evaluate the effect of molecular processing and site on community compositions using the function adonis2 (vegan) with Jaccard dissimilarities (presence/absence) for metazoan and Bray-Curtis dissimilarities for prokaryotes and microeukaryotes. The rationale behind this choice is that metazoans are multicellular organisms of extremely varying numbers of cells, organelles, or ribosomal repeats in their genomes, and can also be detected through a diversity of remains. The number of reads can thus not be expected to reflect the abundance of detected OTUs. Significance was evaluated via marginal effects of terms using 10,000 permutations with site as a blocking factor. Pairwise post hoc comparisons were performed via the pairwiseAdonis package, with site as a blocking factor. Differences between samples were visualized via principal coordinates analysis (PCoA) based on the abovementioned dissimilarities.
Finally, taxonomic compositions in terms of cluster and read abundance were compared between molecular processing methods. In order to compare accurately phylum-level taxonomic compositions, datasets were subsampled to clusters having a minimum hit identity of 86% for rRNA loci and 80% for COI. These values were chosen as they represent approximate minimum identity for reliable phylum assignment (Stefanni et al., 2018).

High-Throughput Sequencing Results
A total of 70 million 18S-V1 reads, 61 million COI reads, 30 million 18S-V4 reads, and 45 million 16S-V4V5 reads were obtained from four Illumina HiSeq runs of pooled amplicon libraries built from triplicate PCR replicates of 75 sediment samples, two mock communities (for 18S-V1 and COI), three extraction blanks, and two to four PCR negative controls (Supplementary Table S6). One to seven sediment samples failed amplification in each dataset. These were always coming from the same sampling sites (MDW-ST117 and MDW-ST38) and predominantly comprised RNA samples (Supplementary Table S6). After bioinformatic processing, read numbers were reduced to 44 million for 18S-V1, 45 million for COI, 16 million for 18S-V4, and 24 million for 16S-V4V5 (Supplementary Table S6). For eukaryote markers, fewer reads were retained in negative controls (2-64%) than in true or mock samples (49-83%), while the opposite was observed for prokaryotes with 16S-V4V5 (62% of reads retained in control samples against 49-57% in true samples). Negative control samples (extraction and PCR blanks) contained 0.001-0.6% of total processed reads compared to 1.3-1.5% in true samples. DNA extracts obtained from the joint DNA/RNA protocol based on the 2-g kit produced fewer eukaryotic reads than DNA extracts from the 10-g kit, while similar yields were obtained for prokaryotes. RNA extracts produced more reads than DNA extracts with the ribosomal loci, while they produced fewer reads with the mitochondrial COI locus (Supplementary Table S6).

Alpha Diversity Between Processing Methods
Rarefaction curves showed that a plateau was reached for all samples, suggesting an overall sequencing depth adequate to capture the diversity present (Supplementary Figure S1). Processing methods significantly affected the number of recovered eukaryote and prokaryote clusters, and significant variability among sites was detected for 18S-V1 for homogeneity and 18S-V4 (Table 1 and Supplementary Figure S2).
Molecular processing designed to remove small DNA fragments (i.e., size selection of DNA to remove fragment <1,000 bp and EtOH reconcentration) did not significantly affect recovered cluster numbers obtained from eDNA extracted from 10 g of sediment for any of the loci investigated (Figure 1 and Table 1; Tukey's HSD multiple comparisons tests, p > 0.9).
Extracts based on the 2-g kit resulted in more variability, reflected by greater standard errors in mean recovered cluster numbers (15-26% of the mean for eukaryotes, 7-9% for prokaryotes) than in DNA extracts based on 10 g of sediment (8-11% for eukaryotes, 3-6% for prokaryotes).

Effect of Molecular Processing Methods on Beta-Diversity Patterns
PERMANOVA showed that although site was the main source of variation among samples (accounting for 20-57% of variability), significant differences existed among molecular methods in terms of community structure for all loci investigated over and above any variation due to site ( Table 1). Pairwise comparisons indicated no significant effect of small DNA fragment removal on revealed community composition (Table 1), and high and significant correlations in Mantel tests (r: 0.92-1.0, p = 0.001) confirmed the minor effect of size selection and EtOH reconcentration. Based on these results, the size-selected and EtOH-reconcentrated DNA data were removed from further analyses, and community structures of the DNA 10-g extracts were compared with those derived from co-extracted DNA/RNA using the 2-g kit.
Pairwise comparisons showed significant differences in community structures between RNA and DNA for all markers analyzed ( Table 1). Ordinations confirmed the predominant effect of site as the first two PCoA axes mostly resolved spatial effects (Supplementary Figure S4) but also revealed that communities detected by RNA differed from those detected by DNA (co-extracted DNA and DNA 10-g), the level of differentiation varying among sites (Figure 3).
Pairwise comparisons also indicated significant differences in community structure between DNA extracts from the 2-g and 10g kits (Table 1) possibly due to higher variability among replicate cores in the DNA 2-g method as seen in ordinations (Figure 3).

Extraction Kit Versus Nature of Nucleic Acid
PERMANOVA of the dataset containing DNA 10-g, DNA 2-g, and RNA 2-g extracts confirmed that site was the predominant effect, explaining ∼20% of variation for metazoans, 33% of variation for micro-eukaryotes, and 54% of variation for prokaryotes. The analysis also indicated that the differences observed between processing methods were predominantly due to the type of nucleic acid rather than the kit used for extraction. Nucleic acid TABLE 1 | Changes in cluster richness and community structures with molecular processing method (DNA 10 g: DNA extracts from ∼10 g of sediment with the PowerMax Soil kit; DNA/RNA 2 g: DNA/RNA extracts from ∼2 g of sediment with the RNeasy PowerSoil kit) and site for the four studied genes.

Locus
Cluster richness Community differentiation ANODEVs were performed on mixed models with negative binomial distributions using rarefied datasets. PERMANOVAs were calculated on normalized datasets by permuting 10,000 times with Site as a blocking factor using Jaccard dissimilarities for 18S-V1 and COI and Bray-Curtis dissimilarities for 18S-V4 and 16S. Significant p-values are in bold. For pairwise comparisons, DNA 10 g comprises all processing methods based on DNA extracted from ∼10 g of sediment, and significance codes are ***p < 0.001; **p < 0.01; *p < 0.05. nature (DNA versus RNA) led to significant differences among assemblages for all loci, while DNA extraction kit resulted in significant differences only for 18S-V1 and 18S-V4 (Supplementary Table S7). This supported observations in relative taxonomic compositions, which were more similar between samples based on DNA (Figure 4), a pattern consistent across cores within each site (Supplementary Figure S5). Expectedly, when looking at read numbers, resolved taxonomic structures were also more similar among DNA-based methods (Supplementary Figure S6). Comparing read and cluster abundances revealed that relative taxonomic compositions based on read numbers (Supplementary Figure S6) were comparable to those based on cluster numbers (Figure 4) for microeukaryotes and prokaryotes and confirmed that this was not the case for metazoans. FIGURE 2 | Mean number of metazoan operational taxonomic units (OTUs) (COI, 18S-V1), protist amplicon sequence variants (ASVs) (18S-V4), and prokaryote ASVs (16S) detected per sample for each of the five processing methods (DNA 10 g: crude DNA extracts from ∼10 g of sediment with the PowerMax Soil kit; DNA 10 g EtOH rec. EtOH reconcentrated 10 g DNA extracts; DNA 10 g S-S: size-selected 10 g DNA extracts; DNA/RNA 2 g: crude DNA/RNA extracts from ∼2 g of sediment with the RNeasy PowerSoil kit). Cluster numbers were calculated on rarefied datasets. Error bars represent standard errors.

DISCUSSION
The aim of this study was to evaluate different molecular methods in order to select the most appropriate eDNA metabarcoding protocol to inventory contemporary deep-sea communities, with the lowest possible bias due to aDNA.
Using RNA rather than DNA to inventory contemporaneous communities has been suggested as a means of avoiding the bias due to long-term persistence of DNA in marine sediments. Indeed, RNA is only produced by living organisms and is thought to quickly degrade when released in the environment due to spontaneous hydrolysis and the abundance of RNases (Torti et al., 2015). Expectedly, in our COI dataset, RNA resulted in fewer OTUs (Figure 1) and detected fewer phyla (Figure 2) than co-extracted DNA. Contrastingly, for ribosomal loci, RNA detected higher cluster numbers than co-extracted DNA (Figure 1), resulting in more clusters per sample for most of the taxonomic groups detected (Figure 2). In these joint datasets, 45-63% of clusters were unique to RNA (Supplementary Figure S2). These unique clusters were not singleton clusters as only up to 2.2% of them had fewer than three reads, even if 5-28% had fewer than 10 reads (data not shown). Although proportions vary strongly among investigations, other studies using ribosomal loci have also reported increased recovery of OTUs in RNA datasets as well as considerable amounts of unshared OTUs between joint RNA and DNA data (Guardiola et al., 2016;Laroche et al., 2017, and references therein).
This difference observed here between COI and ribosomal loci is likely related to the nature of the targeted RNA molecule. The rapid hydrolysis of RNA mostly applies to random coils (like messenger RNA), while helical conformations (including most types of RNA, such as ribosomal RNA, transfer RNA, viral genomic RNA, or ribozymes) are less prone to hydrolysis by water molecules (Torti et al., 2015). The degradation of rRNA is thus likely to be much slower than that of messenger RNA (mRNA), which, combined with decreased digestion by RNases due to adsorption onto sediment particles (Torti et al., 2015), makes long-term persistence of rRNA possible and observed in sediments and even in fossils Cristescu, 2019). Finally, the great abundance of RNA over DNA in living organisms (e.g., 20.5% versus 3.1% in Escherichia coli) may also favor its persistence in the environment. This is especially true for rRNA, which is represented in a cell's RNA pool as many times FIGURE 3 | Principal coordinates analysis (PCoA) ordinations showing community differences between RNA and DNA molecular processing methods using either RNA/DNA extracted jointly from ∼2 g of sediment (RNA 2 g/DNA 2 g) or DNA extracted from ∼10 g of sediment (DNA 10 g) in five deep-sea sites using four barcode markers targeting metazoans (COI, 18S-V1), micro-eukaryotes (18S-V4), and prokaryotes (16S). PCoAs were calculated using Jaccard dissimilarities for metazoans and Bray-Curtis dissimilarities for unicellular organisms. Inserts show pairwise PCoAs.
as there are ribosomes, while only being present in a few copies (10-150) in the genome (Torti et al., 2015).
While RNA has been reported as an effective way to depict the active community compartment (Baldrian et al., 2012;Lejzerowicz et al., 2013b;Pawlowski et al., 2014), variation in activity levels between taxonomic groups as well as differences in life histories, life strategies, and non-growth activities may confound this interpretation and generate taxonomic bias (Blazewicz et al., 2013). Instead, DNA/RNA ratios might reflect different genomic architectures (variation in rDNA copy number) among taxonomic groups rather than different relative activities (Massana et al., 2015). Thus, eRNA data need to be interpreted with caution, as some molecular clusters could be overrepresented due to increased cellular activities . This could explain the higher cluster numbers detected here for ribosomal loci with eRNA compared to eDNA for several taxa (Figure 2).
Moreover, many of the unique RNA ASVs/OTUs may be artifacts from the reverse transcription of RNA to cDNA, a process known to generate errors that are difficult to measure and detect in bioinformatic analyses  but highlighted by the greater amounts of chimeras detected in RNA extracts with ribosomal loci (Supplementary Table S6). This overestimation of RNA-based data will affect non-clustered data more than clustered datasets, in line with the results observed here for microbial ASVs and metazoan OTUs.
In terms of beta diversity patterns, although RNA and DNA detected significantly different communities (Table 1), DNA and RNA samples resolved similar spatial configurations, with samples clustering by site (Figure 3). This is consistent with Guardiola et al. (2016), who also reported similar patterns of ecological differentiation between DNA and RNA in deep-sea sites, although both datasets resolved different communities. Although the comparative study performed here targeted only the first 1 cm layer of sediment, the comparable results obtained by Guardiola et al. (2016) on 5 cm suggest that these findings may be expanded to deeper layers of sediments. However, spatial variation was more pronounced with DNA samples for eukaryotes, which is congruent with Laroche et al. (2017), who suggested that eDNA may be more reliable for assessing differences in community composition.
Thus, due to its suspected persistence in the environment and the unknown but potentially additional sources of bias suspected here, using eRNA for metabarcoding of deep-sea sediments does FIGURE 4 | Patterns of relative cluster abundance resolved by metabarcoding of sediment RNA and DNA from five deep-sea sites using either RNA/DNA extracted jointly from ∼2 g of sediment (RNA 2 g/DNA 2 g) or DNA extracted from ∼10 g of sediment (DNA 10 g) and using four barcode markers targeting metazoans (A: COI, 18S-V1), micro-eukaryotes (B: 18S-V4), and prokaryotes (B: 16S). Values were calculated on balanced datasets. not seem to effectively address the problem of aDNA, and even less so for ribosomal loci. Other studies suggested that a more efficient way to deal with aDNA may be to use joint RNA and DNA datasets and trim for shared OTUs Pochon et al., 2017). This is however particularly stringent (given the low shared OTU proportions observed in this and other studies) and may result in a substantial number of false negatives. With COI, while mRNA may be more effectively targeting living organisms, the approach remains confronted with the taxonomic bias mentioned above, combined with higher in vitro lability of mRNA, making it more challenging to work with (highlighted by the increased failure of RNA extracts in this study; Supplementary Table S6).
Removing small DNA fragments via size selection (removing fragments < 1,000 bp) or EtOH reconcentration did not affect recovered cluster numbers in any of the biological compartments investigated (Figure 1). The methods also did not result in any significant difference in community structures ( Table 1), suggesting that small, likely ancient, DNA fragments have a negligible impact on biodiversity inventories produced through eDNA metabarcoding. This finding is in line with results from the deep-sea (Guardiola et al., 2016;Ramírez et al., 2018) and various other habitats (Lennon et al., 2018), which showed no evidence that spatial patterns were blurred by "dead" DNA persistence, and suggested a minimal effect of extracellular DNA on estimates of taxonomic and phylogenetic diversity.
None of the methods evaluated in the present study removes DNA not enclosed in living cells (e.g., DNA in organelles, DNA from dead cells. . .). It is still unclear how long DNA can remain intracellular after cell death or within organelles. Future research quantifying the rate at which "dead" intracellular DNA becomes extracellular and degraded, and investigation of deeper layers of sediment, will be valuable to estimate the potential bias of archived intracellular DNA in eDNA metabarcoding inventories of extant communities. However, there is increasing evidence that DNA from non-living cells is mostly contemporary (Lennon et al., 2018). This ability to detect extant taxa that were not present in the sample at the time of collection highlights the capacity of eDNA metabarcoding to detect local presence of organisms even from their remains or excretions, and even with a small amount of environmental material.
It remains to be elucidated whether more cost-and timeeffective extraction protocols specifically targeting extracellular DNA offer similar ecological resolution as total DNA kits. This is suggested to be the case for terrestrials soils (Taberlet et al., 2012b;Zinger et al., 2016), although authors have highlighted that conclusions from these studies should be interpreted with caution as results might be influenced by actively released and ancient DNA (Nagler et al., 2018). The only available study testing this in the deep-sea showed that richness patterns were strikingly different in several metazoan phyla between extracellular DNA and total DNA. The authors suggested this to be the result of activity bias: sponges and cnidarians were overrepresented in the extracellular DNA pool because they continuously expel DNA, while nematodes were underrepresented as their cuticles shield DNA (Guardiola et al., 2016). As this comparison was performed on samples collected in two consecutive years, differences observed may partly result from temporal variation. However, another study of shallow and mesobenthic macroinvertebrates showed that targeting solely the extracellular eDNA compartment of marine sediments led to the detection of more than 100 taxa fewer than bulk metabarcoding or morphology, suggesting that extracellular DNA may not be adequate for marine sediments (Aylagas et al., 2016).
Larger amounts of sediment (≥10 g) allowed detecting significantly more eukaryotic clusters. This was not true for prokaryotes, for which both ∼2 and ∼10 g of sediment detected similar numbers of ASVs (Table 1 and Figure 1). It may be suggested that in the joint RNA/DNA kit, DNA elution occurring after RNA elution induces partial DNA loss. However, such effect would be expected to equally affect eu-and prokaryotes, which was not the case here, supporting the fact that the quantity of the starting material significantly affects results for eukaryotes. The importance of adjusting the amount of starting material to the biological compartment investigated has already been documented (Creer et al., 2016;Dopheide et al., 2019), and this study confirms that while 2-5 g of deep-sea sediment may be enough to capture prokaryote diversity, microbial eukaryotes and metazoans are more effectively surveyed with larger sediment volumes.
Finally, the ∼2-g protocols were generally associated with higher variability among replicate cores for all loci investigated (Figures 1, 3). This variability increases confidence intervals, reduces statistical power, and increases the risk of not identifying differences among communities, and thus impacts in EIA studies (Type II errors). Small-scale (centimeters to meters) patchiness has often been reported in the deep-sea (Grassle and Maciolek, 1992;Smith and Snelgrove, 2002;Lejzerowicz et al., 2014). While technical (PCR) replicates allow increasing taxon detection probability (decrease false positives), this within-site variability can only be mitigated by collecting more biological replicates per sampling station and using a sufficiently high amount of starting material to extract nucleic acids.

DATA AVAILABILITY STATEMENT
The data for this work has been submitted to the European Nucleotide Archive (ENA) under the following project: PRJEB33873. Please refer to the sample metadata in Supplementary Material to download samples from ENA.

AUTHOR CONTRIBUTIONS
MB, CL-H, and SA-H designed the study. MB, JP, and CL-H carried out the laboratory work. MB and BT performed the bioinformatic analyses. MB, BT, and NH performed the statistical analyses. MB and SA-H wrote the manuscript. All authors contributed to the final manuscript.

FUNDING
This work was part of the "Pourquoi Pas les Abysses?" project funded by Ifremer and the project eDNAbyss (AP2016-228) funded by France Génomique (ANR-10-INBS-09) and Genoscope-CEA. This work also received funding from the European Union's Horizon 2020 Research and Innovation Program under grant agreement no. 678760 (ATLAS), and the 'Investissements d'Avenir' program OCEANOMICS (ANR-11-BTBR-0008; NH and CV). This output reflects only the authors' view, and the European Union cannot be held responsible for any use that may be made of the information contained therein.

ACKNOWLEDGMENTS
We wish to thank Laure Quintric, Patrick Durand, Caroline Belser, and Stéphane Pesant for bioinformatic and archiving support, as well as Jan Pawlowski and Eva Ramirez Llodra for useful advice on this work. We also wish to express our gratitude to the crew, participants, and mission chiefs of the MarMine cruise (Eva Ramirez Llodra, project 247626/O30, funded by the Research Council of Norway and associated industrial partners) and the MEDWAVES cruise supported by the ATLAS project and the Spanish Ministry of Economy, Industry, and Competitivity (Covadonga Orejas and all the crew from the Sarmiento de Gamboa), as well as to all the people who helped collecting samples (Perregrino Cambeiro, Juancho Movilla, Maria Rakka, Joana Boavida, Anna Addamo). This manuscript has been released as a Pre-Print at https://doi.org/10.1101/836080.
Frontiers in Marine Science | www.frontiersin.org