To compare simultaneous bacterial and fungal DNA extraction effectiveness, we used two types of fecal material, (i) germ-free mice feces, for their absence of viable microbiota and very low background microbial DNA load originating from diet and bedding sterilized by irradiation (Fontaine et al., 2015; Schwarzer et al., 2017), confirmed also in our experiments (Supplementary Table S1), and (ii) a single human stool stock representing a natural sample type control to ensure an identical microbial composition in all samples. Germ-free feces were inoculated with known quantities of clinical bacterial (Enterococcus faecalis) and fungal (Candida albicans and Aspergillus fumigatus) strains to evaluate the methods’ yields using real-time PCR. The strains were chosen as typical representatives of gram-positive, yeast and filamentous microorganisms, respectively. Human feces serving as vehicle controls were analyzed to compare DNA extraction methods in terms of integrity, purity, DNA quantity and microbial composition. Use of the single human stool sample was intended to avoid the effect of inter-individual variability observed elsewhere (Huseyin et al., 2017). The detailed study design is shown in Figure 1.

Germ-free mice fecal samples were obtained from the Laboratory of Gnotobiology, Institute of Microbiology, Academy of Sciences of the Czech Republic. This laboratory is well-established in germ-free animal research (Kozakova et al., 2016; Schwarzer et al., 2016). A single human stool sample was donated by a healthy volunteer who signed informed consent in accordance with the Ethic Committee of the Centre for Cardiovascular Surgery and Transplantation (Protocol no. 201603).

Immediately after delivery, two grams of mice feces were homogenized using TissueRuptor^® (Qiagen, Germany) in 18 mL of autoclaved PBS, divided into 160 μl aliquots (n = 85) and frozen at -80°C until used. Then, the stool aliquots were separated into groups for two independent sets of DNA extraction. For the first set, the aliquots (n = 30) were simultaneously spiked with 20 μl of two different C. albicans concentrations (6.9 × 10⁴ or 6.9 × 10⁵ cells/ml; low or high concentration, respectively) and 20 μl of E. faecalis concentrations (7.7 × 10⁵ or 7.7 × 10⁷ cells/ml), resulting in a 200 μl total sample volume (Figure 1). For the second set, the aliquots (n = 30) were spiked with 20 μl of two different amounts of A. fumigatus (10⁶ or 10⁸ cells/ml approximately) and 20 μl of sterile water (B. Braun Medical, Inc., Germany) to preserve an equal sample volume. Stool aliquots spiked with 40 μl sterile water (n = 25) were included in both groups as the baseline microbial DNA load controls.

In addition, we used 200 mg of human feces, serving as a natural sample type control for the third independent DNA extraction set. The stool was homogenized in 4 mL autoclaved PBS, divided into 200 μl aliquots (n = 15), and frozen at -80°C until used. Dilutions were performed to make samples homogeneous and pipettable to ensure exactly the same volume (the same amount of microbiota) in all samples and thus, identical conditions for each extraction method tested.

Five different fecal DNA extraction methods were evaluated in this study: QIAamp DNA Stool Mini Kit with a pre-treatment step (QIA; Qiagen, Germany), PureLink^TM Microbiome DNA Purification Kit (PL; Thermo Fisher Scientific, United States), ZR Fecal DNA MiniPrep^TM Kit (ZR; Zymo Research, United States), NucleoSpin^® DNA Stool Kit (NS; MACHEREY-NAGEL GmbH & Co., KG, United Kingdom) and non-commercial protocol Q, recommended by the International Human Microbiome Consortium (abbreviated as IHMS). All protocol procedures were performed according to the manufacturer’s instructions with minor differences (detailed protocols are shown in the Supplementary Material). Three independent sets of DNA extraction were performed (Figure 1), and all experimental samples were processed in triplicates for the reproducibility evaluation. Samples with sterile water were used as a negative reagent and process contamination controls. Altogether, 115 samples and reagent controls were extracted in all sets.

DNA concentration was determined fluorometrically on the Qubit^® 3.0 Fluorometer (Thermo Fisher Scientific, United States) using the Qubit^TM dsDNA HS Assay Kit. DNA purity was determined via 260/280 and 260/230 ratios measured on the NanoDrop 1000 (Thermo Fisher Scientific, United States). DNA integrity was determined by 0.6% agarose gel electrophoresis and visualized (Supplementary Figure S1).

Commercial real-time PCR kit for specific E. faecalis detection (PrimerDesign^TM genesig, United Kingdom) was used to evaluate the bacterial yield using the ABI 7500 fast Real-Time PCR System (Applied Biosystems, United States).

