Benchmarking of DNA Extraction Methods for Metagenomics, 16S rRNA Databases and Assessment of the Skimmed Milk Flocculation for Bacterial Enrichment for Metagenomics from Water Samples

1. Introduction

Freshwater and marine ecosystems are vital for the human survival and evolution and are frequently being studied through targeted or shotgun metagenomics. Metagenomics studies of these environments involve extracting genomic DNA (gDNA) from samples and analyzing the diverse organisms present, including bacteria, archaea, eukaryotes and viruses and/or their underlying functional characteristics.

Regardless of methods, study designs and goals, samples of various water volumes are collected, transported to laboratories and optimally processed immediately. The initial step of most bacterial metagenomics studies on water samples involves extracting the biomass (also called enrichment) from the water content. Most often, water samples include small particles such as soil, rocks, plant residues which could be extracted as well due to microorganisms and/or extracellular DNA adhere to those particles or floating free as planktonic forms. In the case of larger sample debris, additional steps are implemented to remove them prior to the biomass extraction. Three main techniques have been developed and used to varying extents for biomass extraction, namely vacuum filtration (VF), ultracentrifugation (UC) and skimmed milk flocculation (SMF). VF is mainly used for characterization of bacterial and fungal communities and is by far considered the gold standard. UC is almost exclusively applied for viral and phage metagenomics, whereas SMF was originally developed to concentrate viruses from coastal waters [1].

Filtration is a straightforward size exclusion or inclusion method used in metagenomics to separate microbial cells and sample debris from the water by vacuum as pressure-driven factor. However, there are numerous factors influencing the filterability during vacuum filtration such as the membrane protein binding affinity, surface charge, hydrophobicity, pore size and structure, and roughness. Additionally, cell size, shape, flexibility, charge and hydrophobicity also influence the potential of bacteria to flow though the filter, and lastly, the particulate nature of the sample suspension itself [2]. For those reasons, research in filtration methods is still evolving. On the other hand, SMF is used mainly for concentrating viruses such as SARS-CoV-2 [3], and others [1,2,3,4,5,6]. SMF is straightforward, affordable, quick to perform, and requires simple laboratory equipment. The flocculation is achieved at 3.5-4.0 pH, where the casein proteins (net positive charge) interact with viral particles from the water samples, which carry a net negative charge due to their functional surface groups like carboxylates and phosphates. The reaction solution is usually agitated at low rpm for approx. 2h, during which through electrostatic interactions, flocs (virus-protein complexes or aggregates) are formed and settle out of the solution either naturally or facilitated through centrifugation at 3500 x g for 30 min [7]. Previously, SMF has been used for simultaneous concentration and quantification of water-borne viruses, bacteria and protozoa [8] mainly for water control purposes and microbial risk assessment studies. However, there is a literature gap regarding its applicability in bacterial microbiome studies. Therefore, we used both 16S rRNA and shotgun metagenomics data to observe microbial differences in the taxonomic composition to assess the applicability of SMF for bacterial metagenomics.

Another critical step is the gDNA extraction from filters or precipitates. Numerous studies have assessed the impact of various DNA extraction protocols on the quantitative analysis of bacterial biomass [9,10,11,12,13], highlighting the necessity for a standardized protocol yielding reproducible results to facilitate cross-study comparisons, especially in research areas with future health diagnostic perspectives such as human gut microbiota profiling. For instance, according to the largest current comparative study based on shotgun sequencing evaluating the bacterial extraction performance of 21 faecal DNA extraction protocols [9], the Protocol Q, which is a slightly modified version of Qiagen’s QIAamp DNA Stool Mini Kit, has been proposed as the standard protocol providing the best results for bacterial DNA extraction from human faeces. To date, there is no equivalent DNA extraction protocol for aqueous environmental samples. In addition, due to limited availability or costs, not all studies could afford to use the same protocol, especially commercial ones. For this, we optimised an in-house protocol that involves elements from other published protocol and compared it with four other commercial DNA extraction kits on Mock community standard samples to inform future study designs.

