How much is enough? Effects of technical and biological replication on metabarcoding dietary analysis

Abstract DNA metabarcoding is increasingly used in dietary studies to estimate diversity, composition and frequency of occurrence of prey items. However, few studies have assessed how technical and biological replication affect the accuracy of diet estimates. This study addresses these issues using the European free‐tailed bat Tadarida teniotis, involving high‐throughput sequencing of a small fragment of the COI gene in 15 separate faecal pellets and a 15‐pellet pool per each of 20 bats. We investigated how diet descriptors were affected by variability among (a) individuals, (b) pellets of each individual and (c) PCRs of each pellet. In addition, we investigated the impact of (d) analysing separate pellets vs. pellet pools. We found that diet diversity estimates increased steadily with the number of pellets analysed per individual, with seven pellets required to detect ~80% of prey species. Most variation in diet composition was associated with differences among individual bats, followed by pellets per individual and PCRs per pellet. The accuracy of frequency of occurrence estimates increased with the number of pellets analysed per bat, with the highest error rates recorded for prey consumed infrequently by many individuals. Pools provided poor estimates of diet diversity and frequency of occurrence, which were comparable to analysing a single pellet per individual, and consistently missed the less common prey items. Overall, our results stress that maximizing biological replication is critical in dietary metabarcoding studies and emphasize that analysing several samples per individual rather than pooled samples produce more accurate results.

multitaxa samples, that is metabarcoding, has greatly renewed the interest in dietary studies, particularly due to the high taxonomic resolution offered by this approach (e.g., De Barba et al., 2014;Kartzinel & Pringle, 2015;Lopes et al., 2015). This has been especially relevant to species whose diet is particularly difficult to study, either due to their secretive behaviour (e.g., Shehzad et al., 2012;Soininen et al., 2015) or due to difficulties to identify prey in dietary remains such as stomach contents, regurgitates and scats (e.g., Arrizabalaga-Escudero et al., 2015;Kaunisto, Roslin, Sääksjärvi, & Vesterinen, 2017;Mollot et al., 2014). However, despite its increasingly widespread use, uncertainties and potential biases associated with the quantification of diets based on metabarcoding are still not well understood, requiring a detailed enquiry on how results are affected by different methodological options (Alberdi, Aizpurua, Gilbert, & Bohmann, 2017;Nielsen et al., 2017).
Diet studies aim to answer three main types of question about animal populations: (a) dietary diversity, generally the number of different prey species consumed; (b) dietary composition, that is the identity of the prey species consumed; and (c) the contribution of each prey species to the diet, quantified as the proportion in numbers, biomass or energetic content (e.g., Baker, Buckland, & Sheaves, 2014;Klare, Kamler, & MacDonald, 2011;Whitaker, McCracken, & Siemers, 2009). Surprisingly, there is a significant knowledge gap on the ability of metabarcoding-based studies to provide accurate estimates of dietary descriptors, particularly under field conditions and involving species with diverse diets (Nielsen et al., 2017). Despite this paucity of quantitative studies, researchers often recognize that metabarcoding can be strongly influenced by numerous factors, which should be accounted for in dietary studies. For instance, dietary descriptors can be strongly affected by amplification bias due to unequal primer binding, which leads to systematic over-or underestimation of the importance of some prey types relative to others (Alberdi et al., 2017;Clarke, Soubrier, Weyrich, & Cooper, 2014).
An important aspect often missed in metabarcoding dietary studies is the impact of both technical and biological replication on final results. Technical replication, that is the number of extractions and PCRs carried out on each sampling unit, is important because both extractions and PCRs have a random component, and a given prey item may be missed in some replicates even if it was present in the original sample. These false negatives are expected particularly if an item's DNA is scarce or if there is a negative primer bias Pansu et al., 2015;Willerslev et al., 2014). Biological replication, that is the number of sampling units analysed per species, including for instance the number of individuals or the number of samples per individual, is important because the number of prey species detected tends to increase with the number of samples analysed. Lack of sufficient biological replication can be detected by either rarefaction or asymptotic species richness estimators, which identify sample sizes as being too small to characterize the biodiversity in a sample (Gotelli & Colwell, 2001). Likewise, the precision of frequency of occurrence estimates is low when biological replication is low, and it varies with the prevalence of the prey items, and thus, a poor description of diet may occur at low sample sizes as a mere consequence of binomial sampling (Trites & Joy, 2005). These problems are worse when there is high variation in diet composition among individuals according for instance to gender, age or individual preferences (e.g., Mata et al., 2016;Pagani-Núñez, Valls, & Senar, 2015;Pleguezuelos & Fahd, 2004), and there may also be intraindividual variations due for instance to temporal changes in prey availability (Burgar et al., 2014;Clare, Symondson, & Fenton, 2014b;Clare et al., 2014a).
Here, we address the impacts of technical and biological replication on the results of metabarcoding dietary analysis, focusing on the European free-tailed bat (Tadarida teniotis). This species was considered suitable because previous studies (Mata et al., 2016;Rydell & Arlettaz, 1994) have shown that it is a specialist predator of moths (Lepidoptera) and thus may be less affected by problems of primer bias than species feeding on a wider range of taxonomic groups. Furthermore, moths are well represented in reference barcode databases, which reduces problems due to unidentified MOTUs. Finally, metabarcoding dietary studies have often focused on bats (e.g., Arrizabalaga-Escudero et al., 2015;Hope et al., 2014;Razgour et al., 2011), thus making it possible to evaluate the implications of our results in the context of widely used replication options. In this study, we evaluate how variability among (a) individual bats, (b) faecal pellets of each bat and (c) PCRs of each pellet affect estimates of diet diversity and composition and on the frequency of occurrence of the prey items. Also, we tested the effects of analysing pools of samples vs. separate samples per individual, as these two variants are often used in dietary studies (e.g., pools: Burgar et al., 2014;Clare et al., 2014aClare et al., , 2014bKrauel, Brown, Westbrook, & McCracken, 2018;individuals: Hope et al., 2014;Mata et al., 2016;Vesterinen, Lilley, Laine, & Wahlberg, 2013). Our results were used to analyse the level of replication required to obtain accurate descriptions of predator diets using metabarcoding.

