Methods

Sampling design

Sampling took place in the Iberian Peninsula, focusing on two sympatric regions in the north where both Myotis escalerai and Myotis crypticus are found (La Rioja-Soria and southern Cantabria), and two allopatric regions: the south (Andalucia: Jaen and Granada), where onlyM. escalerai is found, and the north Atlantic coast (northern Cantabria), where only M. crypticus is found. Additionally, we sampled a single sympatric swarming site in Catalunya, where the two species use the same cave during the autumn mating season (Fig. 1). Within each region, 9-24 locations were sampled based on suitable habitat and accessibility, using monofilament mist nets and a harp trap placed over water sources, forest paths and cave entrances. The sampling period extended from June to September 2017, for a total of 68 sampling nights (Supplementary Table S1 for list of sampling locations). Captured bats were kept in individual cotton bags for up to 1 hour. We collected faecal samples from the cotton bags for diet analysis, and biopsy punches (3 mm) from the wing membrane of the bats to confirm species identification. Dropping samples and wing biopsies were stored in absolute and 70% ethanol, respectively. Bat sampling was carried out under local permits and ethical approval from the University of Southampton (study ID: 26627).
We sampled the arthropod community in bat sampling locations using vegetation sweeping (Barclay, 1991; Swift & Racey, 2002) to assess bat prey selectivity relative to prey availability (Jones, 1990; Kunz, 2009). We chose vegetation sweeping because of the expected low flight and gleaning behaviour of the species based on their morphology and echolocation calls, and the foraging behaviour of the morphologically similar better studied member of the cryptic species complex,Myotis nattereri (de Jong, 1995). During bat sampling nights, we set linear transects in each habitat type in the sampling location and swept the vegetation along each transect. After 10 sweeps, we collected the capture into a plastic bag and moved five steps further without sampling to increase spatial representativeness. Each sampling unit of 10 sweeps and five steps forward was repeated 5-10 times until the capture size was considered representative (> 100 individuals). Transect length ranged between 30 and 80 m. Arthropod specimens captured were separated from vegetation remains in the field and stored in 70% ethanol.

DNA extractions and species confirmation

Bat species identity was confirmed in the Estación Biológica de Doñana Laboratory of Molecular Ecology (LEM, EBD-CSIC, Spain). DNA was extracted from wing biopsy punches through precipitation with isopropanol. Part of the hypervariable region of the mtDNA control region was amplified using the primers CSBC-F 5’-CCTCTTAAATAAGACATCTCGATGG-3’ (Wilkinson & Chapman, 1992) and HV2-Mna-R 5’-ATGCGTGCGTGTGTAATGTC-3’ (Garcia-Mudarra, Ibanez, & Juste, 2020). Species specific differential amplification patterns for this primer set were used to confirm species identity through gel electrophoresis (Garcia-Mudarra et al., 2020).
DNA was extracted from all bat dropping samples using the Qiagen DNA stool mini kit, following the protocol in Zeale et al. (2011). A total of 43 sweeping samples from 23 locations with at least three bat dropping samples were selected (Supplementary Table S1). From those, all arthropod individuals (N = 8366) were first identified morphologically to taxonomic order. Subsequently, whole specimens, if smaller than a drosophila, or a specimen part (leg or head) if larger, were separated out, dried and pooled together for DNA extraction. Arthropod DNA was extracted using the NucleoSpin DNA Insect kit with up to 35 mg of sample dry weight in each tube. Larger samples were split into several tubes. The following modifications were applied to the kit extraction protocol: In steep 2, vortex for 20 minutes in the MN Bead Tube Holder on a Vortex-Genie at maximum speed; after steep 3, pipette 550uL of clean supernatant in to a new 2mL Eppendorf, centrifuge again at 13,000RPM for 2 minutes and continue with steep 4; in steep 6, centrifuge for 3 minutes; in steep 7, add 50uL of ddH20 and incubate for 3 minutes.

High Throughput Sequencing

Both dropping and sweeping samples were sequenced in the Bart’s and the London Genome Centre, London, UK. DNA extracts were checked for quality and concentration on a TapeStation D1000. Two sets of primers were used together in order to reduce primer taxonomic bias (Alberdi, Aizpurua, Gilbert, & Bohmann, 2018), especially given the high diversity of prey types expected in the diet, ZBJ (Zeale, Butlin, Barker, Lees, & Jones, 2011) Forward: ZBJ-ArtF1c 5’-AGATATTGGAACWTTATATTTTATTTTTGG-3’ and Reverse: ZBJ-ArtR2c 5’-WACTAATCAATTWCCAAATCCTCC-3’, and ANML (Jusino et al., 2019) Forward: LCO1490 5’-GGTCAACAAATCATAAAGATATTGG-3’ and Reverse:  CO1-CFMRa 5’-GGWACTAATCAATTTCCAAATCC-3’. For the ZBJ amplicon each 15μl PCR reaction used 7.5μl of Multiplex PCR mastermix (QIAGEN, Germany), 0.25 μL of each primer (10 μM), 5 μL H2O and 2 μL template DNA. Negative and positive controls were included in PCR reactions and later sequencing. The thermal cycling protocol was as follows: 95 ̵̊C for 15 min, 34 cycles of 94 ̵̊C for 40s, 40 ̵̊C for 1 min, 72 ̵̊C for 30s, followed by a final extension of 72 ̵̊C for 5 min. ANML regions were amplified in 15ul reactions following published protocols (Jusino et al., 2019). All products were visualised on a 1.5% agarose gel. Products were tagged using Fluidigm barcodes and checked on a TapeStation D1000 before pooling and sequencing on an Illumina MiSeq using paired end (2 x 250 bp) chemistry (Illumina, San Francisco, USA). We used two technical PCR replicates to reduce biases associated with PCR stochasticity. This led to each sample being sequenced four times (combination of two primer and two PCR replicates).

