Pollen metabarcoding reveals broad and species-specific resource use by urban bees

Field sampling

Floral cover assessment, pollinator sampling and pollinator observations were carried out across 10 pollinator planting strip sites managed by Bournemouth Borough Council in July 2016 (Tables S1 and S2; Fig. S1). All data was collected between 9 am and 5 pm on dry days where the wind speed was less than 5 on the Beaufort scale and the temperature was above 15 °C.

At each site three 1 m² quadrats were placed by selecting the patches with the highest density of open flowers. The total floral cover as a percentage of each plant species was measured in these quadrats (vegetative growth was not recorded). At each of the three quadrats, a single 10 min pollinator count was carried out by recording all pollinators to enter the quadrat. Pollinators were identified by observation only, without netting. Honey bees and bumblebees were identified to species level, while other bees were identified to family level where possible and recorded as ‘other solitary bees’ if not.

On the same day as observational data was collected, pollinators were sampled whilst they were visiting flowers within wildflower strips. As abundances were not being measured, sampling continued until 15 individuals had been caught, regardless of how long this took. Pollen was non-destructively collected from pollinators by confining pollinators in sterile microcentrifuge tubes: 1.5 ml (for A. mellifera and Bombus spp.) or 0.2 ml (for Halictidae spp., hoverflies and beetles). Pollinators were contained in the tubes for 5 min each, in a cool place, allowing pollen to be deposited on the tube walls by insect movement, and then released. Pollinators were captured whilst feeding, and the species of flower visited was also recorded. Bombus lucorum and B. terrestris could not be distinguished in the non-destructive field during sampling, so they have been grouped under the name B. terrestris in this study. Similarly, individuals from the Halicidae could not be reliably identified in the field so individuals are classified at the family level.

DNA extraction

DNA extraction was carried out following the method described by Hawkins et al. (2015). Pollen was resuspended in 400 μl of buffer AP1 from a DNeasy Plant Mini Kit (Qiagen, Venlo, Netherlands), to which 80 μl of proteinase K (one mg ml⁻¹; Thermo Fisher Scientific, Waltham, MA, USA) was added alongside one μl of RNase A (100 mg ml⁻¹; Qiagen, Venlo, Netherlands). Next, the samples were disrupted by shaking for four minutes at a speed of 30 1/s in a Retsch MM200 bead mill with custom adapter. Subsequent steps of the DNeasy Plant Mini Kit were carried out following manufacturer’s instructions, with the omission of the QIAshredder column. Following extraction, DNA was stored at −20 °C.

To test whether the non-destructive method of sampling provided the same information as destructive sampling, 10 honeybees (Apis mellifera) were collected from outside hives in Dorchester (latitude: 50.719°, longitude: −2.419°) and six bumblebees (Bombus terrestris) individuals were collected on Bournemouth University campus (latitude: 50.741°, longitude: −1.894°). Collection was carried out on the 11th of November 2016 for honeybees, and 13th of February 2017 for bumblebees. Each of these individuals was both destructively and non-destructively sampled in a paired test: first, pollen was non-destructively collected from each individual as detailed under ‘Field Sampling’, and then the whole insect was transferred to a fresh tube and frozen at −20 °C prior to further processing. Pollen was ‘washed’ from each insect following the method employed by Lucas et al. (2018a, 2018b): one ml of 1% sodium dodecyl sulphate (SDS) and 2% poly-vinyl pyrrolidinone (PVP) solution in water was added to the insect in each tube. Tubes were vigorously shaken by hand for 1 min, allowed to stand at room temperature for 5 min, and finally shaken vigorously by hand for a further 20 s. Next, the insect was removed and the tube containing pollen and the SDS-PVP solution was centrifuged at 13,000 rpm. Finally, the supernatant was removed and discarded, and DNA extraction was carried out as described above. Each extraction was tested by PCR amplification.

To prevent contamination, all DNA extractions were carried out in a laminar flow hood. Prior to each DNA extraction, surfaces within the hood were cleaned with 10% bleach followed by 95% ethanol, then all reagents and tools were placed within the hood and irradiated with UV light for at least 15 min. The hood was UV irradiated for 1 h every night. A negative control was included with each batch of extractions.

Library preparation and sequencing

