COI Metabarcoding Provides Insights into the Highly Diverse Diet of a Generalist Salamander, Salamandra salamandra (Caudata: Salamandridae)

DNA metabarcoding has proven to be an accessible, cost-effective, and non-invasive tool for dietary analysis of predators in situ. Although DNA metabarcoding provides numerous benefits in characterizing diet—such as detecting prey animals that are difficult to visually identify—this method has seen limited application in amphibian species. Here, we used DNA metabarcoding to characterize the diet of fire salamanders (Salamandra salamandra) (Linnaeus, 1758) in three distinct regions across the northwestern Iberian Peninsula. To test the efficiency of COI-based metabarcoding in determining salamanders’ diet diversity, we compared our COI-based results with results from traditional diet studies from neighboring and distant populations, as well as with recent findings obtained in a DNA metabarcoding study using 18S. Two COI primers were used in combination to investigate the potential impact of primer bias in prey detection. Our COI metabarcoding approach increased taxonomic resolution and supported a generalist diet in S. salamandra. Between primers, there were no significant differences in the diversity and richness of prey detected. We observed differences in the prevalence of prey identified between sampling regions both in our study and in other studies of S. salamandra diet. This COI metabarcoding study provides recommendations and resources for subsequent research using DNA metabarcoding to study amphibian diets.


