Rank-based characterization of pollen assemblages collected by honey bees using a multi-locus metabarcoding approach.

PREMISE OF THE STUDY: Difficulties inherent in microscopic pollen identification have resulted in limited implementation for large-scale studies. Metabarcoding, a relatively novel approach, could make pollen analysis less onerous; however, improved understanding of the quantitative capacity of various plant metabarcode regions and primer sets is needed to ensure that such applications are accurate and precise. METHODS AND
RESULTS: We applied metabarcoding, targeting the ITS2, matK, and rbcL loci, to characterize six samples of pollen collected by honey bees, Apis mellifera. Additionally, samples were analyzed by light microscopy. We found significant rank-based associations between the relative abundance of pollen types within our samples as inferred by the two methods.
CONCLUSIONS: Our findings suggest metabarcoding data from plastid loci, as opposed to the ribosomal locus, are more reliable for quantitative characterization of pollen assemblages. Furthermore, multilocus metabarcoding of pollen may be more reliable than single-locus analyses, underscoring the need for discovering novel barcodes and barcode combinations optimized for molecular palynology.

Quantitative identification of pollen by taxonomic origin is important for applications in pollination biology and conservation (Kearns and Inouye, 1993; Wilson et al., 2010; Forcone et al., 2011; Girard et al., 2012; Cusser and Goodell, 2013), authentication of apicultural products (Louveaux et al., 1978; Jones and Bryant, 1992; Dimou and Thrasyvoulou, 2007), and allergy-related airborne pollen monitoring (Longhi et al., 2009; Kraaijeveld et al., 2015). Traditionally, pollen analysis has been accomplished using microscopic palynology, a technique involving the discrimination of pollen types by morphology (Erdtman, 1943). Due to the expertise required and difficulties associated with accurately distinguishing and identifying pollen from morphologically similar taxa, this technique has been difficult to implement on a large scale. Thus, the development and improvement of novel techniques for pollen analysis is an area of current interest (Keller et al., 2015; Kraaijeveld et al., 2015; Richardson et al., 2015).

The application of DNA barcoding to pollen analysis displays promise as an efficient and reliable approach. Similar to comparing morphological features of unknown pollen to those of reference pollen from voucher specimens, DNA sequences of unknown origin can be compared to sequences from voucher specimens. Initial applications of this approach to pollen analysis have used capillary sequencing technology (Longhi et al., 2009; Wilson et al., 2010; Galimberti et al., 2014); however, advances in the accuracy and length capabilities of next-generation sequencing provide researchers a practical, high-throughput alternative (Keller et al., 2015; Kraaijeveld et al., 2015; Richardson et al., 2015).

Recently, next-generation sequencing was used to characterize the botanical origins of bee-collected pollen using the ribosomal intergenic ITS2 locus (Richardson et al., 2015). This target locus was chosen because previous studies suggested that plastids are rarely incorporated into pollen (Reboud and Zeyl, 1994; Mogensen, 1996; Azhagiri and Maliga, 2007). However, evidence from more recent studies suggests that pollen plastids may be common (Tang et al., 2009), enabling pollen metabarcoding of plastid loci (Galimberti et al., 2014; Kraaijeveld et al., 2015). Although the approach using ITS2 was successful in identifying pollen (Richardson et al., 2015), it suffered from two limitations: (1) the method failed to detect certain prominent taxa identified microscopically and (2) while the method generated a useful taxonomic list, the relative abundance of different pollen types could not be inferred from the sequence data. Here, we present an improvement in pollen metabarcoding by targeting the plastid loci matK and rbcL, in addition to the ribosomal ITS2 locus, to characterize polyfloral samples of pollen collected by honey bees. In addition, we compare our metabarcoding results with results from microscopic analysis to evaluate the range of taxa detected and the capacity for quantitative inference of rank order pollen type abundance using a multilocus metabarcoding approach.

METHODS AND RESULTS

Sample collection and homogenization

During spring 2014, bee-collected pollen samples were collected at six apiaries, all greater than 15 km apart, in west-central Ohio. The latitude and longitude of each apiary is provided in Appendix 1, and apiaries are herein denoted as A, B, C, D, E, and F. Using Sundance I bottom-mounted pollen traps (Ross Rounds, Albany, New York, USA), we collected four samples from each site from 5–11 May, sampling every other day. After collection, samples were pooled by site before homogenization. A 10% subsample (by weight) was taken from each pooled sample, mixed in 50% ethanol, and stirred for 25 min using a magnetic stir plate. Using Buchner funnel vacuum filtration (Whatman grade 1; Sigma-Aldrich, St. Louis, Missouri, USA), we separated the homogenized pollen from the solvent and transferred it to a flow hood to air dry at room temperature.