| Study design
This study was based on the dietary metabarcoding analysis of 20 European free-tailed bats (Tadarida teniotis), using both a 15-pellet pool and 15 separate pellets per bat, and three PCR replicates per each pool and pellet ( Figure 1). The number of individuals analysed is within or close to the range used in previous studies investigating for instance trophic structure in bird and bat assemblages (Burgar et al., 2014;Crisol-Martínez, Moreno-Moyano, Wormington, Brown, & Stanley, 2016;Emrich, Clare, Symondson, Koenig, & Fenton, 2014;Razgour et al., 2011;Sedlock, Krüger, & Clare, 2014). The number of pellets analysed separately for each individual is much larger than that of previous studies, which analysed either a single pellet or a pool of pellets per bat. The number of PCRs per sample is within the range (two to four) of recent studies using multiple PCRs (Biffi et al., 2017;Galan et al., 2018), although the large majority of dietary studies has been based on a single PCR per sample (e.g., Burgar et al., 2014;Crisol-Martínez et al., 2016;Emrich et al., 2014;Razgour et al., 2011;Sedlock et al., 2014). Metabarcoding was carried out separately for each combination of bat × pellet (or pool) × PCR, yielding 960 sampling units, for which we recorded the presence/absence of each prey species. To investigate the effects of pellet sample size on the results of dietary studies, we selected randomly one PCR replicate per pellet (320 sampling units) and quantified how increasing the number of pellets analysed affected estimates of both diet diversity and the frequency of occurrence (FO) of the most important prey species. Also, we compared diet diversity and FO estimates for separate pellets and pooled samples. Finally, we used the overall sample to quantify the contribution of variation among individual bats, pellets and PCR replicates to variation in diet composition.

| Bat pellet sampling
European free-tailed bats (Tadarida teniotis) were mist-netted at their roosts in five bridges located in northeast Portugal (N41°09′-42°00′), in April-October 2012 and 2013, under an ongoing monitoring programme (Amorim, Mata, Beja, & Rebelo, 2015). Individual bats were placed in clean cotton bags, from where guano pellets were collected. We recorded gender, age (juveniles vs. adults) and sampling date of each individual. Pellets were stored in tubes containing silica gel and refrigerated at 4°C until DNA extraction. Pellets from a subset of 143 individuals were used in a previous study to describe the diet of European free-tailed bats (Mata et al., 2016), while for this study, we selected the pellets from a different subset of 20 individuals that had left more than 30 guano pellets in the same capture event.

| Molecular analysis
We extracted DNA from each sample using the Stool DNA Isolation where one pellet and a pool were selected per individual and sequenced at "low coverage" (0.1%) and "high coverage" (1.5%). The actual coverages achieved are provided in Supporting Information Table S1.