In the case of targeted bacterial metagenomics, after sequencing and the initial quality control of the sequencing data, 16S rRNA databases (DB) are utilised to infer bacterial taxonomic composition. While de novo clustering approach is often the choice for studies for initial observation it is not directly influenced by 16S databases, closed or open reference clustering approaches utilise 16S rRNA DBs for clustering and are still preferred widely. As a result, 16S rRNA DBs influence the final taxonomic composition. Several studies benchmarked their performance previously [13,14,15,16], but they undergo constant updates which prompts additional evaluation. Herein, we present a small-scale easy-to-perform benchmarking of six well-established as well as more recently published 16S rRNA DBs on the same set of Mock samples extracted by the five different DNA protocols to evaluate the extent of bias introduced by the choice of 16S rRNA DBs.

This study showed results from comparison of different techniques spanning the three most common critical steps in each metagenomics study: the applicability of SMF for the use of bacterial metagenomics on real samples; the evaluation of a novel DNA extraction method on Mock community standard and on real samples through comparative analysis with other commercially available DNA extraction protocols; and comparison of 16S rRNA databases for taxonomic assignments.

2. Results

First, we evaluated the bias introduced by different DNA extraction protocols and the 16S rRNA DBs on the 16S metagenomics Mock community datasets (n=8). In total, we scored factors such as DNA yield, A_260/280nm, A_260/230nm, Species-level Taxon accuracy rate (TAR), Genus-level Taxon detection rate (TDR), Measurement Integrity Quotient score (MIQ), percent of reads that mapped to the reference sequences, failed quality filter, failed to merge, or are chimeric (Supplementary Table S6). The in-house method yielded the highest amount of DNA, whereas highest purity was achieved with the EurX and EZNA kits. EurX kit yielded an overall good results, compared to the other mock samples and at the time of assessment was considered for the shotgun studies.

2.1. MIQ Score Variability in DNA Extraction Protocols

The tool miqScore16SPublic (2.6.) provided by Zymo was used to establish the bias introduced by DNA kits. The complete reports are available in Suppl. File S1, and only the radar plots were provided in Figure 1. Results indicate that the samples isolated by the FastSpin Soil kit showed the least bias scoring 88 MIQ, followed by EurX, and the in-house method, whereas the worst MIQ scores were in Zymo kit and its datasets with the varying annealing temperatures.

2.2. Bias Introduced by 16S rRNA DBs and/or DNA Extraction Kits on Closed-Reference Clustered Mock Samples

Herein we assessed to what extent the DB could bias the taxonomic composition. For this, closed-reference 99% OTU clustering was applied by using each DB reference sequences. Taxonomic composition plots of all mock samples are visualized on Figure 2. The DBs GTDB-full and Silva generated more OTUs compared to the remaining DBs. This could be considered both beneficial and negative depending on the study purpose. However, in this case the additional OTU Citrobacter_B in GTDB-full that should not be present in the mock sample is skewing the results, resulting in lower MIQ scores (Figure 2 and Figure 4).

We compared the 16S rRNA DBs and the DNA extraction protocols by MIQ scores in parallel and the results are shown in Figure 3. First, 16S rRNA DBs were compared (Figure 3A), and the best theoretical compositions with the least bias were yielded by GG 13_8, followed by GSR and GTDB-214.1-full. Interestingly, GreenGenes database yielded an OTU at family level Enterobacteriaceae, failing to differentiate Escherichia coli and Salmonella enterica. Unfortunately, there was no adequate approach to separate this 99%-clustered OTU for the MIQ calculator to precisely calculate the bias for each microorganism. Therefore, for GG_13.8 specifically, we considered both E. coli and S. enterica as a single organism, and adjusted the reference theoretical composition in the MIQ score script. As a result, all the MIQ scores for GG_13.8 DB are elevated compared to the remaining DBs and should be interpreted with caution. We then compared different DNA extraction kits (Figure 3B), regardless of DB used. Except for ZymoBIOMICs with primer annealing temperature 55°C all the Zymo variants performed worse compared to the others. FastSpin kit yielded the highest MIQs scores (82.6 on average) with all DBs. Second best results were scored by our in-house protocol.

Interestingly, all the 16S DBs except GG_13.8 yielded comparable results with each DNA extraction kits ranging from 59-71 in the worst performing sample (Zymo-62C) up to 80-83 in the best sample (FastSpin). These results suggest that the DNA extraction kit had the largest impact on the final mock composition. The bias in mock samples could not be compensated with the use of a better performing 16S DB (as in the case of both repetitions Zymo-62C in Figure 3C) and vice versa – a sample treated with good performing DNA extraction kit yields comparable taxonomic compositions resembling the theoretically expected one regardless of the 16S DB used.