Pollen identification and quantification by microscopy

We mixed 100 mg of the dried, homogenized pollen sample from each site in 0.5 mL of water and mounted five separate smears onto microscope slides in basic fuchsin jelly (Kearns and Inouye, 1993). We then counted and identified approximately 1000 pollen grains per slide for each pooled sample under a compound microscope at 400–1000× magnification. The voucher specimens used for pollen identification are listed in Richardson et al. (2015). A total of approximately 5000 grains were analyzed per sample. Due to the difficulty in distinguishing some related plant taxa (e.g., within Rosaceae [Moore et al., 1991]), we chose to limit microscopic identification to the family level. The total number of grains of pollen from each plant family, summed from each of the five slides, is available in Appendix S1.

Pollen identification by metabarcoding

After drying our homogenized samples, we freed DNA from 50 mg of pollen per sample using bead-beater pulverization (Mini-BeadBeater-1; BioSpec Products, Bartlesville, Oklahoma, USA) (Simel et al., 1997). Each sample was placed in a 2.0-mL microcentrifuge tube with 600 μL of lysis buffer from the QIAGEN DNeasy Plant Mini Kit (QIAGEN, Venlo, Limburg, Netherlands). Zirconium/silica beads (0.5 mm diameter) were added until the total contents of each tube reached 1.5 mL, and the sample was pulverized for 2 min. Then, 300 μL of deionized water was transferred to each tube and mixed with the contents and a 300-μL portion of the resulting lysate mix was transferred to a sterile 1.5-mL microcentrifuge tube. DNA was extracted using the QIAGEN DNeasy Plant Mini Kit (QIAGEN), and the ribosomal ITS2 and plastid matK and rbcL loci were amplified in separate PCR reactions. Amplification was conducted using previously published primer sets (Fay et al., 1997; Cuénoud et al., 2002; Chen et al., 2010) and the Phusion High-Fidelity PCR Kit (New England Biolabs, Ipswich, Massachusetts, USA) in a Mastercycler ep Gradient PCR machine (Eppendorf AG, Hamburg, Germany). Primer sequences, reagents, and PCR conditions for each barcoding locus are presented in Appendix 2. The ITS2, matK, and rbcL amplicons were subsequently purified using the PureLink PCR Purification kit (Life Technologies, Carlsbad, California, USA). At this point, 500 ng of purified PCR product for each locus was indexed independently using the NEBNext Ultra DNA Library Prep Kit for Illumina and NEBNext Multiplex Oligos for Illumina (New England Biolabs). Multiplexed samples were purified before being pooled (Agencourt AMPure XP; Beckman Coulter, Brea, California, USA). A final nine-cycle library amplification step was performed and samples were analyzed on a Qubit 2.0 fluorometer (Life Technologies) and an Agilent 2100 Bioanalyzer (DNA 1000 kit; Agilent Technologies, Santa Clara, California, USA) to ensure sample quality before sequencing. Paired-end sequencing was performed with the Illumina MiSeq platform using the TruSeq LT assay (600 cycles). Sequence data are available from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (accession code SRP055937).

