The evolution and changing ecology of the African hominid oral microbiome.

The oral microbiome plays key roles in human biology, health, and disease, but little is known about the global diversity, variation, or evolution of this microbial community. To better understand the evolution and changing ecology of the human oral microbiome, we analyzed 124 dental biofilm metagenomes from humans, including Neanderthals and Late Pleistocene to present-day modern humans, chimpanzees, and gorillas, as well as New World howler monkeys for comparison. We find that a core microbiome of primarily biofilm structural taxa has been maintained throughout African hominid evolution, and these microbial groups are also shared with howler monkeys, suggesting that they have been important oral members since before the catarrhine-platyrrhine split ca. 40 Mya. However, community structure and individual microbial phylogenies do not closely reflect host relationships, and the dental biofilms of Homo and chimpanzees are distinguished by major taxonomic and functional differences. Reconstructing oral metagenomes from up to 100 thousand years ago, we show that the microbial profiles of both Neanderthals and modern humans are highly similar, sharing functional adaptations in nutrient metabolism. These include an apparent Homo-specific acquisition of salivary amylase-binding capability by oral streptococci, suggesting microbial coadaptation with host diet. We additionally find evidence of shared genetic diversity in the oral bacteria of Neanderthal and Upper Paleolithic modern humans that is not observed in later modern human populations. Differences in the oral microbiomes of African hominids provide insights into human evolution, the ancestral state of the human microbiome, and a temporal framework for understanding microbial health and disease.

The oral cavity is colonized by one of the most diverse sets of microbial communities of the human body, currently estimated at over 600 prevalent taxa (1). Dental diseases, such as caries and periodontitis, remain health burdens in all human populations despite hygiene interventions (2, 3), and oral microbes are often implicated in extraoral inflammatory diseases (4, 5). To date, most oral microbiome research has focused on clinical samples obtained from industrialized populations that have daily oral hygiene routines and access to antibiotics (1, 6), but far less is known about the global diversity of the oral microbiome, especially from diverse past and present nonindustrialized societies (7). The oral cavity contains at least six distinct habitats, but dental biofilms, including both supra- and subgingival dental plaque, are among the most diverse and clinically important (1, 6, 8). During life, these dental biofilms naturally and repeatedly calcify, forming dental calculus (tooth tartar) (9), a robust, long-term record of the oral microbiome (10). Archaeological dental calculus has been shown to preserve authentic oral bacterial metagenomes in a wide range of historic and prehistoric populations and up to 50 thousand years ago (ka) (10–13). As such, dental calculus presents an opportunity to directly investigate the evolution of the hominid microbiome and to reconstruct ancestral states of the modern human oral microbiome. In addition, because research has shown that evolutionary traits, diet, and cultural behaviors shape modern human microbiome structure and function at other body sites, such as the gut and skin microbiomes (14–18), investigating ancient oral metagenomes has the potential to reveal valuable information about major events in modern human evolution and prehistory, such as predicted dietary changes during the speciation of Homo (19–21) and the direct interaction of Neanderthals and modern humans during the Late Pleistocene (22).

To better understand the evolutionary ecology of the African hominid microbiome, we generated and analyzed 109 dental calculus metagenomes from present-day modern humans (n = 8), gorillas (Gorilla, n = 29), chimpanzees (Pan, n = 20), Neanderthals (n = 13), and two groups of archaeological modern humans associated with major lifestyle transitions (preagricultural, n = 20; preantibiotic, n = 14), as well as New World howler monkeys (n = 5) for comparison (). To account for potential sampling biases, we analyzed multiple subspecies and populations of each African great ape genus, which were obtained from C20th or C21st-collected museum collections, and for modern humans we sampled multiple populations from both Africa and Europe. To this, we added previously published microbiome data from chimpanzees (n = 1) (13), Neanderthals (n = 4) (13), and present-day modern humans (n = 10) (23), for a total dataset of 124 individuals (Fig. 1, , and Dataset S1). We also generated eight new radiocarbon dates for archaeological individuals, for a total of 44 directly or indirectly dated ancient individuals in this study (Dataset S1).

Fig. 1.

Sample locations and oral microbiome authentication of ancient dental calculus. (A) Sample locations. (B) PCoA comparing euclidean distances of microbial genera of well-preserved ancient and present-day dental calculus to environmental proxy controls (degraded archaeological bone) and present-day dental plaque and feces. Ancient dental calculus is distinct from gut and archaeological bone but overlaps with present-day dental plaque. (C) Representative DNA damage patterns for Neanderthals and ancient and present-day modern humans for four oral-specific bacterial species. The Neanderthal and upper Paleolithic modern human individuals show expected damage patterns consistent with authentic aDNA, whereas the present-day individual does not. See also .

