Bumblebees were collected in Cuxhaven (Germany: Lower Saxony) in proximity to the coastline (700 m) between the coordinates 53°51'36.3"N, 8°35'58.4"E and 53°52'55.4"N, 8°37'49.0"E, including the protected area Cuxhavener Küstenheiden (WDPA-ID: 329318, protected since 2004). The sampling sites varied between disturbed ground (suburbia), wet pasture, sandy scrubland and heath. Gardens, farmland (including pastures) and broadleaf forest were within 1 km radius. Sampling took place on the 2019-07-31 between 11 am and 4 pm in sunny weather. Permission to collect specimens and to enter the protected areas was granted by the municipal administration of Cuxhaven (Fachbereich 4: Naturschutzbehörde und Landwirtschaft). For sampling, we focused on female worker bees with visible pollen loads from as many different species as possible. We occasionally caught male bumblebees, which were determined after sampling, but were kept and analyzed alongside the female worker bees. To our best knowledge, no queen bees were caught. Bumblebees were directly caught in new, clean 50 ml centrifuge tubes (one specimen per tube) and then subsequently frozen overnight. We kept tubes upright, however, due to transport conditions, we cannot exclude minimal transfer of corbicula pollen to the bumblebee’s body. Pollen samples were taken within 24h thereafter by removing the pollen from the corbicula (if available) and by dabbing the bumblebee’s body with a toothpick covered in glycerol gelatin (except the hind leg). In 2019, we collected 120 specimens, of which 117 specimens were included in the following analysis.

Historic pollen samples from Cuxhaven were retrieved from the bumblebee collection of the Zoological Museum Hamburg (ZMH, Museum der Natur, Leibnitz Institute for the Analysis of Biodiversity Change, Hamburg). Female worker bumblebees were screened for pollen packages at their hind leg. We focused our efforts on samples collected by Rainer Wagner (Wagner 1969). The pollen packages were removed by precision tweezers (Dumont T5036). If the sample quantity was sufficient, one of the pollen packages was not removed for subsequent analysis. Pollen taken from the body hairs of historic bumblebees did not result in any PCR amplicons, regardless of the marker used (data not shown). In total, we sampled 81 historical bumblebees. These samples represent the subset of bumblebees caught by Wagner (1969), which were deposited in the museum and had visible pollen packages at their hind legs.

Bumblebees were identified independently by COI barcoding with DNA extracted from one hind leg, after pollen removal, and morphologically without disparities. Following the manufacturers’ protocol, the BioSprint 96 DNA Blood Kit (Qiagen) was used for automated DNA extraction with a Biosprint 96 (Thermo-Fisher). PCR targeting the 658 bp mitochondrial COI barcoding fragment was done using the primers HCO2198-JJ and LCO1490-JJ (Astrin and Stüben 2008). PCRs were set up in volumes of 20 μl (2 μl DNA template, 0.8 μl of each primer (10 pmol/μl), 10 µl PCR Multiplex Mastermix (Qiagen, Hilden, Germany) and 6.4 µl H₂O). PCR was performed by two different cycling regimes: 15 cycles (denaturation for 35 s at 94 °C, step-down annealing at for 90 s at 55 °C with –1 °C per cycle, extension for 90 s at 72 °C) immediately followed by 25 cycles (denaturation for 25 s at 94 °C, annealing for 90 s at 45 °C, extension for 90 s at 72 °C). Sanger sequencing with the same primers was outsourced to BGI (Hongkong, China). Sequences were assembled and analyzed by Geneious 7.1.8 (Kearse et al. 2012). Voucher specimens are accessible (deposited in Zoological Museum Hamburg: ZMH844169–ZMH844289).

To avoid contamination, we treated surfaces, labware and plastic equipment with UV light and additionally with DNA AWAY (ROTH X996.2, Karlsruhe, Germany) before use. DNA extractions and the setup of PCRs were performed in a sterile flow cabinet. The DNA extraction and PCR protocols were tested for various parameters and robustness (Kolter and Gemeinholzer (2021a), Suppl. material 1). Approximately 20% of all sequenced PCRs were blanks (PCR template negative controls). Due to the potential risk of cross-contamination into low-DNA samples (1968/69), the decision was made to exclude positive extraction controls from the experiment.

