High-resolution Microbiome Analysis of Host-Rich Samples Using 2bRAD-M Without Host Depletion

doi:10.21203/rs.3.rs-4967509/v1

Download PDF

Article

High-resolution Microbiome Analysis of Host-Rich Samples Using 2bRAD-M Without Host Depletion

https://doi.org/10.21203/rs.3.rs-4967509/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Characterizing human microbiota in samples is crucial for understanding host-microbe interactions but remains challenged by high host DNA contamination. Current depletion strategies are limited by DNA loss and require immediate processing. In this paper, we introduce 2bRAD-M, a reduced metagenomic sequencing method that offers an efficient and integrated analysis of the host-microbe ecosystem, overcoming issues of high host contamination.

Results

Validated on mock samples with > 90% human DNA, 2bRAD-M achieved over 93% in AUPR and L2 similarity metrics. In real saliva samples, it captured diurnal and host-specific microbiome patterns which matched whole metagenome sequencing (WMS) with only 5–10% sequencing effort. In an early childhood caries (ECC) study, 2bRAD-M identified key bacterial indicators and distinguished ECC from healthy subjects with an AUC of 0.92, demonstrating its utility in high host DNA contexts.

Conclusions

2bRAD-M represents a significant technical advancement in microbiome research, particularly in tackling high host DNA contamination. Its ability to provide high-resolution microbial profiles without the need for immediate host depletion processing offers a practical and efficient alternative to existing methods.

Biological sciences/Biological techniques/Sequencing/Next generation sequencing

Biological sciences/Microbiology/Communities/Microbiome

Biological sciences/Microbiology/Communities/Metagenomics

Fig. 1. Schematic illustration of the benchmarking workflow in HoC samples.

(a) HoC mock communities were created by mixing the standard MSA1002 (synthetic community) with 90% and 99% human DNA. Sequencing was performed using 2bRAD-M, 16S, and WMS methods. Obtained profiles were evaluated for microbial identification and abundance estimation using AUPR and L2 similarity, respectively. AUPR: Area under the precision-recall curve, a metric that measures the average performance of a classification model in terms of precision and recall across all possible thresholds, particularly suitable for assessing performance on imbalanced datasets; L2 similarity: a metric that measures the closeness of two data points in space by calculating the inverse of the Euclidean distance between them.

(b) For diurnal saliva samples, we employed 2bRAD-M, 16S, and WMS, with WMS results considered as the gold standard. The acquired feature tables underwent L2 similarity analysis for abundance estimation, followed by diversity analysis. Rarefaction analysis was applied to WMS and 2bRAD-M profiles. Furthermore, the temporal dynamics of oscillating species were examined.

(c) ECC saliva samples were analyzed using 2bRAD-M, short- and long-read 16S sequencing methods. Diagnostic models targeting ECC were developed employing Random Forest algorithms on sequencing datasets. The efficacy of these models was evaluated with ROC curves. Indicative species for the diagnostic model were visualized on a scatter plot, highlighting their importance scores and corresponding AUC values.

Sequencing-based microbiome methodologies have notably expanded our understanding of microbial communities across various anatomical sites of the human body, enabling a progressively nuanced comprehension of host-microbe interactions in human health and diseases. However, compared to the extensively studied gut microbiota, research on microbial communities in high host contamination (HoC) niches, such as saliva, which typically yield > 90% human genome-aligned reads [1], are inadequately understood, implying the pressing need to address the limitations of conventional metagenomic sequencing strategies.

The primary sequencing methods used to establish the taxonomic composition of the microbiome mainly fall into two approaches: (1) targeted sequencing of phylogenetic "marker genes" (such as 16S rRNA sequencing) and (2) whole metagenomic shotgun sequencing (WMS). Although marker-gene-based sequencing is widely employed for characterizing microbial populations in various environmental and host systems, it suffers from several limitations, including ⅰ) limited taxonomic resolution; ⅱ) primer bias affecting bacterial clade representation; ⅲ) PCR bias leading to identification and abundance inaccuracies; and ⅳ) off-target amplification in HoC samples, where the host DNA sequences can be falsely assigned to bacterial taxa [2]. WMS overcomes these hurdles because it analyzes the total DNA content of a sample and does not depend on target-specific primers. However, WMS requires a substantial amount of starting DNA (at least 20–50 ng) and extensive sequencing effort to achieve adequate coverage of microbial genomes in HoC samples [3]. Increasing sequencing depth could improve performance, but is limited by several factors: i) higher costs due to increased human reads; ii) difficulty in detecting the low-abundance yet phenotype-associated microbes in HoC samples; iii) the prohibitive computational costs of ultra-deep sequencing for many labs.

Several strategies have been developed to reduce host DNA in samples, including pre-extraction methods that lyse human cells followed by host DNA degradation and post-extraction methods that separate microbial DNA based on methylation differences [4]. Pre-extraction techniques, such as lyPMA and MEM, necessitate a dual-step approach by exploiting the physical distinctions between prokaryotic microbial cells and eukaryotic host cells or utilizing selective lysis agents like saponin, complemented by enzymatic or chemical DNA degradation methods (e.g., Benzonase nuclease, PMA) [1, 5]. These methods, however, confront challenges in microbial recovery neutrality and preventing DNA loss during processing, hindering their efficacy in low biomass samples [6]. Notably, pre-extraction methods require fresh samples and immediate treatment after sample collection for efficient host depletion, while they may fail to treat frozen samples, a more widely accepted and feasible sample type for large-scale clinical microbiome studies [5]. Regarding post-extraction methodologies, DNA-binding proteins (e.g., methyl-CpG binding protein) and methylation-sensitive endonucleases (e.g., MspJI) are employed to selectively enrich microbial sequences and deplete host DNA respectively [7, 8]. The efficacy of these strategies hinges significantly on the methylation status of the target genomes, which may result in unequal recovery of microbial reads, thus skewing the representation of different microbial lineages [9]; Furthermore, the complete exclusion of host DNA from samples may not be entirely favorable, given that the retention of human DNA yields invaluable insights into the intricate dynamics of microbiome evolution driven by host factors [10].

To address these challenges, we devised a pioneering approach: 2bRAD-M [11]. This reduced genome sequencing strategy employs Type IIB restriction enzymes to digest DNA from organisms into iso-length DNA fragments, typically around 32 bp depending on the enzyme selection, before sequencing [12]. Utilizing iso-length DNA fragments can significantly reduce PCR amplification bias, thereby enhancing the accuracy of taxonomic profiling. This approach is particularly beneficial for samples with a high abundance of host DNA, which typically require numerous PCR cycles to achieve adequate microbial DNA amplification.

We previously assessed the profiling performance of 2bRAD-M on real fecal, skin and FFPE tissue samples [11], yet its performance on saliva samples remained unvalidated. Multiple studies have demonstrated that the salivary microbiome could be used to detect major oral diseases, i.e., periodontal disease and dental caries [13], and was related to a wide range of systemic diseases, including diabetes, cancers and HIV infections [14–16], and even neurodegenerative diseases such as Parkinson's disease [17]. Due to the ease of collection, cost-effectiveness, and non-invasiveness of saliva samples, utilizing salivary microbiome for studying and diagnosing both oral and general diseases has become an increasingly preferred technique [18]. However, the high abundance of host DNA in saliva samples can obscure microbial signals, reducing the sensitivity of microbial detection.

In this study, we first employed mock communities to simulate high host background conditions for assessing the microbial identification and abundance estimation capabilities of 2bRAD-M in comparison to traditional sequencing methods, namely WMS and 16S rRNA sequencing. Subsequently, we evaluated the performance of these three methodologies with real saliva samples, highlighting the high concordance between 2bRAD-M and the current gold standard method, WMS, and identifying several species whose relative abundance fluctuated diurnally. Finally, we validated the ability of 2bRAD-M, alongside 16S short-read and long-read sequencing in distinguishing ECC patients and healthy individuals. Utilizing 2bRAD-M, we identified several key species that significantly enhanced the efficacy of the ECC predictive model.

