Benchmarking protocols for the metagenomic analysis of stream biofilm viromes

Sampling

We sampled benthic biofilms from three streams (Switzerland) draining catchments differing in altitude and land use (Table 1). The Vallon de Nant (VDN) catchment is pristine with vegetation dominated by alpine forests and meadows. The Veveyse (VEV) catchment is characterized by mixed deciduous forests, but also agricultural and urban land use. The Senoge (SNG) catchment is clearly impacted by agricultural land use. During winter, benthic biofilms were randomly collected from stones (5–15 cm in diameter) using sterile brushes and Milli-Q water. Depending on biofilm thickness, we scraped biofilms from cobbles with a total surface area ranging from 1.3 to 2.4 m² into 10 L Milli-Q water. Slurries were transported on ice to the laboratory pending further processing.

Biofilm properties

For bacterial cell counting, samples were fixed with 3.7% formaldehyde (final concentration) and stored at 4 °C. Bacterial cells were disintegrated from the biofilm matrix using 0.25 mM tetrasodium pyrophosphate in combination with rigorous shaking (1 h) and sonication (Bransonic Sonifier 450, Branson, MO, USA) on ice (1 min) (Velji & Albright, 1993). As previously established (Besemer et al., 2009), cells were stained using Syto13 and counted on a flow cytometer (NovoCyte, ACEA Biosciences, San Diego, CA, USA). Chlorophyll a content was determined spectrophotometrically after over-night extraction in 90% EtOH (Lorenzen, 1967). Extracellular polymeric substances (EPS) were extracted from biofilm slurries using 50 mM EDTA and shaking (1 h) (Battin et al., 2003). Carbohydrates were precipitated in 70% EtOH (−20 °C, 48 h) and measured spectrophotometrically as glucose equivalents (DuBois et al., 1956). Proteins were determined according to Lowry et al. (1951).

Transmission electron microscopy

A first confirmation of the presence of VLP in biofilms was obtained from TEM. For this, 5 µL of unprocessed sample was adsorbed to a glow-discharged carbon-coated copper grid (Canemco & Marivac, Gore, QC, Canada), washed with deionized water and stained with 5 µL of 2% uranyl acetate. TEM observations were made on an F20 electron microscope (ThermoFisher, Hillsboro, OR, USA) operated at 200 kV and equipped with a 4,098 × 4,098 pixel camera (CETA, ThermoFisher, Waltham, MA, USA). Magnification ranged from 10,000 to 29,000×, using a defocus range of −1.5 to −2.5 µm.

Protocols for the extraction of VLP

In order to establish an optimized stream biofilm sample-to-sequence pipeline, we explored a variety of protocols for the concentration, extraction and purification of VLP (Fig. 1). The pipeline consists of three main parts: concentration, extraction and purification. The concentration step is required to obtain sufficient genetic material for nucleic acid extraction. Extraction aims to liberate VLP from the biofilm matrix, while purification aims to reduce the amount of contaminant nucleic acids and cellular debris. The order—first volume reduction, then chemical detachment—was chosen because chemicals used for extraction of VLP may damage tangential flow filters. We tested all possible combinations of protocols to identify the most promising laboratory pipeline for the preparation of viral metagenomes.

First, we homogenized the biofilm slurries by manual shaking and split samples into aliquots (1 L). Aliquots were centrifuged at 100×g (15 min) (5810R; Eppendorf, Hamburg, Germany) to remove large sediment particles and organisms contained within the biofilm slurry. We recovered the supernatant for downstream analyses.

Enumerating VLPs using epifluorescence microscopy

