Metabarcoding monitoring analysis: the pros and cons of using co-extracted environmental DNA and RNA data to assess offshore oil production impacts on benthic communities

Field sampling

The study area was located 80 km off the west coast of New Zealand’s North Island in the South Taranaki Bight (latitude: −39.97295°, longitude: 173.2997°; Fig. 1A). In this region, the seabed has a fairly uniform muddy texture with water depth ranging between 100 and 125 m. Biological samples (bacterial, micro-eukaryotic and macro-faunal assemblages) and physico-chemical data were collected along transects radiating to the East–West axis of one Floating Production Storage and Offloading (FPSO) platform. The FPSO has been in activity since 2009 and has been processing and storing the well fluids coming from a wellhead platform (WHP) situated 1.3 km North-West. The main discharges from the FPSO are production waters (formation water and additives) and operational discharges (sewage, deck drainage, etc.; (McConnell et al., 2014)). The FPSO site is characterized by a prevalence of anthropogenic debris (e.g., paint chips, construction material, etc.) and coarser, hard-packed and dark-coloured sediments in its vicinity (≤250 m), possibly indicating lower oxygenation level, elevated barium concentrations radiating away from the production water and drilling discharge point, and differences in the structure of macro-faunal assemblages on the east transect, between the near field station (≤250 m) and far-field stations (>250 m; Johnston, Elvines & Newton, 2015).

The sampling methodology was based on the Offshore Taranaki Environmental Monitoring Protocol (OTEMP; Johnston et al., 2014), which consisted of a distance-graded sample station allocation. Sampling stations were overlaid along the major flow axis with the FPSO located at the center (Fig. 1B). The East–West axis constitutes the main trajectory along which the deposition of drilling mud and cuttings occurs, with the strongest currents usually flowing in a westward direction at a mean speed of 0.12 ms⁻¹ (McConnell et al., 2014). Six stations were sampled at approximately 250, 500 and 1,000 m from the FPSO.

Sediment samples were collected in triplicate at 6 stations for a total of 18 samples, on the 15th and 16th of March 2015 (Table S1). Samples were collected with a modified stainless steel double van Veen grab as described in Laroche et al. (2016). To avoid creating a bow wake effect on the seafloor, the grab sampler was deployed at a constant rate of 0.3 m/s. Recovery of the sampler was similarly processed to avoid displacing superficial sediments. Once retrieved, the top sediment layer was inspected and samples only taken when the surface was undisturbed. The full content of the first grab compartment was sieved (500 µm) and preserved in 70% ethanol for analysis of macro-faunal communities. Subsamples (2 g) of undisturbed surface sediment (approximately 1 cm depth) were collected from the second grab for eDNA/eRNA metabarcoding. These samples were placed in Life Guard™ Soil Preservation Solution (5 mL; MoBio, New York, NY, USA) to preserve microbial and environmental RNA in stasis, using disposable gloves and spatulas. Samples were stored on ice during transportation to the laboratory and kept frozen (−20 °C) until further processing. The second grab compartment was subsampled for physicochemical analyses of sediment as described in Johnston, Elvines & Newton (2015).

Laboratory analyses

Data analysis and statistics

Samples were grouped into near-field (≤250 m of the FPSO) and far-field (>250 m) stations from each side of the FPSO. To visualize the effect of trimming by singletons and by shared OTUs, the number of reads and OTUs in each dataset were compared. The OTUs removed from the shared OTUs trimming method were further investigated by displaying bar plots of the relative abundance of the ten most important classes for bacteria and phyla for eukaryotes using the R package ‘phyloseq’ (McMurdie & Holmes, 2013), and by performing a Mann–Whitney U test with 999 permutations on these taxa to detect significant abundance differences between the near- and far-field stations.

Prior to further analysis, read counts were normalized among samples. The large discrepancy in read abundance among samples (up to five times) prevented the use of typical normalization techniques such as variance stabilization (R package DESeq2; Love, Huber & Anders, 2014) or cumulative sum scaling (R package metagenomeSeq; Paulson et al., 2013), and a subsampling-based normalization strategy was employed following the method described in Aguirre de Carcer et al. (2011). Briefly, multiple rarefactions (100 times) based on the median number of reads among all samples were performed, followed by a rarefaction to the lowest number of reads among all samples. Since a large number of singletons may represent PCR or sequencing artefacts (Majaneva et al., 2015; Brown et al., 2015), these were removed prior to statistical analysis.

