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.