Metagenomics of the modern and historical human oral microbiome with phylogenetic studies on Streptococcus mutans and Streptococcus sobrinus

Abstract

We have recently developed bioinformatic tools to accurately assign metagenomic sequence reads to microbial taxa: SPARSE for probabilistic, taxonomic classification of sequence reads; EToKi for assembling and polishing genomes from short-read sequences; and GrapeTree, a graphic visualizer of genetic distances between large numbers of genomes. Together, these methods support comparative analyses of genomes from ancient skeletons and modern humans. Here, we illustrate these capabilities with 784 samples from historical dental calculus, modern saliva and modern dental plaque. The analyses revealed 1591 microbial species within the oral microbiome. We anticipated that the oral complexes of Socransky et al., which were defined in 1998, would predominate among taxa whose frequencies differed by source. However, although some species discriminated between sources, we could not confirm the existence of the complexes. The results also illustrate further functionality of our pipelines with two species that are associated with dental caries, Streptococcus mutans and Streptococcus sobrinus. They were rare in historical dental calculus but common in modern plaque, and even more common in saliva. Reconstructed draft genomes of these two species from metagenomic samples in which they were abundant were combined with modern public genomes to provide a detailed overview of their core genomic diversity. This article is part of the theme issue 'Insights into health and disease from ancient biomolecules'.

Keywords: ancient DNA; dental calculus; dental plaque; genomic reconstruction; metagenomes; saliva.

Figure 1.

Source specificity of the percentage of species composition in 784 oral metagenomes according to SPARSE. (a) X-Y plot of the first three components from a UMAP (Uniform Manifold Approximation and Projection) [56] dimensional reduction of taxon abundances. (b) Neighbour-joining (FastMe2; [57]) hierarchical clustering based on the Euclidean distances between pairs of metagenomes. Euclidean p-distances were calculated between each pair as the square root of the sum of the squared pairwise differences in the percentage of reads assigned by SPARSE to each microbial taxon. Nodes whose cluster location was inconsistent with the UMAP clustering in (b) are highlighted with black perimeters. Tree visualization: GrapeTree [38].

Figure 2.

Average percentage abundance (left axis) of bacterial species by source for the 40 most discriminating species according to Support Vector Machine analysis. The relative abundances for each of the three sources are indicated by mini-histograms for each species; error bars indicate standard deviations. Species are sorted in descending order by predominant source and then by SVM weight (squared coefficient) in the optimal model. Species belonging to oral complexes are indicated by oral complex-specific shapes and colours. Key legend: source colours used in the mini-histograms and symbol for SVM weight. Asterisk, species designations assigned by RefSeq to single genomes which have not (yet) been confirmed by taxonomists. S. mitis is separated into multiple ANI95% clusters, two of which (s8897; s126097 (electronic supplementary material, table S3)) are among the predominant taxa associated with saliva.

Figure 3.

Average percentage abundances in 784 metagenomes by oral source (key legend) of 28 species from six oral complexes described by Socransky et al. [26]. The oral sources are indicated by three mini-histogram bars for each species. Species are ordered from left to right by oral complex, whose colour designation is indicated at the top. Within each oral complex, the species order is by decreasing total abundance.

Figure 4.

Neighbour-joining (FastMe2; [57]) hierarchical clustering based on the Euclidean distances between pairs of 245 microbial species whose percentage abundance was greater than 2% in at least one metagenome. Members of the six oral complexes [26] are highlighted by coloured species names, whose colours indicate their oral complex membership. These species do not cluster by oral complex, but by other unnamed groupings, four of which are highlighted in grey. An expanded version of the same tree including all species labels is available in electronic supplementary material, figure S2. Branch length distance scale bar is next to the distance of 0.1.

Figure 5.

Numbers of microbial taxa by source. (a) Rarefaction curves of numbers of species by source, with 95% confidence estimates (shadow). Inset data indicate median numbers of species per sample by source, as well as the total numbers for all sources. Rarefactions were performed with the program script called SPARSE_curve.py, using 1000 randomized permutations of the order of samples. (b) Binned histograms of number of species by percentage of samples. The data for this plot was also calculated with SPARSE_curve.py. (c) Venn diagram of overlapping presence of taxa (≥0.0001% abundance) for the three oral sources.

Figure 6.

Reconstruction of pseudo-MAGs (metagenomic assembled genomes) of S. mutans and S. sobrinus from oral metagenomes. (a and c) Numbers of oral samples by source binned by the percentage of reads specific to S. mutans (a) and S. sobrinus (c). (b and d) Numbers of oral samples by source with an average coverage of at least 1x. The data are binned by the predicted read coverage against a reference genome of S. mutans (UA159; (b)) and S. sobrinus (NCTC12279; (d)). (e and f) Read coverage (dots; left) and percentage of the reference genome that was unmasked (≥3 reads; ≥70% consistency) (histogram; right) in S. mutans (e) and S. sobrinus (f). Ordered by decreasing percentage unmasking.

Figure 7.

Maximum-likelihood phylogenies of S. mutans and S. sobrinus genomes. (a) A RaxML [37] tree of 226 genomes of S. mutans (RefSeq: 195; pseudo-MAGs: 31) plus one genome of S. troglodytae as an out-group. The tree was based on 181 321 non-repetitive SNPs in 1.73 Mb. (b) A RaxML tree of 61 genomes of S. sobrinus (RefSeq: 46; pseudo-MAGs: 15) plus six S. downei genomes as an out-group. The tree was based on 160 863 non-repetitive SNPs in 1.13 Mb. Pseudo-MAGs are highlighted by thick black perimeters. Visualization with GrapeTree [38]. Branches with a genetic distance of greater than 0.1 were shortened for clarity and are shown as dashed lines. Legend: numbers of strains by country of origin for both trees.