The fungal yield was evaluated using a broad-range real-time PCR assay as described previously (Nemcova et al., 2015) using the Rotor-Gene 6000 (Corbett Research, Austria). Briefly, amplification of the ITS2 region using primers UNF1 (5′-GCATCGATGAAGAACGCAGC-3′) and UNF2 (5′-TTGATATGCTTAAGTTCAGCGG-3′) was performed with 2 × SensiFAST HRM mix (Bioline, United Kingdom) and RNase-Free Water (Qiagen, Germany).

In both real-time PCR assays, all the analyzed samples were performed in duplicates with serially diluted calibration curves. Sterile water served as no template control. Pure culture extracts served as positive controls for species verification. Altogether, 394 real-time PCR reactions for germ-free fecal samples and 60 real-time PCR reactions for human fecal samples (controls included) were performed.

Statistical analyses were conducted and visualized in R (3.4.2). One-way analysis of variance (ANOVA, function aov in R) followed by Tukey’s multiple comparison test (function Tukey HSD in R) was used to evaluate the DNA yield and real-time PCR data. p-values lower than 0.05 were considered significant. Log transformed data were used for fold change analysis.

Microbial communities were profiled by 16S rDNA and ITS1 rDNA amplicon sequencing using the Illumina MiSeq sequencing platform. Primers with unique barcode sequences for PCR amplification over the bacterial 16S rDNA gene’s V3/V4 region were designed as described previously (Kubasova et al., 2017). Primer pairs ITS1F/ITS2, recommended by the Earth Microbiome Project¹, were used with unique barcode sequences designed in this study, to amplify over the fungal internal transcribed spacer region 1 (ITS1) of the rRNA operon (Supplementary Table S2).

16S Library was constructed according to the “16S Metagenomic Sequencing Library Preparation protocol” recommended by Illumina, with minor differences. Briefly, PCR was performed with 2 × KAPA HiFi HotStart ReadyMix (Kapa Biosystems, Inc., United States) under the following conditions: initial denaturation at 95°C for 15 min, followed by 30 cycles consisting of denaturation at 95°C for 40 s, annealing at 55°C for 45 s and extension at 72°C for 60 s, with a final extension step at 72°C for 5 min. From each reaction, 5 μl were analyzed by 2% agarose gel electrophoresis. PCR products were purified using AMPure XP magnetic beads (Beckman Coulter Inc., United States) diluted into an equimolar concentration and pooled according to their unique barcode sequence, which enables multiplexing. Only samples with different barcode sequences were pooled together. Next, Illumina dual-index barcodes were added to the pooled PCR products with the Nextera XT Index Kit (Illumina, United States). Indexed PCR products were purified and pooled into the equimolar concentration prior to paired-end sequencing with MiSeq Reagent Kit v3 (600 – cycle) (Illumina), following the manufacturer’s directions.

The ITS1 Library was constructed in a similar way to the 16S Library, with the following modifications: PCR was performed with HotStarTaq Master Mix (Qiagen, Germany), 4 mM MgCl₂ (Thermo Fisher Scientific, United States), RNase-Free Water under the following conditions: initial denaturation at 95°C for 15 min, followed by 40 cycles consisting of denaturation at 95°C for 30 s, annealing at 56°C for 45 s and extension at 72°C for 60 s, with a final extension step at 72°C for 10 min. In the first clean up step, we increased the bead amounts to 35 μl per sample to match the input volume. Then the purified PCR products were diluted to an equimolar concentration and samples with different barcode sequences were pooled together. Next, Illumina dual-index barcodes were added to the pools with the Nextera XT Index Kit (Illumina, United States). Indexed PCR pools were purified and pooled into the equimolar concentration prior to paired-end sequencing with MiSeq Reagent Kit v3 (600 – cycle) (Illumina) following the manufacturer’s directions.

For both libraries, raw sequence data analysis was carried out using QIIME (v. 1.9.1.) pipeline (Caporaso et al., 2010). 16S and ITS1 read pairs were demultiplexed based on the unique barcode sequence and then merged using the default QIIME script (join_paired_end.py). Chimeric sequences were identified using VSEARCH (v. 2.6.1) (Rognes et al., 2016) with Greengenes reference database (v. 4feb2011) for bacteria and UCHIME (v. 7.2) reference dataset for fungi. Sequences were clustered into OTUs at 97% threshold using VSEARCH 2.6.1 de novo. A set of sequences representing OTUs was created, and taxonomy was assigned (using script: assign_taxonomy.py) to each sequence using the Greengenes database (v. gg_13_8_otus) and Uclust (v. 1.2.22q) (Edgar, 2010) for bacteria, and using BLAST and UNITE (v. 7.2)² for fungi (input sequences were searched for against a BLAST database of pre-assigned reference sequences from UNITE). This process resulted in the OTU table in BIOM format with the singletons discarded. The bacterial OTU table was filtered for OTUs, with the number of sequences less than 0.005% of the total number of sequences (Bokulich et al., 2013). The fungal OTU table was not further filtered. PyNAST (v. 1.2.2.) was used to align representative sequences to build a phylogenetic tree using FastTree (v. 2.1.3) in bacterial analysis. The non-phylogenetic metric (i.e., Bray–Curtis dissimilarity distance) was calculated for fungi due to the inapplicability of phylogenetic-based metrics (i.e., Weighted/Unweighted UniFrac distance) for ITS1 sequence analysis (Halwachs et al., 2017). Alpha- and beta-diversity calculations were performed and visualized with QIIME script core_diversity_analyses.py.

