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.