Figure 2. Taxa bar plots of the eight mock samples analyzed with closed-reference clustering with each 16S rRNA Database separately and de novo clustered datasets identified with GTDB-full. A database-specific legends with taxa identification (lowest rank) were provided on the right side. For easier interpretation the resulting MIQ score was provided right next to the bars. Higher MIQ score represents taxonomic composition closer to the theoretically expected one (bottom).

Figure 2. Taxa bar plots of the eight mock samples analyzed with closed-reference clustering with each 16S rRNA Database separately and de novo clustered datasets identified with GTDB-full. A database-specific legends with taxa identification (lowest rank) were provided on the right side. For easier interpretation the resulting MIQ score was provided right next to the bars. Higher MIQ score represents taxonomic composition closer to the theoretically expected one (bottom).

Figure 3. Performance of DNA extraction kits and 16S DBs by MIQ score. A. Comparison of 16S DB. B. Comparison of DNA extraction kits/annealing temperatures. Green is better and red is worse. C. Bar graph with the MIQ scores and average root mean square error (RMSE) values.

Figure 3. Performance of DNA extraction kits and 16S DBs by MIQ score. A. Comparison of 16S DB. B. Comparison of DNA extraction kits/annealing temperatures. Green is better and red is worse. C. Bar graph with the MIQ scores and average root mean square error (RMSE) values.

Figure 4. Comparison of taxonomic composition of SMF and VF. (A) Duplicate datasets from the amplicon metagenomics samples from River Perlovska per DNA Extraction kit. (B) Datasets from the shotgun metagenomics samples from River Iskar. Sample pairs (VF and SMF treated) were ordered by sampling date and location (M - villa area Mechka, D - Dragushinovo village). Differential abundance analysis of the (C) amplicon and (D) shotgun samples with all VF samples pooled as a reference compared to the SMF samples. Only taxa with p value < 0.5 and with log fold change (LFC) ≥ 1.0 (C) and 0.5 (D) were presented.

Figure 4. Comparison of taxonomic composition of SMF and VF. (A) Duplicate datasets from the amplicon metagenomics samples from River Perlovska per DNA Extraction kit. (B) Datasets from the shotgun metagenomics samples from River Iskar. Sample pairs (VF and SMF treated) were ordered by sampling date and location (M - villa area Mechka, D - Dragushinovo village). Differential abundance analysis of the (C) amplicon and (D) shotgun samples with all VF samples pooled as a reference compared to the SMF samples. Only taxa with p value < 0.5 and with log fold change (LFC) ≥ 1.0 (C) and 0.5 (D) were presented.

2.3. Taxa Identification Efficiency of 16S DBs on de novo Clustered Mock Samples

While taxonomic composition is the first criterion to assess the bias in mock samples, taxon accuracy rating (TAR) and taxon detection rate (TDR) are two important factors which essentially, evaluate how a 16S rRNA DB performs the taxonomical identification. Although, we presented the results as values for each DB, one should note that they are heavily influenced by the amplified region, the choice of OTU vs. ASV approach and should not be considered entirely as drawbacks to the 16S DBs.

The taxon accuracy rating (TAR) and taxon detection rate (TDR) results were calculated with data from 99% de novo clustered OTUs. The 99% OTU-generated results tend to underrepresent Enterococcus faecalis by splitting it into two separate OTUs g_Enterococcus and s_Enterococcus-faecalis with similar relative frequency values in all extraction kits, thus skewing the total microbial composition. For a small Mock community standard with only 8 bacterial strains, results did not vary drastically, and also TAR and TDR scores at each taxonomic level were identical. At species-level the accurate identification varied from 1/8 bacteria for Silva, 2/8 for GSR, and up to 4/8 for Ezbio, while for genus-level was between 6/8 to 8/8 (Figure 2 and Suppl. Table S7). GreenGenes 13_8 and GTDB 214.1 were designed for genus-level identification, therefore species-level resolution was not possible and not discussed. Only one case of misclassification was recorded, namely Listeria monocytogenes identified as Listeria ivanovii by Ezbio. Underclassifications by two taxonomic levels were observed for Salmonella enterica in GTDB and Silva, for Escherichia coli in GreenGenes and GSR. Additionally, P. aeruginosa was underclassified by two ranks by GreenGenes (Suppl. Figure S7).

