Amplicon Sequence Variants Artificially Split Bacterial Genomes into Separate Clusters

116S rRNA gene sequencing is a powerful technique for describing and comparing microbial communities. Efforts to link 16S rRNA gene sequences to taxonomic levels based on distance thresholds date to at least the 1990s. The distance-based threshold that was developed and is now widely used was based on DNA-DNA hybridization approaches that are not as precise as genome sequencing. Instead, genome sequencing technologies have suggested that the widely used 3% distance threshold to operationally define bacterial taxa is too coarse. As an alternative to operational taxonomic units (OTUs), amplicon sequencing variants (ASVs) have been proposed as a way to adopt the thresholds suggested by genome sequencing to microbial community analysis using 16S rRNA gene sequences. It is widely understood that individual bacterial genomes often have multiple 16S rRNA genes that are not identical and that a 16S rRNA gene sequence could be found with different versions of the sequence in different genomes. This could lead to the problem that ASVs and using too fine a threshold to identify OTUs could split a single genome into multiple clusters. Proponents of ASVs minimize concerns that most bacterial genomes have more than one copy of the rrn operon and that those copies are not identical. Conversely, using too broad of a threshold to define OTUs could cluster multiple bacterial species into the same OTU. An example of both is seen in the comparison of Staphylococcus aureus (NCTC 8325) and Staphylococcus epidermidis (ATCC 12228), where each genome has 5 copies of the 16S rRNA gene. Each of the 10 copies of the 16S rRNA gene in these two genomes is distinct, and they represent 10 ASVs. Conversely, if the copies were clustered using a 3% distance threshold, then all 10 ASVs would cluster into the same OTU. The goal of this study was to quantify the trade-off of splitting a single genome into multiple clusters and the risk of clustering different bacterial species into the same cluster when using ASVs and various OTU definitions.

2To investigate the variation in the number of copies of the 16S rRNA gene per genome and the intragenomic variation among copies of the 16S rRNA gene, I obtained 16S rRNA sequences from the rrn copy number database (rrnDB). Among the 5,972 species represented in the rrnDB, there were 20,427 genomes. The median rrn copy number per species ranged between 1 (e.g., Mycobacterium tuberculosis) and 19 (Metabacillus litoralis). As the rrn copy number for a genome increased, the number of variants of the 16S rRNA gene in each genome also increased. On average, there were 0.58 variants per copy of the full-length 16S rRNA gene and averages of 0.32, 0.25, and 0.27 variants when considering the V3-V4, V4, and V4-V5 regions of the gene, respectively. Although a species tended to have a consistent number of 16S rRNA gene copies per genome, the number of total variants increased with the number of genomes that were sampled (see Fig. S1 in the supplemental material). For example, the 271 genome accessions of Mycobacterium tuberculosis in the rrnDB each had 1 copy of the gene per genome. However, across these accessions, there were 17 versions of the gene. An Escherichia coli genome typically had 7 copies of the 16S rRNA gene, with a median of 5 distinct full-length ASVs per genome (interquartile range of between 3 and 6). Across the 1,390 E. coli genomes in the rrnDB, there were 1,402 versions of the gene. These observations highlight the risk of selecting a threshold for defining clusters that is too narrow because it is possible to split a single genome into multiple clusters.

3A method to avoid splitting a single genome into multiple clusters is to cluster 16S rRNA gene sequences together based on their distances between each other. Therefore, I assessed the impact of the distance threshold used to define clusters of 16S rRNA genes on the propensity to split a genome into separate clusters. To control for uneven representation of genomes across species, I randomly selected one genome from each species and repeated each randomization 100 times. I observed that as the rrn copy number increased, the distance threshold required to reduce the ASVs in each genome to a single OTU increased (Fig. 1). Among species with 7 copies of the rrn operon (e.g., E. coli), a distance threshold of 5.25% was required to reduce full-length ASVs into a single OTU for 95% of the species. Similarly, thresholds of 5.25, 2.50, and 3.75% were required for the V3-V4, V4, and V4-V5 regions, respectively. But if a 3% distance threshold was used, then ASVs from genomes containing fewer than 6, 6, 8, and 6 copies of the rrn operon would reliably be clustered into a single OTU for ASVs from the V1-V9, V3-V4, V4, and V4-V5 regions, respectively. Consequently, these results demonstrate that broad thresholds must be used to avoid splitting different operons from the same genome into separate clusters.

4At broad thresholds, 16S rRNA gene sequences from multiple species could be clustered into the same ASV or OTU. I again randomly selected one genome from each species to control for uneven representation of genomes across species. For this analysis, I measured the percentages of ASVs and OTUs that contained 16S rRNA gene sequences from multiple species (Fig. 2). Without using distance-based thresholds, 4.1% of the ASVs contained sequences from multiple species when considering full-length sequences, and 10.9, 16.2, and 13.1% contained sequences from multiple species when considering the V3-V4, V4, and V4-V5 regions, respectively. At the commonly used 3% threshold for defining OTUs, 27.4% of the OTUs contained 16S rRNA gene sequences from multiple species when considering full-length sequences, and 31.7, 34.3, and 34.8% contained sequences from multiple species when considering the V3-V4, V4, and V4-V5 regions, respectively. Although the actual fractions of ASVs and OTUs that contain sequences from multiple species are dependent on the taxonomic composition of the sequences in the rrnDB, this analysis highlights the trade-offs of using distance-based thresholds.