Using normalized data kept after the two different trimming techniques, observed OTUs and the Shannon diversity index were analyzed for each sample. To test for the effect of the FPSO on the bacterial and eukaryotic alpha-diversity of eDNA/eRNA datasets, a nonparametric t-test using 999 permutations was performed between near and far-field stations. To determine whether the eDNA or eRNA dataset contained the most abundant number of OTUs significantly affected by the FPSO activity, a non-parametric t-test with 999 permutations was conducted on OTUs present in at least three samples, between the near- and far-field station groups. Differences in beta-diversity between and within the near- and far-field stations were tested using an Adonis and PERMDISP tests with 999 permutations on Bray-Curtis distance matrices, using square root transformed data. Non-metric multidimensional scaling plots were obtained using the R package ‘phyloseq’ (McMurdie & Holmes, 2013). Correlations between bacterial and eukaryotic eDNA/eRNA OTUs and macro-fauna beta-diversities, and with environmental variables were investigated using a Mantel test based on 999 permutations. The tests were performed on square root transformed data using Bray-Curtis distance matrices for biological data, and an Euclidean distance matrix for environmental data. For the later, only variables showing a significant change between the near-field and far-field stations were used (i.e., barium and AFDW; Johnston, Elvines & Newton, 2015), and square root transformed prior to the distance matrix computation. Except for the bar plots and nMDS, all analyses mentioned above were computed using the QIIME pipeline (Caporaso et al., 2010).

High-throughput sequencing output

Each of the 18 sediment samples collected in this study were co-extracted for both eDNA and eRNA and PCR-amplified for total bacteria and eukaryotes, resulting in the total production of 72 PCR-amplicons (36 for bacteria and 36 for eukaryotes; Table S1). The no template PCR controls remained negative in all PCR trials. All 72 amplicons were sequenced on an Illumina™ MiSeq instrument, yielding a total of 2,134,982 and 1,040,624 eDNA and eRNA sequence reads, respectively. The water blank negative controls (ddH₂O sample) contained a total of 524 raw reads (130 for bacteria and 494 for eukaryotes) in the eDNA library, and 405 raw reads (234 for bacteria and 171 for eukaryotes) in the eRNA library, of which none remained after quality filtering. The positive controls were assigned at 99.9% to the three target taxa. After removal of unmatched and low-quality reads, a total of 883,155 reads and 33,995 OTUs were obtained for bacteria (eDNA: 25,100; eRNA: 22,211), and a total of 768,374 reads and 11,854 OTUs (eDNA: 7,045; eRNA: 8,059) for eukaryotes. Removal of samples with insufficient diversity coverage led to the removal of one bacterial eDNA and eRNA sample (FE250b), and one eukaryotic eDNA (FE250b) and eRNA (FE250c) sample (Table S1).

Effect of trimming methodologies

For the bacterial dataset, the removal of singletons (i.e., OTUs present only once across the entire eDNA/eRNA datasets) resulted in a 1% loss of total reads and 37% loss of OTUs (Table 1). Similar results were observed for eukaryotes where 0.4% of reads and 29% of OTUs were removed. Trimming OTUs that did not occur in both eDNA/eRNA amplicons of the same sample resulted in a total reduction in OTUs of 78%, and total read reduction of 27%. For eukaryote data, trimming by shared OTUs had a larger impact on OTU and read removal, discarding an additional 51% of OTUs for a total reduction of 80%, and reducing the total read abundance by 28%. The mean number of reads (counts) per removed OTUs was of 3.1 and 8.6 for bacteria and eukaryotes respectively.

