Pellets of proof: First glimpse of the dietary composition of adult odonates as revealed by metabarcoding of feces

Abstract Recent advances in molecular techniques allow us to resolve the diet of unstudied taxa. Odonates are potentially important top‐down regulators of many insects. Yet, to date, our knowledge of odonate prey use is based mainly on limited observations of odonates catching or eating their prey. In this study, we examine the potential use of metabarcoding in establishing the diet of three adult odonate species (Lestes sponsa, Enallagma cyathigerum, and Sympetrum danae) at a site in southwestern Finland. To this purpose, we compared three different methods for extracting DNA from fecal samples: the Macherey‐Nagel Nucleospin XS kit, a traditional salt extraction, and the Zymo Research Fecal Microprep kit. From these extracts, we amplified group‐specific mitochondrial markers (COI and 16S rRNA) from altogether 72 odonate individuals, and compared them to comprehensive reference libraries. The three odonate species show major overlap in diet, with no significant differences between individuals of different size and/or gender, reflecting opportunistic foraging of adult odonates. Of a total of 41 different prey species detected, the most frequently consumed ones were Diptera, with additional records of six other orders. Based on our data, the best DNA extraction method is the traditional salt extraction, as it provides the most information on prey content while also being the most economical. To our knowledge, this is the first study to resolve the species‐level diet of adult odonates. Armed with the appropriate methodological caveats, we are ready to examine the ecological role of odonates in both terrestrial and aquatic food webs, and in transferring subsidies between these two realms.

Insects in the order Odonata, including both dragonflies and damselflies, are a globally diverse group of insects with around 5,900 species described to date. Odonates are an important top predators in many aquatic and riparian ecosystems, thus representing both the aquatic and aerial environments and they are fairly well-known due to decades of research (Corbet, 1999). Indeed, odonates have long been model organisms for ecological research, and form a highly promising taxon for future genomics focused research (e.g., Bybee et al., 2016;Córdoba-Aguilar, 2009). Their role in global ecosystems is likely large, as odonates are common, large and remain predators throughout their life cycle (Askew, 2004). The extended larval stage is spent in water, and both diversity and biomass can be very high (e.g., Corbet, 1999;McCauley et al., 2008). Upon hatching, the adult odonate transfers to the terrestrial realm, thus moving biomass and energy from one habitat to another, and contributing to the predation pressure in a new environment.
However, despite that the role of the odonates in both aquatic and terrestrial ecosystems is likely to be large, the prey use of odonates is poorly documented, especially at adult stage of the life cycle. This is due to multiple constraints: for example, their prey species are usually small and thus hard to identify by observing them in mid-air.
Furthermore, visual prey observations of hunting odonates are likely to be biased toward large prey species. After successful hunting, odonates chew their prey thoroughly, which makes morphological identification of prey remnants from feces practically impossible. Fortunately, molecular tools-involving DNA extraction from predator remains, amplification of prey DNA by PCR, sequencing and identification through comparison to reference sequences (Clarke, Czechowski, Soubrier, Stevens, & Cooper, 2014;King, Read, Traugott, & Symondson, 2008;Pompanon et al., 2012;Roslin & Majaneva, 2016)-are not restricted by these aforementioned obstacles, and for the first time we are able to shed light on precise dietary composition of the odonates.
Understanding the prey use of adult odonates is particularly important, as environmental conditions, for example, food shortage during the adult stage can reduce life span and fecundity, thereby reducing lifetime egg production and leading to numerical effects at the egg and larval stages (Stoks & Cordoba-Aguilar, 2012).
The majority of recent theory about predator-prey interactions is based on the assumption that size is a key factor structuring these interactions (Brose, 2010;Schneider, Scheu, & Brose, 2012). In fact, a number of studies have shown that a significant portion of structural information within food webs can be predicted from body size alone (Stouffer, Rezende, & Amaral, 2011;Williams & Martinez, 2000). This is particularly prominent in aquatic systems where predation is largely limited by the size of the predator's gape-such that the larger a consumer is, the larger its gape and the larger its prey (Brose et al., 2006;Morgan, 1989). While large prey may require too much energy to capture, handle and consume, prey that are too small are not worth the energy invested to capture them (Svanback, Quevedo, Olsson, & Eklov, 2015). This should result in a unimodal relationship between predator and prey body size (Brose, 2010;Woodward, Ebenman, Emmerson et al., 2005;Woodward, Speirs, & Hildrew, 2005); in other words, predators of different size are expected to target prey within a different range, the mode of which should be higher with increasing predator size (Williams & Martinez, 2000). Diet generality may also increase with body size, allowing larger predators to exploit a wider range of prey (Gilljam et al., 2011).
To evaluate the potential for molecular techniques to describe the diet of odonates, we target three sympatric odonate species. Drawing on a comprehensive DNA barcode library of potential prey (the Finnish Barcode of Life, FinBOL; www.finbol.org), we use next-generation sequencing techniques (DNA extraction followed by PCR and Illumina MiSeq sequencing) to answer the following questions: (1) How do methodological choices (DNA extraction techniques and choice of markers) affect our perception of prey use? (2) What prey taxa do these adult odonate predators feed on? and (3) do co-occurring odonate species and sexes of varying size differ in their prey use? We predict, firstly, that methodological choices, especially the selection of PCR primers, will affect the results, with more variable gene regions resolving more prey taxa. Secondly, we expect odonates of different species and sex to differ in size, and this size variation to reflect into prey choice, with larger odonate predators consuming larger prey.