5The results of this analysis demonstrate that there is a significant risk of splitting a single genome into multiple clusters if using ASVs or too fine of a threshold to define OTUs. An ongoing problem for amplicon-based studies is defining a meaningful taxonomic unit. Since there is no consensus for a biological definition of a bacterial species, microbiologists must accept that how bacterial species are named is biased and that taxonomic rules are not applied in a consistent manner (e.g., see references and ). This makes it impossible to fit a distance threshold that matches a set of species names. Furthermore, the 16S rRNA gene does not evolve at the same rate across all bacterial lineages, which limits the biological interpretation of a common OTU definition. A distance-based definition of a taxonomic unit based on 16S rRNA gene or full-genome sequences is operational and not necessarily grounded in biological theory. One benefit of a distance-based OTU definition is the ability to mask residual sequencing error. The analysis in this study was conducted using ideal sequences from assembled genomes, whereas sequences generated in microbiome studies would harbor PCR and sequencing errors. These errors would only exacerbate the inflated number of ASVs. There is general agreement in bacterial systematics that to classify an organism to a bacterial species, phenotypic and genome sequence data are needed. A short sequence from a bacterial genome simply cannot differentiate between species. Moreover, it is difficult to defend a clustering threshold that would split a single genome into multiple taxonomic units. It is not biologically plausible to entertain the possibility that different rrn operons from the same genome would have different ecologies. Individual bacteria are defined at the cellular or chromosomal level and not at the gene level. One could argue that, in practice, communities are compared on a relative rather than an absolute basis. However, communities harboring populations that tend to have more copies of the rrn operon would appear to have higher richness and diversity than those with fewer copies purely due to the propensity for populations with more rrn operons to generate more ASVs. Although there are multiple reasons why proponents favor ASVs, the significant risk of artificially splitting genomes into separate clusters is too high to warrant their use.

6The 16S rRNA gene sequences used in this study were obtained from the rrnDB (https://rrndb.umms.med.umich.edu) (version 5.7, released 18 January 2021). At the time of submission, this was the most current version of the database. The rrnDB obtained the curated 16S rRNA gene sequences from the KEGG database, which ultimately obtained them from the NCBI nonredundant RefSeq database. The rrnDB provided downloadable versions of the sequences with their taxonomy as determined using the naive Bayesian classifier trained on the RDP reference taxonomy. For some genomes, this resulted in multiple classifications since a genome’s 16S rRNA gene sequences were not identical. Instead, I mapped the RefSeq accession number for each genome in the database to obtain a single taxonomy for each genome. Because strain names were not consistently given to genomes across bacterial species, I disregarded the strain-level designations.

7The full-length 16S rRNA gene sequences were aligned to a SILVA reference alignment of the 16S rRNA gene (v. 138) using the mothur software package (v. 1.44.2). Regions of the 16S rRNA gene were selected because of their use in the microbial ecology literature. Full-length sequences corresponded to E. coli strain K-12 substrain MG1655 (GenBank accession number NC_000913) positions 28 through 1491, the V4 region corresponded to positions 534 through 786, V3-V4 corresponded to positions 358 through 786, and V4-V5 corresponded to positions 534 through 908. The positions between these coordinates reflect the fragments that would be amplified using commonly used PCR primers.

8Pairwise distances between sequences were calculated using the dist.seqs command from mothur. The OptiClust algorithm, as implemented in mothur, was used to assign 16S rRNA gene sequences to OTUs. Distance thresholds of between 0.25 and 10.00% in increments of 0.25 percentage points were used to assign sequences to OTUs.

9Because of the uneven distribution of genome sequences across species, I randomly selected one genome from each species for the analysis of splitting genomes and clustering ASVs from different species (Fig. 1 and 2). The random selection was repeated 100 times. Analyses based on this randomization reported the median of the 100 randomizations. The interquartile range between randomizations was less than 0.0024. Because the range was so small, the confidence intervals were narrower than the thickness of the lines in Fig. 1 and 2 and were not included.

10The code to perform the analysis in this article and its history are available as a git-based version control repository at GitHub (https://github.com/SchlossLab/Schloss_rrnAnalysis_mSphere_2021). The analysis can be regenerated using a GNU Make-based workflow that made use of built-in bash tools (v. 3.2.57), mothur (v. 1.44.2), and R (v. 4.1.0). Within R, I used the tidyverse (v. 1.3.1), data.table (v. 1.14.0), Rcpp (v. 1.0.6), furrr (v. 0.2.2), here (v. 1.0.1), and rmarkdown (v. 2.8) packages. The conception and development of this analysis are available as a playlist on the Riffomonas YouTube channel (https://youtube.com/playlist?list=PLmNrK_nkqBpL7m_tyWdQgdyurerttCsPY).

11I used “ASV” to denote the cluster of true 16S rRNA gene sequences that were identical to each other and “OTU” to denote the product of distance-based clustering of sequences. Although ASVs represent a type of operational definition of a taxonomic unit and can be thought of as an OTU formed using a distance of zero, proponents of the ASV approach prefer to avoid the term OTU given the long history of OTUs being formed by distance-based clustering (https://github.com/benjjneb/dada2/issues/62 [accessed 26 February 2021]). For this reason, when an ASV split a genome into different units, those units were called clusters rather than OTUs.