Sequences were analyzed using an alignment-based approach. All computation was performed at the Ohio Supercomputer Center on a 12-core HP Intel Xeon X5650 machine with 48 GB of RAM. Reads were first trimmed by quality using Trimmomatic (v0.32; Bolger et al., 2014) with Phred scale 33 quality thresholds of 20 for both the 5′ and 3′ ends of each read. Reads less than 50 bp in length were discarded. Reads were then dereplicated to minimize PCR amplification bias and converted to FASTA format using the FASTX-Toolkit (version 0.0.13; http://hannonlab.cshl.edu/fastx_toolkit/). Next, reads were aligned against reference ITS2, matK, and rbcL plant sequences downloaded from NCBI GenBank on 23 September 2014. Reference libraries were constrained to only include plant species known to be present in Ohio and surrounding states based on the USDA Plants Database (http://plants.usda.gov/). Reference libraries are available in FASTA format in Appendices S2, S3, and S4. Venn diagrams showing the completeness of each of the reference libraries, at both the genus and species level, are presented in Appendix S5. Complete lists of the genera and species represented in each reference library are presented in Appendix S6. Alignment was performed using the BLASTN algorithm (version 2.2.29+; Altschul et al., 1997). Alignment quality-control thresholds were set as follows: E-value cutoff 1e-150, number of alignments 1, output format 0, number of descriptions 1. An additional setting, percent identity threshold, was used and its value differed between loci. For ITS2, we used a percent identity threshold of 95%, as in Richardson et al. (2015). However, given the relatively low sequence divergence between species at the matK and rbcL loci, we used a stringent setting of 99% identity. Following BLAST, we used MEGAN 5 (version 5.1.5; Huson et al., 2011) to taxonomically summarize our results with the following settings: min support 1, min score 50.0, max expected 1e-150, top percent 100.0, min complexity 0.00, min support percent 0.0 (off), paired end mode. Complete family-level and genus-level metabarcoding results are summarized in Appendix S7 and Appendix S8, respectively.

Analysis of results

After sequencing and implementing quality control, we obtained from 78,975 to 224,428 forward reads and 134,133 to 557,713 reverse reads across all 18 amplicon libraries. The median number of reads per locus was 258,987, 194,856, and 134,183 for ITS2, rbcL, and matK, respectively. In total, these reads had best hits to plant species from 49 families. To limit the potential for false identification, we limited our analysis using a consensus-based approach, counting only families found in more than one of the three amplicon libraries for each sample. Consensus lists of the families detected and their relative abundance in each sample are provided in Appendix S9. Using this approach, we confidently detected 25 plant families across the six sites. Using microscopy, 25 plant families were identified, six of which (Asparagaceae, Elaeagnaceae, Hamamelidaceae, Lamiaceae, Magnoliaceae, and Poaceae) were not identified by the metabarcoding consensus analysis. Although these families were detected microscopically, they were present at very low abundance, never constituting more than 0.5% of the 5000 counted grains in any sample.

To test the ability to infer the rank order abundance of different pollen types from the metabarcoding data, we conducted Spearman’s rank-based correlation between the number of mate-paired read alignments and the number of pollen grains per plant family for each locus individually as well as for the mean of the rbcL and matK loci, excluding the ITS2 data. We chose to exclude ITS2 because data from this locus exhibited poor quantitative capacity in a prior study (Richardson et al., 2015). Lastly, we calculated R coefficients for families detected across at least five of the six samples to determine which families were over- or under-represented in the metabarcoding analysis relative to microscopic analysis. The R coefficient is used in authenticating honey provenance (Bryant and Jones, 2001). In the context of this paper, the R coefficient is the quotient, for a particular taxon, of the relative abundance as inferred by metabarcoding and the relative abundance as inferred by microscopy. We conducted this analysis on rbcL data because this locus exhibited a broad scope of detection and was the only single locus to produce significant rank-based correlations when compared to the microscopy data.

Pollen from the families Rosaceae (commonly species of Malus Mill., Crataegus L., Amelanchier Medik., Prunus L., and other cultivated relatives) and Salicaceae (predominantly Salix L. spp.) comprised over 65% of our samples (Fig. 1). Pollen from plants in the Asteraceae (Taraxacum officinale F. H. Wigg.) and Oleaceae (Fraxinus L. spp.) were also abundant. Using Spearman’s rank-based correlation, we found moderate to strong associations between the rank order abundance of pollen types within our samples as inferred by the molecular and microscopic approaches. For the rbcL locus, ρ values ranged from 0.536 to 0.939, and the associations were significant for five out of six samples (Table 1). For the mean of rbcL and matK, the associations were significant across all samples and ρ values ranged from 0.570 to 0.939 (Fig. 2). When matK and ITS2 were analyzed separately, associations between the molecular and microscopic relative abundances were not significant for any sample (Table 1). In our analysis of average R coefficients, we found that certain families were consistently over- or under-represented in the molecular results relative to the microscopic results (Table 2). In particular, the average R coefficients for Brassicaceae, Caprifoliaceae, and Salicaceae were under-represented in the molecular data by greater than threefold relative to the microscopy data, while Fabaceae and Fagaceae were over-represented by greater than threefold (Table 2).

Fig. 1.

Proportional taxonomic abundances of each sample as estimated by microscopy. Unidentified pollen grains and taxa present at less than or equal to 1% of the sample are grouped as “Other.”

Table 1.

Spearman’s rank-based correlation between the total number of grains per plant taxon as determined by microscopy and the number of mate-paired aligned reads per plant taxon as determined by metabarcoding. Numbers indicate Spearman’s ρ, with the P value for each test in parentheses.

Sampling locationa	ITS2	matK	rbcL
A	0.381 (0.360)	0.469 (0.203)	0.587 (0.049)
B	−0.130 (0.658)	0.515 (0.072)	0.764 (0.001)
C	0.238 (0.582)	0.750 (0.066)	0.939 (<0.001)
D	0.204 (0.504)	0.262 (0.536)	0.575 (0.028)
E	0.067 (0.854)	0.483 (0.194)	0.536 (0.073)
F	−0.005 (0.989)	0.617 (0.086)	0.762 (0.002)

Letters indicate the apiary associated with each sample.

Fig. 2.

Rank-transformed taxonomic abundance as estimated by the mean number of rbcL and matK metabarcoding reads (y-axis) and the number of grains estimated microscopically (x-axis). Spearman’s ρ and P values are provided. Part letters indicate the apiary associated with each sample.

Table 2.

Average R coefficients for taxa present in at least five of the six samples.

Family	Average R	SD
Rosaceae	1.507	0.868
Asteraceae	2.922	2.744
Brassicaceae	0.190	0.116
Fagaceae	5.972	4.734
Oleaceae	2.674	5.186
Caprifoliaceae	0.091	0.017
Fabaceae	7.509	9.631
Salicaceae	0.097	0.083
Aceraceae	0.618	0.893

Note: SD = standard deviation.

CONCLUSIONS

We employed multilocus metabarcoding alongside traditional microscopic pollen identification, with the latter being considered the current standard of practice. Using a consensus-based approach, we found significant rank-based correlations between rbcL sequence abundance and microscopically examined pollen grain abundance for five of six samples. However, using the mean of rbcL and matK sequence abundance, we found significant associations between metabarcoding and microscopic results across all sites. This suggests that while the rbcL locus may be quantitatively useful, the simultaneous use of multiple loci may improve quantitative measurement of pollen abundance.

Our multilocus, consensus-based method exhibits promise as a powerful approach to pollen identification using metabarcoding. While no significant associations were found between matK sequence abundance and microscopy data, significant associations were found across all samples when matK sequence abundance was averaged with rbcL abundance. The poor performance of matK when used individually may be a result of incomplete universality displayed by the matK primer set (Chen et al., 2010). Despite its discriminatory power as a rapidly evolving plastidial coding region (Hilu and Liang, 1997), our data suggest the matK primer set used here may not be ideal for characterizing diverse pollen samples and may only be useful for supplementing data from other loci through average- or median-based analyses. Performing such analyses could enable researchers to both broaden the scope of detectable taxa and increase the quantitative capacity of metabarcoding efforts. Using one primer set to coamplify a genetic region across taxonomically diverse samples can be problematic, because priming site sequence divergence may hinder or prevent amplification for some taxa, potentially leading to under-representation or even nondetection in the metabarcoding sequence data. Employing a suite of primers enables researchers to overcome this limitation.

An additional metabarcoding issue involves minimizing the potential for false-positive identifications. Across a diverse sample, it can be expected that some closely related taxa exhibit little sequence divergence at a particular locus. Employing multiple loci in conjunction with consensus-based analysis limits the potential for false-positive identifications as the probability of the same false-positive identification occurring across multiple independent loci is decreased relative to the probability for a single locus. Lastly, the completeness of the reference database is crucial for the successful application of metabarcoding. Although none of the libraries used here were entirely complete with respect to Ohio taxa, a large majority of the known species were represented (Appendix S5).

Future research into different bioinformatic analyses, such as classifier-based analysis as opposed to alignment-based analysis, is warranted. The current alignment-based approach does not provide confidence estimates for individual sequence to taxon assignments. Classifier-based approaches are commonly used in microbial ecology, where they have been designed for the analysis of ribosomal amplicon libraries (Wang et al., 2007). Keller et al. (2015) successfully applied a classifier-based approach to ribosomal amplicons originating from pollen DNA, but to our knowledge, this approach has never been applied to nonribosomal loci, such as matK or rbcL. Successful application of this approach may enable researchers to better understand the confidence of taxonomic assignments on a read-by-read basis as well as across taxonomic ranks.

Although significant associations were found between the microscopic and molecular method, the presence of outliers cannot be overlooked (Fig. 2). Our analysis of family-specific R coefficients shows that some families were consistently over- or under-represented in the molecular results when compared to the microscopic results (Table 2), suggesting that, in addition to stochastic sampling error, some systemic mechanism may bias results. Such systemic biases could be attributable to aspects of pollen plastid biology, such as taxon-specific rates of plastid incorporation or relationships between average pollen grain volume and plastid abundance. To our knowledge, no studies have directly addressed such basic questions of plastid biology within pollen tissue. Alternatively, these biases may be the result of decreased amplification efficiency for certain plant families, resulting in nondetection or underestimation of abundance.

Unless validated, pollen metabarcoding data should be questioned in terms of its capacity for quantitative inference. Such validation requires comparison between the results of novel molecular approaches and the standard method of microscopic palynology. Contemporary studies have applied this approach, using statistical tests including Pearson’s product moment correlation, Spearman’s rank-based correlation, and generalized linear modeling to determine the quantitative capacity of metabarcoding techniques (Keller et al., 2015; Kraaijeveld et al., 2015; Richardson et al., 2015). One conclusion consistent with the analysis presented here, and that presented by Kraaijeveld et al. (2015), is that low-copy-number plastid loci provide generally quantitative results. However, studies disagree on the quantitative capacity of the repetitive ITS2 locus for metabarcoding (Keller et al., 2015; Richardson et al., 2015).

We employed multilocus metabarcoding to characterize the taxonomic composition of polyfloral honey bee–collected pollen. Requiring only minute quantities of pollen, our approach can easily be applied to studying the foraging habits of individual honey bees, as well as other ecologically and economically important pollinators, such as solitary bees and bumble bees. Our results suggest that sequencing plastid loci produces semiquantitative results. Furthermore, our results support the use of a multilocus, consensus-based approach over single-locus barcoding. Further research is needed to validate these trends across a larger sample of plant taxa and additional barcode loci.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Appendix 1.

GPS coordinates of each sampling location.

Sampling location	Latitude (°N)	Longitude (°W)
A	40.09	83.39
B	40.05	84.15
C	39.96	83.43
D	39.99	83.59
E	39.91	84.00
F	39.86	83.66

Table A1.

PCR mastermix reagents.

Reagent	Volume (μL)	Final concentration
Nuclease-free water	30.5
5× Phusion HF	10	1×
10 mM dNTPs	1	200 μM
10 μM Forward primer	2.5	0.5 μM
10 μM Reverse primer	2.5	0.5 μM
2 U/μL Phusion DNA polymerase	0.5	0.4 U/20 μL PCR
Template DNA	3	150 ng (total amount)
Total	50

Table A2.

Primer sequences and PCR conditions used for the amplification of plant ITS2, matK, and rbcL loci.

Loci	Primer sequences (5′–3′)
ITS2	F: ATGCGATACTTGGTGTGAAT
	R: GACGCTTCTCCAGACTACAAT
matK	F: CGATCTATTCATTCAATATTTC
	R: TCTAGCACACGAAAGTCGAAGT
rbcL	F: ATGTCACCACAAACAGAAAC
	R: TCGCATGTACCTGCAGTAGC

Table A3.

PCR conditions for ITS2 locus.

Step	Cycles	Temperature (°C)	Time
Initial denaturation	1	98	30 s
Denaturation	30	98	10 s
Annealing	30	59	30 s
Extension	30	72	30 s
Final extension	1	72	10 min
Hold	1	4	∞

Table A4.

PCR conditions for matK locus.

Step	Cycles	Temperature (°C)	Time
Initial denaturation	1	98	30 s
Denaturation	30	98	10 s
Annealing	30	58	30 s
Extension	30	72	30 s
Final extension	1	72	10 min
Hold	1	4	∞

Table A5.

PCR conditions for rbcL locus.

Step	Cycles	Temperature (°C)	Time
Initial denaturation	1	98	30 s
Denaturation	30	98	10 s
Annealing	30	57	30 s
Extension	30	72	30 s
Final extension	1	72	10 min
Hold	1	4	∞

End Prep Enzyme Mix (green)	3.0 μL
End Repair Reaction Buffer (10×) (green)	6.5 μL
Fragmented DNA	55.5 μL
Total volume	65 μL

Blunt/TA Ligase Master Mix (red)	15 μL
NEBNext Adaptor for Illumina* (red)	2.5 μL
Ligation Enhancer (red)	1 μL
Total volume	83.5 μL

*The NEBNext Adaptor is provided in the Multiplex Oligos for Illumina kit.

Adapter-ligated DNA fragments	23 μL
NEBNext High Fidelity 2× PCR Master Mix (blue)	25 μL
Index primer* (blue)	1 μL
Universal PCR primer* (blue)	1 μL
Total volume	50 μL

*The primers are provided in NEBNext Multiplex Oligos for Illumina.

Step	Cycles	Temperature (°C)	Time
Initial denaturation	1	98	30 s
Denaturation	9	98	10 s
Annealing	9	65	30 s
Extension	9	72	30 s
Final extension	1	72	5 min
Hold	1	4	∞