Microbial community structure is important to characterize because it can influence many ecosystem processes (Graham et al., 2016). For assessments of overall microbial community composition, the 16S rRNA gene is typically used because it is highly conserved across prokaryotes and generally not subject to horizontal gene transfer. Given the ubiquity of 16S rRNA assessments within microbiome research, the bioinformatics pipelines used to process raw amplicon sequence data were developed with a focus on 16S rRNA gene sequence data specifically. Over the last decade, the approaches used to preprocess sequences, cluster unique sequences into OTUs, and assign taxonomic classification have been continually expanding and improving. However, side-by-side comparisons of 16S rRNA datasets resulting from pipelines which differ in only one or two key steps have demonstrated that upstream processing decisions (e.g., clustering of sequences) can influence conclusions about differential abundance, composition, taxonomic identity, and richness and diversity measures (Chen et al., 2013; Nguyen et al., 2016; López-García et al., 2018).

As a complement to overall community composition, studies focusing on specific microbially mediated ecological processes often characterize functional groups relevant to the process being investigated. These functional groups fall along a diversity spectrum ranging from processes which are performed by a comparatively narrow selection of taxa, or phylogenetically “narrow” processes, to those which can be performed by a large variety of taxa, or “broad” processes (Schimel and Schaeffer, 2012). The same bioinformatics tools developed for 16S rRNA gene sequences are used to process diagnostic gene sequences for these types of functional groups as well. Though efforts have been made to identify discrepancies among results generated from different processing methods using 16S rRNA gene sequence datasets, the impact of processing choices on downstream community analyses for functional genes of varying diversity has not yet been explored.

For both 16S rRNA and functional genes, amplicon sequence data are often generated via Illumina sequencing in the form of millions of paired-end reads. Generally, primers for amplicon sequencing are designed to generate forward and reverse reads which overlap and can be assembled into continuous sequences, or contigs. The first step in processing raw sequence data for downstream OTU clustering is therefore merging the two reads and filtering the resulting contigs for quality control. Modern sequencing instruments return a quality score for each base in each sequence, and this quality score, together with agreement between bases in the overlapping portion of paired reads, is used to filter out low-quality contigs. These steps can be achieved using many common bioinformatics software packages, including Usearch (Edgar, 2013), FLASH (Magoè and Salzberg, 2011), Mothur (Schloss et al., 2009), or QIIME (Caporaso et al., 2010). The implementation of these steps is similar among all packages and allows the user to provide arguments to customize quality cutoffs as desired. The end result of this stage of preprocessing is a list of all unique sequences which passed quality filtering.

After quality screening, unique sequences are assigned to an operational taxonomic unit (OTU), which is most often achieved by clustering sequences according to some similarity percentage. The purpose of clustering is twofold: clustering sequences together by similarity helps to eliminate erroneous sequences formed during the PCR preamplification step carried out prior to sequencing, as each of these erroneous sequences should deviate from one another by only a few bases, thus reducing diversity to true biological diversity (Hugerth and Andersson, 2017). In addition, collapsing the full breadth of sequence diversity into groups within some percentage similarity of one another reduces the total number of “variables” in downstream analyses, making them more computationally tractable.

The similarity threshold typically chosen is 97% (Gevers et al., 2005), as similarity levels lower than this in the 16S rRNA gene region are considered unlikely to be derived from the same species and unlikely to achieve 70% DNA-DNA hybridization at the genome level, a previously common metric for determining bacterial species assignment (Gevers et al., 2005). This, however, can only evaluate the gene region considered, and does not necessarily reflect 97% similarity across the full length of the gene. This is also an arbitrary cutoff, as individual taxa may possess 16S rRNA genes that are more than 97% similar but still represent ecologically distinct clades based on the remainder of their genome content (Fox et al., 1992; Gevers et al., 2005). As sequencing throughput and quality has increased in recent years, a 99% similarity cutoff for inclusion in an OTU has become increasingly common (Hugerth and Andersson, 2017). Acceptable percent similarity cutoffs for OTUs generated from functional gene sequences have not yet been established, and 97% is still typically used, regardless of the diversity within the gene’s phylogeny.

