Optimised biomolecular extraction for metagenomic analysis of microbial biofilms from high-mountain streams

Sample origin & collection

DNA extraction protocols were benchmarked using three different GFS sediments from the Swiss Alps: Corbassière (CBS, collection date: 13.11.2018), 2444 m above sea level (m a.s.l) and Val Ferret (FE), 1995 m a.s.l at the glacier snout (up site, FEU, collection date: 23.10.2018) and one kilometer further downstream (down site, FED, collection date: 24.01.2019). Sampling was always performed later in the morning, before noon. Sediments (0.25 to 3.15 mm) were collected using two flame-sterilized metal sieves with a mesh size of 3.15 mm and 0.250 mm respectively. CBS differs from FEU and FED in terms of bedrock geology, with clastic sedimentary limestone dominating the catchment of CBS and brecchiated gneiss dominating in FEU and FED. Sediments generally contain more organic material further downstream from the glacier, which may inhibit DNA extraction. Sediments (0.25 to 3.15 mm) were collected using flame-sterilized sampling equipment. Wet sediments were transferred into 10 ml sterile, DNA/DNase-free tubes and immediately flash-frozen in liquid nitrogen in the field. Samples were transferred to the laboratory and kept at −80 °C until analysis. All necessary measures were taken to ensure contamination-free sampling.

DNA extraction methods

Four different DNA extraction methods were applied to the samples. The key characteristics of the different methods are summarized in Table 1. Method-1, -2 and -4 were non-commercial protocols differing primarily in the lysis step (bead-beating and lysis buffer composition; Roume et al., 2013; Lever et al., 2015; Sambrook & Russell, 2006) while method-3 was a modified protocol of the DNeasy PowerMax Soil Kit (Cat.No. 12988-10) provided by Qiagen (based on communication exchanged with the manufacturer). Due to the very low microbial abundance, additional precautions were taken to establish contamination-free conditions, including daily decontamination of equipment/areas with bleach, using DNA/DNase-free glassware and plasticware, reagents and chemicals. Additionally, “germ-free” sediment is not a viable option and is hard to remove any and all microorganisms from sediments. So, extraction blanks, i.e., tubes without any sample, were used as controls, which underwent the same extraction protocols along with the other samples. Post-DNA recovery, we assessed whether any of the eluted samples had DNA via both NanoDrop and Qubit, and found them to be below detectable levels. Additionally, the PCR reactions during library preparation did not yield any product confirming serving as a contamination-check. The input weight of sediment ranged from 0.5–5 g and are described in the respective protocols.

Method-1 was based on a previously established method (Griffiths, Whiteley & Bailey, 2000). Introduced modifications concerned primarily the step of mechanical lysis and DNA precipitation that was rendered more stringent to improve the recovery of small amounts of DNA. Sample cell lysis was achieved by adding 0.5 g of sediment into a lysing matrix E tube with beads of variable diameter provided by the manufacturer (MP Biomedicals, SKU 116914050), 500 µl CTAB buffer (5% CTAB, 120 mM KPO4, pH 8.0) and 500 µl of phenol/chloroform/isoamyl alcohol (ratio 25:24:1). Samples were loaded on a Precellys beater for 45 s at 5.500 r/s. DNA was extracted once more with chloroform/isoamyl alcohol (24:1) and precipitated with 2 vol PEG-6000, 15 µg/ml linear polyacrylamide (LPA) and 2 h incubation on ice (Supplemental Information).

Method-2 was an adaptation to alpine stream sediments of the modular method for DNA extraction previously published (Lever et al., 2015). The appropriate modules of the method, based on the nature of our samples, were put together in our protocol without further modification. Samples were prepared by mixing 5 g of sediment, 10–20% of 0.1 mm zirconium beads and 1 ml of 100 mM dNTP solution. Cell lysis was achieved with 5 ml lysis buffer (30 mM Tris-HCl, 30 mM EDTA, 1% Triton X-100, 800 mM guanidium hydrochloride, pH 10.0) and incubation at 50 °C for 1 h with gentle agitation (Hybridization oven, Labnet Problot L6). The supernatant was extracted once with chloroform/isoamyl alcohol (24:1) and DNA was precipitated with 10_260∕280g/ml LPA, 0.2 vol 0.5 M NaCl, 2.5 vol ethanol and 2 h incubation at RT in the dark (Supplemental Information). The input weight of 5 g sediment was a modification from previously established protocol and the subsequent reagent volumes were adjusted accordingly.