| Bioinformatics and prey identification
We used OBITools (Boyer et al., 2016) for general sequence processing. Briefly, paired-end reads were aligned and assigned to samples, barcodes and primers were removed, and finally, sequences were collapsed into haplotypes. Singletons were removed, as well as sequences smaller than 155 bp and longer than 159 bp. The remaining haplotypes went through "obiclean," a method that allows the removal of haplotypes differing 1 bp from each other, if one has a higher read count than the other in every sample. From each PCR, we further removed haplotypes representing less than 1% of the total number of reads and those containing stop codons. We then compared the haplotypes retained against known sequences within the BOLD database (www.boldsystems.org) and unpublished sequences of arthropods collected in northern Portugal. Haplotypes that best matched the same species were collapsed into a single taxon unit. For the haplotypes for which only family, order or classlevel identification was possible, a neighbour-joining tree was built with all haplotypes to cluster similar sequences (>98% similarity) into distinct taxa (e.g., Cerambycidae haplotypes with divergences above 98% among them were clustered into Cerambycidae 1, Cerambycidae 2 and so on). Although this approach may artificially increase the number of taxa present in some cases, it was taken to avoid removing from further analysis taxa that are less represented on BOLD and for which genus or species-level identification is often not possible.

| Data analysis
We analysed how pellet sample size affected estimates of diet diversity by building species accumulation curves per individual, as a function of the number of pellets analysed (Colwell & Coddington, 1994).
We used both the actual number of species recorded and the Chao2 estimator of species richness (Chao & Chiu, 2016). We then averaged estimates for each pellet sample size across the 20 bats analysed, to produce a mean species accumulation curve per individual.
Estimates along this curve were compared to richness estimates obtained from the analysis of a pellet pool per individual. To evaluate the effects of sequencing depth, we tested for the difference in species richness in estimates based on one pellet and on a pool of 15 pellets, both at low and at high coverage. We used generalized mixed linear models (GLMM) with logit link and binomial errors, specifying individual bats as the random component, to test whether the probability of detecting a given prey item in pools was related to its frequency of occurrence in the sample of separate pellets (FOpel).
The contribution of biological and technical replication to variation in diet composition was analysed using PerMANOVA (Anderson, 2001). Specifically, we modelled the contribution of three independent components: (a)  FOpel was used to investigate whether error rates tended to be systematically lower (or higher) in prey that was frequently consumed by particular individuals, although not necessarily at the population level. We also used beta regression to estimate whether the error rates of FO estimates in pools varied in relation to FOtot and FOpel.
Simulations were implemented in the R script described in Supplementary Material, while beta regression was carried out using the "BETAREG" package (Cribari-Neto & Zeileis, 2010).