Clustering approaches can be reference-based or de novo. The former uses a reference taxonomic database to classify sequences into taxonomic bins based on known taxonomy, while the latter allows the data to “speak for themselves” by assigning sequences to clusters based on similarity alone (Schloss and Westcott, 2011). Reference-based clustering can be either closed reference, wherein sequences are mapped to their best possible match within a database, and those that do not match sufficiently are discarded, or open reference, where those sequences which do not match to the reference are then clustered using the de novo approach. In de novo clustering, approaches may be hierarchical (based on single, average, or complete linkage) (Schloss and Handelsman, 2005) or heuristic in strategy. Single-linkage hierarchical approaches place a sequence into a cluster if it has a similarity above some threshold to at least one other sequence in the cluster, while complete linkage conversely requires a sequence to have a similarity above the threshold to all others in the cluster; average linkage requires that the average similarity between a sequence and all others be above the threshold (Hugerth and Andersson, 2017). Average-linkage de novo clustering has been demonstrated to produce higher quality OTUs based on the Matthew’s Correlation Coefficient (MCC), a metric for describing the ratio of False Positives (FP), False Negatives (FN), True Positives (TP), and True Negatives (TN) commonly used in machine learning control theory (Westcott and Schloss, 2015). However, because it is computationally expensive to run all-against-all comparisons on datasets containing millions of reads, heuristic approaches were developed. Among these are the UPARSE algorithm implemented via USEARCH (Edgar, 2010), which approximates average-linkage approaches by comparing a sequence to only one centroid sequence within each cluster. The choice between true hierarchical clustering and heuristic clustering therefore represents a tradeoff between computational speed and accuracy.

The distribution of distances between sequences in clusters will differ depending on the clustering approach used, even if using the same similarity cutoff (Hugerth and Andersson, 2017). Among commonly used modern tools, hierarchical clustering is available through the “Opti-MCC” method implemented in Mothur (Westcott and Schloss, 2017), which is now included as the default clustering method. Heuristic approaches are available through USEARCH (Edgar, 2013), QIIME (Caporaso et al., 2010), which runs Uclust in the background, and VSEARCH (Rognes et al., 2016), an open-source alternative to USEARCH. Among these heuristic options, abundance-based greedy clustering (AGC) is often the default implementation. The AGC clustering algorithm begins with the most abundant unique sequences in the dataset and begins building OTUs from these. This operates under the assumption that the most abundant sequences are more likely to be biologically “accurate,” and do not represent sequencing errors. That the AGC method is “greedy” means that as it works through the list of sequences, it places a sequence with the first match it finds that meets the percent similarity threshold—it does not continue looking to see if there is an even closer match, and once a decision is made, it cannot be changed afterward. The implication of this is that resulting clusters are less accurate, and thus often leads to fewer OTUs and more dissimilar sequences assigned to each OTU, despite taking far less time to compute. In addition, this approach may introduce spurious correlations between samples or treatments which are overrepresented in the most abundant sequences. Conversely, Opti-MCC, a hierarchical clustering method, uses an iterative approach which repeatedly reevaluates the clusters formed until the MCC (ratio of FP, FN, TP, TN) (Westcott and Schloss, 2015) is optimized. This method begins with each unique sequence as its own OTU, and checks whether combining each pair of OTUs will improve the MCC—if it does, they are combined. This progresses until no further combinations remain that will improve the MCC. Unsurprisingly, this iterative approach requires much more time to execute, but the results are optimized clusters which more accurately group sequences according to percent similarity. Therefore, the choice of clustering approach can affect downstream analyses and conclusions regarding microbial community diversity and structure due to its direct effect on assignment of sequences to OTUs.

To further complicate matters, many common multivariate analyses used in microbiome studies are sensitive to uneven count data (unequal number of sequence reads per sample), thus requiring normalization prior to analysis. Historically, this has been achieved through rarefaction, which involves randomly subsampling each sample’s reads down to an even depth. This, however, has been demonstrated to reduce statistical power, and may dramatically change the conclusions drawn from downstream multivariate analyses like PERMANOVA (McMurdie and Holmes, 2014). Comparing analytical results from repeated rarefying trials has been previously suggested (Navas-Molina et al., 2013), but this is time consuming and is not generally practiced. Therefore, most published results are obtained from analyses performed on a rarefied OTU table which was produced by randomly subsampling the original OTU table a single time. When the entire burden of proof for downstream analyses rests on this single subsample, consistency among random samples becomes critically important. Otherwise, we allow chance to determine whether or not our subsample contains enough statistical power to reject the null hypothesis. For differential abundance analyses focusing on individual OTUs, we are able to sidestep these pitfalls of rarefaction by using the negative binomial mixed model implementation in R package DESeq2 (Anders and Huber, 2010); however, this approach still depends on the accuracy of the upstream OTU clustering. Ultimately, regardless of the analysis methods used downstream, the accuracy of OTU assignment will play a role in interpreting DNA sequence data.