A section of the rbcL gene was amplified and prepared for sequencing following the protocol of De Vere et al. (2017), adapted from that described by Illumina for the V4 region of 16S rRNA genes in bacteria (Illumina, 2013). The rbcL gene was chosen because a complete rbcL database has been created for native plants within Wales (De Vere et al., 2012), containing the majority of plants found in the UK as a whole. This protocol involves two PCR amplification steps: one to amplify the region of interest, and a second to add index and adapter sequences for sequencing. Following each PCR, samples are purified using AMPure beads. The final step involves library quantification, normalisation, and pooling. Only samples which produced a visible band in the first PCR step were carried through to further library preparation steps.

First, the rbcL gene was amplified using primers described in Table S3. Each reaction was at a final volume of 20 μl, and contained two μl of template DNA, Phusion High-Fidelity Master Mix at 1X concentration (New England Biolabs, Ipswich, MA, USA) and primers at 0.2 mM each. Thermal cycling conditions were as follows: 95 °C for 2 min; 35 cycles of 95 °C for 30 s, 50 °C for 1 min 30 s, 72 °C for 40 s (40 cycles); and 72 °C for 5 min. PCR clean-up was carried out using AMPure XP beads (Beckman Coulter, Brea, CA, USA) according to manufacturer’s instructions. The second stage PCR was carried out at a final volume of 25 μl, with each reaction containing 12.5 μl Phusion High-Fidelity Master Mix, 2.5 μl Nexera XT Index Primer 1 (N7XX), 2.5 μl Nextera XT Index Primer 2 (S5XX), five μl water and 2.5 μl of purified product from the first PCR. Thermal cycling conditions were as follows: 95 °C for 3 min; eight cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s; and 72 °C for 5 min. A second PCR clean-up was carried out as described above. The product was quantified using a Qubit Fluorometer and Qubit dsDNA HS Assay Kit according to manufacturer’s instructions, and pooled at equal concentrations to generate the final library pool. Prior to sequencing, the library was again quantified by Qubit and adjusted to 10 nM concentration with 0.1M Tris-HCl/0.01% Tween 20 solution, prior to denaturing and loading onto an Illumina MiSeq (Illumina, San Diego, CA, USA) following manufacturer’s instructions.

Two negative controls were included in the sequence run: one containing a randomly chosen DNA extraction negative control, and one containing purified water.

Data analysis

Sequencing data was analysed using a workflow previously described by De Vere et al. (2017) and available at https://github.com/colford/nbgw-plant-illumina-pipeline. Adapters and low-quality bases were trimmed using Trimmomatic (Bolger, Lohse & Usadel, 2014), then paired-end reads were merged using FLASH (Magoc & Salzberg, 2011). Singleton reads and merged sequences less than 450 bp in length were removed. Next, megablast (McGinnis & Madden, 2004) was used to search unique sequences against a custom BLAST database which consisted of all sequences from the Barcode Wales project (De Vere et al., 2012) alongside selected other sequences downloaded from GenBank (Benson et al., 2012). Results were manually filtered to remove plants that do not occur in the UK, based on Stace (2010), and Cubey & Merrick (2014).

All further analyses were carried out in R (R Computing Team, 2017). Rarefaction curves were generated using the R package ‘vegan’ (Oksanen et al., 2017), and rank-abundance curves in BiodiversityR. Bipartite pollinator–plant networks were drawn in R package ‘bipartite’ (Dormann, Gruber & Fruend, 2008). To test for significant differences in the number of reads and genus diversity between species and sites, generalized linear models were fitted with poisson (where no overdispersion was detected) or quasipoisson distributions using function ‘glm’. Post hoc Tukey comparisons were carried out using package ‘lsmeans’ (Lenth, 2016).

Comparison of destructive and non-destructive sampling

To compare destructive and non-destructive sampling methodologies, DNA was collected using both methods from ten A. mellifera and six B. terrestris individuals. A single A. mellifera individual and four B. terrestris individuals produced visible bands following PCR amplification, and were thus sequenced. For destructive sampling, two A. mellifera and six B. terrestris individuals produced visible bands.

Overall, five genera were detected across all individuals collected for this part of the study: Camellia, Erica, Hedera, Ulex, and Viburnum. Hedera was the sole genus detected on A. mellifera, while Camellia, Erica, Ulex, and Viburnum were detected on B. terrestris. A minimum of 98% of genera in all pairs of destructively and non-destructively collected samples were shared between both samples (Table 2), and overall community composition was broadly similar between methods (Fig. 1). In all but one case, non-destructive and destructive sampling detected the same ‘major’ plant genera (defined as genera making up more than 5% of reads): the only exception to this was bee ‘Bt3’, where Ulex was detected in both samples but made up less than 1% of reads in the sample collected using non-destructive methods.

