Bats play an important role in pest control, seed dispersal and pollination (Kunz et al. 2011, Baroja et al. 2019), but are threatened by the loss of foraging habitat and insect declines (Tiede et al. 2020). Determining the prey spectrum and identifying suitable foraging areas is thus key to making conservation decisions in habitats occupied by insectivorous bat species. In addition, dietary variation in bats has been shown to be dependent on landscape and agricultural practices (Aizpurua et al. 2018); therefore, changes in land use can lead to the loss of foraging habitat, as well as source habitats suitable for the whole life cycle of prey species. The importance of extending bat conservation areas beyond directly-used hunting grounds to include prey source habitats has been highlighted (Arrizabalaga-Escudero et al. 2015) and particular attention should be paid to cover all habitats for arthropod prey species that have several life stages, as many arthropods are known to show ontogenetic habitat shifts (Carr et al. 2020, Kirse et al. 2021Kirse et al. 2021, Carr et al. 2020This is also highly relevant for bat protection under the EU Habitats Directive (92/43/EEC), where protected sites should include not only roosts, but habitats necessary for the entire life cycle of bats.

N. leisleri is a bat of 13-18 g and a wingspan of 26-32 cm (Dietz et al. 2009). This bat is a small member of the genus Nyctalus with a western Palaearctic range (Europe and north-west Africa) with scattered records in the eastern Palaearctic (Pakistan, Afghanistan, the Himalayas) (Juste and Paunović 2016). In late summer and autumn, large parts of the European population migrate south to spend the winter under milder conditions in southern France, Spain and Italy (Boston et al. 2021). Mating happens during late summer on or prior to the beginning of migration. After hibernation and spring-migration, the bats return to their summer habitats and the females form nursery colonies, mainly in tree cavities (Boston et al. 2021), but also in buildings behind wall sidings. They have long, narrow wings adapted for fast flying speeds and for catching insects during flight in the open airspace above the canopy and water bodies, as well as near street-lights and forest edges (Boston et al. 2021). The Conservation status on the IUCN Red List, as well as the European Red List is LC: Least Concern (European Commission et al. 2007, Juste and Paunović 2016).

Current threats to N. leisleri include: 1) the reduction in insect abundance due to increased pesticide use; 2) changes in land use leading to the disappearance of fallow land, permanent grassland, hedges and margins, causing the loss of insect-rich habitats; 3) habitat loss due to the draining of wetlands and water bodies in forests and open countryside; 4) habitat degradation through reduction of natural or semi-natural forests; 5) the loss of old trees with high roost potential; 6) renovation work on buildings leading to loss of roosts and roosting opportunities and 7) wind-energy development due to direct collision with rotor blades especially during migration (Boston et al. 2021).