The DNA extraction buffer was modified from Sellers et al. (2018), while the extraction protocol was modified from the NucleoMag DNA extraction kit (REF744400.4 Macherey-Nagel, Düren, Germany). Samples were homogenized by bead milling in 1.5 ml tubes with six 1.5 mm stainless steel (grade SAE 316L) balls for 2.5 minutes at 30 Hz (MM400 Retsch, Haan, Germany) and DNA was subsequently isolated with a downscaled magnetic bead extraction protocol (Suppl. material 2). DNA was eluted by adding 35 µl of pre-warmed buffered H₂O (5 mM Tris/HCl pH 8.5).

The ITS2 PCR protocol and sequencing strategy was previously described in Kolter and Gemeinholzer (2021a). Identical protocols were applied to ITS1, using the primers ITS-2plR1 and ITS-u1 (Cheng et al. 2016), modified by adding a TruSeq HT (Illumina) primer tail to enable tagging by a subsequent PCR (Suppl. material 3). The trnL-P6 primers g and h by Taberlet et al. (2007) were modified in the same way. Due to their age, samples from 1968/69 were not amplified with the ITS1/2 markers, as initial trials revealed little success (data not shown). PCR was performed in individually tagged triplicates (Suppl. material 3). The paired-end 300 bp MiSeq sequencing run was performed on 1 ½ flow cells by LGC Genomics (Berlin, Germany).

Reference databases for ITS1/2 and trnL-P6 were generated from GenBank data files and filtered subsequently by a custom R script (Suppl. material 4). Reference database alignments were done with MAFFT and ITS regions were extracted using ITSx (Bengtsson-Palme et al. 2013; Katoh and Standley 2016).

The custom bioinformatic pipeline to analyze the MiSeq data used the R packages dada2 (featuring the DADA2 algorithm), vegan and ShortRead (Morgan et al. 2009; Callahan et al. 2016; R Core Team 2021; Oksanen et al. 2022). In contrast to ITS1/2, the read length of 300 bp was sufficient to always cover the whole trnL-P6 amplicon (< 160bp). To maximize the number of reads available for downstream analysis, the reverse (R2) sequencing reads were eliminated from the trnL-P6 workflow (Suppl. material 5). Amplicon sequence variant (ASV) generation by DADA2 was followed by multiple filtering steps. (Suppl. material 5). ASVs which could not be found in two out of three PCR replicates per sample were removed; this has been shown to minimize false positives and false negatives (Yang et al. 2021). The taxonomic assignment by SINTAX (Edgar 2016) was manually checked for plausibility. To address background noise-type contamination and more sporadic contamination types, a two-step process was implemented. First, ASVs whose sum of read count across all blanks surpasses 10% of their total read count across all samples were eliminated. Second, background contamination was mitigated by calculating the 90^th percentile of read counts across all blank samples for each ASV, respectively, and subtracting this value from the read counts of the corresponding ASV in all other samples. Any resulting negative values were adjusted to zero. ASVs only found in blanks were ignored (Suppl. material 5). Sequentially, read numbers of ASVs which were identified to the same taxonomic entity by SINTAX were conflated. The read numbers were transformed into presence / absence data. Pollination network graphs were generated by the R package bipartite (Dormann et al. 2009). The bioinformatic R pipeline, including a sample file, is available alongside the used reference databases (Suppl. material 6) and the ASV sequence data (Suppl. material 7).

For the comparison of the efficacy of the ITS1/2 makers, we assessed the number of detected taxa per sample after rarefaction, but excluded any subsequent filtering steps to minimize pipeline bias (Suppl. material 5). This was only applied to 2019 samples where both, ITS1/2, were successfully amplified. The Jaccard similarity was calculated for each of the taxonomic levels (family, genus, species) to compare the taxon lists recovered by the ITS1/2 and trnL-P6 loop. We ultimately decided to limit our analysis to presence / absence data because there currently is no consensus whether trnL-P6 read numbers show a significant correlation with pollen counts in the eDNA sample (Deagle et al. 2019; Baksay et al. 2020; Polling et al. 2021).

