Intrinsic challenges in ancient microbiome reconstruction using 16S rRNA gene amplification

Abstract

To date, characterization of ancient oral (dental calculus) and gut (coprolite) microbiota has been primarily accomplished through a metataxonomic approach involving targeted amplification of one or more variable regions in the 16S rRNA gene. Specifically, the V3 region (E. coli 341-534) of this gene has been suggested as an excellent candidate for ancient DNA amplification and microbial community reconstruction. However, in practice this metataxonomic approach often produces highly skewed taxonomic frequency data. In this study, we use non-targeted (shotgun metagenomics) sequencing methods to better understand skewed microbial profiles observed in four ancient dental calculus specimens previously analyzed by amplicon sequencing. Through comparisons of microbial taxonomic counts from paired amplicon (V3 U341F/534R) and shotgun sequencing datasets, we demonstrate that extensive length polymorphisms in the V3 region are a consistent and major cause of differential amplification leading to taxonomic bias in ancient microbiome reconstructions based on amplicon sequencing. We conclude that systematic amplification bias confounds attempts to accurately reconstruct microbiome taxonomic profiles from 16S rRNA V3 amplicon data generated using universal primers. Because in silico analysis indicates that alternative 16S rRNA hypervariable regions will present similar challenges, we advocate for the use of a shotgun metagenomics approach in ancient microbiome reconstructions.

Figure 1. Unusual microbiome profiles observed in 16S rRNA gene V3 amplicon data from archaeological dental calculus.

Relative abundance charts summarizing: (a) Frequency of oral-associated genera in dental calculus and control samples. The Dutch, UK, Nepalese, and Spanish calculus samples show a greater proportion of oral-associated genera compared to the St. Helena and pre-Columbian Caribbean samples. (b) Contribution of oral, gut, and environmental sources to microbiome composition estimated by Bayesian source tracking. The oral microbiome (saliva, supragingival plaque, subgingival plaque) is a major source (>50%) in only a small proportion of archaeological dental calculus samples (10%), and an oral source is not indicated for more than a quarter of samples (26%). Laboratory controls (osteologist hands and osteology lab bench surfaces) and extraction blanks are largely consistent with a skin microbiome source and unassigned contaminants. (c) Frequency of microbial phyla inferred from V3 amplicon sequencing. The taxonomic profile reveals an unusual and non-biological pattern of exceptionally high Euryarchaeota levels in the Dutch, UK, and some Caribbean dental calculus samples. All Euryarchaeota OTUs are assigned to the genus Methanobrevibacter, the only prevalent genus of Archaea in the oral cavity. Methanobrevibacter is typically found at low frequencies (<0.5%) in healthy human dental plaque, but in archaeological samples it may reach frequencies >60%, as seen here. Starred samples (*) were also analyzed using shotgun metagenome sequencing. Sites are ordered from left to right by increasing thermal age (see Table 1). Figure data is available in Supplementary Data 1.

Figure 2. Length distribution box plots of aDNA extracted from archaeological dental calculus and calculated V3 and V4 16S rRNA amplicon lengths for microbes in the SILVA SSU 111 database.

As expected for aDNA, the genetic material within dental calculus is highly fragmented to median lengths <100 bp: 214C, 75 bp; 37C, 77 bp; F1948C, 80 bp; 454C, 91 bp. This is significantly shorter than the median lengths of the 16S rRNA V3 (183 bp) and V4 (292 bp) amplicon targets. The number of read pairs comprising each box plot are as follows: 214C, 20,355; 37C, 17,962; F1948C, 26,517; 454C, 15,736; SILVA V3, 651,163; SILVA V4, 649,660.

Figure 3. Simplified 16S small subunit ribosomal RNA secondary structure.

Secondary structure of Escherichia coli (J01695) 16S rRNA (main panel). Amplicon targets (primer inclusive) for the third (V3, 341F/529R) and fourth (V4, 515F/806R) variable regions are highlighted in pink and blue, respectively. Overlapping V3/V4 target sequences are highlighted in purple. Although widely used in ecological studies, the V4 region is impractical for aDNA research because of its long length (approx. 292 bp, primer inclusive). The V3 region is considerably shorter, but comparative sequence analysis (a–d) reveals that the V3 region exhibits extensive length polymorphisms (arrows) in archaeal (e.g., Methanobrevibacter oralis) and bacterial (e.g., Corynebacterium diphtheria, Streptococcus pyogenes) taxa, with predicted V3 amplicon lengths ranging from 150–194 bp when queried against the SILVA SSU 111 16S rRNA database (d). By contrast, the V4 region is relatively length invariant (e–h), ranging from 290–295 bp (h). 16S rRNA secondary structure has been adapted from Comparative RNA Web Site and Project.

Figure 4. Heatmap of 16S rRNA V3 amplicon lengths reveals high variability but broad taxonomic patterns.

(a) 16S rRNA V3 sequence data was analyzed in silico for 36,634 OTUs belonging to 31 representative oral microbiome genera from 9 major microbial phyla: Euryarchaeota (yellow), TM7 (orange), Chlorflexi (red), Actinobacteria (green), Firmicutes (blue), Fusobacteria (purple), Bacteroidetes (gray), Proteobacteria (pink), Spirochaetes (brown). (b) Log fold changes in genus frequency within archaeological dental calculus when comparing data generated by targeted V3 U341F/534R amplicon sequencing to non-targeted shotgun metagenomics data. Methanobrevibacter, Anaerolinea, and TM7, which have very short predicted V3 U341F/534R amplicon lengths, strongly overamplify compared to frequency data obtained from non-targeted shotgun metagenomic sequencing. Most other taxa underamplify, especially genera with very long predicted V3 U341F/534R amplicon lengths, such as Treponema and members of Proteobacteria, and Bacteroidetes.

Figure 5. Fold changes in taxon frequency between 16S rRNA V3 U341F/534R amplicon and shotgun metagenome data.

Genera with relatively short median amplicon lengths (Methanobrevibacter, Anaerolineae G-1, TM7) are overrepresented in the 16S rRNA V3 U341F/534R amplicon dataset, while genera with relatively long median amplicon lengths (Treponema, Neisseria) are strongly underrepresented.

Figure 6. Predicted effect of thermal age on reconstructed taxonomic frequencies of selected oral bacteria from V3 U341F/534R amplicon data.

Using a random DNA degradation model, the relative abundance of taxa presented in Table 2 is modeled at different thermal ages, corresponding to the thermal ages of the archaeological sites in this study. The probability of chain scission (λ) is estimated using temperatures estimated using previously published kinetic parameters. Starting taxonomic frequencies were taken from mean HOMD values, and taxa with a frequency of <0.01% in the HOMD are not shown. Together, the taxa shown account for 25% of the HOMD human oral microbiome.