| Study species
To evaluate the potential for molecular, DNA-based techniques based on locus-specific amplification of gene regions to describe the diet of odonates, we target three odonate species, the northern bluet Enallagma cyathigerum (Charpentier, 1840) (Coenagrionidae), common spreadwing Lestes sponsa (Hansemann, 1823) (Lestidae), and black darter Sympetrum danae (Sulzer, 1776) (Libellulidae) into this study (images of the species in Fig. 4). The target species were chosen to represent locally common dragonfly and damselfly species that are phylogenetically divergent, and have different life history strategies while sharing the same habitat and overlapping phenology (Corbet, 1999;Dijkstra, 2006;Dijkstra & Kalkman, 2012). Enallagma cyathigerum, a Coenagrionidae species, overwinters as a larva and develops into the adult stage later than most damselfly species in Finland. Contrary to the majority of damselfly species in Finland, E. cyathigerum forages commonly on open areas, also above watersheds. The second focal species, L. sponsa belonging to Lestidae, overwinters at the egg stage in Finland, and develops rather fast through the larval stage during the summer. The adults are among the most common damselfly species flying in July-August. Lestes species hunt commonly near or inside dense vegetation. The third target species, S. danae belonging to Libellulidae, overwinters as an egg, develops quickly in the spring, and then hatches mainly in July. Sympetrum danae is by far the largest and strongest flyer of the focal species, foraging in open areas mainly by chasing its prey in mid-air.

| Study site and sample collection
Odonate samples were collected with aerial sweep nets. To remove variation between different foraging habitats and available prey, all our study samples were collected from one location in South West The age of focal odonate individuals was determined by the stiffness of their wings and the coloration of their bodies (as described for genus Calopteryx in Plaistow & Siva-Jothy, 1996). Only sexually mature individuals were included in this study. After being caught, sample individuals were placed into individual plastic Sarstedt 10-ml tubes with a piece of moist paper towel added to avoid the dehydration of animals. Individuals were kept in the containers for 24 hr to allow complete defecation. All the fecal material (typically some 1-4 individual fragments of irregular size and shape) produced during this time was regarded as one sample. As a proxy for body size, we measured the length of hind wings and calculated average hind wing length for all the individuals. This metric has been previously shown to correlate well with body size (but see Schneider et al., 2012), and it is fairly easy to measure precisely. Thereafter, the feces were collected into Eppendorf tubes and frozen at −20°C until further analysis.

| Procedures for prevention of contamination
To minimize the risk contamination, we tried to adhere to the principles of ancient DNA processing as far as possible in the current laboratory. All the extraction steps were carried out in carefully cleaned laboratory space, using purified pipettes with filter tips. All the PCR's were carried out in a separate room, and no amplified DNA was transferred back to the pre-PCR facilities. Negative controls containing all but template DNA were carried out for each PCR assay.

