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

**Adam J. D. Marques 1,2,3,\*, Vanessa A. Mata 1,2 and Guillermo Velo-Antón 1,2,4,\***


**Abstract:** 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.

**Keywords:** COI; diet; DNA metabarcoding; prey; salamanders

#### **1. 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–18]. Fecal content inspection is less invasive, but introduces a bias favoring hard-bodied prey species that are not fully digest [13].

**Citation:** Marques, A.J.D.; Mata, V.A.; Velo-Antón, G. COI Metabarcoding Provides Insights into the Highly Diverse Diet of a Generalist Salamander, *Salamandra salamandra* (Caudata: Salamandridae). *Diversity* **2022**, *14*, 89. https://doi.org/ 10.3390/d14020089

Academic Editor: Sebastiano Salvidio

Received: 31 December 2021 Accepted: 24 January 2022 Published: 28 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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–33].

#### **2. 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

purified with AMPure XP beads (Beckman Coulter), eluted to 25 μL, and pooled into equimolar concentrations per fragment. The pooled libraries were quantified with qPCR, normalized to 4 nM, and sequenced in an Illumina MiSeq with an expected coverage of 20,000 reads per sample.

**Figure 1.** Maps include the species distribution of *Salamandra salamandra*, our study region in the northwest Iberian Peninsula, and photos illustrating the habitats found within the three sampled regions of Ons, Morrazo, and Porto.

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].

#### **3. Results**

#### *3.1. 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.

#### *3.2. 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 Lumbricidae wherein 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).


*3.3. 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).

**Figure 2.** Species diversity extrapolated as a factor of sample coverage for (**a**) different CO1 fragments; (**b**) sampling regions. Points indicate the observed richness, while error bars indicate the 95% confidence intervals of the response curves.

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.

**Figure 3.** Variance in the prey taxa summarized by the first two eigenvectors of a PCoA for (**a**) different CO1 fragments; (**b**) sampling regions. Ellipses capture the distance of points from the centroid within one standard deviation.

#### **4. Discussion**

#### *4.1. 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–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.

#### *4.2. 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].

#### *4.3. 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–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.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/d14020089/s1, Figure S1: Rarefaction comparison between fragments; Figure S2: Rarefaction comparison between regions; Table S1: Extended list of OTUs after sequence annotation and filtration.

**Author Contributions:** Conceptualization, G.V.-A. and V.A.M.; methodology, V.A.M.; validation, A.J.D.M. and V.A.M.; formal analysis, A.J.D.M.; investigation, G.V.-A. and A.J.D.M.; resources, G.V.- A. and V.A.M.; data curation, G.V.-A., A.J.D.M. and V.A.M.; writing—original draft preparation, A.J.D.M.; writing—review and editing, G.V.-A. and V.A.M.; visualization, A.J.D.M.; supervision, G.V.-A. and V.A.M.; project administration, G.V.-A.; funding acquisition, G.V.-A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Fundação para a Ciência e a Tecnologia, grant numbers PTDC/BIA-EVL/28475/2017 and PTDC/BIA-CBI/2278/2020. V.A.M. was funded by the FCT through the program 'Stimulus of Scientific Employment, Individual Support—3rd Edition' (2020.02547.CEECIND). G.V.-A. was supported by the FCT (CEECIND/00937/2018), and recently by a Ramón y Cajal research grant (Ref. RYC-2019-026959-I/AEI/10.13039/501100011033). This work was co-funded by the project NORTE-01-0246-FEDER-000063, supported by Norte Portugal Regional Operational Programme (NORTE2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF).

**Institutional Review Board Statement:** Fieldwork for collecting individuals and fecal samples was carried out with the corresponding permits from the regional and national administrations (Xunta de Galicia, Ref. EB-031/2021; Portugal Ref. 544/2021/CAPT). Sampling procedures were carried out following the Guidelines for Use of Live Amphibians and Reptiles in Field and Laboratory Research, 2nd Edition, revised by the Herpetological Animal Care and Use Committee (HACC) of the American Society of Ichthyologists and Herpetologists, 2004.

**Data Availability Statement:** The data presented in this study are openly available in Zenodo at https://zenodo.org/badge/latestdoi/429739169 (accessed on 30 December 2021).