Here, we investigate the structure, function, and core microbial members of the human oral microbiome within an evolutionary framework, seeking to determine whether a core microbiome can be defined for each African hominid group, whether the core is phylogenetically coherent, and whether some members of the core are specific to certain host groups. We test whether the oral microbiome of hominids reflects host phylogeny, finding that African hominid oral microbiota are distinguished by major taxonomic and functional differences that only weakly reflect host relationships and are likely influenced by other physiological, dietary, or behavioral factors. We compare the microbial profiles of Neanderthals and modern humans and, contrary to expectations (12, 13), find a high consistency of oral microbiome structure within Homo, regardless of geography, time period, or diet/lifestyle. We detect the persistence of shared genetic diversity in core taxa between Neanderthals and Upper Paleolithic humans prior to 14 ka, supporting a growing body of evidence for earlier admixture and interaction in Ice Age Europe (24, 25). Finally, we explore possible implications of our findings on Homo-associated encephalization (19, 26) and the role of dietary starch in human evolution (20, 21) by investigating the evolutionary history of amylase-binding capability by oral streptococci. We find that amylase binding is an apparent Homo-specific trait, suggestive of microbial coadaptation to starch-rich diets early in human evolution.

Results

Preservation of Oral Microbiota in Dental Calculus.

Authenticating ancient DNA (aDNA) preservation is a necessary and essential step for all paleogenomic studies. However, these methods have been underdeveloped for ancient microbiomes. Here, we apply a multistep procedure of both conventional and new methods to evaluate and validate oral microbiome preservation in our dataset (). First, we applied a reference-based metagenomic binning of reads to the National Center for Biotechnology Information (NCBI) nucleotide (nt) database (27) (Dataset S2) and then developed and applied a method to assess the decay of the cumulative percentage of known oral taxa in samples compared to a panel of oral and nonoral reference metagenomes (). This allowed us to remove samples that did not exhibit a taxonomic composition consistent with an oral origin. We then cross validated these results using SourceTracker (28) () and inspection by principal coordinate analysis (PCoA, ). To samples exhibiting good oral microbiome preservation (Fig. 1), we then applied the R package decontam (29) to detect and remove putative laboratory and environmental contaminants prior to downstream analysis (). Next, we examined each dataset and confirmed the presence of DNA damage characteristics of ancient samples, including short fragment lengths and elevated levels of cytosine to thymine deamination (Fig. 1 and ) (30). Finally, to reduce potentially spurious assignments for compositional analysis, we removed low-abundance taxa using thresholds optimized at different taxonomic levels (). The resulting 89 well-preserved dental calculus datasets consist of samples ranging from the present day up to 100 ka.

The Core African Hominid Oral Microbiome.

We performed PCoA on our dataset of well-preserved samples and found considerable overlap in the microbial composition of African hominid dental calculus, as well as howler monkeys (Fig. 1), suggesting the existence of a core microbiome that has been maintained for more than 8 My, based on fossil and molecular evidence of host divergence among African hominids (31, 32), and possibly since before the catarrhine–platyrrhine split ca. 40 Mya (33, 34). At the same time, small but significant differences were indicated by Permutational Multivariate Analysis of Variance (PERMANOVA) (35) at both the microbial genus and species levels between each hominid genus (100 bootstrap replicates, ɑ = 0.05; genus: F = 5.22 ± 1.42, df = 3, R2 = 0.27 ± 0.05, P = 0.001; species: F = 6.67 ± 2.52, df = 3, R2 = 0.32 ± 0.07, P = 0.001; ), and this pattern remained robust after controlling for unequal sample sizes ().

Dental plaque biofilms in humans form by the microbial succession of early, bridging, and late colonizers (36), and in contrast to the gut, which has high interindividual variability at the microbial phylum level (37) and is sensitive to subsistence changes over short and long timescales (15, 38, 39), oral microbial communities have been found to be more stable and consistent, particularly at the genus level (40–42), and even when challenged by antibiotics (43). Because of this, we sought to begin to define the African hominid core oral microbiome as a group and for each genus separately. For a microbial taxon to be considered “core” (44, 45), we required it to be present in at least two-thirds of the populations making up a given host genus, counting as present only those populations in which it is found in at least half of individuals to account for variation in preservation (). We then calculated the intersection of each core microbial genus (Fig. 2) and species (Fig. 2) across all host taxa (Dataset S3). Most “core” taxa are shared across all three African hominid genera (Gorilla, Pan, and Homo) and howler monkeys, whereas fewer are “core” only to African hominids (Gorilla, Pan, and Homo), Pan and Homo, or Homo (Fig. 2 and ). Despite smaller sample sizes than studies of present-day microbiomes, bootstrapping analysis to assess consistency of calculations supported most core microbiome assignments, with lower values possibly indicating taxa influenced by factors such as biofilm maturity (). This suggests a high degree of genus-level microbial taxonomic conservation during African hominid, and possibly broader primate, host evolution and speciation.