| DNA extraction using three different methods
Our dataset consisted of total 72 samples: 24 fecal samples for each three study species as equally distributed among females and males.
The fecal material was not pretreated in any specific way prior to extraction, and the amount was so minimal (approximately 1 × 1 mm) that it was not practical to weigh the samples. The total set of samples was divided into three subgroups each consisting of 24 samples with equal representation of females and males of each study species (each subgroup contained four males and four females per species). One group was processed using ZR Fecal DNA MicroPrep (hereafter abbreviated as ZR; product nr D6012, Zymo Research, Irvine, California, U.S.A.), the second group using NucleoSpin ® Tissue XS Kit (abbreviation NS; product nr 740901, Macherey-Nagel, Düren, Germany), and the third group with a traditional salt extraction method (abbreviation SE) (see Appendix S1 for detailed salt extraction protocol applied: Aljanabi & Martinez, 1997;Pilipenko, Salmela, & Vesterinen, 2012).
We did not measure DNA concentrations, but expected them to be rather low due to the small amount of sample. Moreover, as we amplified both predator and prey DNA, the total DNA concentration as such would not be informative in any case. Thus, we used 1 μl of template DNA regardless of potential differences in the DNA concentrations.

| PCR and Illumina library construction
PCRs were prepared using the protocol of Clarke et al. (2014), with slight modifications related to the different indexing scheme as identified in Vesterinen et al. (2016). For this study, we used dual indexing designed for Illumina sequencing platform, and thus, both forward and reverse primers were tagged with different linkers, unique barcodes and Illumina-compatible adapters (Shokralla et al., 2015). All the individual samples were tagged with a unique index combination.
We chose to include the most common mitochondrial markers used for molecular identification of animals: cytochrome oxidase subunit I (hereafter abbreviated as COI) and 16S ribosomal RNA (16S) (Hebert, Cywinska, Ball, & DeWaard, 2003;Yang et al., 2014). The choice of COI region is natural, as most DNA barcoding to date has been carried out using this gene, resulting in millions of reference sequences available in the BOLD database (ref). The 16S region is the gene region second most commonly employed in DNA metabarcoding, and as this region is more conserved than COI, 16S primers usually amplify a larger set of taxa (Clarke et al., 2014). To amplify suitable fragments of approximately same lengths, we applied two primer sets-COI: primers ZBJ-ArtF1c and ZBJ-ArtR2c after Zeale, Butlin, Barker, Lees, and Jones (2011) and 16S: primers Ins16S-1F and Ins16S-1Rshort after (Clarke et al., 2014).
For this study, the PCR setup was further optimized as follows: for a reaction volume of 10 μl, we mixed 3.4 μl distilled water,

| Sequencing output analysis and OTU identification
The sequencing run yielded 637,087 quality-controlled paired-end reads. The reads separated by each original sample were uploaded to CSC servers (IT Center for Science, www.csc.fi) for trimming and further analysis. Trimming and quality control of the sequences were carried out as follows. Paired-end reads were merged and trimmed for quality using USEARCH version 9 (Edgar, 2010). Primers were removed using cutadapt version 1.11 (Martin, 2017). The reads were then collapsed into unique sequences (singletons removed), chimeras were removed, and reads were clustered into OTUs and mapped back to the original trimmed reads to establish the total number of reads in each sample using USEARCH version 9. Zero length OTUs do not practically differ from traditional clustering of OTUs, but the UNOISE algorithm performs better in removing chimeras, PhiX sequences and Illumina artifacts (Edgar & Flyvbjerg, 2015). Finally, our dataset consisted of 11,793 (COI) and 48,985 (16S) reads which were assigned to species. The OTUs were identified to species or-when specieslevel determination could not be achieved-to higher taxa using BLAST (Altschul, Gish, Miller, Myers, & Lipman, 1990) and the Python script package "bold-retriever," version 1.0.0 (Vesterinen et al., 2016).
Nearly all reads could be identified to at least order level and were thus retained for further analyzes. Rest of the reads were discarded: about 1% of COI reads were identified as Bacteria and plants; of 16S reads less than 1% matched human and other mammalian DNA. A detailed description of the bioinformatics applied is available from the authors upon a request.
Data on taxon-specific size (body length of the prey taxa) were then extracted from literature or pictures from the BOLD database.
In other words, prey taxa as per the taxonomic assignment of our sequences were used to estimate the body size range of the prey consumed, allowing the later testing of our explicit predictions regarding relationships between predator and prey size. From these tests, prey identified as "Hemiptera sp." was omitted, as the size range within this compound taxon is too large to be informative.