2.4. Evaluation of Skimmed Milk Flocculation

SMF and VF were compared only on real samples as the Zymo Mock standards are not designed to be pre-treated prior to DNA extraction as cells are stored in DNA/RNA Shield^TM and are partially lysed. The relative taxonomic composition of both types of datasets were presented in Figure 4A and Figure 5B. Full length figures of the DAA results are available in Supplementary Figure S2 and S3. In all sample pairs (a pair being VF and SMF treated) noticeable differences were observed by the separation of SMF- and VF-treated samples in the PCoA plots (Supplementary Figure S4).

In the 16S rRNA amplicon datasets (River Perlovska), genus Polaromonas, Agitococcus and an unidentified OTU were observed to be overrepresented in the SMF treated samples. In Green Genes and the GSR DBs the Polaromonas was identified to species-level as P. naphtalenivorans, while the unidentified OTU was either left at domain Bacteria, order Bacteroidales, or identified as family Saprospiraceae by EzbioCloud (both the free DB v2018 and their website non-free latest DB v2023.08.23). All three taxa and a few more were also shown to be significantly enriched (p < 0.5) with log fold change (LFC) of 2.0 and above as shown by the DAA in Figure 4C.

Similar trend was observed in the shotgun metagenomics datasets (River Iskar) but with different enriched taxa. Not all enriched taxa were visible on the bars except for the genera Streptococcus (brown) and Lactococcus (pale yellow) in SMF-S2, SMF-S3 and SMF-S4 in Figure 4B. Sankey plots were provided for better visualization of all taxa in Supplementary Figure S1. According to the DAA the lactic acid bacteria members Lactococcus, Leuconostoc, Streptococcus, Enterococcus, and Lactobacilus were significantly enriched (p < 0.5) with LFC between 1.75 - 5.0 in the SMF-treated samples. Lastly, the genus Macrococcus was also significantly overrepresented and Bracken species-level hits are mostly Macrococcus caseolyticus. Interestingly, genus Streptococcus was also detected in the amplicon metagenomics samples but was not overrepresented. The complete non-filtered taxonomy tables are available as supplementary Table S7, S8 and S9.

3. Discussion

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn. This is the first study to evaluate the applicability of SMF in bacterial metagenomics. Despite the small scale of the study design, with the use of DAA analysis it was clearly shown that SMF has skewed the taxonomic composition of real water samples, therefore rendering the SMF technique inapplicable for bacterial enrichment in metagenomics. Interestingly, not all taxa were altered, rather only specific ones. According to a previous study on the concentration of specific species such as Escherichia coli and Helicobacter pylori SMF has provided comparable results as expected [8] and could be used for detection of specific pathogens or indicator species.

Skimmed milk consists mainly of protein (casein and whey) and lactose in addition to other nutrients and minerals that could act as growing factors for bacteria and is usually supplied in culture media. In this regard, it was previously shown that the genus Polaromonas was enriched in dairy products removal tanks [17], while other studies have identified it as the third most abundant genus in the mixed-species dairy biofilm on a biofilter [18]. Genus Agitococcus was also significantly enriched in SMF-treated samples and while no species level identification was achieved, Agitococcus lubricus, a species described in 1981 was tested positive for skim milk proteolysis [19], and members of the genus Agitococcus would likely be positive as well. Unfortunately, the largest significantly enriched OTU namely the f_Saprospiraceae (EzbioCloud) or o_Bacteroidales (Green Genes) or d_Bacteria (GTDB-full) remained unidentified and no conclusion about its SMF utilization could be drawn. Similarly, in the shotgun datasets from River Iskar, all of the enriched lactic acid bacteria are usually found in decomposing plants and milk products, which produce lactic acid as the major metabolic end product of carbohydrate fermentation by utilizing the lactose from the skimmed milk. The acidification of samples (pH=3.5) during the SMF which facilitates the flocculation process might be advantageous to their replication. Lastly, according to the Bracken reports, most of the Macrococcus read hits were Macrococcus caseolyticus, which is again shown to efficiently hydrolyse casein and is a natural component of the secondary microflora in cheeses and sausages [20,21]. In conclusion, the skimmed milk was highly likely metabolized during the 2 h. incubation protocol resulting in replication of specific taxa and skewing significantly the final taxonomic composition of SMF-treated samples.

