Table for five, please: Dietary partitioning in boreal bats

Abstract Differences in diet can explain resource partitioning in apparently similar, sympatric species. Here, we analyzed 1,252 fecal droppings from five species (Eptesicus nilssonii, Myotis brandtii, M. daubentonii, M. mystacinus, and Plecotus auritus) to reveal their dietary niches using fecal DNA metabarcoding. We identified nearly 550 prey species in 13 arthropod orders. Two main orders (Diptera and Lepidoptera) formed the majority of the diet for all species, constituting roughly 80%–90% of the diet. All five species had different dietary assemblages. We also found significant differences in the size of prey species between the bat species. Our results on diet composition remain mostly unchanged when using either read counts as a proxy for quantitative diet or presence–absence data, indicating a strong biological pattern. We conclude that although bats share major components in their ecology (nocturnal life style, insectivory, and echolocation), species differ in feeding behavior, suggesting bats may have distinctive evolutionary strategies. Diet analysis helps illuminate life history traits of various species, adding to sparse ecological knowledge, which can be utilized in conservation planning.


| INTRODUC TI ON
Coexistence of sympatric species is facilitated by differences in the use of resources, that is, resource partitioning (Schoener, 1974). Resource partitioning occurs in several dimensions, with regard to resources.
Ultimately, the sum of these dimensions constitutes the ecological niche of an organism, that is, the set of biotic and abiotic conditions in which a species can persist (Holt, 2009). This includes both the distribution of a species and its interactions with other species, but also factors relevant to the fine-scale distribution of species (e.g., microhabitats), their biotic interactions as well as their diet (Wiens et al., 2010).
With a notable adaptive radiation in their evolutionary history, and over 1,300 known species worldwide (Fenton & Simmons, 2015), bats have an important role in supporting global ecosystems through their dietary preferences. This is evidenced primarily through the consumption of nocturnal insects and dispersal of nutrients, pollen, and seeds (Patterson, Willig, & Stevens, 2003). Research on the feeding behavior of species is essential to understanding ecosystem function and the impacts of pollution, habitat destruction, and global climate change (Boyles & Storm, 2007;Kunz, Braun de Torrez, Bauer, Lobova, & Fleming, 2011;Vesterinen, 2015;Vesterinen et al., 2016). Furthermore, establishing factors influencing the extinction risk of bats is essential for their conservation, because they help identify endangered species and provide the basis for conservation (Safi & Kerth, 2004). However, these factors may be difficult to discern between species of bats, of which many appear to share portions of their ecological niches, such as habitat and apparently diet. imals (Roslin, Majaneva, & Clare, 2016;Vesterinen et al., 2016;Vesterinen, Lilley, Laine, & Wahlberg, 2013). Vesper bats within the same feeding guild appear to share a great proportion of their diet (Roswag, Becker, & Encarnação, 2018). Because insectivorous bats opportunistically consume prey that may be periodically abundant (Vesterinen et al., 2013), this leads to significant temporal changes in the diet (Vesterinen et al., 2016), but could additionally result in a large overlap in dietary niches, suggesting resource partitioning occurs in other ecological dimensions.
Here, we unravel the resource partitioning of five resident vesper bats in southwestern Finland through deep dietary analysis, including prey species identification, an estimate for prey body size and temporal changes in diet using fecal DNA barcoding. At high northern latitudes, the distribution of bats is constrained by extreme environmental demands and prey availability is more seasonal than elsewhere in their range Shively, Barboza, Doak, & Jung, 2017). The ranges of these five species (Eptesicus nilssonii [Keyserling & Bläsius, 1839], Myotis daubentonii [Kuhl, 1817], M. mystacinus [Kuhl, 1817], M. brandtii [Eversmann, 1845], and Plecotus auritus [Linnaeus, 1758]) show considerable overlap, suggesting that trophic resource partitioning is important in supporting the species in Fennoscandia. We expect to see clear guild-specific segregation in diet between the three different feeding guilds presented by our species, trawling (M. daubentonii), gleaning (P. auritus), and aerial hawking (Figure 1;M. brandtii,M. mystacinus,and E. nilssonii), and that we will see at least a partial dietary overlap among the members of the aerial hawkers. Because of the opportunistic foraging behavior of insectivorous bats (Vesterinen et al., 2013), we also predict significant temporal changes in diet throughout the sampling season (but see Vesterinen et al., 2016). Finally, we predict a positive correlation between predator and prey size, which could be due to the negative correlation between bat size and echolocation frequency, hindering the ability to detect small prey items (Brigham, 1991). To the best of our knowledge, of the species studied here, molecular data on diet exist only for M. daubentonii (Galan et al., 2018;Krüger, Clare, Greif, et al., 2014;Krüger, Clare, Symondson, Keišs, & Pētersons, 2014;Vesterinen et al., 2013Vesterinen et al., , 2016 ), although the dietary contents of all species have previously been described through morphological analysis of fecal remains (Rydell, 1986;Vaughan, 1997).
F I G U R E 1 One of the study species, Myotis brandtii, foraging in its natural environment near the study area in southwestern Finland. M. brandtii catches its prey mainly in flight in an open or semi-open environment. The current study is the first ever published molecular analysis of its diet: Geometrid and tortricid moths constituted half of its diet, while mosquitos, midges, and flies formed another large part of the menu, approximately one-third. Photograph credits: Mr. Risto Lindstedt 2 | MATERIAL S AND ME THODS