**Acknowledgments:** We thank Lucia Alarcón-Ríos and Iria Pazos for assistance during fieldwork, and the staff of the *Parque Nacional Marítimo-Terrestre das Illas Atlánticas de Galicia* for providing lodging and facilitating the trips to the island.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Ecological Observations on Hybrid Populations of European Plethodontid Salamanders, Genus** *Speleomantes*

**Enrico Lunghi 1,2,3,4,\*, Fabio Cianferoni 2,5, Stefano Merilli 6, Yahui Zhao 1,\*, Raoul Manenti 7,8, Gentile Francesco Ficetola 7,9 and Claudia Corti <sup>2</sup>**


**Abstract:** *Speleomantes* are the only plethodontid salamanders present in Europe. Multiple studies have been performed to investigate the trophic niche of the eight *Speleomantes* species, but none of these studies included hybrid populations. For the first time, we studied the trophic niche of five *Speleomantes* hybrid populations. Each population was surveyed twice in 2020, and stomach flushing was performed on each captured salamander; stomach flushing is a harmless technique that allows stomach contents to be inspected. We also assessed the potential divergence in size and body condition between natural and introduced hybrids, and their parental species. Previously collected data on *Speleomantes* were included to increase the robustness of these analyses. In only 33 out of 134 sampled hybrid *Speleomantes* we recognized 81 items belonging to 11 prey categories. The frequency of empty stomachs was higher in females and individuals from natural hybrid populations, whereas the largest number of prey was consumed by males. We compared the total length and body condition of 685 adult salamanders belonging to three types of hybrids and three parental (sub)species. Three group of salamanders (one hybrid and two parental species) showed significantly larger size, whereas no difference in body condition was observed. This study provided novel ecological information on *Speleomantes* hybrid populations. We also provided insights into the potential divergence between hybrids and parental species in terms of size and body condition. We discuss our findings, and formulate several hypotheses that should be tested in the future.

**Keywords:** *Speleomantes*; *Hydromantes*; trophic niche; body condition; cave biology; biospeleology; parental species; diet; size; capture-mark-recapture

#### **1. Introduction**

European cave salamanders of the genus *Speleomantes* are the only plethodontids living in the Europe, and are almost all endemic to Italy [1]. Seven of the eight *Speleomantes* species (*S. ambrosii*, *S. italicus*, *S. flavus*, *S. supramontis*, *S. imperialis*, *S. sarrabusensis*, *S. genei*) live exclusively in Italy, while the range of one species (*S. strinatii*) extends to part of French Provence [1]. Each species is distributed in a well-defined area, and no range overlap exists; *Speleomantes* distribution is likely shaped by geomorphology [2,3]. For *S.*

**Citation:** Lunghi, E.; Cianferoni, F.; Merilli, S.; Zhao, Y.; Manenti, R.; Ficetola, G.F.; Corti, C. Ecological Observations on Hybrid Populations of European Plethodontid Salamanders, Genus *Speleomantes*. *Diversity* **2021**, *13*, 285. https:// doi.org/10.3390/d13070285

Academic Editors: Salvidio Sebastiano and Michael Wink

Received: 25 May 2021 Accepted: 17 June 2021 Published: 23 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

*ambrosii*, the river Magra marks the separation of the two allopatric subspecies: the western *S. ambrosii ambrosii* and the eastern *S. a. bianchii* [1,4]. Although some phenotypic variability can be observed among *Speleomantes* [5–7], a valid method to phenotypically distinguish among species/subspecies is still lacking, thus their identification is based mostly on geography [1,3]. The southern distribution of the latter reaches the norther limit of *S. italicus*, creating a narrow contact zone in which viable hybrid populations are found [1,8] (see Figure 1 of Ref. [8]) These natural hybrid populations exhibit some genetic divergence, which is mainly influenced by the relative abundance of individuals of each parental species: in the northern part of the range there are populations of *S. a. bianchii* introgressed with *S. italicus*, whereas in the southern area, the opposite occurs [8]. However, hybrids do not show a clear divergence of phenotypic characters from their parental species; therefore, hybrids can also only be recognized on the basis of their geographic distribution [8,9]. In addition to this natural hybrid zone, hybrids between *S. ambrosii* and *S. italicus* are also related to a human-mediated translocation. In 1983, for scientific purposes individuals of both *S. italicus* and *S. a. ambrosii* were introduced into a natural cave in southern Tuscany, outside the natural range of *Speleomantes* [1,10]. Thirty years after its introduction, the population was genetically characterized and none of the individuals had a pure genotype. For 77% of individuals, the majority of alleles (>75%) matched alleles specific to *S. a. ambrosii*; for 6% of individuals, the majority of alleles matched those of *S. italicus*; and for 16% of individuals, alleles of *S. a. ambrosii* and of *S. italicus* were recombined [11]. However, the lack of ecological information on this population prevents us from evaluating the potential difference between these non-natural hybrids and the hybrids living in the natural hybrid zone, and to assess the trophic relationships between these populations and the local fauna [8,12]. Here we provide the first assessment of the size, body condition, and diet of both natural and introduced hybrid populations of *Speleomantes*.