| Data analysis
To compare the size of odonates (hind wing length), we used ANOVA to model body size (predator hind wing length average) as a function of predator, sex and predator×sex. An equivalent model was fitted to data on prey size (prey body length) and to number of prey taxa detected (count of prey items in each sample).
In terms of prey used, we characterized the frequency of each prey taxon by its presence/absence at the level of individual odonate droppings. This approach was chosen as with PCR-based approaches, the number of reads has been shown to carry little information about the original quantity of template DNA (Deagle & Tollit, 2007;Pompanon et al., 2012). Frequencies were calculated for each odonate species, for males and females and for different extraction methods.
To compare the effect of the gene region amplified, we calculated the number of prey items found in each sample separately for COI and 16S primers. We then used a Kruskal-Wallis Analysis of Variance procedure to compare the performance of each primer set in terms of the frequency distribution of prey items detected (Kruskal & Wallis, 1952). Likewise, we compared the total number of predator and prey reads among individual extraction methods. To further evaluate the performance of each DNA extraction method, we calculated the performance rate (as a percentage) by dividing the number of samples that produced at least some (prey or predator) reads with the number of samples used in the study. For this performance analysis of each extraction method, we used the read numbers remaining after quality control (see above). These performance metrics were calculated separately for each DNA extraction methods, odonate predator species, as well as for males and females within species.
We used ANOVA to model the number of samples that produced taxonomically assignable reads (as explained above) as a function of the different DNA extraction methods. For each extraction method, we also compared the ratio of predator versus prey reads. To visualize the trophic interactions structures resolved by the molecular data, we used package bipartite (Dormann, Fründ, Blüthgen, & Gruber, 2009) implemented in program R (R Core Team 2012). Semi-quantitative webs were constructed for each odonate predator species, using proportional frequencies as explained above. To study the effects of body size, sex and predator species on variation in prey species composition, we conducted a permutational multivariate analysis of variance (PERMANOVA (Anderson, 2001) for presence/absence data, using 999 random permutations to assess statistical significance.
PERMANOVA analysis was carried out using software R with package "vegan" (Oksanen et al., 2017). To visually compare odonate predator diets at the prey family level, we measured the proportion of the most frequent families (Diptera: Chironomidae, Sciaridae) consumed by each predator species.

| How does the choice of DNA extraction method affect the results?
Different extraction methods returned proportionally similar amounts of predator and prey reads (Fig. 1), although a very different amount of absolute target reads. The salt extraction method produced the highest number of reads (32,947), and also the highest number of prey taxa (27 distinct prey taxa). However, although commercial kits did not produce as many reads as did salt extraction, Nucleospin extraction resulted in almost as many prey taxa (26 prey species), as did the salt extraction method. Zymo Research performed the worse in both metrics, allowing the identification of only 11 different prey taxa.
Success rates for each extraction method varied from 50% to 96% (Table 1). In terms of prey item detection rate per sample, we found no significant differences between different DNA extraction methods (Fig. 2a). Overall, COI primers produced a much lower number of predator reads than did 16S primers, but generally, the pattern of prey detection was similar between these two markers, as most of the samples only contained one prey item, while just a few contained more than four distinct prey species (Fig. 2b). The highest success rate was found when using salt extraction and 16S primers (95.8%), and the lowest success rate was observed with the Zymo Research kit and COI primers (50%; Table 1). On average, samples from females and males yielded similar success rates for both markers (COI: 61.1% vs. 66.7%; 16S: 94.4% vs. 97.2%; Table 1). For COI, most OTUs (33/48) and for 16S, less than half (15/37), were successfully identified to the species level or to an unequivocal higher taxonomic level (genus, family or order). In terms of reads (not OTUs), 87% of the trimmed COI reads and 98% of trimmed 16S reads offered a match in a database (BOLD and GenBank for COI; GenBank for 16S) and were retained for subsequent analysis. We found statistical differences between the DNA extraction methods and also between genetic markers in terms of how many prey items was retrieved per approach ( Labeled raw reads, read counts, and OTU data are available in the Dryad Digital Repository: https://doi.org/10.5061/dryad.5n92p.

| What prey species do odonates feed on and is there variation in the diet between odonate species and sexes?
Altogether, we found DNA from 41 different prey taxa, representing 25 different families and seven orders (Appendix S2). Of these, 10 could only be assigned to family or higher taxa. Overall, the three odonate species differed in size, S. danae being significantly larger than the other two species (Fig. 3a; F 2, 69 =227.7, p = <.0001), with no difference among the two sexes ( Fig. 3a; Sex F 1, 70 =1.36, p = .25; Sex × Predator F 2, 69 = 1.02, p = .36). Nonetheless, the size of the prey consumed did not differ in size among either odonate species (Fig. 3b; In terms of exact prey composition, different odonate species shared many prey species (Fig. 4). This dietary similarity was further confirmed by the multivariate analysis: after accounting for the effect of the body size, no significant differences remained between predators or sexes (ADONIS: R 2 = 0.04, p = .08; Table 2). The most common prey taxa were found in dipteran families Chironomidae (midges) and Sciaridae (dark-winged fungus gnats), and when diet was compared at the family level, then the diet of the three predator species was indeed strikingly similar (Fig. 5).

| DISCUSSION
To our knowledge, this is the first study to shed light on the complete species-level diet of adult odonates. With the help of an extensive national barcode library, we were able to identify over forty prey taxa from the fecal samples. Of the 41 distinct prey identified, 28 were assigned to at least the genus level, and the rest (13) to a family (with one taxon only to the order) level.

| Prey taxa consumed
Among the three odonate species studied here (E. cyathigerum the most commonly consumed prey order was Diptera. This group is undoubtedly one of the most abundant prey types available in the wet habitat where the study was conducted. It is also reported to be the most common prey taxon in previous odonate studies based on visual prey identification and sticky traps (e.g., Baird & May, 1997). Within Diptera, the most frequently observed taxa were families Chironomidae and Sciaridae. Indeed, these are some of the most abundant and diverse insect families in Finland, with approximately 700 and 340 species, respectively (Paasivirta, 2012(Paasivirta, , 2014Vilkamaa, 2014). This also concur with previous study conducted with sticky traps in Japan, which found that at least 80% of individual prey items available for dragonfly Mnais pruinosa were small Diptera (Higashi, Nomakuchi, Maeda, & Yasuda, 1979).
We did not find statistical differences in the diet of the three focal odonates. Although the S. danae was significantly larger than other predator species, no significant difference was found in terms of prey size, detected prey items per sample, or prey assemblage.
S. danae diet largely overlapped with that of the other predators, and included only two prey items unique to this species. In  F I G U R E 2 Prey identification success with different extraction methods or genetic markers. (a) Shown is the number of prey items identified per sample using each of the three DNA extraction methods, with (b) an equivalent graph for the two markers used: COI versus 16S. The zero prey items mean that the sample did not produce any sequences after bioinformatic pipeline. There was no significant difference in the frequency of prey detection between methods

(a) (b)
contrast, the other two predators (L. sponsa and E. cyathigerum) had nine and 11 unique prey species, respectively. At the level of prey families, the dietary patterns were highly similar between all three predators. Thus, all three species will likely exhibit an opportunistic hunting behavior, with slight differences in the exact prey species consumed reflecting chance events associated with the prey taxa encountered.

| Similarities to other air-borne insectivores
Interestingly, the prey assortment detected in odonates was similar to that of bats living in similar habitats in the same region (Vesterinen F I G U R E 4 A semi-quantitative food web of the odonate predator species and their prey, combining data from all extraction methods and both markers. The blocks in the upper row represent predators in each web and the blocks in the lower row the prey species. A line connecting a predator with a prey represents a detected predation event, and the thickness of the line represents the proportional frequency of each predation event. The web was drawn using method "cca" which minimizes the cross-links between predators in R package "bipartite" (Dormann et al., 2009 T A B L E 2 Permutational multivariate analysis of variance using Sørensen dissimilarity matrix of presence or absence of prey species in each sample. Terms were added sequentially to the model meaning that the significance of each term is evaluated against the background of terms above it. Predator body size was measured as the average hind wing length of each individual, with factor Predator referring to the species of dragonfly (three levels)   et al., 2013, 2016). This similarity likely reflects joint features in habitat selection and foraging strategies, but also the high general availability of the main prey taxa (Diptera). Of the other diurnal species living in the same areas, the diet of birds is actually less explored by comparable techniques. Yet, results to date suggest that although birds undoubtedly consume large quantities of dipteran insects, they also catch more butterflies and moths (Lepidoptera), most probably at the larval stage (E. Vesterinen unpublished data). Thus, our methods offer a promising tool for assessing ecological similarities and dissimilarities among predator groups for which data have previously been hard to come by, and allow us to finally start mapping out the ecological significance of odonates. While the three focal odonate species differ in size and foraging tactics, there were no differences in the size of the prey they consumed. These results are rather surprising, as we a priori expected larger odonates to hunt larger prey species.

| Source of DNA
One of the issues complicating the interpretation of molecular information derived from food web studies is the source of DNA: whether it originates directly from the prey of focal predators or does it derive from lower steps in the food chain, a phenomenon commonly referred to as secondary predation (Boyer, Cruickshank, & Wratten, 2015;Sheppard et al., 2005). As the odonates are among the top predators of the insect world, they might consume many other predatory species, resulting in secondary predation.
Furthermore, in case of parasite(oid)ism, it is possible that remnants of host species' DNA could end up in these parasites and further to their predators. In this study, nonetheless, practically all of the prey items seemed to be species which are either herbivorous or do not feed as adults. Thus, the risk of false positives in the prey species list seems low in the current study. One exception for this in our results could be Trombidiformes, which include water mite species that commonly parasitize aquatic insects (e.g., Di Sabatino, Martin, Gerecke, & Cicolani, 2002). As our odonate species hunt close to water bodies, it is highly likely that they have consumed other insects that have been parasitized by the members of the order Trombidiformes and that this DNA is consequently represented in our results mainly via secondary predation. Another question is, whether the prey was caught by active hunting or scavenging for example from spider webs, behavior reported for helicopter damselflies by Ingley, Bybee, Tennessen, Whiting, and Branham (2012). However, even in such a case, it is not to be taken as a fault in the results or methods, but instead a challenge to be tackled by other methods, such as complementary direct observations.
Needless to say, contamination is also a real risk in any study dealing with tiny amounts of degraded DNA. Especially, when DNA is amplified, the unwanted DNA originating from contamination would amplify alongside resulting with false positives. In this study, we followed the procedures from our earlier works to cutting the risks to minimum (Vesterinen et al., 2013(Vesterinen et al., , 2016Wirta et al., 2015). The greatest of caution needs to be taken not to introduce any sources of contaminating material to the laboratory while handling the samples.
The amplification has to be done in a separate room (post-PCR), and no amplified DNA should be taken back to the pre-PCR facilities. The inclusion of negative control samples is a standard nowadays, but also positive "mock community" samples are used increasingly (Beng et al., 2016). The idea of mock communities is to add a sample containing a known mixture of potentially expected DNA and using the information from final data to interpret the molecular data quality. This is something to be built on in the future, although there is no easy way of standardizing the mock community approach between different studies. Despite the lack of positive controls with known DNA mixtures, we trust that we have succeeded in preventing the contamination and that what was found was actually eaten.  (allowing the reconstruction of longer sequences) may be amplified from the same samples. For the current cost-effective library construction protocols, several complementary primers can be added without significantly increasing the total costs. PCR-free methods offer another alternative (see Roslin & Majaneva, 2016), but offer additional problems and will be challenging for diet studies dealing with the current tiny amounts and degraded quality of prey DNA (Paula et al., 2015).

CONFLICT OF INTEREST
None declared.

AUTHORS CONTRIBUTION
Kari M Kaunisto involved in original idea, field work, and writing the manuscript. Tomas Roslin involved in writing the manuscript and statistics.
Ilari E Sääksjärvi: involved in writing the manuscript. Eero J Vesterinen: involved in laboratory analysis, statistics, and writing the manuscript.