A reference cytochrome c oxidase subunit I database curated for hierarchical classification of arthropod metabarcoding data

Introduction

With the increasing availability of high-throughput DNA sequencing, scientists with a wide diversity of backgrounds and interests are increasingly utilizing this technology to achieve a variety of goals. One growing area of interest involves the use of metabarcoding, or amplicon sequencing, for biomonitoring, biodiversity assessment and community composition inference (Yu et al., 2012; Guardiola et al., 2015; Richardson et al., 2015). Using universal primers designed to amplify conserved genomic regions across a broad diversity of taxonomic groups of interest, researchers are afforded the opportunity to survey biological communities at previously unprecedented scales. While such advancements hold great promise for improving our knowledge of the biological world, they also represent new challenges to the scientific community.

Given that bioinformatic methods for taxonomic inference of metabarcoding sequence data are relatively new, the development, validation and refinement of appropriate analytical methods is ongoing. Relatively few studies have characterized the strengths and weaknesses of different bioinformatic sequence classification protocols (Porter & Golding, 2012; Bengtsson-Palme et al., 2015; Peabody et al., 2015; Somervuo et al., 2016; Richardson, Bengtsson-Palme & Johnson, 2017). Further, researchers continue to utilize a diversity of methods to draw taxonomic inferences from amplicon sequence data. Relative to alignment-based nearest-neighbor and lowest common ancestor-type classification approaches, methods involving hierarchical classification of DNA sequences are popular as they are often designed to estimate the probabilistic confidence of taxonomic inferences at each taxonomic rank. However, studies explicitly examining the accuracy of classification confidence estimates are rare (Somervuo et al., 2016).

When performing hierarchical classification, the construction, curation and uniform taxonomic annotation of the reference sequence database is an important methodological consideration. Database quality can affect classification performance in numerous ways. For example, artifacts within the taxonomic identifiers of a reference database can represent artificial diversity and the inclusion of sequence data adjacent to the exact barcoding locus of interest likely display sequence composition that is unrepresentative of the barcoding locus. Lastly, sequence redundancy within reference databases increases computational resource use and is particularly problematic for classification software programs that classify sequences based on a set number of top alignments. In general, such database artifacts have the potential to bias model selection and confidence estimation both with k-mer style classifiers such as UTAX, SINTAX and the RDP Naïve Bayesian Classifier (Wang et al., 2007; Edgar, 2015; Edgar, 2016) and alignment-based classification approaches such as Metaxa2 and Megan (Huson et al., 2011; Bengtsson-Palme et al., 2015). Thus, it is important to identify and manage reference sequence database artifacts during curation for optimal downstream classification performance.

The use of molecular barcoding and metabarcoding in arthropod community assessment and gut content analysis has gained popularity in recent years (Corse et al., 2010; Yu et al., 2012; Mollot et al., 2014; Elbrecht & Leese, 2017). However, as with other non-microbial taxonomic groups of interest, few researchers have developed hierarchical DNA sequence classification techniques for arthropods (Porter et al., 2014; Tang et al., 2015; Somervuo et al., 2017). Here, we detail the construction, curation and evaluation of a database designed for hierarchical classification of amplicon sequences belonging to a 157 bp COI locus commonly used for arthropod metabarcoding (Zeale et al., 2011). This work will serve as both a resource for those conducting experiments using arthropod metabarcoding and as a template for future work curating and evaluating hierarchical sequence classification databases.

Methods

Results

Following curation and extraction, we obtained 199,206 reference amplicon sequences belonging to 51,416 arthropod species. Over 90 percent of the references were between 142 and 149 bp in length, with a minimum reference sequence length of 94 bp. For the final database creation and classifier training procedure, many reference amplicons were shorter than the 157 bp region of interest due to the incompleteness of some reference sequences and the trimming of taxonomically uninformative ends during Metaxa2 training. Prior to this step, 82 percent of the sequences were between 150 and 157 bp in length following the original extraction and these longer sequences can be found at the GitHub repository associated with this work. The taxonomic representativeness of the database across different arthropod classes and insect orders, including the number of families, genera and species in each, are presented in Tables 2 and 3.