These major considerations and shortcomings related to OTU clustering for 16S rRNA gene sequences, the most commonly sequenced gene for which all of these methods were developed, do not even begin to address methodological considerations for the quagmire of diversity found within common functional genes of interest. For example, the suite of functional genes commonly evaluated in microbiome studies concerning nitrogen (N) cycling range from phylogenetically narrow (e.g., N-fixation, nitrification) to broad [e.g., denitrification, dissimilatory nitrate reduction to ammonium (DNRA)]. However, a method which can accurately characterize sequences harboring comparatively little diversity (which can be binned into fewer OTUs) may not accurately resolve the diversity within a gene comprised of a large number of clades. The diagnostic gene for N-fixation, nifH, segregates into only 3 distinct clades (Raymond et al., 2004), and thus an algorithm like AGC can be expected perform sufficiently because each of these clades are likely to be represented within the most abundant sequences. In contrast, the diagnostic gene for DNRA, nrfA, is comprised of 18 distinct clades (Welsh et al., 2014). In this case, an abundance-based method which preferentially assigns sequences to the largest clusters runs the risk of pulling sequences that might otherwise represent their own smaller clusters into the larger initial clusters. Because AGC is also “greedy,” this means that there is no reassessment afterward to reconcile this error. A true hierarchical clustering method, however, will be more likely to detect these smaller clusters of similarity and assign them to their own OTU as clustering proceeds. Since functional gene community analyses are typically presented together to illustrate treatment effects on specific functional groups, these differences in accuracy may introduce biases in the conclusions that may be drawn from the resulting OTU datasets.

In this study, we evaluated the performance of two clustering methods, one heuristic (AGC) and one hierarchical (Opti-MCC), on a functional gene comprised of few clades (nifH) and a gene comprised of many clades (nrfA). We predicted that both the AGC and Opti-MCC methods would perform similarly on nifH sequences, but that Opti-MCC would produce more statistical power than AGC when applied to nrfA sequences due to its greater ability to characterize the diversity within the gene. We processed the same two sets of raw Illumina sequences in Mothur, once using the Opti-MCC clustering algorithm at the standard 97% similarity cutoff, and twice using the AGC algorithm: Once at 97% similarity, and again at 98%, to evaluate if tightening the similarity threshold was able to better resolve greater diversity. Each resulting OTU table was then subjected to 10 repeated rarefaction trials, which generated 10 rarefied OTU tables from each raw OTU table. The ability to detect community differences using common multivariate methods was assessed for each rarefied dataset and results were aggregated for comparison between methods. Alpha diversity metrics were also assessed for each unrarefied and rarefied dataset. Finally, differential abundance analyses were conducted on each original unrarefied OTU table to evaluate performance independent of rarefaction, as well as performance in taxonomic classification. Based on these comparisons, we conclude with recommendations about appropriate use cases for each method, and when priority should be placed on clustering accuracy vs. computational speed.

The processing necessary to convert raw Illumina amplicon sequence data to a usable OTU table can be achieved through any number of approaches within the continuously growing toolbox of bioinformatics methods. However, the choices made at each step during upstream processing are liable to impact downstream results, and the ways these choices influence different types of data are only beginning to be explored. In this study, we conducted a comparison of OTU tables generated from pipelines which differed only in their clustering method on each of two functional genes: nifH, the diagnostic gene for N-fixation, which is comprised of only 3 clades, and nrfA, the diagnostic gene for DNRA, which is comprised of at least 18 clades. AGC, a quicker heuristic method, performed comparably to the slower, more accurate hierarchical Opti-MCC method for the phylogenetically narrow gene, nifH. For nrfA, however, which harbors a much greater diversity than does nifH, the AGC method did not perform as well as the Opti-MCC method, producing an OTU table which resulted in unreliable downstream community analyses post-rarefaction and evidence of reduced granularity in classifying OTU taxonomy. Our results demonstrate how clustering methods optimized for 16S rRNA phylogenies may perform differently depending on the diversity and lineages of the functional gene sequences being processed.