**Figure 1.** Map indicating the sampled populations. Blue labels indicate natural hybrid populations, whereas yellow indicates those introduced outside the natural distribution of *Speleomantes*. The coordinates of the surveyed sites are not provided to increase species protection [13].

#### **2. Materials and Methods**

In 2020 we surveyed five hybrid populations of *Speleomantes* that inhabit natural caves (Figure 1); each population was surveyed twice, once before (29 June–14 July) and after (3–10 September) the aestivation period [1].

Three populations are located within the natural hybrid zone between *S. italicus* and *S. ambrosii bianchii* occurring in north-western Tuscany (Province of Lucca); two populations included introgressed *S. italicus* with >10% of *S. ambrosii* alleles, whereas the other is a *S. ambrosii* population introgressed with >10% of *S. italicus* alleles [8]. The natural hybrid populations were selected following the genetic characterization of Ruggi, Cimmaruta, Forti and Nascetti [9]. The two other hybrid populations are found in southern Tuscany (Province of Siena), where cave salamanders are not native and have been introduced [1].

One of these two populations was discovered during this study. In the surroundings of the cave where *Speleomantes* were introduced [10], some authors (CC, EL, SM) have found and explored another cave about 53 m in a straight line from the previous one; because several *Speleomantes* were observed during the exploration, the individuals from this second cave were considered to be an additional population [14]. Surveys and animal handling were performed by taking measures to avoid the spread of pathogens (i.e., using disposable gloves and disinfecting boots and equipment after each survey). During each survey, all captured salamanders were placed in a perforated plastic box (60 × 40 × 20 cm). After salamanders were captured, we recorded the following data in sequence: individual sex (males were recognized by the presence of the distinctive mental gland, females were salamanders with SVL ≥ 40 mm without the mental gland, and remaining salamanders were considered juveniles) [1]; weight (using a digital scale, 0.01 g); a photo was taken from the dorsal view of individuals positioned along a reference card [6]; the harmless stomach flushing technique was used to evaluate *Speleomantes* foraging activity [15]. All salamanders were released at their collection points. The photos were analyzed with ImageJ software to measure the total length (TL) of the salamanders and estimate the snoutvent length (SVL) [16]. Stomach contents were preserved in a 75% ethanol solution and subsequently observed under an optical microscope, where prey items were recognized and counted following [17]. When we were unable to recognize any item at the order level, we considered the content as "unidentifiable"; when no remains were found, the stomach was considered "empty". The dorsal pattern of salamanders was used for individual recognition [18].

We used binomial Generalized Linear Mixed Models (GLMMs) to assess the potential effects that the considered variables may have on the stomach condition. Four individuals were captured twice, so we only used data from their first capture. We used stomach condition (empty/full) as the dependent variable, and the independent variables were salamander sex (male, female, juvenile) and hybrid identity (natural = *S. a. bianchii* × *S. italicus*; introduced = *S. a. ambrosii* × *S. italicus*); because the frequency of empty stomach changes over time [19], the survey month was added, together with the identity of the cave, as a random variable. Similarly, we used a GLMM to assess the potential correlation between the number of recognized prey and the independent variables mentioned above; cave and individual identity were still used as random variables.