Fig. 2.

Core oral microbiome of African hominids shows a deep evolutionary conservation of biofilm structure. UpSet plots showing the number of microbial genera (A) and species (B) core to host groups and group combinations. (C) Core taxa of the human oral microbiome (inclusive of all African hominid and howler monkey ranks). Human biofilm spatial organization based on refs. 8 and 100. Taxa are colored by the broadest host group for which they are core. “Other” taxa are those that fall into paraphyletic host groupings (e.g., Alouatta:Homo). Dashed lines separate the biofilm into basal, intermediate, and peripheral regions (100). Taxa with unknown spatial location are marked with an asterisk (*); taxa core to Homo with any combination of other host genera at the species level but not at the genus level are marked with a dagger (†). Reference Dataset S3 for additional information.

Core taxa at both the genus and species levels include well-known members of each stage of plaque biofilm formation (8, 36), including the early colonizers Streptococcus and Actinomyces, the bridging taxa Fusobacterium and Corynebacterium, and the late colonizers Porphyromonas and Treponema, although the latter two are “core” to only chimpanzees and Homo (Fig. 2). Major periopathogens, bacteria associated with periodontal disease, are found among the different host core combinations, and, focusing on Porphyromonas and Tannerella specifically because of their clinical significance today, we find that their major virulence factors are shared across multiple primates and thus are not specific to modern humans (). The presence of periopathogens within the core microbiome supports the hypothesis that they are not pathogens in a conventional sense but rather that their pathogenic character in present-day humans may be related to an imbalance between the biofilm and the host, as has been suggested by recent ecological studies (46). Although some of the African hominid “core” taxa are periopathogens or their close relatives, most core members are known today to play important structural and functional roles in the formation and maturation of plaque, implying deep coevolutionary relationships between these taxa and their hosts.

African Hominid Oral Microbiome Structure Shows a Weak Relationship with Host Phylogeny.

Hierarchical clustering shows that calculus metagenomes tend to cluster by host genus, confirming intragroup similarity, but these relationships exhibit differences from host phylogeny (Fig. 3). We find, for example, that howler monkeys and gorillas fall together in a single clade and a subset of Homo clusters with chimpanzees. With respect to the latter, available metadata do not provide any clear associations with factors such as geography, time period, or disease to explain this pattern (). Overall, gorillas and howler monkeys are characterized by a wide diversity of aerobic and facultatively anaerobic taxa, while chimpanzees have higher levels of obligately anaerobic taxa, including many putative periopathogens (e.g., Porphyromonas gingivalis, Treponema denticola, Tannerella forsythia, Filifactor alocis, and Fretibacterium fastidosum). Neanderthals consistently fall within the diversity of modern humans. Homo is notable for its high abundance of Streptococcus spp., while this genus is found at substantially lower levels in Pan.

Fig. 3.

African hominid dental calculus microbiomes cluster by host genus and other factors. Hierarchical clustering of howler monkeys, chimpanzees, gorillas, Neanderthals, and ancient and present-day modern humans based on species-level prokaryotic taxonomic assignments. Bacterial oxygen tolerance is associated with biofilm maturation stage in modern humans, and colored names indicate species corresponding to Socransky complexes (111) (reference for a summary). Microaerophilic is defined based on the BacDive database and is roughly synonymous to facultative anaerobe. The tree is schematic, and bifurcations are shown until all host genera are represented. Microbial species names are collapsed to genus level. Species and sample names can be located in .

Many of the taxa identified in human and nonhuman primate dental calculus are poorly characterized, making further exploration difficult. Indeed, several species within the human core genera remain unnamed (Ottowia sp. oral taxon 894, Olsenella sp. oral taxon 807) or understudied (Pseudopropionibacterium propionicum, F. fastidiosum) and some even lack genus designations ([Eubacterium] minutum, TM7x, Anaerolinaceae bacterium oral taxon 439). Their absence from most discussions of the modern human oral microbiome points to a major gap in current oral microbiology research, and targeted investigation of these species is needed to identify their functional and structural roles within plaque biofilms (47–49). Host genus patterns in community structure may be influenced by differences in salivary flow or composition (50), as well as differences in diet texture, quality, and nutrient content (51) (). We also investigated microbial community structure within modern humans, but in contrast to previous studies (13), we found no difference among broad dietary patterns or time periods (). These findings accord with the results of modern oral microbiome studies, which also show minimal, if any, broad and sustained compositional changes in response to diet (e.g., refs. 41, 52, and 53). The relative stability of the oral microbiome may be due in part to the extensive community interdependencies (54, 55) that have developed within these biofilms to metabolize complex host salivary glycoproteins, which are the major nutrient source for most members of the oral microbiota (56). This is in contrast to studies demonstrating strong associations between diet and taxonomic/functional composition in modern gut microbiomes (57, 58).