Bats prey upon a wide variety of arthropod species of various sizes, diurnal and nocturnal and flying or non-flying, but many studies find that Lepidoptera, Diptera and Coleoptera represent the dominant prey orders (Alberdi et al. 2012, Arrizabalaga-Escudero et al. 2015, Baroja et al. 2019, Tiede et al. 2020, Alberdi et al. 2020. Leisler's bat is an insectivorous aerial hawker known to catch insects in flight in the open air space above the forest canopy and close to forest edges, some of which are caught in swarms (Waters et al. 1995, Kaňuch et al. 2005). Radio tracking shows that N. leisleri commutes to foraging sites up to 13.4 km away from the roost (Shiel et al. 2006a). Based on visual analysis of taxonomically-informative remains found in faecal pellets of N. leisleri, the most frequently encountered prey were from the insect orders Lepidoptera, Diptera, Coleoptera and Trichoptera (Beck 1995, Kaňuch et al. 2005, Waters et al. 2006). The presence of Trichoptera in the diet indicates that N. leisleri hunts over water bodies (Beck 1995, Vaughan 1997, Shiel et al. 2006b). However, soft-bodied prey species may be underestimated using this method and taxonomically-important parts of prey species may be missing (Beck 1995, Alberdi et al. 2012). While most dietary studies focus on the analysis of faecal pellets, no high-resolution molecular studies of N. leisleri diet exist to date. High resolution dietary analyses in bats can help answer a wide variety of ecological questions, such as the relationship between dietary niche breadth and spatial distribution (Alberdi et al. 2020) or help shed light on bat foraging ecology (Alberdi et al. 2012).

Direct observation of feeding is generally very difficult in nocturnal bats and visual identification of prey arthropod remains in bat faeces does not generally result in taxonomic identification of the prey below order or family level (Alberdi et al. 2012). DNA metabarcoding has revolutionised the field of dietary analysis, revealing much higher prey diversity than previously recorded through morphological analysis in many taxa, ranging from bats (Clare et al. 2009, Zeale et al. 2010, Bohmann et al. 2011, Tiede et al. 2020), to fish (Leray et al. 2013, Jakubavičiūtė et al. 2017, Bourlat et al. 2021) and even invertebrates (Waldner and Traugott 2012, Sint et al. 2015). In general, primers covering a wide range of taxa are used for gut content analysis, where shorter fragments of 100 - 250 bp of the mitochondrial cytochrome c oxidase I (COI) gene are generally sufficient to provide taxonomic resolution at the species level (Meusnier et al. 2008, Zeale et al. 2010). This is advantageous for amplification from dietary remains, where DNA is expected to be highly degraded. Molecular methods enable prey identification up to species level, which can be very useful for addressing questions relating to specific features in the prey species, such as wingspan, whether they are tympanate species or whether they are nocturnal or diurnal species. In addition, foraging habitats can be inferred from the consumed prey species using finer scale taxonomic resolution of prey allowing the broadening of the scope of ecological studies on bats (Alberdi et al. 2012), such as their role in providing important ecosystem services (e.g. pest control) in agricultural landscapes (Baroja et al. 2019).

In this study, we sampled bat droppings (guano) at the roost of N. leisleri in a natural forest reserve in North Rhine-Westphalia in Germany during March to September 2017. There were three defined objectives to our study. First, to provide a high-resolution analysis of arthropod prey species and seasonal trends in the insectivorous bat species N. leisleri. Second, to compare the performance of two fragments of varying lengths for COI, the 313 bp ‘mini barcode’ (mlCOIintF combined with dgHCO2198, hereafter COImldg) (Meyer 2003, Leray et al. 2013) and the 157 bp fragment (ZBJ-ArtF1c combined with ZBJ-ArtR2c, herafter COIArt) (Zeale et al. 2010) and a 110 bp region of the mitochondrial 16S gene (IN16STK-1FW combined with IN16STK-1Rv, herafter 16S) (Kartzinel and Pringle 2015) in the identification of arthropod species from bat faeces. Third, to identify the ecological requirements for some of the prey species identified and, based on this, make recommendations for the conservation of N. leisleri prey habitats.

The data underlying this study have been submitted to the NCBI SRA archive under accession number PRJNA752700.

Research site

Sampling was carried out in the EU Natura 2000 site 'Waldreservat Kottenforst' (DE5308303), located near Bonn, Germany between 180 and 200 m above sea level. The forest has an area of 2450 ha and is dominated by sub-atlantic and medio-atlantic oak (Quercus robur, Quercus petraea) and oak-hornbeam (Quercus sp., Carpinus betulus) forest, partially with varying admixture of beech (Fagus sylvatica). Hydromorphic soils with high water tables provide numerous small water bodies. The forest is managed for wood production and includes "wilderness areas", corresponding to unmanaged stands. To the north and east, the forest borders urban areas of the city of Bonn and the highly urbanised and industrialised Rhine valley. To the north-west, the Kottenforst is connected to the Waldville forested area and west to southeast, the forest borders agricultural areas.

Capture and radio-tracking of bats

Bats were radio-tagged between 2014 and 2016 to find roost trees within the Natura 2000 site, in order to align management decisions with nature conservation goals. Bats were caught with mist nets at small water bodies and at potential foraging sites. Suitable animals (female, not pregnant, no injuries, minimum average weight) were equipped with transmitters (Telemetrie-Service Dessau) weighing between 0.3-0.5 g (Permission number: RSK 67.1-1.03.20-18/14-M). Tags were attached between the scapulae using surgical glue. Roost trees were tracked down the day after tagging and checked for the presence of bats for the next 10 - 14 days until transmitters fell off or the transmitter battery was consumed. For this study, bat guano was sampled at the main roost tree of Nyctalus leisleri, which was identified in 2014 in a woodpecker cavity at a height of 9 m. This roost tree was used constantly during the summer months for several consecutive years.

Bat guano sampling

A novel type of guano trap was installed beneath the roost entrance (Fig. 1). The guano trap is a lightweight rectangular frame of PVC-pipes (25 mm in diameter) with a mosquito-net (mesh width 1.4 mm), attached to the trunk of the roost tree at 3.5 - 4 m height, approximately 3 to 6 m below the roost entrance. The catchment-area of the trap is 2.2 m².

Returning bats show pre-dawn swarming behaviour at occupied roost-trees (Naďo and Kaňuch 2013, Stanton 2016). The bats fly in close circles around the roost tree, land and leave in close proximity to the roost entrance (‘touch and go’) (Kaňuch 2007, Naďo and Kaňuch 2013) and stick guano pellets on the trunk close to the roost entrance. It is assumed that dawn swarming is part of group decision-making and day roost selection of the colony (Schöner et al. 2010, Chaverri et al. 2018). Since bats have a rather low digestive efficiency (Barclay et al. 1991) and a short retention time of prey remains in the digestive system (Roswag et al. 2012), they frequently drop faeces pellets which can be caught by the guano trap.

The trap was checked after nights when swarming was likely to happen. It was assumed that good conditions for swarming were warm nights, with no wind and no rain in the second half of the night. During unfavourable weather conditions, the trap was checked on a regular basis every 2-3 days to remove leaves, small twigs and other debris. Pellets were collected from the net, stored in 15 ml sterile sampling tubes and dried with silica gel and/or stored in 2-propanol. This sampling method is non-invasive and bats do not have to be caught or disturbed to collect faeces for dietary analyses. However, pellets collected can originate from different individuals and possibly even different bat species, due to interspecific swarming behaviour at the roost. Therefore, species identity of the bats was checked using both COI and 16S primers upon library sequencing and data analysis. After species identity check, nine samples confirmed to be from N. leisleri were included for further analyses. All samples of bat faeces collected and analysed in this study are detailed in Suppl. material 1.

DNA extraction and amplicon library preparation from bat faeces

DNA was extracted from bat guano pellets using the Zymo Quick-DNA™ Fecal/Soil Microbe Midiprep kit, following the manufacturer’s instructions. Guano pellets stored in ethanol were first dried and approximately 40 mg of guano were subsampled from the pellet pool for DNA extraction. All samples (stored in ethanol and silica) were extracted in three replicates including a negative control consisting of sterile water. DNA concentration was measured using the Quantus™ Fluorometer with the QuantiFluor® dsDNA System (Promega). All samples were diluted to 2 ng/µl.

PCR amplification was performed using a 2-step PCR approach. The first PCR was carried out in a total volume of 15 µl per replicate, using 7.5 µl of Q5 Hot Start High ‐ Fidelity 2X Master Mix (NEB), 0.5 µl of each primer (10 µM), 0.5 µl of Bovine Serum Albumin (Thermo Fisher Scientific), 5 µl Sigma H₂O and 1 µl of DNA. PCR1 conditions involved denaturation at 98°C for 2 min, followed by 20 cycles at 98°C for 40 sec, 50°C for 40 sec and 72°C for 30 sec and a final extension step at 72°C for 3 min. DNA extraction negative controls and PCR negative controls (water) were included for every PCR reaction.

PCR1 products were purified using the HT ExoSAP-IT^TM (Thermo Fisher Scientific), with 4 µl ExoSAP for 15 µl PCR 1 product, following the manufacturer’s protocol.

In a second PCR step, the Illumina index adaptors were attached to the purified PCR1 product, which was split into two tubes, each with 7 µl of PCR1 product. Amplifications were carried out in a total volume of 25 µl with 12.5 µl of Q5 Hot Start High ‐ Fidelity 2X Master Mix (NEB), 1.2 µl of each primer (10 µM), 1 µl of Bovine Serum Albumin (Thermo Fisher Scientific), 2 µl Sigma H₂O and 7 µl PCR 1 product. PCR2 conditions involved a denaturation at 98°C for 2 min, followed by 20 cycles at 98°C for 40 sec, 55°C for 30 sec and 72°C for 30 sec and a final extension step at 72°C for 3 min.

All replicate PCR2 products were pooled, visualised by electrophoresis on a 2% agarose gel (120 V 20 min, 150 V 40 min, 450 mA, 150 W) and purified with the QIAquick gel extraction kit (Qiagen). All purified PCR products were then diluted to the same concentration (3 ng/µl) and pooled into two amplicon libraries. One library comprised the 313 bp COI fragment and the second library the 157 bp COI and 110 bp 16S fragments.

Sequencing

The purified amplicon library pools were sequenced on four runs on the Illumina MiSeq platform (2 x 300 bp) using the v.2 Chemistry at the Centre for Genomic Research (CGR, Liverpool University).

Bioinformatic methods

Data sequenced at the Centre for Genomic Research (Liverpool, UK) had already undergone a first quality check. The raw fastq files were trimmed for the presence of Illumina adapter sequences using Cutadapt version 1.2.1 (Martin 2011). Sequences were further trimmed using Sickle version 1.200 (Joshi and Fass 2011) with a minimum window quality score of 20. Reads shorter than 20 bp after trimming were removed. Only sequences passing this first quality check were available for download from the CGR server. The downloaded sequences were checked for the presence of the three primer pairs using Cutadapt version 2.10 with Python 3.6.10 (Martin 2011) with the following settings: maximum error rate (-e): 0.1, minimum overlap (-O): 20, minimum sequence length (-m): 150. Each primer-pair dataset was analysed separately. Only sequences with both forward and reverse primers were retained for further analysis. The primers were removed from the sequences before being uploaded to the QIIME2 pipeline (Bolyen et al. 2019). For denoising using Dada2 (Callahan et al. 2016), sequences were truncated to the following lengths: forward and reverse reads of COI mIdg to 175 bp and 170 bp, respectively; forward and reverse reads of COIArt to 216 bp and 169 bp, respectively; forward and reverse reads of 16S to 126 bp and 126 bp, respectively.

Depending on marker, two different reference databases were used. COI sequences were blasted against the German Barcode Of Life (GBOL) database, downloaded from (https://doi.org/10.20363/gbol-20210128) on 29 January 2020 using the following settings: (a) 'query coverage high-scoring sequence pair percent' (-qcov_hsp_perc) was set to 90, meaning that a sequence was reported as match when 90% of the query formed an alignment with an entry of the reference file; (b) minimum percent identity (-perc_identity) was set to 97, requiring the reference and query sequence to match by at least 97% to be reported as a match. The format of the output file was customised using the –outfmt settings ‘6 qseqid sseqid pident’. Taxonomic assignment with the GBOL database yielded 36 arthropod species for mldg and 241 arthropod species for COIArt in the nine guano samples from N. leisleri (Suppl. material 2).

The mitochondrial 16S sequences were blasted against a customised 16S reference database downloaded from NCBI GenBank on (29 December 2020). The following search parameters were applied:16S[All Fields] AND (animals[filter] AND is_nuccore[filter] AND mitochondrion[filter] AND ("100"[SLEN] : "1000"[SLEN])). Taxonomic assignment with the GenBank database using a 97% blastID yielded 119 arthropod species for the nine guano samples from N. leisleri (Suppl. material 2).

For the ecological analyses, ASV tables converted to a presence/absence matrix were uploaded into R studio (version 1.4.1106; R version 4.0.4.). For statistical analysis, nine guano samples, assigned uniquely to N. leisleri with a 100% Blast match, were analysed (KF01-01, KF01-02, KF01-03, KF01-06, KF01-07, KF01-08, KF01-09, KF01-10, KF01-11). Venn Diagrams were prepared using the package VennDiagram (version 1.6.20) (Chen and Boutros 2011). Assessed community composition depending on markers was visualised using the R packages ggplot2 (version 3.3.3., Wickham (2016)) and RColorBrewer (version 1.1-2, Neuwirth (2011)). For visualisation of N. leisleri diet over time, the R packages ggplot2 (version 3.3.3.) and ggpubr (version 0.4.0.Kassambara (2018)) were used. As it has previously been shown that diet analyses based on presence/absence data are very conservative and sometimes overestimate the food consumed in small quantities, we show in parallel analyses based on relative read abundance (RRA) calculated using the formula of Deagle et al. (2018). RRA was calculated for all arthropod and lepidopteran taxa (Suppl. material 3, Suppl. material 4). For calculation of the RRA per sample, the number of reads assigned to each arthropod taxon within a sample was divided by the sum of the of reads for all arthropod taxa in that sample and multiplied by 100 ((number_of_reads_per_taxon_and_sample / total_number_of_reads_per_sample)*100). For calculation of the total RRA in all samples combined, the sum of the reads assigned to a taxon across samples was divided by the sum of the reads for all arthropod taxa in the dataset (all samples combined) and multiplied by 100 ((Number_of_reads_per_taxon / total_number_of_reads)*100). Corresponding R code for figures 2, 3 and 4 can be found in the supplementary materials (Suppl. material 5).

Taxonomy assignment for the 16S mitochondrial marker was carried out against the NCBI database. This database is more incomplete than the GBOL database for the arthropods, especially arthropod species from Germany. These assignments are more likely to represent the best available match and are, therefore, likely biased at the species level. Based on this, we excluded taxonomic assignments with 16S from the analysis based on RRA.

Molecular identification of species occupying the roosts

All bat species identification was confirmed using both COI and 16S primers upon library sequencing and data analysis, since guano provides a non-invasive source of DNA that includes information from the bat as well as dietary items, parasites and pathogens (Swift et al. 2018). When taxonomy assignment for each guano sample retrieved exclusively N. leisleri with 100% BLAST match, we could confirm that the pellet originates from N. leisleri and the sample was included for further analysis (Table 1, samples marked in bold).

Denoising with Dada2 yielded 1519 ASVs (amplicon sequence variants) for COI mIdg, 1107 ASVs for COIArt and 565 ASVs for 16S for the samples included in our analysis (KF01-01, KF01-02, KF01-03, KF01-06, KF01-07, KF01-08, KF01-09, KF01-10, KF01-11). Taxonomic assignment with the GBOL database yielded 36 arthropod species for mIdg and 241 species for COIArt for the nine guano samples (Suppl. material 2). Taxonomic assignment with the GenBank database using a 97% blastID yielded 119 arthropod species for 16S (Suppl. material 2).

Molecular identification of species occupying the roosts

High-resolution analysis of arthropod prey species in N. leisleri guano and comparison of different markers

Seasonal trends in prey consumed

Most abundant species in the bat guano

Guano samples provide a non-invasive source of DNA that includes information from the bat, but also dietary items, parasites, and pathogens (Swift et al. 2018). In this study, our analyses confirmed the presence of five different bat species (Myotis nattereri, Nyctalus leisleri, Plecotus auritus, Myotis mystacinus and Myotis bechsteinii), with multiple species sometimes found at the same roost. The latter can occur since bats eavesdrop on other bat species to find roosts in forests (Jones 2008, Ruczyński et al. 2009). In particular, Schöner et al. 2010 showed that M. bechsteinii, M. nattereri and P. auritus can approach bat boxes with played-back bat calls. We assume that bat calls, even from other species, were a cue for N. leisleri to approach and find possible roosts within our study area.

Dietary analyses and seasonal trends

Ecological requirements of the prey species

Recommendations for the management of N. leisleri

In this study, we show that metabarcoding has the capacity to improve the quality and resolution of ecological data, such as diet and prey data, which can be a turning point for the success of habitat and conservation management measures. From the N. leisleri prey data obtained, we derive a set of key recommendations for N. leisleri habitat and conservation management:

• Preserve rivers, streams and ponds (Ephemeroptera, Trichoptera, Chironomidae).
• Preserve a varied landscape with drier and wetter herbaceous meadows and forests, as well as wetlands around the bat habitat.
• Avoid extensive mowing of meadows. Avoid complete mowing at once, keep or install rotational systems which leave part of the grassland or forest borders unmown every year to support overwintering and survival of above-ground immature insect stages of many Lepidoptera and Diptera.
• Avoid simultaneous yearly cutting of forest track margins with tall herb stands.
• Avoid spraying nettles (Urtica) with herbicides. This is an important larval host plant for many of the Lepidoptera species that serve as food for N. leisleri.
• Preserve trees with woodpecker cavities where N. leisleri roosts, as well as standing dead wood.
• Many of the target prey species (Lepidoptera, Ephemeroptera) are attracted by artificial light sources. Use insect friendly street lighting, to reduce the impact on the insects and the bats (Rowse et al. 2015).
• Monitor changes in prey phenology. Climate induced phenological shifts could affect prey availability at the time of highest energy requirement for the bats (birth and lactation).

This study was funded by the German Federal Ministry of Education and Research, through the German Barcode of Life project (GBOLII, FKZ01LI1501) and the European Commission Life Forests-waterworlds project (LIFE13 NAT/DE/000147). We thank the Biological Station Bonn-Rhein-Erft, Karina Jungmann, Heinz Schumacher and Rolf Mörtter for help with the fieldwork and Lisa Apfelbacher for the illustration.

SJB and VGF designed and supervised the project. MK designed the guano traps and MK, RM and AS conducted the fieldwork. KL conducted the molecular lab work. HG and AK conducted the bioinformatic analysis and SJB, AK, VGF, MK and ME the data interpretation. SJB wrote the paper with contributions from VGF, AK, KL, MK, ME. All authors approved the final manuscript.

The authors have declared that no competing interests exist.