The practicality of using 2bRAD-M for detecting the HoC microbiome through in vitro simulation

The performance in microbial identification and abundance estimation of 2bRAD-M was benchmarked against 16S rRNA sequencing and WMS sequencing methods. The conventional 16S method targeting the V4-V5 region was applied to generate genus-level taxonomic profiles with QIIME2 [19] based on the Silva database [20]. Additionally, a recently developed 16S rDNA sequencing protocol called the "5R 16S method" was employed to predict taxonomic profiles at the species level [21]. For WMS sequencing, widely recognized bioinformatic tools like MetaPhlAn4 [22] and mOTUs2 [23] were utilized to taxonomically derive profiles down to genus and species levels (Fig. 2 and Table. S1 in Additional File). To critically assess the performance of profiling tools, the precision-recall curve was used as the main indicator for microbial identification [24]. It graphically represents precision and recall scores at distinct abundance thresholds, offering a comprehensive evaluation of classification performance. By comparing the generated taxonomic profiles with the ground truth, the area under the precision-recall curve (AUPR) was calculated to yield a singular metric for consolidating precision and recall scores and effectively evaluating the presence/absence patterns of microbes [25]. Abundance estimation was also assessed based on L2 similarity for a more thorough analysis.

To assess the efficacy of the aforementioned methods in profiling microbial communities amidst substantial interference from host DNA, we employed a mock microbial community comprising evenly mixed DNA from 20 bacterial species (across 18 genera) [26]. This composite served as the basis for preparing stocks spiked with 90% and 99% human DNA respectively. For each type of stock, two technical replicates were generated and subjected to 2bRAD-M, 16S rRNA sequencing, and WMS sequencing methods. Subsequently, the resultant profiles by each methodology, along with the AUPR and L2 similarity metrics, were then compared with the ground truth across the taxonomic ranks, encompassing the genus and species levels.

2bRAD-M exhibited robust microbial identification and abundance estimation capabilities, coupled with commendable technical reproducibility across replicates. In the scenario where host DNA comprises 90% of the sample, the AUPR and L2 similarity scores for 2bRAD-M exhibited exemplary performance, attaining notably high percentages at both the genus and species levels. Conversely, 16S rRNA sequencing yielded relatively lower AUPR and L2 similarity scores across taxonomic levels, whereas WMS demonstrated metric values akin to those of 2bRAD-M. Upon reaching a host DNA proportion of 99% in the sample, the AUPR and L2 similarity scores for 2bRAD-M significantly outstripped those of 16S rRNA sequencing at both the genus and species levels, underscoring its superior capability in delineating the taxonomic composition of the microbial community, even in the presence of high host DNA contamination. Although WMS exhibited AUPR values comparable to those of 2bRAD-M, its L2 similarity values were lower. Notably, under the HoC scenarios, outstanding challenges have been widely observed in experiments using 16S rRNA sequencing, characterized by a pronounced false positive issue (Fig. 2a-b). The elevated false positive rate and diminished accuracy in the abundance estimation correlated with the augmentation of host DNA within mock samples.

Fig. 2. 2bRAD-M revealed highly accurate microbiota composition in HoC mock samples.

In both (a) and (b), the left panel shows the taxonomic profiling results of MSA1002 mixed with 90% human DNA, while results from MSA1002 mixed with 99% human DNA are shown in the right panel. In the stacked bar plots, the colored bars indicate bacteria taxa that can be found in the ground truth, whereas the white bars indicate false positives identified in the profiling results of a given method. AUPR measures the microbial identification accuracy, while L2 similarity measures the similarity between the ground truth and predicted taxonomic profiles.

(a) Genus-level profiling results of the three sequencing methods (16S, 2bRAD-M, WMS) on the same DNA mock community (MSA 1002).

(b) Species-level profiling results of the three metagenomics methods.

Benchmarking the taxonomic profiling performance of 2bRAD-M with diurnal saliva samples

Although the salivary microbiome has been reported to be significantly associated with various oral and systematic diseases, the underlying mechanisms remain unclear. Elucidating the dynamical changes of the salivary microbiota could contribute to a deeper understanding of its relationship with human diseases. To scrutinize the intricate temporal dynamics of salivary microbiota as well as the performance of 2bRAD-M on actual HoC samples, a cohort of eight participants was enlisted to provide saliva samples at four distinct time points: 9 AM, 11 AM, 1 PM, and 5 PM, resulting in a total of 32 samples (hereafter referred to as diurnal saliva samples). Each sample was partitioned into three aliquots for subsequent analysis via WMS, 2BRAD-M, and 16S rRNA sequencing, respectively. (Fig. 1b).

Considering that the metagenomic profilers employed in our study rely on disparate reference databases, the potential for this discrepancy to introduce confounding variables during the evaluation of classification performance across different methodologies becomes apparent. Hence, it is imperative to standardize the database before initiating comparative analyses. 16S data analysis hinges upon the utilization of 16S rRNA databases, such as Greengenes and the SILVA database. In contrast, 2bRAD-M relies on the GTDB; WMS data analysis is conducted using MetaPhlAn4, which provides a script facilitating the conversion of profiles to GTDB-based outcomes. It is noteworthy that the SILVA and GTDB databases do not share a common data universe for microbiome investigations. To facilitate a more rigorous comparison of technical disparities among sequencing methodologies, we undertook the training of a taxonomic classifier utilizing the GTDB within the QIIME2 framework (refer to the code provided in Additional File 2). In light of the potential influence of sequencing depth-related inaccuracies of WMS, the conventional gold standard, with authentic samples, the comparative analysis did not incorporate the AUPR values. Rather, the evaluation centered on the computation of Pearson correlation coefficients and L2 similarity metrics across different methodologies.

After standardization, to comprehensively evaluate the profiling performance of 2bRAD-M with diurnal saliva samples, we compared the taxonomic profiles generated by 2bRAD-M, whole metagenome shotgun sequencing (WMS), and 16S rRNA sequencing at the genus level (Fig. 3a-c). The taxonomic profiles obtained from 2bRAD-M and WMS exhibited high concordance, as evidenced by an average Pearson correlation coefficient (R) of 92.7% and an average L2 similarity score (L2) of 94.7% (average Pearson correlation R = 92.7% and average L2 similarity L2 = 94.7%) (Fig. 3a; Fig. S1 in Additional File 1). Conversely, a notable disparity is observed when comparing these taxonomic profiles to those generated using 16S rRNA sequencing: the average R and L2 are 87.4% and 83.9% between 2bRAD-M and 16S, as well as 82.5% and 85.1% between WMS and 16S (Fig. 3b and 3c; Fig. S2-S3 in Additional File 1).

Fig. 3. The correlation of microbial abundance at the genus or species level between each pair of three sequencing methods (2bRAD-M, 16S, and WMS) in diurnal saliva samples.

Within each group, three out of the total 32 samples were selected to illustrate the correlation between different methodologies. "Average R" and "Average L2" labeled in the figure represent the respective average indices for all 32 samples. The gray box on top of each scatter plot displays the sample ID. Points on the coordinate axis represent unique features identified by the corresponding method, while shared features are indicated by light blue points in the white area of the plot.

(a-c) Comparative analysis of genus-level taxonomic profiling results obtained from 2bRAD-M, WMS and 16S rRNA sequencing. The profiles of 2bRAD-M and WMS are highly concordant.