Method-3 has been previously applied successfully on sand and clay soils (Hale & Crowley, 2015). In this protocol, the standard lysis capacity of the DNeasy PowerMax Soil Kit (Qiagen, Cat. No. 112988-10) was modified and enhanced by the addition of phenol:chloroform:isoamyl alcohol along with PowerBeads (kit provided) and C1 solution to 5 g of sediment. Then, the manufacturer-suggested sequence of treatments and rinses with the standard buffers of the kit were followed to reach elution of extracted DNA from silica columns with 6 ml of elution buffer. Further concentration of extracted DNA was carried out with the addition of 240 µl 5M NaCl, 2.5 vol ethanol and 10µg/ml LPA (Supplemental Information). LPA was an additional modification to the original protocol for improved DNA recovery.

Method-4 involved chemical and enzymatic treatment of samples according (Green & Sambrook, 2017) with minor modifications. Five g of sample was mixed with 10 ml of lysis buffer, incorporating the SDS as well, (0.1 M Tris-HCl pH 7.5, 0.05 M EDTA pH 8, 1.25% SDS) and 10 µl RNase A (100 mg/ml). Then sediment was vortexed for 15 s and incubated at 37 °C for 1 h in a hybridization oven. 100 µl Proteinase K (20 mg/ml) was added in a subsequent step and the mixture was incubated for 10 min at 70 °C. Samples were extracted once with phenol/chloroform/isoamyl alcohol (ratio 25:24:1) and supernatants were extracted subsequently with chloroform/isoamyl alcohol (24:1). More stringent DNA precipitation conditions were applied with the addition of 10 µg/ml LPA, and overnight incubation at −20 °C (Supplemental Information).

All DNA extracts were suspended in 100 µl of DNA/DNase–free water (ThermoFisher Scientific). Due to the inadequacy of DNA obtained from Method-1 given the 0.5g input sediment weight, we scaled the extraction to 5 g starting weight prior to sequencing. Extracted DNA was thereafter stored at −20 °C until further use. Due to the low DNA yields, it was necessary to use Qubit dsDNA HS kit (Invitrogen), a fluorescent DNA quantification method with high sensitivity. Quality assessment, with Nanodrop and DNA visualization on 0.8% agarose gel containing GelRed nucleic acid stain, was possible only for DNA extracted with method-4 and for DNA concentrations higher than 0.5 ng/ul. All samples yielded sufficient DNA, i.e., 50 ng (total input), for metagenomic sequencing and subsequent analyses. Additionally, a commercially-available microbial mock community (ZymoBIOMICS, Cat.No. D6300) was extracted using Method-4 and used for subsequent sequencing.

DNA sequencing

50 ng of DNA from all samples were subjected to random shotgun sequencing. The sequencing libraries were prepared using the NEBNext Ultra II FS DNA Library Prep kit for Illumina (Cat.No. E7805) using the protocol provided with the kit. The libraries were prepared considering 350 base pairs (bp) average insert size. Qubit (Invitrogen) was used to quantify the prepared libraries while their quality was assessed on a Bioanalyzer (Agilent). We used the NextSeq500 (Illumina) instrument to perform the sequencing using 2 ×150 bp read length at the Luxembourg Centre for Systems Biomedicine Sequencing Platform.

Genome reconstruction and metagenomic data processing

Paired sequences (i.e., forward and reverse) were processed using the Integrated Meta-omic Pipeline (IMP) (Narayanasamy et al., 2016). The metagenomic workflow encompasses pre-processing (read quality filtering and trimming), assembly, and genome reconstruction in a reproducible manner. The adapter sequences were trimmed in the pre-processing step including the removal of human reads. Thereafter, de novo assembly was performed using the MEGAHIT (version 2.0) assembler (Li et al., 2015). Default IMP parameters were retained for all samples. Subsequently, we used MetaBAT2 (Kang et al., 2019) and MaxBin2 (Wu, Simmons & Singer, 2016) for binning in addition to an in-house binning methodology previously described (Heintz-Buschart et al., 2017). The latter method initially ignores the ribosomal RNA sequences in kmer profiles based on VizBin embedding clusters (Laczny et al., 2015). In this context, VizBin utilises density-based non-hierarchical clustering algorithms and depth of coverage for genome reconstructions. Subsequently we obtained a non-redundant set of metagenome-assembled genomes (MAGs) using DASTool (Sieber et al., 2018) with a score threshold of 0.7 for downstream analyses. The abundance of MAGs in each sample was determined by mapping the reads to the reconstructed genomes using BWA-MEM (Li, 2013), taking the average coverage across all contigs. Diversity measures from metagenomic sequencing were assessed by determining the abundance-weighted average coverage of all the reads to identify the number of non-redundant read sets (Rodriguez-R & Konstantinidis, 2014).