Introduction
Diet studies are fundamental to understanding species' dietary habits [1], food webs [2], and trophic niches [3,4], which are key traits in many ecological processes and for the conservation and management of species and ecosystems. DNA metabarcoding as a means of dietary analysis has been used in many taxonomic groups [5,6], but has been underutilized in amphibians [7]-particularly in salamanders. Salamanders serve an important role as mesofaunal predators [8], often comprising a large portion of ecosystem biomass [9], and have low energy requirements, making them potential energy sinks in ecosystems [10]. Moreover, salamanders may also exert a top-down effect on invertebrate community composition and nutrient cycling [7,11,12], which makes studying the diets of salamander species especially relevant for understanding their role in ecosystem functioning.
Visual inspection of stomach or fecal contents is a useful but inconsistent means of diet characterization in salamanders [13,14]. Stomach contents provide insight into recent consumption [3,15,16], and while stomach-flushing avoids sample mortality, it is an invasive approach to diet analysis [14,[16][17][18]. Fecal content inspection is less invasive, but introduces a bias favoring hard-bodied prey species that are not fully digest [13].
DNA metabarcoding can help to identify prey species that are consumed over a longer period of time, avoids the detection bias against soft-bodied prey, and requires less taxonomic training in prey identification [19]. Similar to visual inspections, DNA metabarcoding can also indicate the preferential diets of salamanders and their role as predators in communities [7,17,20,21], and can inform us about invertebrate biodiversity. To our knowledge, only one study has applied DNA metabarcoding to investigate the diet of adult salamander species. Specifically, Wang et al. [22] used the 18S ribosomal RNA (18S) region to characterize the diets of adult fire salamanders (Salamandra salamandra) (Linnaeus, 1758) by collecting fecal samples from three Belgian forests. While 18S has proven useful to detect potential prey items, the use of more informative (i.e., variable) DNA fragments such as the cytochrome c oxidase I (COI) benefits from a large reference database that is supported by the Barcode of Life Data System (BOLD) [23,24], and the relatively high variability of the region allows for high-resolution taxonomic assignment [5,[25][26][27].
This study aims to provide an update on the diet of S. salamandra while evaluating the use of COI metabarcoding as an efficient, non-invasive method for diet studies, and comparing it to previous works [15], as well as to gain better insights into the diets and functional roles of salamanders as generalist terrestrial predators [28]. Fecal samples were collected from salamanders across the northwestern Iberian Peninsula. To determine whether a significant difference in prey detection could be attributable to primer bias [29], we compared the performance of two COI primers. Technical considerations include the evaluation of (1) sampling effectiveness to capture prey species, and (2) the usefulness of COI primers as barcodes [7,19,[30][31][32][33].

Materials and Methods
Fecal samples were collected from three distinct regions across the northwestern Iberian Peninsula, including the extended metropolitan area of Porto, a forested area across the Morrazo Peninsula, and the island of Ons, which is mostly dominated by bushes ( Figure 1). Nocturnal sampling was conducted in Morrazo and Porto in the spring of 2021, and in Ons during November of 2020, coinciding with the highest annual activity peaks for the species and under suitable climatic conditions (i.e., rainy nights and temperatures of 10-20 • C). Up to 20 individuals from each site were collected and placed in sterilized individual or group containers. All individuals were returned to their original sampling sites.
To mitigate primer bias, two COI-specific primer pairs-fwh1 (fwhF1 5 -YTC HAC WAA YCA YAA RGA YAT YGG-3, fwhR1 5 -ART CAR TTW CCR AAH CCH CC-3 ) [34] and LerayXT (jgHCO2198 5 -TAI ACY TCI GGR TGI CCR AAR AAY CA-3 , 185 bp;  mlCOIintF-XT 5 -GGW ACW RGW TGR ACW ITI TAY CCY CC-3 ; 313 bp) [27,35]-were used. Salamandra salamandra COI sequences from NCBI (KX094979 & GQ380404) were used to design blocking primers for fwh1 (Ssal_fwhF1-blk 5 -CAA AGA CAT TGG CAC CCT CTA CCT AAT TTT TGG [SpC3]-3 ) and LerayXT (Ssal_mlCO1intF-blk 5 -GAA CAG TCT ACC CCC CCC TTG CCG GAA ATC TGG [SpC3]-3 ). Initial PCR mixes comprised 10 µL of Qiagen Multiplex PCR Master Mix, 0.3 µL or 0.4 µL of 10 mM target primer (fwh1 and LerayXT, respectively), 8.0 µL of 10 mM blocking primer, 2 µL of DNA template, and enough water for a final volume of 25 µL [36]. Thermocycling conditions included an initial denaturing step at 95 • C for 15 min, followed by 40 cycles of denaturing at 95 • C for 30 s, annealing at 45 • C for 60 s, extension at 72 • C for 30 s, and a final extension at 72 • C for 10 min. Amplification success was verified by running 2 µL of PCR product on a 2% agarose gel. Successfully amplified PCR product was diluted at 1:3 of the initial concentration, in order to reduce the primer dimer during indexing. Illumina indices were then annealed to the PCR product with a PCR composed of 7 µL of KAPA Taq ReadyMix, 1.4 µL of Nextera index [37], 2.8 µL of DNA template, and enough water for a final volume of 14 µL. Thermocycling conditions followed an initial denaturing step at 95 • C for 15 min, followed by 8 cycles of denaturing at 95 • C for 30 s, annealing at 55 • C for 30 s, extension at 72 • C for 30 s, and a final extension at 72 • C for 10 min. Indexed PCR products were  Paired reads were aligned with PEAR [38], and successfully assembled reads went through 'ngsfilter' from OBITools [39] to remove primer sequences and annotate sample information. Trimmed reads were then collapsed into unique sequence variants using 'obiuniq' and denoised with '-cluster-unoise' from VSEARH [40], using default parameters, except for minimum sequence length, which was set as 150 bp for fwh1 and 300 bp for LerayXT. Resulting zero-radius operational taxonomic units (zOTUs) went through chimera removal with '-uchime3_denovo', and were clustered at 99% identity [41] with '-cluster_size'. Finally, reads were mapped to the remaining OTUs with 99% identity using '-usearch_global'. To further remove potential nuclear mitochondrial copies (NUMTs) and surviving PCR and sequencing errors, the R package LULU [42] was used with the default parameters. Extraction and PCR negatives were used to correct for contamination. The maximum number of reads of any OTU identified in either extraction or PCR negative was subtracted from the number of reads observed of that OTU in each sample. OTUs were assigned to a taxon using BOLDIGGER v.1.2.5 [43]. OTUs with a minimum of 90% similarity to a taxon included in the phyla Annelida, Arthropoda, or Mollusca were retained as plausible invertebrate prey [22]. Samples with less than 100 reads assigned to dietary OTUs were discarded, as were OTUs comprising less than 1% of the total dietary reads per sample, so as to avoid errors from tag jumping or overrepresentation of rare prey [44]. Prey items were identified at the genus level, as assignment accuracy at the species level was often missing or undefined in the reference database.
To estimate and compare sampling completeness for each region and fragment, as well as prey richness, we used rarefaction curves based on Hill numbers using the 'iNEXT' function from the iNEXT package in R [45,46]. Prey occurrences for each site and for each fragment were converted into incidence frequencies using 'incfreq', and then sample coverage and prey richness were calculated. Sample coverage gives us the proportion of the diet composed of prey species already sampled, and is considered a better reference than sample size to compare species richness among differently sampled groups [47]. To compare the prey composition among samples of different regions and fragments, we calculated a pairwise distance matrix using the Jaccard dissimilarity indices using 'vegdist' available in the R package vegan [48] to quantify the differences between regions and fragments based on prey occurrence. This matrix was then tested using a permutational multivariate analysis of variance (PERMANOVA) with the Jaccard method and 1000 permutations using 'adonis'. One of the assumptions of PERMANOVA is that there are no differences in dispersion among groups; thus, we further conducted a beta dispersion test using 'betadisper' to confirm this homogeneity. Test results were summarized and displayed via principal coordinates analysis (PCoA). Finally, to assess which prey items were significantly contributing to differences between regions and fragments, we conducted a similarity percentage test using 'simper' with 1000 permutations [49].

Sample Collection and Sequence Amplification
A total of 50 individual fecal pellets were extracted (Morrazo = 32, Porto = 8, Ons = 10), with two replicate extractions from a single pellet collected in Porto, and three extraction negatives. Fragments were successfully amplified in 30 samples with fwh1 (Morrazo = 20, Ons = 10) and 38 samples with LerayXT (Morrazo = 21, Porto = 8, Ons = 9), each including the extraction negatives as well as a PCR negative. Post-filtering, we retained dietary reads from a total of 35 individuals, with 25 samples sequenced at either fragment (83% and 66% success rates for fwh1 and LerayXT, respectively) and 15 samples sequenced at both. Each extraction replicate identified three prey taxa, of which two were common to both replicates, while those prey found in only one replicate comprised less than 5% of the total dietary reads.

Diet Characterization
Across both fragments, a total of 95 unique OTUs were retained, corresponding to 58 prey taxa (Table 1). Two families of Annelida were identified-Almidae and Lumbricidaewherein one and five genera were identified, respectively. Included among these, Lumbricus (present in 49% of samples) was the most common Annelida prey. Arthropoda was the most diverse phylum, comprising 6 known (1 unknown) classes, 18 orders, 32 known (one unknown) families, and 30 known (9 unknown) genera. Among all of these families, no more than two genera were identified. Arthropoda included some of the most common prey-namely, millipedes (Diplopoda), and in particular the genera Polydesmus (present among 51% of all samples), Glomeris (31%), and Ommatoiulus (26%). Mollusca prey comprised only Gastropoda-namely, the orders Pulmonata and Stylommatophora, the former corresponding to a single genus, Cochlicella, and the latter comprising 11 families and 12 genera. The second most common prey overall were roundback slugs of the genus Arion (49%). While in some instances species-level resolution was available for some prey (e.g., all OTUs assigned to the genus Glomeris were also identified as the species G. occidentalis), in many cases, taxonomic assignments were unresolved at the species level (e.g., of the nine OTUs assigned to the genus Arion, seven different published sequences were identified as matches, but all lacked species-level designations). To avoid potentially inflating the prey richness as an artifact of unresolved taxonomic assignments, we opted instead to use the genus-level resolution, which was available for the majority of OTUs identified. Table 1. The frequency of occurrence of prey taxa observed among samples in each region, and in total. Where an OTU at the genus-level resolution could not be identified, the next highest taxonomic resolution (e.g., family, order, etc.) is provided. Significant differences in pairwise comparisons of average abundance between regions are shown in bold with an asterisk (p < 0.05).

Species Richness
When comparing samples sequenced at both fragments, we observed a near-identical sample coverage of 65% and 62% for fwh1 and LerayXT, respectively, with fwh1 detecting a higher number of prey species (39) than LerayXT (25) (Figure 2a). However, when comparing both fragments at similar levels of sampling completeness, the estimated prey richness did not differ significantly between the two fragments (overlapping 95% confidence intervals of rarefaction curves; Figure 2a). A two-sample t-test confirmed that fwh1 produced a higher number of total dietary reads (10,811 ± 9896) post-filtering compared to LerayXT (1802 ± 3015), but no differences in the total number of filtered reads or the ratio of dietary reads to filtered reads (t = 7.0748; p < 0.001).  Sample coverage was 85%, 72%, and 60% for Morrazo, Porto, and Ons, respectively (Figure 2b). Observed prey richness was 29, 25, and 12 for Morrazo, Ons, and Porto, respectively. Rarefaction curves showed that estimated prey richness in Porto was lower than in Morrazo and Ons, while the latter exhibited similar levels of prey richness. The PERMANOVA model showed no differences in the composition of prey identified by either fragment (Figure 3a), but significant differences between regions (p < 0.001). However, the beta dispersal test suggests that the significant differences in prey composition observed in the PERMANOVA between sites may be inflated due to the lack of homogeneity in variance across groups (PCoA; Figure 2b; p = 0.01967). Notable differences in prey prevalence between regions include the absence of millipedes among samples from Ons. This coincides with an increased diversity of soft-bodied prey from among Ons samples, including several gastropod genera-Cochlicella, Oestophora, and Milax-found to be significantly more common, and the only instance of earthworms from the family Alma. Annelida was notably absent among samples from Porto. Several genera of arthropods, however, were significantly more common in Porto than in other locations, including Polydesmida: Oxidus, Diptera: Eupeodes, Hemiptera: Chaitophorus, Lepidoptera: Peridroma, and Isopoda: Armadillidium and Eluma, although we should note that the low sample coverage from Porto may inflate the significance of this observation.

COI Metabarcoding for Salamanders' Diet Characterization
Concordant with previous characterizations of S. salamandra diet, our results from the DNA metabarcoding of COI suggest a prevalence of low-mobility terrestrial detritivores, including millipedes (Diplopoda), roundback slugs (Arionidae), and earthworms (Lumbricidae) [15][16][17]50]. While each of these prey taxa are similarly prevalent in the salamander diet as a whole, regional differences such as the absence of Annelids in Porto and Diplopoda in Ons were observed. Comparing the results of COI metabarcoding to diet inspections by Bas et al. [15] in S. s. gallaica (northwest Iberia), we observed a clear difference in the types of prey detected. The most common prey identified by Bas et al. [15] were Insecta, with far fewer soft-bodied prey compared to the findings of this study. This was expected, as visual inspection may favor the detection of hard-bodied prey that are more slowly digested compared to soft-bodied prey, which may become visually unrecognizable several days after consumption [13]. We observed this anecdotally in the prevalence of Coleoptera identified by Bas et al. [15], of which 16 genera were identified, as compared to the four genera identified in this study, with both studies identifying Nalassus as the most common prey among the class Insecta. Conversely, 18S metabarcoding by Wang et al. [22] found that the most common prey were Gastropoda, which far exceeded other prey in prevalence. Gastropoda have been reported as a common prey elsewhere [16,50], as well as in our results. However, in many of these cases, often only a few prey taxa could be identified, because of either digestion or low taxonomic resolution. Our study identified 58 different prey taxa from three populations at variable taxonomic resolutions-fewer than the 76 prey taxa identified by gastric inspection from a higher number of individuals (N = 72), localities (N = 10), and a wider environmental and ecological survey [15], but nearly threefold more than from previous taxonomic assignment using DNA metabarcoding of 18S across three forest populations in Belgium [22]. Thus, DNA metabarcoding of COI increases taxonomic resolution and provides a cost-effective and expedient method for characterizing the diets of salamanders. Indeed, the increased taxonomic resolution provided by COI suggests that salamanders are able to utilize a variety of Gastropoda prey. This was most evident among our samples from Ons, which often comprised soft-bodied gastropods and a complete absence of Diplopoda-possibly as a result of differences in region, season, or habitat compared to our other sites. For instance, seasonal variation was found to influence prey richness in the diets of salamanders, with the greatest prey richness reported in autumn [4,28]. This seems a probable explanation for the observed prey richness among samples from Ons sampled during autumn. Relatively few arachnid (Arachnida) and centipede (Chilopoda) taxa were identified among our samples, although they have been reported in the diets of salamanders [15,22], suggesting an absence either among our samples or in local abundance.

COI Primers as Barcodes
While studies should always strive to include variable fragments in order to account for primer biases [29,30], our results suggest no discernible difference in the results gathered by using either fwh1 or LerayXT. Greater species richness among fwh1 sequences compared to LerayXT was unexpected given previous comparisons of COI primer performances in literature [51], although sample degradation may favor the shorter fragment. The absence of any clear differences in the prey composition between these primers, however, casts doubt on whether LerayXT underperformed, as inequalities in read output and sample coverage discourage definitive conclusions. Further sequencing and more distributive sampling will be necessary for verification. The number of dietary reads generated by LerayXT was significantly lower than the number generated by fwh1, even with no discernible difference in the prey being identified by either primer. During sequencing, the smaller fragment-in this case, fwh1-will usually be favored [52]; however, a comparison between COI primers found that fwh1 has a higher likelihood of mismatch between the primer and the template, potentially identifying fewer prey species [51]. Despite expectations, LerayXT identified fewer prey species, with a possible explanation being among the unfiltered reads, as 32% of all OTUs and 26% of all reads were 352 bp-longer than the target fragment length-and either unassigned or identified as Flavobacterium. When present, these OTUs were the most abundant reads among a subset of individuals, and may represent an instance of nuclear mitochondrial DNA (NUMT) resulting from transposition of COI into the nuclear genome [53]. Pseudogenes such as these may evade blocking primers and dominate the amplification reaction. No OTUs were assigned to Salamandra salamandra, indicating high efficiency of the blocking primers; however, considering the large size and repetitive nature of the salamander genome, pseudogenes are to be expected [54].

Dietary Variation across Regions
Although differences in prey prevalence were observed between regions, overlap in the prey composition should deter us from drawing any premature conclusions about diet preferences or prey abundance. Instead, we can refer to these preliminary results as a starting point for subsequent studies. Based on our extended rarefaction results (Supplementary Materials Figure S2), we also recommend that future studies aim for a minimum sample size of 20 sample units per site for sample coverage. The differences in prey taxa observed between Ons and the mainland regions-primarily the prevalence of soft-bodied prey such as land snails and slugs (Stylommatophora) and segmented worms (Alma) that were otherwise undetected in the mainland samples-is of particular interest. We might have expected islands to have lower alpha diversity than the mainland [55]; however, samples from Ons were temporally distinct from those of Morrazo and Porto, and must be resampled in the same temporal period in order to control for the known effects that seasonal variation has on prey richness [4,19,50]. Additional observations, such as the relative absence of Diplopoda in the diet of Ons samples, may indicate a scarcity of this common prey, driving prey diversification [56]. In Porto, conversely, we anticipated higher prey richness by taxa-namely, Isopoda and Limacidae-able to utilize anthropogenic spaces [57]. Although prey richness was low when compared to other regions, there was a prevalence of pill woodlice (Armadillidiidae), which may be of interest for studies in ecotoxicology, as terrestrial Isopoda often serve as model organisms in soil ecotoxicology [58]. However, we must also consider that these samples were not sequenced at fwh1, due to poor amplification, and have fewer overall reads to compare (t = 2.1898; p < 0.05). Previous studies investigating prey consumption in S. salamandra detected dietary differences between sexes [22], ages [59], seasons [50], and populations [15]. Follow-up studies should also consider comprehensive sampling of distinct habitats throughout the species' range, as well as the remarkable intraspecific differentiation in reproductive modes [60], head shape [61], and behavioral strategies [62][63][64], both between and within subspecies of S. salamandra. The inclusion of these variables may help to elucidate the factors that contribute to the dietary variation observed in this study.