Bioinformatics

Sequencing runs were merged using USEARCH (Edgar, 2010) and primers and adaptors removed using cutadapt (Martin, 2011). Sequences were processed on the mBRAVE platform (http://www.boldsystems.org/bin) (Ratnasingham & Hebert, 2007) setting the following parameters: Min QV = 0 qv, Min Length = 100 bp, Max Bases with Low QV (<20) = 75%, Max Bases with Ultra Low QV (<10) = 75%, ID Distance Threshold = 1.5%, Exclude from OTU Threshold = 3%, Minimum OTU Size = 1, OTU Threshold = 2%. Sequences were compared with the BOLD reference libraries SYS-CRLINSECTA and SYS-CRLNONINSECTARTH to established Barcode Index Numbers (BINs). BINs are a type of Operational Taxonomic Unit (OTU) integrated in the BOLD system with advantages over traditional OTUs, such as being unique and stable (Ratnasingham & Hebert, 2013). We used BIN identity as a proxy of taxonomic prey item unit.
After obtaining BIN (prey item) composition per sample and run, we removed singletons, i.e. BINs that only had a single read per run and sample, because they are likely to be PCR or sequencing errors (Alberdi et al., 2018). We established the threshold for the minimum number of reads per sample to retain an identification based on comparing composition similarity between molecular and morphological ID from sweeping samples and any prey present in blanks. Additionally, we controlled for potential contamination during extractions by removing the BINs present in extraction blanks from other samples in the same extraction run if present with a high number of reads (>100 in any blank sample) or with less than 10 times more reads in samples than in blank. We controlled for contamination from the sequencing process by removing BINs present in sequencing blanks from all the samples from the same primer following the same criteria as above.
We used two alternative approaches to combine data from PCR replicates. In the first, the additive criteria, BIN composition from both PCR replicates of each sample were added together. In the second, the conservative approach, only BINs that appeared in both PCR replicates were considered (Alberdi et al., 2018). Under this second criteria only samples in which the four runs contained a minimum number of reads (>100) were considered given that a failed sequencing run in a sample would lead to a null composition for both of the PCR replicates of a primer (52 samples, Supplementary Table S3). Finally, we combined taxa recovered from both primers to obtain the prey composition per dropping sample for downstream ecological analysis. Duplicated BINs in the same sample coming from different PCR replicates or primers were removed. A flow chart describing the methods is shown in Supplementary Fig. S1. Ecological results from both approaches were very similar, thus we present results based on the additive approach (See Supplementary Fig. S9 for diet based on the conservative approach).

Characterising the diets of the two bat species

The contribution of different elements to the diet for a set of samples was quantified using weighted Percent of Occurrence (wPOO), which measures the relative occurrence of diet elements (prey items/OTUs/BINs) in a set of samples considering first their relative proportion per sample (Deagle et al., 2019). For example, a prey item found in a sample with 9 other prey items will be interpreted to contribute in the diet 1/10 of what it would if it was the only prey item present. Contributions to diet based on the two other commonly used metrics: Percent of Occurrence (POO) and Relative Read Abundance (RRA) (Deagle et al., 2019), are shown in supplementary materials (Supplementary Fig. S2). We tested for differences in the number of BINs per sample between bat species for each of the orders that constitute at least 10% of the diet of either one of the bat species (Araneae, Diptera, Lepidoptera, and Hemiptera), using negative binomial generalised linear models (GLMs; in R) to fit data structure based on the distribution of model residuals. We measured order level and prey species level (BIN-level) diet composition overlap between bat species using Pianka’s measure of niche overlap (Ojk) (R package: EcoSimR (Gotelli, Hart, & Ellison, 2015)). We tested with an ANOSIM test (R package vegan (Oksanen et al., 2019)) whether Jaccard distance in BIN composition was greater between than within bat species. The ANOSIM statistic R is based on the difference of mean ranks between and within groups, with a range between -1 and +1. A value of zero indicates that the group does not explain compositional differences. We visualised ordination of samples depending on their BIN composition with Non-Metric Multidimensional Scaling (NMDS, R package vegan: Oksanen et al., 2019). We calculated Levins’ (1968) standardised measure of niche breadth (BA) at the prey species (BIN) level for each bat species.

Functional diet assessment

Prey items were classified based on the literature (outlined in Supplementary Table S2) and an expert entomological taxonomist into three functional categories: non-volant, not actively-volant, nocturnally volant. Categorisation depended on their mobility and type of activity, reflecting their likelihood of being captured by gleaning or aerial hawking (Supplementary Data file S1, Supplementary Table S2). The categorisation was done at family or finer taxonomic level by checking the literature for data on daily activity patterns of each family and presence in nocturnal light traps (Supplementary Table S2 for criteria used). The non-volant category included wingless arthropod groups (Araneae, Isopoda and wingless insects such as some members of Blattodea, Orthoptera). The not actively-volant category included those able to fly but unlikely to have been captured by the bat through aerial hawking because they are not active fliers, either at night (diurnal Diptera), or not active fliers in general (e.g. Hemiptera, some Blattodea, Orthoptera and Coleoptera). The nocturnally volant category comprised arthropods with aerial and nocturnal activity, and therefore likely to be captured by aerial hawking (e.g. non-Ropalocera Lepidoptera, nocturnal Diptera, Neuroptera, Ephemeroptera, Trichoptera). This classification represents the likelihood of being captured by gleaning or aerial hawking rather than direct inference of the capture mode because nocturnally active aerial prey can also be captured by gleaning when resting on vegetation and not active nocturnal fliers could also be captured in the air (e.g. ballooning in spiders). Once all prey items were classified into functional groups, we obtained the functional diet of both bat species using weighted percent of occurrence (wPOO), and compared the percentage of not nocturnally volant (including both non-volant and not actively volant categories) per sample between bat species using a linear model.

Trophic niche overlap in allopatry vs sympatry across spatial scales

Locations from Andalucía (Mediterranean climate) and northern Cantabria (Atlantic climate) were classified as regionally allopatric. Locations from La Rioja and southern Cantabria (climatically Mediterranean to sub-Atlantic) as regionally sympatric (based on data from Razgour et al., (2019) and EBD records). At the fine-scale within the sympatric regions, we classified locations as locally sympatric or allopatric depending on whether they were within 3 km of records of the other species based on a conservative estimation of the home-range diameter of the better studied cryptic congener M. nattereri (Boye & Dietz, 2005). The swarming location in Catalunya was removed from the fine-scale analysis because bats gather in swarming sites from distances of up to 60 km from their colonies for the purpose of breeding rather than foraging (Rivers, Butlin, & Altringham, 2005), and therefore it is unclear whether those individuals forage in sympatric areas (Supplementary Table S1 for sampling locations and their broad and fine-scale sympatry category).
To identify differential use of certain prey orders and functional groups, we tested separately for allopatric and sympatric locations whether 1) the number of BINs per sample for each of the main arthropod orders, and 2) the percentage of not nocturnally volant functional groups differed between bat species. We used negative binomial zero inflated GLMs and a linear model respectively. We run separate models for broad and fine spatial scales given that both fine scale sympatry and allopatry treatments are within regional sympatry. Cases where resource (prey order or functional group) use was different between bat species when sympatric but not when allopatric were regarded as evidence of resource partitioning. We measured prey species (BIN) level niche overlap (Ojk) between bat species in sympatry and in allopatry, and tested, using null models (R package ecosimR) whether overlap was lower or higher than random in sympatry versus allopatry. We tested whether Ojk differed between sympatric and allopatric locations by pooling the diet composition of each bat species per location and measuring Ojk between pairs of locations. At the regional scale we used a Gaussian Hurdle model due to the high density of zeroes in overlap values. At the local scale we used a linear model with log transformed values of Ojk to meet assumptions of normal distribution. All statistical analysis was carried out in R (R core team, 2019).

Prey consumption relative to availability

For each location, we quantified the relative availability of each arthropod order and functional group using weighted percent of occurrence (wPOO) after pooling together sweeping samples from the different habitats. Similarly, we obtained bat diet composition (wPOO) per location of each arthropod order and functional group by pooling diet composition of all individual bats. Then, we subtracted from the bat diet wPOO the prey availability wPOO to obtain prey arthropod and functional group selection per location. A higher proportion of a given arthropod order in the diet than in sweeping samples indicates the bats may be preferentially consuming this resource, based on prey availability at the sampled strata. This analysis is based on assumptions of sampling representativeness (see discussion for detailed overview of methodological limitations).

Testing primer performance and representation of the DNA metabarcoding approach

For each arthropod order we described the number and proportion of BINs identified by each primer. Morphological identification of the arthropod communities allowed us to compare the performance of the primers and metabarcoding approaches. We compared the presence of orders in each sweeping sample based on molecular and morphological identification to determine whether metabarcoding offers a good estimation of arthropod community composition.