Benchmarking a targeted 16S ribosomal RNA gene enrichment approach to reconstruct ancient microbial communities

16S rRNA gene RNA probe design

Full-length 16S rRNA genes were obtained from all species identified from a previous ancient dental calculus study (Weyrich et al., 2017) from the Ribosome Database Project (RDP) (Cole et al., 2014). To further increase diversity, 16S rRNA genes from species in the Human Oral Microbiome Database (HOMD) (Chen et al., 2010) that were not identified by Weyrich et al. (2017) were also included. This added 285 sequences, which resulted in 570 full-length 16S rRNA genes to be used for probe design (FigShare DOI: 10.25909/5cc11894b0cc2). RepeatMasker (Smit, 2004) removed simple and low-complexity repeats from the sequences. A total of 80 bp (base pair) RNA baits with 20 bp probe spacing and 4× tiling density were synthesized by Arbor Biosciences (formerly MyBaits). Baits were collapsed if they had fewer than 11 mismatches between each other, yielding a total of 19,634 80 bp RNA baits. While the probe design was based on 16S rRNA genes from ancient and modern oral taxa, the redundancy built into the probes should also capture microbial diversity that was not used as input for the probe sequences.

Hybridization enrichment and DNA sequencing

A total of four samples that were previously published in Weyrich et al. (2017) using unenriched, whole-genome alignment approach (UeWGA) methods were selected for hybridization enrichment (Table S1). In brief, each library was created by amplifying existing library DNA in four 25 µL PCR reactions (each containing: 13.625 µL dH₂O, 2.5 µL 10X AmpliTaq Gold Buffer, 2.5 µL of 25 mM MgCl₂, 0.625 µL of 10 mM dNTPs, 2.5 µL of 10 µM forward and reverse primer, 0.25 µL AmpliTaq Gold, 3 µL template DNA) with the following PCR conditions: denaturing at 94 °C for 12 min before 13 cycles of (30 s denaturing at 94 °C, 30 s annealing at 60 °C, 45 s extension at 72 °C) followed by a final extension of 10 min at 72 °C (Eisenhofer, 2018). PCR amplifications were pooled and then cleaned with AMPure XP beads to reach 100 ng of DNA for input into hybridization enrichment. A modified version three of the MyBaits protocol was used for hybridization enrichment, whereby the input RNA bait concentration was reduced to 25% of the recommended amount, and custom oligonucleotides were utilized to block the P5/P7 library adapters. To test the efficiency of hybridization at different temperatures, four samples were enriched at 55 °C and 65 °C for 40 h. Based on these results, the rest of the samples were enriched at 65 °C for 40 h. Following hybridization enrichment capture, the streptavidin beads were washed three times using Wash Buffer 2 (MyBaits v3 manual) and resuspended in a PCR mastermix (as above) before amplification with 13 cycles of PCR. Amplified libraries were cleaned with AMPure, quantified with an Agilent TapeStation, pooled at equimolar concentrations, and sequenced on an Illumina HiSeq X Ten platform (2 × 150 bp). Sequencing data for the samples processed with the hybridization capture method can be found on the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under the accession number: PRJNA1031168.

Sequence data processing and 16S rRNA gene enrichment analysis