Statistical analyses were performed using online Calypso software (version 8.82) (Zakrzewski et al., 2017). Data were normalized using cumulative-sum scaling (CSS) and log2 transformation. Beta-diversity calculations were visualized using the principal coordinate analysis plots (PCoA), based on unweighted UniFrac, weighted UniFrac (bacterial data) and Bray–Curtis (fungal data) distances, and compared using the nonparametric analysis of similarities (ANOSIM) test. Moreover, the nonparametric Kruskal–Wallis and generalized linear model (GLM) tests implemented in ALDEx2 test were applied to detect differences in taxa abundances at genera level. These tests were performed with only the higher abundance taxa (>0.01% of total). The p-values were adjusted according to the Benjamini–Hochberg procedure. p-values lower than 0.05 were considered significant. To characterize method specific composition profiles, “core microbiome” measurements were performed.

To compare the DNA extraction methods, an artificial system was constructed. The homogenized germ-free mice stool aliquots were spiked with verified clinical bacterial and/or fungal cultures, each in two defined concentrations differing in one to two orders of magnitude. Non-spiked stool aliquots were used as a baseline microbial load, and sterile water was used as a blank control. All extractions were performed in triplicates and compared using bacterial- and fungal-specific real-time PCR, resulting in 394 reactions.

All methods varied in terms of technical reproducibility, but the variability among replicates was considerably lower than among the same samples extracted by different methods. Extraction from the low-concentrated (LC) and high-concentrated (HC) communities showed similar patterns (Figure 2). However, microbial DNA recovery was significantly influenced by the DNA extraction method.

A real-time PCR analysis of the extracted bacterial DNA revealed a similar amount of extracted DNA in IHMS, PL, QIA and ZR methods, but a significantly lower DNA yield for the NS method (LC and HC level; p < 0.001 and p < 0.05, respectively) (Figure 2A). Moreover, real-time PCR analysis indicated moderate PCR inhibition by extraction with the ZR method (Supplementary Figure S2).

Extracting fungal DNA yielded variable results between the methods used (Figure 2B,C) and between C. albicans and A. fumigatus, selected as typical yeast and mold representatives. For both species, the IHMS and ZR methods provided the highest DNA yields consistently (Figure 2B,C). The significantly lower yields (p < 0.001) were produced by QIA and PL methods in both concentration levels. The QIA method produced higher DNA yields for C. albicans (Figure 2B) and the PL method was more successful in A. fumigatus extraction (Figure 2C). The NS method results were inconsistent and were species and load-dependent (Figure 2B,C). It is also notable to mention that the ZR method’s profile was not observed as being inhibited in any fungal assays.

During E. faecalis DNA detection, a positive signal was occasionally captured in one of the baseline controls’ real-time PCR duplicates (once in both, n = 5/15) across methods, with values lower than 1 copy (Supplementary Table S1). Therefore, it was considered as being close to the detection limit, and thus negligible. No E. faecalis DNA signal was captured in blank controls. During fungal DNA detection, a positive signal was captured in all baseline controls (n = 25). However, only the ZR method significantly differed (p < 0.05) from every other method, resulting in a higher fungal DNA yield. In addition, fungal DNA was also detected in the blank controls (n = 8/10) with values lower than 10 copies for four methods (QIA, PL, NS, and IHMS) and 100 copies for the ZR method (Supplementary Table S1).

Using all five extraction methods, genomic DNA was successfully isolated from the human feces serving as natural sample type controls. While the performed technical triplicates were reproducible, the DNA yields and quality of DNA extracted differed between methods (Table 1). The highest DNA yields were obtained using the IHMS and ZR methods, while NS, QIA and PL methods resulted in five times lower amounts (p < 0.05). In addition, DNA extracted with the IHMS, PL, or NS methods were not contaminated with RNA, protein, or organic contents, while ZR extraction showed the lowest DNA purity (Table 1).