Figure 2 summarizes the total number of bacterial and eukaryotic OTUs that were either shared between both eDNA and eRNA data of the whole datasets, shared between eDNA and eRNA from the same sample, or those that were restricted to one dataset type following the two distinct data treatments. The percentage of shared OTUs between the eDNA and eRNA bacterial datasets increased from 22% to 36% after singleton removal (Table 1). For eukaryotes, removing singletons improved the percentage of shared OTUs in the eDNA and eRNA datasets from 20 to 28%. The percentage of unshared OTUs were highest in the eDNA dataset for bacteria, while in the eukaryotic dataset most of the unshared OTUs were in the eRNA.

After trimming by singletons, the bacterial datasets (kept OTUs) were mainly composed of Gammaproteobacteria (eDNA 57.2%, eRNA 77.9%), Deltaproteobacteria (eDNA 12%, eRNA 8%) and Alphaproteobacteria (eDNA 8.5%, eRNA 4.8%; Table S3). Metazoa composed the large majority of eukaryote OTUs (eDNA 84.6%, eRNA 85.5%), followed by Dinophyta (6.2%) and Radiolaria (2.8%) in eDNA, and by Cercozoa (7.4%) and Cilliophora (3.5%) in eRNA (Table S3).

Trimming by shared OTUs slightly changed the relative abundance of each taxa of the kept data, but did not change their decreasing order (Table S4 ). The majority of bacterial OTUs removed from trimming by shared OTUs were Proteobacteria belonging to Deltaproteobacteria (eDNA 19%, eRNA 30.4%) and Gammaproteobacteria (eDNA 22.7%, eRNA 28.2%; Fig. 3, Table S3). In eDNA, the taxonomic groups Planctomycetia (9.6%), Alphaproteobacteria (7.1%) and Verrucomicrobiae (5.7%) also accounted for a substantial part of the trimmed OTUs. Most of the eukaryote OTUs removed were Metazoa (eDNA 60.5%, eRNA 42.8%), followed by Dinophyta (13.9%) and Apicomplexa (10.1%) in eDNA, and by Ciliophora (19.7%) and Cercozoa (19.7%) in eRNA (Table S3). Among all these groups, only eDNA Verrucomicrobiae (T-stat = 2.39, P-value = 0.03; Table S6) and eRNA Planctomycetia (T-stat = −2.78, P-value = 0.03; Table S6) expressed a significant relation with distance from the platform, both being more abundant in the near-field stations.

Comparison of bacteria and eukaryotes eDNA and eRNA datasets

Significant OTUs, beta-diversity and correlation to environmental variables

The non-parametric t-test analysis identified OTUs significantly associated with the near- or far-field station groups, by comparing their frequencies across samples groups. For both taxa, these OTUs were more abundant in eDNA than eRNA (Table S7). Trimming by shared OTUs reduced the number of significant OTUs by 52 (eDNA) and 26 (eRNA) in bacteria, and by 99 (eDNA) and 86 (eRNA) in eukaryotes (Table S7).

Non-metric multidimensional scaling plots (Fig. 4) demonstrated that near- and far-field stations tended to cluster separately for both taxa. Our distance-based test for homogeneity of multivariate dispersion (Permdisp) indicated that communities within the near- and far-field stations groups were homogeneous in all datasets (Table S8). However, significant variance between replicates were observed for bacterial eDNA (p-value = 0.03) and eRNA (p-value = 0.04) data of the trimmed by singleton dataset, and for bacterial eDNA (p-value = 0.03) of the trimmed by shared OTUs dataset (Table S9).

Table 3:

Adonis and Mantel tests.

Analysis of the strength and statistical significance of sample groupings among datasets, and correlations between beta-diversity matrices of micro (bacteria and eukaryotes) and macro-fauna, and with distance matrix of selected environmental variables.

Datasets			Adonis		Mantel test
			Near-field VS Far-field		Correlation with macro-fauna		Correlation with env. variables		Correlation DNA/RNA
			R2	p-value	r	p-value	r	p-value	r	p-value
Bacteria	Trimmed by singletons	eDNA	0.07	0.01	−0.53	<0.01	0.52	<0.01	−0.05	0.80
		eRNA	0.07	0.21	−0.06	0.78	0.04	0.82
	Trimmed by shared OTUs	eDNA	0.10	<0.01	−0.51	<0.01	0.58	<0.01	0.31	0.08
		eRNA	0.07	0.01	−0.08	0.73	0.08	0.63