Of the 99 sampled bumblebees from 2019 (99 body swap samples + 18 corbicula pollen samples), 91 produced molecular data for all three plant barcoding markers. This dataset was used for foraging preference analysis and included the following bumblebee species, identified visually and via DNA barcoding: 40 B. terrestris, 24 B. lapidarius, 23 B. pascuorum, 3 B. pratorum and 1 B. cryptarum (Appendices 1, 2). Of the 81 corbicula pollen samples sampled from museum specimens (1968/69), we successfully generated NGS data of the trnL-P6 loop from 65 samples (39 B. pascuorum, 10 B. veteranus, 7 B. distinguendus, 5 B. hortorum, 3 B. lucorum, 1 B. muscorum).

To assess the technical performance of ITS1/2 and the trnL-P6 loop, we minimized the filter steps to increase sensitivity (Tables 1, 2, Appendix 4). As a result, we detected 37 plant families, 96 genera, and 61 species in the 99 bumblebee pollen samples from 2019 (Tables 1, 2, Appendix 4). The barcoding markers identified distinct sets of plant taxa. Differences in taxonomic resolution between the ITS1/2 markers were mainly observed at the genus and species levels (Appendix 4). Although the Jaccard similarity at the genus level was 0.84 (Table 2), the ITS1/2 markers can be considered roughly equivalent, given that most variation at the genus level involved plant taxa found in only one specimen respectively, such as Papaver or Clematis (Appendix 4). However, it is worth mentioning that the ITS1 marker detected 15 plant genera not detected by ITS2, while ITS2 only identified three genera not detected by ITS1 (Appendix 5).

In addition to variation in the identified taxa, the ITS1/2 markers also exhibited discrepancies in the sum of presence detections of a taxon across all samples (Table 2). Aggregating the presence detections for all plant taxa revealed that the ITS1 marker outperformed the ITS2 marker by approximately 20% in terms of identifications at the genus and species levels (714 vs. 562) (Table 1). This trend persisted even when excluding the taxa exclusively detected by ITS1 from the count.

The trnL-P6 marker detected the most plant families (n=35), compared to ITS1 (n=32) and ITS2 (n=31), but fewer plant taxa at genus level (trnL-P6 n=42, ITS1 n=77, ITS2 n=65) and species levels (trnL-P6 n=13, ITS1 n=52, ITS2 n=45) (Table 1). Considering the detection counts, the trnL-P6 loop detected more individual family signals (n=637) as the ITS1/2 markers (454 and 436, respectively) across all samples (Table 1). However, ITS1/2 outperformed the trnL-P6 on the individual genus and species level detections (Table 1). The families Asteraceae, Malvaceae and Rosaceae were most affected by the low trnL-P6 resolution (Appendices 4, 5). The low Jaccard similarity of detections between the trnL-P6 loop and ITS1 or ITS2 on genus and species level can be primarily attributed to the low resolution of the trnL-P6 loop (Table 1, 2). In Asteraceae, for example, the resolution of the trnL-P6 loop is restricted to family level (except Achillea), while the ITS1/2 markers successfully distinguished between 23 genera (Appendix 4).

We also calculated the similarity between the ITS1 and ITS2 data of the same sample individually (1 vs. 1) instead of comparing the whole sample pool (all ITS1 vs. all ITS2). Excluding genera exclusively being detected in either ITS1 or ITS2, the Jaccard similarity of bumblebee samples from 2019 (male and female, body and corbicula) between ITS1 and ITS2, per specimen, is 0.59 (n=110) (Suppl. material 7). The similarity between the markers rises asymptotically to 0.81 at a cutoff rate of 400 reads and reaches saturation at a similarity index of 0.88 at a cutoff rate of 870 reads (Suppl. material 7).

Data used to analyze the plant-pollinator interactions was filtered more strictly to exclude rare and infrequent plant signals that may be of minor importance in terms of their impact on the bumblebee colonies nutritional status (Table 3; Fig. 1, Appendices 1–3). The plant taxa identification in ITS1 (n=42) and ITS2 (n=45) were at least on the genus level (Appendices 1, 2). In trnL-P6, the highest taxonomic rank of the identified plant taxa (n=47) was either on family (n=12) or genus level (n=35) (Appendix 3). In detail, successful detections at genus or family level in the trnL-P6 loop, which were not detected by ITS1 and ITS2, can be partly attributed to Atriplex, Alnus, Robinia, Rumex, Acer, Pinus or Salicaceae (Appendix 4). In body samples (n=2) of B. pratorum the additionally detected taxa by the trnL-P6 loop in comparison to ITS1 and ITS2 were Convolvulaceae, Cucurbitaceae, Rosaceae, Ericaceae: Calluna, Fabaceae: Lotus, Lamiaceae: Mentha and Onagraceae: Oenothera (Appendix 3).