Many of the differences observed between results from the two methods can be attributed to the strategy each uses to assign sequences to clusters. The Opti-MCC method initializes each unique sequence as its own OTU, and proceeds to combine these based on all-against-all similarity comparisons for the sequences in each cluster (Westcott and Schloss, 2017). The metric used to optimize cluster assignment, the MCC, represents a balance of not only true positives, but also false positives, true negatives, and false negatives. In contrast, the AGC method approximates this process by assigning sequences to clusters based on similarity to an averaged centroid sequence, beginning with the most abundant OTUs (He et al., 2015). As a “greedy” algorithm, it places a sequence with the first OTU it finds whose centroid sequence is within the chosen similarity cutoff. Since these comparisons begin by considering the most abundant OTUs, there is a risk of placing a sequence with an abundant OTU when it would be more accurately classified in a smaller OTU which was not considered for comparison. In other words, OTU tables generated using AGC may have a higher rate of false positives within the most abundant OTUs, and a higher rate of false negatives in rarer OTUs. This is manifested in the differing performance of the AGC method on genes of differing diversity, as the full breadth of diversity is more likely to be represented among the most abundant OTUs for a phylogenetically narrow gene compared to a broad one.

The clearest demonstration of this bias toward the most abundant OTUs in AGC can be seen in the differing level of consistency among PERMANOVA results for the full OTU table vs. only the 100 most abundant OTUs. While both AGC and Opti-MCC methods produced consistent PERMANOVA results for analyses of all OTUs and only the top 100 OTUs for nifH, they differed in performance for nrfA. The Opti-MCC method yielded reliable results using the full OTU table, but many pairwise differences became undetectable when analyzing only the top 100 OTUs. This indicates that there is enough explanatory power in the data beyond of the top 100 OTUs that differences become difficult to detect when the less abundant OTUs are excluded. For the AGC method, however, we observed the opposite: for many PERMANOVA results which were inconsistent among rarefaction trials when the full dataset was analyzed, results became more reliable when only the top 100 OTUs were considered. These results demonstrate that explanatory power is concentrated more heavily in the most abundant OTUs when clusters are assembled using AGC.

Although the Opti-MCC method failed to produce consistent PERMANOVA results using only the top 100 OTUs, the taxonomic classification within the top 100 OTUs was more granular than that achieved from AGC. This speaks to the accuracy of the original cluster assignment, as fewer false positives result in a more precise representative sequence which can be better resolved against a reference database. When the AGC method assigns sequences preferentially to abundant OTUs without searching for a better fit among less abundant OTUs, the result is “fuzzier” clusters with a wider variety of sequences, resulting in an increasingly generic representative sequence. This may lead to poorer granularity in taxonomic classification. While we can circumvent some of the drawbacks of AGC-clustered OTUs by focusing only on the most abundant OTUs, taxonomy assignments may still suffer.

In contrast to the variable results obtained using PERMANOVA, the results from the same datasets were strikingly consistent when analyzed using ANOSIM: across all clustering methods and rarefaction trials, ANOSIM results consistently agreed in both statistical significance and approximate effect size. ANOSIM is a rank-based omnibus test which detects differences in several aspects of underlying data structure (Anderson and Walsh, 2013). Therefore, the observed consistency in results among clustering methods indicates that the agricultural treatment differences present in the sequence data were preserved in both the AGC and Opti-MCC datasets—and still differ with the same quantifiable magnitude—but that these differences are simply being structured differently within the OTU tables. In the case of the AGC-constructed nrfA OTUs, this structure made it difficult to detect these treatment differences via PERMANOVA when assessing the full OTU table. In addition to consistent ANOSIM results, processed data from both methods shared many similarities in terms of alpha diversity metrics. Results generated from each method followed the same treatment patterns and exhibited comparable effect sizes, for both the unrarefied OTU tables as well as each independently rarefied OTU table. While the outcomes of some analyses may differ between methods, these results demonstrate that the two can still produce very similar results for others.

