Ethics Statement

The Daubenton's bats were sampled with active sampling licenses from the Finnish National Animal Experiment Board (ESLH-YM-2007-01055) and Centre for Economic Development, Transport and the Environment (LOS-2007-L-182-254). Care was taken to handle and trap bats according to ethical guidelines presented by the National Animal Experiment Board.

The study area and its characteristics

The sampling took place on the mainland and islands in the northern part of the Archipelago Sea, situated in southwestern Finland (Figure 1). The area has a large river outlet of fresh water and the salinity is very low (typically less than 5 psu) [35], which affects the aquatic as well as terrestrial species composition of the area, including insects with aquatic life stages. Due to the nutrient runoff from rivers, the Archipelago Sea is eutrophic [36]. The coastline in the northern part of the Archipelago Sea is covered by a wide (up to approx. 30 m) bed of Common reed (Phragmites australis (Cav.) Trin. ex Steud). Uncluttered areas in the vicinity of the reed beds are primary feeding areas for Daubenton’s bats in the area.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. The study area, sampling sites, and regions.

River (Nautelankoski rapid), North (Korjaustelakka and Tamminiemi), West (Pakinainen), Central (Kukonpää), East (Rauhala and Brinkhall), and South (Seili and Hemsundet). Land is shaded grey and water is white.

https://doi.org/10.1371/journal.pone.0082168.g001

Bat trapping and fecal sampling

Bats were caught with a combination of mist nets and harp trap. A Sussex Autobat acoustic lure, which plays species-specific ultrasound social calls, was placed in the center of harp trap to attract the bats [37]. This multi-trap combination was placed across the flying corridor of bats commuting between roosts and foraging areas.

All fecal sampling was conducted in the field. Upon catching, bats were kept in clean cloth bags until handling for identification and sampling. Caught bats were identified to species and sex and banded. Age was determined from bone structures, especially from the joint of phalanges. The bones of juveniles are not fully ossified. With the onset of flight, the phalanges are almost completely ossified, but growth plates near the joints are apparent as a light cartilaginous gap. Fecal pellets were collected either from laundered, individual handling bags or directly from the bat while handling. In the first case (bags) the fecal droppings were collected using disposable tweezers immediately when the bat was removed from the bag. In the second case (direct collecting) the fecal droppings were collected immediately upon excretion from the surroundings of anus using disposable tweezers. The pellets were then stored in 1.5 ml, 94 % ethanol-filled Eppendorf tubes, which were individually labeled and subsequently stored at +4. See Table 1 for detailed information of the sampled bat individuals.

DNA extraction and PCR setup