The sequencing data was converted into the fastq format using Illumina’s bcl2fastq software and were demultiplexed using AdapterRemoval2 based on unique P5/P7 barcode combinations (Schubert, Lindgreen & Orlando, 2016). We reduced the computational time of downstream analyses, seqtk (https://github.com/lh3/seqtk) was used to randomly subsample each 16S rRNA enriched library to 100,000 reads. SortMeRNA (Kopylova, Noé & Touzet, 2012) was used with the SILVA bacterial and archaeal 16S rRNA databases (SILVA 119, 95% clustered) to identify 16S rRNA gene fragments from samples. We constructed a BLAST database (Altschul et al., 1990) using the SILVA 128 NR99 16S rRNA reference database (Quast et al., 2012) to assign taxonomic identifications for the putative 16S rRNA gene fragments. The putative 16S rRNA gene fragments were aligned to the database using the default parameters, with the following exceptions: –evalue = 0.01 for added stringency and –outfmt = 0 to allow for import into MEGAN. We tested the impact of missing reference sequences on taxonomic classification by creating a filtered SILVA database and removed reference sequences of taxa from of our simulated dataset using the Filterbyname.sh script from BBtools (https://jgi.doe.gov/data-and-tools/bbtools/). We used grep to create a list of taxa to remove from the reference sequences in the SILVA database using Filterbyname.sh (Table S2). Lastly, BLAST outputs were then imported into MEGAN CE (v. 6.24.15) (Huson et al., 2016) using the “Import from BLAST” option with the synonyms mapping file (SSURef_NR99_128_tax_silva_to_NCBI_synonyms.map.gz) obtained from the MEGAN community website: (https://software-ab.cs.uni-tuebingen.de/download/megan6/welcome.html).

Generating simulated dataset with ancient DNA characteristics

To test if our 16S rRNA gene fragment analysis pipeline could accurately reconstruct microbial communities, we constructed a simulated dataset using Gargammel (Renaud et al., 2017). 16S rRNA genes from 19 phylogenetically diverse prokaryotic species were obtained from the SILVA 128 NR99 16S rRNA database (Quast et al., 2012) (Table S2) and randomly fragmented to create 100,000 16S rRNA gene fragments fitting a log-normal aDNA fragment-length distribution (Fig. S1) (gargammel -n 100,000 --loc 4 --scale 0.3). Varying levels of cytosine deamination (aDNA damage) were simulated using deamSim from gargammel on our simulated dataset to create three different levels of simulated single-stranded overhang deamination: 10% (low) deamination (-damage 0.03, 0.25, 0.01, 0.1); empirical (moderate) deamination from a published mapDamage profile (Jónsson et al., 2013; Olalde et al., 2014); and 50% (high) deamination (-damage 0.03, 0.25, 0.01, 0.5). Lastly, the resulting simulated 16S rRNA gene datasets were mapped against the SILVA 128 NR99 database and the filtered SILVA database using BLAST, as described above.

Taxonomic classification of dental calculus data processed with the UeWGA pipeline

Sequencing data bioinformatically processed with the UeWGA pipeline come dental calculus samples that were previously published (NCBI SRA database under BioProject PRJNA685265) and steps for this pipeline are described in Eisenhofer & Weyrich (2019). In brief, reads were aligned against a NCBI microbial reference database containing 47,713 genomes from bacteria and archaea using MALTn (Herbig et al., 2016) with default parameters and outputting BLAST text files. The LCA parameters used for the UeWGA data were: bitscore = 50, E-value = 0.01, minsupp = 0.1 (i.e., a taxonomic assignment requires at least 0.1% of reads to pass), and the weighted LCA algorithm (80%), as suggested in Huson et al. (2016). The resulting BLAST text files were converted into RMA6 files using the blast2rma script in MEGAN. Operational Taxonomic Units (OTUs) resulting from a 16S rRNA amplicon approach on the same samples were also obtained from a previous study (Weyrich et al., 2017). For the 16S samples, OTUs were clustered using UCLUST at 97% and taxonomically identified using the Greengenes 12.10 database.

Comparison of 16S rRNA gene enrichment, 16S rRNA amplicon datasets, and UeWGA datasets

LCA parameters used for the 16S enrichment method were the same as the ones in the UeWGA dataset except, except that the naïve LCA algorithm (80%) was used because false-positive taxonomic assignments were identified in the simulated 16S rRNA gene fragments with the weighted algorithm. 16S rRNA amplicon data was directly imported into MEGAN (v. 6.24.15) (Huson et al., 2016).

Statistical analyses

Feature-level (i.e., sequences assigned to the lowest taxonomic rank) and genus-level assignments in MEGAN were exported as BIOM tables and then imported into QIIME2 (v.2021.11) (Caporaso et al., 2012). Samples with less than 1,000 counts were excluded from downstream analyses. For the alpha and beta diversity analyses, unfiltered data sets were compared to ascertain the differential impacts of contaminant taxa. For comparisons of specific taxa across methods, genera and species identified in the extraction blank controls were removed from biological samples in MEGAN. When calculating beta diversity using Bray Curtis, the counts of all samples were rarefied to the sampling depth for the sample with the lowest number of feature-level (9,980 for A12017 EuroHG2) and genus-level assignments (9,731 for A13208_AfrSF2). Alpha diversity measures were calculated based on the number of ‘observed features’ and ‘observed genera’ using the alpha-group-significance plug-in (Kruskal-Wallis test). Phyloseq was also used to perform a rarefaction curve analysis using different alpha diversity indexes (Shannon, Chao1, and total observe features/genera). Bray Curtis distance matrices at the feature and genus level were computed with the qiime diversity core-metrics plugin and used for the Principal Coordinate Analysis (PCoA) and PERMANOVA tests in QIIME2. Aitchison distances were computed with the qiime deicode rpca plugin and used for the PCoAs visualizations and PERMANOVA tests in QIIME2.

Assessing taxonomic assignment of simulated ancient 16S rRNA gene fragments

We created a simulated dataset containing 19 phylogenetically diverse prokaryotic species to test how accurately we could classify short 16S rRNA gene fragments (Table S2). The 16S rRNA gene from each species was selected and randomly fragmented to fit a commonly observed aDNA fragment length distribution (i.e., log normal with a mode of 50 bp; Fig. S1). DNA fragments were classified taxonomically by mapping DNA fragments to the SILVA 16S rRNA database (NR 99, release 128) using BLAST nucleotide alignment, following by applying the default lowest common ancestor (LCA) algorithm in MEGAN6 (v. 6.24.15). We found that 98.4% of the total DNA fragments were assigned taxonomically (Fig. 2A), and of these, 53.2% were assigned to the genus level (Fig. S2), and 3.5% to the species level (Table 1; Fig. S2). The resulting genus-level taxonomic composition matched that of the input sequences with one exception: Yersinia 16S rRNA gene fragments could not be assigned to the genus Yersinia and were pushed up to the order Enterobacterales (Fig. 2B). The 1.6% of fragments that were not aligned were extremely short (between 15–30 bp) (Fig. S1), suggesting that accurate classifications are limited to at least 30 bp for 16S rRNA fragments, as previously observed for whole-genome alignments (Eisenhofer & Weyrich, 2019).

To test if cytosine deamination (a characteristic of aDNA) influences taxonomic classification of short 16S rRNA gene fragments, we simulated three levels of sequence deamination on our simulated dataset representing low (10%), moderate (simulated from a real mapDamage profile from an ancient specimen; –30% (Olalde et al., 2014); and high (50%). We simulated moderate damage with the La Brana individual because the damage profile for this sample has been benchmarked in previous simulation studies and represents a typical ancient sample (Oliva et al., 2021). We found that deamination did not have a notable impact on taxonomic classification, as there was no loss of taxonomic assignments or misclassifications (Figs. 2E–2G). Only a 1.5% loss of reads assigned taxonomically was observed when assessing the highest level of deamination (50%) compared to the non-deaminated simulated dataset (Table 1). Overall, these findings suggest that deamination does not meaningfully hinder the ability to classify ancient 16S DNA fragments and that the method developed here could be used to examine taxonomic diversity even in highly, degraded ancient samples. Our results are consistent with prior research, indicating that deamination has only a marginal effect on taxonomic classification for data generated with shotgun sequencing (Eisenhofer & Weyrich, 2019; Velsko et al., 2018; Warinner et al., 2017).

To test the impact of missing reference sequences on our ability to taxonomically classify 16S rRNA gene fragments, we performed a species and genus exclusion experiment whereby we removed reference sequences corresponding to species or genera present within the simulated dataset from the SILVA database (Table S2), and then repeated BLAST alignments using this modified database. Sequences originating from excluded species could be classified to their respective genera (Fig. 2C, Table 1). Similarly, removal of all reference sequences attributed to the Streptococcus genus (which accounts for 25% of the simulated dataset) resulted in a 12% reduction in assignments to the genus-level (53.2% vs. 41.2%), and an associated increase in sequences at higher taxonomic ranks (Fig. 2D, Table 1). Importantly, the total number of aligned reads did not substantially decrease (97.7% vs. 98.4% for non-exclusion) (Fig. S2), suggesting that sequence conservation and phylogenetic signal within the 16S gene allowed for placement of these sequences higher in the taxonomy, rather than discarding them outright.

Optimizing hybridization enrichment efficiency

We designed RNA baits to capture 16S rRNA gene fragments from a diverse range of microbial taxa. We tested whether the proportion of on-target sequences and non-target sequences impacted by differing hybridization temperatures, specifically 55 °C and 65 °C. We found that enrichment at 65 °C yielded the highest on-target enrichment of 16S rRNA gene fragments (average 53.42% putative 16S rRNA gene fragments). This was more than a two-fold increase over the 55 °C enrichment (average 26.17% putative 16S rRNA gene fragments), and a 334-fold increase over the unenriched shotgun libraries (average 0.16% putative 16S rRNA gene fragments) (Table S3). We found that increasing hybridization temperature selected for longer mean DNA fragment lengths (unenriched = 56 bp, 55 °C = 69 bp, 65 °C = 79 bp), but did not find changes in the mean GC content between temperatures (Table S4). This enrichment of longer DNA fragments by higher temperatures could be explained by shorter DNA fragments not having sufficient Watson-Crick hydrogen bonding sites to remain attached to probes at higher temperatures, resulting in the preferential binding and enrichment of longer DNA fragments.

Hybridization temperature minimally influences diversity

Given that hybridization temperature was found to select for longer DNA sequences, we next tested to see if this influences the alpha and beta diversity of samples. Due to constraints in the alignment speed of BLAST, each sample from each treatment (55 °C, 65 °C) was randomly subsampled to 100,000 reads. After alignment and classification using the SILVA database, genus and feature level assignments were exported separately from MEGAN6 (v.6.24.15) into QIIME2 (v.2021.11). The alpha diversity between the 55 °C and 65 °C groups were insignificant at the feature level (Kruskal-Wallis; p = 0.77) and genus level (Kruskal-Wallis; p = 1) (Table S5). Beta diversity between the two groups were also similar (PERMANOVA; p > 0.05) when using either Bray-Curtis and Aitchison distances at the genus and feature levels (Figs. 3A and 3B; Table 2). Regarding taxonomic assignments specific to hybridization temperature, taxonomic bar plots representing the sum of samples per treatment were nearly identical in both abundance and diversity (Fig. S3). Both enrichment temperatures allowed for the detection of Corynebacterium (a genus commonly found in modern human dental plaque) and was not identified in unenriched samples. The 65 °C treatment was also the only one to detect Pseudoramibacter (an oral taxon previously identified in ancient dental calculus and is present in modern oral studies (Eisenhofer & Weyrich, 2019; Siqueira & Rôças, 2003)) (Fig. S4). Overall, these results suggest that a hybridization enrichment can be used to identify microbial taxa within dental calculus and that temperature of 65 °C may provide a mild improvement on the recovery of 16S rRNA gene fragments. Therefore, a hybridization temperature of 65 °C was used for subsequent enrichments.

Differences in diversity are observed between 16S rRNA enrichment, whole-genome alignment, and 16S rRNA amplification methods

We next sought to examine if the diversity and composition of communities obtained from 16S rRNA gene enrichment were distinct from the UeWGA and 16S rRNA amplified strategies. We compared 16S rRNA enrichments on samples with previously published UeWGA and 16S rRNA amplification data (Eisenhofer & Weyrich, 2019; Weyrich et al., 2017). We enriched 11 additional dental calculus samples from Weyrich et al. (2017) at 65 °C, bringing the total number of samples to 15 (including four samples from the previous section).

Alpha-rarefaction curves indicated that sampling completeness for both feature-level and genus-level was reached at a sequencing depth of 1,000 reads (Shannon’s Diversity) and >3,000 sequences (observed taxa) for 16S enriched, 16S rRNA amplicon, and UeWGA, (Fig. S5). Group significance for alpha diversity was calculated using the Kruskal-Wallis test for both observed features and genera BIOM tables. Results indicate that UeWGA recovered more taxa at the feature and genus level than both the 16S hybridization enrichment and amplicon strategies (Figs. 4A and 4B). The alpha diversity among the three treatments were significantly different when using feature (p = 0.00000009; H = 32.372) and genus-level assignments (p = 0.002; H = 16.923) for each of the metrics examined (q < 0.05; Table S6). Each pairwise comparisons among the three groups were also significant at the feature and genus level (q < 0.05; Table S6) across all three groups, suggesting that distinct microbial profiles were obtained using either of the three methods. This observation corroborates with the findings of previous studies (Weyrich et al., 2017; Ziesemer et al., 2015). The median number of observed features in the 16S rRNA amplicon data were less than both the 16S enriched and UeWGA approaches, likely because the amplicons were too long to recover most of the taxa. It is interesting that the enrichment and the 16S amplicon methods recovered similar number of genera from the samples, since the enrichment method only recovers oral microbes. Nevertheless, the results make it clear that the UeWGA method identifies the greatest number of taxa and had less range of variation than either method, demonstrating that it should still be the preferred choice when surveying the microbial diversity of samples.

When comparing compositions across the methods, PCoA analysis of either Bray Curtis or Aitchison distances resulted in distinct clustering based on library type rather than the individual (PERMANOVA; p = 0.001; Figs. 5A–5D), although the 16S enriched and UeWGA datasets clustered more closely together than the 16S amplicon data sets on Axis 1. This supports previous findings highlighting the strong contaminant biases present in ancient 16S rRNA gene amplified libraries (Weyrich et al., 2017). Biplots indicate that distinct clustering of the samples generated by 16S amplicon sequencing had higher abundances of contaminant taxa at the feature and genus level (e.g., Acientobacter, Comamonas, Comamonadaceae, and Eubacteriales) (Figs. 5B and 5D), which has already been previously noted (Ziesemer et al., 2015).

We also explored differences between enriched 16S and UeWGA methods using a presence/absence approach. Taxa identified in the UeWGA data included Bacteroides, Microbacter, Parabacteriodes, Desulfoplanes, Pseudomonas, Merismopedia, Blautia, Acholeplasma, and Candidatus Phytoplasma (Fig. S5). These taxa are not commonly found in the oral cavity (Chen et al., 2010) and are likely present because of either laboratory or reagent contamination. It is possible that they were identified here because the minimum percentage needed for taxonomic assignment (minimum support percent 0.1%) for the reads of the UeWGA data is lower (i.e., the minimum support percent for 3,000 assignments is (3,000*0.001 = 3)) than the 16S enriched data (minimum support percent for 30,000 (30,000*0.001 = 30)).

Specific taxonomic differences between 16S rRNA enriched and UeWGA data sets

Given the known biases that 16S rRNA amplification methods have when applied to aDNA datasets, we removed samples generated from this method and proceeded to compare the compositional differences between the 16S enrichment and UeWGA methods (Table S7). A pairwise comparison based on Bray Curtis and Aitchison distances indicated that the microbial compositions were significantly different at the genus and feature level (q = <0.05; Table 2). As expected, the UeWGA method identified more species (average of 64.8% of assigned reads at the species level versus 1.8% for the 16S enrichment method). Importantly, the UeWGA method could align 59% of reads, whereas only 0.8% could be aligned for the 16S enrichment method (Table S7). We next tested if there were differences in genus-level assignments between the two methods. The 16S rRNA gene enrichment method classified 20 genera that were not present in the UeWGA method, although these assignments each had a mean abundance <1% (Table 3). The UeWGA method had 18 genus identifications absent in the 16S rRNA gene enrichment method, and all but three of these each had a mean abundance of <1% (Table 3). Comparing these assignments to taxa found in the Human Oral Microbiome Database (HOMD), 9/20 were putatively oral for the 16S enrichment method, while 14/18 for the UeWGA method. Notably, the 16S enrichment method did not detect the common oral genus Tannerella, while the UeWGA not only identified this genus but also had a mean abundance of 3.42% across samples. Additionally, the putative oral genus Olsenella was not detected in the 16S enrichment method but had a mean abundance of 7.85% in the UeWGA data set. These findings suggest that rare or low abundant taxa are likely contributing to the minimal differences in composition between the two methods. As such, the 16S enrichment strategy missed key oral taxa that were identified in the UeWGA approach, meaning future should investigate how to improve probe designs for dental calculus samples to minimize the loss of important taxa.

Comparison of cost estimates for enriched and shotgun methodologies

Given these observable compositional differences, we compared the costs associated with implementing 16S hybridization capture and UeWGA, given equal sequencing depth. When the sequencing data for each method was subsampled to the same depth, the 16S rRNA gene enrichment method assigned an average of 21.1% of sequences a taxonomic classification to the genus and species level, while the UeWGA method assigned 34.5% of the sequences to a similar level. This represents a 13.4% increase in the number of reads assigned at the genus or species level for UeWGA over the 16S rRNA gene enrichment method, despite the increase in the percentage of reads assigned for the 16S rRNA gene enrichment method (56% for 16S rRNA gene enrichment, 41% for UeWGA). This suggests that while the 16S enrichment method can align and assign a higher percentage of sequences, these assignments are not as specific as the UeWGA method, and given the extra costs associated with the 16S rRNA gene enrichment method, the UeWGA method is currently the most cost-effective means of classifying genus or species level taxonomy in highly degraded ancient samples.