Variation between bee species

The sequenced samples comprised 21 B. terrestris individuals, one individual each of B. hypnorum and B. lapidarius, two B. pascuorum individuals, nine A. mellifera individuals, and seven belonging to family Halictidae. Both hoverflies (Diptera: Syrphidae) and beetles (Coleoptera) were collected, but did not yield useable DNA. The numbers of pollinators collected and sequenced at each site are described in Table 1.

The three most abundant pollinator taxa in the sequencing dataset (A. mellifera, B. terrestris and Halictidae spp.; Table 1) were chosen for intra-species and intra-site comparisons. Reads from B. terrestris belonged to the broadest range of plant genera, followed by Halictidae and then A. mellifera (Table 3). Of the plant genera which individual bees were collected from, metabarcoding detected five of seven plant genera visited by B. terrestris, two of the four genera visited by A. mellifera, and four of the seven genera visited by Halictidae spp. Additionally, metabarcoding detected a number of plant genera which were not present in wildflower plots (Table 3).

There were significant differences in genus richness of pollen found on individual insects of different species (resid. dev. = 36.1, df = 27, p = 0.02) and at different sites (resid. dev. = 47.1, df = 34, p = 0.03). In particular, pollen from B. terrestris individuals contained significantly more plant genera than that from A. mellifera (z = 2.5, p = 0.04). The mean number of plant genera per individual (±S.E.) was 4.2 (±0.4) for B. terrestris, 2.7 (±0.5) for A. mellifera, and 4.1 (±1.5) for Halictidae. The minimum number of genera detected on an individual was one for all bee species, with a maximum of seven for B. terrestris, five for A. mellifera, and 12 for Halictidae.

Use of sown wildflower strips by pollinators

Metabarcoding of pollen collected from the bodies of bees detected pollen from a wide variety of plant families, many of which were not present in the wildflower plots (Fig. 2; Table 3).

In order to allow comparison between the observational and metabarcoding datasets, only data for bees which were sequenced is shown in Fig. 2A. Bees were collected from the flowers of nine different plant genera, all of which were flowering within wildflower strips at the time of sampling. The largest numbers of insects were collected on Phacelia, followed by Centaurea and Chrysanthemum (Fig. 2A). Conversely, while both Phacelia and Chrysanthemum were abundant in DNA metabarcoding data, Papaver made up a larger proportion of metabarcoding reads than of available floral resources (Fig. 2B) and Centaurea was not detected by DNA metabarcoding. In addition, a number of plant genera were detected on bee bodies, but were not detected in floral surveys: particularly abundant were Ligustrum, Rosa, and Achillea (Fig. 2B). These genera were often found on only a small subset of bees.

Plant genera present in wildflower plots accounted for the majority of reads (Fig. 3B), making up 69% of reads overall. Particularly well represented plant genera were Papaver and Phacelia. Conversely, the three members of the Asteraceae (Anthemis, Centaurea, and Chrysanthemum) found in sown wildflower plots were typically poorly represented amongst sequencing reads relative to the proportion of flowers available at each site: an exception to this was site 29, which was dominated by Chrysanthemum (Fig. 3) in the sequence reads. For three individuals, all pollen was assigned to a single genus that was not present in the wildflower strips: instead, all sequences were assigned to Achillea, Sonchus, or Rubus. Other insects appeared to mix genera found within the wildflower plots with genera from external sources.

Comparison of destructive and non-destructive sampling

Previous studies have employed destructive sampling methods to collect pollen from individual insects (Pornon et al., 2016; Lucas et al., 2018a, 2018b), in which bees are killed prior to pollen removal. While effective, this approach may be problematic under certain circumstances, for example, sampling of rare non-colony-forming pollinators, such as threatened species of hoverfly (Ball & Morris, 2014) or the sampling of bumblebee workers early in the season, which has lasting effects on colony reproductive success (Muller & Schmid-Hempel, 1992). The non-destructive sampling method employed in the current study involved trapping individual pollinators in Eppendorf tubes for long enough for pollen to fall off, then releasing the insect to continue foraging. However, using this approach the number of individuals yielding sufficient pollen for analysis was considerably lower, particularly for hoverflies and solitary bees, meaning that the utility of non-destructive sampling methods may be limited unless refinements can be made to improve insect contact with the collection tube.

