Impact of preservation method and storage period on ribosomal metabarcoding of marine microbes: Implications for remote automated samplings.

Automated sampling technologies can enhance the temporal and spatial resolution of marine microbial observations, particularly in remote and inaccessible areas. A critical aspect of automated microbiome sampling is the preservation of nucleic acids over long-term autosampler deployments. Understanding the impact of preservation method on microbial metabarcoding is essential for implementing genomic observatories into existing infrastructure, and for establishing best practices for the regional and global synthesis of data. The present study evaluates the effect of two preservatives commonly used in autosampler deployments (mercuric chloride and formalin) and two extraction kits (PowerWater and NucleoSpin) on amplicon sequencing of 16S and 18S rRNA gene over 50 weeks of sample storage. Our results suggest the combination of mercuric chloride preservation and PowerWater extraction as most adequate for 16S and 18S rRNA gene amplicon-sequencing from the same seawater sample. This approach provides consistent information on species richness, diversity and community composition in comparison to control samples (nonfixed, filtered and frozen) when stored up to 50 weeks at in situ temperature. Preservation affects the recovery of certain taxa, with specific OTUs becoming overrepresented (SAR11 and diatoms) or underrepresented (Colwellia and pico-eukaryotes) after preservation. In case eukaryotic sequence information is the sole target, formalin preservation and NucleoSpin extraction performed best. Our study contributes to the design of long-term autonomous microbial observations in remote ocean areas, allowing cross-comparison of microbiome dynamics across sampling devices (e.g., water and particle samplers) and marine realms.

Introduction

Microbial communities have fundamental ecological and biogeochemical roles in nutrient recycling and carbon sequestration (Jørgensen and Boetius, 2007; Fuhrman et al., 2015). Understanding the consequences of global change for marine ecosystems requires a robust assessment of microbial community dynamics over temporal and spatial scales (Sunagawa et al., 2015; Buttigieg et al., 2018). Automated sampling devices attached to observational platforms, e.g., ocean moorings, enable time-series observations of microbial dynamics (Herfort et al., 2016; Zhang et al., 2019, 2021). Typically combined with physical and chemical sensors, automated samplers are of particular value in remote and inaccessible areas, such as seasonally ice-covered environments (Liu et al., 2020; von Appen et al., 2021; Wietz et al., 2021; Ramondenc et al., 2022). For instance, autonomous sediment traps allow linking particle flux with microbial diversity over extended periods, based on microscopic counts (Salter et al., 2007, 2012, 2014; Nöthig et al., 2020; Zúñiga et al., 2021) and DNA sequencing (Metfies et al., 2017; Bachy et al., 2022; Valencia et al., 2022).

There is a growing toolbox and increasing application of automated water and particle sampling approaches (Supplementary Table S1). As in situ molecular analysis is still an emerging technology (Moore et al., 2021) and beyond the resource capacity of many observing programs, automated samplers mostly perform in situ preservation of sample material (Yamahara et al., 2019; Lindsay, 2021; Truelove et al., 2022). In situ preservation intends to minimize signal modification over the extended duration of device deployment and laboratory processing. Formalin and mercuric chloride are commonly used to preserve sinking particles in long-term monitoring programs (Lee et al., 1992; UNESCO-IOC, 1994; Bauerfeind et al., 2009; Lampitt et al., 2010; Fischer et al., 2016). Although these chemicals originally aimed to preserve tissues, particles and cells for bulk biogeochemical analyses, recent studies have demonstrated that downstream molecular analyses are feasible with both mercuric chloride- (Metfies et al., 2017) and formalin-fixed (Boeuf et al., 2019) sediment trap samples. Likewise, preservation with mercuric chloride (Liu et al., 2020; Wietz et al., 2021) and formalin (Stern et al., 2015) allows ribosomal metabarcoding of microbes in autonomously collected seawater. Also the nucleic acid stabilizers RNAlater and DNAgard can preserve environmental DNA (Gray et al., 2013; Rachel and Gieg, 2020), however requiring frozen storage in stabilizer solution or the concentration of microbial biomass on filters (Ottesen et al., 2011). Both reagents have been tested as preservative in automated microbial samplings (Boeuf et al., 2019; Formel et al., 2021; Poff et al., 2021), but can lead to DNA loss (Renshaw et al., 2015) and are likely unsuitable in remote regions where samples cannot be frozen immediately. Hence, although automated technologies – in particular comparative sampling across different regions – offer exciting perspectives, preservation method and storage time are challenging factors for microbial diversity studies (Sherr and Sherr, 1993; Rissanen et al., 2010; Metfies et al., 2017; Spens et al., 2017; Sano et al., 2020; Pratte and Kellogg, 2021).