Evolutionary Histories of Oral Microbial Species Reflect Homo Interactions.

We next examined the phylogenies of individual microbial taxa to determine if host evolutionary relationships are reflected at the microbial genome level. To improve genome coverage and reduce potential noise from DNA damage, we selected a representative subset of well-preserved calculus samples across all host genera (n = 19; and Dataset S1) and constructed uracil-DNA glycosylase-treated (UDG) libraries to remove deaminated cytosines (59), which we then deeply sequenced and analyzed together with a subset of four of the present-day modern humans. Genome-level sequence reconstruction from diversity-rich ancient microbiomes is challenging due to both the highly fragmented nature of aDNA and the low relative abundance of each species, which makes strain separation difficult (). Furthermore, a lack of sufficient reference genomes for many commensal oral microbes increases mismapping and noise artifacts when identifying single-nucleotide polymorphisms (SNPs) (). Nevertheless, despite these challenges and using representative genomes from core taxa, we were able to reconstruct phylogenetic trees with high bootstrap support on internal nodes for eight oral bacteria (Fig. 4 and ).

Fig. 4.

African hominid oral taxa cluster phylogenetically by host genus. Selected neighbor-joining SNP-based phylogenetic cladograms of representative core oral microbiome genomes from deep-sequenced calculus metagenomes (). Actinomyces and Tannerella trees are rooted on the branch leading to howler monkeys (Alouatta, blue), Fretibacterium tree is midpoint rooted. Positions refer to non-N nucleotide calls in the alignment. Node values represent node support out of 100 bootstrap replicates. Asterisk (*) represents the Upper Paleolithic individual from El Mirón (EMN001), which consistently falls near Neanderthal individuals. The remaining eight trees, with tip labels, are provided in .

As with compositional analysis, reconstructed genome-level sequences tend to cluster with those from the same host genus but do not closely reflect host phylogeny (Fig. 4 and ). Overall, genome-level sequences reconstructed from gorillas and chimpanzees fall closer to each other than do those of chimpanzees and Homo. Biases from the use of modern human-derived microbial reference genomes may in part contribute to this pattern, but microbial exchange due to overlapping territorial ranges of gorillas and chimpanzees throughout their evolution may also be a contributing factor. Within Homo, Neanderthals consistently group together, indicating shared within-species microbial diversity. However, we also note that the Upper Paleolithic individual from El Mirón in Iberia (18.6 ka) clusters in all trees with Neanderthals, rather than with other Pleistocene hunter-gatherers of the African Later Stone Age or more recent Holocene-era European or African populations. Recently published human genomic data including this individual has revealed that its associated genetic ancestry component was largely displaced across Europe after 14 ka (24, 60) during postglacial warming. Turning to our low-coverage metagenomic datasets, we assessed additional European Upper Paleolithic and Mesolithic groups () and found that they show a similar pattern (albeit at lower resolution), with the oral taxa of individuals dated to before 14 ka mostly falling with Neanderthals and those after 14 ka mostly clustering with present-day modern humans (24, 60). This pattern suggests that the reconstructed oral bacterial genomes from El Mirón reflect a standing microbial diversity in Homo that was present in Europe during the Middle and Upper Paleolithic, but which was later replaced following subsequent migrations of modern human populations from elsewhere. Because oral microbiota are primarily inherited through caregivers (61, 62), additional sampling and ultradeep sequencing of Paleolithic European and Asian dental calculus may prove informative about the poorly understood interaction dynamics between archaic and modern humans.

Homo-Specific Shifts in Oral Biofilm Are Linked to Dietary Starch Availability.

The metabolic potential of a microbial community, which is inferred from its total gene content, can offer insights into biofilm ecology and function that cannot be understood from taxonomy alone. To better characterize the metabolic and functional differences among hominid oral biofilms, we compared the gene content of dental calculus metagenomes from well-preserved samples of the larger sequencing dataset using two different methods of functional classification, HUMAnN2 (63) and AADDER (64), and found moderate concordance in overall results. Principal Components Analysis (PCA) of the protein-level functional assignments cluster host genera distinctly with a high degree of separation between hosts and functional content (Fig. 5 and ), whereas we observe only a moderate degree of separation in the taxonomic PCoA (Fig. 1), suggesting that gene content of the taxa shared by hominids is more host-specific than taxonomic assignments, a pattern that has also been seen for other microbial systems (65). The genes that drive separation of Homo from nonhuman primates consistently relate to carbohydrate processing (), are much more abundant in Homo, and largely derive from Streptococcus (), something also observed in primate gut microbiomes (66). We therefore investigated the distribution of Streptococcus across our samples (Fig. 5) using a classification system based on biochemical characteristics and genetic relatedness (67). We find that Streptococcus species belonging to the Mitis, Sanguinis, and Salivarius groups are dominant in Homo, while these same groups are effectively absent in chimpanzees, and nonhuman primates in general are characterized by much higher proportions of Streptococcus species in the Anginosus, Mutans, and Pyogenic groups (Fig. 5).