(d-f) Comparative analysis of species-level taxonomic profiling results of 2bRAD-M and WMS based on overall divergence, sequencing method disparity and computational pipeline variance. 2bRAD-M and WMS are highly similar in sequencing methods, but significant differences were observed in their corresponding computational pipelines. (WMS)PMA: analyzing WMS data using MetaPhlAn4; (WMS)2bRAD-M: analyzing WMS data using the 2bRAD-M computational pipeline; (2b)2bRAD-M: analyzing 2bRAD-M sequencing data with 2bRAD-M computational pipeline.

We next scrutinize the taxonomic profiling efficacy of 2bRAD-M at the species level. Given the inherent limitations in accurately classifying species using 16S data, we opted to exclusively benchmark 2bRAD-M alongside WMS (Fig. 3d-f). Initially, the two approaches were quite dissimilar (average R = 52.0% and average L2 = 69.9%) and both methods identified substantial numbers of unique species absent in the other (Fig. 3d; Fig. S4 in Additional File 1). Although comparing the shared species could enhance the average Pearson correlation to 71.1%, the outcome remained suboptimal (Fig. S5 in Additional File 1). We postulated that the disparity observed in the taxonomic analysis outcomes were attributed to variations in both the sequencing methodologies and the subsequent computational pipelines. The combination of variances from these two factors might have contributed to the high heterogeneity evident in the overall taxonomic profiles.

To assess the disparity in sequencing methodologies, we contrasted the species-level profiling results generated from WMS data processed within the computational framework of 2bRAD-M, denoted as the (WMS)2bRAD-M combination, with those from 2bRAD-M sequencing data processed through the 2bRAD-M pipeline, referred to as the (2b)2bRAD-M combination. Remarkably, the paired results exhibited high degrees of similarity, with an average R = 96.7% and L2 = 96.7% (Fig. 3e, Fig. S6 in Additional File 1). While certain unique species were detected by either data processing combination, the shared species identified by both combinations represented 98.65% and 99.76% of all reads in the 2bRAD-M data and the WMS data, respectively. These results suggest that taxa identified as rare in the two combinations concerned only a minor fraction of the total biomass. Additionally, correlation analysis conducted on the shared species revealed marginally enhanced metrics in the taxonomic profiling outcomes, with both average R and average L2 = 96.9% (Fig. S7 in Additional File 1).

Furthermore, the discrepancy of computational pipelines was compared. We juxtaposed the taxonomic profiles derived from WMS data analyzed with MetaPhlAn4, denoted as (WMS)MPA combination, with those obtained from the (WMS)2bRAD-M combination, revealing a notable discrepancy (Fig. 3f, Fig. S8 in Additional File 1). These findings corroborated a substantial concordance between the sequencing methodologies of 2bRAD-M and WMS, with disparities in taxonomic profiles primarily attributed to different computational pipelines.

2bRAD-M exhibits a significant congruence with WMS when conducting diversity analyses

To substantiate the capability of 2bRAD-M in yielding biologically equivalent results to other methodologies, particularly WMS, we conducted simultaneous analyses of these samples using 2bRAD-M, 16S, and WMS methodologies. All saliva samples (n = 32) were stratified into the gingivitis group (n = 16) and the healthy group (n = 16) according to the oral health status of the participants, with four individuals diagnosed with gingivitis and four without. Both 2bRAD-M and WMS exhibited significant differences between the two groups in alpha diversity, as evidenced by Shannon diversities of 7.36 and 12.28, respectively (p < 0.05). In contrast, 16S data did not reveal any significant difference between the two groups (p = 0.1) (Fig. 4a). We further conducted beta diversity analysis using four distance matrices (Jaccard distance, Bray-Curtis dissimilarity, unweighted UniFrac distance, and weighted UniFrac distance) for the taxonomic profiles generated by each sequencing method. Based on the pseudo-F statistic from the Adonis test, all methods distinguished between gingivitis and health. Remarkably, the statistical power of 2BRAD-M closely approximated that of WMS, and significantly surpassed that of the 16S method (Fig. S9 in Additional File 1). Collectively, both 2bRAD-M and WMS methods demonstrated high efficacy in discerning microbial diversity variances between gingivitis and health states. Conversely, the 16S method exhibited slightly lower sensitivity to inter-group differences due to limited resolution.

To evaluate the technical consistency across methodologies, we performed the Mantel test to assess the similarity of the corresponding distance matrices derived from different sequencing methods (Fig. 4b). We observed a robust correlation (Pearson r = 0.97) between the genus-level Bray-Curtis dissimilarity matrices for 2bRAD-M and WMS. However, the similarity between the 16S and WMS matrices was notably lower (Pearson r = 0.49). At the species level, the Bray-Curtis and weighted UniFrac distance matrices derived from 2bRAD-M and WMS demonstrated strong consistency, with outstanding Pearson correlations of 0.97 and 0.95, respectively. However, the performance of the 16S method was far less satisfactory, as indicated by similarity values of 0.39 and 0.25 to WMS for the respective distance matrices. These suggest 16S yields less accurate abundance estimates for HoC samples, resulting in less correlated distance matrices compared to WMS and 2bRAD-M. Evidently, 2bRAD-M and WMS exhibited high concordance in beta diversity analysis in terms of Bray-Curtis and weighted UniFrac distances. While 2bRAD-M and WMS might have identified unique species, the shared species identified by both methods were predominant and accounted for the majority of the sequencing reads. The identity and relative abundance of dominant species identified by the two methods exhibited an ultra-high degree of similarity. Collectively, 2bRAD-M technology demonstrates both biological and methodological robustness comparable to those of WMS when sequencing HoC samples, represented by saliva specimens in this study.

Fig. 4 2bRAD-M has a high correlation with whole metagenomic sequencing in saliva specimens.

(a) Comparison of alpha diversities (in Shannon index) between health and gingivitis-afflicted samples using 2bRAD-M, WMS and 16S rRNA sequencing. Kruskal-Wallis test was employed to determine the statistical significance of the observed differences.

(b) Mantel test for the correlation between the corresponding distance matrices derived from different methods, with each point denoting the distance between two samples. The distance matrices obtained from 2bRAD-M, whether based on non-phylogenetic distance (Bray-Curtis dissimilarity) or phylogenetic distance (Weighted UniFrac distance), are highly similar to those obtained from WMS.

(c) Rarefaction analysis of 2bRAD-M and WMS datasets. In the graph, each colored line represents a sample. From each sample, three subsamples were randomly sampled at several fixed sequencing depths (if the set subsampling volume exceeded the total data volume of the sample, subsampling would be halted). The average alpha and beta diversity indices for each of the three subsamples at each sequencing depth were calculated and compared to the pre-rarefaction results. The suggested minimal sequencing depth for 2bRAD-M is defined as the depth at which the similarity between the index obtained through resampling 1M and the pre-rarefaction result in 2bRAD-M does not improve by more than 1%. In WMS, it is defined as the depth at which the similarity after resampling 10M does not improve by more than 1%.

2bRAD-M can capture similar biological signals to WMS with just 5-10% of the sequencing effort from host-rich samples