In the present study, we examined how preservation and DNA extraction methods affect molecular microbial analyses after long-term storage of seawater samples. Specifically, we addressed DNA yields, PCR amplification efficiency and microbiome composition after sample storage for 10, 28, and 50 weeks (0°C) to mimic long-term autosampler deployments. The approach was chosen to match deployment conditions of autonomous samplers in polar waters, which are installed on moorings and typically serviced only once per year (e.g., von Appen et al., 2021). We evaluate the consistency of 16S and 18S rRNA sequence information obtained from samples after different periods of post-sampling storage. We focus on formalin and mercuric chloride as they are widely used preservatives (Supplementary Table S1) and functionally different, particularly with respect to long-term storage at in situ temperatures. Furthermore, we aimed to assess how results from freshly preserved samples align with those from legacy samples, and indeed allow decadal-scale characterization of ecosystem dynamics. Our results have implications for microbial time-series collected with automated samplers, both regarding short-term methodological aspects and long-term archiving of biodiversity information.

Materials and methods

Experimental design and sampling regime

Approx. 6 l of surface seawater were collected at the pier on Helgoland Island in the German Bight (54° 10′ 58.3″N, 7° 53′ 19.9″E) on March 30, 2017. The water sample was kept at 4°C in the dark for ~ 35 days, then well mixed and split into 40 ml subsamples. Five subsamples were directly filtered as reference. The following preservatives were added to four sets of five replicate subsamples: (i) saturated mercuric chloride (HgCl2) solution (0.15% w/v final concentration per sample), (ii) 20% formalin (1.8% v/v final concentration per sample), (iii) RNAlater (1% final concentration per sample), and (iv) DNAgard (1% final concentration per sample). Preserved 40 ml subsamples were stored in the dark at 0°C to mimic conditions during high-latitude mooring deployments. After 10, 28, and 50 weeks, respectively (hereafter referred to as 10w, 28w, 50w), five replicates per preservation method were subjected to DNA extraction with two different kits after filtering each 20 ml onto Isopore membrane filters (Millipore, Burlington, MA, United States; 0.2 μm pore size, 47 mm diameter). Filters were stored frozen at −20°C for the same amount of time until DNA extraction with the NucleoSpin II (NS; Macherey-Nagel, Germany) or PowerWater (PW; QIAGEN, Germany) kit following the manufacturers’ protocols. Filters from formalin-preserved samples were subjected to additional rinsing steps before DNA extraction following Bucklin and Allen (2004). DNA extracts were quantified using a Nanodrop 1000 photometer (Thermo Fisher Scientific, Germany) and stored frozen until library preparation.

Amplicon sequencing

