How reliable is metabarcoding for pollen identification? An evaluation of different taxonomic assignment strategies by cross-validation

Barcodes: rbcL vs. ITS2

We chose to work with databases dedicated to the ITS2 and rbcL barcodes because they are commonly used and recommended for plant/pollen metabarcoding (e.g.: Bell et al., 2022; Lowe et al., 2022). All available sequences were retrieved from the National Center for Biotechnology Information (NCBI) and cleaned with the DB4Q2 pipeline (Dubois et al., 2022). In summary, sequence and taxonomy data was retrieved from NCBI and reformatted, low-quality sequences were discarded, a step of dereplication was carried out to reduce redundancy, two additional filters were applied to remove sequences suspected to have a fungal origin or to be misidentified, and a last optional step consisted in the database restriction, i.e. extracting from reference sequences only the portion amplified by a primer set of interest. The resulting databases cover a minimum of 81% (ITS2) to 87% (rbcL) of the species present in the wild in our study case area (Belgium) and approximately 95% of the genus.

Region of the barcode gene considered: General, Restricted, Restricted-General

Since the reference sequences were downloaded from NCBI without a stringent length filter, they displayed variable lengths: inter-quartile range between 547 and 688 pb for ITS2 (min = 109 pb, max = 14,531 pb) and between 551 and 1,331 pb for rbcL (min = 100 pb, max = 2,550 pb). However, with a metabarcoding approach, only a fraction of these sequences is amplified depending on the chosen primers. For example, the primer sets ITS-S2F (5′-ATGCGATACTTGGTGTGAAT-3′)/ITS4R (5′-TCCTCCGCTTATTGATATGC-3′) and rbcLaF (5′-ATGTCACCACAAACAGAGACTAAAGC-3′)/rbcLr506 (5′-AGGGGACGACCATACTTGTTCA-3′) that have been used in our case study, amplify fragments with median lengths of 343 and 458 pb respectively, with little variation.

Using the General DB containing the full-length sequences could artificially increase the accuracy of the taxonomic predictions. Indeed, with longer sequences more information is available to distinguish closely related species, whereas with a metabarcoding approach the typical amplicon length is rather short. Alternative databases were therefore built for each barcode restricted to the region amplified by the ITS-S2F/ITS4R and rbcLaF/rbcLr506primers (Restricted DB). The restricted databases contained, however, a much lower number of sequences: $\sim$22,000 and $\sim$15,000 sequences for ITS2 and rbcL respectively while the general databases contained $\sim$170,000 and $\sim$135,000 sequences respectively. Also, the restricted DB were not an exact subsample of the general DB because the clipping step appears before the cleaning step in our pipeline. In some cases, the full sequence was discarded because it did not pass the quality checks, while the sequence restricted to the amplified region passed the quality check and was conserved in the database. The drop in the number of sequences can be explained by two main reasons: (i) the full-length sequences only partially overlap the target region and/or (ii) the primer hybridization site of the full-length sequence has too many SNPs compared to the primer sequences provided for the database clipping.

In these two approaches, the Restricted or General DB were cross-validated against themselves. A third approach, “Restricted-General”, was also tested, in which the cross-validation involved blasting the restricted sequences against the general databases.

Species origin for the target sequence: Foreign sp. vs. Local sp.

To evaluate how the taxonomic assignment performs for species of different origin, the target sequences were classified into two categories: Local species vs. Foreign species. The Local sp. category is considering the flora susceptible to be found in our area of interest (Belgium) in a wide sense ( $\sim$2,000 species): it contains native species but also crops and non-native plants present in the wild. All other species from the the database were considered as “Foreign species” even if they might be planted in urbanized areas. These exotic plants can, however, be an important source of pollen for bees (Casanelles-Abella et al., 2022) and are of particular interest in our study case.

Database geographic area: Local DB vs. World DB

Reference databases containing only the local species (Local DB) were created, to compare with the result obtained when all sequences were kept in the database, i.e. corresponding to species from the entire world (World DB). We can then blast Local sp. sequences against the Local DB or World DB and compare the quality of taxonomic assignments obtained when we blast Foreign sp. against the same databases.