The effectiveness of metagenomic sequencing demonstrates a robust correlation with sequencing depth [27]. To evaluate the minimal sequencing depth for elucidating microbial composition and studying community ecology in 2bRAD-M and WMS sequencing, rarefaction analysis was performed on datasets from both methods across a range of increasing sequencing depths. We sequenced 32 saliva samples with an average read count of 20 million per sample using 2bRAD-M. In comparison, WMS sequencing produced an average of 207 million reads per sample, representing a tenfold increase in read count compared to 2bRAD-M. Subsequently, we bootstrapped these samples down to 0.1, 1, 2, 5, 7.5, 10, 12.5, 15, 17.5, 20, and 22.5 million reads for 2bRAD-M (bootstrap was not conducted if the sequencing depth did not meet the corresponding threshold, totaling 256 samples). WMS data were subsampled to 1, 10, 20, 50, 75, 100, 125, 150, 175, 200, and 225 million sequences per sample (totaling 308 samples). We aimed to determine the sequencing depth necessary for accurate quantification of crucial ecological indicators, including alpha diversity (Shannon entropy) and beta diversity (Bray-Curtis dissimilarity and UniFrac distance), in both 2bRAD-M and WMS datasets. Notably, with a sequencing depth of 5 million reads per sample, corresponding to about 360 megabases (Mb) of sequencing data, the rarefaction curve of 2bRAD-M reaches saturation and attains Shannon Entropy, Bray-Curtis distance, and weighted UniFrac Distance values that are highly congruent with those of the original samples (similarity > 99%). Conversely, WMS data requires 50 million reads to generate accurate diversity estimation comparable with its deepest read depth (Fig. 4c). Given that the reads in the 2bRAD-M dataset are 32bp long, while those in the WMS dataset are 150bp long, our findings indicate that 2bRAD-M needs merely 5-10% of the sequencing data required by WMS to generate the species-resolved taxonomic profile with comparable accuracy, thereby decreasing costs by 50 to 100 times. These cost reductions are especially significant in studies requiring deep sequencing, such as revealing true microbial composition or detecting scarce but ecologically important species in intricate microbial communities with high host DNA. Therefore, 2bRAD-M is particularly well-suited for analyzing HoC samples.

2bRAD-M captured the subtle temporal changes in the saliva microbial communities

Studies have reported the diurnal oscillation patterns of the oral microbiome [28, 29], but further investigation is needed to elucidate the dynamics of microbial communities at the high taxonomic resolution over shorter time frames, such as during the daytime. To benchmark the capacity of 2bRAD-M in capturing subtle oscillations within the microbiome of HoC samples, we partitioned the diurnal saliva dataset into four cohorts based on the sampling time and computed the Bray-Curtis dissimilarity for each participant's salivary microbiome at 11 AM, 1 PM, and 5 PM compared to their baseline (9 AM). Both the analysis of 2bRAD-M and WMS unveiled marked alterations in microbial communities between the 9 AM-5 PM and 9 AM-11 AM time combinations (Fig. 5a, p = 0.02). This observation aligns with previous literature, indicating the detectability of rhythmic microbial community variations over extended periods [28]. In contrast, the results obtained from 16S sequencing showed no significant changes in microbial communities at different time points compared to the baseline (p = 0.20). These suggest that 2bRAD-M exhibits high sensitivity in capturing the temporal oscillations in microbial communities immersed in host DNA.

Using either the 2bRAD-M or the WMS method, we clustered the microbiotas using the principal coordinate analysis (PCoA) plot and further identified which individual species were responsible for the overall diurnal pattern of saliva microbiota. Samples from two times (9 AM and 5 PM) were differentiated and connected by gray lines to indicate if they belong to the same host (Fig. 5b). The taxonomic profiles were projected onto the first principal coordinate (PC1) and showed a gradient-like heterogeneity of microbial community over the daytime period. Changes in PC1 within the host were associated with microbiome alterations in beta diversity from 9 AM to 5 PM. Therefore, PC1 can serve as the primary descriptor for quantitatively measuring the development of microorganisms during daytime oscillation.

Next, we identified individual microbes responsible for changes in the microbiome along PC1 or compositional change using both WMS and 2bRAD-M methods. While we identified 58 PC1-associated bacterial species using 2bRAD-M, 19 were identified using WMS, where eight species were commonly found in both profiles. In the PCoA plot, the 9-AM samples consistently appear on the right side of those at 5-PM ones (Fig. 5b). Six out of 58 species identified by 2bRAD-M negatively correlated with PC1 (i.e., enriched at 5 PM), the remaining 52 species increased in abundance with PC1 (i.e., enriched at 9 AM). Likewise, 11 out of 19 species identified by WMS positively correlated with PC1. The relative abundance of eight shared species between the two methods commonly increases along PC1 (i.e., enriched at 9 AM). We further confirmed the declining trend in the relative abundance of these eight species from 9 AM to 5 PM for each host (Fig. S10 in Additional File 1; Spearman rho > 0.7, FDR q < 0.2). These shared species included Filifactor alocis, Desulfobulbus oralis, Eubacterium_M brachy, Tannerella_forsythia, Treponema_B denticola, Porphyromonas endodontalis, F0058 sp000163695, and Campylobacter_A rectus. In previous literature, the genus Porphyromonas has been shown to exhibit a notable decreasing trend in the saliva of one or multiple individuals [28]. Consistent with these observations, we noted a similar trend for Porphyromonas endodontalis in our study (Fig. 5c). Therefore, 2bRAD-M provides a promising avenue for studying the subtle dynamics of the microbiome in HoC samples and holds potential applications in assessing the relationship between this variation and the host's physiological state.

Fig. 5 The diurnal pattern of saliva microbiota profiled and compared by three sequencing methods.

(a) The box plot shows the Bray-Curtis distance between each pair of samples collected from the same hosts at different time points (i.e., 9 AM VS 5 PM, 9 AM VS 11 AM). Distance values were calculated from the taxonomic profiles produced by 2bRAD-M, WMS and 16S rRNA sequencing. To ensure visual clarity, time point combinations with non-significant differences (e.g., 9 AM vs 1 PM, 11 AM vs 1 PM, 11 AM vs 5 Pm and 1 PM vs 5 PM) are not depicted in the plot. The Wilcoxon signed-rank tests were conducted to indicate the significance between the 9 AM-5 PM combination and the 9 AM-11 AM combination.

(b) Gradient-like changes in microbial beta diversity over two sampling times, 9 AM (red) and 5 PM (blue). PCoA plot was used to visualize the Bray-Curtis dissimilarity matrices of samples collected at the two time points. Gray lines connected samples from the same host.

(c) The relative abundance of Porphyromonas endodontalis increased along the PC1 axis. Each graph includes a color gradient in the bottom left corner indicating relative abundance, ranging from low (dark blue) to high (brown). Rho represents the Spearman correlation coefficient, and p.adj represents the p-value adjusted by Bonferroni correction.

2bRAD-M demonstrates a stronger discriminatory power in classifying early childhood caries (ECC) than conventional metagenomics

Understanding the human microbiome is key to enhancing disease classification, treatment, and prevention by leveraging insights from microbiota profiles. The use of oral microbiota to predict disease states demonstrated in studies utilizing saliva and plaque samples for diseases like Early Childhood Caries (ECC) [30] and gingivitis [31], is a burgeoning area of research. ECC, affecting nearly half of children worldwide, results in considerable social and economic burdens [32]. It causes permanent dental damage, increasing the risk of further decay and tooth loss over a child's life. Therefore, timely diagnosis and prevention of ECC are of utmost clinical significance. Recent studies have shown the effectiveness of using species-level taxonomic profiles generated from 16S short-read sequencing of saliva and plaque to predict ECC [30]. However, models based on saliva alone were less effective (AUC=0.68), potentially due to the dynamic nature of saliva as a sink for other oral microbiome habitats, adding variability to its microbial profile [33]. Furthermore, implementing more refined sequencing techniques like 2bRAD-M may enhance model performance by providing finer species-level or even strain-level resolution.

To evaluate the performance of 2bRAD-M against conventional sequencing methods in diagnosing ECC, we initiated our study by selecting 19 high-quality samples from each of the ECC and healthy control groups in the aforementioned study. These selections were based on the integrity and completeness of data from previous 16S short-read sequencing efforts, which specifically targeted the V3V4 regions. Subsequently, these samples were subjected to both 2bRAD-M and 16S long-read sequencing. For a consistent and fair comparison across all three methods, we used the GTDB release 202 as our reference database.