DNA extraction is a critical step in a metagenomics workflow and it is known to be influenced by numerous parameters, which are challenging to evaluate comprehensively. In addition, the choice of the DNA extraction method strongly affect the detection and composition of bacterial communities [22,23]. Moreover, in-house methods and protocols and or commercial products are constantly being developed and used widely, therefore cross-study comparison is difficult and either updated benchmarking studies are required or standardization. While, we developed a good performing DNA extraction protocol, further improvements are needed to match the performance of the EurX kit used here. Nevertheless, budget study design could implement the in-house method. Surprisingly, EZNA Universal Pathogen kit is not designed for metagenomics and certainly is not optimized to extract DNA equally from Gram positive and Gram negative as the cell wall of Gram-positive bacteria is composed of a thick layer of peptidoglycan. Nevertheless, it performed more or less equally to other kits such as ZymoBIOMICs, EurX without implementing a bead-beating step, which is presently widely adopted and recommended to facilitate balanced lysis [24].

Next, our meticulously optimized 16S rRNA amplicon library sequencing protocol has shown good results, yielding MIQ scores ≥ 80 with most of DNA Extraction kits and the in-house method, which according to the Zymo miqScore16SPublic calculator were considered good and validates the real sample results. Nevertheless, benchmarking of sequencing datasets from simultaneous 16S amplicon libraries generating protocols or commercial kits are required to fully evaluate the applicability of the 16S protocol.

In regard to the 16S DBs comparison, we aimed to present the most sample- and primer pairs-specific taxonomic identification by first truncating the reference sequences to the primer regions with the use of in-silico PCR and then building a classifier. This way, the detection and the accuracy of identification of each 16S DB have been adjusted specifically to the used primer pair, thus allowing for normalized comparison. TAR and TDR values were not as informative as initially perceived mainly due to the small number of bacteria included in the Microbial community standard. 16S DB results presented here should be interpret together with the amplified region as identification is also heavily influenced by the amplified region. TAR of resulting OTUs clustered at threshold equal or below 99% usually suffered from identification bias [25] and as anticipated no 8/8 TAR was achieved with OTUs. While most reliable identification is usually achieved with amplicon sequence variants (ASV), a large portion of studies still rely on OTU clustering [26,27].

A few limitations of the study could be listed. The SMF protocol applied in this study was originally optimized for virus concentration. In the literature, SMF protocols adopted or adjusted for bacteria are lacking and additional pre-treatment steps could be implemented to inhibit bacteria growth. Ideally, microbial community standards with higher number of bacteria (20+) would present more insightful results compared to the 8 bacteria community standard used here.

4. Materials and Methods

Both 16S and shotgun metagenomics sets were included in the current study. Chronologically, 16S rRNA amplicon sequencing was performed first. Initially, a novel in-house method for DNA extraction for metagenomics was optimized. Then, Mock Community Standard together with real water samples were used to evaluate and compare the new protocol for metagenomics to other commercial DNA extraction kits followed by a small-scale 16S rRNA databases comparison. Next, the applicability of SMF for both 16S amplicon and shotgun metagenomics was analysed on real water samples only and separately from the comparative studies on mock samples. For the shotgun samples, only one DNA extraction protocol was applied.

4.1. Samples

Two different sets of water samples were collected and processed separately. For shotgun metagenomic sequencing four composite water samples (1 litre each) were collected in pairs along the River Iskar from the two locations (42.367698, 23.555463 – Dragushinovo village and 42.431095, 23.531900 – villa area “Mechkata”) with an automatic sampler over the period of 24 h to avoid day-night bias on 3rd Nov 2022 and 17th Nov 2022. They were transferred to the laboratory within 6 h and immediately processed. For 16S rRNA amplicon sequencing one non-composite water sample (1 litre) was collected in a sterile HDPE plastic container from a small urban River Perlovska at location 42.692164, 23.343892 and transported within 30 min to the laboratory. All samples were divided into two equal parts (500 ml) for SMF and VF steps as described in Figure 5. SMF and VF was applied on both set of samples to analyse the applicability of SMF in both 16S and shotgun bacterial metagenomics.