Eukaryotes	Trimmed by singletons	eDNA	0.08	0.04	−0.35	0.05	0.52	<0.01	0.47	0.02
		eRNA	0.12	<0.01	−0.14	0.48	0.24	0.09
	Trimmed by shared OTUs	eDNA	0.10	0.02	<0.01	0.99	0.32	0.02	0.89	<0.01
		eRNA	0.11	<0.01	−0.09	0.67	0.18	0.30

DOI: 10.7717/peerj.3347/table-3

Notes:

eDNA

environmental DNA

eRNA

environmental RNA

Significant p-values are in bold.

When comparing the correlation of bacterial eDNA/eRNA beta-diversity matrices between near- and far-field stations, eDNA provided stronger correlations and/or lower p-values (Table 3). This trend was also observed with correlations to the distance matrices of macro-fauna and environmental variables (Table 3). For eukaryotes however, this tendency was not always observed as eRNA showed slightly stronger correlations and more significant results between the near- and far-field stations, but weaker correlations with macro-faunal assemblages and with environmental variables. Overall, correlations between eDNA and eRNA beta-diversities were stronger for the trimmed by shared OTUs datasets, and only significantly positive among the eukaryote datasets (Table 3). For bacteria, the trimmed by singletons datasets demonstrated particularly weak correspondence (r =  − 0.05; p-value = 0.8) between eDNA and eRNA.

High-throughput sequencing outputs

In our datasets, the number of bacterial eDNA OTUs (Table 1) was slightly higher than eRNA. Environmental DNA is more resistant to posthumous degradation, and may persist several years in marine sediment (Dell’Anno, 2005). Analysis of the eukaryote data showed a differing pattern with more eRNA OTUs compared to eDNA. The proportion of OTUs found only in eDNA or eRNA datasets varies substantially among previous investigations. For example, Dowle et al. (2015) found more bacterial eRNA than eDNA OTUs in their datasets, and Hu et al. (2016) and Laroche et al. (2016) obtained more protist eRNA than eDNA OTUs, while Pawlowski et al. (2014) and Pochon et al. (2015) detected more protist eDNA than eRNA OTUs. It is possible that many of the unshared eRNA OTUs represent PCR artefacts, especially from the reverse transcription of RNA to cDNA, and sequencing errors (Egge et al., 2013; Ficetola et al., 2015). The RNA conversion to complementary DNA (cDNA) requires the use of a reverse transcriptase which, as opposed to DNA polymerase, lack proof reading activity, creating point mutations in some of the cDNA sequences (Svarovskaia et al., 2003; Houseley & Tollervey, 2010). Reverse transcriptase also frequently ‘jump’ from one template to another during the transcription process. Known as template-switching, this action can lead to two different outcomes: (1) production of chimeric cDNA sequences from intermolecular template switching, and (2) creation of shortened isoform sequences from intramolecular template switching (Cocquet et al., 2006); the latter being more problematic as it is unlikely to be detected during bioinformatics analysis. Finally, the use of random hexamer primers for cDNA synthesis is known to introduce nucleotide bias at the beginning of the 5′-end of sequence reads (Hansen, Brenner & Dudoit, 2010). Because the pre-mentioned artefacts are subsequently augmented by the gene-specific PCR amplification, they would become integrated into the HTS library and potentially not be removed during bioinformatics filtering processes. The incertitude linked to potential intrinsic error rates from reverse transcription limits accurate comparison of eDNA and eRNA datasets. The use of a RNA control (e.g., synthetic oligomers) and technical (PCR) replicates could help identifying these artefacts and improve concordance between the eDNA and eRNA profiles. Trimming OTUs by their shared presence in eDNA and eRNA represents an alternative way of discarding reverse transcription artefacts.