In our study, we performed diversity analyses and developed machine learning models to compare the effectiveness of three sequencing methods in differentiating ECC from healthy oral microbiomes. Our analyses revealed that both 2bRAD-M and 16S short-read sequencing effectively distinguish between the ECC-affected and healthy groups. Conversely, the long-read sequencing method was not as successful in making this distinction. Specifically, the results of alpha diversity analysis indicate a significant difference in the Shannon Index between the two groups for both 2bRAD-M and 16S short-read sequencing, as confirmed by the Kruskal-Wallis test (p < 0.05) (Fig. 6a). Notably, the Shannon diversity was higher when using 2bRAD-M, suggesting it may detect a broader array of species compared to the 16S method. We then performed beta diversity analysis to examine ecological difference between oral microbiotas of ECC and healthy groups. 2bRAD-M outperforms other methods, evidenced by significant between-sample disparities using both unweighted UniFrac distances and Bray-Curtis dissimilarities (Fig. 6b). 16S short-read data performed moderately well, distinguishing the two groups under Bray-Curtis dissimilarities but not in unweighted UniFrac distances. In contrast, 16S long-read sequencing data exhibits poor performance probably due to its high sequence error rate, failing to discern ECC from health. These findings highlight the superior performance of 2bRAD-M in capturing the complexity of oral microbiomes, emphasizing its potential for enhancing ECC diagnostics.

The efficacy of microbial abundance profiles as diagnostic models for ECC was evaluated using the random forest algorithm. We measured the models’ effectiveness by their area under the receiver operating characteristic curve (AUC), which quantifies their ability to differentiate between ECC and ECC-free cases. The 2bRAD-M sequencing approach yielded the most accurate results (AUC=0.92), surpassing the 16S short-read sequencing (AUC=0.83) and significantly outperforming the 16S long-read sequencing, which showed minimal discrimination (AUC=0.52) (Fig. 6c). Both the 2bRAD-M and 16S short-read sequencing models pinpointed Streptococcus mutans as a key indicator for ECC. Notably, the 2bRAD-M model identified additional biomarkers, underscoring its capability to reveal more nuanced disease state indicators (Fig. 6d). Within the predictive architecture of the 2bRAD-M model, an intriguing phenomenon was observed: confining the predictive input to merely the quartet of most discriminative bacterial species culminated in an optimal model performance. This apex of predictive accuracy was substantiated by a rigorous 10-fold cross-validation approach, underscoring the outstanding discriminative power of these selected taxa (Fig. 6e). Alongside Streptococcus mutans, this cohort includes Propionibacterium acidifaciens, Mitsuokella sp000469545, and Bifidobacterium dentium. Notably, the genus Mitsuokella, in concert with the other three species, has been empirically corroborated as exhibiting a pronounced association with pathogenicity in dental caries, as evidenced by a litany of scholarly inquiries [34-37] (Fig. 6F). These findings demonstrate that the occurrence of dental caries is not solely associated with a single species but rather with multiple species or even a complex sub-community. The high-precision profiling results provided by 2bRAD-M are crucial for identifying key species related to ECC and for constructing highly performant machine learning models.

Fig. 6. Sequencing of ECC saliva samples using 2bRAD-M, long-read 16S sequencing, and short-read 16S sequencing, with diversity analysis and disease classification results based on species-level profiles.

Diversity analysis encompasses (a) the alpha diversity indices and (b) the beta diversity disparities discerned between healthy (blue) and ECC-afflicted (red) samples. The Kruskal-Wallis test was used to evaluate alpha diversity differences, and PERMANOVA was used for beta diversity. The H and F values represent the test statistic, quantitatively indicating the magnitude of the differences.

(c) Classification performance of Random Forest model using species profiles of the three sequencing strategies, assessed by AUC (area under ROC curve).

(d) The scatter plot displays the relationship between the importance scores of different microbial species and their respective AUC values. Red dots indicate an increase in ECC, blue dots signify a decrease, and grey dots represent a neutral status. The size of the dots represents the mean abundance of these species.

(e) Relationship between the numbers of variables included in the 2bRAD-M model and the corresponding predictive performance (the error bar of each dot denotes the standard deviation).

(f) The four most discriminant species in the predictive model were shown using boxplot. The left half of the graph shows the log10-transforme relative abundance of each species in the healthy or ECC group. The utility of each taxon as a potential ECC marker is assessed by AUC in the right half of the graph.

Our study underscored the robust capabilities of 2bRAD-M sequencing in the precise delineation of microbial communities within HoC samples. We implemented 2bRAD-M sequencing, benchmarking its performance against conventional metagenomic methods using both mock and real saliva specimens. The outcomes accentuated the adeptness of 2bRAD-M in the precise microbial identification and abundance estimation within samples containing more than 90% host DNA. In an examination of salivary microbiome oscillations, 2bRAD-M demonstrated not only its ability to yield profiles congruent with those obtained via WMS with reduced sequencing effort but also its adeptness in discerning minor temporal alterations in microbial communities. A cohort analysis concerning ECC, efficaciously distinguished the discriminative microbiota of individuals afflicted with ECC from healthy counterparts, underscoring its utility in identifying critical microbial indicators in HoC samples. The study concludes that 2bRAD-M represents a high-performance, cost-effective sequencing methodology for the analysis of microbial communities in HoC samples, offering profound insights into the intricate dynamics of host-microbe interactions.

Host depletion strategies remain problematic for addressing the challenges in analyzing HoC samples, highlighting the importance of the rational development of 2bRAD-M in this study. Current host depletion strategies generally exploit structural or genomic differences between host and microbial cells, though each approach comes with its own set of limitations [4]. Pre-extraction methods employ selective cell lysis buffers to process host cells in the sample, followed by the degradation of host DNA using DNase or propidium monoazide (PMA). Research has shown that saponins minimally disrupt microbial cells [6], and benzonase nuclease is deemed effective for host DNA removal due to its high specificity and compatibility with various lysis agents [38]. PMA, offering a lower-cost and simpler alternative to traditional enzymatic degradation, eliminates the need for DNA-damaging wash steps [39]. Nonetheless, systematic evaluations indicate that PMA's efficacy is considerably constrained by factors like initial biomass, sample type, and chemical environment, potentially distorting the true profile of microbial communities in complex samples [40]. Moreover, almost all pre-extraction methods, including those modified schemes based on PMA and other newly proposed methods like MEM, require immediate processing of fresh samples to prevent multiple freeze-thaw cycles that significantly damage microbial integrity in clinical samples [5]. These protocols also involve multiple lysis and centrifugation steps, inevitably leading to DNA loss of both host and microbial cells, thus reducing their effectiveness in processing HoC samples [6, 41].

Post-extraction methods aim to separate microbial DNA from host contamination, for instance, by employing selective amplification with DNA probes that specifically bind to microbial/host DNA [42], or by capitalizing on the differences in DNA methylation between eukaryotes and prokaryotes through the use of DNA-binding proteins or methylation-sensitive restriction enzymes to enrich microbial DNA [43]. The application conditions for post-extraction methods are even more stringent. One limitation of using DNA probes is the prerequisite knowledge of the microbial DNA sequences to be enriched. Methods employing hybridization of probes to host DNA are unsuitable for HoC samples due to the efficiency being severely hampered by the vast and complex nature of host DNA. Moreover, only a few microbial methylomes are well-defined, and some prokaryotic methylation patterns are highly similar to those of eukaryotes, severely restricting the application of methylation-based separation methods [44]. In contrast, the 2bRAD-M method offers several advantages over these approaches, including being pretreatment-free and requiring minimal processing while simultaneously producing species-level bacterial, archaeal, and fungal profiles [11]. Unlike host depletion strategies, 2bRAD-M also retains host DNA information in the profile, which has the potential to facilitate studies on host-microbe interaction in clinical samples.