Fig. 5.

Metabolic function and Streptococcus amylase-binding gene content is distinct between African hominid oral microbiomes. (A) PCA of microbial gene functions (SEED classification) clusters well-preserved samples by host genus (PERMANOVA R2 = 0.345). Homo is functionally distinct from nonhuman African hominids and howler monkeys, particularly with respect to carbohydrate metabolism (). (B) Bar plot of proportion of alignments to different Streptococcus groups show differences between host genera. Color of squares below bars corresponds to legend in C. Amylase-binding activity has been observed among members of the Sanguinis, Mitis, and Salivarius groups (68). (C) Ratios of reads aligning to amylase-binding-protein annotated sequences versus a genus-wide Streptococcus “superreference” show higher values in Homo than nonhuman primates, based on a deep-sequenced subset of samples and four present-day modern humans. Note the ratio on the y-axes of abpA and abpB are scaled differently.

The Mitis, Sanguinis, and Salivarius groups are notable for their ability to express amylase-binding proteins to capture salivary ɑ-amylase (68, 69), which they use for their own nutrient acquisition from dietary starch, as well as dental adhesion (70, 71). Amylase-binding protein genes (e.g., abpA and abpB) share no homology but rather confer a similar phenotype through convergent evolution, and they are found almost exclusively in oral Streptococcus species (68). Alpha-amylase is the most abundant enzyme in modern human saliva and modern humans express it at higher levels than any other hominid (50, 72). In contrast to most other nonhuman primates, modern humans exhibit high salivary ɑ-amylase (AMY1) copy number variation, with a reported range of up to 30 diploid copies (16, 73, 74). This copy number expansion is estimated to have occurred along the modern human lineage after the divergence from Neanderthals in the Middle Pleistocene (75, 76). It has been argued this increase relates to dietary shifts during the evolutionary history of modern humans and specifically to an increased reliance on starch-rich foods (20, 21).

We next calculated the ratio of reads aligning to abpA and abpB sequences compared to all Streptococcus reads in the deep-sequenced dataset. We find that abpA and abpB reads are nearly absent in the nonhuman groups but are prevalent and significantly more abundant in Homo (Mann–Whitney U test Homo versus non-Homo: abpB, α = 0.05, U = 128, P = < 0.001, 95% CI = 0.686 to 0.851; abpA, α = 0.05, U = 112, P = < 0.001, 95% CI = 0.398 to 0.861). In particular, abpB is present in all deeply sequenced Homo individuals, and abpA is especially prevalent in modern humans (Fig. 5). This suggests that oral streptococci evolved in association with changes in host diet and supports an early importance of starch-rich foods in Homo evolution.

Discussion

Commensal microbes of the oral microbiome represent an underutilized and independent source of information about host evolutionary and ecological differences (15, 77). With generation times orders-of-magnitude shorter than their hosts and the ability to acquire new functions through horizontal gene transfer across distantly related groups, microbes are a particularly dynamic and temporally resolved system for understanding human evolution. After applying a rigorous strategy to identify, decontaminate, and authenticate well-preserved dental calculus specimens up to 100 ka, we identify a core group of 10 bacterial genera within the African hominid primate oral microbiome that are also shared with howler monkeys, suggesting that these microbial groups have played a key role in oral biofilms since before the catarrhine–platyrrhine split ca. 40 Mya (33, 34). Today, these core taxa are primarily involved in providing structural support within the dental plaque biofilm, and their study holds promise for understanding biofilm growth and maturation in the ancestral human microbiome (78, 79). Identifying the role of such taxa is critical for the successful long-term treatment, prevention, and control of dysbiotic biofilms, such as those found in dental and periodontal diseases (80).

Further, we identify 27 genus-level members of the Homo core oral microbiome, and these include many well-known and clinically relevant taxa, such as Streptococcus and the periopathogens P. gingivalis, T. forsythia, and T. denticola; however, nearly all of these are also core microbiome members of other African hominids. Only Veillonella parvula, a commensal species known to have a synergistic relationship with the cariopathogen Streptococcus mutans (81), is primarily found in humans. Surprisingly, not all members of the core Homo oral microbiome are well-known—three have no genus designation and several lack species names, revealing a major gap in oral microbiology research that in part relates to the difficulties in growing and propagating these microbes.