Libraries were prepared according to the standard instructions of the 16S Metagenomic Sequencing Library Preparation protocol (Illumina, San Diego, CA, United States). The V4 region of eukaryotic 18S rRNA genes was amplified using PCR primers 528F (5′-GCGGTAATTCCAGCTCCAA-3′; Elwood et al., 1985) and 964iR (5′-ACTTTCGTTCTTGATYRR-3′; Balzano et al., 2015). The V4-5 region of bacterial and archaeal 16S rRNA genes was amplified using primers 515F (5′-GTGYCAGCMGCCGCGGTAA-3′) and 926R (5′-CCGYCAATTYMTTTRAGTTT-3′; Parada et al., 2016). All PCRs had a final volume of 25 μl and contained 12.5 μl KAPA HiFi HotStart ReadyMix (Roche, Basel, Switzerland), 2.5 μl of each primer (1 μM) and 2.5 μl template. Amplification included initial denaturation (95°C, 3 min) followed by 25 cycles of denaturation (95°C, 30 s), annealing (55°C, 30 s), and extension (72°C, 30 s) with a single final extension (72°C, 5 min). 18S rRNA PCR products were gel-purified using the AMPure XP PCR purification kit (Beckman Coulter, Pasadena, CA, United States) according to the manufacturer’s protocol. All PCR products were quantified using a Quantus Fluorometer (Promega, Madison, WI, United States). Indices and sequencing adapters were attached via PCRs (final volume 50 μl), each containing 25 μl of KAPA HiFi HotStart ReadyMix (Roche), 5 μl of each Nextera XT Index Primer [1 μmol/l], 5 μl template (~5 ng DNA in total) and 10 μl PCR grade water. Amplification included initial denaturation (95°C, 3 min) followed by 8 cycles of denaturation (95°C, 30 s), annealing (55°C, 30 s), and extension (72°C, 30 s) with a single final extension (72°C, 5 min). 18S rRNA libraries were gel-purified using the AMPure XP PCR purification kit (Beckman Coulter). All libraries were quantified using a Quantus fluorometer (Promega) and sequenced using MiSeq and the MiSeq Reagent Kit V3 (2 × 300 bp) according to the manufacturer’s protocol (Illumina).

Processing and analysis of amplicon reads

Reads were processed using Trimmomatic v0.38 (Bolger et al., 2014) by scanning each sequence from the 5′ to 3′ end, trimming the 3′ end if average Phred Q-score of < 8 in a sliding window of 3 bp. Paired ends were merged using VSEARCH v2.3.0 (Rognes et al., 2016), discarding pairs with <50 bp overlap and > 5 mismatches in the overlapping segment. To guarantee identical orientation, sequences were filtered so forward sequences occur before reverse complement sequences. If sequences did not match this pattern, their reverse complement was also scanned using cutadapt v1.17 (Martin, 2011), requiring minimum overlaps of 17 and 13 bp for forward and reverse primer sequences respectively, and only one mismatch. Primer sequences were truncated, and sequences feature-filtered using VSEARCH. Sequences were discarded if (i) < 300 bp or > 550 bp, (ii) containing ambiguous bases (assigned as RYSWKMBDHVN per IUPAC nomenclature), or (iii) having an expected error (sum of all base error probabilities) > 0.25. Each sample was independently dereplicated, and the abundance of each sequence added to the sequence header. Chimeras were sample-wise predicted de novo by VSEARCH with default settings and removed. Subsequently, only samples with at least 10,000 sequences were used. Cleaned sample files were pooled and dereplicated in total, keeping amplicon abundances in the sequence headers. The pooled file was used as input for OTU clustering with SWARM v2.2.2 (Mahé et al., 2014), using the most abundant amplicon of an OTU as representative for annotation. Sequences were annotated with the default classifier implemented in mothur v1.38.1 using the Protist Ribosomal database v4.11.1 (Guillou et al., 2013) and the Silva v132 database (Quast et al., 2013) for 18S and 16S rRNA amplicons respectively, with a confidence cut-off of 80. One representative sequence was used to annotate the full OTU cluster, discarding singletons as well as OTUs with < 0.005% relative abundance. Statistical evaluation was carried out with R v.4.1.1 in RStudio using packages phyloseq, ampvis2, iNEXT, vegan, ape, tidyverse and scico (McMurdie and Holmes, 2013; Oksanen et al., 2013; Hsieh et al., 2016; Andersen et al., 2018; Paradis and Schliep, 2019; Wickham et al., 2019; Crameri, 2021). As our 16S rRNA dataset contained almost no archaeal sequences, 16S results are only referred to as “bacteria”. Relative abundances were Hellinger-transformed (the square root of the relative abundance per OTU and sample), an ecologically relevant transformation to correct for the compositionality of amplicon sequence data (Legendre and Gallagher, 2001).