The benchmarking analysis conducted on the mock samples revealed that the 2bRAD-M method achieved results in microbial identification and abundance estimation that are closest to the ground truth. The 16S rRNA sequencing exhibited pronounced false positive rates, erroneously identifying a substantial array of microbial taxa absent from the empirical ground truth. Notably, even the deployment of the advanced 5R 16S technique fell short in providing dependable species-level profiles against a backdrop of high host DNA content. This phenomenon may be attributed to "off-target amplification" when sequencing HoC samples with the 16S method, resulting from marker gene primers binding to the host DNA (e.g., the mitochondrial DNA in humans). It was reported that amplified human or mouse DNA sequences were erroneously clustered and assigned to bacterial taxa [2]. In our study, the relative abundance of E. flexneri significantly increases with the augmented proportion of host DNA, potentially implying a proclivity for amplification of the Escherichia genus.

Regarding ECC saliva samples in our study, the long-read 16S rRNA sequencing failed to differentiate between disease and health. A recent report highlighted that full-length 16S rRNA gene sequencing with PacBio, compared to Illumina short-read sequencing, enhanced the taxonomic resolution in human microbiome samples [45]. Employing data pre-processing and quality filtering methods previously described (i.e., DADA2), we noted that the data had an average passing rate of only 38.4%. This low passing rate may be a critical factor affecting its performance, despite long-read 16S sequencing identifying a greater number of species (191 vs. 146 species) compared to short-read sequencing. Given the emergent nature of the long-read 16S sequencing technology, it is imperative to conduct more comprehensive benchmark for a better understanding and improvement of its technical performance in real-world applications.

We acknowledge the limitations of this study. The validation of 2bRAD-M's performance in HoC samples was primarily based on saliva specimens. Further experimental validation of 2bRAD-M's efficacy in analyzing microbiotas collected from other anatomical sites with high host DNA content, such as the cervicovaginal region and intestinal biopsy tissues, should be conducted. Additionally, since 2bRAD-M is a reference-dependent profiler, its performance is limited by database composition. While we attempted to utilize GTDB to analyze 2bRAD-M data and 16S data types, bioinformatic profiling tools analyzing WMS data typically utilize other reference databases, leading to potential discrepancies when aligning their profiles to GTDB-based outcomes. Given 2bRAD-M can also take WMS as input, we have benchmarked the performance of different sequencing methods using GTDB with 2bRAD-M in mock samples, minimizing technical differences in comparisons.

Sample preparation and sequencing

Mock samples

The mock sample employed for evaluating the efficacy of the 2bRAD-M sequencing approach, in comparison with other sequencing techniques, was a 20-Species Uniform Genomic Mix, designated MSA 1002 (ATCC). The selected mock comprises genomic DNA prepared from fully sequenced, characterized, and authenticated ATCC Genuine Cultures. These cultures were judiciously chosen to represent a spectrum of pertinent phenotypic and genotypic characteristics, including Gram staining reaction, guanine-cytosine (GC) content, genome magnitude, and spore-forming capability. Subsequently, experimental conditions were meticulously devised to mimic the human microbiome samples of high complexity (HoC), formulated to contain either 90% or 99% human DNA, with two replicates generated for each variant.

Diurnal saliva samples

For this experiment, 16 participants were randomly recruited from the Qingdao branch of Shanghai OE Biotech Co., Ltd. All sample collections were performed by the same professional dentist to minimize exogenous contamination in the samples. To this end, participants were instructed not to eat or drink for two hours before the sampling time. The exclusion criteria for participants included the use of antibiotics or medication in the previous month. Samples from eight participants, for whom the clean data accounted for more than 70% of the raw data according to the 2bRAD-M analysis, were selected for further Whole Metagenome Sequencing (WMS) and 16S rRNA sequencing.

ECC saliva samples

Between March and June 2019, an oral epidemiological survey was conducted among 5-year-old children across 13 randomly selected kindergartens in Qingdao, China. This study received ethical approval from the Ethics Committee of Qingdao Municipal Hospital of Stomatology, and informed consent was obtained from all participants or their guardians, detailing oral examination procedures, sample collection, environmental factor gathering, and the use of data in subsequent research and publications. Informed consent was obtained from the participants. All participants gave written consent for their clinical details to be published. A total of 63 children participated in the study, comprising 34 males and 29 females. Based on the decayed, missing, and filled surfaces (dmfs) index, participants were classified into two groups: a high-caries group (n = 32) with severe early childhood caries (dmfs ≥ 6) and a healthy group (n = 31) with no caries (dmfs = 0). Eligibility criteria included no antibiotic use within three months prior to the examination, a dentition consisting entirely of primary teeth with at least 20 teeth present, normal oral mucosa in color and texture, no congenital malformations or systemic diseases, and no orthodontic appliances or accessories worn.

Sample collection

Participants were instructed to tilt their heads downward, slightly open their mouths, and rest their lower lips against the opening of a 15ml sterile centrifuge tube, facilitating the collection of unstimulated saliva. The collection period lasted approximately 5 minutes, during which the participants were asked to directly expectorate saliva into the collection tube. Approximately 5 ml of saliva was collected from each participant per session.

DNA extraction, whole metagenomic shotgun sequencing, 16S rRNA short-read and long-read sequencing

Sample collection tubes were thawed on ice and then centrifuged thoroughly. Adhering to the manufacturer's guidelines with minor adjustments, bacterial genomic DNA from saliva samples was extracted using the Tissue and Blood DNA Isolation kit (Qiagen).

For the mock samples, we attempted sequencing of the 16S rRNA gene's highly variable regions V2V3, V4V5, and five regions from V2 to V8 (termed 5R 16S) and compared the results. The 5R 16S rRNA sequencing approach amplifies regions covering approximately 68% of the full length of the 16S sequence, which significantly enhances the coverage and resolution of bacterial species detection, making it particularly suitable for low-biomass microbial sample analysis.

For the diurnal saliva samples, considering cost issues and the current state of 16S rRNA technology application, we employed the most commonly used sequencing approach by amplifying the V3V4 highly variable regions of the 16S rRNA gene using primer pair (338F/806R). PCR amplification for each sample was performed in triplicate, using 15 µl of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 µM of both forward and reverse primers, and approximately 10 ng of template DNA.

For long-read 16S rRNA sequencing, total DNA was extracted from the ECC saliva samples. Using the full-length 16S primers (i.e., 27F and 1492R), barcode-tagged specific primers were synthesized for PCR amplification. The PCR products were then purified, quantified, and normalized to prepare the sequencing library (SMRT Bell). Once the library was constructed, it underwent quality control checks. Libraries that passed QC were sequenced using the PacBio Sequel system. The PacBio Sequel generates data in BAM format, which is then exported as CCS files using the SMRT Link analysis software. The data for different samples, identified by barcode sequences, are converted into a fastq format for further analysis. All sequences were processed according to the standard QIIME2 (2023.2) pipeline for preprocessing and downstream bioinformatics analysis. The clustering of ASVs was performed using the DADA2 plugin provided by QIIME2.

Both mock community samples (N = 6) and diurnal saliva samples (N = 32) were analyzed with whole metagenomics sequencing (WMS) using the MGI DNBSEQ-T7RS platform. Clean reads were produced and taxonomically characterized using MetaPhlAn4 and Bracken2 with standard settings.