A quantitative count of the pollination network links (=connection) reveals that the markers show a different trend (Table 3). The median of detected taxa (statistical median of number of network links across specimens of the same species) in the trnL-P6 marker is always higher or equal to the other markers (Table 3). In contrast to most of the ITS1 or ITS2 evaluations (except B. pascuorum in ITS2), trnL-P6 recovered more taxa in the corbicula samples than the body samples (Table 3). Despite the high abundance of some taxa, e.g. Scorzoneroides, Tanacetum, Lotus and Lythrum, detected in more than 10 specimens in at least one bumblebee species by ITS1/2 (Appendices 1–3), the median number of plant taxa on individual female worker bumblebees was 2.5 to 4 for ITS1 and ITS2 (where n > 3 bumblebees, Table 3). The old corbicula samples (1968/69) did not recover less taxa when compared to the 2019 samples (where n > 3 bumblebees, Table 3). Plant taxa were found in a median of two to three bumblebee body samples (where n > 3 bumblebees, Table 3).

Figure 1.

Plant-pollinator inference network of historic and recent bumblebee specimen by trnL-P6 pollen metabarcoding. Corbicula pollen was sampled from bumblebees caught in 1968/69 (blue) and bumblebees caught in 2019 (orange). In addition, body pollen was sampled from all bumblebees caught in 2019 (yellow). The width of the colored bars reflects the sum of unique interactions (i) with plant taxa for all specimen (n) of the respective sample type (color). The width of the black bars reflects the total number of bumblebee specimens (s) in which the respective plant taxa has been found. Plant taxa were ordered to minimize network overlap. To avoid clutter, plant taxa found in less than two samples and bumblebee species represented by less than two specimens were omitted. Taxa represented by family names showed insufficient resolution in the trnL-P6 maker (e.g., Asteraceae).

85% (ITS1) and 75% (ITS2) of the plant taxa, at genus level, found in the corbicula samples (n=18) from female bumblebees caught in 2019 (B. terrestris, B. lapidarius and B. pascuorum) can also be detected in the respective body samples of the same specimen (Suppl. material 7). On family level the overlap rises to 91% for ITS1 and 78% for ITS2 (Suppl. material 7). Due to high similarity indices between corbicula and body samples of the same individual in ITS1/2, the corbicula samples can be broadly viewed as a subset of the body samples. For the trnL-P6 loop, where a higher number of plant taxa has been detected in the corbicula samples, 75% of the plant taxa detected in body samples have been found in the corbicula samples of the same specimen (Suppl. material 7). At family level, the overlap increases to 79% (Suppl. material 7).

In general, the most often visited taxa in 2019 were also detected in the 1968/69 samples and vice versa, albeit at a different frequency (Fig. 1). In detail, each of the well-covered plant families (s > 20) was found on each of the well-represented bumblebee species (n > 20) included in this study (Fig. 1). However, trends of plant-pollinator interactions in the 2019 samples versus the samples from 1968/69 using the trnL-P6 loop plant barcode marker are different. In detail, the ratio of recent / historic interactions per plant family (s > 20) varies with Asteraceae, Ericaceae and Lythraceae showing, relatively, more recent interactions and Fabaceae and Oleaceae showing the smallest ratio of recent interactions (Fig. 1). Phaseolus and Lathyrus (both Fabaceae), with one exception, can only be found on 1968/69 bumblebees. The difference in visitation pattern for B. pascuorum is especially striking, as we could rarely find Vicia and Trifolium plant-pollinator interactions in recent data, compared to historic B. pascuorum specimens (Fig. 1). Bombus pascuorum visited the same plant taxa (with at least 3 interactions in historic and recent data) in historical and recent collections, except for Phaseolus, which was additionally detected in the 2019 survey (Appendix 3).

