Enhancing Pollinator Habitat in Pacific Northwest Croplands Using DNA Metabarcoding Techniques

Study Sites

The field component of the project was conducted in 2018 and 2019 at three locations in eastern Oregon, US: Threemile Canyon Farms (45.7513° N, 119.9376° W), the United States Forest Service (USFS) Starkey Experimental Forest and Range (45.2332° N, 118.5511° W), and The Nature Conservancy’s (TNC) Zumwalt Prairie Preserve (45.5559° N, 116.9587° W). These locations were chosen with the goal of examining plant-bee associations in different land use types (i.e., agroecosystem, riparian forest, and grassland respectively). The three study sites occur within 350 km of each other.

Threemile Canyon Farms (Threemile) is a 37,636 ha industrial, pivot-irrigated farm located in Morrow County (elevation 100-300 m). Threemile has 9,308 ha in conservation area, 15,985 ha of irrigated conventionally managed cropland, and 6,178 ha of irrigated organically managed cropland. Uncultivated habitat includes the conservation area, which consists of arid grassland and shrub-steppe, as well as field margins, which are dominated by non-native vegetation. The Starkey Experimental Forest and Range (Starkey) is located in Union County (elevation 1,130-1,500 m); this long-term research site established in 1940 is grazed by cattle under a standard USFS grazing permit, and also supports herds of deer (Odocoileus spp.) and elk (Cervus canadensis). TNC’s Zumwalt Prairie Preserve (Zumwalt) is a 13,269 ha remnant bunchgrass prairie in Wallowa County (elevation 1,100-1,700 m). The Zumwalt Prairie has been used as summer pasture for horse, sheep, and cattle for over 100 years, but the majority of the area remains dominated by native plant species.

Bee Sampling

Bees were sampled during peak foraging hours (0900-1800) once a month from each site during June, July, and August of each year. Each bee was caught directly in a glass vial, if possible, or with an insect net and placed in an individual glass vial. All vials and nets were placed in 10% bleach for at least one minute and dried prior to sampling. Nets were replaced with clean, bleached nets after each use to prevent pollen contamination among samples. No killing agent was used in the vials, and bees were frozen at the end of each field day. Vials were labeled with the time, date, location, and flower species that the bee was foraging on when it was collected. The flower species noted for each bee was considered the bee foraging observation for that specimen. Each bee was given a unique sample ID so that it could later be associated with its specific pollen load. Bees were pinned, labeled, sexed and identified to the lowest taxonomic level possible, usually species.

Pollen Isolation

Bees were washed in the glass vial in which they were collected in the field. Sterilized water was added to the vial, and the vial was shaken vigorously until all visible pollen was removed from the bee. The pollen solution was pipetted from the vial and transferred to a 50 mL centrifuge tube. If pollen was still visible on the bee or the walls of the vial, it was rinsed with additional sterilized water, shaken vigorously, and the solution was pipetted from the vial and transferred to the same 50 mL centrifuge tube. The bees were then stored in a freezer until they could be pinned and identified. Each tube containing an individual bee’s pollen load was centrifuged at 2,000 rpm for 2 min, and the supernatant was discarded. The pollen pellet was resuspended in 1 mL of water, and the solution was transferred to a 1.5 mL screw cap microcentrifuge tube. The screw cap tubes were centrifuged at 12,500 rpm for 30 s and the supernatant was discarded. The pollen pellet was washed in 1 mL of 100% ethanol, and the supernatant was discarded. The pollen pellet was then dried for 30 min in an Eppendorf vacufuge. Dried pollen pellets were stored at -20°C until further processing.

DNA Extraction and PCR

DNA extractions took place using the Macherey-Nagel Nucleospin Food kit (Macherey-Nagel, Bethlehem, Pennsylvania, USA). We followed the “isolation of genomic DNA from honey or pollen” supplementary protocol. We added 1 mm glass beads to each 1.5 mL screw cap microcentrifuge tube and homogenized the samples in a Mini-BeadBeater-24 (Biospec Products Bartlesville, Oklahoma, USA). We included negative controls with each round of DNA extraction using sterilized water instead of pollen.