2bRAD-M sequencing

The 2bRAD-M library preparation, with minor modifications from the original method [12], initiated by digesting genomic DNA (1 pg to 200 ng) with 4 U BcgI (NEB) in a 15-µl reaction at 37°C for three hours. 5 µl of the product was assessed on a 1% agarose gel for complete digestion. The ligation reaction, conducted at 16°C overnight, consisted of a 20-µl reaction mixture comprising enzyme-digested DNA, 2 µl of 1× T4 ligase buffer, 1 µl of adaptors, and 800 U of T4 DNA ligase (NEB).

Ligation products underwent PCR amplification in a 40-µl reaction including 7 µl of DNA, 0.1 µM of each Illumina-specific primer (Primer1 and Primer2), 0.3 mM dNTPs, 1× Phusion HF buffer, and 0.4 U Phusion DNA polymerase (NEB). The PCR, using a DNA Engine Tetrad 2 thermal cycler (Bio-Rad), involved 16–28 cycles at 98°C for 5 s, 60°C for 20 s, and 72°C for 10 s, with a final extension at 72°C for 10 min. The desired band (~ 100 bp) was excised from an 8% polyacrylamide gel, and DNA was eluted into nuclease-free water for 12 hours at 4°C. Barcode introduction was performed in a secondary 40-µl PCR containing 50 ng of extracted DNA, 0.2 µM of each barcode-bearing primer (Primer1 and Primer3), 0.6 mM dNTPs, 1× Phusion HF buffer, and 0.8 U Phusion DNA polymerase for seven cycles. The amplified products were purified using the QIAquick PCR Purification Kit (Qiagen) and sequenced on the Illumina HiSeq platform.

Data analysis

Evaluation Metrics: Precision, Recall, AUPR, L2 Similarity, and Pearson Coefficient

To gauge overall performance, we employed precision and recall for the accuracy of microbial identification and L2 distance to assess the precision of abundance estimations in samples. Precision is determined by the ratio of correctly identified species to the overall species identified by the technique, and recall measures the proportion of correctly identified species relative to all species present in the sample. Adjusting the threshold for classifying a species as a true positive varies the precision and recall outcomes. Plotting these metrics, with recall on the x-axis and precision on the y-axis, creates the Precision-Recall (PR) curve. The area under this curve, known as the Area Under the Precision-Recall Curve (AUPR), quantifies the average performance of the model across various decision thresholds in terms of precision and recall. To assess the precision of the abundance profiles, the L2 distance was calculated between the abundance profile of the ground truth and those produced by various metagenomic sequencing methods at specific taxonomic levels (like species or genus). For a clearer performance comparison, L2 similarity was determined by one minus the L2 distance. Furthermore, in the Pearson correlation analysis, when a species was uniquely identified by one method, its abundance was set to zero in the computation of the Pearson correlation coefficient, facilitating a more direct comparison of benchmarking results.

Diagnosis model of ECC samples

A random forest model was trained to identify the disease status based on taxonomic profiles and subsequently evaluated by the area under the Receiver Operating Characteristic (ROC) curve. The model was implemented using default parameters from the R package "randomForest" (ntree = 5000, errortype="oob", n.oob = 10, nfolds = 10). The performance of the microbiome-based model was further assessed through 10-fold cross-validation using species profiles. According to the "rfcv" function in the random forest package, the top four ranked important taxa facilitated a considerably accurate classification of ECC status. A comprehensive method description is provided in the supplementary information.

Ethics approval and consent to participate
All participants expressed their consent to partake in the study and signed the corresponding informed consent forms. This study was approved by the Ethics Committee of Qingdao Municipal Hospital (certificate number: 2022KQYX031), and conducted by the ethical principles established in the Declaration of Helsinki (2008).

Consent for publication
This manuscript includes individual person’s data in case reports, images, and detailed descriptions. Consent for publication was obtained from all participants involved in the study. In the case of children, permission was obtained from their parents or legal guardians. All participants or their guardians have signed informed consent forms that permit data publication. These consent forms are archived and can be provided upon request.

Competing Interests

Jiang Liu, Lisha Zhou, Qianya Niu, and Xiuping Wang are affiliated with Qingdao OE Biotechnology Company. The other authors state that the research was carried out without any commercial or financial relationships that could be interpreted as a potential conflict of interest.

Author Contribution

Y. J and S. H conceived the study. S. H, F. T, and X. W supervised the experimental design. J. L, L. Z, and Q. N conducted the experiments. Y. J, J. L, and Y. Z analyzed the data. Y. J, E. K, and SW. H wrote the manuscript with support from Y. L, Z. X, T. D, and S. H. All authors approved the final manuscript.

Acknowledgements

This study was funded by Health and Medical Research Fund (No. 10212276). The funder played no role in study design, data collection, analysis and interpretation of data, or the writing of this manuscript.

Data Availability

The sequencing data of 2bRAD-M, WMS, next-generation 16S and third-generation sequencing for Mock samples (ATCC MSA 1002 combined with 90% and 99% human DNA) and real saliva samples (diurnal saliva samples and ECC saliva samples) were stored at Sequence Read Archive (SRA): https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1131785. The 2bRAD-M computational pipeline and its associated database files are publicly accessible on GitHub at https://github.com/shihuang047/2bRAD-M. The 2bRAD-M computational pipeline is distributed under the MIT license. All source data and codes for the generation of figures and tables in the manuscript can be downloaded from: https://github.com/jiangys30/HoC_project.git.