Preliminary sequence analyses showed that only HgCl2 and formalin performed well in our experimental design (Supplementary Figure S1). The nucleic acid stabilizers RNAlater and DNAgard were originally tested, since being used in some automated sampling approaches (Supplementary Table S1). However, as nucleic acid stabilizers are not designed for long-term sample storage without freezing, we omitted results from RNAlater and DNAgard from further analysis.

Data and code availability

The entire workflow from raw sequence processing to statistical evaluation is available at https://github.com/matthiaswietz/MicroPreserve. Sequence data have been deposited in the European Nucleotide Archive (ENA) under accession number PRJEB43307, using the data brokerage service of the German Federation for Biological Data (GFBio) in compliance with MIxS standards (Yilmaz et al., 2011).

Results and discussion

We evaluated microbial community composition in seawater samples following two different preservation methods, based on poisoning (HgCl2) and fixation by protein cross-links (formalin). The concentrations of HgCl2 and formalin, common preservatives to study water column biogeochemistry and microbiology, were at the higher end of the range typically used, aiming at the observation of the strongest preservative effect expected. HgCl2 (0.15% w/v) and formalin (1.8% v/v) concentrations correspond to those used in particle traps (Bauerfeind et al., 2009; Lampitt et al., 2010). HgCl2 concentrations in water sampler deployments can be tenfold lower (von Appen et al., 2021; Wietz et al., 2021) as biomass in seawater is commonly lower compared to particles.

DNA yields and PCR amplification

Preservation with HgCl2 resulted in a higher proportion of successful DNA extractions compared to formalin (Table 1; Supplementary Figure S2) and approx. tenfold higher yields, despite pre-treatment of formalin-preserved samples (Bucklin and Allen, 2004). DNA-protein cross-linking through formalin may explain lower success rates and DNA yields. For both formalin and HgCl2, approximately twofold higher DNA yields were observed with PW extraction, likely corresponding to the combined bead-beating and enzymatic lysis compared to only chemical lysis with NS extraction (Yuan et al., 2015). In general, preservation decreased DNA yields two to fourfold compared to non-preserved controls. The impact of preservatives on DNA yields was observed at the earliest experimental time-point (10w), with no significant further decreases over the experimental period (Supplementary Figure S2). Hence, the chemical effect of preservatives is the major determinant of DNA yields, without further impact of prolonged storage, at least for up to 50 weeks. Independent of extraction kit, PCR amplification failures were ~ 50% for formalin compared to < 10% for HgCl2 (Table 1). Although formalin-preserved samples allowed DNA extraction and amplification in several cases, our results hence support that formalin can impede downstream molecular analyses (Hoffman et al., 2015; Reid et al., 2017).

Table 1

DNA yields and successful PCRs after preservation in comparison to the unpreserved reference, when extracted with either PowerWater (PW) or NucleoSpin (NS).

	Extractions with detectable DNA yield/sample number	DNA yield [ng μl−1]	Successful PCRs (16S/18S rRNA)
Reference; PowerWater	5/5	0.85 ± 0.2	5/5
Reference; NucleoSpin	5/5	0.2 ± 0.1	5/5
Mercuric chloride; PowerWater	15/15	0.24 ± 0.17	15/14
Mercuric chloride; NucleoSpin	13/15	0.03 ± 0.03	12/15
Formalin, PowerWater	5/15	0.2 ± 0.08	5/2
Formalin; NucleoSpin	11/15	0.01 ± 0.02	9/14

In some cases, PCR was successful despite NanoDrop did not detect DNA, probably related to detection sensitivity of the instrument. Reference: directly filtered environmental sample without preservation, immediately frozen at –20°C.

Microbial community composition

We obtained a mean of 34,000 and 62,000 chimera-filtered 16S and 18S rRNA amplicon reads, respectively (Supplementary Table S2). Principal coordinates analysis revealed clear clustering of both eukaryotic and bacterial communities by preservation method (Figure 1; PERMANOVA, p < 0.01), with little effect of storage time or extraction kit. Hence, in line with DNA extraction and PCR results, preservation method is the major determinant of ribosomal metabarcoding results under the specific microbial community and storage conditions tested in this study.

Figure 1