Taxonomic classification for metagenomic operational taxonomic units

We used the trimmed and pre-processed reads from the IMP workflow to determine the microbial abundance and taxonomic profiles based on the mOTU (v2) tool (Milanese et al., 2019). Based on the updated marker genes in the mOTU2 database including those from the TARA Oceans Study (Sunagawa et al., 2015) and recently generated MAGs (Tully, Graham & Heidelberg, 2018), taxonomic profiling was performed on our sequence datasets. We used a minimum alignment length of 140 bp to determine the relative abundances of the mOTUs, including the normalisation of read counts to the gene length, also accounting for the base coverage of the genes. Additionally, we used CheckM (Parks et al., 2015) to assess completeness and contamination. Subsequently, taxonomy for MAGs recovered after the redundancy analyses from DASTool was determined using the GTDB (Genome Taxonomy Database) toolkit (gtdb-tk) (Parks et al., 2018).

Data analysis

All figures for the DNA concentrations, library preparation, assembly metrics and supplementary figures were generated using GraphPad Prism (v8.3.0). Taxonomical assessment and diversity measures were created using version 3.6 of the R statistical software package (R Core Team, 2013). DESeq2 (Love, Huber & Anders, 2014) with FDR-adjustments for multiple testing were used to assess significant differences in the MAG abundances. The genomic cluster figure for the mock community was obtained as an output from the IMP metagenomic pipeline.

Phenol-chloroform-based extraction method results in higher DNA yields

To ensure native sequencing, by minimizing the number of PCR (polymerase chain reaction) cycles within the library preparation protocols, we tested four protocols for biomolecular extraction, with an aim of acquiring large quantities (>50 ng) DNA from glacier-fed stream benthic sediments. The four methods tested were selected because of their wide applicability on related environmental samples (Method-1 & -2) (Griffiths, Whiteley & Bailey, 2000; Lever et al., 2015; Tatti et al., 2016) and their improved chances of higher yields (Method-3; Qiagen communication). Since method-4 is considered the gold standard of DNA extraction in biomedical sciences (Dairawan & Shett, 2020) and bacterial cultures (Green & Sambrook, 2017), it was included in our study. The four protocols are largely based on the same principles, viz. sample preparation, cell lysis, purification, precipitation and washing (Table 1). From preliminary tests, it became apparent that a small-scale approach (Method-1; 0.5 g input sediment) did not yield sufficient amounts of DNA for metagenomics due to, on average, limited microbial biomass in the samples. Thus, all protocols (aside from Qiagen’s - already produced for maxi scale) were scaled up to 5 g of input sediment and a co-precipitant, like linear polyacrylamide, was included in all precipitation steps. This was essential for the quantitative recovery of the small amounts of extracted DNA from high solution volumes (6–10 ml).

Overall, we found that extractions using the commercial kit from Qiagen (method-3) yielded increased total DNA as compared to a commonly used protocol (method-1; Fig. 1A). Furthermore, method-3 was similar in terms of DNA yield when compared to a generalized protocol (method-2) previously proposed (Lever et al., 2015) (Fig. 1B). On the other hand, the phenol-chloroform based extraction protocol (method-4) was tested against both methods 2 and 3, using sediment samples collected from the three different glacier floodplain streams (CBS, FEU, FED) from Switzerland. Method-1 was omitted from these tests due to insufficient DNA concentrations in the preliminary extractions. We found that for all three GFS, the phenol-chloroform extraction yielded the highest DNA concentrations. In some cases, and notably samples with low cell abundance, we even obtained one order of magnitude more DNA (Fig. 1C).

Quality assessment of these DNA extracts with Nanodrop showed OD260/280 ratios between ∼1.4 and ∼1.6. Agarose gel electrophoresis revealed a high-molecular weight band with no apparent shearing, smearing or residual RNA, indicative of high-quality DNA (Fig. 2A). A secondary effect appearing in certain samples, but without any perceived consequences in the quality of extracted DNA whatsoever, was the development of a pink-red color of varying intensities with the addition of phenol:chloroform:isoamyl alcohol (Fig. 2B). This was pH dependent since samples were decolorized with the addition of sodium acetate pH 5.2 in the precipitation step. This could possibly be due to a ferric-chloride-phenol compound formed when chloride and phenol constituents of the protocol interact reversibly with Fe+3 ions contained in certain samples depending on local geology (Banerjee & Haldar, 1950). Similar coloration has been previously reported (Lever et al., 2015).