To be clear: the Local species and Foreign species labels concern the sequences that are being blasted, while Local DB and World DB concern the databases against which the sequences are blasted. The Local DB contains all the sequences of the Local sp. while the World DB contains the Local sp. + Foreign sp. sequences.

Cross-validation strategy: 10-fold CV vs. leaked CV

We tested two different ways to perform the cross-validation (CV). The 10-fold CV is the most conventional approach: the sequences are divided into 10 groups (=folds), each containing a random sample of 10% of the sequences. The sequences from each fold are then blasted against the remaining 90% of the sequences from the other folds. With the other approach named here (leaked CV), the sequences were blasted against the full database without excluding the target sequences. These two approaches are testing different hypotheses. With leaked CV, we test the ability of blast to match the exact sequence when it is present in the database among many other sequences. With 10-fold CV we test how blast performs to find alternative sequences that are taxonomically close to the true target sequence.

Also, we did not perform a full cross-validation for all the sequences. We kept all sequences corresponding to Local species ( $\sim$32,000 sequences across ITS2 and rbcL databases) as well as a random sample of 10% of the sequences corresponding to “Foreign species” from each database ( $\sim$25,000 sequences across ITS2 and rbcL databases). We also excluded the sequences with imprecise identifications (i.e. not identified to the species level), doubtful names and hybrids taxa. In the end, $\sim$31,000 ITS2 and $\sim$25,000 rbcL sequences were blasted against the general databases and $\sim$4,800 ITS2 and $\sim$3,500 rbcL sequences against the restricted databases.

For the Restricted vs. general 10-fold CV approach, the restricted sequences were split into 10 folds containing 10% of the target restricted sequences. Before blasting each fold against the general database, the full-length sequences corresponding to the restricted sequences present in the fold were removed from the reference database.

Descriptive statistics

Blast taxonomic assignment methods: TopHit, TopHitPlus, TopN, TopNplus

Four different methods were used to assign a taxon from the top 10 blast hits, two are based mainly on the highest bit score and the two others on the consensus scores. Both the identity and consensus scores are computed for each of the four methods and for each taxonomic level because they can be useful to assess the quality of the taxonomic assignments.

• i)

  With TopHit, the taxon with the best Bit score is selected (default blast output if we keep only one match).

• ii)

  With TopHitPlus, the taxon with the best bit score is also retained, but the potential bit score ties are broken by choosing the taxon with the highest consensus score among the ties.

• iii)

  With TopN, the taxon with the highest consensus score among the top 10 hits is selected at each taxonomic level. With this approach, the chosen species might belong to a different genus than the taxon chosen at the genus level for example.

