Introduction

Traditionally, diets of free-ranging animals have been determined by direct observation of feeding bouts or food remains, and by microhistological inspection of gut contents or faeces. Even if these approaches have provided much of the currently available dietary information on wildlife, they have some limitations: either cannot be applied to elusive animals [1], are limited by ethical reasons [2,3], or the results heavily depend upon the researchers’ skills [4–6] and the remains left by preys (reviewed in [7]). In order to overcome these difficulties, researchers have innovatively adopted a wide array of molecular approaches with varying success [1, 8–12]. Nonetheless, they have not been capable of determining the diets’ components at the species level.

Conversely, DNA based dietary studies allow the examination of the range and diversity of prey taken by generalist predators/consumers [1, 13], the identification of bulk samples even within highly degraded samples such as faeces, gut contents or regurgitates [13], the processing of DNA from many different consumed species [14] and, using non-invasive procedures, also the diet characterization of elusive species [15, 16]. Consequently, DNA metabarcoding and High Throughput Sequencing (HTS) techniques have recently become common-use in dietary studies (e.g. [16–25]), offering great detection ability and identification of consumed prey and plants to the species level. However, achieved results may be biased due to the sequencing technology [26], the genetic marker choice, the performance of the primers’ and PCR amplification, and laboratory workflows or bioinformatic analyses (e.g. [27–29]).

Cytochrome oxidase gene subunit I (COI) has become the most commonly used marker region in DNA metabarcoding diet studies (e.g. [19,24]) because it has the most extensive information in genomic databases (BOLD, GenBank). Besides, it includes very short fragments of DNA–"mini-barcodes” (representative fragments of COI)–easily recoverable from degraded samples (e.g. [24, 30–32]). Amplification of such fragments is currently carried out using different primer sets with varying success. Most studies on the trophic ecology of insectivorous bats have relied on the primers proposed by Zeale [19] (e.g. [17, 33–40], but see [41, 42]). The Zeale primers have proven to be highly successful to amplify DNA from Diptera and Lepidoptera [43, 44], and therefore have been widely used. Nevertheless, primer choice is a critical step in metabarcoding studies as the resulting list of prey may be unwittingly biased [17] due to the particular taxonomic affinity of the selected primers (e.g. [28, 29, 45]).

Alberdi et al. [45] highlighted that the use of multiple primers targeting the same taxonomic group reduces the effect of each primer sets’ biases and increments the taxonomic coverage, obtaining a more complete view of the diet of the predator. Accordingly, Esnaola et al. [46] showed that five primer sets targeting sections of varying lengths within the COI region performed differently when amplifying faecal DNA of the generalist insectivorous Pyrenean desman Galemys pyrenaicus. Even if most of the primers used were able to identify the most common arthropod prey taxa consumed, the differences regarding less abundant prey groups were considerable, and hence the diet composition depended on the chosen primer sets. Esnaola et al. [46] found that the combination of two primer sets was the most successful, namely the Zeale primers mentioned above [24] and a second primer set modified following Gillet et al. [31], targeting a shorter 133 bp mini-COI sequence. These two primer sets have different length and degeneration levels, and allegedly best reveal the prey range of generalist predators [46].

We aim to apply the aforementioned pair of primer sets to the molecular diet analysis of an ensemble of horseshoe bats (family Rhinolophidae) composed by Rhinolophus hipposideros (Bechstein, 1800), R. ferrumequinum (Schreber, 1774) and R. euryale Blasius, 1853. These three species have the greatest distribution range and broadest overlap in Europe. So far, the diet of R. hipposideros and R. ferrumequinum have been analysed using only Gillet primers by Galan et al. [42], while that of R. euryale has been characterized either with Zeale primers [33,34] or with Gillet primers [42]. Based on those studies we expect that R. euryale will mainly prey upon Lepidoptera, R. ferrumequinum upon Lepidoptera, Diptera and Coleoptera and finally R. hipposideros upon Lepidoptera, Diptera and Neuroptera. In accordance with previous research [27,46], we expect that the results of each primer set will be different and complementary, coming with a more precise view of their diet, since the identification of prey species will be primer-dependent. Secondly, we want to evaluate the combination of the aforementioned two primers sets to characterize the diet overlap of the three horseshoe bat species, which show varying preferences for moths. Andreas et al. [47] studied the niche partitioning of the aforementioned horseshoe bats based on microscopic identification and showed that their trophic niche overlap was considerably low. Thus, we aim to see how the results on the overlap of their trophic niches reflect the choice of the primer set.