To the best of our knowledge, this is the first pollen metabarcoding study of historic bumblebee pollen samples dating back ~50 years and the only bumblebee metabarcoding study reporting results from multiple endangered Bombus species. Existing pollen metabarcoding studies are primarily focused on B. terrestris (Wilkinson et al. 2017; Biella et al. 2019; Potter et al. 2019; Bänsch et al. 2020; Piko et al. 2021; Bontšutšnaja et al. 2021). This can possibly be explained by the difficulties associated in locating rare bumblebee species in the field. In the study of Beyer et al. (2020), only 1.2% of the bumblebees caught belonged to an endangered bumblebee species (Westrich et al. 2011). This underlines the usefulness of historic natural history collections in reconstructing trophic interactions between species, such as pollinator-plant interactions, especially when rare species are included (Scheper et al. 2014). However, it is important to recognize that plant-pollinator interactions must be supplemented with knowledge about the effects of shifts in flowering time, climate change, parasites or diseases, pesticides, and competition before conservation measures can be taken. (Goulson et al. 2015; Miller-Struttmann et al. 2015; Marshall et al. 2018; Soroye et al. 2020).

Our data on shifts in flower visitations revealed trends that differ between current and historical bumblebee specimens. The bumblebee species caught in 2019 are commonly reported to be present in urban environments and display a highly generalist foraging behavior (Banaszak-Cibicka and Żmihorski 2012; Zajdel et al. 2019; Sikora et al. 2020). This is different from two species only found in the historic samples (B. distinguendus and B. muscorum), which have been reported to be amongst the bumblebees with the narrowest dietary breadth (Wood et al. 2021). On the other hand, B. lapidarius, B. terrestris and B. pascuorum, in combination, have been reported in other areas with high population density, indicating their status as hemerophiles (Vray et al. 2019). Although we could show their generalist foraging strategy, most plant taxa could only been found on a small subset of all analyzed samples (Table 3). This indicates that a higher sampling depth is required to fully characterize the foraging habits. We can hypothesize that the greater frequency of Calluna flower visitation is due to the establishment of a strictly protected heathland that was not yet a formal protected area in 1969 (Wagner 1969). The higher visitation frequency of Lythrum could be a result of an increased amount of drainage channels, which provide an ideal habitat, due to intensified agriculture (Wagner 1969; Dierssen 1997). While difficult to predict, a shift in flowering time could also have played a role in floral availability (Rafferty et al. 2020), which, however, was not tested here.

Our analysis of recent and historic plant-pollinator interactions revealed distinct floral visits (Fig. 1). While the overlap of analyzed species is restricted to B. pascuorum, it can be inferred that Fabaceae are an important part of bumblebee diet (Wood et al. 2021), particularly for B. lapidarius (Hülsmann et al. 2015). Therefore, if available, all accessible floral resources from the Fabaceae family would have been used by B. pascuorum and B. lapidarius, as observed in Lotus (Fig. 1). In contrast to our results, bumblebee diversity was more extensive in the Cuxhaven area in 1968/69 (Wagner 1969). An important consideration for our results is whether the species is truly absent from the area, or we just could not find it. Therefore, we will discuss probable causes for its absence in our 2019 sampling effort.

B. ombus lucorum is reported to be a generalist forager (Wood et al. 2021), which aligns with our data (Fig. 1, Appendix 3). Reports of B. lucorum in the area are inconsistent. It has been classified as an abundant species in the area (Witt 2016), however, other studies report relatively low catch rates of 1.6% (Persson et al. 2015). In 1968/69 it accounted for 4.5% of the bumblebees captured (Wagner 1969). Overall, we conclude that B. lucorum was likely missed in 2019 due to insufficient sampling. Another possible explanation for the absence of B. lucorum in our samples is that the number of hot days (>25 °C) was higher than in 1968/69 (DWD 2022) and shortened the colony’s lifecycle (Maebe et al. 2021).