Altogether, we sampled 37 fecal pellets from 34 bat individuals for fecal DNA extractions (Table S1). We took one fecal pellet at a time from each individual to extract DNA. In order to avoid contamination, fecal pellets were processed in the Herbarium laboratory at the University of Turku, where only plant specimens have been processed so far. Total fecal DNA was extracted using QIAamp DNA Stool Mini Kit (Qiagen, catalog number 51504). Negative control extractions containing all the chemicals but no fecal pellets were performed alongside each batch of extractions to monitor for contamination of the extraction chemicals. The fecal DNA concentrations ranged from less than 1 to over 20 nanograms per microliter. The PCR reactions were also prepared in the insect-DNA-free laboratory and thermocycling was later performed in the TEGlab facilities (Laboratory of Genetics, University of Turku (for a review of methods, see 4). The following PCR conditions were applied to our DNA templates: 1x KAPA HiFi Hotstart ReadyMix (Kapa Biosystems, catalog number KK2602), 10 µM each dNTP, 8 % (final concentration) BSA, 0.5 µM each primer (ZBJ1c and ZBJ2c after [3]), and 1–2 µL fecal DNA extraction. PCR grade water was added up to 25 µL reaction volume. Negative control reactions were performed alongside each batch of PCR to monitor for contamination of the chemicals used. The thermal profile included a 5 min initial denaturation at 95 °C followed by 40 cycles of a 30 s denaturation at 95 °C, a 30 s annealing at 50 °C, and a 30 s elongation at 72 °C. Final elongation was conducted at 72 °C for 5 minutes. The PCR success was visually inspected under UV light using a 2 % agarose gel stained with Midori Green Advanced (Nippon Genetics Europe, catalog number MG04). The successful reactions were purified using NucleoMag 96 PCR Kit (Macherey-Nagel, catalog number 744100.1). The final elution volume was 25 µL.

Ion Torrent PGM library preparation

We employed a time and cost effect method to pool individually tagged DNA samples into two sequencing libraries [38]. The first sequencing library included 15 individuals, and the second 23 individuals. Three individuals (band numbers 2015, 2026, and 2038) were included in both libraries. Libraries were prepared according to [38] with the following modifications to their protocol to make it more suitable for Ion Torrent PGM (steps as in [38]):

Step 1: The adapter mix was prepared using multiplex identifier (MID) tagged adapters. Adapters were designed in order to include an individual MID and to match Ion Torrent specific priming sites. The molar concentration of the adapter mix was 20 µM for each adapter. We used a 200 bp long PCR product as positive control DNA template through steps 1 to 26.

Step 4: We performed the blunt-end reaction with half of the volume, i.e., we added 10 µL of blunt-end master mix to 25 µL of sample. The next purification using SPRI beads was adjusted accordingly to ensure a template:bead ratio of 1:1.8. In the post-ligation purification steps template:bead ratio was 1:1.

Step 20: To verify adapter ligation success, the SPRI-purified MID-tagged positive control library and a subset of MID-tagged library aliquots were separated side by side with non-MID-tagged templates on 2% agarose gel stained with Midori Green Advanced electrophoresis (45 min, 95 V) and visualized using UV light.

Step 26: For a subset of samples, we measured DNA concentration and fragment size distribution using Bioanalyzer (Agilent Technologies, Inc.) and DNA concentration using Qubit Fluorometer (Life Technologies Corporation). The results from Qubit were consistent with Bioanalyzer results. We omitted Bioanalyzer run for each individual samples and quantified the libraries using the Qubit Fluorometer. Individual amplicon libraries were subsequently mixed in equimolar ratios.

The stock library was re-amplified in two separate reactions with following setup: 5 µL pooled library stock was added to a master mix consisting of 5 U of Herculase II polymerase (Agilent technologies, catalog number 600677), 1x Herculase II reaction buffer, 25 mM each dNTP, 10 µM each primer, and added PCR grade water up to 50 µL. We designed the re-amplification primers based on the sequence of the IonTorrent specific priming sites and thus generated millions of MID-tagged copies including binding sites necessary for subsequent sequencing utilizing IonTorrent technology. The thermocycling profile included a 30 s denaturation at 98 °C followed by 15 cycles consisting of a 20 s denaturation at 98 °C, a 30 s annealing at 64 °C, and a 30 s elongation at 72 °C. Final elongation was conducted at 72 °C for 5 minutes.

To clean the re-amplified library, size-selection was done by separating the entire library using 2% Size-Select Agarose E-Gel and E-Gel Electrophoresis System (Life Technologies, catalog numbers G6610-02 and G6500) following the manufacturer’s instructions. The DNA concentration of amplified library pool stock was again measured with Qubit Fluorometer. The library pool stock was then diluted to a final concentration of 26 pM. For template preparation, an 18 µL aliquot of the library dilution (approximately 2.8 x 10⁸ molecules) was transferred into the sequencing reaction set up.

Emulsion PCR and Ion Torrent PGM Sequencing were carried out on two 314 chips (Life Sciences, catalog number 4462923) according to the manufacture’s protocol (Publication Part Number: 4471974 Rev. C).

Processing of Ion Torrent PGM reads

The resulting reads were binned and renamed by MID, i.e., the original individual samples using the software Geneious Pro 6.1 [39]. Then the original arthropod-specific primers (ZBJ-1c and ZBJ-2c; [3]) were trimmed off and all reads were trimmed for poor quality parts using 0.05 error probability limit and then the reads shorter than 100 bp were excluded using Geneious Pro 6.1. Thus, only reads containing both the Ion Torrent adapter with a MID and the original primer were passed on for further analysis. Subsequent analyses were carried out using super computers at the IT Center for Science (Espoo, Finland, www.csc.fi). Sequences were assigned to species using BLASTN 2.2.25+ algorithm [40,41] against two databases, built the following way: DB-1) a public database consisting of all the arthropod COI Barcode of Life Data systems [42] and Genbank [34], and DB-2) all the sequences as in database one plus sequences from National Finnish Barcode of Life (FinBOL) campaign and sequences previously generated by ourselves [43]. We blasted all individuals separately against DB-2 and we also blasted all the query sequences against both databases using super-computing servers at Finnish Grid Infrastructure (FGI, http://www.csc.fi/english/collaboration/projects/fgi).

Analyzing the BLAST results

The BLAST output was imported into the program MEGAN [44] to get an overview of the results. We used default parameters, except for following settings: Min support=1, Min score=100, Min complexity=0. We utilized a built-in comparison tool in MEGAN for database comparison. As recommended by [8,9], we discarded all hits below 98% identity. We also discarded hits with e-value over 1e-20. We further examined the hits one by one to see if there were reads matching more than one species in the database. In these cases we followed these steps:

1. 1. Select the hit with the longest alignment length,
2. 2. Select a match to an adult specimen over juvenile,
3. 3. Check the database entry for errors and omit wrong entry

Following these simple steps we were able to reliably select species identity for 23253 reads. We calculated prey species accumulation curve [45] based on samples (= individual bat fecal pellets) to analyze the total dietary species richness using program EstimateS (version 9) using 1000 runs for sample order randomization ([46], see also [47] and references there on). The expected species curve confidence limits where calculated by rarefaction after [48]. Non-parametric Chao 1 (using classic formula setting) [49] and ICE [50,51] richness estimators were calculated and plotted against individual samples to see the behavior of traditional estimators for bats as diversity sampling units. The effective number of species (true diversity of order 1) was calculated after [52]. The dominance-diversity relationship of dietary species was analyzed by plotting species occurrences per species against species rank (from most abundant to least abundant). The differences in species diversity between regions could not be statistically analyzed due to small sample size and biased sampling. To visually demonstrate the prey species assemblages we drew a species cloud based on the identification and frequency of all the prey species at each region.

Overview of Ion Torrent PGM data

Two sequencing runs produced 12.1 Mb of raw data of which 8.67 Mb showed a quality of Q20 or higher (see Supporting information S1 and S2 for details). The raw reads are deposited in the Sequence Read Archive (SRA). The SRA accession codes for each individual are listed in Table 1. Total number of reads was 110142, of which 101266 reads (91.9%) included perfectly matching MID. The rest of the reads (8.1%) either included MID-tag with too many mismatches to be reliably identified or did not include a MID at all. In the latter case, the read usually appeared to be an artifact (for example primer-dimer), that is, not a real DNA sequence, which was confirmed by BLASTing these sequences against Genbank’s nucleotide database. After trimming away reads shorter than 100 bp we ended up with 41575 high quality reads. All control reactions proved to be free of contamination.

Comparison of BLAST results from public and local databases

We compared the blast results for the selected public databases and the local database first at order level. Most of the orders were equally represented in both blast outputs with Diptera being the most frequent order (Figure 2. A). For the order Hemiptera there were more hits in the public database whereas for the order Psocoptera there were more hits in the local database. We then compared genus and species level hits for the main prey group, the family Chironomidae of the order Diptera (Figure 2.B and C). Some of the genera and species are totally missing from the public database, or in some cases the sequences are available, but with incomplete identification information (Figure 2.B and C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Comparison of results from public and private databases.

The bar chart comparison of public databases (Genbank, BOLD) and public+private databases (Genbank, BOLD + FinBOL sequences) in identification of reads originating from Daubenton’s bats fecal samples: A. Comparing results at order level for all arthropods. B. Comparing results at genus level for family Chironomidae. C. Comparing results at species level for family Chironomidae.

https://doi.org/10.1371/journal.pone.0082168.g002

BLAST results

The BLAST to the local database returned matches to 40769 (98%) of the query reads. Of these, 23253 query reads were qualified for further analysis, that is, at least 98% identity to database sequences and e-value < 1E-20. Altogether, 128 prey species from the Daubenton’s bats fecal DNA were identified, of which 79 with 100 % similarity to the reference database (Table S1). The effective number of species for the whole dataset was 99. The reads matched database sequences from three arthropod classes: most of the hits were assigned to class Insecta, and a small proportion to Arachnida and Collembola (Table 2). More than 70 % of all the species belonged to two proportionally most abundant orders: Diptera and Lepidoptera (Table 2).

The insect order Diptera was found in 94% of samples, Trichoptera 69%, and Lepidoptera 63% (Figure 3). More than every fifth (21 %) dietary species belonged to one single dipteran family Chironomidae (Figure 4), although there were a total of 49 families identified (Table S1). The most common prey family was the dipteran family Chironomidae, which was found in 89% of all the samples (Table 3). Other frequent families included the caddisfly family Leptoceridae, the moth family Tortricidae and another dipteran family, Tipulidae (Table 3). The most common prey species were non-biting midges (Diptera: Chironomidae) Microtendipes pedellus (50%), Glyptotendipes cauliginellus (44%), and Procladius ferrugineus (41%). The number of different prey species identified in each sample varied from two to twenty four. On average, one bat individual consumed 10 prey species.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. The most frequent arthropod orders in the Daubenton’s bats diet.

The abundance of all the arthropod orders occurring in the Daubenton’s bats diet. The bar chart shows the frequency of the order in all the sampled bat individuals (n = 34).

https://doi.org/10.1371/journal.pone.0082168.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Most abundant families in the Daubenton’s bats diet.

Pie chart shows proportions of the dietary families (n=49) obtained from Daubenton’s bat fecal DNA extractions amplified using arthropod-specific mtDNA .

https://doi.org/10.1371/journal.pone.0082168.g004

Prey species diversity

The species accumulation curve shows, that when resampled, every bat individual adds more species to the prey species list (Figure 5). The number of uniques (species occurring only once) and duplicates (species occurring only twice) are also high. Both non-parametric diversity estimators (Chao1 and ICE) predict overall species diversity higher than observed. The prey species rank abundance plot clearly indicates, that there are a few very dominant species in the diet, and many species occurring only rarely (Figure 6). The species identification cloud indicates that insects in the order Diptera are among the most consumed species in all regions (Figure 7). A trichopteran species, the mass-emerging Oecetis ochraea, is the most frequent single species only at a single location (South). The built-in comparison tool used for database comparison in MEGAN uses NCBI taxonomy (available at http://www.ncbi.nlm.nih.gov/taxonomy) as such, and for this reason our final prey species (Table S1) list differs from the chironomid genera and species list used in database comparison analysis (Figure 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. The prey species accumulation curve based on individual bat samples.

The y-axis shows number of samples (t) accumulated. S_expected, expected number of prey species in t pooled samples with 95% confidence intervals. S_obs, observed number of prey species in t pooled samples. Uniques and duplicates, species observed only once or twice (respectively) in all the pooled samples. Chao1, incidence-based (presence/absence) non-parametric species richness estimator. ICE, Incidence coverage-based non-parametric species richness estimator. See text for references.

https://doi.org/10.1371/journal.pone.0082168.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Dietary species’ dominance-diversity relationship.

Log₁₀ percentage of species occurrences are plotted against species rank (from most abundant to least abundant).

https://doi.org/10.1371/journal.pone.0082168.g006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Daubenton’s bats prey species by frequency and taxonomic order.

Species identification cloud of prey species identified from Daubenton’s bats fecal DNA samples on six sampling regions. Font size represents the proportional occurrence of each prey species in all the samples sequenced. See Figure 1 for regions on the map.

https://doi.org/10.1371/journal.pone.0082168.g007

Molecular methods for dietary studies

We analysed 34 fecal pellets collected straight from bats while handling them. The library preparation protocol [38] was slightly modified in order to analyze fecal samples. The two sequencing runs did not produce as much data as we expected, but the data are nevertheless useful in demonstrating the method. Based on the preliminary results obtained in this study, we can conclude that Daubenton’s bats in the Archipelago Sea feed mainly on dipteran insects, especially of the family Chironomidae. The present conclusion is in agreement with previous studies on the diet of the Daubenton’s bat, in which almost all fecal pellets contained traces of chironomids [10,13,14,25]. The second most frequent prey groups are caddisflies and small moths, which were found to be equally frequent. This finding contradicts earlier studies, which seemed to find that moths are only marginal prey items for Daubenton’s bat, occurrences varying from very low to medium [26–28,30–32]. This may be due to variation in sampling habitats and different analysis methods, which vary from anecdotal evidence or direct observation of hunting bats [29], prey remain analysis [28], morphological analysis of fecal remains [10,13,14,26,27,30–32], and analysis of stomach contents [53,54]. It seems that while morphological identification of highly degraded remains yields similar results as DNA-based methods, the resolution of identification is not as good [3,12]. Based on the above we come to conclusion that modern sequencing technology, here Ion Torrent PGM, is suitable for accurate dietary studies.

Importance of local DNA barcode library for molecular approach

In comparison to numerous earlier studies on insectivorous bat diets [3,9–14,17,25–33,53–56], the application of novel molecular techniques and an extensive DNA barcoding reference library in this study helped to truly identify dietary objects at the species level. Furthermore, with additions to the library, we can also apply the method to other bat species. As an improvement for future studies, more individuals should be sampled throughout the field season for temporal variation from several sites, and technical replicates such as several fecal pellets from a single individual and PCR replicates need to be analyzed.

Daubenton’s bats dietary summary and taxonomic report of prey groups

Overall prey species richness for Daubenton’s bat was high: 128 identified species in 49 families and 10 arthropod orders. Although high, the number of prey species is in accordance with earlier studies on Daubenton’s bats diet [10,14,25]. All of the main prey species are aquatic or related to water habitat [57]. Many of the other prey species are also, for example many moths occurring only once or twice, using shore vegetation as a food source in their larval stage. The dietary composition of main prey groups was found to be similar across all three years, making the results more robust. Although the species richness of total diet was high, the effective number of species (true diversity) using species occurrence frequencies was lower (98) because of the few dominant species and many rare species. According to the resampling of the individual bat diets, the overall species richness is likely to be higher than observed, but there is little doubt that the observed main dietary groups are already revealed by this study. Overall, there were many arthropod groups and many species present in the diet. Here we will take an overview of taxonomic groups, and discuss the commonness and importance of various arthropod groups in Daubenton’s bats diet.