To assess the relative viral extraction efficiency of the various protocols, we counted VLP under an epifluorescence microscopy (AxioImager Z2, Zeiss, Germany) as described by Patel et al. (2007). Because of high background noise owing to the biofilm matrix constituents, samples could only be counted after the purification steps. It should be noted that pyrophosphate concentrations >10 mM tend to interfere with VLP counting (Danovaro et al., 2001); however, we used five mM pyrophosphate during VLP liberation from biofilms. Subsamples were fixed with formaldehyde, stained with 0.5 µL of 1,000 × SYBR Gold (Molecular Probes; ThermoScientific, Waltham, MA, USA) and incubated at room temperature in the dark (30 min). After incubation, each sample was filtered onto a 0.02 μm pore size membrane filter (Anodisc; Whatman, Maidstone, UK). The filters were mounted on glass slides with a drop of VECTASHIELD antifade mounting medium (Vector Laboratories, Burlingame, CA, USA). VLP were visualized using blue light (488 nm) excitation and green (512 nm) emission. For each sample, 15–20 randomly selected images were acquired with a camera (Axiocam 506 mono, Zeiss, Germany) mounted onto the microscope. VLPs were discriminated from bacteria by size (0.015–0.2 µm) and enumerated using a custom script in Fiji (Schindelin et al., 2012).

Nucleic acid extraction

To further assess the efficiency of the tested protocols to obtain viral nucleic acids for metagenome sequencing, we extracted DNA and quantified its concentration. Following Sambrook, Fritsch & Maniatis (1989), 0.1 volume of TE buffer, 0.01 volume of 0.5 M EDTA (pH 8) and 1 volume of formamide was added to purified viral samples, and incubated at room temperature (30 min). Then, two volumes of cold 100% EtOH were added and incubated for 30 min (4 °C). Samples were centrifuged at 17,000×g (4 °C, 20 min) and pellets washed twice with 70% cold EtOH. Pellets were air-dried and resuspended in 567 µL of TE buffer (10 mM Tris and one mM EDTA (pH 8.0)). Thirty µL of warm 10% SDS and three µL of proteinase K (20 mg/mL) were added and samples were incubated for 1 h at 37 °C. Next, 100 µL of five M NaCl and 80 µL of warm hexadecyltrimethylammonium bromide (CTAB) were added and incubated for 10 min (65 °C) (Doyle & Doyle, 1987). DNA was extracted in a series of chloroform, phenol:chloroform:isoamyl alcohol (25:24:1) and chloroform treatments, with centrifugation at 16,000×g (10 min) for phase separation. Finally, DNA was precipitated overnight in isopropanol (−20 °C). DNA was concentrated by centrifugation at 16,000×g (4°C) for 20 min, and the pellet washed twice with cold 70% EtOH, air-dried, resuspended in 50 μL nuclease-free H₂O and stored (−20 °C). DNA concentration was measured using Qubit and the dsDNA high-sensitivity kit according to the manufacturer’s instructions (Life Technologies, Carlsbad, CA, USA).

Library construction and sequencing

Based on the epifluorescence microscopy counts of VLP and the DNA concentration from purified biofilm samples, we selected samples processed with the best performing pipeline for metagenome sequencing. For sequencing library construction, DNA was sheared with an S2 focused ultrasonicator (Covaris, Woburn, MA, USA) to achieve a target size of DNA fragments of around 350 bp. We opted for the ACCEL-NGS® 1S PLUS DNA library kit (Swift Biosciences, Ann Arbor, MI, USA) which allows low quantities of both single- and double-stranded DNA as input (Roux et al., 2016). Library construction and multiplexing was performed following the manufacturer’s instructions for DNA inputs (<1 ng/μL) and 20 cycles of indexing PCR. Paired-end sequencing (2 × 300 bp) was performed on a MiSeq System (Illumina, San Diego, CA, USA) at the Lausanne Genomic Technologies Facilities. Raw sequences have been submitted to the European Nucleotide Archive under accession number PRJEB33548.

Bioinformatic analyses