Figure 5. Graphical representation of the SMF and VF water treatment steps. In brief, the sample was divided into equal parts. For SMF, the sample was adjusted to 0.05% skimmed milk, with pH 3.5-4.0, agitated at slow speed for 2 h, aliquoted and centrifuged at 3500 x g for 30 min. For VF, the sample was further divided into 100 ml portions and filtered through 0.2 µm nylon filters, aseptically cut into smaller pieces and added to DNA extraction tubes. A pair of two tubes with pellets / filters were resuspended in lysis buffer from each kit / protocol.

Figure 5. Graphical representation of the SMF and VF water treatment steps. In brief, the sample was divided into equal parts. For SMF, the sample was adjusted to 0.05% skimmed milk, with pH 3.5-4.0, agitated at slow speed for 2 h, aliquoted and centrifuged at 3500 x g for 30 min. For VF, the sample was further divided into 100 ml portions and filtered through 0.2 µm nylon filters, aseptically cut into smaller pieces and added to DNA extraction tubes. A pair of two tubes with pellets / filters were resuspended in lysis buffer from each kit / protocol.

4.2. Skim Milk Flocculation

For SMF a previously described protocol was used [7]. In brief, 5% skim milk (HiMedia Laboratories) was autoclaved for 15 minutes at 115°C, 18 psi. Five ml of the 5% pre-flocculated skim milk solution was added to 500 ml sample to achieve 0.05% final skim milk concentration. The sample pH was adjusted to 3.5-4.0 with 1M HCl, placed on a horizontal shaker and agitated at 200 rpm for 2 h at room temperature. It was then distributed into 50 ml conical tubes, and centrifuged at 3500 x g for 30 min at 4⁰C. The supernatant was decanted and tubes left upside down to drain residual water for 5 min. Pellets were subjected to DNA extraction (see 2.4.).

4.3. Vacuum Filtration

The samples were subjected to VF through 47 mm diameter 0.2 µm pore size nylon filters (Cytiva, Whatman™), using Lafil 400 - LF 30 Filtration System (Rocker). The sample from River Perlovska (16S amplicon metagenomics) was divided among five filters (100 ml each) and subjected to the five DNA extraction kits (2.4.), while for River Iskar samples (shotgun metagenomics), the entire volume of 500 ml was filtered through one filter and extracted with only one DNA extraction kit (see 2.4). Regardless of the sample sets, each filter was cut sterilely in two equal halves to be further extracted in pairs. Each half was further cut into smaller pieces for better homogenization and directly added to extraction tubes (see 2.4.).

4.4. DNA Extraction

The four DNA extraction kits and the in-house protocol, were evaluated in parallel only on the Perlovska River and Mock samples in duplicates. The following kits used were by adhering to the manufacturer instructions: 1) E.Z.N.A. Universal Pathogen Kit (OMEGA Bio-Tek, Inc., USA); 2) ZymoBIOMICS DNA Miniprep Kit (Zymo Research, USA); 3) FastDNA Spin Kit for Soil (MP Biomedicals, LLC, Germany); and 4) Environmental DNA & RNA Purification Kit (EURx Sp. z o.o., Poland). The in-house protocol was developed by combining elements from other published protocols [28] with slight modifications and optimised locally. Complete detailed protocol and reagents are provided in Suppl. Tables S1. The Microbial Community Standard used herein was (ZymoBIOMICS, Zymo Research, cat. D6300), from which we obtained five DNA samples. The Perlovska River water sample was split into SMF-treated and VF-treated portions, and each was subjected to five DNA extraction methods, resulting in ten DNA samples (Figure 5). Next, DNA from the Iskar River samples were extracted with Environmental DNA & RNA Purification Kit (EURx Sp. z o.o., Poland).

4.5. 16S rRNA and “Shotgun” Metagenomics

16S rRNA metagenomics was applied on Perlovska river samples and Mock community standard samples. The 16S rRNA V3-V4 region was amplified with previously published primer pairs Pro341F and Pro805R [29] and with a rigorously optimized 16S rRNA amplification protocol described in Suppl. Tables S2÷S5. Sequencing was carried out on Illumina MiSeq V3 (2x300bp). Additionally, for the Zymo miniprep DNA kit we used the same DNA sample in three additional amplification reactions with primer annealing temperature gradient (55℃, 58.5℃, 62℃), therefore resulting in 8 Mock Community Standard and ten Perlovska River datasets.

