Role of Bats in Controlling Agricultural Pests

Study Species and Area

We examined the diet of Big Brown (Eptesicus fuscus) and Red (Lasiurus borealis) in a typical agroecosystem of the United States Corn Belt. These species are two of the most abundant species in much of the Midwest United States and therefore the most likely to provided substantial ecological services to agriculture.

These species also believed to represent diverging dietary preferences. Previous morphological and DNA analysis of Big Brown bats show a fondness for beetles (Order Coleoptera) (Clare et al., 2014b). Conversely, Red bats are believed to primarily consume Lepidopteran species (moths) (Clare et al., 2009).

Fieldwork was conducted at two sites in Lancaster County Nebraska (Fig. 1). The Lancaster County (219,210 ha) landscape is dominated by row crop agriculture (56% of land area). Grassland, which is primarily used for haying, is the second most common land cover type (22%). Developed area in Lancaster county is concentrated around the city of Lincoln, NE and accounts for 13% of landcover (U.S. Geological Survey Gap Analysis Program, 2016).

Agricultural in Lancaster county is dominated by corn and soybean production. In 2016, 56,858 ha corn (25.9% of land area) and 54,511 ha soybeans (24.9% of land area) were planted in the county. Additionally, 19,667 ha (8.9% of land area) of alfalfa was harvested. Corn genetically modified to resist pests (bt corn) comprises 80% corn planted in NE (USDA National Agriculture Statistics Service). Bt corn in Nebraska is primarily used to provide resistance against the European corn borer and Corn root worm, however additional bt corn strains target corn earworm, Fall armyworm. Common stalk borer, black cutworm, and true armyworm are also available (Peterson et al., 2018).

Insect Pest Species

A list of 62 insect pests (Table s1) that damage corn, soybean, and alfalfa in the area was compiled using information from the University of Nebraska Extension office (Wright et al., 2012b, 2012a, 2009a, 2009b, 2009c, 2009d).

Guano collection

Bats were captured June-September 2016 using mist nets systems 2.6-10.4 m high. All captured bats were identified, weighed, and measured. Bat capture and handling was conducted in accordance to guidelines of American Society of Mammologist (Sikes and Gannon, 2011) and in accordance with University of Nebraska IACUC protocol (Project 1556). Bats were kept in cloth bags for up to 1 hour to deposit guano samples and released on site. Fecal pellets were stored in silica and frozen at the end of each sampling night.

DNA extraction and PCR amplification

Equal amounts of all fecal samples were extracted using a Qiagen QIAamp DNA Stool Mini Kit following manufacturers protocol. DNA was amplified using two different sets of primers amplifying a portion of COI. We used a common set of primers, ZBJ (ZBJ-ArtF1c and ZBJ-ArtR2c) (Zeale et al., 2011), and a recently designed set, ANML (LCO1490/CO1-CFMRa) (Jusino et al., 2019). Both sets of primers included Illumina TruSeq sequences creating fusion primers (Glenn et al., 2019b; Kieran et al., 2017).

Amplicon library construction was performed in two steps, first initial replicate PCRs using fusion primers was performed for each primer set using Promega Go Taq (Promega). Reaction consisted of 3 μl Buffer, 0.3 μl of 10 mM dNTPs, 0.1 μl Taq, 7.16 μl molecular grade water, 0.24 μl 10mg/mL BSA, 0.6 μl of 5μM forward primer, 0.6 μl of 5 μM reverse primer, and 3 μl of DNA at the following thermocycler conditions: 94°C for 2 min, followed by 5 cycles at 94°C for 1 min, 45°C for 1.5 min, 72°C for 1.5 min, then 35 cycles at  94°C for 1 min, 50°C for 1.5 min, 72°C for 1 min and a final extension at 72°C for 5 min. After verification on a 1.5% agarose gel, replicates were combined separately for each set of primers (ZBJ and ANML). Amplicons for each primer set were purified using SPRI-beads (Thermo-Scientific, Waltham, MA) in a 1.7:1 ratio and reconstituted in 20 μl TLE.

Second step PCR used iTru primers(Glenn et al., 2019a) with Illumina TruSeqHT compatible 8 nt indexed primers, which when used combinatorially provide 96 unique index combinations. We performed replicate 25 μl reaction of KAPA HiFi HotStart Kits using 2.5 μl of 5x Buffer, 0.25 μl of 10 mM dNTPs, 0.1 μl HotStart, 4.15 μl molecular grade water, 1 μl of 5 μM forward primer, 1 μl 5 μM reverse primer, and 3 μl of 1^st round PCR product using the following thermocycler conditions: 98°C for 3 min, followed by 6 cycles at 98°C for 30 s, 60°C for 30 s, 72°C for 30 s and a final extension at 72°C for 5 min. Product was pooled in equal volume, cleaned with SPRI-beads in a 1.7:1 ratio, and pooled with other uniquely indexed samples for sequencing on an Illumina MiSeq using a PE300 kit (Illumina, San Diego, CA) at the Georgia Genomics and Bioinformatics Core (http://dna.uga.edu).

Bioinformatic analyses

Sequencing data outer indices were demultiplexed using bcl2fastq (Illumina v1.8.4) to identify the amplicon pool (ZBJ, ANML). Each amplicon pool data was then demultiplexed by internal indices and primers removed using Mr. Demuxy v1.2.0 (https://pypi.python.org/pypi/Mr_Demuxy/1.2.0). Reads were then imported into Geneious v10.1 (Kearse et al. 2012, Biomatters Limited, NJ), set as paired-reads with an expected insert size of 150bp, and quality trimmed with a score of 0.001. Paired-end sequences were then merged using FLASH v1.2.9 plugin (Magoc and Salzberg, 2011) excluding reads below 80bp in length. Data were then exported as FASTQ files and imported into Qiime2 v2018.11 (Bolyen et al., 2018).

Full records of all Insecta COI sequences available on the Barcode of Life Database (BOLD) from the United States, Canada, and Mexico were downloaded on 24 April 2019. We used custom code in R (cite) to format the data for the taxonomy and reference sequence database files used by Qiime. Additional reference sequences for taxa of interest not found on BOLD were downloaded from Genbank as FASTA files and added to our custom taxonomy and reference sequence database files. Complete reference sequence database was duplicated and imported into Geneious where sequences were trimmed independently using ANML and Zeale primers to limit references to amplicon region of interest.

 Custom taxonomy and trimmed reference sequence databases formatted for Qiime were imported into QIIME2 v2018.11 along with preprocessed FASTQ files from Geneious. Chimeric sequences were removed, and de novo clustered at 99% identity using vsearch (Rognes et al., 2016). We then used QIIME2’s feature-classifier plugin with the classify-consensus-blast to identify sequences against out custom reference database at 97% and 99% identify with minimum consensus threshold of 60% and 80%. These resulting feature tables were collapsed at the species level and exported in BIOM format for further analysis in R. We tested other analysis methods using DADA2 (Callahan et al., 2016) to denoise and remove chimeras, followed by clustering with vsearch, and using classify-sklearn to identify sequences.

Statistical Analysis

We first examined the species composition of diets using heat trees (Citation). We developed a heat tree of subfamilies consumed, and individual heat trees for the 4 most dominate orders consumed. Additionally, we quantified the number of families, subfamilies, and genera in each sample. Finally, we classified each prey species as pest or non-pest and determined the number of pest subfamilies and genera present in each sample. Furthermore, we calculated the proportion of total reads attributable to pests and non-pests in each sample as a crude measure of relative abundance.