Analyzing the number of sequences per species in the input reference sequence data, we observed a heavily right-skewed distribution, with a median of 2 and a mean of 11.1 sequences per species (Fig. 1A). Further, 32.0 percent of species were represented by 5 or more sequences and 40 species, including Bemisia tabaci and Delia platura, were represented by between 1,000 and 9,736 entries. After conducting amplicon sequence extraction using the Metaxa2 database builder tool, we were able to extract the COI region of interest from 80.8 percent of the input sequences, which corresponded to 73.4 percent of the genus-level diversity found in the original input data. Following sequence dereplication, removal of sequences with three or more ambiguous base calls and taxonomic lineage curation, our final database contained approximately 13 percent of the input extracted sequences, which represented 98.2 percent of the genus-level richness of the input extracted reference amplicon sequences (Figs. 1B and 1C).

Regressing classification outcome against the Metaxa2 reliability score yielded a similar best fit model for both the 80 bp and full length test sequence datasets (Fig. 2). For both regressions, the probability of sequence mis-assignment was below 10 percent for reliability scores above 70. Thus, for our evaluations, we chose a reliability score of 68, which corresponded to family-level error probabilities of approximately 11.3 percent and 9.5 percent for 80 bp and full length sequences, respectively.

In evaluating the performance of our classification database when analyzed with Metaxa2 at a reliability score cutoff of 68, we found consistently low proportions of false positives across all ranks, though the proportion of true positive, true negative and false negatives varied more considerably from kingdom to species (Fig. 3A). Further, while the Metaxa2 false discovery rate increased with higher resolution ranks (Fig. 3B), it was generally low, never exceeding 5 percent at the genus level. Interestingly, Metaxa2 displayed low variance in the proportion of sequences assigned when classifying 80 bp sequences relative to full length sequences of 147 bp in median length. Overall, the proportion of sequences assigned was greater than 90 percent through the order level for both full length and half length sequences. Beyond the order level, this statistic decreased to 53 and 56 percent at the species level for half length and full length sequences, respectively. Conversely, the proportion of false positives varied more strongly by sequence length and was greatest at higher-resolution taxonomic levels. At the species level, 1.97 percent of 80 bp sequences were misclassified, compared to only 1.13 percent for full length sequences. At the order level, the percent of sequences misclassified was 0.59 and 0.65 for 80 bp and full length sequences, respectively. As measured in errors per assignment, the classification false discovery rate was similarly highest at the species level, with 7.3 and 6.3 percent of assignments being incorrect for 80 bp and full length sequences, respectively. False discovery rates decreased to 1.2 and 0.7 percent of assignments being incorrect at the order level for 80 bp and full length sequences, respectively.

During our evaluation of taxonomic overclassifiction, we found between 3,141 and 3,202 sequence test cases belonging to species not represented in the corresponding training data and between 612 and 630 sequence test cases belonging to genera not represented in the corresponding training data across the three iterations of training and testing data. At the species level, the proportion of these cases which were overclassified was roughly equal for full length and half length sequences, with 5.4 and 5.2 percent being overclassified (Fig. 4A). With respect to genus level classification performance on these species overclassification test cases, approximately 31 and 37 percent of test cases were classified correctly as true positives or true negatives for half length and full length sequences. Genus level false positive proportions for these sequence cases were 6.2 and 5.1 percent (Fig. 4B). For the genus level overclassification cases, the difference in overclassification rates by sequence length was slightly larger, with approximately 8.9 and 10.9 percent of full length and half length test cases being overclassified (Fig. 4C). Family level performance on genus level overclassification cases was slightly lower with the proportion of true positive and true negative identifications summing to 25 and 27 percent for half length and full length sequences, with corresponding false positive proportions of 8.8 and 8.0