Shotgun metagenomics was applied on Iskar River samples and a single Mock community standard (used for pre-sequencing quality control). NGS libraries were constructed Illumina DNA Prep kit (Illumina, Inc., San Diego, California, U.S.) with 50 ng genomic DNA input. The libraries were pooled and sequenced on NextSeq 550 with V2.5 Mid Output kit (2 x 150 bp) (Illumina, San Diego, California, U.S.).

4.6. Bioinformatic Analysis

For DNA extractions and 16S rRNA DBs comparisons with Mock samples we used the Qiime2 platform [30] and miqScore16SPublic tool (https://github.com/Zymo-Research/miqScore16SPublic, accessed on 22.07.2024) to determine which DNA extraction kit yields microbial composition of mock samples closest to the “ground truth”. Raw reads were automatically demultiplexed and trimmed from adapters in BaseSpace. Next, cutadapt v4.6 was used to remove the 16S rRNA primers and trim low quality bases (<20 Q at 5’end, and <15 Q at 3’end). We used six 16S rRNA DBs and trained classifiers (n=6) locally. First, full-length 16S DBs (Silva 138-99%-OTU, GSR-99%-OTU, GTDB-214.1-99%-OTU-less, GTDB-214.1-99%-OTU-full, GreenGenes-13_8-99%-OTU, EzBio-Cloud-v.2018) were downloaded. Reference sequences were extracted based on in-silico PCR with the 16S primer pairs, and used as input for the Naïve-bayes algorithm within sk-learn to train personalized classifiers with default parameters, used later for taxonomic assignments with confidence index at 0.7 (default). Both de novo and closed-reference clustered (for each 16S rRNA DB) OTUs were obtained at 99% for each mock dataset. De novo clustered OTUs provided an equal starting basis to efficiently compare the identification accuracy and detection limits of each 16S rRNA DBs, as the taxonomic composition of all resulting OTUs is the same regardless of the DB. For this we used the “evaluate-composition” plugin [31] in Qiime2 to compare the expected and observed features, the taxon detection and taxon accuracy rates.

On the other hand, closed-reference OTUs are clustered according to each 16S rRNA DB and thus influence the final taxonomic composition. To evaluate the composition the miqScore16SPublic tool was used, as it calculates a MIQ score (ranging from 0 to 100), a numeric value to assess the bias between the observed and the theoretically expected microbial profile. However, it was designed to function only with raw reads and performs amplicon sequence variants (ASV) instead of OTUs. And while we used to as is and considered the obtained MIQ scores as reference, it could not be used with OTU tables obtained from our closed-reference clustered OTUs. For this we used the miqScore16SPublic intrinsic formula for MIQ score calculation in a separate Python script (https://github.com/maddne/MIQ-calc-from-OTU-tables), which accepts OTU tables as input regardless of their origin.

Next, for comparison of SMF and VF we used differential abundance analysis (DAA). Taxonomic assignment to raw reads was performed with Kraken2 v2.1.2 [32] with the PlusPF (Standard plus Refeq protozoa & fungi), built in January 2023, followed by Bracken v2.8 [33] as described in protocol [34]. The resulting OTU count tables were imported into Qiime2 platform, filtered from low count taxa (n=20 in at least 3 samples), and used for differential abundance analysis with ANCOM-BC2 plugin [35].

5. Conclusion

In this study, we systematically evaluated the effectiveness of SMF for bacterial metagenomics, introduced and benchmarked a novel in-house DNA extraction method against four commercial kits, and assessed the performance of six 16S rRNA databases for taxonomic identification.

The findings reveal that skimmed milk flocculation is unsuitable for bacterial microbiome studies as it significantly alters the microbial composition due to the proliferation of lactic acid or casein utilizing bacteria, leading to an increased relative abundance compared to the traditional vacuum filtration method. Our in-house DNA extraction protocol demonstrated competitive performance, particularly in comparison to the commercial kits, which were rigorously optimized for minimal bias. This in-house protocol provides a cost-effective alternative for researchers with limited access to commercial kits, offering robust and reliable results for metagenomic studies.

Lastly, the evaluation of 16S rRNA databases showed that while there are variances in taxonomic assignments, the choice of DNA extraction protocol has a more pronounced impact on the microbial composition than the choice of 16S rRNA database. This underscores the importance of selecting an appropriate DNA extraction method to minimize biases in metagenomic studies. This comprehensive benchmarking study offers critical insights for designing future water metagenomic studies, emphasizing the importance of method selection at various stages to ensure accurate and reliable microbial community profiling.