Materials and methods

DNA extraction, PCR amplification, library preparation and sequencing

Individual faecal samples of 10–40 mg were used for DNA extraction with the DNeasy PowerSoil Kit (Qiagen, Valencia, CA), following the manufacturer steps. Extracted DNA was PCR-amplified twice using to different primer sets, targeting different mini-COI segments of the mitochondrial DNA cytochrome c oxidase subunit I barcode region (COI): Zeale primers (ZBJ-ArtF1c and ZBJ-ArtR2c) [24] were used to amplify a 157 bp section, and Gillet primers (modified LepF1 and EPT-long-univR, following [31] to amplify another 133bp section. Both amplifications were performed using QIAGEN Multiplex PCR Kit (Qiagen Iberia, S.L. Madrid) in 25 μl PCR reactions. Each reaction contained 2.5 μl Buffer 10X, 1.5 μl MgCl₂ 50mM, 0.5 μl nNTPs 25mM and 0.125 μl of taq polymerase. In the case of Zeale primers, 0.6 μl of each primer (forward and reverse), 17.175 μl deionised water and 2 μl sample DNA were added. With Gillet primers, 0.75 μl of each primer, 14.875 μl deionised water and 4 μl sample DNA were added. Each primer set had its own PCR program, modified from the reference to the used reactive. Thermocycler conditions for Zeale primers were: 95°C– 15 min; 50 cycles of 94°C– 30 sec, 52°C– 30 sec, 72°C– 30 sec; 72°C– 6 min (modified from [51]). For Gillet primers we used: 95°C– 15 min; 40 cycles of 94°C– 30 sec, 45°C– 45 sec, 72°C– 30 sec; 72°C– 10 min [31]. For the library preparation, each sample was tagged with a unique combination of Multiplex Identifier primers (MID) [52]. PCR outputs were sequenced by Ion Torrent sequencing platform, one run making above one million reads.

Results

Diet of horseshoe bats

We successfully extracted and amplified DNA from faeces of 24 R. ferrumequinum individuals, 18 R. euryale and 31 R. hipposideros, obtaining above one million sequence reads (Table 1). 309 OTUs were then built and 135 of them were assigned to potential prey species consumed by bats.

Download:

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

Table 1. Results obtained in the different steps of bioinformatic analyses with each of the primer sets.

https://doi.org/10.1371/journal.pone.0220081.t001

We identified 62 prey species of R. ferrumequinum: 34 lepidopterans, 17 dipterans, 7 coleopterans, 2 neuropterans and 1 trichopteran. Lepidoptera and Diptera were the most frequently consumed, followed by Coleoptera (Tables A, B and C in S1 Supporting Information File). Among the most frequently consumed species Pharmacis fusconebulosa (Hepialidae; FOO = 46%) prevailed among Lepidoptera and Rhipidia maculata (Limoniidae; FOO = 71%) and Tipula maxima (Tipulidae; FOO = 29%) among Diptera. Within Coleoptera, the most consumed were the elaterid Stenagostus rhombeus (FOO = 46%) and the cerambycids Arhopalus rusticus and Prionus coriarius (FOO = 42% and 17% respectively), completed with scarabeids Aphodius sp. and Serica brunnea (FOO = 29% and 17% respectively).

Out of the 81 prey species identified in faeces of R. euryale, 61 were lepidopterans, 10 dipterans, 5 neuropterans, 2 ephemeropterans, 1 trichopteran, 1 hemipteran and 1 hymenopteran. Lepidoptera was the most frequently consumed order followed by Diptera (Tables A, B and C in S1 Supporting Information File), although the FOO of most of them was less than 20%. The exceptions were the noctuids Capsula sparganii, Cosmia trapezina and Lycophotia porphyrea (FOO > 28%), the geometrid Idaea biselata (FOO = 56%) and the limonid Austrolimnophila ochracea (FOO = 94%).

Finally, 73 prey species were identified for R. hipposideros, including 33 lepidopterans, 28 dipterans, 3 hemipterans, 3 neuropterans, 2 coleopterans, 1 hymenopteran, 1 trichopteran, 1 spider and 1 psocopteran. Among them, Diptera were the most frequently consumed, followed in descending order by Lepidoptera, Ephemeroptera, Trichoptera and Neuroptera (Tables A, B and C in S1 Supporting Information File). Among Diptera, lesser horseshoe bats mostly preyed upon limonids Austrolimnophila ochracea (FOO = 100%), Rhipidia maculata (FOO = 58%), Neolimonia dumetorum (FOO = 52%), Limonia nubeculosa (FOO = 29%) and Dicranomyia modesta (FOO = 26%), and the tipulid Tipula helvola (FOO = 55%). They also consistently preyed upon the neuropterans Hemerobious humulinus and Wesmaelius nervosus (FOO > 19%). Noteworthy, the occurrence of most moth species was below 3 with a maximum of 8 occurrences of the autostichid Anania hortulata.

As a whole, 12 prey species have been identified in the faeces of the three predators: 6 were lepidopterans (Acronicta rumicis, Cyclophora punctaria, Idaea degeneraria, Anaplectoides prasina, Noctua sp. and Udea ferrugalis), 5 dipterans (Rhipidia maculata, Austrolimnophila ochracea, Limonia nubeculosa, Neolimonia dumetorum and Tipula helvola) and one neuropteran (Hemerobius humulinus).

Performance of primers

Gillet primers yielded the highest numbers of reads, whereas Zeale ones got the highest numbers of either primary OTUs, positively identified OTUs, occurrences of prey and prey species identified (Table 1). Moreover, some of the OTUs built from Gillet primers were identified as belonging to algae and mammal species (4.41% of the total taxa), and so they must be considered as environmental pollution instead of "potential prey" consumed by bats.

We first tested that there was not significant geographical effect of the two sampling sites in the diet (F_(1,68) = 1.031; R² = 0.132; p = 0.37). Therefore, the location variable was not considered in further analyses. The difference between species diets is significant for the whole data set (F_(2,135) = 8.277; R² = 0.092; p = 0.001), but also for the results obtained with each of the primer sets by their own (Gillet:F_(2,70) = 10.466; R² = 0.230; p = 0.001; Zeale: F_(2,65) = 4,772; R² = 0.128; p = 0.001). The primer choice significantly affected the resulting diet composition (F_(1,135) = 14.438; R² = 0.082; p = 0.001; Fig 1). Consequently, the sum of the partial results enlarged the entire prey species list. Of the total species identified as potential prey 40% have been identified with Gillet’s set and 74.8% with Zeale’s, i.e., only 21 out of the 135 (15.5%) potential prey species have been amplified by both primer sets. Anyway, we can see that both primer and bat species affect the list of consumed prey, with a slightly higher explanation of the variation in the case of the bat species. The interaction of primer and species also shows a significant difference among the results, even if it explains less variation than primers and species on their own (F_(2,135) = 4.841; R² = 0.055; p = 0.001). The complete lists of potential prey species identified with each primer set is are included in Tables A, B and C; sequences of all the OTUs built are available in Table D, all of them in S1 Supporting Information File.

Download:

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

Fig 1. NMDS ordination of samples.

Stress = 0.1997; k = 2; non-metric fit R² = 0.96. Dots represent prey species and colours different primer sets (Red: Zeale; Green: Gillet). More distant dots indicate more different prey composition of samples. Individual bat samples are represented as grey triangles.

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

There were big qualitative and quantitative differences among primers at a broader taxonomic level as well (F_(1,135) = 61.157; R² = 0.239; p = 0.001) (Fig 2, Tables 2 and 3, and Table A, in S1 Supporting Information File). Thus, Zeale primers were able to identify five orders of potential prey–namely Lepidoptera, Diptera, Neuroptera, Hemiptera and Psocoptera–whereas OTUs yielded from Gillet primers were assigned to species of fifteen orders: namely Lepidoptera, Diptera, Coleoptera, Trichoptera, Neuroptera, Ephemeroptera, Hemiptera, Hymenoptera and Araneae for prey species, as well as Mucorales, Artiodactyla, Primates and Chiroptera for environmental DNA. Besides, Zeale primers yielded more occurrences than Gillet ones (Table 1). The three predator species were identified with Gillet primers in all samples but in 2 R. ferrumequinum.

Download:

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

Fig 2. Results of the three bats’ diet obtained with each primer and combining both primers.

Results are represented as percentages of occurrences (POO) (2a: R. ferrumequinum, 2b: R. euryale, 2c: R. hipposideros). “Others” comprise the orders with lesser frequencies: Araneae, hemiptera, hymenoptera, psocoptera and trichoptera. GIL: Gillet; ZEA: Zeale; COMB: Combination of both primer sets.

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

Download:

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

Table 2. Species identified in faeces and their occurrences with each primer set (Zeale’s and Gillet’s), and combining results, arranged by prey orders.

https://doi.org/10.1371/journal.pone.0220081.t002

Download:

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

Table 3. Main orders of prey consumed identified in faeces of the three species of horseshoe bats.

https://doi.org/10.1371/journal.pone.0220081.t003

The interspecific overlap of the diet obtained with Zeale primer sets is not significantly different than the expected by chance (O_jk = 0.18, P = 0.076). R. ferrumequinum and R. hipposideros show the highest overlap (O_jk = 0.34), followed by R. euryale and R. ferrumequinum (O_jk = 0.12), and R. euryale and R. hipposideros (O_jk = 0.08). On the contrary, the overlap based on Gillet primers was significantly higher than expected by chance (O_jk = 0.50, P = 0.002) with the least overlap between R. euryale and R. ferrumequinum (O_jk = 0.25) and the highest overlap between the other two species pairs (R. euryale—R. hipposideros: O_jk = 0.74; R. ferrumequinum—R. hipposideros: O_jk = 0.52). In any case, the effect size is still generally large and the p-value generally small. Similarly, when results of both primer sets are combined the overlap was higher than expected by chance (O_jk = 0.34, P = 0.001). Again, the least overlap was showed by R. euryale and R. ferrumequinum (O_jk = 0.15), while the other two couples show a higher overlap (R. euryale—R. hipposideros: O_jk = 0.38; R. ferrumequinum—R. hipposideros: O_jk = 0.50).

Discussion

Our results confirm the relevance of combining complementary primers to describe the diet of generalist insectivorous bats with amplicon metabarcoding techniques. In general, each pair of primers revealed a subset of the prey composition, with a small fraction of the species being detected by both of them. As a result, the interplay between the primer taxonomic affinity and dietary composition of the bat affected the niche overlap among the three horseshoe bats pictured by each primer.

Due to their more generalist character [31], Gillet primers amplify and identify a higher number of different orders, showing a more diverse diet composition. This generalist character, though, doesn’t cover a full representation or important prey orders such as Lepidoptera. Moreover, their high amplification success comes with the impossibility of identifying a substantial fraction of the amplified DNA (Table 1). On the contrary, the higher selectivity of the Zeale primer set for Lepidoptera and some Diptera [43, 44] might elicit the underestimatimation of other groups of consumed prey, such as Coleoptera, Ephemeroptera, Hymenoptera or Orthoptera [44]. On the positive side, maybe due to their lesser degeneration level, a higher proportion of the OTUs got with Zeale primers were assigned to know taxa (68%, Table 1), providing a deep coverage of lepidopteran prey species.

As both primer sets used in our study amplify regions of the same well-represented COI marker region, the final prey list did not depend of the availability of model species’ sequences in the databases. Moreover, both primers target mini-COI segments of similar size, short enough to be present in faecal samples after digestion [24,30–32]. In fact, the slightly higher amount of primary OTUs yielded by Gillet primers is consistent with the fact that this primer set amplifies moderately shorter fragments than Zeale ones (133 vs 157 bp, respectively). Nevertheless, the more degenerated Gillet primers could have amplified more DNA fragments, generating more OTUs but with a lower assignment to prey taxa whereas Zeale primers produced less OTUs but with a higher assignment to taxa, consistently with their lower degree of degeneration.

The latest molecular study carried out by Galan et al. [31] described the diet of R. ferrumequinum as mainly consisting in Lepidoptera and Diptera, whereas morphological studies have described Coleoptera and Lepidoptera as the most important prey orders [47,62–64]. In our study Diptera occurs in 96% of the samples, closely followed by Lepidoptera (79%) and Coleoptera (71%). Nevertheless, lepidopterans were the main prey order detected with Zeale primers, followed by dipterans, whereas with Gillet ones coleopterans and dipterans prevailed, lepidopterans falling down to a modest third place. In fact, coleopterans were only amplified by Gillet primers and some of the most important lepidopteran families (namely Geometridae, Noctuidae and Totricidae) were disclosed by Zeale. We report the family Geometridae and frequently occurring species such as Rhipidia maculata (Diptera, Limoniidae), Serica brunnea (Scarabaeidae) and Pharmacis fusconebulosa (Lepidoptera, Hepialidae), for the first time among prey of R. ferrumequinum.

R. euryale has been widely considered a moth specialist [62] and, according to Koselj [65] and Dietz [66], lepidopterans make up 90% of its diet. In our study, separate molecular studies performed with the two primers showed a narrow specialization level of R. euryale for lepidopterans, but seasonally complemented by ephemeropterans, hemipterans, hymenopterans and trichopterans [33,42]. In fact, ephemeropterans had been previously reported as prey of R. euryale in North Africa [67]. Noteworthy, lepidopterans are almost the only preyed order if using Zeale, whereas Gillet gives similar importance to lepidopterans and dipterans. Four out of the five most frequently occurring lepidopterans—namely Capsula sparganii (Noctuidae), Udea ferrugalis (Crambidae), Lycopohotia porphyrea (Noctuidae), Scoparia sp. (Crambidae) where solely amplified by Zeale. The two most preyed species Capsula sparganii and Austrolimnophila ochracea have not been described before in the diet of R. euryale.

R. hipposideros is known to prey mostly upon Diptera Nematocera, followed by Lepidoptera and Neuroptera [68, 69]. In our study, the combined use of both primer sets overall confirms the diet composition depicted in previous studies [42, 63, 67–70], even if the family choices within dipterans and neuropterans differ. When only Gillet primers were used, though, prevalence of dipterans (mainly limonids) inflated, while that of lepidopterans and neuropterans (hemerobids) deflated. For Zeale, instead, dipterans and lepidopterans appeared almost in the same frequencies, closely followed by neuropterans. This results agree with Andreas et al. [47] who reported a highly prevalence of Lepidoptera in the pellets. Some of the most important families within Diptera reported by morphological studies [71], namely Tipulidae, Empididae, Muscidae and Culicidae are also represented within the most frequent prey species in the current study, Empididae only amplified by Gillet and Muscidae only by Zeale. Two of the most frequent limonid prey species—namely Neolimonia dumetorum and Limonia nubeculosa—had been previously reported by Galan et al. [42]. Conversely, some other frequent limonids such as Austrolimnophila ochracea, Rhipidia maculata, Dicranomyia modesta or along with other frequent prey species—Hemerobius humulinus (Neuroptera, Hemerobiidae), Wesmaelius nervosus (Neuroptera, Hemerobiidae), or Pseudatemelia josephinae (Lepidoptera)—had not been reported before.

Noteworthy, most of the molecular diet studies carried out on bats exclusively with Zeale primers not surprisingly have concluded that moths or/and Diptera were their main prey: e. g. Barbastella barbastellus, Plecotus macrobullaris, Chalinolobus gouldii, Vespadelus regulus, Nyctophilus gouldi, Eptesicus nilssonii, Myotis brandtii, M. daubentonii, M. mystacinus and Plecotus auritus [17, 24, 35, 39]. In some of the studies [35] previously known prey species–such as Coleoptera, Hymenoptera, Isoptera and Trichoptera–were lacking. Notwithstanding the well-settled importance of moths and dipterans as prey of insectivorous bats (e.g. [63,3]), the reliability of trophic scenarios depicted so far only with Zeale primers is still to be ascertained. Therefore, new studies using combination of primers are highly advisable in order to acquire the fullest dietary view whether to confirm the results obtained with Zeale.

Furthermore, the strong primer bias reported herein cast doubts on the results of previous studies comparing the trophic niche overlap between sibling bat species, carried out exclusively with a single primer set (Zeale). For example, a study comparing the niches of R. euryale–R. mehelyi [34] showed a high degree of diet overlap. Razgour et al. [51] obtained similar results for Plecotus austriacus and P. auritus. Some other studies have also analysed the diets of sympatric bat species [40] based on Zeale primers. Even though the diet overlaps these studies reported cannot be denied, other primers may well unveil additional consumed prey and higher levels of resource partitioning among the species pairs.

Last but not least, previous studies have shown that Gillet primers are useful to identify predators’ DNA [42,46]. In this study we identified almost all the faecal samples for their predator, in except from two R. ferrumequinum samples. Galan et al. [42] argued that a mismatch (T/C) at the 30-end of the reverse primer could be at the origin of their higher rates of amplification failure for some bat species, including R. ferrumequinum. We also identified DNA remains indicating unexpected interactions, including secondary predation events. Thus, we found one R. euryale faecal sample containing Bos taurus sequences, likely traces of bovine animal excrements coming from the common housefly (Musca domestica). Lichtheimia ramosa was identified in R. euryale and R. hipposideros. This is a fungus living in soil and vegetable wastes that infects both animals and humans. These results must be considered with caution, though, because field contamination cannot be fully discarded.

Conclusion

On the one hand, the present study shows that the combination of primer sets with different degeneration degrees that amplify different sub-regions of a specific marker allows identifying a broader and more complete prey spectrum for generalist predators like insectivorous bats. The complementarity of the results yielded by both primer sets lie at the species level, since very few prey species’ sequences were amplified by both primers. For instance, thanks to the use of both primer sets we were able to reveal that R. euryale, though considered a moth specialist, complemented its diet with very diverse prey. On the other hand, our results stress the constraints of the PCR-based metabarcoding diet studies, due to the biases of the many methodological procedures and steps involved in them. Due to biases involving false positives, false negatives and varying affinities to amplify different sequences, we must be extremely cautious when drawing any conclusion from the results gathered, and even when comparing results of different studies. This strong bias when amplifying prey sequences will yield erroneous or at least partial pictures of the trophic requirements of the consumers and resource partitioning among them. In this context, the use of complementary primers improves any assessment of species trophic spectrum and resource partitioning. In this case, we have seen that among the three studied horseshoe bats resource partitioning exists, mostly between R. euryale and R. ferrumequinum.

Supporting information

Acknowledgments

We thank the colleagues who helped us during the fieldwork, and also the Sequencing and Genotyping Unit—Genomic Facilities—SGIker (UPV/EHU/ERDF, EU) for technical and human support provided by them. We are also very thankful to editor and two anonymous reviewers, whose suggestions notably improved our previous manuscript.