Evaluating classification performance for each order resulted in unsurprising outcomes across taxonomic ranks, from family to species (Fig. 5). Generally, the proportion of sequences assigned decreased and false discovery rate increased with increasing taxonomic resolution. Between orders, there was noteworthy variation in performance. For example, the proportion of sequences assigned was highest among trichopteran sequences and lowest among lepidopteran sequences, with approximately 94 and 50 percent of sequences being assigned to the family rank for each order, respectively. Further, while the proportion of sequences assigned was similarly lowest for lepidopteran sequences at the genus and species levels, wherein 39 and 21 percent of sequences were assigned for each respective rank, the highest proportion of sequences assigned at these ranks was not observed with trichopteran sequences. Instead, genus and species level proportion assigned was highest among Sessilia sequences, with 86 and 78 percent of sequences being assigned to each respective rank.

Discussion

While species-specific PCR and immunohistochemistry-based methods have been useful in documenting arthropod food webs (Stuart & Greenstone, 1990; Symondson, 2002; Weber et al., 2006; Blubaugh et al., 2016), the narrow species-by-species nature of such approaches has limited their utility for answering large-scale or open-ended ecological questions. With the increasing availability of high-throughput sequencing, arthropod metabarcoding will continue to become more broadly applicable to scientific questions spanning a diversity of research areas. The development of improved methods for drawing maximal inferences from sequence data is an important area for further methodological research. In creating a highly curated COI reference amplicon sequence database and evaluating its performance when used with the Metaxa2 taxonomic classifier, we have developed a new method to aid researchers in the analysis of arthropod metabarcoding data.

Though predictions vary greatly, researchers have estimated the species richness of arthropods to be between 2.5 to 3.7 million (Hamilton et al., 2010). Further, according to the literature review of Porter et al. (2014), 72,618 insect genera have been described to date. Thus, the 51,416 species and 17,039 genera represented in our database account for only a small fraction of arthropod biodiversity. The limited representativeness of currently available, high quality reference sequence amplicons for the COI region highlights the need for continued efforts to catalogue arthropod biodiversity with molecular techniques. Despite this current limitation, the combination of molecular gut content analysis with high-throughput sequencing is a promising path toward investigating arthropod trophic ecology and biodiversity monitoring with greater sensitivity and accuracy relative to alternate approaches.

The results of our inventory of sequences per species and genus level richness at various stages in the database curation process revealed that our amplicon extraction procedure was highly sensitive, extracting and trimming approximately 81 percent of the input sequences down to the 157 bp region of interest. Further, approximately 87 percent of these extracted sequences represented sequence redundancies and were removed during dereplication. As mentioned previously, the trimming of sequence residues adjacent to the barcode of interest and removal of redundant sequences not only makes computational analysis less resource intensive, it can also improve classification performance. For k-mer style classifiers, extraneous sequence residues can bias model selection during classifier training, while abundant sequence duplicates can result in an overwhelming number of identical top hit alignments for alignment-based classifiers.

Overall, the best fit local regression models summarizing the relationship between the Metaxa2 reliability score and the probability of classification error were useful in that the likelihood of misclassification was always less than what would be expected based on the reliability score. For example, a reliability score of 90 corresponded to only a 3.3 percent probability of family-level misclassification for full length sequences. We selected a reliability score of 68 for subsequent analysis as this provided a balanced trade-off between sensitivity and accuracy. Using this reliability score we observed minimal false positive rates and overall proportions of misclassification when comparing our results to those of similar studies (Porter et al., 2014; Bengtsson-Palme et al., 2015; Edgar, 2016; Richardson, Bengtsson-Palme & Johnson, 2017). Given that the family level probability of error was only 9.5 percent at a reliability score of 68, a lower reliability score threshold may be justifiable for certain research situations. However, further testing should be conducted to ensure that the relationship between reliability score and classification confidence is similar across taxonomic ranks and between different DNA barcodes.

With respect to sensitivity using a reliability score of 68, our results were highly dependent upon the rank being analyzed, with sensitivities, as measured by the total proportion of sequences classified, above 60 percent only being achieved at the family and order ranks. To some degree, these sensitivity estimates reflect the large degree of database incompleteness at the genus and species ranks, wherein approximately 3.6 and 25.5 percent of unclassified sequences were true negatives. However, to our knowledge, no other studies have reported classification sensitivity data for this COI amplicon locus. This makes it difficult to ascertain if Metaxa2 –with an average of approximately 44.7 percent of false negative assignments at the genus and species ranks –exhibits relatively low sensitivity or if this locus is limited in discriminatory power. The short length of the 157 bp COI amplicon region relative to other barcoding regions such as the ITS and 18S rRNA regions (Hugerth et al., 2014; Wang et al., 2015) could be a cause of such limited discriminatory power.