The number of reads and GC content of each sample before and after quality control were calculated using a custom python script. For virome classification, we followed the recommendations to assemble viral contiguous sequences (contigs) according to Roux et al. (2019). First, BBDuk (v35.79) was used to remove Illumina adapters, filtering and trimming (trimq = 12). Next, reads with >93% similarity to a human reference genome were discarded (using BBmap). The ACCEL-NGS^® 1S Plus library preparation kit includes a low complexity adaptase tail, which was clipped (10 bases) according to the manufacturer’s instructions. We then used the error correction capability of Tadpole (v. 37.76) to correct for sequencing errors (mode = correct ecc = t prefilter = 2). Prior to assembly, we de-duplicated our datasets using clumpify (v37.76) with parameters set such that identical reads were identified and only one copy was retained (dedupe subs = 0). Finally, we used the SPAdes assembler (v3.13.0, Bankevich et al., 2012) in single-cell (—sc) mode, with error correction disabled (—only-assembler) and kmers set to 21, 33, 55, 77, 99, 127 to assemble contigs. We co-assembled the paired-end reads from both sequencing runs for each sample individually. To obtain an overview of the potential contaminant sequences (e.g., human, bacterial and phiX), we uploaded the quality-trimmed and de-replicated reads (forward orientation only) to the web interface of taxonomer (Flygare et al., 2016). Taxonomer is an ultrafast taxonomy assigner, which assigns and classifies reads to human, bacterial, viral, phage, fungal, phix, ambiguous (i.e., reads fit to more than one bin) and “unknown” bins. For identification and classification of viral contigs, we used MetaPhinder (v2.1, Jurtz et al., 2016) through the web interface hosted at the Center for Genomic Epidemiology at the Danish Technical University (DTU). We mapped the reads to the contigs using BWA-MEM (http://bio-bwa.sourceforge.net/) using default setting (e.g., a mismatch penalty of four and a minimal alignment score of 30) and then counted the number of reads mapping each contig using samtools (Li et al., 2009). The average number of reads per contig were adjusted by contig length and normalized by library size to obtain the average percent of reads mapping to contigs. In order to specifically target ssDNA contigs (Trubl et al., 2019), we matched the contigs to the Viral_rep and Phage_F domains of the PFAM database using hmmsearch (HMMER v3, Eddy, 2011) with cutoff scores ≥50 and e-values ≤ 0.001.

Biofilm properties

Bacterial abundance ranged between 4.1 × 10¹¹ cells m⁻² at the lowest stream (SNG) and 2.3 × 10⁹ cells m⁻² at the uppermost stream (VDN; Table 1). This pattern was mirrored by an increase in chlorophyl a content and the protein and carbohydrate concentrations in EPS that were higher in the lowland than the high-altitude stream (Table 1). Despite these differences in biofilm properties, no differences in extraction efficiency among the different protocols were detected. In fact, the average number of VLPs extracted by the various protocols were strongly correlated among the different samples (pairwise Spearman’s r_s ranging between 0.94 and 0.98).

Transmission electron microscopy

Direct electron microscopic observations of raw biofilm samples revealed the presence and morphological diversity of virions, including tailed bacteriophages and lemon-shaped viruses, in biofilms from all three streams (Fig. 2). Polyhedral, spherical and filamentous VLPs were also observed, which may include untailed bacteriophages or viruses infecting eukaryotes. Some amorphous structures may also be interpreted as membrane vesicles.

Virome extraction and purification

In total, we evaluated 16 protocols to concentrate, extract and purify viruses from benthic biofilms from three streams (Fig. 1). Based on the number of VLPs and DNA yields retained at the end of each protocol, we observed significant differences in relative extraction efficiencies among protocols. Across all samples, average VLP counts and DNA yields were correlated (Spearman’s r_s = 0.74), suggesting conformity among these two means of evaluation. The combination of TFF for concentration, tetrasodium pyrophosphate and sonication for extraction, and sucrose gradient centrifugation for purification resulted in the highest VLP counts in all three samples (two-way ANOVA, p < 0.01; Fig. 3). This pipeline also yielded the highest DNA concentration (Fig. 4). Protocols involving PEG precipitation generally resulted in a lower recovery of VLPs (on average only 13.6% compared to the best performing pipeline) and DNA yields below detection limit. This may be attributable to the formation of a visible, viscous layer upon addition of PEG to the biofilm samples. The dissociation of VLPs from the biofilm matrix was most effective using tetrasodium pyrophosphate and sonication (two-way ANOVA, p < 0.01). Protocols based on TFF and using DTT or chloroform extracted on average 25.0% and 33.2% less VLPs than protocols using TFF followed by tetrasodium pyrophosphate treatment and sonication (Fig. 3). Samples without any physicochemical treatment to extract VLP from biofilms yielded on average 54.8% less VLPs than the tetrasodium pyrophosphate and sonication treatment. To further purify viruses, two discontinuous gradient formulations with sucrose and CsCl were tested. On average, 1.9 times more VLPs were retained by ultracentrifugation using the sucrose gradient than using the CsCl gradient in samples concentrated using TFF and treated with tetrasodium pyrophosphate and sonication (paired-T-Test, p< 0.01). This is also reflected in DNA yields, which reached 1.22 ng µL⁻¹ in VDN, 1.31 ng µL⁻¹ in VEV and 18.7 ng µL⁻¹ in SNG using TFF, pyrophosphate and sonication followed by sucrose gradient centrifugation. DNA yields were on average 3 times higher using ultracentrifugation in the sucrose compared to the CsCl gradient (paired T-Test, p= 0.03) and on average 1.5, 1.7 and 1.9 times higher using tetrasodium pyrophosphate and sonication compared to DTT, no physico-chemical detachment and chloroform, respectively.

Given that stream biofilms contain abundant prokaryotic, eukaryotic and extracellular DNA, it is crucial to verify the absence of DNA potentially contaminating the samples. Negative PCR results from samples treated with DNase I confirmed the absence of DNA contamination from prokaryotic cells.

Stream biofilm viromes

From the two sequencing runs we obtained 24388096, 22459218 and 23804190 paired-end reads from SNG, VEV and VDN, respectively (Table 2). After quality control, error correction and deduplication, on average 97.5% of the reads remained. Initial screening using taxonomer showed that human contaminant reads accounted for 0.13–0.25% of the reads, while bacterial contaminant reads accounted for 1.47–8.11% of the reads.

We obtained 3,698, 11,323 and 13,591 contigs from de-novo assembly of quality-controlled and deduplicated reads from SNG, VEV and VDN, respectively. The largest contigs were 9,493, 46,665 and 54,492 bp in SNG, VEV and VDN, respectively. Hmmsearch against the PFAM databases did not yield ssDNA contigs in the three viromes. Between 726 and 2613 contigs were classified as of viral origin in the three viromes (Fig. 5). In all samples, the majority of contigs were identified as not further classified Siphoviridae (28.3–53.3%), followed by Myoviridae (9.6–18.0%) and Podoviridae (4.7–5.2%). Contigs classified as T4virus, Cp220virus, Kayvirus and P12024virus were common in all three biofilm viromes. However, contigs classified as Twortvirus, Phicbkvirus, Coopervirus and L5virus were only detected in viromes obtained from VEV and VDN but not in SNG. Across the three samples, on average 554.6 reads mapped to contigs classified as Caudovirales, averaging 63.4, 63.8 and 89.1% bp accounting for contig size in SNG, VEV and VDN, respectively. On average, 153.8 reads mapped to Podoviridae (range 2.2–30.8% of bp) and 115.1 reads mapped to Myoviridae (range 0.4–6.1% of bp), whereas only 29.4 reads mapped to Siphoviridae (range 1.8–17.7% of bp). There was substantial variation in the number of reads mapped to contigs in the different samples. For instance, more reads mapped on average to unclassified bacterial viruses and Myoviridae in SNG (311.1 and 339.7 reads accounting for 9.3 and 6.1% of bp) than in VDN (2.6 and 2.8 reads accounting for 0.6 and 0.8% bp) or VEV (36.7 and 2.8 reads accounting for 1.5 and 0.4% bp).