• iv)

  With TopNPlus, the sequences which do not fulfill some predefined requirements are first eliminated, then the consensus score is computed and the taxon with the highest consensus score is computed at each taxonomic level. In this study, we discarded the matched sequences with an alignment length < 100 base pairs and with an E score > 10⁻¹⁰. The identity score had also to be at least 97% for the species level, 90% for the genus level and 80% for the family level. Also the difference of identity score between the best hit and the matched sequence had to be <1% for the species level, 10% for the genus level and 20% for the family level (approach inspired by e.g. Milla et al. (2022) and https://github.com/Joseph7e/Assign-Taxonomy-with-BLAST).

Note that the rationale behind the three first methods is to assign a taxon at each taxonomic level (e.g. species, genus, family) and then to evaluate the reliability of the prediction in a second stage to choose a reasonable taxonomic level. These three methods always provide a predicted taxonomy, with the rare exception where blast finds no match (representing <2% of the sequences for particular cases). With TopNPlus, we must choose a priori a series of criteria that will prevent the algorithm from providing any prediction at a given taxonomic level if these criteria are not met. Sometimes none of the top 10 sequences matches the requirements. In that case, no taxon is attributed.

Comparing four blast taxonomic assignment approaches

The methods using the most frequent taxon among the top 10 blast hits (TopN, TopNPlus) never outperformed the methods using the best blast hit (TopHit, TopHitPlus) and, in many situations, they performed worse (Figs. 1 and 2).

The method we called TopHitPlus outperformed the simpler TopHit approach in several circumstances, particularly for rbcL at species level (up to 13.6% more correct identifications at species level when restricted sequences from local species were blasted against a world general database). The maximum benefit for ITS2 reached 3.9% better prediction for TopHitPlus relative to TopHit. The TopHit approach never outperformed TopHitPlus or with a very limited improvement (maximum 0.3% of better assignments). The rbcL barcode is less specific and ties between the highest bit scores are more frequent than for ITS2. In this case, preferring the most frequent taxon (as performed by TopHitPlus) provides a slight advantage.

For the remainder of this article, we will concentrate on the results of this TopHitPlus approach.

Leaked CV vs. 10-fold CV

The accuracy estimates obtained with leaked cross validation (CV) were better than with the 10-fold CV approach (Fig. 3). However, even though in that case, the true sequence is present in the reference database, blast had still difficulties retrieving the correct species level identification among all the other sequences in 4–14% of the cases for ITS2 and in 16–53% of the cases for rbcL (depending on the strategy used). The lower accuracy for rbcL is most likely due to the lower specificity of this barcode which makes the distinction between closely related species more difficult.

As expected, 10-fold CV (left panels on Fig. 3) accuracy values were clearly worse than for leaked CV. In real biological samples, we can expect that the exact target sequence will sometimes be present in the reference database and sometimes not. So the accuracy in real biological samples is probably intermediate between the results obtained with 10-fold CV and leaked CV. In an unpublished dataset based on 354 real pollen samples, we observed a 100% identity match for only 40.3% of 12,971 ITS2 ASVs and 64% of 13,132 rbcL ASVs. So in practice the exact sequence does not seem to be always present in the reference database (or at least not matched). In the remainder of this article, we will therefore concentrate on the 10-fold CV results which represent a more realistic scenario.

Restricted vs. general databases

Using 10-fold CV, the accuracy obtained with the database restricted to the region amplified by a specific primer was systematically worse than with the general database containing the full-length, untrimmed sequences (see left panels on Fig. 3, yellow lines). This is probably because these restricted databases contain much fewer sequences after trimming, cleaning and dereplication. So, once the true sequence is removed from the reference database, blast struggles to find alternative sequences taxonomically close to the true target sequence. Other studies on different barcode genes have also reported a decrease in taxonomic accuracy when using reference databases with sequences trimmed to a specific primer pair, i.e. on COI gene for Arthropoda and Chordata (Robeson et al., 2021) or on fungal ITS gene (QIIME2 Documentation, 2022).

When restricted sequences were blasted against a general database, the accuracy was similar to the one obtained when general databases were cross-validated against themselves (Fig. 3, left panels, red and brown lines). We expected that the longer sequences in the general vs. general cross-validation would artificially increase the accuracy thanks to the additional DNA portions providing more information to distinguish closely related species. However, this is not what was observed from the data. When we use the general database, blast probably struggles to find alternative sequences (once the true sequence has been removed by the cross-validation process) for many very long sequences.

In the remainder of this article, we will therefore focus on the results for the general database CV as a more general case with more tested sequences.

Local and foreign species vs. local or world databases

Blasting the sequences of local species (i.e. species present in Belgium in our study case) against a local reference database (i.e. containing only sequences of Belgian species) instead of a world database (i.e. with all NCBI sequences) increased the species-level accuracies by $\sim$12% for both rbcL and ITS2 (Fig. 4). For ITS2, the accuracy rose from $\sim$63% to $\sim$75% when we used a world database or a local database respectively. For rbcL, the accuracy rose from $\sim$46% to $\sim$58%. Despite this improvement, the error rate at species-level is still rather high (25% for ITS2 and 42% for rbcL). Of course, when sequences of foreign species were blasted against a local database, the accuracy at the species level was 0%, while when using a world database, the accuracy increases weakly to 41% for ITS2 and 20% for rbcL. So, the identification of foreign species at species level was always poor even when a world database was used. This is probably due to the fact that this foreign flora is not as well represented in the databases as the local flora considered in this study case (i.e. Northern Europe flora). So once the true sequence was removed from the database, blast struggled to match alternative sequences that are close to the real sequence.

However, at higher taxonomic levels, using a world database is the best option in all cases. Indeed, for local species, using a local or world database made little difference at the genus level ( $\sim$94% accuracy for ITS2 and $\sim$ 90% for rbcL; Fig. 4). For foreign species, however, using a world database caused a substantial improvement ( $\sim$93% accuracy for ITS2 and $\sim$81% for rbcL). If a local database is used instead, these accuracies for foreign species dropped to $\sim$41% for ITS and $\sim$32% for rbcL. At the family level, the accuracy was $\sim$99% when a world database was used, for both foreign and local species.

Generalized linear models (GLMs)

As expected, the probability of correct identification increased when both the identity and consensus scores increased at all taxonomic levels (Fig. 5). Local and foreign species followed a similar trend with, on average, a higher probability of correct identification for the local species. However, there was no perfect separation line between correct and incorrect identification based on the combination of these two scores. Even with a 100% identity score, many of the identifications remained wrong at the species level. Even when the consensus score is as low as 10%, 20–40% of the identifications remained correct at species level. This means that as soon as we set a consensus or identity threshold under which we do not trust the identification, we will also discard a lot of “correct” identifications.

Classification trees

Classification trees can help define simple rules based on the identity and consensus scores. Instead of a GLM, we fit a classification tree, with the correct/incorrect identification as response and the consensus and identity scores as predictors at each taxonomic level, for each barcode and World vs. Local databases. An example is displayed on Fig. 6 for rbcL at species level (see Supplements section 4.3 for the trees corresponding to other cases). The nature of this recursive partitioning algorithm makes it particularly well suited to “detect” best thresholds (rather than linear increases). A visual examination of figures like Fig. 5 with univariate binomial GLMs is also helpful in making better-informed decisions on the best thresholds.

The rules “proposed” by the tree must always be adapted depending on the type of errors one favours. For example, here, we have established a set of rules quite stringent at the species level. The aim was to eliminate as many false positives as possible since we know that the percentage of correct identification is rather low at the species level. On the contrary, much less stringent rules were chosen at the genus level. The error percentage would already be minimal if we accepted all the identifications at the genus level (particularly for ITS2). So, the aim was to eliminate only the most obviously doubtful identifications.

Based on the results from the classification trees (Fig. 6) and the univariate GLM graphs (Fig. 5), the following rules were decided: at the genus level, we used only the World DB results and discarded the taxonomic assignments with Consensus < 40% and Identity < 98.5% for ITS2 or Identity < 99.7% for rbcL. At the species level, for ITS2, with a Local DB, we discarded all predictions with Identity < 99%. Then, for the sequences whose identification was discarded with the local database approach, we considered the results from the World DB and discarded all sequences with a consensus < 40% and an Identity < 99.9% or a consensus < 20% (even if the identity was 100%). For rbcL, with a Local DB, we discarded all predictions with Identity < 99.7% or with an Identity <= 100% and a consensus < 50%. Then, we also used the World DB and discarded all sequences with a consensus < 30% or a consensus < 60% and an identity < 100% (NB: such an approach combining Local and Global databases was for example used by Bell et al. (2020) or Casanelles-Abella et al. (2022)).

Reduction of prediction errors after applying the rules based on identity and consensus scores

Figure 7 shows what happened when these rules were applied in our study case. The percentage of sequences wrongly identified at the species level dropped from 37–80% to 5–17%. However, this comes at a cost: the sequences without identification rose from 0% to 17–92%. For foreign species with rbcL for example (worst case), the percentage of prediction error was 80% (on a World database) before filtering out low quality assignments. After applying these rules, the percentage of errors dropped to 5.5% but the percentage of unassigned sequences skyrocketed at 91.5%, with only 2.9% of the sequences correctly identified. Even for the best case (local species with ITS2), after filtering, $\sim$17% of sequences were unassigned, $\sim$68% were correctly identified, and $\sim$17% displayed wrong identifications. This means that, even after filtering, 20% of the assignments (17/(17 + 68)) were still incorrect at the species level. At the genus level, the gains were marginal (the raw predictions were already good enough). The most important change concerned the assignment of foreign species with rbcL: the error rate dropped from $\sim$19% without rules to 12% after applying the rules. In conclusion, in our study case, while the identification accuracy at the genus level was satisfactory, it seemed difficult to remove the wrong identifications at species level without also removing most of the correct ones, even by carefully choosing the consensus and identity thresholds.