| Study species
Of the 13 species of bats occurring in Finland, the species sampled here represent the most common and accessible (Myotis daubentonii, Eptesicus nilssonii, M. brandtii, M. mystacinus, and Plecotus auritus).

| Field sampling
Fecal pellets were collected between April and July 2014 (Table 1) from day roosts of five species of bats in southwestern Finland, and all these roosts were in buildings within approximately 60 km of each other ( Figure 2f). The pellets were collected by placing a clean paper sheet under the roosting bats the day before the collection, and collecting the droppings the next day. The collection was repeated for two or three consecutive days within a period of two weeks. Pellets were stored in RNA later at −20°C until laboratory analysis.

| Laboratory work
We aimed to pool 25 droppings (from the same roost and same time point) into each sample to maximize the number of droppings without the need to analyze hundreds of fecal pellets individually. Only four samples included less than 25 droppings, and for these, we pooled every available pellet for the given time point per roost. We focused sampling on roosts inhabited by a single species, and likewise, we intended to pool pellets from a single species into a single pooled sample. In total, we initially sampled 1,252 fecal pellets from the five bat species in this study ( we wanted to allow comparison of our results with those of other studies using the same primers Kaunisto, Roslin, Sääksjärvi, & Vesterinen, 2017;Koskinen et al., 2018;Krüger, Clare, Greif, et al., 2014;Vesterinen et al., 2013Vesterinen et al., , 2016Wirta et al., 2015;Eitzinger et al., 2018). The PCR and library construction closely followed Kaunisto et al. (2017), TA B L E 1 Information on the sampling details and characteristics of the field and molecular data. Time/roost sampling points per bat species denote how many times per roost the species was sampled: M. daubentonii was sampled from only a single roost (NAU; see Figure 2 for locations of the roost sites in the current study), E. nilssonii was sampled separately from two roosting sites (SJÄ, SSA), M. mystacinus and P. auritus were sampled from the same roost (LAI), and M. brandtii was sampled at two locations (RUI), one of which was shared by P. auritus (ROT). We found no statistical differences between samples from different bat species in the total reads, total prey species richness, or the average number of prey in each pellet     control of the sequences were conducted according to Kaunisto et al. (2017). Consequently, paired-end reads were merged (SFF: ~90% reads successfully merged; ZBJ: ~85%) and trimmed for quality using program USEARCH with "fastq_maxee_rate" algorithm with threshold 1 (Edgar, 2010). Primers were removed using python program cutadapt (SFF: ~99% reads passed; ZBJ: ~96%) (Martin, 2011). We then dereplicated reads using USEARCH "fastx_uniques" algorithm with option "minuniquesize 2", and then, we applied USEARCH UNOISE3 algorithm to cluster these unique reads into ZOTUs (zeroradius operational taxonomical units; Edgar, 2016). In short, UNOISE algorithm allows the simultaneous a) detection and removal of chimeras (PCR artifacts where two fragments of different origin bind together), point errors (substitutions due to incorrect base calls and gaps due to omitted or spurious base calls), and b) results in ZOTUs (zero-radius OTUs) that are superior to conventional 97% OTUs for most purposes, because they provide the maximum possible biological resolution given the data available (Edgar, 2016). Finally, reads