Principal Coordinates Analysis of Bray–Curtis dissimilarities of bacterial (left) and eukaryotic (right) communities after preservation and storage over different intervals. DNA extraction kits are not visually separated, as the influence of fixation significantly outweighs that of extraction.

Bacterial communities

Preservation significantly influenced bacterial community composition compared to the unpreserved reference (PERMANOVA, p < 0.001). However, differences to the unpreserved control were minor, with a taxonomic distance of ~ 0.2 particularly for HgCl2 samples and little change over time. In addition to preservation method, an effect of the extraction kit was observed. While HgCl2 + PW, HgCl2 + NS and formalin + NS performed comparably for bacterial communities, communities obtained from formalin + PW clustered separately (Supplementary Figure S3). The inverse Simpson index, considering both evenness and richness to determine alpha-diversity, was elevated after HgCl2 preservation (Figure 2; Kruskal–Wallis with Dunn’s post hoc test, p = 0.04). This concurred with higher relative abundances of planctomycetes, Deltaproteobacteria, and Actinobacteria (Figure 3), indicating that preservation can overestimate the rare biosphere. Among the major classes, preservation influenced the representation of alphaproteobacterial and gammaproteobacterial abundances (Figure 3), mainly relating to SAR11 clade Ia (higher) and Colwellia (lower abundances) respectively (Figure 4A). These taxa are at the lower and higher size spectrum of pelagic marine bacteria, respectively (Bowman, 2014; Giovannoni, 2017), indicating that preservation might favor smaller bacterial cells. Alternatively, cell wall structure and glycosylation (Dadon-Pilosof et al., 2017) might influence preservation efficiency. Compositionality effects can amplify such observations, but can be alleviated by normalizing relative abundances (Legendre and Gallagher, 2001; Weiss et al., 2017). Indeed, Hellinger-transformed relative abundances provided a more even picture of community structure (Figure 4B), with smaller differences for Colwellia while identifying highest variability for Amphritea (Gammaproteobacteria: Oceanospirillales). Previous studies have identified seasonal microbial dynamics in polar waters based on HgCl2 + PW preserved, autonomously collected samples (Liu et al., 2020; Wietz et al., 2021). Our results indicate that detection of Colwellia in such samples (Wietz et al., 2021) represented a true ecological finding, supported by stable OTU numbers from Alpha- and Gammaproteobacteria in HgCl2 + PW samples (Supplementary Figure S4).

Figure 2

Bacterial species richness and inverse Simpson index by preservation, storage time, and DNA extraction. The number of samples per group is shown in Supplementary Table S2.

Figure 3

Relative abundances of major bacterial classes (average of all replicates per sampling event) by preservation, storage time, and DNA extraction. PW: PowerWater, NS: NucleoSpin.

Figure 4

Relative abundances (A) and Hellinger-transformed relative abundances after HgCl2 preservation (B) of major bacterial genera (average of all replicates per sampling event) in relation to preservation, storage time, and DNA extraction. PW: PowerWater, NS: NucleoSpin, uc: unclassified.

Eukaryotic communities

For eukaryotes, PW extraction of formalin-preserved samples largely failed (Table 1). Hence, we restricted assessment of method performance to formalin + NS and HgCl2 samples. Preservation significantly influenced eukaryotic community composition compared to the unpreserved reference (PERMANOVA, p < 0.001), albeit with minor differences to the unpreserved control (maximum taxonomic dissimilarities of ~0.3) comparable to bacteria. Formalin + NS, HgCl2 + NS and HgCl2 + PW performed similarly, providing comparable composition and diversity patterns compared to the reference (Figures 1, 5; Supplementary Figure S5). HgCl2 and formalin resulted in higher proportions of Bacillariophyta (i.e., diatoms) in comparison to the unpreserved reference. In addition, Filosa-Imbricatea and unclassified stramenopiles were overrepresented in HgCl2 + NS (Figure 6). The total number of OTUs detected within stramenopile groups was lower after both HgCl2 and formalin preservation (Supplementary Figure S4), contributing to an overall lower species richness compared to the reference (Figure 5). The relative abundances of picoplankton classes Picozoa, MAST and Choanoflagellata were most similar between the reference and formalin + NS (Figure 6). As opposed to the overrepresentation of smaller bacterial cells, HgCl2 preservation favored larger-size eukaryotes such as centric diatoms, with higher abundances of especially unclassified Mediophyceae compared to the reference (Figure 7A). If resources allow, additional microscopy (Metfies et al., 2017), quantitative PCR or flow cytometry are advised to assess the effect of preservation on cell numbers and/or size classes. As for bacteria, Hellinger-transformed data provided a more even picture of community structure (Figure 7B).

