Diet diversity and environment determine the intestinal microbiome and bacterial pathogen load of fire salamanders

Diverse communities of symbiotic microbes inhabit the digestive systems of vertebrates and play a crucial role in animal health, and host diet plays a major role in shaping the composition and diversity of these communities. Here, we characterized diet and gut microbiome of fire salamander populations from three Belgian forests. We carried out DNA metabarcoding on fecal samples, targeting eukaryotic 18S rRNA of potential dietary prey items, and bacterial 16S rRNA of the concomitant gut microbiome. Our results demonstrated an abundance of soft-bodied prey in the diet of fire salamanders, and a significant difference in the diet composition between males and females. This sex-dependent effect on diet was also reflected in the gut microbiome diversity, which is higher in males than female animals. Proximity to human activities was associated with increased intestinal pathogen loads. Collectively, the data supports a relationship between diet, environment and intestinal microbiome in fire salamanders, with potential health implications.


Results
Diet analysis based on fecal DNA metabarcoding. A total of 20 prey taxa were identified, belonging to the phyla Mollusca, Annelida and Arthropoda.
Fire salamander ingested prey biomass was determined by the relative number of sub-OTU reads for each prey category. We used this as a proxy for the relative amount of ingested prey within all salamander fecal samples. More than 70% of the sub-OTU reads were identified as gastropods, followed by millipedes at 10.8%, centipedes at 5.8%, and soil mites at 4.4%. All other taxa were rare in fire salamander diet. Diet taxa prevalence (presence/ absence) across all salamanders revealed gastropods to be the most prevalent taxonomic class at 40.6%. Both millipedes and annelids (earth worms) were found at 14.5%. Centipedes, true flies and crustaceans made up 5.8% of salamander diet. Relative abundance of prey sub-OTUs differed significantly between females (n = 12) and males (n = 21) (RDA F = 4.58, df = 2.28, p = 0.011). Prey ingested by females was made up largely of gastropods (73.3%), with millipedes making up most of the remaining prey (16.7%). Males, on the other hand, had a more balanced diet, with a large portion of their ingested prey being divided between gastropods (36.4%), millipedes (21.9%), soil mites (16.7%), and centipedes (11.6%) (Fig. 1a). Prevalence of prey taxa also differed significantly between both sexes (χ 2 1,12 = 35.12, p < 0.001) (Fig. 1b). Fire salamander diets were largely made up of the same general taxa, and the relative abundances of prey reads did not differ between the forests of HG (Heilig Geestgoed) (total n = 11, male n = 5, female n = 5, juvenile n = 1), M (Makegem) (total n = 17, male n = 9, female n = 5, juvenile n = 3) and S (Smetledebos) (total n = 9, male n = 7, female n = 2, juvenile n = 0) (RDA F = 0.50, df = 2.28, p = 0.738). Gastropods and earthworms were found to have a similar prevalence in salamander diets between forests (Supplementary Figure S1). Differences were seen in millipedes, where values fluctuated between locations, although not significantly (χ 2 2,12 = 4.35, p = 0.113). Many of the Insecta taxa were not part of the salamander diet in all three forests (Supplementary Figure S1).
Based on redundancy analysis (RDA), the gut microbiome composition was significantly affected by location (RDA F 2,36 = 4.07, p = 0.001) but not by sex (RDA F 2,36 = 0.85, p = 0.720). Alpha diversity of the gut microbiome varies between locations (Fig. 3b), with significant differences for ASV richness (χ 2 = 7.52, df = 2, p = 0.023), Chao1 (χ 2 = 9.18, df = 2, p = 0.010) and Shannon Index (χ 2 = 11.09, df = 2, p = 0.004). We found that forest HG has significant lower alpha diversity than the other two forests (Fig. 3b). When comparing between sexes, alpha diversity of the gut microbiome is significantly lower for females than males: ASV richness (W = 43, p = 0.013), Chao1 (W = 48, p = 0.021), Shannon Index (W = 38, p = 0.005). The dissimilarity in microbiome composition was analyzed with PERMANOVA on Bray-Curtis distance, where no difference could be seen comparing between sexes (Fig. 2c, Pseudo-F 2,29 = 0.9, p = 0.556), but when comparing between locations, the three forests were significantly different from each other (Fig. 2d www.nature.com/scientificreports/ We included pathogens of the genera Flavobacterium, Chryseobacterium, Sphingobacterium, Aeromonas, Citrobacter, Yersinia, Acinetobacter and Stenotrophomonas as a proxy for pathogen load, given their known involvement in amphibian pathology [45][46][47][48][49][50][51][52][53][54][55] . The presence of potential pathogenic amphibian bacteria in the fire salamander gut microbiome was analysed at genus level (Supplementary Figure S2). We focused on eight genera well known to contain pathogenic bacterial taxa ( Table 1). The total relative abundance of the eight pathogenic bacterial genera makes up less than two percent of the fire salamander gut microbiomal communities in forest M (0.5%) and S (1.4%), but more than 13% in forest HG. The relative abundance of the pathogen load is significantly higher in forest HG, compared to forest M (p < 0.0001). Similarly, a tendency towards increased pathogen load was observed when comparing forests HG and S (p = 0.8649). Comparison of beta diversity of the pathogen load between forests was analysed with PERMANOVA on Bray-Curtis distance and revealed a significant difference between locations (Supplementary Figure S3

Diet and gut microbiome. A Mantel test revealed that dietary differences between individual salamanders
are not similar to the microbiome differences (r = − 0.05, p = 0.747). Exploring the structure between diet and microbiome we found a co-inertia coefficient of 0.42, indicating that the two categories varied independently (p = 0.721). Furthermore, no correlation was observed between the alpha diversity (Chao1 index) of diet and gut microbiome (Supplementary Figure S4  www.nature.com/scientificreports/

Discussion
To estimate the accuracy of the DNA barcoding technique, validation studies have been conducted using feeding experiments of captive animals, where exact diet inputs were compared to fecal DNA sequencing outputs [56][57][58] . Whilst the species of prey could be successfully identified in these studies, the proportion of detected DNA varied and was not completely on par with the dietary proportions. Amounts of DNA in the fecal matter need to be exactly proportionate to ingested prey biomass, which is not always the case, due to biological and technical biases, such as different speeds of digestion, size of ingested prey, possible presence of multicopy genes, DNA degradation in fecal samples and availability of DNA reference sequences of potential prey in public databases [59][60][61] . Nonetheless, DNA sequencing of fecal matter is a more viable technique to identify soft-tissued, easily digestible prey 59 , and able to provide higher resolution than conventional stomach content analysis 62,63 . This technique has been applied successfully in unravelling diets for many species, including birds 44,56 , mammals 59,64,65 , fish 66,67 and reptiles 68 , but so far has not been used for amphibians.
In this study, we characterise the diet of fire salamanders by means of DNA metabarcoding of fecal samples. Our data identifies gastropods as the most prevalent prey of fire salamanders, and reveals the diet differences between sexes. As reported previously, fire salamanders are opportunistic and generalist predators that focus on slow-moving and soft prey items [69][70][71][72] . Our findings are in line with previous publications of fire salamander diet compositions, where Gastropoda was also found to be the most important prey 71 . Other amphibian species have significantly different diet compositions. Monte Albo cave salamander, for example, mostly consume Hymenoptera, Ambrosi's cave salamander feed mostly on Arachnida and Strinati's cave salamander mostly consume Diplopoda [73][74][75] . This likely explains why our results are in line only with previous studies on diet conducted specifically on fire salamanders. The sex-dependent change in diet has been thoroughly investigated in a number of amphibian species 73,75-78 , but no difference of diet composition has been found. Our study observed for the first www.nature.com/scientificreports/ time diet differences between sexes in amphibians, with female fire salamanders consuming more gastropods and millipeds. While in other vertebrates, such as reptiles, this is often linked to sexual dimorphism 79-81 and/ www.nature.com/scientificreports/ or the utilisation of different microhabitats 82 , we attribute the observed difference to different activity patterns.
In North American plethodontid salamanders it has been found that females are less active than males [83][84][85] , resulting in them eating more slow-moving organisms, such as millipedes and gastropods 86 . Similarly, during the breeding season in autumn, female fire salamanders are less active outside the shelters than males 82 , which could explain the diet difference between sexes. Interestingly, the alpha diversity of the gut microbiome was also found to be significantly different between sexes, with diversity of species being lower in female than in male animals. Looking at this lower diversity of the gut microbiome in female animals, especially with regards to the sex-dependent diet, these findings suggest that decreased prey variation results in lower gut microbiome diversity. Previous studies on gut microbiome compositions between sexes in amphibians have yielded mixed results. One study on Chinese concave-eared frogs observed a significant difference in gut microbiome composition at family and genus levels between sexes 87 . In cane toads a study elucidated that the difference in variation of gut microbiome communities is mostly due to the factor of sex, but did not statistically quantify the effect on microbiome diversity 88 . In rice field frogs, however, no statistical difference of gut microbiome composition could be found between sexes 89 . Therefore, differences of gut microbiome diversity between sexes depends on the respective species, and is most likely subject to individual behaviour, prey selection, biology and environmental factors.
The alpha diversity of the fecal micriobiome was also found to be significantly different between locations, with a lower diversity of species found at forest HG, which is located adjacent to farmland, unlike the other two forests in this study. Forest HG also stood out, in that it had a remarkably higher pathogen load, compared to the other forests. We assume that the higher pathogen load coincided with the close vicinity to farmland. This hypothesis is supported by the relatively high amount of Proteobacteria, and the presence of members of the phylum Elusimicrobia in these salamanders, which has been associated with exposure to fertilizers and pesticides in farmland frogs 90 . Generally, a substantial proportion of Proteobacteria is typical for amphibian gut microbiomes, but more so in larvae than adults 17,91,92 . An additional explanation is that, in contrast to the other two sites, the salamanders from this forest are exposed to raw sewage from the neighbouring houses, which may also explain increased numbers of Yersiniaceae, notably of the genus Yersinia, which has been previously reported as an indicator of fecal pollution in fish 93 . Yersinia has also been reported in amphibian species, such as Necturus maculosus 94 , Rana pipiens 95 and Rana clamitans 45 . This, however, is the first time Yersinia has been found in gut microbiome of amphibians from a farmland adjacent habitat, which may prossibly be linked to sewage and fecal pollution. Therefore, the salamander gut microbiome is likely associated with land use and/or pollution, which may have consequences for salamander physiology and health.
Furthermore, we observed no correlation between diet alpha diversity and gut microbiome alpha diversity in fire salamanders. However, previous studies on the correlation between microbiome richness and dietary richness did not yield consistent results. A positive correlation between diet diversity and gut microbome diversity was found in humans 96 , black howler monkeys 97 , kudus 98 and predatory insects 99 . In contrast, a negative correlation was observed in fish, where a study showed that feeding on mixed diets resulted in a lower gut microbiome diversity, compared to pure diets 11 . Moreover, no correlation between diet diversity and gut microbiome diversity was observed in numerous mammalian herbivore species, such as pika, elephant and camel 98 . Diet-microbiome correlation is thus difficult to distill into a general rule across the animal kingdom. A previous study looked specifically at the diet and gut microbiome correlation, where a positive correlation was found in tadpoles of Malagasy frogs 17 . We were unable to confirm this in our study. So for now it appears that variety of diet does not always correlate with the variety and biodiversity of the gut microbiome in amphibians.
Similar studies of the gut microbiome have been conducted in fire salamander larvae 20 , where they found Proteobacteria, Firmicutes and Bacteroidetes to be the most abundant phyla. This is in line with our findings, where we have found Bacteroidetes, Firmicutes and Proteobacteria to be the most abundant phyla. Exception is the phylum Actinobacteriota, which was found to be more abundant in both the pond larvae (3.9%) and stream larvae (6.6%), than in our adult fire salamanders (0.1%). Nonetheless, these findings could indicate that the most abundant phyla in the gut microbiome composition do not undergo major changes throughout the different life stages of fire salamanders. The aforementioned publication also found that gut microbiome composition Table 1. Mean relative abundance and level of significance of potentially pathogenic bacterial genera in the fire salamander gut microbiome between different forests (HG = forest Heilig Geestgoed, M = forest Makegem and S = forest Smetledebos). www.nature.com/scientificreports/ changes depending on habitat-specific diet in ponds versus streams. We were unable to see such stark contrasts of habitat-dependent diet in our adult animals, since prey was similar across different forest habitats. To conclude, this is the first study to investigate the correlation between the diversity of diet and gut microbiome of adult fire salamanders in Belgian forests, using high-throughput DNA metabarcoding techniques. We show that diet composition is driven by sex, and influences microbiome composition in the fire salamander intestine. However, no correlation was observed between diet diversity and gut microbiome diversity.

Methods
This study was carried out in compliance with the ARRIVE guidelines and all methods were carried out in accordance with relevant guidelines and regulations. Sample collection. Fire salamanders were collected from each of these three forests. They were kept individually without food in sterile boxes (16 × 11 × 5 cm) with moist towel, air holes and hiding places for 1 to 3 days at 15 °C at the Faculty of Veterinary Medicine, University of Ghent in Merelbeke and then returned to their exact same locality in the forest. All animals were weighed (to 0.1 g) and measured (to 0.1 cm). Scaled mass index (SMI) of body condition ( M i ) was calculated for all salamanders: . Animals were checked and the boxes were cleaned each day and any fecal samples were collected in Eppendorf tubes and frozen at − 20 °C until DNA extraction. In total, 60 fecal samples were collected from 49 individual fire salamanders. 13 individual fire salamanders were collected from forest HG (7 males, 5 females and 1 juvenile), 21 individuals from forest M (11 males, 6 females and 4 juveniles) and 15 individuals from forest S (9 males, 4 females and 2 juveniles). Diet. Sequencing methods were used as previously described 101 . The V9 regions of nuclear 18S rRNA gene were amplified using modified primers 1391F (5′-GTA CAC ACC GCC CGTC-3′) and EukBr (5′-TGA TCC TTC  TGC AGG TTC ACC TAC -3′) 102 . We performed index PCR with unique combinations of indexed forward and reverse primers for each sample using Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific). For quantification of amplicon DNA concentration, we roughly assessed intensity of their signal on agarose gels and then added 2 to 6 µl to a pooled library. This library was gel-purified by cutting out the band of the correct amplicon size, and subsequently cleaned with Qiagen MinElute Kit. DNA concentration of the purified library was determined with a Qubit 2.0 fluorometer. The library was sequenced on the Illumina MiSeq platform using MiSeq Reagent Kit v3 for 300 cycles in both directions, including 10% phiX. Sequences were processed with MacQIIME v1.9.1 103 , filtering the forward reads as previously described 101 . Quality filtered sequences were clustered into sub-operational taxonomic units (sub-OTUs) using the deblur workflow 103,104 ; (https:// github. com/ bioco re/ deblur). Within this workflow, all sequences were trimmed to 150 bp and subsampled to 13,000 reads per sample. Sub-OTUs clusters with less than 10 reads were removed. Sequences of sub-OTUs were taxonomically identified through BLAST searches. Only potential salamander prey items were selected from the 18S data, discarding sub-OTU reads from taxa such as Bacteria, Fungi and Protists. Prey taxa were identified into the main taxa: Gastropoda, Diplopoda, Chilopoda, Arachnida (Acari and Opiliones), Crustacea, Annelida (Lumbricidae), and six insect clades: Hymenoptera, Hemiptera, Dipera, Coleoptera, Blattodea and Collembola (Supplementary  Table S2).
Gut microbiome. The V3-V4 hypervariable region of the 16S rRNA gene was amplified using gene-specific primers S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21 105 . The PCR amplification procedures were performed according to a previous study 106 . The final barcoded libraries were combined to an equimolar 10 nM pool and sequenced with 30% PhiX spike-in using the Illumina MiSeq v3 technology (2 × 300 bp, paired-end) at the Oklahoma Medical Research Centre (Oklahoma City, USA). Demultiplexing of the amplicon dataset and deletion of the barcodes was done by the sequencing provider. Quality of the raw sequence data was checked with the FastQC quality-control tool (Babraham Bioinformatics, Cambridge, United Kingdom; http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/), after which the sequences were trimmed, quality-filtered and www.nature.com/scientificreports/ dereplicated using the DADA2 algorithm (v1.14) within R 107 . An initial amplicon sequence variant (ASV) table was constructed before chimaeras were identified using the removeBimeraDenovo function. Finally, taxonomy was assigned using DADA2's native naïve Bayesian classifier against the Silva database (v138) 108 .
To select the appropriate subsampling depth, alpha rarefaction curves were generated ( Supplementary Figure S5). One sample (sample MF12) was unsufficiently sequenced, and was therefore excluded from the final ASV table. Amplicon sequence variants were subsampled (rarefied) by random sampling to a depth of 5047 ASVs per samples (depth of the least sequenced sample).
As our sequencing samples are high bacterial biomass faecal samples, negative controls (buffer control) were included in our sequencing runs, to guard against reagent contamination, as this is only a problem in low bacterial biomass samples. We rarefied to 1300 reads per sample for 18S and 5047 ASVs for 16S and eliminated all sub-OTUs with fewer than 10 reads overall. 43 fecal samples from 37 individuals were used for subsequent diet statistic analysis and 41 fecal samples from 35 individuals were used for gut microbiome analysis. The mean abundance of 16S and 18S data of multiple samples from same individuals was calculated.
Statistical analysis. Sequence reads per fecal DNA sub-OTU found in the samples of each salamander were used as approximate proxy for ingested biomass. We calculated the prevalence of each prey taxon in the salamander diet, determined by the number of samples in which a prey category was found present or absent.
To evaluate alpha diversity (in-samples diversity) of salamander diets, OTU richness and Chao2 were calculated from prey prevalence data, and for bacteria ASV richness, Chao1 (estimated ASV richness) and Shannon indexes (estimated community diversity) were calculated from relative abundance data using the fossil package v 3.2.5 109 . To compare alpha diversity between location and sexes, a Kruskal Wallis One-Way Analysis of Variance was performed between locations and a Wilcoxon rank sum test was performed between sexes (due to non-normal distribution in both cases). To evaluate the differences of the prey prevalence data between different locations and sexes, Pearson's Chi squared was calculated. Differences in microbial relative abundance between locations were assessed at the phylum, family and genus level, and differences between sexes were assessed at the phylum and family level, using the DESeq2 algorithm from the phyloseq package (v 1.30) 110 . Differentiallly abundant taxa were identified through applying the negative binomial Wald test with p-values corrected for multiple hypothesis testing using the Benjamini-Hochberg method 111 . The fire salamander gut microbiomes contained at least eight potentially pathogenic genera (Flavobacterium, Chryseobacterium, Sphingobacterium, Aeromonas, Citrobacter, Yersinia, Acinetobacter, Stenotrophomonas). The differences in relative abundances of pathogen load between locations were assessed using SPSS (IBM SPSS Statistics for Windows, Version 26.0. Armonk, NY, USA), with Kruskal-Wallis analysis followed by Dunn's multiple comparison tests (significance values adjusted by the Bonferroni correction for multiple tests).
Redundancy analyses (RDA) were performed to examine the effect of location and sex on the relative abundance of ingested biomass, as well as to further quantify the effect of location, sex, SMI and diet on gut bacteria. Spearman's Rank correlation was used to test whether the alpha diversity of individual diet and gut bacteria was significantly related to salamander SMI. Additionally, correlations were used to examine the relationship of alpha diversity of diet on gut bacteria.
To explore the beta diversity, principal coordinates analysis was performed, based on the Jaccard distance 112 , to compare community dissimilarities with presence or absence of OTUs for diet. We performed the principal coordinates analysis based on Bray-Curtis distance, to compare the community dissimilarities of relative abundance of ASVs, for both gut microbiome and pathogen load. We used PERMANOVA to test if the divergent species composition differed significantly between locations and sexes. Analyses were done with the package vegan (v 2.4) 113 and the functions vegdist and adonis 2 in R (v 3.4).
To measure the correlation between individual salamander diet and the fecal microbiome, Mantel tests were performed on the Jaccard and Bray-Curtis distance matrices using R. To further explore the inter-relationship between the diet and the bacterial community matrices, a co-inertia analysis was run using the ade4 package 114 in R.