While tightening the similarity cutoff for the AGC method resulted in many more OTUs, it still did not produce the same consistency among rarefied analytical results as the Opti-MCC method. While some pairwise comparisons improved in reliability, e.g., two trials resulting in 4/10 significant differences detected at 97% similarity improved to 7/10 and 9/10 significant under 98% similarity, others did not improve. In addition, the large number of OTUs rendered detection of differences within the top 100 OTUs impossible, and consideration of only the top 100 OTUs was no longer appropriate for determining differential abundance among treatments. Therefore, although using a higher similarity threshold marginally improved performance, it introduced additional issues, such as the need to reevaluate and identify appropriate cutoffs for differential abundance analyses. This clearly demonstrates the importance of prioritizing quality of OTU cluster formation over quantity of OTUs, as more is not necessarily better.

This study aimed to explicitly compare the impact of two OTU clustering algorithms on downstream analyses while holding constant all other aspects of the bioinformatics pipeline. While this approach allowed us to directly quantify the impact of clustering algorithm alone, this also limits the comparison to clustering approaches that may be implemented using the same pipeline and otherwise identical steps. However, an increasing number of researchers are beginning to turn to the use of amplicon sequence variants (ASVs) in lieu of OTUs, an approach which “denoises” the unprocessed sequence reads by clustering them into biologically meaningful groups independently of a predefined level of similarity (Hugerth and Andersson, 2017). A popular implementation of this approach is through the package DADA2 (Callahan et al., 2016), which begins by initializing clusters based on amplicon abundance (where, as in the AGC heuristic, more common sequences are assumed to be more biologically accurate) and sequence distance from other reads. The quality scores assigned to each base by the sequencing platform are then used to build an error model to “correct” reads by assigning low frequency reads to higher frequency reads from which they may have been derived via sequencing errors (Hugerth and Andersson, 2017). While this approach circumvents some of the pitfalls of OTU clustering at a fixed threshold, it is important to note that such denoising algorithms are still performing clustering operations and many of the same considerations regarding how performance may vary between functional genes of contrasting diversity may still apply.

Functional genes are increasingly assessed in microbiome research because they can shed light on the ways that our experimental treatments impact different functional group communities and the ecological processes they perform. However, in some cases, these functional genes represent a vast range of diversity and divergent lineages, with some comprised of only a few clades, and others an order of magnitude more. This study demonstrated some of the differences that can manifest when applying differing clustering algorithms to sequences from genes on opposite ends of the diversity spectrum. While the AGC method performed well on a gene with little diversity, it resulted in unreliable analytical outcomes for some tests when applied to a more diverse gene. This is especially problematic when we consider that functional gene community analyses are typically presented together as a means to demonstrate which functional groups are differentially impacted by a treatment. If we choose a method which obstructs detection of treatment differences in diverse genes, we introduce biased Type II errors by failing to identify treatment effects in some genes but not in others. Surprisingly, despite rarefying to a very low read depth for nifH, Type II errors were not introduced, and treatment differences were detected consistently among methods. This emphasizes the importance of considering gene characteristics when selecting methods, as more Type II errors were introduced simply by clustering nrfA sequences with the AGC method than were introduced by rarefying to a low read depth for nifH.

Ultimately, the reliability of our downstream analyses and subsequent conclusions is only as good as our upstream processing. For questions which can be answered using tests agnostic to clustering method (e.g., ANOSIM), or for genes of relatively low phylogenetic diversity (e.g., nifH), most upstream processing methods should lead to similar conclusions from downstream analyses. For studies involving more diverse genes, however, care should be exercised to ensure accurate clustering for all genes. Most importantly, we must continually reassess the performance of our preferred bioinformatics tools as technology continues advancing and more sophisticated methods become available.

Amplicon sequence data for 16S rRNA genes and N-cycling functional genes are available for download on the NCBI SRA database via the BioProject accession number: PRJNA752786 (https://www.ncbi.nlm.nih.gov/sra/PRJNA752786).

SE conceived of the study, conducted the study, and wrote the manuscript. RS, AK, and WY provided guidance and edited the manuscript. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.