| RESULTS
Metabarcoding of free-tailed bat faecal pellets detected 153 taxa from nine insect orders, of which 65.4% were Lepidoptera (Supporting Information Table S2). Most taxa (77.1%), including 95% of the Lepidoptera, were unambiguously assigned to a single species or to a group of two or three closely related species within the same genus. The seven species with the highest frequencies of occurrence in a single pellet (low: 6.3 ± 3.9; high: 6.2 ± 3.7). The GLMM indicated that the probability of detecting a given prey item in a pool was strongly related to its frequency of occurrence in the diet esti- Information Figure S1).
PerMANOVA showed that variation in species composition among sampling units was significantly affected by variation among individuals, pellets within individuals and PCRs within pellets (Table 1). However, the highest variation in the identity of species  Table S3). The error rates always declined with the number of pellets analysed per individual, but for a given sample size, the error rates tended to be higher for species with high FOtot (i.e., prey items consumed frequently by the population) and that they tended to be lower for species that had higher FOpel (i.e., prey items consumed frequently by particular individuals) (Figure 4) that FOpel was the main factor affecting variation in the error rate of pool FO estimates across prey items ( Figure 5; Supporting Information Table S4).   Although our results are based on a single case study that may be affected by some idiosyncrasies and limitations, this is unlikely to affect the generality of our conclusions to a significant extent. One possibility is that our results were largely driven by the particular species studied, as it consumes a wide range of different prey items (Mata et al., 2016;this study), and thus, it may require higher levels of replication than species with less diverse diets. Although diverse diets may indeed be more difficult to estimate (Nielsen et al., 2017), there are many species such as insectivore bats and birds that feed on a very wide range of taxa and thus may be as prone to insufficient biological replication as European free-tailed bats. Another limitation is that we did not have information on the "true" diet, against which our metabarcoding results could be compared. Previous field studies have circumvented this problem by comparing metabarcoding results with those from visual or stable isotope analysis (Nielsen et al., 2017), but this is not without problems, because all methods have their own errors and biases. Therefore, these comparisons do not show which method is closer to the "truth," but only whether different methods provide consistent results. In these circumstances, we believe that our approach of assessing how estimates of diet descriptors vary with replication levels is warranted, although further research is needed on the extent to which the method provides accurate estimates of what is actually eaten by free-ranging animals. | 171 Finally, our study was based on the analysis of just 20 bats, with all pellets of each bat collected in the same night, and thus, it might be argued that our own study had insufficient biological replication.
Although this sample size is comparable to that of previous studies, we recognize that it may be insufficient to describe in detail the diet of European free-tailed bats. However, it highlights the difficulties of accurately estimating what 20 individuals have eaten during a single night, thereby emphasizing the challenges of inferring diets for entire populations over long time frames. This problem is not restricted to DNA metabarcoding studies of diet, however, except that their increased sensitivity of detection will make real biological variation in diet more detectable. Population diet is an inherently complicated ecological trait to characterize by any methodology, and this has been noted in the past for many dietary studies using different methodologies (Nielsen et al., 2017).
Our results support the view that technical replication affects the estimates of diet descriptors (e.g., Alberdi et al., 2017;Willerslev et al., 2014), although its impact was much lower than that of biological replication. Although there was variation among PCR replicates in the composition of prey items, this was about 13 times lower than variation among pellets of the same individual bat and about 200 times lower than variation among bats.
The low variation among PCR replicates suggests that prey DNA concentration was high and its degradation was low in bat faecal pellets, which are factors known to affect the amount of false positives and negatives and thus technical reproducibility in metabarcoding studies . In contrast to PCR replicates, the magnitude of variation among individuals was particularly striking, suggesting that different individuals fed on different prey items. Reasons for this are unknown, but they may be related to the effects of season, gender, age or foraging habitat. Random factors may also have played a major role, related to haphazard encounters between each foraging bat and a particular set of prey items in the night when pellets were collected. Variation among pellets of the same individual is also noteworthy, with the accuracy of diet diversity and frequency of occurrence estimates increasing markedly with the number of pellets analysed. These results seem surprising, because it might be expected that different pellets collected in the same time from a single individual would be representative of a single meal consumed in that night, thereby leading to low variability in dietary information among pellets. However, bats have an extremely rapid digestion and a high passage rate of food through the digestive tract (Staliński, 1994), and thus, differences in pellet content within indi- is not entirely clear as it seems somewhat counterintuitive, because the DNA from species in individual pellets should also be present in a mix of the same pellets. However, common species in a mix will become proportionally more abundant, and rare species, which appear in low quantities in just a few pellets, will show an even smaller proportion. Therefore, the most likely explanation is that low abundance templates are not detected because of competition during PCR with proportionally more abundant templates. It is also possible that during DNA extraction, pooled samples might saturate the spin column and only the most common species get eluted. Nevertheless, the error in frequency of occurrence estimates is still slightly higher for pools even for common species. This is because pools seem to detect mostly what is highly abundant within individuals, meaning that the analysis of a single pellet is as likely to detect abundant species as is the analysis of a pool. It should be noted, however, that pooling may still be a necessary step when the initial DNA template is too low for extraction and amplification, although results need to be interpreted carefully given the errors associated with sample pooling revealed in our study.
Taken together, our results have important implications for the design of metabarcoding dietary studies, emphasizing the prominent role of biological replication to obtain robust estimates of diet diversity and composition, and the frequency of occurrence of prey items.
In particular, the high variability reported here both among and within individuals points out that large numbers of individuals and sufficiently large numbers of samples per individual need to be analysed if the true diversity of the population's diet is to be recovered.
Determination of sufficient levels of biological replication in general, however, will depend on the particular scientific questions being asked and the dietary characteristics of the species being studied.
For instance, although in conventional studies of bat diets, it is generally agreed that 20-50 samples should be analysed for each ecological group under study (e.g., species, site, season, gender and age; Whitaker et al., 2009), this may or may not be sufficient dependent on the levels of variability within groups and the actual differences in the value of diet descriptors among groups. Larger sample sizes may thus be needed to detect differences in trophic niche between two species showing high intraspecific dietary heterogeneity due to gender, age or seasonal effects than between adult males and females of the same species on a given season, for example. On the other hand, smaller sample sizes may be more acceptable in studies aiming to provide broad descriptions of dietary patterns in diverse species communities, particularly when these include species that are rare or otherwise difficult to study than when testing specific hypothesis in community ecology requiring precise dietary estimates.
Therefore, scoping studies may need to be conducted before embarking in full-scale projects, using power analysis to estimate the levels of biological replication required to detect a given effect size at a predefined probability level (Ferry & Cailliet, 1996). When this is impractical, researchers may need to take a precautionary approach and try to maximize the number of samples analysed, which is increasingly feasible due to the ever lower costs of high-throughput DNA sequencing.

ACKNOWLEDGMENTS
The study was funded by Fundação para Ciência e Tecnologia (