Our data confirms that B. distinguendus has a preference for Fabaceae (Diekötter et al. 2006; Witt 2016). However, we also observed the use of alternative pollen sources (Vray et al. 2017; Phelan et al. 2021), such as Rosa and Calluna (Appendix 3). In contrast to reports between 1909 and 1969 (Alfken 1913; Wagner 1969; Rasmont et al. 2015), B. distinguendus is assumed to be absent from most parts of Lower Saxony nowadays (Sprichardt 2010; Rasmont et al. 2015; Witt 2016), which may explain the absence in our study (Fig. 1). A similar trend of habitat destruction or extinction was also found in multiple other countries (Goulson et al. 2005; Charman et al. 2010; Dupont et al. 2011; Bommarco et al. 2012; Rasmont et al. 2015; Rollin et al. 2020). Common recommendations to protect B. distinguendus include modifying mowing practices of open grassland areas and converting pastures into species-rich grassland, provided suitable nesting sites are available (Charman et al. 2010; Witt 2016; Phelan et al. 2021).

The decline of B. hortorum, which is generally not considered a rare species in Germany, between 1959–1962 and 1968–1969 in the Cuxhaven area, can be attributed to land use changes (Wagner 1969). Although it was present in 2007 and 2009 in the Cuxhaven area, it was only observed in 2 out of 85 sampling days (Sprichardt 2010). We were unable to capture any B. hortorum individuals in 2019. Bombus hortorum is known to show an extreme dietary preference for Trifolium (Goulson et al. 2008; Kleijn and Raemakers 2008), and likely also other Fabaceae species. Interactions between bumblebees and Trifolium were underrepresented in our 2019 data, compared to the data from 1968–1969 (Fig. 1). This possibly hints at the fact that the land use changes were accompanied by a change in floral availability. While some studies have not demonstrated a dietary focus on Fabaceae (Wood et al. 2021), it is challenging to compare results across dietary studies as many parameters are not being controlled for (i.e., floral availability).

It has been hypothesized that B. veteranus is closely associated with Fabaceae (Goulson et al. 2008; Wood et al. 2021), and our data confirms this finding (Appendix 3). However, another study showed B. veteranus to be tightly associated with Cardueae (Vray et al. 2017), which was not reflected in our data (Fig. 1). Our data also revealed interactions with plants associated with anthropogenic influence, such as Ligustrum (Appendix 3). This suggests that the habitat of B. veteranus may be threatened by nearby land use changes, especially urbanization. Bombus veteranus is considered a rare species in Germany (Westrich et al. 2011) and has experienced rapid decline in Belgium (Rollin et al. 2020). In the Cuxhaven area, its last record is from 1969 (Wagner 1969). In Lower Saxony only one large population is known (Witt 2016).

In summary, the detected floral interactions of bumblebees caught in 2019 have shifted away from many Fabaceae genera. This is important for three reasons: 1) Host plant availability has been identified as the main driving factor of wild bee decline (n = 57 species), with Fabaceae dependent species (29 out of 57) showing the highest decline (Scheper et al. 2014). 2) Herbaceous Fabaceae, such as Vicia and Trifolium, are of principal importance for many rare bumblebees (Bäckman and Tiainen 2002; Goulson et al. 2005; Goulson et al. 2008; Kleijn and Raemakers 2008; Dupont et al. 2011; Timberlake et al. 2019; Wood et al. 2019; Sikora et al. 2020; Wood et al. 2021). And finally, 3) Fabaceae has been found to be the most effective plant family in mitigating negative effects of urbanization, promoting their usefulness also in urban landscapes (Hülsmann et al. 2015).

We tested three genetic markers, including the trnL-P6 loop, positioned between the trnL (UAA) exon 1 and the trnL (UAA) exon 2 (Taberlet et al. 2007), which has been utilized in multiple pollen and honey metabarcoding studies (Chiara et al. 2021; Milla et al. 2021; Polling et al. 2021). Compared to studies using pollen microscopy, it shows worse taxonomic resolution in the Asteraceae family and approximately the same taxonomic resolution in other plant families (Wood et al. 2019; Wood et al. 2021).