Figure 5

Eukaryotic species richness and inverse Simpson index by preservation, storage time, and DNA extraction. The number of samples per group is shown in Supplementary Table S2.

Figure 6

Relative abundances of major eukaryotic classes (average of all replicates per sampling event) by preservation, storage time, and DNA extraction. PW: PowerWater, NS: NucleoSpin, uc: unclassified.

Figure 7

Relative abundances (A) and Hellinger-transformed relative abundances after HgCl2 preservation (B) of major eukaryotic genera (average of all replicates per sampling event) by preservation, storage time, and DNA extraction. PW: PowerWater, NS: NucleoSpin, uc: unclassified.

A comparison of results from the five technical replicates per treatment and time point allowed assessing the variability introduced by extraction and sequencing. Community structure in technical replicates were was highly reproducible for bacteria, but varied more for eukaryotes (Supplementary Figure S6). This observation potentially corresponds to disproportional distribution of large eukaryotic cells in some replicates, suggesting sample volumes should be maximized whenever possible. However, volumes and replicate numbers often need to be balanced with the desired temporal resolution, which can be challenging in remote locations relying on autonomous sampling.

Conclusion

Understanding the ecological and biogeochemical roles of marine microbes substantially benefits from automated sampling in long-term ecological observatories. We herein assessed the combined effects of preservation, DNA extraction and storage time on ribosomal metabarcoding of bacterial and eukaryotic communities. These insights inform the design of automated microbial observation in remote waters, which rely on in situ preservation and ex situ extraction after extended storage between sample collection and retrieval of the sampler. We present four major conclusions:

HgCl2 + PW provided the best representation of bacterial diversity and composition, even after 1 year of storage. Despite altering some patterns observed in the original community, abundances of the major taxa were overall reproducible and differences restricted to only few taxa.

Formalin + NS performed best for eukaryotes, despite low DNA yields. Although logistically demanding in (automated) field studies, sampling volumes should be as large as possible to maximize the robustness of analyses. Hellinger or centered-log ratio transformations can counteract the inherent compositionality of amplicon data and provide a more reasonable picture of microbial dynamics.

For parallel assessment of bacteria and eukaryotes, we recommend HgCl2 + PW, as this provides good 16S and reasonable 18S rRNA sequence information from single DNA extracts. Our results indicate that the HgCl2 shortcomings in eukaryotes outweigh the formalin shortcomings in bacteria, indicating HgCl2 as most suitable for observatories aiming to study both groups based on DNA from the same samples. Nonetheless, individual time-series should perform similar benchmark studies, as the respective strengths and weaknesses might differ at other in situ temperatures and for other microbial communities.

In order to minimize bias, we recommend that the choice of preservation should also consider potentially desired comparisons with other sites, as well as other samples from the same observatory. For instance, in case of the FRAM observatory of the Alfred Wegener Institute, the use of HgCl2 + PW facilitates cross-comparability with metabarcoding of sinking particles from sediment traps, including decade-old legacy samples that are treated similarly.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ebi.ac.uk/ena, PRJEB43307.

Author contributions

MW and KM analyzed the data and wrote the paper. CW performed the experimental work. CB, FJ, IS, and AB co-designed the study and contributed to writing. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by institutional funds of the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, funds provided within the framework of the Helmholtz infrastructure program Frontiers in Arctic Marine Monitoring (FRAM), the Horizon 2020 project AtlantOS (grant ID 633211), and the ERC AdvGrant project ABYSS (grant ID 294757) to AB.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.