Overall, where sufficient DNA was obtained for amplification with both methods, both destructive and non-destructive sampling of pollen gave broadly similar results. In all sample pairs, the vast majority of reads belonged to genera present in both samples. Additionally, in all but one case, the number of plant genera detected was the same in both destructively and non-destructively collected samples once rare genera were excluded. The proportion of reads assigned to each species was broadly similar between both methods, although some differences did exist: in particular, where multiple genera were present in a sample Ulex typically made up a greater proportion of reads when destructive sampling was used rather than non-destructive. It is possible that this is due to differences in the effectiveness of non-destructive pollen sampling in removing different pollen species, for example due to species-specific differences in pollen attachment (Pacini & Hesse, 2005). An additional factor to consider when comparing destructive and non-destructive sampling methodologies is that it is not possible to separate pollen stored in the scopa (pollen storage structure) from pollen found on the rest of the body. Pollen load analysis is commonly carried out on isolated scopal pollen in order to only include deliberately collected pollen rather than pollen accidentally collected when foraging for nectar (Wood, Holland & Goulson, 2016), or alternatively may be carried out with the scopa removed in order to only include pollen available for pollination (Pornon et al., 2016).

Variation between bee species

Each bee species carried a unique range of plant genera (Table 3), similar to previous studies which have identified differences in floral preferences between pollinator species (Geslin et al., 2013; Kells, Holland & Goulson, 2001; Leonhardt & Blüthgen, 2012). The largest number of plant genera were detected on B. terrestris and the smallest number on Halictidae spp., although it should be noted that this pattern reflects the number of individuals sampled for each species. There were also differences in the choice of flowers used by each bee species: in observational data, Halictidae spp. were primarily observed feeding on Asteraceae such as Chrysanthemum and Anthemis, while B. terrestris and A. mellifera were observed to use Phacelia. Metabarcoding data detected similar species-specific patterns of plant use: Chrysanthemum was detected on the highest proportion of Halictidae individuals, Papaver on B. terrestris, and Phacelia on B. terrestris and A. mellifera. This pattern reflects differences in tongue length: B. terrestris has a longer tongue than A. mellifera, which in turn has a longer tongue than Halictidae spp. (Prys-Jones, 1982; Williams, 1997). Halictidae individuals were found to feed primarily on open, brush and composite flowers, whose pollen and nectar is easier for short-tongued pollinators to access than tube-shaped flowers, for example, those of Phacelia (Inouye, 1980; Willmer, 2011). While there were very few observations and opportunities to sample mid- or long-tongued length bumblebees, the two B. pascuorum individuals (a longer-tongued species) which were caught with sufficient pollen loads for analysis were both found to be feeding predominantly on the tube-shaped flowers of the genera Phacelia and Echium.

The number of genera detected on individual bees by metabarcoding was unexpectedly high, particularly for A. mellifera. Earlier work suggests that A. mellifera foragers typically concentrate on a single plant species within a given foraging trip, while Bombus spp. and Halictidae feed on multiple species per foraging trip (Beil, Horn & Schwabe, 2008; reviewed by Grant, 1950; Leonhardt & Blüthgen, 2012; Wood, Holland & Goulson, 2016; but c.f. Brodschneider et al., 2018). Conversely, in the current study A. mellifera individuals carried pollen from an average of 2.7 plant genera. However, in many samples diversity was increased by plant genera which made up only a small proportion of reads: while metabarcoding data is only semi-quantitative (Elbrecht & Leese, 2015; Richardson et al., 2015a), the low abundance of reads from these genera raises the possibility that they do not represent genuine food plants, and are instead a result of pollen grains left over from previous foraging trips or transferred from the bodies of other bees in the hive or nest (Free & Williams, 1972). Alternatively, these plants may represent secondary, or ‘minor’, food sources, which may become more important as some pollinators exhibit more generalist foraging strategies in urban habitats (Geslin et al., 2013).

Use of sown wildflower strips by pollinators