Next, to verify fungal DNA presence in the human stool DNA extracts, fungal-specific real-time PCR analysis was performed (Supplementary Figure S3). All extractions were performed in triplicates, in the same way as fungal assays in the murine system. For human feces, a similar “yield pattern” was observed to the assay with a high C. albicans community concentration, which indicates a yeast dominance in the examined human stool samples. Blank controls were also quantified, resulting in mean values 3.3 logs lower than the samples (Supplementary Table S1).

In general, all evaluated DNA extraction methods produced a sufficient DNA quantity and quality for subsequent human fecal sample NGS analysis.

The pooled 15 samples (5 method; 3 replicates) returned a total of 384579 16S rDNA gene sequence reads after raw sequence filtration (see section “Materials and Methods”) and chimera removal. These reads (range 71124–83383 per method; range 20938–30270 per replicate) were distributed into 60 bacterial taxa assigned at genus level (Supplementary Table S3). Only 21 unassigned reads (0.005% of the total) passed filter criteria. Forty-nine taxa were detected in all extracted samples, and thus represented a “core microbiome.” Among the other eleven taxa, no unique method-associated taxon was observed (Supplementary Table S3). The blank controls, also sequenced, contained a very low number (n = 1018) of sequencing reads. However, Acinetobacter reads (n = 943) detected in the ZR method’s blank control, were also reflected as increased relative abundance of Acinetobacter in the ZR method’s microbial composition profile (1.2% in ZR vs. 0.003% in others). The complete taxa list detected in the blank controls is shown in Supplementary Table S4.

Bacterial alpha-diversity was similar between the IHMS, NS, PL, and ZR methods, while the QIA method showed the lowest rarefaction curve (Supplementary Figure S4A). Principal coordinate analysis (PcoA) of beta-diversity revealed observable clusters according to the method, using both the unweighted (Figure 3A) and weighted (Figure 3B) UniFrac distance. Moreover, analyses of similarity (ANOSIM), performed for both UniFrac distance matrices, confirmed significant variability among the used extractions (p < 0.001; Figure 3A,B).

To compare taxa abundance differences at genus level between methods, we performed Kruskal–Wallis and GLM tests with adjusted p-values according to the Benjamini–Hochberg procedure. We found that the examined taxa (n = 52; relative abundance > 0.01% of total) varied greatly between extraction methods. We detected significant differences (p < 0.05) in 19/52 taxa applying the Kruskal–Wallis test and 42/52 taxa applying the GLM test (Supplementary Table S5). The ten most abundant bacterial taxa significantly differing between methods are shown in Figure 4.

The same 15 pooled samples (5 method; 3 replicates) were analyzed for their ITS1 target. Chimeric and singleton sequence removal resulted in a total 483224 sequence reads. These reads (range 73630–119305 per method; range 23218–44101 per replicate) were distributed into 72 fungal taxa assigned at genus level. Only 0.06% of the total reads (n = 294) were not assigned. Two taxa (i.e., unidentified Dipodascaceae and Helotiales) with an average relative abundance of 99% were dominant in all extracted samples and the other 70 taxa, (altogether <1% of total abundance) varied between methods. Among them, a majority (n = 47) of the taxa were consistent between replicates (Supplementary Figure S5), while the remaining 23 taxa were occasionally detected in one of the methods’ replicates and were considered environmental/kit contamination. Moreover, most of these rare taxa (36/72) were also present in the blank controls (Supplementary Table S4), although positive signals detected by real-time PCR were several logs beyond the analyzed samples’ signals (Supplementary Table S1).

Contrary to bacteria, the fungal “core microbiome” constituted of only five taxa (Supplementary Table S3) and the most of the recovered taxa (n = 27; out of 47) were uniquely detected by only one method. This variable fraction was probably responsible for the observed alpha-diversity differences (Supplementary Figure S4B).

Fungal beta-diversity analysis showed similar results to the bacterial samples (Figure 3). PCoA based on non-phylogenetic Bray–Curtis dissimilarity revealed distinct clusters according to each method (Figure 3C), and subsequent ANOSIM analysis confirmed lower variability within methods than between them (Figure 3C).

Similar to bacteria, we performed Kruskal–Wallis and GLM tests with adjusted p-values according to the Benjamini–Hochberg procedure to assess the differences between each method’s taxa abundances. Only taxa at genus level with a higher relative abundance than 0.01% of total (n = 16) were tested. We detected significant differences (p < 0.05) in 15/16 taxa by applying a GLM test, although there were no significant differences after applying a Kruskal–Wallis test (Supplementary Table S5).