As expected, analyzing cases of overclassification in our data revealed that sequences from taxa lacking representation in the database are far more likely to be misclassified relative to sequences from well-represented taxa. This is supported by an approximately 10 percent probability of genus level overclassification for sequences from unrepresented genera relative to a 1 to 2 percent probability, depending on sequence length, for all sequence test cases. Interestingly, the genus level overclassification rate was approximately double that observed at the species level. This seems counterintuitive but is expected in light of the discussion put forth by Edgar (2018), wherein the percent identity difference from the closest reference sequence match is considered one of the most important predictors of classification performance. For a species-level overclassification case, the closest reference sequence match in the corresponding database is a sequence from a congeneric species in the best case scenario. Since the closest reference sequence match for a genus level overclassification case is a sequence from a confamilial species at best, the average percent identity difference from the closest reference is greater for a genus level overclassification case than for a species level overclassification case. Thus, for overclassification cases in particular, higher levels of error are expected at more inclusive taxonomic ranks.

While the genus level overclassification estimates we observed are not desirable, they are considerably lower than similar estimates for the RDP classifier, which range from 21.3 percent to 67.8 percent depending on the database, locus analyzed and cross-validation approach used (Edgar, 2016; Richardson, Bengtsson-Palme & Johnson, 2017). Further, the observed degree of genus level overclassification using Metaxa2 with our COI database was similar to or less than that of the recently developed SINTAX classifier (Edgar, 2016; Edgar, 2018). Interestingly, the recent analysis of Edgar (2018), resulted in a similar estimate of the Metaxa2 genus level overclassification rate while revealing a considerably higher overclassification rate for the majority of other classifiers tested. Though this analysis found Metaxa2 to be relatively less sensitive than other classifiers, a weakness of the work was that it only tested Metaxa2 using the default reliability score of 80, while testing multiple confidence thresholds for alternate classifiers, such as SINTAX and RDP. Given that our analysis has revealed the Metaxa2 default reliability score to be too conservative—at least for the COI locus—such results are difficult to interpret. In general, such comparisons across studies should be approached with caution as multiple factors complicate the interpretation of classification performance, such as locus discriminatory power, database completeness and the choice of evaluation metrics used. Ultimately, direct comparisons of classification methods using standardized loci and databases are needed to more rigorously compare performance.

With respect to Metaxa2 classification of full length relative to half length amplicon sequences, we observed surprisingly small differences in performance. Consistently, the proportion of misclassified sequences was greater for half length sequences. When considering error and sensitivity together, the false discovery rate or errors per assignment for full length sequences was consistently less than or equal to that achieved during the classification of half length sequences. Lastly, when considering the relationship between the Metaxa2 reliability score and the probability of classification error at the family level, we noted highly similar local polynomial regression models of error probability for both full length and half length sequences.

Conclusions

Here, we assembled a highly curated database of arthropod COI reference amplicon sequences, trained a recently developed hierarchical DNA sequence classifier using the database and conducted extensive in silico performance evaluations on the resulting classification pipeline. Overall, we found a high degree of sequence redundancy within the initial, uncurated dataset, highlighting the importance of effective sequence dereplication during the creation of databases designed for metabarcoding analysis. Further, the limited representativeness of the database with respect to arthropod biodiversity indicates that additional sequencing effort is needed to further improve the performance of arthropod metabarcoding techniques. Though the performance evaluations presented in this work were conducted on a large corpus of available biological data, the results are not necessarily directly transferable to all experimental settings. For example, variations in sequence error profiles and taxonomic distributions among datasets are potential confounding factors. Despite these limitations, this work provides researchers with a new resource for arthropod COI sequence analysis and novel data for gauging the strengths and limitations of different approaches to arthropod metabarcoding.