Extraction method affects library preparation efficiency

The DNA extractions based on method-3 and using phenol-chloroform methods were subsequently subjected to library preparation for high-throughput whole genome shotgun sequencing. Despite similar quality of DNA across both methods (∼1.4–1.6 OD260/280), library preparation using the modified commercial kit did not yield any successful libraries (Fig. 3). To assess if any impurities or inhibitors hampered library preparation we tested two clean-up methods for the DNA extracted from the commercial kit: (1) ethanol precipitation and (2) magnetic-bead based clean-up. For magnetic bead clean-up the SPRIselect beads (Beckman Coulter, 23318) were used according to the manufacturer’s protocol. We found that the magnetic-bead method leads to a complete loss of sample (i.e., undetectable DNA quantity via Qubit analyses) during the process, especially if starting with a low input DNA concentration. Although we lost six out of twelve samples using the magnetic-bead clean-up, we achieved 100% efficiency as indicated by a concentration of greater than 0.5 ng/ul after library preparation quantified by Qubit, with the remaining six samples. On the other hand, ∼20% of the samples cleaned via ethanol precipitation failed library preparation. Contradictory to these methods, DNA extracted using the phenol-chloroform based method (method-4) yielded 100% efficiency in terms of library preparation without any additional clean-up (Fig. 3). Additionally, we found that the distribution of the total yield after library preparation using the phenol-chloroform method was more uniform across samples compared to the other methods (Fig. 3).

Whole genome shotgun assembly unaffected by extraction methods

Extraction methods for whole genome shotgun sequencing may affect the sequencing itself, including the quality and assembly of the reads downstream. To assess this, we used the libraries prepared as described above (Fig. 3), and performed whole genome shotgun sequencing on an Illumina NextSeq500. The average quality across all three methods based on short-read sequencing was Q30 after trimming the leading and trailing sequences (described in Methods). We assessed several assembly metrics including the average length of contigs (N50), largest alignment, total aligned length and coverage. We did not find any significant differences among any of these measures across all three methods (Figs. 4A–4C, 4E). Using a diversity index metric, we however found a more uniform distribution across all samples prepared using method-4, albeit no significant differences to the commercial kit-based extraction and library preparation (Fig. 4D).

Extraction methods influence metagenomic profiles

It is well established that extraction methods (Mackenzie, Brett & Taylor, 2015) and library preparation (Bower et al., 2015) protocols affect the taxonomic profiles and genomes recovered after high-throughput sequencing. We determined if the preparation methods affected the overall diversity of taxa recovered and found that phenol-chloroform and the magnetic-bead clean-up methods demonstrated similar levels of diversity (Shannon) as compared to samples precipitated using ethanol (Fig. 5A). Overall, the community profiles of the ethanol precipitation-based method were highly diverse (Fig. 5B). Interestingly, the genomes recovered and their abundances were similar in the phenol-chloroform and magnetic-bead methods as well (Fig. 5C). However, we observed a significant increase (p < 0.001, FDR-adjusted p-value) in the abundance of a Ralstonia genome when prepared with the ethanol precipitation protocol (Fig. S1). Additionally, we found that the number of genomes recovered using the phenol-chloroform was more consistent with previously reported 16S rRNA gene sequencing profiles for GFS from Austria (Wilhelm et al., 2013; Besemer et al., 2012; Wilhelm et al., 2014). Simultaneously, we used an approach to identify metagenomic operational taxonomic units (mOTUs) and found that the phenol-chloroform and magnetic-bead methods showed similar profiles of mOTUs compared to that of ethanol precipitation (Fig. 5D).

Efficiency of phenol-chloroform extraction on a mock community including eukaryotes

To determine whether the phenol-chloroform extraction method is biased against eukaryotes, we used a commercially-available mock community (ZymoBIOMICS Microbial Community Standard #D6300) to assess bias and errors. After sequencing, we recovered high quality (>90% completion, <5% contamination) bacterial genomes (Fig. 6A). Additionally, the abundance of the microbial genomes, including one of the eukaryotes - Saccharomyces cerevisiae, were similar to the expected levels in the mock community (Fig. 6B). On the other hand, the protocol enabled the identification and partial recovery of the Cryptococcus neoformans genome, albeit at lower levels possibly due to increased melanisation of the cell wall (Grossman & Casadevall, 2017) affecting lysis and subsequent extraction.