Focusing on oral microbiome evolution within Homo, we reconstruct authentic oral metagenomes of Neanderthals dating up to 100 ka and modern humans dating up to 30 ka, finding a high degree of similarity in microbial community structure, while also documenting indications of strain-level differences within core taxa. Interestingly, we find that Neanderthal-associated strain-level sequence variants are consistently present in Upper Paleolithic Europeans but not afterward, which accords with a described modern human genomic turnover around 14 ka (24, 60). Comparing human and nonhuman primates, we show that within Streptococcus, amylase-binding groups play a central role in the oral biofilms of Homo, likely aided by both their enhanced ability to colonize the dentition and their exclusive access to dietary starches. These Streptococcus groups and abpB are a general feature of Homo, suggesting that starch-rich foods, possibly modified by cooking (20) (), first became important early in Homo evolution prior to the split between Neanderthal and modern human lineages more than 600 ka (82, 83), a finding with potential implications for the energetics of Homo-associated encephalization (19–21, 26). Subsequent copy number expansion of AMY1 in the modern human genome and the rise of abpA in oral streptococci may signal an even greater reliance on starch-rich foods by modern humans.

Further research on the evolution of abpA, abpB, and other amylase-binding proteins, including phylogenetic reconstruction and demographic modeling, promises to refine questions regarding biofilm formation and the nature and timing of dietary change in Homo. In addition, future research on non-African hominids (orangutans) and additional catarrhines, in particular cercopithecenes with high or unusual salivary amylase expression, such as gelada and hamadryas baboons (73, 84, 85), may yield further insights into the diverse evolutionary trajectories of primate oral microbiomes in response to habitat and dietary change. In addition, it is clear that more research on core genera is urgently needed, as many of the highly conserved and potentially key structural taxa in hominid oral biofilms are understudied and even lack formal names. Furthermore, future sequencing projects focusing on within-species genomic diversity will be critical to understanding microbiome evolution and coadaptation within the human lineage. This study demonstrates that integrating evolutionary studies of the modern human microbiome with wild primate and ancient Homo metagenomic data provides valuable insights into the ancestral states of the human oral microbiome, the nature of microbial–host relationships, and major events in the evolution of modern humans and Neanderthals.

Materials and Methods

Materials.

Our sampling strategy aimed to collect dental calculus from a minimum of two independent populations, each consisting of at least five individuals, for each host genus and modern human lifestyle group (excepting Alouatta) ( and Dataset S1). Dental calculus was sampled from twentieth-century skeletal remains of wild Alouatta (A. palliata), Gorilla (G. berengei beringei; G. berengei graueri; G. gorilla gorilla), and Pan (P. troglodytes schweinfurthii; P. troglodytes ellioti; P. troglodytes verus) and from archaeological Neanderthals and modern humans using established protocols (DOIs: 10.17504/protocols.io.7vrhn56 and 10.17504/protocols.io.7hphj5n). Although many present-day human dental plaque datasets are publicly available, they have been shown to not be directly comparable to dental calculus (23), and consequently we generated dental calculus data for present-day humans. The study of deidentified present-day dental calculus was approved by the Institutional Review Board for Human Research Participant Protection at the University of Oklahoma (IRB no. 4543). All samples were collected under informed consent during routine dental cleaning procedures by practicing dental odontologists. For additional sample context descriptions and additional ethical approval information, reference .

Laboratory Methods.

For all museum and field station samples, we performed DNA extraction in dedicated cleanroom facilities using a protocol optimized for the recovery of degraded and fragmentary DNA (86). Present-day calculus was extracted as previously described (23). For all samples, DNA was built into dual-indexed Illumina libraries (87) and shotgun sequenced. In addition, a subset of samples were separately subjected to UDG treatment (88), followed by deep sequencing. Negative controls were included in all extraction and library construction batches. Sequencing was performed on either Illumina NextSeq. 500 or HiSeq. 4000 platforms. For details, reference and protocols.io under DOI: 10.17504/protocols.io.bq7wmzpe.

Data Processing and Quality Filtering.