In this study, bacterial and eukaryote RNAs were treated using the same protocol during the reverse transcription process. However, eukaryote cDNA required more PCR cycles (35 instead of 30 for bacteria) to amplify sufficient material, potentially producing more false transcripts and explaining the differences observed between the two taxonomic groups (Fig. 2). Another possibility could be the presence of rare but highly active eukaryote OTUs which were not detected in the eDNA dataset because of their low abundance. Particularly high activity and growth from some of the rare biota has been previously reported for both bacteria (Campbell et al., 2011; Hugoni et al., 2013) and microbial eukaryotes (Logares et al., 2014).

Effect of trimming methodologies

Trimming OTUs not simultaneously present in both eDNA and eRNA assemblages from the same samples resulted in twice as many OTUs being removed compared to discarding singletons only, and removed 27-fold (bacteria) and 70-fold (eukaryotes) more reads. However, further investigation suggested that it mostly resulted in the removal of legacy DNA, DNA deposited from the water column, and technical artefacts. This improved the sensitivity and accuracy of our environmental metabarcoding analysis (Tables 2 and 3) and enabled the detection of differences between the effects of the production platform on near- and far-field stations. In particular, it allowed us to disclose the divergence between beta-diversities of near-field and far-field stations, and improved the correlations between the beta-diversity and environmental variables distance matrices (Table 3).

Removal of single sequence detection (singletons) to reduce PCR and sequencing artefacts is a controversial practice. Although many of these detections are true errors (Parada, Needham & Fuhrman, 2016), some may constitute real OTU occurrences. In studies not particularly interested by the rare biota, the trade-off of producing false negatives to remove false positives is usually accepted (Majaneva et al., 2015; Bálint et al., 2016). Systematic inclusion or exclusion of singletons is arbitrary and not necessarily adapted to all situations. For example, the likelihood of retaining false positives increases with the sequencing depth of HTS libraries (Ficetola, Taberlet & Coissac, 2016). Some researchers have therefore suggested the use of a read detection threshold based on sequencing depth (Bálint et al., 2016). However, Ficetola, Taberlet & Coissac (2016) and Lahoz-Monfort, Guillera-Arroita & Tingley (2016) strongly discouraged systematic removal of taxa based on an abundance threshold, unless high detection probability of real OTUs and sufficient biological or technical replicates are available to identify false positives. Instead, they propose the application of site occupancy-detection modelling (SODM) methods which take into account false positive and negative detection probabilities, and which have proven to work well with few replicate samples (Schmidt et al., 2013; Lahoz-Monfort, Guillera-Arroita & Tingley, 2016).

In our study, we used co-extracted DNA and RNA to segregate real OTUs from artefacts. This approach allows for an informed elimination of artefactual OTUs while also discarding data from legacy DNA. However, trimming data by shared OTUs between co-extracted eDNA and eRNA material is particularly stringent and may result in the removal of a substantial number of false negatives. As an alternative, biological or technical replicates could be used in parallel to discriminate genuine and false detections. For example, biological replicates could be pooled together prior to trimming by shared OTUs, and subsequently demultiplexed to allow evaluation of within station variability. It is worth mentioning that the application of unsupervised oligotyping methodologies such as the one proposed by Eren et al. (2015), could further improve sensitivity of environmental metabarcoding monitoring studies by discriminating closely related taxa, and will likely gain in popularity as tools are being developed.

Operational taxonomic units not shared between both eDNA and eRNA datasets

Among the bacterial datasets, most of the removed OTUs (i.e., those not simultaneously found in both eDNA and eRNA) belonged to Gammaproteobacteria and Deltaproteobacteria, which also constituted the main classes of the singleton trimmed dataset. This suggests that most OTUs removed by the shared OTUs trimming method may have arisen from PCR or sequencing errors originating from highly abundant unique sequences which subsequently associated to these classes (Huse et al., 2010; Quince et al., 2011). Another possibility could be that these abundant taxa have more polymorphisms. In such cases, removal of these rare variants would produce false negatives. Nonetheless, their small contribution to weighted (abundance-based) analyses, as demonstrated in the present study, would unlikely affect results.