| Bioinformatics and prey list construction
were mapped back to the original trimmed reads to establish the total number of reads in each sample using USEARCH "otutab" algorithm.
After processing, our datasets for this study consisted of 5,449,755 prey reads (produced with primers ZBJ-ArtF1c and ZBJ-ArtR2c) and 1,452,602 bat reads (produced with primers SFF-145f and SFF-351r).
The remaining reads (roughly 30% of total output of the sequencing run; ZBJ: 2,618,342 + SFF: 721,684) were used in another study.
We used the following strict criteria for including prey species in the data: (a) Sequence similarity with the reference sequence had to be at least 98% for the ZOTU to be given any (even higher taxa) assignation, and (b) at least ten reads of the final assigned prey species were required to be present in the final data. We assigned the ZOTUs to species as accurately as possible, utilizing a large reference sequence collection orchestrated by the Finnish Barcode of Life campaign (FinBOL: www.finbol.org) and BOLD database (Ratnasingham & Hebert, 2007), and confirmed that all the prey species were actually recorded from (southern) Finland. After the above trimming, we were able to identify and retain 93% of all the prey reads. To account for the even distribution of reads into separate samples, we used ANOVA to test samples from different bat species for differences in the total reads per sample, total prey species richness per sample, and the average number of prey in each pellet (prey richness divided by the number of pooled pellets). The reads originating from bats in the second dataset were used to confirm the bat species identity. unclear identification or no data on size, we ended up with 1,553 distinct individuals from the bat banding data.

| Data analysis
Traditionally, the read count (or read abundance) data produced in metabarcoding studies are directly transformed into presence/absence data, considered to be more cautious and less biased than using read counts. However, the latest opinion on the field seems to suggest that using normalized read abundance data could be even less biased than mere converting to p/a data (Deagle et al., 2018; see also Vesterinen, 2015;Vesterinen et al., 2016). For this reason, we chose to use relative read abundance (RRA: calculated as the proportion of reads per each prey item in each sample). To make the comparison to earlier studies possible, we also prepared the secondary set of analysis using p/a data or more precisely the modified frequency of occurrence (MFO) data throughout the analysis. MFO was calculated as the proportion of occurrences of each prey taxa in each sample scaled to 100% across all prey items (see Deagle et al. (2018) for the terminology and further discussion on the topic).
To begin our data analysis, we calculated prey species accumulation curves to account for sampling adequacy (Colwell & Coddington, 1994). We used R package "iNEXT" to resample the prey reads and frequencies for each bat species and plotted these against accumulated prey species richness (Hsieh, Ma, & Chao, 2016;R Core Team, 2013).
In order to unfold the trophic interactions resolved by the DNA analysis, we used package bipartite (Dormann, Gruber, & Fründ, 2008) implemented in program R to draw interaction webs for each bat predator species using both RRA and MFO data. For those two cases, where two different bat species were observed in the same roost, we constructed additional webs to analyze the diet between separate samples in each location using RRA data. To further estimate patterns among the dietary assemblages of the five species, we used principal coordinates analysis (PCoA) based on Bray-Curtis dissimilarity (Jaccard similarity for presence/absence data) between samples (Davis, 2002;Podani & Miklós, 2002).
Then, to study the effects of predator species and temporal variation (as week number) on variation in prey species composition in each sample, we conducted a permutational multivariate analysis of variance (with Bray-Curtis for RRA and Jaccard for presence/ absence data), using 9,999 random permutations to evaluate statistical significance (Anderson, 2001)(PERMANOVA; Anderson, 2001).
Analysis of variance was carried out using "adonis" in software R with package "vegan" (Oksanen et al., 2013). Variation was further dissembled using pairwise analysis of variance with package "pairwise.adonis" between all bat species using Bonferroni correction for Kruskal-Wallis analysis of variance (nonparametric ANOVA) procedure to compare body size (FA length) as a function of predator size using command "kruskal.test" in R (Kruskal & Wallis, 1952). To further study the difference between bat species pairs, we applied the Tukey and Kramer (Nemenyi) test with Tukey-Dist approximation for independent samples with R package "PMCMR" (Pohlert, 2014;Sach, 1997, pp. 395-397, 662-664). The same tests were applied to test prey size (wingspan or thorax length as explained above) differences between the bat species.