We used the dual-indexing pollen DNA metabarcoding strategy. We used the second internal transcribed spacer (ITS2) primers ITS S2F and ITS4R and existing universal ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit (rbcL) primers, rbcL2 and rbcLaR. Illumina overhang adapter sequences were added to each of the primers. The polymerase chain reaction (PCR) was conducted as a multiplex, targeting both sequence regions for amplification in one reaction. The PCR reactions contained 10 μL of 5X Green GoTaq® Reaction Buffer, 0.3 uL of 10 mM dNTPs, 1 uL of each primer, 0.2 uL of GoTaq® DNA Polymerase, 33.5 uL of water, and 2 uL of the template DNA for a total volume of 50 uL per reaction. PCR cycles included an initial period of heat activation for 3 min at 95°C. This was followed by 35 cycles of 30 s at 95°C, 30 s at 55°C, and 1 min at 72°C. There was a final extension of 10 min at 72°C and then the samples were held at 10°C until further processing. A negative control was included in each round of PCR, where 2 uL of water were used instead of template DNA. A 5 uL sample of each PCR product was electrophoresed in a 2% agarose gel, stained with GelRed^TM (Biotium Inc, Fremont, CA), and visualized under UV light to confirm the presence of the appropriately-sized amplicons. The PCR products were purified using the Promega Wizard SV Gel and PCR Clean-Up System. DNA was quantified for each sample using a Nanodrop 2000 Spectrophotometer (Thermo Scientific), and the samples were diluted with sterilized water in 96-well plates. The samples were indexed, pooled, and sequenced at the Center for Genome Research and Biocomputing at Oregon State University in a standard flow cell v3 300bp paired-end full service run of the Illumina MiSeq instrument.

Bioinformatics- Read Quality Filtering and Denoising

Raw reads were retrieved across all samples for ITS2 and rbcL, respectively. The open-source QIIME2 pipeline (version 2019-04) was used to analyze species composition of pollen loads. Initially, the QIIME2 DADA2 plugin was used to filter read quality and denoise reads. Sequence reads with a median sequencing quality score below 22 and noisy reads were filtered out. Contigs were created by joining paired-end reads, and duplicate sequences and sequences with chimera were removed. The resulting feature table containing optimized sequences and amplicon sequence variants (ASV) was used for additional analyses (e.g., taxonomic classification).

Bioinformatics- Reference Database Construction

Lists of all plant species known to occur at the three study sites have been developed over multiple years by botanists conducting field work at each site. A regional plant list was developed that included the unique species of each study site. Reference sequences of the ITS2 region and rbcL gene were extracted from publicly available data). A python package, NCBI-Companion (https://github.com/lixiaopi1985/NCBI_Companion), was used to query the NCBI nucleotide database to obtain sequences of species that were not present in published datasets. ITSx and MetaCurator software were used to detect and extract ITS2 and rbcL regions in the NCBI queried reference sequences, respectively. The final “regional” ITS2 and rbcL database was then split into three smaller “local” databases for ITS2 and rbcL that contained only species known to occur at each individual site.

Bioinformatics- Training Classifier and Taxonomy Classification

The Naïve Bayes algorithm provided by the QIIME2 feature-classifier plugin was used to train a classifier for taxonomic assignment based on the ITS2 and rbcL databases. Initially, the default settings were applied with high accuracy and low recall. In one case, taxonomic assignment of the classifier resulted in most of the samples being unassigned, so the kmer length was increased and confidence threshold reduced (high recall kmer = 32; confidence threshold = 0.6) to increase taxonomic assignments. The classify-sklearn tool of the QIIME2 feature-classifier plugin was used to assign taxonomy to each ASV.

Sequence Count Removal Threshold

Using the average and standard deviation of the number of sequencing reads in all negative control samples, we added 1.645 standard deviations to the average and used this number as our sequence count removal threshold. This allowed us to account for 95% of possible “background noise” detected in the pollen samples. We set the sequence count to zero for any taxonomic assignment whose read count fell below the sequence count removal threshold.

Network Analyses

Bipartite plant-pollinator networks were created, and network statistics were calculated using the bipartite package in R. Plant-pollinator networks were based on bee foraging observations and plant species assignments derived from DNA metabarcoding data using reference libraries of regional (MB-RDB) and local (MB-LDB) plant species. For each location, data from all sampling periods were combined into one network for each detection method (i.e., bee foraging observations, MB-RDB, and MB-LDB). Analysis of variance (ANOVA), blocked by location, was used to compare four parameters of pollinator networks: connectance (the realized proportion of possible links), H2’ (a frequency based index that increases with greater ecological specialization), linkage density (the mean number of links per species within the network), and bee generality (the mean number of plant species visited by each pollinator species within the network). Network parameter values were scaled to the highest value for a given parameters at a given site prior to analysis to account for differences in location. Tukey’s honest significant difference test (HSD) was used to compare means among detection methods. Analyses were conducted using R, version 4.0.0.