Comparing the trnL-P6 loop with ITS1/2 revealed two important findings. First, in accordance with Milla et al. (2021), we observed that the trnL-P6 loop is able to capture greater diversity at the family level than the ITS1/2 marker (Table 1). This increased diversity can be attributed to the trnL-P6 loop’s ability to amplify degraded DNA due to its shorter length (Valentini et al. 2009). There are multiple possible explanations for the additionally detected taxa in trnL-P6 compared to ITS1/2 in samples of 2019. Plant material attached could be accidentally deposited in the corbicula during combing alongside with pollen. Regurgitated nectar, used to fixate pollen on the corbicula (Thorp 2000), which also contains ingested pollen (Owen et al. 2013), could also add traces of degraded DNA to the pollen package. This could also explain why, in contrast to the ITS1/2 marker, the trnL-P6 loop was able to identify more plant taxa in the corbicula samples, compared to the body samples (Table 3).

Second, supported by Polling et al. (2021), in contrast to Milla et al. (2021), our results show that the ITS1/2 markers recover a higher number of taxa at species and genus rank (Table 1). After the manual curation of SINTAX results, the trnL-P6 loop identified only ~60% of the genera found by ITS1/2 in recent samples (Table 1).

These results demonstrate that sequencing shorter DNA fragments, such as the trnL-P6 loop, alongside with longer DNA fragments, such as the ITS1/2 barcode marker, will yield different insights and are well worth exploring. Unfortunately, we could not find comparable studies in literature and further controlled experiments are required to understand the differential detections of ITS1/2 and trnL-P6. In summary, the trnL-P6 loop was generally able to recover a higher taxonomic breadth, while the ITS1/2 maker was generally able to recover a higher taxonomic depth.

Our study shows that the ITS1 possesses favorable attributes, compared to the ITS2 marker. This can be demonstrated by more detected taxa on species and genus level, as well as the overall higher number of taxa in all samples (Table 1). The higher taxonomic resolution aligns with previous results (Wang et al. 2015; Kolter and Gemeinholzer 2021b). The increased richness of taxa per sample suggests a more even amplification profile of mixed samples, potentially explained by the more conserved flanking regions of ITS1 (Wang et al. 2015; Kolter and Gemeinholzer 2021a). Another possible explanation could be the, on average, lower GC content of ITS1 compared to ITS2 (Wang et al. 2015), resulting in less stable secondary template structures during PCR. Our study, due to improvements in methodology (Cheng et al. 2016; Kolter and Gemeinholzer 2021a), contrasts the findings of Chen et al. (2010), who excluded ITS1 as a barcode candidate due to amplification problems. The overall performance of ITS1 also contradicts the study of Gous et al. (2019), which, however, used primers which were also designed for fungal amplification (White et al. 1990; Kolter and Gemeinholzer 2021a). Subsequently, their findings are possibly a result of preferential amplification of fungal DNA and must be verified with plant-specific primers.

One disadvantage of ITS1 is the presence of extremely long ITS1 sequences in certain Gymnosperms (Cheng et al. 2016), which currently exceeds the technical limits of the Illumina sequencing platform (2×300 bp). ITS1 has been used sparingly in pollen metabarcoding studies (Pornon et al. 2016; Gous et al. 2019; Baksay et al. 2020; Gous et al. 2021), and further investigations, ideally including mock communities, are required, before it can be recommended over ITS2. In this context, it is important to mention that the number and quality of taxa recovered from any eDNA sample by metabarcoding depend heavily on the metabarcoding pipeline used (Pauvert et al. 2019). Finally, we conclude that the application of ITS1 in pollen metabarcoding studies needs more comparative studies but shows promise.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This DFG project (project number: 352447832) was part of the SPP 1991: Taxon-Omics: New approaches for discovering and naming biodiversity (project number: 313688472).

Andreas Kolter: Writing - original draft; Writing - review and editing; Data curation; Formal analysis; Investigation; Methodology; Validation; Visualization. Martin Husemann: ; Writing - original draft; Writing - review and editing; Investigation; Project administration;

Resources. Lars Podsiadlowski: Writing - review and editing; Investigation; Resources.Birgit Gemeinholzer: Conceptualization; Writing - review and editing; Writing - original draft; Funding acquisition; Project administration; Resources; Supervision.

Andreas Kolter https://orcid.org/0000-0001-9327-2053

Martin Husemann https://orcid.org/0000-0001-5536-6681

Birgit Gemeinholzer https://orcid.org/0000-0002-9145-9284

All of the data that support the findings of this study are available in the main text or Supplementary Information. Raw sequence reads are deposited in NCBI BioProject PRJNA841517.