| General aspects of the diet and the study
Altogether, we identified 547 distinct prey species in 13 arthropod orders (  Figure 3a), although when using presence/absence data, the curves did not seem to reach the plateau yet ( Figure 3b). Nevertheless, we kept M. mystacinus in all the analysis, but interpret the results with relevant caution.

| Dietary patterns of the studied bats
The quantitative prey assemblages (RRA) seem to be very different for all the bat species, as revealed by the bipartite analysis ( Figure 4a). However, when using frequencies (MFO), these patterns are not that clear (Figure 4b). In the current study, different bat species were mainly sampled in different roosts, but luckily prey use does not seem to be vastly related to the roost site, as can be seen from the bipartite analysis from the two sites where two different bat species were sampled from the same roost (Figure 5a,b). The prey use patterns were further illustrated in the PCoA ordinations: Both RRA and presence/absence data ordinations grouped the bat species according to their respective feeding guilds based on differences in the prey species assemblages (Figure 6a,b). In the RRA plotting, first coordinate explained 10.5% and the second coordinate 7.5% of the variation in the data (Figure 6a), and in the plot using presence/absence data, the first and the second coordinates explained 15% and 9.9% of the variation (Figure 6b) both data types a large part of the variation remained unexplained.
Altogether 44 common prey species were shared by all the bat species, and 90 more equally common prey species were shared by four bat species (

| Dietary patterns in the feeding guilds
The feeding guilds are also easily separated by looking the diet at the prey family level (here using percentages from relative read abundance data, but approximately the same ratios can be drawn from the presence-absence data; Table 2): The trawling species (M. daubentonii) predominantly consumes a single prey family, Chironomidae (45.8% of all the reads), which is a highly abundant and species-rich family in southwestern Finland (Lilley, Ruokolainen, Vesterinen, Paasivirta, & Norrdahl, 2012;Paasivirta, 2012Paasivirta, , 2014, but constrained to the vicinity of aquatic environment, where the bat collects its prey from the water surface (Nilsson, 1997). The gleaner (P. auritus) relies on the plentiful moth family Noctuidae (57.2%), which is either caught in flight or from surfaces on vegetation, as some of the prey species are mainly diurnal (Silvonen et al., 2014). The other largely consumed prey family for P. auritus was the coleopteran family Carabidae (18.7%), which is most probably foraged from the ground. The third guild, hawkers, consists of three bat species (E. nilssonii, M. brandtii, and M. mystacinus), which all have distinct prey family spectrum.

| Temporal aspects and predator-prey size analysis
The strong assorting patterns of different bat species seen in plotwebs and PCoA were confirmed when comparing all bat species' diet's together in the analysis of variance (Table 3: Predator: RRA data, df = 4, R 2 = 0.12, p = 0.0001; PA data, df = 4, R 2 = 0.05, p = 0.0033). Despite the limited temporal span of the sampling for each bat (   When the prey assemblages were analyzed separately in pairwise PERMANOVA between species, the diet was significantly different in all compared pairs, except those with M. mystacinus, which was present in the sample with only one sample (    (Table 6). On average, P. auritus consumed the largest prey (Figure 7b,c; Table 6), whereas M. brandtii consumed the smallest prey (Figure 7b,c; Table 6).

| D ISCUSS I ON
Co-occurring species with a relatively short active season offer an excellent setup for the study of dietary strategies. Here, we identified 547 prey species in the diet of five common and abundant boreal vespertilionid bat species. All species fed mainly on two insect orders (Diptera or Lepidoptera), which undoubtedly are among the most available dietary groups (with Coleoptera) in terms of species richness (Erwin, 1982;Stork, 2018) and probably for biomass, although reliable biomass estimates are lacking. The three feeding guilds (trawlers, hawkers, and gleaners) are clearly separated by diet in the data. Moreover, the dietary composition between all bat species differed significantly, a pattern that persisted throughout the results. This pattern was strong enough to be observed in all the interpretations of the molecular data (presence/absence, frequencies, and read count data analysis). The sampling week did not explain the diet for any bat species, but we found differences in average prey size consumed by the bat species, and a positive correlation between bat species size and size of prey, although with a fine marginal.
In concordance with dietary studies on insectivorous bats, we also revealed a high frequency of lepidopteran and Dipteran species in the diets of the sampled species Vesterinen et al., 2016). In fact, combined, these two orders constitute the majority of all predation records in the whole study, regardless of the data type (read counts, frequency, or presence/absence). Especially, P. auritus appears to utilize lepidopteran prey species to a higher degree compared to the other species, although rather surprisingly, ~20% of the diet (in terms of relative read abundance) of P. auritus appears to consist of Coleoptera, particularly ground beetles. All other invertebrate orders are less relied on, although Trichoptera and Neuroptera constitute a small part of the diet in some species. This is expected, seeing as these orders include mass-emerging species, such as Oecetis ochraea (Trichoptera, Leptoceridae), or species which are active and available as prey throughout the season, such as Brachyderes incanus (Coleoptera, Curculionidae), or otherwise very common and abundant species, such as Chrysoperla carnea (Neuroptera), are all found in this study (Vesterinen et al., 2013(Vesterinen et al., , 2016. This primarily highlights the huge biomass and species diversity found in Lepidoptera and Diptera, but secondly, also further establishes the importance of these orders to bat species diversity. Because of the huge biomass of insects worldwide, there are numerous predators in addition to bats, such as fish, birds and even predatory insects, consuming these, and other arthropods as their primary food source (fish: Jakubavičiūtė, Bergström, Eklöf, Haenel, & Bourlat, 2017;dragonflies: Kaunisto et al., 2017;birds, spiders: Wirta et al., 2015). Surprisingly, the prey order-level similarity between different predator taxa is surprisingly high when comparing our results to the aforementioned studies, especially between bats and other flying insectivores.  Vesterinen et al., 2013). This suggests that the diets of our study species could be determined by the abundance and availability of insect prey instead of any particular predator-specific characteristic. In fact, it has previously been reported that bat diet responds to local insect population fluctuations Clare, Barber, Sweeney, Hebert, & Fenton, 2011;Sedlock, Krüger, & Clare, 2014;Vesterinen et al., 2016). Razgour et al. (2011) reported temporal shifts in the proportional frequencies of Lepidoptera and Diptera prey of P. auritus.
We found no evidence of shift in these frequencies in our P. auritus samples. At the latitude where our study was conducted, there is only a two-to four-week difference between the highest abundance peaks for Diptera and Lepidoptera, and furthermore, it may be that even during the period of low abundance, there are still more than enough prey items available for bats (Vesterinen et al., 2016).
Diet comparisons between sympatric bat species using molecular methods are still relatively scarce, but often show considerable overlap in diet, even at the lower taxon level (Krüger, Clare, Greif, et al., 2014;Salinas-Ramos et al., 2015;Ware, 2016). Most studies focus on either closely related species, or species which share a feeding guild, such as the two trawling bats (M. daubentonii and M. dasycneme) in a study by Krüger, Clare, Greif, et al., 2014; Plecotus auritus, the species considered a gleaner and moth specialist, showed a marked difference in PCoA ordination compared to the other two groups. We also discovered a significant difference in the size of prey consumed, with the larger P. auritus consuming larger prey species, whereas the smaller bat, M. brandtii, consumed smaller prey items. This is not surprising as it is generally accepted that the echolocation used by aerial insectivorous bats renders smaller prey items unavailable to larger bats (Brigham, 1991;Waters, Rydell, & Jones, 1995). Additionally, P. auritus, among other members of the genus, possesses a suite of morphological characters (low wingloading, large pinna, low-frequency hearing), which allow them to use both acoustic gleaning and aerial-hawking foraging strategies to capture prey (Coles, Guppy, Anderson, & Schlegel, 1989;Norberg & Rayner, 1987). It is possible that some noctuid prey individuals have been foraged as larvae, as the flight peak of most noctuid prey in the current study is later than the sampling period.(Finnish Biodiversity Information Facility/FinBIF. https://tun.fi/HBF.31668; accessed 2018-08-26). These strategies permit the genus to occupy a specialized feeding niche within European bat assemblages (Roswag et al., 2018 (Dietz, Nill, & Helversen, 2009). This is resource partitioning that could be further dissected by looking at isotopic niches, for instance, to give a complementary scenery to dietary ecology besides DNA-based analysis (Schmidt, Mosbacher, Vesterinen, Roslin, & Michelsen, 2018). Another option would be to increase sampling effort to obtain an even more robust overview of the main prey items. Information on the identified major dietary taxa could then be used to deduct the main foraging habitat, as presented by Alberdi, Garin, Aizpurua, and Aihartza (2012).
The molecular work carried out in this analysis not only highlights the deep insight offered by metabarcoding, but also underlines the dynamic and complementary nature of DNA-based analysis. Based on our earlier field work, we had chosen species-specific roosting sites for the diet analysis of five bat species, to obtain an equal sampling effort. However, when confirming the fecal "donor" by the means of metabarcoding, we noticed some discrepancies between the field data and confirmed data, that is, our M. mystacinus roost was confirmed as an E. nilssonii roost. In future, the molecular confirmation of noninvasively collected samples should be a standard approach, either by traditional Sanger sequencing or cost-effective next-generation sequencing (NGS), depending on the number of samples and the predator and prey species. Also, the importance of a comprehensive reference library (Mutanen et al., 2012;Pentinsaari, Hebert, & Mutanen, 2014;Pilipenko, Salmela, & Vesterinen, 2012), which allows the correct and reliable identification of most prey items, needs to be pointed out once more. This offers the possibility of deeper ecological dietary studies, such as prey size analysis (Pentinsaari et al., 2014). While some prey items had not been described with a scientific species-level name in this study, a reliable estimate of their size could be inferred using the so-called barcode index numbers (BIN; Ratnasingham & Hebert, 2013) to trace the images for measurements. This emphasizes the significance of public and easy-accessible reference library systems, such as BOLD (Ratnasingham & Hebert, 2007). Although some studies still rely on OTUs (operational taxonomical units) instead of biological species, we highlight the importance of actual prey species determination, which allows a deeper and more robust insight into dietary ecology.
The main drawbacks of the molecular methods are the highly challenging interpretations of the quantitative aspects of the diet, that is, are the most frequently consumed prey items also the most important in terms of biomass and energy gain? While the current practice in many molecular ecological dietary studies using metabarcoding appears to mostly rely on frequency of occurrence (but see Vesterinen et al., 2016), the read counts may actually hold some important quantitative information (Deagle et al., 2018). Here, we tested our data using both frequency of occurrence and read count data and found no major differences in the outcome of the analysis, or more importantly, in the interpretation of the results. This suggests our data have strong ecological message that holds despite the methodological approach used.
Our study supports the existence of dietary flexibility in generalist bats and dietary niche overlapping, especially in bats of the same feeding guild in a highly seasonal ecosystem (Roswag et al., 2018). In fact, it could be the flexibility in feeding strategies which allows species to sustain populations in arctic and subarctic regions . Additionally, a great proportion of niche differentiation most likely also occurs outside the diet dimension where an almost infinite number of possible axes exist for competing species in the n-dimensional niche hyper-volume (Hutchinson, 1957). Even minor differences in a number of different axes can result in a substantial overall difference (Privitera et al., 2008). Clearly, the next logical step is to utilize deep dietary analysis, alongside other ecological (LIDAR: light detection and ranging method, etc.) and behavioral (GPS-tracking) datasets to begin to understand niche realization and resource partitioning in species to a far higher accuracy than has been available to date.

AUTH O R CO NTR I B UTI O N S
EJV and TML designed the study, collected the data, and wrote the first version of manuscript. ASB collected samples in the field and gathered prey species measurements and the map data. AIEP and EJV conducted the molecular work and data analysis. All authors contributed to the final version of the manuscript.

DATA ACCE SS I B I LIT Y
Labeled raw reads and OTUs are available in the Dryad Digital