We also performed a comprehensive analysis by combining these data with previously published datasets [5,8,20] to evaluate differences in size and body condition between hybrids and parental species. In these previous studies, the threshold to discriminate between adults and juvenile was set on the base of the smallest TL measured among males, which was 69 mm; in this study the smallest male had a SVL of 44 mm and TL of 69 mm. For each site, only the survey with the highest number of measured salamanders was considered, with the exception of hybrid data, in which pattern recognition [18] allowed the inclusion of individuals captured for the first time during the second survey. We used a Linear Mixed Model (LMM) to evaluate the potential differences in adult size (TL) between different groups of salamanders (*S. a. ambrosii*, *S. a. bianchii*, *S. italicus*, *S. a. ambrosii* × *S. italicus*, *S. a. bianchii* × *S. italicus*). A Shapiro–Wilk test showed a non-normal distribution of data related to *Speleomantes* size (log-transformed TL, W = 0.98, *p* < 0.001); however, LMM is appropriate for the analysis of non-normal distributed data [21]. The log-transformed TL was used as a dependent variable, and the group of salamanders as the independent variable. Considering the natural divergence in maximum size occurring between males and females [1], the sex of salamanders was used as a random factor. We used the Residual Index (RI) as a proxy of the body condition of the salamanders; this index provides information on the difference between the observed and the expected body mass [22,23]. To calculate the RI, we first log-transformed weight and TL, and then extracted residuals from the regression analysis for each species/hybrid group, in order to avoid bias due to different size [22,23]. *Speleomantes* body condition peaks during the foraging periods (i.e., when precipitation is higher and temperature relatively cold), and is

poorest during inactivity periods (i.e., when climatic conditions are too hot and dry), when the salamanders mostly consume energy that was previously stored [1,24,25]. Therefore, in this analysis we only used data collected before the aestivation period (June–July). We used an LMM to assess the potential correlation between the RI (dependent variable) and two independent variables: the sex of the salamanders and the species group. Considering that data were collected over different periods, we included the month and year of the survey, in addition to the identity of the cave, as random factors. The significance of GLMM and LMM variables was tested with a likelihood ratio test.

#### **3. Results**

In the hybrid populations, we obtained 138 salamander detections corresponding to 134 individuals; four individuals (one male and three females) were observed twice. The size of the recaptured individuals (both SVL and TL) did not change between the two surveys: the difference between the first and second measurement was <2 mm; this difference is comparable to measurement error [16]. Most of the individuals sampled (90) had an empty stomach. The frequency of empty stomachs was significantly different between sexes (χ<sup>2</sup> = 7.31, df = 2, *p* = 0.026) and type of hybrid (*ß* = 4.41, SE = 1.14, χ<sup>2</sup> = 14.91, df = 1, *p* < 0.001); the frequency of empty stomachs was higher in females and individuals from natural hybrid populations. Eleven individuals had stomach contents in a state of advanced digestion and, therefore, the contents were considered not identifiable. We were able to recognize a total of 81 prey items from 33 individuals; the recognized prey belonged to 11 different categories: Sarcoptiformes (1), Mesostigmata (1), Araneae (4), Pseudoscorpiones (4), Polydesmida (3), Isopoda (8), Hemiptera (1), Hymenoptera (1), Coleoptera (3), Coleoptera\_larva (1), and Diptera (55). Diptera were the most consumed prey: they were observed in 32 individuals, representing 67% of recognized prey. The number of prey consumed was significantly affected by the sex of salamanders (χ<sup>2</sup> = 10.74, df = 2, *p* = 0.005); the largest number of prey was consumed by males.

When we combined the data of hybrid populations with those from previous surveys, we obtained data on the size and body condition of 678 salamanders (additional details in Figure 2). The maximum and average (±SD) size (TL) measured for the salamander groups considered were: *S. a. ambrosii*, females max. 125 and average 91 (±14) mm, and males max. 105 and average 92 (±6) mm; *S. a. bianchii*, females max. 114 and average 89 (±12) mm, and males max. 124 and average 105 (±12) mm; *S. italicus*, females max. 120 and average 94 (±13) mm, and males max. 118 and average 101 (±8) mm; *S. a. bianchii* × *S. italicus*, females max. 119 and average 90 (±15) mm, and males max. 119 and average 95 (±8) mm; *S. a. ambrosii* × *S. italicus*, females max. 134 and average 98 (±20) mm, and males max. 108 and average 103 (±2) mm. The size of the adult salamanders significantly differed between the groups (F4,506 = 4.61, *p* = 0.001); *S. italicus*, *S. a. bianchii*, and *S. a. ambrosii* × *S. italicus* hybrids had the largest size (Figure 2). No significant differences in body condition were observed between sexes (F2,641 = 0.69, *p* = 0.5) or species groups (F4,25 = 0.16, *p* = 0.956).