Marotz, C.A., et al., Improving saliva shotgun metagenomics by chemical host DNA depletion. Microbiome, 2018. 6(1): p. 42.
Bedarf, J.R., et al., Much ado about nothing? Off-target amplification can lead to false-positive bacterial brain microbiome detection in healthy and Parkinson's disease individuals. Microbiome, 2021. 9(1).
Knight, R., et al., Best practices for analysing microbiomes. Nat Rev Microbiol, 2018. 16(7): p. 410–422.
Shi, Y., et al., Metagenomic Sequencing for Microbial DNA in Human Samples: Emerging Technological Advances. Int J Mol Sci, 2022. 23(4).
Liu, F., et al., Systematic Evaluation of the Viable Microbiome in the Human Oral and Gut Samples with Spike-in Gram+/- Bacteria. mSystems, 2023. 8(2): p. e0073822.
Hasan, M.R., et al., Depletion of Human DNA in Spiked Clinical Specimens for Improvement of Sensitivity of Pathogen Detection by Next-Generation Sequencing. J Clin Microbiol, 2016. 54(4): p. 919–27.
Xu, C., et al., The structural basis for selective binding of non-methylated CpG islands by the CFP1 CXXC domain. Nat Commun, 2011. 2: p. 227.
Oyola, S.O., et al., Efficient depletion of host DNA contamination in malaria clinical sequencing. J Clin Microbiol, 2013. 51(3): p. 745–51.
Liu, G.H., et al., Epigenetic Segregation of Microbial Genomes from Complex Samples Using Restriction Endonucleases HpaII and McrB. Plos One, 2016. 11(1).
Groussin, M., F. Mazel, and E.J. Alm, Co-evolution and Co-speciation of Host-Gut Bacteria Systems. Cell Host Microbe, 2020. 28(1): p. 12–22.
Sun, Z., et al., Species-resolved sequencing of low-biomass or degraded microbiomes using 2bRAD-M. Genome Biol, 2022. 23(1): p. 36.
Wang, S., et al., 2b-RAD: a simple and flexible method for genome-wide genotyping. Nat Methods, 2012. 9(8): p. 808–10.
Sanz, M., et al., Role of microbial biofilms in the maintenance of oral health and in the development of dental caries and periodontal diseases. Consensus report of group 1 of the Joint EFP/ORCA workshop on the boundaries between caries and periodontal disease. J Clin Periodontol, 2017. 44 Suppl 18: p. S5-S11.
Sabharwal, A., et al., The salivary microbiome of diabetic and non-diabetic adults with periodontal disease. J Periodontol, 2019. 90(1): p. 26–34.
Kageyama, S., et al., Characteristics of the Salivary Microbiota in Patients With Various Digestive Tract Cancers. Frontiers in Microbiology, 2019. 10.
Lewy, T., et al., Oral Microbiome in HIV-Infected Women: Shifts in the Abundance of Pathogenic and Beneficial Bacteria Are Associated with Aging, HIV Load, CD4 Count, and Antiretroviral Therapy. AIDS Res Hum Retroviruses, 2019. 35(3): p. 276–286.
Wallen, Z.D., et al., Metagenomics of Parkinson's disease implicates the gut microbiome in multiple disease mechanisms. Nat Commun, 2022. 13(1): p. 6958.
Bhattarai, K.R., H.R. Kim, and H.J. Chae, Compliance with Saliva Collection Protocol in Healthy Volunteers: Strategies for Managing Risk and Errors. Int J Med Sci, 2018. 15(8): p. 823–831.
Bolyen, E., et al., Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol, 2019. 37(8): p. 852–857.
Quast, C., et al., The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res, 2013. 41(Database issue): p. D590-6.
Nejman, D., et al., The human tumor microbiome is composed of tumor type-specific intracellular bacteria. Science, 2020. 368(6494): p. 973-+.
Blanco-Míguez, A., et al., Extending and improving metagenomic taxonomic profiling with uncharacterized species using MetaPhlAn 4. Nature Biotechnology, 2023.
Milanese, A., et al., Microbial abundance, activity and population genomic profiling with mOTUs2. Nat Commun, 2019. 10(1): p. 1014.
Ye, S.H., et al., Benchmarking Metagenomics Tools for Taxonomic Classification. Cell, 2019. 178(4): p. 779–794.
Saito, T. and M. Rehmsmeier, The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets. PloS one, 2015. 10(3): p. e0118432.
Amos, G.C.A., et al., Developing standards for the microbiome field. Microbiome, 2020. 8(1): p. 98.
Zhu, Q.Y., et al., Phylogeny-Aware Analysis of Metagenome Community Ecology Based on Matched Reference Genomes while Bypassing Taxonomy. Msystems, 2022. 7(2).
Hu, Y., et al., Diurnal and eating-associated microbial patterns revealed via high-frequency saliva sampling. Genome Res, 2022. 32(6): p. 1112–1123.
Takayasu, L., et al., Circadian oscillations of microbial and functional composition in the human salivary microbiome. DNA Res, 2017. 24(3): p. 261–270.
Teng, F., et al., Prediction of Early Childhood Caries via Spatial-Temporal Variations of Oral Microbiota. Cell Host Microbe, 2015. 18(3): p. 296–306.
Huang, S., et al., Predictive modeling of gingivitis severity and susceptibility via oral microbiota. Isme Journal, 2014. 8(9): p. 1768–1780.
Casamassimo, P.S., et al., Beyond the dmft: the human and economic cost of early childhood caries. J Am Dent Assoc, 2009. 140(6): p. 650–7.
Parisotto, T.M., et al., Early Childhood Caries and Mutans Streptococci: A Systematic Review. Oral Health & Preventive Dentistry, 2010. 8(1): p. 59–70.
Obata, J., et al., Pathogenic mechanisms of cariogenic Propionibacterium acidifaciens. Arch Oral Biol, 2019. 105: p. 46–51.
Skelly, E., et al., Response of Salivary Microbiota to Caries Preventive Treatment in Aboriginal and Torres Strait Islander Children. J Oral Microbiol, 2020. 12(1): p. 1830623.
Becker, M.R., et al., Molecular analysis of bacterial species associated with childhood caries. J Clin Microbiol, 2002. 40(3): p. 1001–9.
Okada, M., et al., PCR detection of Streptococcus mutans and S. sobrinus in dental plaque samples from Japanese pre-school children. J Med Microbiol, 2002. 51(5): p. 443–447.
Nelson, M.T., et al., Human and Extracellular DNA Depletion for Metagenomic Analysis of Complex Clinical Infection Samples Yields Optimized Viable Microbiome Profiles. Cell Rep, 2019. 26(8): p. 2227–2240 e5.
Fittipaldi, M., A. Nocker, and F. Codony, Progress in understanding preferential detection of live cells using viability dyes in combination with DNA amplification. J Microbiol Methods, 2012. 91(2): p. 276–89.
Wang, Y., et al., Whole microbial community viability is not quantitatively reflected by propidium monoazide sequencing approach. Microbiome, 2021. 9(1): p. 17.
Ferretti, P., et al., Experimental metagenomics and ribosomal profiling of the human skin microbiome. Exp Dermatol, 2017. 26(3): p. 211–219.
Smith, M., et al., An In-Solution Hybridisation Method for the Isolation of Pathogen DNA from Human DNA-rich Clinical Samples for Analysis by NGS. Open Genomics J, 2012. 5.
Barnes, H.E., et al., Selective microbial genomic DNA isolation using restriction endonucleases. PLoS One, 2014. 9(10): p. e109061.
Bewick, A.J., et al., Diversity of cytosine methylation across the fungal tree of life. Nat Ecol Evol, 2019. 3(3): p. 479–490.
Buetas, E., et al., Full-length 16S rRNA gene sequencing by PacBio improves taxonomic resolution in human microbiome samples. BMC Genomics, 2024. 25(1): p. 310.

Supplementary S1 is not available with this version.

Competing interest reported. Jiang Liu, Lisha Zhou, Qianya Niu, and Xiuping Wang are affiliated with Qingdao OE Biotechnology Company. The other authors state that the research was carried out without any commercial or financial relationships that could be interpreted as a potential conflict of interest.

Supplementaryfigures.zip

Download PDF

Reviews received at journal
14 Nov, 2024
Reviewers agreed at journal
22 Oct, 2024
Reviewers agreed at journal
08 Sep, 2024
Reviewers invited by journal
08 Sep, 2024
Editor assigned by journal
02 Sep, 2024
Submission checks completed at journal
26 Aug, 2024
First submitted to journal
24 Aug, 2024

You are reading this latest preprint version

High-resolution Microbiome Analysis of Host-Rich Samples Using 2bRAD-M Without Host Depletion

Status:

Version 1

Abstract

Background

Results

Conclusions

Figures

Background

Results

The practicality of using 2bRAD-M for detecting the HoC microbiome through in vitro simulation

Benchmarking the taxonomic profiling performance of 2bRAD-M with diurnal saliva samples

2bRAD-M exhibits a significant congruence with WMS when conducting diversity analyses

2bRAD-M can capture similar biological signals to WMS with just 5-10% of the sequencing effort from host-rich samples

2bRAD-M captured the subtle temporal changes in the saliva microbial communities

2bRAD-M demonstrates a stronger discriminatory power in classifying early childhood caries (ECC) than conventional metagenomics

Discussion

Methods

Sample preparation and sequencing

Mock samples

Diurnal saliva samples

ECC saliva samples

Sample collection

DNA extraction, whole metagenomic shotgun sequencing, 16S rRNA short-read and long-read sequencing

2bRAD-M sequencing

Data analysis

Evaluation Metrics: Precision, Recall, AUPR, L2 Similarity, and Pearson Coefficient

Diagnosis model of ECC samples

Declarations

Competing Interests

Author Contribution

Acknowledgements

Data Availability

References

Supplementary Table

Additional Declarations

Supplementary Files

Status:

Version 1