Most of the eukaryotes which were removed during the trimming procedure were Metazoa, Dinophyta and Apicomplexa in eDNA, and Ciliophora and Cercozoa in eRNA (Fig. 3, Table S5). These phyla are also the most abundant in the datasets trimmed by singletons (Table S3). None of these showed a significant correlation with distance from the platform (Table S6). Just as for bacteria, the important abundance of these phyla in the discarded eDNA/eRNA OTUs is likely a consequence of PCR/sequencing errors originating from the most abundant unique sequences. This is further supported by the fact that the decreasing order of relative abundance of the most important taxa is the same between the datasets of kept OTUs and removed OTUs resulting from the shared OTUs trimming.

Comparison between eDNA and eRNA

When performing alpha-diversity on the trimmed by singletons datasets, differences were observed between eDNA and eRNA, with the latter mainly showing stronger dissimilarities between the near-field and far-field stations. This is not surprising as eDNA can contain genetic information of past and current communities, inflating alpha-diversity and potentially masking the present effects of the platform (Corinaldesi, Danovaro & Dell’Anno, 2005; Dell’Anno, 2005).

Beta-diversity analysis of the singleton trimmed data showed more significant differences between the near- and far-field stations with bacterial eDNA, while only small differences could be observed between eukaryote eDNA and eRNA. Relationships between micro-eukaryotic and the macro-faunal assemblages and environmental variables were also much stronger than with eDNA. This suggest that while eRNA may represent a better proxy for assessing ongoing anthropogenic effects on species diversity, eDNA seems more reliable for assessing effects on community composition. This is especially true for bacteria where no correlation could be observed between eRNA data and macro-fauna, or with environmental variables (Table 2). These results contrast with other similar studies (Pawlowski et al., 2014; Visco et al., 2015; Dowle et al., 2015; Pochon et al., 2015; Laroche et al., 2016; Hu et al., 2016) where eRNA always out-performed eDNA when correlating microbial and macro-faunal assemblages, or when correlating these assemblages with environmental conditions. Nonetheless, these studies used a less stringent OTU filtering method prior to statistical analyses, removing only those not present in both eDNA/eRNA compartment of their whole dataset. Consequently, an important number of OTUs from legacy DNA may have subsisted in their datasets, decreasing the observed sensitivity to local conditions of the eDNA-derived communities. Also, the relationship between eRNA concentrations and activity does not always align. Differences such as life history, life strategies and non-growth activity between organisms play an important role in RNA production (Blazewicz et al., 2013). Ribosomal eRNA such as those targeted in this study (16S and 18S) can have a much wider range of expression (over 3 orders of magnitude) than ribosomal eDNA (ca. 1 order of magnitude; (Fegatella et al., 1998)). Fluctuation in DNA copy number between taxa and in RNA expression may account for the strong variability observed among the bacterial eDNA and eRNA replicates (Table S8). This variability increases confidence intervals, reduces statistical power, and increases the risk of not identifying impacts (Type II errors; Underwood, 1993).

For accurate assessment of false positive and negative OTU occurrences, and sufficient evaluation of within station biological variability, optimal number of replicate samples is often a compromise between available funds and time, and the desired level of precision and sensitivity. Parameters such as targeted taxa, the choice of primer set and the sequencing depth of HTS libraries all play an important part in taxa detection probability, which in turn, dictates the optimal replication level (Ficetola et al., 2015). In studies assessing whole communities, the number of replicates should ideally depend on the taxa with the lowest detection probability. Ficetola et al. (2015) give an example where a taxon with a 0.5 probability would normally require a minimum of eight replicates to ensure detection. For routine monitoring surveys of benthic communities however, where the non-detection of rare biota is unlikely to change the outcomes of studies, the replication level does not have to be as high. In such cases, the number of replicates may be more influenced by the level of heterogeneity found within stations.

An important consideration when using single-time sampling designs over a contamination gradient is that it assumes that natural factors equally influence communities among samples (Wiens & Parker, 1995). This bias can be particularly problematic with small and short-lived organisms like bacteria. In low impacted environments, such as those explored in this study, biological interactions may supplant the effect of environmental factors, increasing spatial variability (Schafer, 2000), and creating stronger discrepancies among biological replicates. In this situation, collecting more biological replicates per stations should be favoured over technical replicates, in order to mitigate the effect of micro-patchiness.