All species of bee included in the study were found to utilise a range of plant genera which were not present in pollinator strips, demonstrating the ability of pollen identification to detect plant–pollinator interactions on a broader spatial scale then could be easily achieved with observation alone, especially given that honeybees, solitary bees and bumblebees can all forage over long distances (Beekman & Ratnieks, 2000; Knight et al., 2005; Zurbuchen et al., 2010). Overall, five of the 21 B. terrestris pollen loads were dominated by non-sown genera (>50% of reads): this is comparable to B. terrestris individuals foraging in sown wildflower strips on arable land (Carvell et al., 2006). Outside of wildflower strips, the majority of plant genera detected contained native members: for example, Rubus and Rosa. Previous work has found that native species are often preferred by pollinators (Corbet et al., 2001; Pardee & Philpott, 2014; Salisbury et al., 2015). Most of the plant genera originating from outside wildflower strips are commonly found in parks and gardens, either cultivated (e.g. Lupinus, Hydrangea, Buddleja, Ligustrum) or as wild plants (e.g. Rubus, Sonchus, Lactuca). Gardens provide a large proportion of urban green space (Loram et al., 2007; Thompson et al., 2003) and contain diverse plant species (Gaston et al., 2005). Alongside previous studies (Matteson & Langellotto, 2009; Osborne et al., 2008), the presence of garden species in the current study highlights the importance of gardens for urban bees. The abundance of species from outside the wildflower strips suggests that while wildflower strips alone are not enough to provision nearby bees adequately, bees are able to flexibly utilise a wide range of plant taxa and urban landscape features in order to obtain adequate floral resources.

Despite the fact that all the insects studied here were collected in urban pollinator strips, there was a relatively weak correspondence between the floral composition of the strips (i.e. which plants were flowering in strips at the time of sampling) and the composition of pollen collected from bees. The plant genera detected in the highest proportion of pollen samples were Chrysanthemum, Papaver, and Phacelia, suggesting that these plants may be valuable contributors to wildflower mixes sown to support urban bees. Each of these genera established well and produced large quantities of flowers at the time of sampling. Phacelia in particular is a common component of wildflower mixes sown in agricultural margins, and is often a significant component of foraging on these margins, particularly for B. terrestris (Carreck & Williams, 1997; Kells, Holland & Goulson, 2001; Pywell et al., 2005). Members of family Asteraceae are the plants most visited by small bees in agricultural margins (Wood, Holland & Goulson, 2017). In the current study, similar patterns were found: A. mellifera and B. terrestris disproportionately carried pollen from Phacelia and Papaver, while Halictidae spp. carried pollen from Chrysanthemum. Therefore, a mixture of these species appears to provide floral resources for the range of Hymenoptera studied. Phacelia and Papaver both contain high levels of protein and essential amino acids (Hanley et al., 2008; Roulston, Cane & Buchmann, 2000; Weiner et al., 2010), although this is not the case for Chrysanthemum (Roulston, Cane & Buchmann, 2000). However, it should be noted that all sampling was carried out in July, and it is likely that the relative contribution of different plant species to bee foraging varies over the course of the season.

In several cases, bees were observed to feed on plant genera within the wildflower strips that were not present in metabarcoding data: examples include Borago, Linaria, Leucanthenum, and Anthemis. Only DNA extractions which produced a band following PCR amplification were processed for sequencing, and this may have biased sequencing towards samples which were taken from pollen-foraging individuals since these individuals carry more pollen grains and thus are more likely to yield adequate DNA for amplification. In A. mellifera, individual bees are specialised for either pollen or nectar collection (Robinson & Page, 1989). Similarly, while Halictidae and B. terrestris females do not show individual specialisation, they may forage exclusively for either nectar or pollen on separate trips (Batra, 1964; Delph & Lively, 1992; Konzmann & Lunau, 2014). In the current study, pollen stored in the scopa was not separated from pollen on the rest of the body, and so it is not possible to distinguish deliberately collected pollen from pollen accidentally collected while foraging for nectar. However, previous studies show that pollen-feeding insects may make different foraging choices to nectar-foraging individuals: in particular, Asteraceae and Boraginaceae (families containing Anthemis, Leucathernum, and Borago, which were represented in observational but not metabarcoding data) are heavily used for nectar but not for pollen by bumblebees (Goulson et al., 2005), although it should be noted that small quantities of each family are found in bumblebee pollen loads (Kleijn & Raemakers, 2008).

However, detection may have been inhibited by methodology: while family Boraginaceae is well-detected by the primers chosen (De Vere et al., 2012), the amplification efficiency of different species in mixed samples is variable (Pornon et al., 2016). Additionally, rbcL may offer only poor discrimination for Asteraceae due to low levels of interspecific divergence (Gao et al., 2010), which were highly represented in wildflower plots. While at least one member of family Asteraceae (Chrysanthemum) was confidently identified to genus level in sequence datasets, a large number of reads could only be assigned to Asteraceae at the family level and may originate from species that were observed in floral surveys but not metabarcoding data (e.g. Anthemis).