Systematic Bias Introduced by Genomic DNA Template Dilution in 16S rRNA Gene-Targeted Microbiota Profiling in Human Stool Homogenates

The genomic DNA input for stool samples utilized for microbiome composition has not been determined. In this study, we determined the reliable threshold level under which conclusions drawn from the data may be compromised. We also determined the type of microbial bias introduced by less-than-ideal genomic input.

has a key role in health and in disease (1,2). The study of microbially diverse populations is routinely performed with 16S rRNA gene library construction and sequencing (3,4), a protocol that requires an adequate microbial DNA concentration for accurate representation of the microbial community.
The effects of multiple variables on gut microbiome representation have been studied, and biases can be introduced through DNA extraction (4)(5)(6), PCR (7)(8)(9)(10), and sequencing (11). For instance, PCR parameters such as the choice of primers and polymerase and the number of reaction cycles affect the richness (total number of species) and evenness (relative contribution of a given species to the total number of species) of the resulting sampling of the bacterial communities (4,(7)(8)(9)(10). The impact of the template genomic DNA (gDNA) concentration in the community representation has been studied in sediment samples, where Wu et al. (8) reported that a dilution of from 20-fold to 200-fold reduced the determined taxon richness without having a significant impact on taxon evenness. In contrast, Biesbroek et al. (12) showed that diluting saliva DNA to below Ͻ1 pg/l affected relative abundances and species representations, with an increased relative abundance of Proteobacteria and a significant decrease in the abundance of Bacteroidetes compared with identical samples with DNA concentrations above this threshold.
Because of the high concentration of native microbes in stool samples, the gDNA concentration typically is not a parameter of concern for evaluating gut microbiomes. However, this scenario can be affected by disease conditions such as diarrhea, recent antibiotic use, or chemotherapy drug use. Regardless, the initial gDNA concentrations, which often are not reported in publications, can potentially introduce relevant bias in study findings. In addition, dilution of stool gDNA templates has become a commonplace method to achieve optimal PCR efficiency and to dilute reaction inhibitors and host tissue DNA. However, the effect of the starting gDNA concentration in stool samples has not been studied. Similarly to the studies that established gDNA thresholds for accurate community representation in sediment and saliva samples, we sought to determine the level of input gDNA concentrations at which the microbial community becomes misrepresented in stool samples. We also studied the impact of different extraction methods and PCR conditions on the resulting microbial community profile.

RESULTS
Effects of technical variations on the overall microbiota profile. The total number of sequence reads was 6,888,783, with a sample range of 11,328 to 216,333. The median number of sequence reads was 70,601, and the median absolute deviation was 57,817. Among the tested variables, gDNA concentration and 16S rRNA gene yield had the strongest effects on the microbial community representation, explaining 22.4% to 38.1% and 20.7% to 42.8% of the overall microbiota variation, respectively, as summarized in Table 1. In contrast, the tested extraction methods (manual PowerSoil [Mo Bio Laboratories, Inc.] and automated Chemagic) and PCR volumes (25 and 50 l) did not significantly affect microbial composition. PCR annealing temperatures (70°C annealing temperature and 69°C touchdown annealing temperature) significantly affected the microbiota representation but explained less than 3.7% of the community variation. Analysis of correlations between the 2 major predictors of the microbial profile (gDNA and 16S rRNA gene concentration) showed a high correlation (Spearman r ϭ 0.83) ( Fig. 1), suggesting that the effect of dilution on the microbial profile is likely mediated by the 16S rRNA gene yield.
Overall diversity in microbiota composition associated with different levels of gDNA input. The levels of diversity in microbial composition observed with different levels of gDNA input were compared by assessing within-sample diversity (␣-diversity) and between-sample diversity (␤-diversity). The ␣-diversity, which was investigated by determining the number of observed operational taxonomic units (OTUs) (measure of sample richness) and the Shannon diversity index values (measure of sample richness and evenness), showed that low gDNA input was associated with lower overall diversity (Fig. 2). In particular, compared with a gDNA input of 1 ng/l, a gDNA input of 8 ϫ 10 Ϫ3 ng/l or lower was associated with a significant decrease in species richness (P Ͻ 0.001), whereas a gDNA input of 1.6 ϫ 10 Ϫ3 ng/l or lower was associated with decreased overall diversity, as measured by the Shannon index (P Ͻ 0.001) ( Table 2). Consistent with the ␣-diversity, the ␤-diversity analyses (Bray-Curtis, unweighted Uni-Frac, and weighted UniFrac distances) showed microbiota dissimilarity for different levels of gDNA input. Although principal-coordinate analyses based on all three ␤-diversity measures showed that a decrease in gDNA input was associated with a clear change in the representation of the microbiota (Fig. 3), the fact that the separation was more evident in the unweighted UniFrac analysis indicates that the major difference in microbiota occurs in the presence of and with an abundance of rare species (13) ( Table 3). On the basis of the unweighted UniFrac distance data, a gDNA input of 4 ϫ 10 Ϫ2 ng/l or lower resulted in a significant difference in the ␤-diversity profile compared with an input of 1 ng/l. This difference was no longer significant if the gDNA input was at least 2 ϫ 10 Ϫ1 ng/l.
Specific microbial lineage bias associated with dilution. Differential abundance analysis was performed to identify differentially abundant taxa at different levels of gDNA input. In particular, using 1.6 ϫ 10 Ϫ3 ng/l of gDNA input as a cutoff, the analysis identified 55 differentially abundant genera (false-discovery-rate [q] value, Ͻ0.01) (Fig. 4a). The magnitude of difference for each of the taxa that were differentially abundant between samples with high gDNA input (Ͼ1.6 ϫ 10 Ϫ3 ng/l) and low gDNA input (Յ1.6 ϫ 10 Ϫ3 ng/l) is indicated with the Ϫlog P value (Fig. 4b). Interestingly, the abundances of several genera of the phylum Proteobacteria, including Pseudomonas, Delftia, Phenylobacterium, and Stenotrophomonas, were highly overestimated in samples with low gDNA input. In contrast, numerous genera from the phylum Firmicutes (Christensenella, Faecalibacterium, Clostridium, Megamonas, Eubacterium) were underrepresented with low gDNA input (Fig. 4b). Moreover, the abundances of taxa of the low-gDNA-input group showed high variability (Fig. 4c). These differentially abundant genera could effectively separate the 2 groups, as revealed by hierarchical clustering based on the abundance profiles (Fig. 4d).

DISCUSSION
In this study, we showed that template dilution and 16S rRNA gene yield had a significant effect on the representation of the microbial composition of stool samples. In particular, gDNA input at or below 8 ϫ 10 Ϫ3 ng/l was associated with significantly lower within-sample diversity (␣-diversity) and divergences in between-sample diversity (␤-diversity). Moreover, template dilution was associated with a change in the microbial lineage, and a gDNA input level of 1.6 ϫ 10 Ϫ3 ng/l or less was associated with overrepresentation of Proteobacteria and underrepresentation of Firmicutes. In particular, of the 22 taxa enriched by dilution, 10, namely, Pseudomonas, Delftia, Stenotrophomonas, Flavobacterium, Pedobacter, Chryseobacterium, Enterobacter, Ralstonia, Pelomonas, and Streptococcus, have been previously identified as common contaminants (14). This indicates that some of the bias introduced may relate to the increase in the contamination signal caused by the use of diluted template. However, the remaining results representing taxa overrepresented by dilution effects are likely unrelated to contaminating sources. A low gDNA concentration also increases the variability of the relative abundances and hence results in a larger measurement error. The tested DNA extraction methods (manual PowerSoil [Mo Bio, Inc.] and automated Chemagic), PCR annealing conditions (70°C and 69°C touchdown), and PCR volumes (25 and 50 l) had no significant impact on the microbiome composition.
Previous studies evaluated the effect of template dilution on the microbial community representation of sediment samples (8) and saliva (12), but their results are not directly comparable to ours because of the high microbial abundance that characterizes stool samples. However, consistent with our findings, dilution of saliva was associated with an increase in the relative abundance of Proteobacteria (12).
To determine the impact of gDNA concentration on the accurate representation of the bacterial community and detection of rare species, which can be key components of a healthy or diseased community, we analyzed the ␣-diversity. Of the 2 ␣-diversity metrics used (number of observed OTU and Shannon index), the effect of reduced gDNA concentration on the number of observed OTUs was the more prominent, suggesting that the strongest effect is on species richness (detection of rare species) rather than on species evenness. Our results suggest that use of a gDNA concentration at or below 8 ϫ 10 Ϫ3 ng/l would result in a significant underestimation of the number of species present and hence would influence the conclusions of the study.
To determine the potential bias introduced by gDNA concentration with respect to the variation between samples, we analyzed the ␤-diversity. Bias introduced at this level can artificially inflate sample divergence, leading to artificial groupings or diminishing the capacity to differentiate between truly divergent cohorts. We evaluated ␤-diversity with 3 metrics (Bray-Curtis, unweighted UniFrac, and weighted UniFrac distances) and showed that the effect of dilution was strongest on the unweighted UniFrac data, which capture the difference in rare and less-abundant species. This finding was consistent with what we observed with ␣-diversity, indicating that the loss in detection of rare species affects not only the within-sample diversity metrics but also the between-sample comparisons, potentially skewing results toward erroneous con- clusions. The microbiota variability caused by different gDNA concentrations is therefore a critical criterion for consideration. If DNA concentrations cannot be normalized among samples, particularly between cases and controls, statistical confounders may be created. More concerning, the gDNA concentration could confound the association of interest if it is also associated with the variable of interest. Unlike the findings of other studies (4,5,8), the effects of DNA extraction and PCR conditions on microbially diverse populations in our study were very limited, with DNA extraction not being a significant parameter and PCR conditions explaining less than 4% of the microbial variation. Both of these conclusions apply only to the variations that we have tested in both of these factors.
The present study had numerous strengths. To our knowledge, this was the first study to evaluate the effect of template dilution on the representation of microbial communities from stool samples. Because dilution has different effects (depending on microbial abundance) (12), these results will allow an improvement in the interpretation of the variation in the microbial community of stool samples. Moreover, because all fecal sample collection methods have been demonstrated to be reproducible, stable, and accurate (15), our results are generalizable to all studies evaluating fecal samples, independently of the collection methods. Another of the strengths of the study is that we used samples from the same homogenate, which eliminated the subject-to-subject variability that occurs in a population-based study and that would markedly affect variability (16). The study was thus powered to detect even very small effects of the technical variables evaluated. We note that the conclusions of this study are limited to the gDNA concentrations and processing conditions tested.
Conclusion. We performed a systematic and comprehensive assessment of the impact of various experimental parameters on the resulting 16S rRNA gene library and showed that microbial community representation can be affected by gDNA concentration and 16S rRNA gene yield. On the basis of these results, we recommend a minimum concentration above 4 ϫ 10 Ϫ2 ng/l gDNA input, and, ideally, Ͼ2 ϫ 10 Ϫ1 ng/l, to achieve an unbiased representation of the microbial community. The use of gDNA input levels of Յ1.6 ϫ 10 Ϫ3 ng/l is not recommended because such levels would introduce taxonomic biases that would result in a definitive misrepresentation of the microbiome.

MATERIALS AND METHODS
Sample collection and homogenization. Stool samples were prospectively collected from healthy volunteers (age, Ն18 years) at Mayo Clinic (Rochester, MN) from 4 September through 7 November 2014. Written informed consent for each subject was obtained following a protocol approved by the Mayo Clinic Institutional Review Board (protocol 13-005217). Participants were excluded from the study if they had received antibiotic or probiotic therapy within the 2 weeks before the collection date, had received chemotherapy, or had a history of pelvic radiation. Details of sample collection are reported elsewhere (15).
For the present study, we used freshly collected stool from 8 patients; samples were combined with sterile phosphate-buffered saline and blended in a Stomacher 400 paddle blender (Seward Laboratory Systems Inc.). Each 150-l aliquot of blended stool was placed in a 2-ml screw-cap tube (Starstedt; catalog no. 72.694.217) and stored at Ϫ80°C. The contents of one PowerBead tube (Mo Bio Laboratories, Inc.; catalog no. 12888-50-PBT) and 60 l of solution C1 (Mo Bio Laboratories, Inc.; catalog no. 12888-50-1) were added to each aliquot. Aliquots were homogenized for 40 s with a FastPrep-24 instrument (MP Biomedicals LLC) at a speed setting of 6.0. Stool homogenates were clarified by centrifugation, pooled (2-ml total volume), and diluted to 3.4 ml with beading solution (Mo Bio Laboratories, Inc.; catalog no. 12888-50-BS). Four 800-l aliquots were prepared from each homogenate pool and stored at Ϫ80°C before extraction testing. The study design is shown in Fig. 5. gDNA Concentration Bias in 16S Analysis permutational multivariate analysis of variance (PERMANOVA; "Adonis" function in the R "vegan" package). Significance was assessed with 1,000 permutations. A distance-based R 2 value was used to quantify the percentage of variability of the microbiota that was explained by various technical factors (28). Ordination plots based on Bray-Curtis, unweighted UniFrac, and weighted UniFrac distances were generated with principal-coordinate analysis as implemented in R ("cmdscale" function in the R "stats" package). To identify bacterial genera associated with gDNA concentration, we conducted differential abundance analysis at the genus level by using the Wilcoxon rank sum test. We filtered out genera with a prevalence of less than 10%. We normalized the count data into relative abundances (proportions) by dividing by the total read count. Taxa with a maximum proportion of less than 0.2% were excluded from testing to reduce the number of the tests. False-discovery-rate control (B-H procedure; "padjust" in standard R packages) was used to correct for multiple testing, and false-discovery-rate-adjusted P values (q values) of less than 0.01 were considered significant.
Data availability. The sequencing data sets are available in the National Center for Biotechnology Information repository (study ID SRP098583).