For detailed descriptions of preprocessing and analysis procedures, including code, reference and external data repository (GitHub repository: https://github.com/jfy133/Hominid_Calculus_Microbiome_Evolution; Archive DOI: 10.5281/zenodo.3740493). Additional ancient (13) and present-day dental calculus (23) data from previous studies were downloaded from the Online Ancient Genome Repository (OAGR) (https://www.oagr.org.au/) and the European Bioinformatics Institute (EBI) European Nucleotide Archive (ENA) (https://www.ebi.ac.uk/ena/) databases, respectively. Comparative metagenomes from present-day modern human microbiome and environmental sources were additionally downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra/). Accession numbers and download instructions for all FASTQ files are provided in . The EAGER pipeline (89) was used to perform initial preprocessing of sequencing data to remove possible modern human DNA sequences that can interfere with taxonomic profiling (due to present-day modern human DNA contamination in microbial reference genomes). We used relaxed bwa aln (90) mapping parameters for aDNA (−n 0.01), and nonhuman reads from replicate samples and libraries were then concatenated per individual. Human-mapped sequences were then poly-G clipped prior to reporting of mapping statistics. Processing statistics are provided in .

Taxonomic Binning and Preservation Assessment.

For taxonomic binning, we used the aDNA-optimized high-throughput aligner MALT (27, 91) together with the NCBI nt database (October 2017; uploaded to Zenodo under DOI: 10.5281/zenodo.4382154) and a custom NCBI RefSeq database (containing bacteria, archaea, and Homo sapiens, October 2018, ) and employed a relaxed percent identity parameter of 85% and a base tail cut off (“minimum support”) of 0.01%. Resulting RMA6 files were loaded into MEGAN6 CE (64) and prokaryotic Operational Taxonomic Unit (OTU) tables were exported (Dataset S2). A comparison of the two databases is provided in . Given the challenges of low preservation and contamination in ancient microbiome studies, we performed a multistep procedure to screen for and remove poorly preserved samples and contaminant OTUs from the non-UDG-treated dataset (). We developed a visualization for the identification of calculus samples with weak oral microbiome signatures (). This procedure involves comparing identified taxa to their previously reported isolation source(s), ranking these taxa from most to least abundant, and tracking the cumulative percentage of oral taxa along this rank (termed here as “decay”). Samples with a low percentage of oral taxa after an initial “burn-in” based on stabilization of curve fluctuation were removed from downstream analysis. Reference for details. We compared this method to results obtained using SourceTracker (28)—which was performed on 16S-mapped reads filtered from shotgun data using EAGER (with comparative present-day modern human and environmental metagenomes as sources), followed by closed-reference clustering using QIIME (92)—and found concordance between the two methods (). We next used the R package decontam (29) to statistically detect putative laboratory and environmental contaminants (as present in negative controls and a set of archaeological bone samples—), which were then removed prior to downstream analysis. To authenticate the remaining OTUs, we utilized the output of MaltExtract (93) in the MaltExtract-Interactive Plotting App (MEx-IPA) tool (DOI: 10.5281/zenodo.3380011), which we developed for rapid visualization of characteristic aDNA patterns, such as cytosine to thymine deamination, short fragment lengths, and edit distance from reference (). After mapping with EAGER to well-known oral taxa, we also validated DNA damage patterns using DamageProfiler (Fig. 1) (94).

Microbial Compositional Analysis.

To remove low-abundance environmental contaminants or spurious hits, we selected a minimum abundance cutoff of 0.07% of alignments for genus-level and 0.04% of alignments for species-level identifications (). We normalized profiles through phylogenetic isometric-log-ratio transformation (95) of the abundance-filtered OTU tables and then performed PCoA on the resulting euclidean distances ( and ). To statistically verify host genus clusters, we used the adonis function from the R package vegan to perform PERMANOVA (35) analysis after controlling from unequal sample sizes (). After removal of poorly preserved samples, oral communities show distinct centroids for each host genus (bootstrapped PERMANOVA, ɑ = 0.05, P = 0.001, pseudo-F = 5.23, R2 = 0.28); Alouatta was excluded due to small sample size. We performed bootstrapped hierarchical clustering (96) on the euclidean distances of centered log ratio–transformed OTU tables and visualized the results in the form of a heatmap (Fig. 3 and ). Sample and taxon clustering was performed with the McQuitty hierarchical clustering algorithm, and taxon blocks within the heatmap were selected by visual inspection. Bootstrap values of sample clusters were estimated through the R package pvclust (96). Species oxygen-tolerance metadata was obtained from the BacDive database (97) via the BacDiveR R package (DOI: 10.5281/zenodo.1308060). For validation of the observations made on the heatmaps, we also performed grouped indicator analysis (98) (). Clustering of human oral microbiomes by variables such as time, geography, and dietary subsistence was assessed using PCoA, PERMANOVA, and hierarchical clustering ().

Core Microbiome Analysis.

Using the contaminant-filtered OTU tables of well-preserved samples, we converted all taxa above the minimum support threshold to a presence/absence profile. Taxa were required to be present in at least half (50%) of the members of a population for it to be considered core to the population and to be present in at least two-thirds (66%) of populations to be considered core to a host group (; reference for parameter experimentation details). We then generated UpSet plots (99) to visualize the microbial intersections of each host group at both the species and genus levels (Fig. 2 ), and we also compared the results between both databases. Further discussion on the exclusion of the common soil genus Mycobacterium from core genera is provided in . Validation of results through smaller sample sizes was carried out by bootstrapping analysis, which was performed by randomly subsampling (with replacement) individuals from each host genus and rerunning the core calculation procedure to 1,000 replicates (). We created a diagram of the core human oral microbiome (Fig. 2) based on published fluorescence in situ hybridization (FISH) images of human dental plaque (8, 100). For species/genera that were not analyzed in these publications, literature searches were performed to find evidence of their localization within plaque based on immunohistochemistry, immunofluorescence, or FISH (). All members of the human core microbiome are shown, including those also shared with other African hominids and howler monkeys. For further details, reference .

Genomic Analysis.

We used EAGER to map (see below for more details) the deep-sequenced UDG-treated dataset and four samples from present-day individuals (Alouatta, 3; Gorilla, 3; Pan, 4; Neanderthal, 3; ancient modern human, 6; present-day modern human, 4; total: 23) against the reference genomes of Tannerella forsythia and Porphyromonas gingivalis (). We used bedtools (101) to calculate the breadth and depth coverage of a set of known virulence factors for these two taxa. To reduce the risk of spurious alignments (e.g., from cross mapping of conserved sequences), we filtered out genes that had a breadth of coverage less than 70% and/or that appeared to have strongly different coverage depths compared to the rest of the genome (reference for more details). The resulting genes were visualized as a heatmap for comparison (). We selected all species-level Streptococcus alignments from the shallow sequenced dataset minimum support filtered NCBI nt–MALT OTU tables and assigned them to one of eight species groups based on the literature (67) (reference for group definitions). We then calculated the fraction of alignments for each species group over all taxonomic alignments for each sample (Fig. 5). To further validate the results, we calculated a similar ratio but based on the mapping of the deep-sequenced dataset against a superreference of 166 Streptococcus genomes (see below). We identified abpA- and abpB-like gene coordinates from the superreference using panX (102), then extracted the number of reads mapping to these annotations and calculated the fraction of these reads over all Streptococcus superreference mapped reads. We then applied a Mann–Whitney U test to test the null hypothesis of no difference between the distributions of ratios of Homo and nonhuman primates, as well as compared these results to a distribution of P values of 100 randomly shuffled group assignments (reference for more details). Reference sequences of abpA and abpB were extracted from Streptococcus genomes in RefSeq and indexed for mapping. All shallow sequencing dataset samples were mapped against all reference strains. For samples with a gene coverage of at least 40% at 1×, a consensus sequence was exported from the Integrative Genome Viewer (IGV) (103). An input file of the consensus sequences and references was generated in BEAUTi and used to run BEAST2 (104) for Bayesian skyline plot analysis. For details, reference .

Microbial Phylogenetics.

We first attempted a competitive-mapping strategy against genus-wide superreferences of identified core taxa (reference ), but this approach yielded only limited results (). We then instead performed phylogenetic reconstruction by mapping the same dataset to a single representative genome for each genus, considered as representing a population of related taxa. To account for challenges with low-coverage ancient data, we called SNPs using MultiVCFAnalyzer and required each SNP call to have a minimum of 2× coverage and a support of ≥70% of reads (). The resulting FASTA alignments were loaded into R. Samples with fewer than 1,000 SNPs were removed, and pairwise distances were calculated based on the JC69 model (105). A bootstrapped neighbor-joining algorithm from the R package ape (106) was applied to the distance matrices with 100 replicates (). Trees were visualized with ggtree (107). Finally, we retained trees where the basal internal nodes had bootstrap supports of ≥70% (). The same procedure was then applied to the shallow sequencing dataset with the additional samples described above in the main text (). To test whether pre-14 ka individuals clustered with Neanderthals due to reference bias, we calculated the median number of positions that were shared between EMN001 and Neanderthals to a histogram of median pairwise comparisons between all modern human individuals ().

Functional and Metabolic Pathway Analysis.

We took two approaches to characterizing the functional profiles of the calculus metagenomes. First, we used HUMANn2 (63) [with MetaPhlAn2 (108) generated taxonomic profiles] to generate functional profiles based on the UniRef90 (109) and ChocoPhlAn (July 2018) (63) databases. Preservation was independently assessed for pathway abundance and Kyoto Encyclopedia of Genes and Genomes (KEGG) ortholog functional profiles, . We compared the functional profiles of well-preserved calculus between host groups using pathway abundance (n = 94) and gene families converted to KEGG orthologs (n = 109) using PCA (). Orthologs with the strongest loadings were visualized with biplots (), and the species from which these orthologs were derived were determined (). The clustering of host genera in PCAs using only orthologs in specific pathways (carbohydrates, amino acids, lipids) was also explored (). For details, reference . Second, we used AADDER (included within MEGAN6 CE) (64) to profile the number of alignments to annotations present in the custom RefSeq database as aligned by MALT (see above). We then used MEGAN6 to export SEED category (110) profiles. Preservation was independently assessed for SEED protein functional profiles, reference . We compared the functional profiles of well-preserved calculus (n = 95) between host groups using proteins but not higher-level pathways (). The proteins with the strongest loadings were visualized using biplots (), and the species from which these proteins were derived were determined (). The clustering of host genera in PCAs using only proteins in specific pathways (carbohydrates, amino acids, lipids) was also explored (). For details, reference .