Metagenome sequencing-based strain-level and functional characterization of supragingival microbiome associated with dental caries in children

Abstract

Studies of the microbiome associated with dental caries have largely relied on 16S rRNA sequence analysis, which is associated with PCR biases, low taxonomic resolution, and inability to accurately study functions. Here, we employed whole metagenome shotgun sequencing, coupled with high-resolution analysis algorithm, to analyze supragingival microbiomes from 30 children with or without dental caries. A total of 726 bacterial strains belonging to 406 species, in addition to 34 bacteriophages were identified. A core bacteriome was identified at the species and strain levels. Species of Prevotella, Veillonella, as yet unnamed Actinomyces, and Atopobium showed strongest association with caries; Streptococcus sp. AS14 and Leptotrichia sp. Oral taxon 225, among others, were overabundant in caries-free. For several species, the association was strain-specific. Furthermore, for some species, e.g. Streptococcus mitis and Streptococcus sanguinis, sister strains showed differential associations. Noteworthy, associations were also identified for phages: Streptococcus phage M102 with caries and Haemophilus phage HP1 with caries-free. Functionally, potentially relevant features were identified including urate, vitamin K2, and polyamine biosynthesis in association with caries; and three deiminases and lactate dehydrogenase with health. The results demonstrate new associations between the microbiome and dental caries at the strain and functional levels that need further investigation.

Keywords: Bacteria; High-throughput nucleotide sequencing; dental caries; metagenome; microbiota.

Figure 1.

The microbiome profile of supragingival plaque. DNA extracted from supragingival samples was shotgun-sequenced on an Ion S5XL sequencer. The generated reads (200 bp) where classified to the strain level using the CosmosID metagenomic pipeline that employs clade-specific k-mers derived from a comprehensive reference genome database (see text for details). (a) bacterial phyla, (b) bacterial species and (c) bacterial strains identified at average relative abundance ≥ 1%. (d) Relative abundance of major groups of bacterial phages detected.

Figure 2.

Species richness and alpha diversity. The taxonomic profiles obtained with compositional analysis of the sequences were used to calculate alpha diversity indices employing standard QIIME scripts. Significance of differences between groups were then sought with the Kruskal-Wallis test. The figure shows a comparison of observed species richness, Chao1 (expected richness), Shannon and Simpson indices between the study groups. * P ≤ 0.05.

Figure 3.

β-Diversity analysis. Inter-sample distances were calculated based on abundance Jaccard index employing standard QIIME scripts. The samples were then clustered using Principle Coordinate Analysis (PCoA) to visualize distances at the species (a) and strain (b) levels. (c) Quantitative representation with boxplots of distances between the samples based on abundance Jaccard index at the species and strain levels. * P ≤ 0.001, Wilcoxon rank sum test.

Figure 4.

Differentially abundant taxa. (a) species and (b) strains that showed significant differences in relative abundance between the caries and caries-free groups, as identified by linear discriminant analysis (LDA) effect size analysis (LEfSe), with ‘caries’ and ‘caries-free’ as classes and ‘early-caries’ and ‘advanced caries’ as subclasses. Note: many unnamed species, e.g. Streptococcus sp_AS14, are represented in the database by only one strain, and in this case the species and strain names are identical.

Figure 5.

Differential association with health and disease at the strain level. Examples of strains belonging to same species but showing different patterns of association with the study groups. ¶Identified as differentially abundant by LEfSe (see Figure 4). For the right panel, significance of differences was sought with the Kruskal-Wallis test. NS: no significant difference; * P ≤ 0.05.

Figure 6.

Differentially enriched microbial features. Quality-filtered sequences were functionally analyzed with HUMAnN2, which involves nucleotide search by mapping to pangenomes in ChocoPhlAn and translated search against Uniref90. Generated gene lists are collapsed into protein families, class families and pathways using Pfam and Metacyc databases. Linear discriminant analysis (LDA) effect size analysis (LEfSe) was then used to identify pathways (a) and enzyme classes (b) with significant differences in relative abundance between the caries and caries-free groups.

Figure 7.

Species contributions to key functions. For each identified enzyme class or pathway, HUMAnN2 calculates species contribution based on results from nucleotide search in which reads are mapped to annotated pangenomes of known species. In the figure are species contributions to six features of potential relevance to dental caries.