MATERIALS AND
METHODS
Study
Area
The study area encompassed c. 3500 hectares of conservation property in
western Montana (www.mpgranch.com; 46°41’ N, 114°00’ W). Historic
management practices include cattle grazing, logging, and agriculture.
Current conservation strategies include restoring native grasslands and
shrublands, primarily through weed control, seeding and planting
efforts, wildlife management, and irrigation. Sampling occurred in
mid-elevation forest/grassland, mid-elevation forest, floodplain forest,
mid-elevation sagebrush, mid-elevation sagebrush/woodland, and
mid-elevation shrubland/grassland plant communities. Elevation ranged
from approximately 970 meters in floodplain areas to around 1650 m in
higher elevation forests.
Sample collection and
processing
We collected fecal samples from NAIs May through September during 2017
and 2018. We captured bats monthly after evening emergence in mist nets
set over dry land, streams, and ponds. We placed bats in individual
paper bags to collect their fecal pellets. Six additional bat species
occur in the study area but were excluded from this study due to low
sample sizes. We collected fresh fecal samples from Flammulated Owls,
Common Poorwills, and Common Nighthawks captured in mist nets on or near
breeding territories. Common Nighthawks and Common Poorwills were also
sampled opportunistically near roads or at known nest or roost sites via
hand nets. We placed all fecal samples in vials containing ethanol in
the field and stored them in a freezer at -20° C until further
processing. We labeled samples by the plant community in which they were
collected, though sampling location does not always equate to plant
community used while foraging. A subsample of Common Nighthawks,
Flammulated Owls, and Common Poorwills were also fitted with GPS
tracking devices to gather data on home and foraging ranges. Telemetry
data for Common Nighthawks indicated foraging ranges up to 400 ha,
whereas Common Poorwills ranged between 0.5 and 3.0 ha and Flammulated
Owls generally foraged in <1.0 ha (Supplementary Table 1). The
daily foraging range of bats sampled varies between < 1 km for
Long-eared Myotis to > 4.4 km for Big Brown Bats (Maxell,
2015). Based on observational and telemetry data, home and foraging
ranges overlapped for all species.
The Canadian Centre for DNA Barcoding (CCDB) performed all DNA
extractions, amplification, and sequencing. DNA extraction and PCR
amplification followed CCDB protocols as described in Moran et al.,
(2019). Samples were incubated overnight in a lysis buffer, concentrated
and dried by centrifugation, and finally eluted using a Tris-HCl elution
buffer. The CCDB also processed negative controls in parallel with
samples to ensure that contamination did not occur. The cytochrome C
oxidase 1 (CO1) region was amplified from each sample using the
arthropod-specific primers, ZBJ-ArtF1c_t1 and ZBJ-ArtR2_t1 (Zeale et
al., 2011), as described previously (Moran et al., 2019; Prosser &
Hebert, 2017). Following amplification, samples were pooled and
purified. The CCDB performed sequencing on an Ion Torrent PGM following
standard protocols (Prosser & Hebert, 2017).
Constructing a DNA barcode library
from local
Arthropoda
In 2017 and 2018, we collected nocturnal insects monthly May-August
using mercury vapor and black lights placed in front of a white sheet
and an aerial flight-intercept trap at sites across our study area. In
2019, we expanded insect sampling to include bulk samples collected
weekly over 13 weeks (May-August) from flight-intercept, pitfall, and
yellow and blue pan traps. We sent samples to the CCDB for sequencing
and identification (deWaard et al., 2019; Ratnasingham & Hebert, 2007).
Technicians at CCDB counted total insect abundance by order, weighed
biomass, and collected tissue samples from several representatives of
each morphospecies. All records are publicly available on BOLD under the
datasets MPG and MPGR with photos of specimens to aid in future
identification. The resulting local arthropod DNA barcode library
consisted of 56,191 Arthropoda specimens collected May-September from
2017 to 2019 at 48 sites within our study area. Nearly all (99.5%) of
the specimen sequences were assigned to order, 92.8% to family, 58.0%
to genus, and 24.4% to species. A total of 52,033 of the sequences
gained Barcode Index Numbers (BINs) in the Barcode of Life Database
(BOLD), comprising 6,080 total unique BINs. This effort added 1,529
previously undocumented arthropod records to BOLD, and represented 38
orders, 383 families, 1,810 genera, and 1,740 total species. Dominant
orders represented in the final database included Diptera, Hymenoptera,
Hemiptera, Lepidoptera, and Coleoptera (Supplementary Figure 1).
Data
analysis
We processed demultiplexed sequences using QIIME2 version 2020.2 (Bolyen
et al., 2018). We removed all primers prior to analysis using the
cutadapt plugin (Martin, 2011) and denoised sequences using the DADA2
denoise-pyro plugin (Callahan et al., 2016). DADA2 is sensitive to
single base-pair differences among sequences and produces unique
“amplicon sequence variants” (ASVs). The median base pair quality
score for all sequences was maintained above 25. Denoised sequences
shorter than 100 bp were removed from analyses. We clustered sequences
into operational taxonomic units (OTUs) based on a 97% sequence
similarity threshold (Vamos et al., 2017), using the VSEARCH plugin
(Rognes et al., 2016). We removed sequences only occurring in a single
sample or that were represented by fewer than 0.001% of sequences to
limit artifactual sequences.
We determined taxonomic assignments using our local DNA barcode library
and the BLAST plugin within QIIME2, with a coverage value of 0.7 and
sequential percent matching identities of 100, 99, 98, and 97%. If
taxonomy could not be assigned to our local database using these
parameters, we used a global COI database compiled from BOLD and GenBank
and a pretrained RDP classifier (Porter & Hajibabaei, 2018; Wang et
al., 2007). The BOLD accession ID associated with each taxonomic
identification is indicated where available. We verified all taxonomic
identifications based on the plausibility that they may occur within or
nearby the study area. All sequences not matching to Arthropoda using
either the global COI database or the local database were removed from
further analyses. We rarefied samples at 500 sequences, which was
sufficient to adequately characterize the species within each sample
(Supplementary Figure 2). In total, 77% of OTUs recovered from NAI
fecal samples matched 97% or greater with locally collected specimens,
while 23% were assigned taxonomy using the RDP classifier.
All statistical analyses were conducted in RStudio Version 1.1.453 using
R version 3.6.0, (R Core Team, 2018). We used linear regression to
assess differences in diversity metrics based on species, plant
community, and collection month. We used a two-way Anova (car package,
Fox and Weisberg, 2011) with a type II sum of squares for unbalanced
data to test significance of diversity metrics. Where significant, the
emmeans package was used for pairwise analyses of diversity metrics
between each species (Lenth et al., 2021). To determine insect taxa
significantly associated with different NAI species or foraging
strategies (i.e., aerial hawking vs. sit-and-wait strategies, or
echolocation vs. visual hunting), we used the indicspecies package
(Cáceres & Legendre, 2009).
To determine resource partitioning among NAI diet composition, we
analyzed both relative read abundance and presence/absence data.
Presence/absence data are considered a more conservative option in
insectivore fecal analyses (Jusino et al., 2019). However,
presence/absence data can also overestimate the importance of prey
consumed in small quantities, and it is generally thought that relative
read abundances provide more accurate population-level data (Deagle et
al., 2019). Even so, we chose to analyze both relative read abundance
and presence/absence data and found similar results. We performed all
compositional comparisons on either Bray-Curtis distances of Hellinger
transformed relative read abundances, or Raup-Crick transformed
presence/absence data using the vegan package (Oksanen et al., 2019). We
assessed differences in diet dispersion (distance from the mean) among
species, as significant differences in group dispersions can lead to
false positives when testing differences among groups. Because we saw no
differences in data dispersion, we used the Permanova function to assess
differences in NAI diet composition at the OTU level. For pairwise
compositional analyses, we used the pairwiseAdonis2 function in the
‘pairwiseAdonis’ package (Arbizu, 2017/2021), and adjusted for multiple
comparisons (Benjamini & Hochberg, 1995). For each NAI, we also
calculated diet turnover among samples using the vegetarian package
(Charney & Record, 2012). Diet turnover within each NAI species was
calculated based on Shannon beta diversity where zero equals no
difference between samples and one represents completely different
samples.