Sonia L. Ghose1,2*
Jonathan A. Eisen1,2,3- 1Genome Center, University of California, Davis, Davis, CA, United States
- 2Department of Evolution and Ecology, University of California, Davis, Davis, CA, United States
- 3Department of Medical Microbiology and Immunology, University of California, Davis, Davis, CA, United States
Introduction: The amphibian skin microbiome is an important line of defense against pathogens including the deadly chytrid fungus, Batrachochytrium dendrobatidis (Bd). Bd is known to preferentially infect ventral skin surfaces and feet of host amphibians, often leaving dorsal surfaces like the back uninfected. Within-individual variation in infection distribution across the skin, therefore, may relate to differences in microbiomes among skin regions. However, microbiome heterogeneity within amphibian individuals remains poorly characterized.
Methods: We utilized 16S rRNA gene amplicon sequencing to compare microbiomes of 10 body regions from nine captive Rana sierrae individuals and their tank environments. These individuals were naive to Bd, allowing us to assess whether microbiomes differed among body regions prior to any impacts that may be caused by infection.
Results: We found that frog skin and tank environments harbored distinct microbial communities. On frog skin, the bacterial families Burkholderiaceae (phylum Proteobacteria) and Rubritaleaceae (phylum Verrucomicrobia) were dominant, driven in large part by relative abundances of undescribed members of these families that were significantly higher on frogs than in their environment. Within individuals, we detected differences between microbiomes of body regions where Bd infection would be expected compared to regions that infrequently experience infection. Notably, putative Bd-inhibitory relative abundance was significantly higher on body regions where Bd infection is often localized.
Discussion: These findings suggest that microbiomes in certain skin regions may be predisposed for interactions with Bd. Further, our results highlight the importance of considering intraindividual heterogeneities, which could provide insights relevant to predicting localized interactions with pathogens.
1 Introduction
Communities of microbes associated with multicellular organisms, also known as microbiomes, can play a significant role in the health and disease of their hosts (Cho and Blaser, 2012; Lee and Hase, 2014; Oever and Ten Netea, 2014; Robinson et al., 2010; Zilber-Rosenberg and Rosenberg, 2008). We are interested here in the skin-associated microbiome of frogs and the roles it may play in frog health. In general, the skin microbiome composition and structure in animals can influence host health by contributing to immune defenses and maintaining skin homeostasis (Sanford and Gallo, 2013). In amphibians, one key role of the skin microbiome is that it can serve as a primary defense mechanism against invading pathogens (Walke and Belden, 2016). The role of the amphibian skin microbiome in pathogen defense has become of great interest recently due to the global spread of the pathogen Batrachochytrium dendrobatidis (Bd), a chytrid fungus that causes the disease chytridiomycosis, which has led to dramatic declines and species extinctions in amphibians around the world (Fisher et al., 2009; Fisher and Garner, 2020; Scheele et al., 2019; Skerratt et al., 2007).
Bd infects keratinized epidermal cells, disrupting host osmoregulation and electrolyte balance, and often leading to mortality (Berger et al., 1998; Berger et al., 1999; Voyles et al., 2009). Interestingly, susceptibility to Bd varies widely among amphibian species, populations, and individuals (Jiménez and Sommer, 2017; Rosenblum et al., 2010). While this variation is influenced by multiple factors, including host genetics and environmental conditions, the skin-associated microbiome may also play a crucial role in determining susceptibility (Becker et al., 2015a). Studies have shown that several amphibian skin microbes can inhibit Bd, and that the structure of the skin microbiome including the number and proportion of known and/or predicted Bd-inhibitory taxa can predict the severity of infection and disease outcomes (Bates et al., 2018; Becker et al., 2015a, 2015b; Bell et al., 2018; Chen et al., 2022; Harris et al., 2009a; Jani et al., 2017; Woodhams et al., 2007b, 2015). These findings have spurred interest in probiotics for amphibians, although results have been mixed (Becker et al., 2009; Harris et al., 2009a, 2009b; Kueneman et al., 2016a; Becker et al., 2011, 2015a, 2021; Knapp et al., 2022; Woodhams et al., 2012, 2020). Probiotic effectiveness often depends on the ability of beneficial bacteria to persist on the skin, which is influenced by the existing microbial community (Becker et al., 2015a; Knapp et al., 2022; Woodhams et al., 2012). This underscores the need for a deeper understanding of skin microbiome complexity and dynamics in amphibians (Becker et al., 2011; Costello et al., 2012; Garner et al., 2016).
Bd infections in frogs are primarily limited to the ventral skin surfaces and toes (Berger et al., 1998; Berger et al., 2005; North and Alford, 2008; Pessier et al., 1999). This pattern of infection may indicate a difference in the microbial communities and niche space available in certain regions of the body, warranting an examination of the microbiomes of different body regions to understand potential regional defenses against Bd. Microbial heterogeneity across the skin could arise from the variability in epithelia among regions, with dorsal and ventral surfaces harboring differences in the types of glands and secretions produced (Berger et al., 2005; Varga et al., 2019). Additionally, it has been shown that bacterially produced compounds can act synergistically with host-produced anti-microbial peptides (AMPs) to inhibit Bd growth (Myers et al., 2012). Thus, the combination of skin architecture and bacterial composition are likely directly relevant to the distribution of Bd infection across the skin.
Evidence indicates that frog skin selects for specific microbes from the environment (Bates et al., 2018; Loudon et al., 2016; Walke et al., 2014), but whether there is selection for different microbes in body regions that are known to be preferentially infected by Bd has not been examined. In fact, few studies have assessed within-individual variability in amphibian microbiomes, although such variation has been documented in humans and other animals (Asangba et al., 2022; Bouslimani et al., 2015; Grice et al., 2009; Krog et al., 2022; Shibagaki et al., 2017; Sugden et al., 2021). Heterogeneity in microbiome structure among body regions has been detected in certain amphibian species (Bataille et al., 2016; Sabino-Pinto et al., 2016; Sanchez et al., 2017), suggesting that for at least some species, different skin regions may harbor distinct microbial communities.
In this study, we utilized high-throughput sequencing of bacterial 16S rRNA gene amplicons to characterize the skin microbiome of captive adult Sierra Nevada yellow-legged frogs (Rana sierrae) within individuals and their tank environments. This species has experienced dramatic population declines due to invasive fish and disease (Vredenburg et al., 2007; Vredenburg et al., 2010). Restoration efforts for this species often involve head-starting, where frogs are reared to adulthood in captivity before being reintroduced into the wild. Captivity is known to alter the amphibian skin microbiome, with several studies finding differences in microbiome structure and diversity between captive and wild individuals across many amphibian species, likely due to environmental and dietary differences (Antwis et al., 2014; Becker et al., 2014; Kueneman et al., 2022; Loudon et al., 2014; Sabino-Pinto et al., 2016). These captivity-induced shifts in the microbiome could impact the success of reintroduction programs, warranting closer attention to the microbiome in captivity (Redford et al., 2012). We were therefore interested in investigating their microbiomes in the captive setting. Further, we were interested in examining microbiomes in Bd-naive captive frogs so that we could determine whether frog body regions harbored differences corresponding to regions where Bd is expected to infect, rather than detecting differences that may have been driven by acute Bd infection.
By examining the skin microbiome in a captive-reared population of R. sierrae, we sought to address the following questions: (1) How do captive R. sierrae skin microbiomes differ from their tank environment microbiome? (2) How much variation is there among microbiomes of different body regions within individuals? and (3) Are there consistent differences in the skin microbiome that correspond to body regions known to be targets of Bd infection? We hypothesized that we would detect differences between frogs and their tank environments (Bataille et al., 2016; Walke et al., 2014) and among body regions (Bataille et al., 2016; Sabino-Pinto et al., 2016; Sanchez et al., 2017). Further, we hypothesized that certain microbes would be differentially abundant between body regions that tend to harbor Bd infections (ventral surfaces and feet) and body regions where infection is often absent (dorsal surfaces like the back).
2 Materials and methods
2.1 Frog population and handling
Frogs sampled for this study were members of a captive population reared to adulthood at the San Francisco Zoo from egg masses collected from a wild population in the Desolation Wilderness (El Dorado County, California; ~2,500 m elevation). This captive population was established as part of a conservation effort that was unrelated to the present study and therefore was not for research purposes alone. Adults, i.e., those with snout–vent length (SVL) ≥ 40 mm, were tagged with 8 mm unique passive integrated transponder (PIT) tags, which allow for differentiation among individuals. Groups of 8–13 frogs were housed in 52 gallon tanks (91.44 cm × 50.80 cm × 43.18 cm) filled with ~30 gallons of tap water that was filtered via reverse osmosis to remove chlorine, amines, and solids, and then supplemented with Kent Marine R/O Right, a proprietary formulation of dissolved solids and electrolytes used to restore natural water chemistry. Tanks contained biological filters to remove toxic nitrogenous compounds. We randomly selected adult frog individuals from which we collected samples and collected samples from tank environments housing those individuals.
2.2 Sample collection
We wore nitrile gloves during sample collection from frogs and surfaces in tanks using sterile synthetic fine tip dry swabs (Medical Wire & Equipment, Corsham, Wiltshire, UK; MW113). Prior to sampling, we rinsed each frog individual with 60 mL of sterile water (Culp et al., 2007; Lauer et al., 2007). From each frog, we collected a separate swab from each body region in the following order: back, outer hindlimbs, snout, vocal sack, ventral abdomen, inner forelimbs, forefeet, inner hindlimbs, hindfeet, and cloaca (Figure 1). Body regions were swabbed by taking 10 strokes to standardize sample material from regions of various sizes, taking care to only swab to the target region. The chosen target body regions were not overlapping to ensure that collected swabs were from distinct skin regions. For limbs and feet, both the left and right were sampled with the same swab. The cloaca was the smallest region sampled, for which we placed the swab under the frog vent and spun it 10 times to collect the most localized sample possible. For each frog, we also recorded the sex, tank identity, and individual identity. Next, we sampled tank environments (Supplementary Figure S1). We collected swabs of surfaces in tanks by taking 40 strokes across each surface. Surfaces sampled included (1) rock perches (affixed to the tank wall above the water surface; two samples per tank collected from the top surface of each perch), (2) underwater rocks (submerged in water; one sample per tank collected from the top surface of the rock), and (3) tank walls (three samples per tank collected from above the water surface). Tank water was sampled by filling a 60 mL syringe, passing the water through a 0.22 μm Sterivex filter (Millipore, Burlington, MA, United States), and repeating this process four times (total water filtered = 240 mL per sample; two filter samples collected per tank; Ellison et al., 2019). All samples were kept on dry ice during collection and transferred to a − 80°C freezer for storage on the same day.
Figure 1. Diagram of Rana sierrae body regions sampled in this study. Body regions sampled are encircled in dashed lines. Ventral surfaces sampled are indicated with purple text and arrows. Dorsal surfaces sampled are indicated with orange text and arrows. Feet surfaces sampled are indicated with blue text and arrows. For limbs and feet, both the right and left were sampled. Forefeet were only sampled on the ventral surface, whereas hindfeet samples were collected from both the ventral and dorsal surfaces. Frogs sampled ranged from 46.4 to 54.4 mm snout-to-vent length (SVL) and from 9.5 to 17.8 g in weight.
In this study, we analyzed samples collected from nine frogs (n = 90 microbiome swabs) and their tank environments (n = 18 microbiome swabs; n = 6 water filters) on July 28, 2015. As we randomly selected frogs to sample, they were not distributed among tanks in a balanced manner: six frogs were co-housed in one tank, two frogs were co-housed in a second tank, and one frog was housed in a third tank.
2.3 DNA extraction
We extracted DNA from microbiome swabs and water filter samples using PowerSoil DNA Isolation Kits (MoBio Laboratories, Carlsbad, CA, United States) with a modified protocol for low-biomass samples discussed with the manufacturer. These kits have been used in several previous amphibian microbiome studies (Chen et al., 2022; Ellison et al., 2019, 2021; Kueneman et al., 2014). Swabs or filters were swirled in the PowerBead tubes and left inside these tubes. Modifications to the manufacturer’s standard protocol included the following: (1) after adding Solution C1 and vortexing to mix, tubes were incubated at 65°C for 10 min; (2) tubes were then secured in a bead beater set to “homogenize,” and bead-beated for a total of 3 min (90 s on, 60 s rest, followed by 90 s on); (3) all centrifugation steps throughout were done for 1 min at 13,000 × g unless otherwise noted below; (4) we combined steps for Solutions C2 and C3 by adding 100 μL of each at once prior to 5 min incubation on ice, (5) after the C2/C3 step, we transferred 700 μL of lysate to a clean collection tube and added 700 μL of Solution C4 and 600 μL of 100% ethanol before loading on the spin filters; (6) before washing the filter with solution C5, we inserted a step to wash with 650 μL of 100% ethanol; (7) After washing with solution C5, we dried the spin column by centrifuging for 2 min at 13,000 × g; (8) We added 60 μL of Solution C6 (heated to 60°C) to the filter membrane, and we allowed this solution to sit on the filter for 5 min before centrifuging into a storage tube. Following DNA extraction, we quantified DNA concentration using a Qubit (Invitrogen, Carlsbad, CA, United States) and the dsDNA High Sensitivity Kit, and stored DNA extracts at −80°C.
2.4 Sequence generation
Sequencing libraries were prepared following the protocol “16S Metagenomic Sequencing Library Preparation” (Part # 15044223 Rev. B, Illumina, Inc., San Diego, CA, United States) with some modifications. Briefly, we PCR amplified the hypervariable V3-V4 region of the bacterial 16S rRNA gene using 341F and 805R primers (Klindworth et al., 2013) with overhang adaptors (forward primer with overhang = 5′ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG; reverse primer with overhang = 5′ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC) from each sample in triplicate, using 4 uL DNA extract per reaction. We pooled PCR products from each sample (75 μL pool) and purified them using magnetic beads (AxyPrep Mag PCR Clean-Up Kit, Axygen, Union City, CA, United States) using 60 μL beads per pool (for a 0.8X ratio of beads to PCR product), eluting in 25 μL of 10 mM Tris pH 8.5. Next, we attached dual indices and Illumina sequencing adaptors in a second round of PCR described in the Illumina protocol. We purified and normalized 25 μL of each index PCR product using SequalPrep Normalization Plate Kits (Invitrogen, Carlsbad, CA, United States), following the manufacturer’s protocol with an extension of the binding step incubation to 2–6 h. We then pooled 10 μL of purified, normalized, and indexed PCR product per sample, and used the Zymo Clean and Concentrator Kit (Zymo Research, Irvine, CA, United States) to increase the DNA concentration of the pool following the manufacturer’s protocol (using a ratio of 5:1 of DNA Binding Buffer:PCR product, and final elution using 200 μL DNA Elution Buffer). We quantified DNA in the final pool using a Qubit (Invitrogen) and sent the pooled libraries to the UC Davis Genome Center DNA Technology Core for sequencing on an Illumina MiSeq (Illumina, Inc., San Diego, CA, United States) with v3 chemistry in 2 × 300 bp run mode.
We sequenced 20 negative control samples in addition to the biological samples. These included six controls for swab sample collection (three dry swabs and three swabs rinsed with sterile water that were placed in collection tubes at the zoo on the day of sample collections and processed in the same way as biological samples), eight blank DNA extraction kit controls (two to three preparations from each of three PowerSoil kits used to extract DNA from biological samples, for which no sample was added to the PowerBead tube but otherwise processed in the same way as biological samples), and six PCR controls (one from each time biological samples were amplified, for which no sample DNA was added to the first PCR step and subsequent processing was the same as for biological samples).
2.5 Sequence processing
To demultiplex the sequence data, we used a modified version of a custom script designed by G. Jospin (https://github.com/gjospin/scripts/blob/master/Demul_trim_prep.pl). Primers were removed using cutadapt v. 3.5 (Martin, 2011) with Python v. 3.9.10 (van Rossum and Drake, 2009), discarding reads for which primers were not present. We processed the resulting sequences using the DADA2 v.1.24.0 (Callahan et al., 2016) workflow in R v. 4.2.1 with RStudio v. 2022.07.0–548 (Posit Team, 2022; R Core Team, 2022). We trimmed forward and reverse reads at 250 base pairs, truncated at the first quality score of 2, and removed them if the expected errors were greater than 4 (this removed 32.8% of sequences). We then merged reads and inferred amplicon sequence variants (ASVs; 7.9% of sequences did not pass through these steps). Next, we identified 1.8% of merged reads to be chimeric and removed them. After chimera removal, samples had a mean read depth of 23,832 with a range of 294 to 98,808 reads. We assigned taxonomy to genus level using the Ribosomal Database Project (RDP) naive Bayesian classifier algorithm and the SILVA high quality ribosomal RNA database v. 132, and species level assignments were made based on exact matching of ASVs to reference strains in the SILVA database (Quast et al., 2012; Wang et al., 2007; Yilmaz et al., 2014). We then assigned unique names to ASVs, beginning with “SV” (sequence variant) followed by a number (e.g., SV1, SV2, etc.). We removed ASVs based on taxonomic classifications that were (1) non-bacterial at the domain level (including Eukaryota, Archaea, and those unclassified to domain), (2) chloroplasts, and (3) mitochondria, which resulted in 2,373 unique ASVs in the dataset.
We used Decontam v. 1.16.0 to identify putative contaminants, implementing the prevalence method with a probability threshold of 0.5 (which identifies sequences that have a higher prevalence in negative controls than in biological samples) and setting the batch argument so that contaminants were identified independently within groups of samples associated with specific negative controls (Davis et al., 2018). We identified contaminants separately for each of the following four control sample groups: (1) dry swabs (n = 3) setting batch by the sample material so as to identify contaminants in swab samples and not filters, (2) swabs rinsed with sterile water (n = 3) setting batch by whether sampling involved rinsing with sterile water so as to identify contaminants from frogs that were rinsed prior to swabbing, (3) extraction kit blanks (n = 8) setting batch by the PowerSoil kit used so as to identify contaminants from each kit separately, and (4) PCR negative controls (n = 6) without specifying batch so as to identify contaminants associated with PCR across all samples. We then compiled a list of putative contaminants identified using each control group (102 unique ASVs) and removed them, leaving 2,271 unique ASVs in the dataset.
There has been an ongoing debate in the literature regarding the validity of rarefying read counts as a sample normalization technique for microbiome data (Cameron et al., 2021; Gloor et al., 2017; McKnight et al., 2019; McMurdie and Holmes, 2014; Weiss et al., 2017). We chose to implement this method for many analyses because we were interested in community level comparisons that can become distorted using other normalization methods (McKnight et al., 2019). Additionally, rarefying was shown to be more effective than other methods at controlling effects of sample library size when sample depths are very uneven (Schloss, 2023; Weiss et al., 2017), which was the case for our dataset (read counts ranged from 2,727 to 93,792 for biological samples). We confirmed that most sample rarefaction curves plateaued at the minimum frog sample depth (i.e., 2,727 reads; Supplementary Figure S1). For all subsequent sequence analyses except for DESeq2 differential abundance testing (see below), samples were rarefied at an even sampling depth of 2,727. As negative control samples had fewer reads, the process of rarefying removed them from the dataset.
After rarefying the dataset, we aligned remaining sequences using DECIPHER v. 2.22.0 (Wright, 2016) and built a maximum likelihood tree with a GTR + Γ(4) + I model using phangorn v. 2.8.1 (Schliep et al., 2017; Schliep, 2011) on the UC Davis Bioinformatics Core High Performance Computing Cluster in R v. 4.1.0 (R Core Team, 2022). We midpoint rooted the tree using phangorn v. 2.9.0 (Schliep et al., 2017; Schliep, 2011).
The resulting dataset analyzed for this study included 1,861 unique ASVs across 114 frog and tank environment samples.
2.6 Microbial sequence analysis and visualization
2.6.1 Alpha diversity analyses
We considered two metrics of within-sample microbial community diversity (i.e., alpha diversity): observed richness (i.e., the number of ASVs in the rarefied dataset) and Shannon diversity, which we chose to use because it incorporates both ASV richness and relative abundance (i.e., richness and evenness). We calculated these metrics using the estimate_richness function in phyloseq v. 1.40.0 (McMurdie and Holmes, 2013). Shapiro–Wilk normality tests for groups of alpha diversity estimates that we sought to compare revealed that estimates for at least one group in each comparison were not normally distributed (p < 0.05), warranting use of nonparametric statistical tests. We therefore implemented Kruskal-Wallis rank sum tests for significant differences in alpha diversity values between metadata groupings of the samples (including frog vs. environmental sample type and frog body region) using the kruskal.test function in base R v. 4.2.1 (R Core Team, 2022). For significant Kruskal-Wallis results (p ≤ 0.05), we performed post hoc Dunn tests with a Benjamini-Hochberg correction to control the false discovery rate (FDR) with multiple comparisons (dunnTest function in FSA v. 0.9.3; Ogle et al., 2022).
2.6.2 Beta diversity analyses
To assess community structure, we compared between-sample diversity (i.e., beta diversity) using three ecological distance metrics: unweighted UniFrac, weighted UniFrac, and Bray-Curtis dissimilarities. Unweighted UniFrac distance is calculated from the community phylogenetic tree as the unique fraction of branch length within a sample community that is not shared with other communities sampled (Lozupone and Knight, 2005). This metric, therefore, can be thought of as measuring distances based on community membership (presence/absence). Weighted UniFrac takes into account the relative abundances (i.e., evenness) of branch lengths in addition to membership, giving more weight to dominant organisms than rare ones (Lozupone et al., 2007). Bray-Curtis also takes into account species richness and relative abundances, but this metric is not informed by phylogeny. We used the ordinate function in phyloseq v. 1.40.0 (McMurdie and Holmes, 2013) to calculate these distances for different subsets of the data, and visualized ordinations using principal coordinate analysis (PCoA).
To test for significant effects of metadata variables on microbial community structure, we ran permutational multivariate analysis of variance (PERMANOVA) using the adonis function in vegan v. 2.6.2 with 9,999 permutations (Anderson, 2001; Oksanen et al., 2022). For each ecological distance metric, we ran a PERMANOVA on the whole dataset (frog samples and environmental samples) to test whether significant variation in microbiome structure was explained by frog and environmental sample types. We then ran sequential PERMANOVAs on the frog samples alone to test for significant effects of frog body regions after controlling for individual and/or tank effects, as well as to test for significant effects of other metadata factors like frog sex (i.e., males versus females). Next, for factors of interest that rejected the null hypothesis in these PERMANOVAs (p ≤ 0.05), we performed post hoc pairwise PERMANOVA tests to identify which levels within factors differed significantly, using the adonis.pair function in EcolUtils v. 0.1 with 9,999 permutations (Salazar, 2022). We corrected p-values for multiple comparisons using the Benjamini-Hochberg procedure.
PERMANOVA tests are sensitive to differences in group dispersion (i.e., within-group variance) and thus significant effects detected by these tests could indicate differences in the average location of groups in ordination space (i.e., group centroids), differences in group dispersion, or some combination of the two (Anderson, 2001). If PERMANOVA tests indicate a significant effect of a factor and group dispersions for that factor do not differ significantly, then we know that centroid location differs for these groups. However, if the group dispersions are significantly different, then our tests are unable to distinguish whether differences are due to dispersion alone or some combination of dispersion and centroid location. We therefore also calculated mean dispersion for factors included in post hoc pairwise PERMANOVAs and tested for significant differences using the betadisper and permutest.betadisper functions in vegan v. 2.6.2 (Oksanen et al., 2022). For significant comparisons (p ≤ 0.05), we implemented post hoc Tukey honest significant differences tests to identify which levels within factors differed significantly in their group dispersion (using the TukeyHSD function in base R, which corrects p-values for the family wise error rate in multiple comparisons).
2.6.3 Analyses of presence and relative abundance of taxa
We calculated the proportion of taxa shared between frog samples and environmental samples in the rarefied dataset. To do this, we first merged samples by frog and environment sample categories (i.e., collapsing read counts within each sample category), and removed any ASVs that had zero counts across frog and environmental samples. Next, we calculated the proportion of ASVs that were present in both frog and environment samples (“shared”). We also calculated the proportion of bacterial families that were shared between frog and environmental samples by collapsing ASVs from the same families using the tax_glom function in phyloseq v. 1.40.0 (specifying NArm = FALSE to include unclassified taxa; McMurdie and Holmes, 2013) prior to merging samples by frog and environment sample categories.
We visualized and compared the mean relative abundance of taxa for frog and environmental samples. We first transformed the rarefied sample counts to relative abundances and merged ASVs from the same families using the tax_glom function in phyloseq (specifying NArm = FALSE to include unclassified taxa). We grouped data by a metadata factor to compare frog samples to the four environmental sample types and by taxonomic ranks and then calculated mean relative abundances. We used microshades v. 1.11 (Dahl et al., 2022) to generate a stacked bar chart that simultaneously displays a higher taxonomic rank (phylum) and a lower taxonomic rank (family). We assigned colors to the five phyla that represented the highest relative abundance across samples and to the highest relative abundance families within those phyla, with remaining taxa represented in “Other” categories.
To determine whether the relative abundance of bacterial families or ASVs differed significantly between frog and environmental sample types, we implemented nonparametric Kruskal-Wallis rank sum tests followed by post hoc Dunn tests when results of Kruskal-Wallis tests were significant (p ≤ 0.05). We corrected p-values from Dunn tests using the Benjamini-Hochberg procedure.
2.6.4 DESeq2 differential relative abundance testing
We used DESeq2 v. 1.36.0 in R on filtered but un-rarified merged read counts to determine which ASVs showed significant log2 fold differences between frog body regions (Love et al., 2014). Based on previous findings that Bd infection is more prominent on ventral surfaces and toes than on the dorsal back surface of many anurans (Berger et al., 1998; Berger et al., 2005; North and Alford, 2008; Pessier et al., 1999), we chose to compare frog back samples to (1) abdomen samples, (2) inner hindlimb samples, (3) hind feet samples, and (4) forefeet samples, to examine whether bacterial community members were differentially associated with these body regions.
First, we used the phyloseq_to_deseq2 function to format the raw read count data from frog samples for DESeq2. Our ASV counts table was sparse, with only two ASVs present across all samples. Therefore, to prevent geometric means and estimated size factors for DESeq2 sample normalization from being influenced solely by these ASVs, we calculated geometric means across samples for each ASV by ignoring samples with zero counts. Then, we estimated size factors for each sample based on the geometric means (applying the estimateSizeFactors function). We filtered out low relative abundance ASVs with 10 or fewer reads total across frog samples (this filtered out 1984 unique ASVs, leaving 287 unique ASVs across the frog samples). We then ran the DESeq function on the dataset for each contrast of interest, identifying ASVs that showed significant log2 fold differences (i.e., that showed differential relative abundance; Benjamini-Hochberg corrected p-values ≤ 0.05). For ASVs that showed differential relative abundance based on DESeq2 normalized counts, we calculated the mean relative abundance across body regions from the rarefied dataset. We confirmed that a nonparametric test was appropriate (Shapiro–Wilk, p < 0.05 for multiple body regions) and implemented Kruskal-Wallis rank sum tests to determine whether mean rarefied relative abundance also differed significantly among body regions.
2.6.5 Bd-inhibitory predictions
We predicted putative Bd-inhibitory function of frog-associated microbial community members. Predictions were based on a database of full length 16S rRNA gene sequences from bacteria that were isolated and assayed for their effects on growth of Batrachochytrium pathogens (Woodhams et al., 2015). This database, which is regularly updated, has been used in several previous studies to predict anti-Bd function from amplicon data (Bletz et al., 2017; Chen et al., 2022; Jiménez et al., 2022; Kueneman et al., 2016a, 2022; Muletz Wolz et al., 2018). We used the strict inhibitory subset of the database (AmphiBac_InhibitoryStrict_2023.2; accessed from https://github.com/AmphiBac/AmphiBac-Database/) which included sequences from 2,056 inhibitory taxa. We used nucleotide BLAST v.2.9.0 to make a multiple sequence alignment with ASV sequences as queries and the inhibitory database as subject sequences, specifying a minimum e-value threshold of 1e−10 (Camacho et al., 2009). We documented the top hit result for cases where query coverage was 100% (i.e., cases where the alignment included the entire query ASV sequence) and the percent identity was ≥ 99% (i.e., cases where the query ASV shared ≥ 99% sequence similarity with an inhibitory database taxon). We defined these ASVs as putatively Bd-inhibitory. We also defined a strict subset of putative Bd-inhibitory taxa which included only ASVs that shared 100% sequence similarity with inhibitory database taxa.
We calculated putative Bd-inhibitory relative abundance in frog samples based on our strict subset of taxa (i.e., those sharing 100% sequence similarity to known inhibitory taxa) by taking the sum of read counts for the identified ASVs and dividing by the total number of reads (i.e., 2,727 reads in the rarefied dataset). To test for significant differences among frog body regions, we implemented nonparametric Kruskal-Wallis rank sum tests after finding that values were not normally distributed for multiple body regions (Shapiro–Wilk, p < 0.05). We then conducted post hoc Dunn tests, correcting p-values using the Benjamini-Hochberg procedure.
3 Results
3.1 Alpha diversity
3.1.1 Frogs vs. tank environment microbiome diversity
Within-sample diversity was significantly lower in frog samples than in tank environment sample types (including rock perch, tank wall, tank water, and underwater rock samples) based on both the observed richness (Kruskal-Wallis, chi-squared = 55.62, df = 4, p < 0.001; Dunn tests, p < 0.05; Supplementary Table S1; Figure 2A) and Shannon diversity (Kruskal-Wallis, chi-squared = 56.77, df = 4, p < 0.001; Dunn tests, p < 0.05; Supplementary Table S1; Figure 2B).
Figure 2. Alpha diversity-based comparisons of sample groupings. Within-sample diversity in terms of observed richness (A,C) and Shannon diversity (B,D). (A,B) Frog and tank environment samples with boxplots and points colored by the type of sample substrate. Panel background shading differentiates frog samples from tank environment samples. (C,D) Frog samples with boxplots and points colored by frog body region. Panel background shading differentiates groups of body regions sampled (dorsal, ventral, and feet). (A–D) Results of Kruskal-Wallis tests are shown. Post hoc Dunn test results are displayed as significance bars where applicable (*p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001).
3.1.2 Variation in frog microbiome diversity by body region
We found that Shannon diversity differed significantly by frog body region, (Kruskal-Wallis, chi-squared = 21.54, df = 9, p = 0.01; Figure 2D), but that observed richness did not (Kruskal-Wallis, chi-squared = 8.67, df = 9, p = 0.47; Figure 2C). Frog forefeet harbored higher Shannon diversity than the abdomen, back, and hind-feet in post hoc comparisons (Dunn tests, p < 0.05; Supplementary Table S2). In other words, while all frog body regions harbored a similar number of microbial community members, the forefeet harbored communities with higher evenness than certain other regions.
3.2 Beta diversity
3.2.1 Frog vs. tank environment microbiome structure
Microbial community structure was significantly different between frog samples and tank environment sample types based on all three ecological distance metrics examined (PERMANOVA, p < 0.001; Figures 3A–C; Supplementary Table S3). Post hoc pairwise PERMANOVA tests revealed that all groups (frog, tank water, tank wall, rock perch, and underwater rock) were significantly different from each other based on all three metrics (PERMANOVA, p < 0.05; Supplementary Table S4). However, the relative importance of community characteristics measured by each metric differed. We found that differences between microbial assemblages on frogs and those from environmental samples explained the highest amount of variation in weighted UniFrac (69%), followed by Bray-Curtis (45%), and finally unweighted UniFrac (29%) (Supplementary Table S3).
Figure 3. Beta diversity-based comparisons of sample groupings. (A–C) Microbial community structure of frog and tank environment samples with points colored and shaped by the type of sample substrate. (D–F) Microbial community structure of frog samples with points colored and shaped by body region. (A,D) PCoA visualizations of unweighted UniFrac distances. (B,E) PCoA visualizations of weighted UniFrac distances. (C,F) PCoA visualizations of Bray-Curtis dissimilarities.
We found group dispersion between frogs and environment sample types did not differ significantly for unweighted UniFrac (betadisper permutest, p = 0.16; Supplementary Table S5), indicating that differences based on phylogenetically informed community membership were due to differences in mean centroid locations of groups in ordination space. There were significant differences in group dispersion, however, for weighted UniFrac and Bray-Curtis (betadisper permutest, p < 0.001; Supplementary Table S5). For weighted UniFrac, pairwise comparisons of group dispersion (for frog, rock perch, tank wall, tank water, and underwater rock sample groupings) revealed that five out of ten comparisons were significantly different, which included two out of the four comparisons with frog samples (TukeyHSD, p < 0.05; Supplementary Table S6). For Bray–Curtis dissimilarity, three out of ten pairwise comparisons of group dispersion were significantly different, and these represented three of the four comparisons with frog samples (TukeyHSD, p < 0.05; Supplementary Tables S6, S7). However, ordination visualizations (Figures 3A–C) showed that frog samples clustered separately from environmental sample types for all three measures of community structure, which is evidence that in cases where differences in dispersion were detected between frog and environment samples, both the group dispersions and centroid locations may have been distinct.
3.2.2 Drivers of frog microbiome structure
Our sequential PERMANOVA models for microbial community structure within frog samples revealed that body region explained the same amount of variation after accounting for tank identity and individual identity as it did when included as the first term (Supplementary Table S7). The interaction terms for these factors were not significant and were dropped from the models. Body region explained the highest amount of variation in Bray-Curtis (32%; p < 0.001), followed by weighted UniFrac (28%; p < 0.001), and explained the least amount of variation in unweighted UniFrac distances (11%; p < 0.05). We also tested whether frog sex (i.e., males versus females) explained variation in the microbiome. However, when this factor was included after individual identity in sequential PERMANOVA models, it was not significant; thus, sex was dropped from the model.
3.2.3 Variation in frog microbiome structure among body regions
We used pairwise PERMANOVA tests to determine which frog body regions drove the variation in community structure explained by this factor. For unweighted UniFrac, none of the 45 pairwise comparisons of body regions were significantly different after correcting p-values for multiple comparisons (PERMANOVA, p > 0.05; Supplementary Table S8; Figure 3D). For weighted UniFrac, nine pairwise comparisons of body regions were significantly different (PERMANOVA, p < 0.05; Supplementary Table S8; Figure 3E), and these included seven significant differences between the back and other body regions, and two significant differences between the snout and other body regions (hindfeet and cloaca were significantly different from both the back and the snout; the abdomen, forefeet, inner forelimbs, inner hindlimbs, and outer hindlimbs were all significantly different from the back). For Bray Curtis, 13 pairwise comparisons of body regions were significantly different (PERMANOVA, p < 0.05; Supplementary Table S8; Figure 3F), including the same seven significant differences with the back and two significant differences with the snout that were identified for weighted UniFrac, as well as four additional significant differences with the snout only identified for Bray-Curtis (the snout also differed from the abdomen, forefeet, inner hindlimbs and outer hindlimbs for Bray-Curtis). The vocal sack was the only region that never differed significantly from other body regions in terms of community structure.
Group dispersion by body region did not differ significantly for any distance metrics (betadisper permutest, p > 0.05; Supplementary Table S9). These results indicate that significant results from pairwise PERMANOVA tests comparing frog body regions represented significant differences in group centroids for community structure and not differences in dispersion among body regions.
3.3 Presence and relative abundance of frog vs. tank environment microbial taxa
To investigate the similarity in presence of taxa between frog samples and environmental samples, we calculated the proportion of taxa shared between frog and environment. We found that 15.4% of ASVs were shared between frog and environmental samples (212 out of 1,378 ASVs). We also looked at the proportion of families shared between frog and environmental samples, which was 38.4% (91 out of 237 bacterial families).
To better understand the composition of microbiomes defining different types of samples, we visualized taxonomic families across all samples. This further revealed that the distribution of taxa was distinct between frog samples and the four types of environmental samples (Figure 4). Frog-associated communities were dominated by the families Burkholderiaceae (phylum Proteobacteria; mean relative abundance of 48.0 ± 2.6% SE across frogs), Rubritaleaceae (phylum Verrucomicrobia; mean relative abundance of 39.5 ± 2.5% SE across frogs), and to a lesser extent the family Pseudomonadaceae (mean relative abundance of 7.3 ± 1.2% SE across frogs) and other families in phylum Proteobacteria. The environment-associated communities were, for the most part, made up of lower relative abundances (mean relative abundance < 17%) of an increased number of families representing several phyla in addition to the Proteobacteria and Verrucomicrobia, including the Acidobacteria, Actinobacteria, Bacteroidetes, Deinococcus-Thermus, and Firmicutes (Figure 4).
Figure 4. Mean relative abundances of bacterial phyla and families. Stacked bar chart for frog, rock perch, tank wall, tank water, and underwater rock samples was generated using microshades v. 1.11 (Dahl et al., 2022). Taxa are colored by phylum, and families within each phylum are colored with unique shades of the associated phylum color. Phyla and families are ordered so that higher mean relative abundance groups appear lower in stacked bars. Blastocatellaceae belong to the phylum Acidobacteria. Clostridiaceae and Bacillaceae belong to the phylum Firmicutes.
Next, we examined whether there were significant differences in relative abundances of taxa of interest between frogs and environmental sample types. Relative abundance of the family Burkholderiaceae was significantly higher on frogs than on the four environment sample types (Kruskal-Wallis, chi-squared = 357.31, df = 4, p < 0.001; Dunn test, p < 0.001; Supplementary Table S10; Figure 4). Within the Burkholderiaceae, one ASV (“SV2,” which was unclassified at the genus level) was dominant on frogs, making up 98.8% of Burkholderiaceae rarefied read counts across frog samples and 31.42% of Burkholderiaceae rarefied read counts across environmental samples. While SV2 was present in 100% of samples, its relative abundance was significantly different between frogs and environmental sample types (Kruskal-Wallis, chi-squared = 56.49, df = 4, p < 0.001), with significantly higher relative abundance on frogs than in the environmental sample types (Dunn tests, p < 0.01; Supplementary Table S11).
The family Rubritaleaceae was also dominated by one ASV present in 100% of samples (“SV1,” which was unclassified at the genus level) that made up 99.9% of Rubritaleaceae rarefied read counts across frog samples and 92.68% of Rubritaleaceae rarefied read counts across environmental samples. While the relative abundance of the family Rubritaleaceae was not significantly different between frogs and environment sample types (Kruskal-Wallis, chi-squared = 3.66, df = 4, p = 0.45; Figure 4), SV1 relative abundance was significantly different between these groups (Kruskal-Wallis, chi-squared = 45.62, df = 4, p < 0.001), with higher relative abundance on frogs than in environmental samples (Dunn tests, p < 0.01; Supplementary Table S11).
3.4 DESeq2 differential relative abundance testing among frog body regions
To identify ASVs that were differentially abundant between body regions known to experience higher Bd infection and pathogenesis (e.g., ventral surfaces and toes) and body regions known to have markedly lower Bd infection (e.g., dorsal surfaces, mainly the back) we implemented DESeq2 analysis on raw read count data. We individually compared DESeq2 normalized counts from the abdomen, the inner hindlimbs, the hindfeet, and the forefeet to the back. The analysis identified one unique ASV, an undescribed member of the family Burkholderiaceae, that showed significant log2 fold higher normalized counts on the abdomen compared to the back (SV56; estimated log2 fold difference of 24.18; Supplementary Table S12). The mean relative abundance of this ASV in the rarefied dataset, however, was not significantly different among body regions (Kruskal Wallis, chi-squared = 5.56, df = 9, p = 0.78; Supplementary Figure S3).
3.5 Bd-inhibitory predictions among frog body regions
We identified 73 frog-associated ASVs that were putatively Bd-inhibitory, based on sharing ≥ 99% sequence identity with taxa known to inhibit Bd growth (Woodhams et al., 2015). Assigned taxonomy for all 73 ASVs is shown in Supplementary Table S13. Our strict subset of putative Bd-inhibitory taxa included 29 of these ASVs that shared 100% sequence identity with inhibitory database taxa (Supplementary Table S13). The dominant undescribed Burkholderiaceae across frog samples (SV2), shared 100% sequence identity with a known Bd-inhibitory taxon, AmphiBac_1576 (Woodhams et al., 2015). Another undescribed Burkholderiaceae that showed significant log2 fold higher normalized relative abundance on the R. sierrae abdomen than on their back (SV56; see above), shared 99.53% sequence identity with the same inhibitory taxon.
Total putative Bd-inhibitory relative abundance for our strict subset of ASVs (i.e., 29 ASVs that had identical sequences to inhibitory database taxa) differed significantly by frog body region (Kruskal-Wallis, chi-squared = 34.19, df = 9, p < 0.001; Figure 5). The back harbored significantly lower relative abundance of putative Bd-inhibitory taxa than seven other body regions (the abdomen, cloaca, forefeet, hindfeet, inner forelimbs, inner hindlimbs, and outer hindlimbs; Dunn tests, p < 0.05; Supplementary Table S14) and the snout harbored significantly lower relative abundance of putative Bd-inhibitory taxa than four other body regions (the abdomen, cloaca, hindfeet and inner hindlimbs; Dunn tests, p < 0.05; Supplementary Table S14). We confirmed that results were similar for ASVs sharing ≥ 99% similarity with inhibitory database taxa.
Figure 5. Relative abundances of putative Bd-inhibitory taxa among frog body regions. Values were calculated based on our strict subset of 29 frog-associated ASVs that shared 100% sequence identity with taxa known to inhibit Bd growth (Woodhams et al., 2015). Boxplots and points are colored by frog body region. Panel background shading differentiates groups of body regions sampled (dorsal, ventral, and feet). Results of Kruskal-Wallis tests are shown. Post hoc Dunn test results are displayed as significance bars where applicable (*p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001).
4 Discussion
Our fine-scale analysis of the skin microbiome identified characteristics that vary between captive frogs and their tank environments, and among body regions within individuals. We found that frogs harbored distinct microbial communities compared to their local tank environment. In addition, there were detectable differences between microbiomes of body regions known to be targets of Bd infection compared to those regions that infrequently experience infection. Further, elevated relative abundances of putatively Bd-inhibitory microbes were localized in body regions where we would expect interactions with Bd to occur. Together, these results help elucidate the captive microbiome of the endangered Sierra Nevada yellow-legged frog, R. sierrae, and provide a basis for predicting microbiome-pathogen interactions.
4.1 Frog skin microbiomes were distinct from their tank environment microbiome and were dominated by fewer organisms
We hypothesized that frog skin microbiomes would be distinct from their surrounding tank environment microbiome, which was supported by our results for within-sample community diversity (i.e., alpha diversity), community structure (i.e., beta diversity), and presence and relative abundance of microbial community members. We found that tank substrates and water harbored significantly higher within-sample diversity than frogs (Figures 2A,B), which agrees with a previous study of wild R. sierrae showing that lake water communities had higher observed richness than frog associated communities (Ellison et al., 2019). However, this finding differed from previous evidence that lake water microbiomes had reduced or equal diversity compared to microbial communities associated with several other species of post-metamorphic amphibians (Bates et al., 2018; Kueneman et al., 2014). The discrepancy between our results here and those of these prior studies has multiple possible explanations. One possibility is that lake water collected by Ellison et al. (2019) and tank water collected here were unusually diverse compared to other environments. Another possibility (and we note, both could be occurring) is that R. sierrae may harbor lower diversity microbiomes than other species. Reduced diversity is linked to clinical signs of chytridiomycosis (Becker and Harris, 2010) while higher community richness has been shown to correlate with host persistence after Bd invasion (Jani et al., 2017). Therefore, if R. sierrae microbiomes harbor lower diversity and richness than other species, this could relate to their high susceptibility to Bd (Vredenburg et al., 2010). Additional studies that directly compare the diversity of R. sierrae microbiomes to other species would be useful here.
We also found that community structure significantly differed between frog-associated and environment-associated microbiomes (Figures 3A–C), which has been previously reported in studies of both captive and wild amphibian populations (Albecker et al., 2019; Bates et al., 2018; Fitzpatrick and Allison, 2014; Jani et al., 2017; Jani and Briggs, 2014; Kueneman et al., 2014; Walke et al., 2014). We note that the ordinations comparing frog and environmental samples showed a characteristic “horseshoe effect,” which could be a consequence of saturation of the distance metrics (Morton et al., 2017). Since distance metrics cannot discriminate between samples that do not share ASVs, this saturation can occur when the dataset has sparse counts of ASVs across samples creating a “band table” (Morton et al., 2017). This effect does not alter our interpretation of the ordinations, and it actually provides further evidence of the high dissimilarity between environment and frog samples. There are several factors that make amphibian skin a unique and complex environment that could lead to such differences. Mucosal secretions, anti-microbial peptides (AMPs), and other secretions produced by the host regulate microbial presence and abundance on the skin, as do other microbes and the anti-microbial metabolites they produce, all of which vary between and within host species (Lillywhite and Licht, 1975; Myers et al., 2012; Tennessen et al., 2009; Walke et al., 2014; Woodhams et al., 2006a, 2006b, 2007a, 2010). By affecting which microbes can exist and persist on the skin, these interacting skin components act as a filter for microbes from the environment.
Previous studies found variable proportions of taxa shared between amphibian-and environment-associated communities, and usually dominant microorganisms on amphibians were different from those in environmental assemblages (Bates et al., 2018; Kueneman et al., 2014; Walke et al., 2014). Further, abundant microorganisms on amphibians have been shown to be rare in the environment (Bates et al., 2018; Kueneman et al., 2014; Walke et al., 2014). Here, we found that only 15.4% of ASVs were shared between frogs and environmental samples. We also looked at the proportion of shared bacterial families between frog and environmental samples, which was higher at 38.4%. We note that the proportion of taxa shared with their environment may be lower for captive frogs than their wild counterparts, as was shown previously (Bataille et al., 2016).
Additionally, in our study, a major difference in the distribution of taxa between frogs and their tank environments was that the two sequence variants found to be dominant on frogs (SV1 in the family Rubritaleaceae and SV2 in the family Burkholderiaceae) both showed significantly lower relative abundance on tank and perch substrates and in tank water. Similarly, abundance-weighted metrics for microbial community structure explained a higher proportion of variation between frogs and environmental sample types than did an unweighted metric, supporting that relative abundances of community members better define differences than community membership alone. These results may suggest that high relative abundances of certain bacteria, including the two dominant taxa on frogs, were selected for by their skin (Loudon et al., 2016; Walke et al., 2014). A caveat of these statistical comparisons is that the data used is compositional (i.e., relative abundances must sum to 100%). Therefore, care must be taken with the interpretation of differences in relative abundances across samples since they do not represent absolute abundances and are standardized to the rarefied read count. For example, the higher relative abundances of these two taxa on frogs than in their environment could have resulted from higher absolute abundances on frogs, but it also could have resulted from reduced abundances of other taxa on frogs that inflated the relative abundance of these two ASVs. Regardless, it is interesting that microbial relative abundances on R. sierrae skin were dominated by only two sequence variants, and this result was consistent with previous studies of both wild and captive amphibians that reported dominance by one or few bacterial strains (Bates et al., 2018; Kueneman et al., 2014, 2016b; Loudon et al., 2014).
The Rubritaleaceae are a little studied family of Gram-negative bacteria in the phylum Verrucomicrobia, containing only five described species isolated from marine animals or marine sediment (Kasai et al., 2007; Rosenberg, 2014; Scheuermayer et al., 2006; Yoon et al., 2007; Yoon et al., 2008). The 16S rRNA genes of these species are very highly conserved and the species are not distinguishable by 16S amplicon analysis (Rosenberg, 2014). This may explain why we were unable to assign taxonomy below the family level for the dominant Rubritaleaceae sequence variant on frogs. Described Rubritaleaceae species are non-motile, obligate aerobes that synthesize carotenoid pigments, resulting in red-colored colonies (Rosenberg, 2014). The production of carotenoid pigments by Rubritaleaceae on frog skin may affect the skin-associated microbial community, as previous studies have shown that dietary carotenoid intake by amphibians increased community richness and shifted community structure of frog skin microbiomes (Antwis et al., 2014; Edwards et al., 2017).
The only previous mention of the Rubritaleaceae in amphibian microbiomes was from a study of captive Rana muscosa, the sister species to R. sierrae, conducted in the same facility at the San Francisco Zoo as our study (Jani et al., 2021). However, the phylum Verrucomicrobia, which includes Rubritaleaceae, has been detected on amphibian skin in several studies based on 16S rRNA gene data (Becker et al., 2014; Belden et al., 2015; Kueneman et al., 2014, 2016b; Longo et al., 2015; Loudon et al., 2014; Sabino-Pinto et al., 2016; Sanchez et al., 2017). Verrucomicrobia is a widely distributed phylum that has been found in various environments including soil, marine and freshwater, and animal intestines (Bergmann et al., 2011; Hugenholtz et al., 1998; Parveen et al., 2013; van Passel et al., 2011a; van Passel et al., 2011b). Although a previous study found that Verrucomicrobia were higher in relative abundance on wild than captive Panamanian golden frogs (Atelopus zeteki; Becker et al., 2014), we hypothesize that the high relative abundance of Rubritaleaceae and Verrucomicrobia observed in the present study may be unique to captivity. This owes to the fact that Verrucomicrobia, though present, were not high in relative abundance in previous studies of wild R. sierrae populations (Jani and Briggs, 2014), even for populations from the same site used to source the San Francisco Zoo population examined here using the same primers for 16S rRNA gene amplification (Ellison et al., 2019, 2021). It is possible that in captivity, these taxa replace other taxa with similar functional abilities on the skin in the wild, but this requires further investigation.
The other dominant amphibian associated sequence variant was a member of the family Burkholderiaceae. This family consists of ecologically, phenotypically, and metabolically diverse Gram-negative bacteria found in soil, water, and in association with plants, animals, and fungi (Coenye, 2014). While some Burkholderiaceae are pathogens to plants and animals including humans (Coenye, 2014), others have been shown to suppress fungal pathogens (Carrión et al., 2018). The dominant Burkholderiaceae sequence variant observed here shared 100% sequence identity with a bacterial isolate from the Bd inhibitory database, suggesting that this bacterium may help suppress Bd proliferation on the skin (AmphiBac_1576/Ranamuscosa-inhibitory_37; Supplementary Table S13; Woodhams et al., 2015).
The order Burkholderiales, which includes the family Burkholderiaceae, has been identified as highly relatively abundant on amphibians in several studies (Bataille et al., 2016; Bates et al., 2018; Kueneman et al., 2014), including studies of R. sierrae (Ellison et al., 2019, 2021). Previous research on wild R. sierrae found Burkholderiaceae on their skin, but it showed lower relative abundance compared to another family in the order, Comamonadaceae (Ellison et al., 2019). Notably, several amphibian microbiome studies report that a single Comamonadaceae sequence variant dominated the community in much the same way as the dominant Burkholderiaceae sequence variant did here (Bates et al., 2018; Kueneman et al., 2014, 2016b). Interestingly, while the taxonomic assignment for this bacterium was to Burkholderiaceae, the Bd inhibitory bacterial isolate with identical amplicon sequence was classified as an undescribed Comamonadaceae in the inhibitory database metadata (Woodhams et al., 2015). This discrepancy illustrates how choice of assignment algorithm and/or taxonomic database can impact such classifications. Therefore, the dominant sequence variant we identified as a Burkholderiaceae is likely to be closely related to the dominant Comamonadaceae found in other amphibian studies, or it may even represent the same bacterium. Despite amplicon sequence similarity to a known Bd-inhibitor, further work is needed to determine whether the specific taxon identified here exhibits anti-Bd function.
4.2 Frog body regions showed spatial variation in the microbiome, which corresponded to expected spatial variation in Bd infection
We tested our hypotheses that we would detect differences in microbiomes among body regions within frog individuals and that some variation would correspond to body regions where Bd infection would be expected to occur. We examined uninfected frogs (see Supplementary material) in order to determine whether they natively harbored such differences, rather than describing Bd-driven impacts to microbiomes. While our findings for within-sample diversity provide some support for differences among body regions, findings for community structure and relative abundance of putative Bd-inhibitory taxa showed stronger support for heterogeneity among body regions and that differences were associated with predicted Bd localization, discussed below.
Within-sample diversity, for the most part, did not differ between frog body regions; the only exceptions were significantly higher relative Shannon diversity on frog forefeet compared to the abdomen, hindfeet, and back (Figures 2C,D). This differs from previous studies of the Bombina orientalis microbiome that found ventral surfaces harbored higher richness and diversity than dorsal surfaces (Bataille et al., 2016; Sabino-Pinto et al., 2016) but is consistent with findings from other amphibian species that showed no such differences (for Bufo japonicus, Cynops pyrrhogaster, Odorrana splendida, and Rana japonica; Sabino-Pinto et al., 2016). Interestingly, there did not appear to be a relationship between the size of each body part and microbiome diversity. The forefeet were one of the smallest body regions sampled, but harbored higher diversity than larger body regions like the back. This suggests that standardizing the number of strokes of each body region was sufficient to control for differences in the area occupied by each body region.
Microbial community structure differed primarily between the back and other body regions. Previous studies also found significant differences in microbiome structure across the skin for two different amphibian species (Bataille et al., 2016; Sabino-Pinto et al., 2016; Sanchez et al., 2017). Though in one study, community structure only differed significantly between dorsal and ventral surfaces in captivity and not in the wild, warranting future investigation of within-individual microbiome heterogeneity of wild R. sierrae (Bataille et al., 2016). Here, while we did not detect differences in unweighted community membership among body regions, we found that the back differed significantly from seven other body regions based on two metrics of relative abundance weighted community structure (Figures 3D–F). This suggests that shifts in relative abundances of community members are more important to differences between body regions than presence/absence of organisms. Supporting this claim, we found that the relative abundance of putative Bd-inhibitory taxa differed significantly in the same comparisons of the back with other body regions (Figure 5; Woodhams et al., 2015). These taxa had significantly higher relative abundance on body regions including the abdomen, inner hind-limbs, and feet compared to the back. In addition, one putatively Bd-inhibitory undescribed member of the Burkholderiaceae showed significant log2 fold higher normalized read counts on the abdomen compared to the back (SV56; Supplementary Tables S12, S13; Woodhams et al., 2015).
Our finding that much of what defines heterogeneity in the skin microbiome are differences between the back and other body regions supports our hypothesis that microbiome variation corresponds to spatial heterogeneity in Bd infection across the skin. Studies have shown that Bd infection occurs most on ventral surfaces, hindfeet, and toes, and is either absent or minimal (i.e., very few Bd sporangia) on dorsal surfaces like the back (Berger et al., 1998; Berger et al., 2005; North and Alford, 2008; Pessier et al., 1999). In addition, it has been shown that the back of an amphibian experiences fewer pathological changes due to Bd infection than do other body surfaces (Berger et al., 2005). The fact that putatively Bd-inhibitory taxa showed higher relative abundance on ventral surfaces and feet compared to the back suggests that they may directly interact with Bd upon infection. Further, they may help to control infection on body regions where Bd is known to infect. Indeed, previous work has shown that higher putative Bd-inhibitory relative abundance correlates with lower Bd loads (Chen et al., 2022). However, isolation and functional characterization of these R. sierrae associated taxa are needed to determine if they would act to inhibit Bd growth in practice. Additionally, as discussed above, the data is compositional, and therefore quantitative analyses of taxa are needed to determine whether differences in relative abundance of putatively inhibitory taxa were driven by differences in their absolute abundances or by differences in abundances of other taxa.
The reasons that both Bd infections and microbiomes differ among body regions may relate to differences in skin architecture of the amphibian host. For example, there are usually larger and more numerous granular glands (also referred to as serous glands) on the back compared to ventral surfaces (Berger et al., 2005; Varga et al., 2019). Granular glands secrete bioactive molecules that assist in host defense, including AMPs (Varga et al., 2019). Such differences in the skin landscape may contribute to lower Bd infection on the back and to the differences in the skin microbiome between the back and other body regions that we observed here.
Our results suggest that where you collect an amphibian skin swab from (i.e., which body regions) will affect the resultant community observed. However, we emphasize that this may not apply to all types of amphibians. The heterogeneity in skin structure, microbiome structure, and Bd localization across the skin are all likely related to the evolved ecology of the amphibian. R. sierrae is a semi-aquatic species that spends much of its time basking at the edges of lakes and streams, keeping ventral surfaces and toes more moist than the back. Bd zoospores require water to disperse, so differences in moisture across the skin due to an amphibian’s lifestyle and ecology may contribute to spatial heterogeneity of Bd infection across the skin. For example, in a fully aquatic amphibian species (Xenopus tropicalus), no differences were detected in Bd infection between dorsal and ventral regions (Parker et al., 2002). Future studies would benefit from comparing differences in the microbiome structure, skin architecture, and spatial heterogeneity in Bd infection across amphibians of differing ecologies to determine whether there are consistent patterns and to elucidate the role that ecology has played in the evolution of such differences.
5 Conclusion
We found that captive R. sierrae microbiomes were distinct from microbiomes of their local environment and that there was variation in microbiomes among different regions of their skin. Some differences across the skin aligned with regions where Bd is known to infect, despite the fact that sampled frogs were naive to Bd. We suggest that understanding microbiome variation within captive individuals could be important to interventions designed to protect amphibians threatened by Bd, which are often applied in the captive setting. For example, differences across the skin could be exploited to focus on altering microbiomes of ventral surfaces and feet that gain higher Bd loads, perhaps by boosting abundances of Bd-inhibitory microbes that occur natively in these skin regions. Additionally, by elucidating microbiome variation within individuals, we can better understand and develop models to predict corresponding variation in Bd intensity. Future studies that examine variation within wild individuals would also be useful for this purpose. Detected natural variation may relate to how susceptible frogs will be to high levels of infection (Ellison et al., 2019; Jani and Briggs, 2014), which could also be an indicator of how healthy and resilient frogs are in the face of other pathogens or environmental stressors.
Data availability statement
The raw 16S rRNA gene amplicon sequence data for this project has been deposited in the National Center for Biotechnology Information Sequence Read Archive under BioProject PRJNA1219149.
Ethics statement
The animal study was approved by the UC Davis Institutional Animal Care and Use Committee (Protocol #18732) and the San Francisco Zoo Research Review Committee. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
SG: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Visualization, Writing – original draft, Writing – review & editing. JE: Conceptualization, Funding acquisition, Resources, Supervision, Writing – review & editing, Methodology.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by grants from the UC Davis Center for Population Biology awarded to SG and the Alfred P. Sloan Foundation to JE. SG was supported in part by the NIH Animal Models of Infectious Diseases Training Program T32 AI060555 Ruth L. Kirschstein National Research Service Award to SG. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Acknowledgments
We thank Jessie Bushell and the San Francisco Zoo for permitting and facilitating sampling. We thank Roland A. Knapp and the Sierra Nevada Aquatic Research Laboratory for assaying our swab samples for Bd. We thank Vance T. Vredenburg for providing advice on experimental design. We thank Cassandra L. Ettinger and John J. Stachowicz for providing advice on analysis and on the manuscript. The work presented is derived from the doctoral dissertation of SG (Ghose, 2024) and was previously made available as a preprint on bioRxiv (Ghose and Eisen, 2025).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The authors declare that no Gen AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2025.1579231/full#supplementary-material
References
Albecker, M. A., Belden, L. K., and McCoy, M. W. (2019). Comparative analysis of anuran amphibian skin microbiomes across inland and coastal wetlands. Microb. Ecol. 78, 348–360. doi: 10.1007/s00248-018-1295-9
Anderson, M. J. (2001). A new method for non parametric multivariate analysis of variance. Austral Ecol. 26, 32–46. doi: 10.1111/j.1442-9993.2001.01070.pp.x
Antwis, R. E., Haworth, R. L., Engelmoer, D. J. P., Ogilvy, V., Fidgett, A. L., and Preziosi, R. F. (2014). Ex situ diet influences the bacterial community associated with the skin of red-eyed tree frogs (Agalychnis callidryas). PLoS One 9:e85563. doi: 10.1371/journal.pone.0085563
Asangba, A. E., Mugisha, L., Rukundo, J., Lewis, R. J., Halajian, A., Cortés-Ortiz, L., et al. (2022). Large comparative analyses of primate body site microbiomes indicate that the oral microbiome is unique among all body sites and conserved among nonhuman primates. Microbiol. Spectr. 10:e0164321. doi: 10.1128/spectrum.01643-21
Bataille, A., Lee-Cruz, L., Tripathi, B., Kim, H., and Waldman, B. (2016). Microbiome variation across amphibian skin regions: implications for chytridiomycosis mitigation efforts. Microb. Ecol. 71, 221–232. doi: 10.1007/s00248-015-0653-0
Bates, K. A., Clare, F. C., O’Hanlon, S., Bosch, J., Brookes, L., Hopkins, K., et al. (2018). Amphibian chytridiomycosis outbreak dynamics are linked with host skin bacterial community structure. Nat. Commun. 9, 693–611. doi: 10.1038/s41467-018-02967-w
Becker, M. H., Brophy, J. A. N., Barrett, K., Bronikowski, E., Evans, M., Glassey, E., et al. (2021). Genetically modifying skin microbe to produce violacein and augmenting microbiome did not defend Panamanian golden frogs from disease. ISME Commun. 1, 57–10. doi: 10.1038/s43705-021-00044-w
Becker, M. H., Brucker, R. M., Schwantes, C. R., Harris, R. N., and Minbiole, K. P. C. (2009). The bacterially produced metabolite violacein is associated with survival of amphibians infected with a lethal fungus. Appl. Environ. Microbiol. 75, 6635–6638. doi: 10.1128/AEM.01294-09
Becker, M. H., and Harris, R. N. (2010). Cutaneous bacteria of the redback salamander prevent morbidity associated with a lethal disease. PLoS One 5, 1–6. doi: 10.1371/journal.pone.0010957
Becker, M. H., Harris, R. N., Minbiole, K. P. C., Schwantes, C. R., Rollins-Smith, L. A., Reinert, L. K., et al. (2011). Towards a better understanding of the use of probiotics for preventing chytridiomycosis in Panamanian golden frogs. Eco Health 8, 501–506. doi: 10.1007/s10393-012-0743-0
Becker, M. H., Richards-Zawacki, C. L., Gratwicke, B., and Belden, L. K. (2014). The effect of captivity on the cutaneous bacterial community of the critically endangered Panamanian golden frog (Atelopus zeteki). Biol. Conserv. 176, 199–206. doi: 10.1016/j.biocon.2014.05.029
Becker, M. H., Walke, J. B., Cikanek, S., Savage, A. E., Mattheus, N., Santiago, C. N., et al. (2015a). Composition of symbiotic bacteria predicts survival in Panamanian golden frogs infected with a lethal fungus. Proc. R. Soc. B Biol. Sci. 282:20142881. doi: 10.1098/rspb.2014.2881
Becker, M. H., Walke, J. B., Murrill, L., Woodhams, D. C., Reinert, L. K., Rollins-Smith, L. A., et al. (2015b). Phylogenetic distribution of symbiotic bacteria from Panamanian amphibians that inhibit growth of the lethal fungal pathogen Batrachochytrium dendrobatidis. Mol. Ecol. 24, 1628–1641. doi: 10.1111/mec.13135
Belden, L. K., Hughey, M. C., Rebollar, E. A., Umile, T. P., Loftus, S. C., Burzynski, E. A., et al. (2015). Panamanian frog species host unique skin bacterial communities. Front. Microbiol. 6:1171. doi: 10.3389/fmicb.2015.01171
Bell, S. C., Garland, S., and Alford, R. A. (2018). Increased numbers of culturable inhibitory bacterial taxa may mitigate the effects of Batrachochytrium dendrobatidis in Australian wet tropics frogs. Front. Microbiol. 9:1604. doi: 10.3389/fmicb.2018.01604
Berger, L., Speare, R., Daszak, P., Green, D. E., Cunningham, A. A., Goggin, C. L., et al. (1998). Chytridiomycosis causes amphibian mortality associated with population declines in the rain forests of Australia and Central America. Proc. Natl. Acad. Sci. USA 95, 9031–9036. doi: 10.1073/pnas.95.15.9031
Berger, L., Speare, R., and Hyatt, A. (1999). Chytrid fungi and amphibian declines: overview, implications and future directions. Canberra: Environment Australia, 23–33.
Berger, L., Speare, R., and Skerratt, L. F. (2005). Distribution of Batrachochytrium dendrobatidis and pathology in the skin of green tree frogs Litoria caerulea with severe chytridiomycosis. Dis. Aquat. Org. 68, 65–70. doi: 10.3354/dao068065
Bergmann, G. T., Bates, S. T., Eilers, K. G., Lauber, C. L., Caporaso, J. G., Walters, W. A., et al. (2011). The under-recognized dominance of Verrucomicrobia in soil bacterial communities. Soil Biol. Biochem. 43, 1450–1455. doi: 10.1016/j.soilbio.2011.03.012
Bletz, M. C., Perl, R. B., Bobowski, B. T., Japke, L. M., Tebbe, C. C., Dohrmann, A. B., et al. (2017). Amphibian skin microbiota exhibits temporal variation in community structure but stability of predicted Bd-inhibitory function. ISME J. 11, 1521–1534. doi: 10.1038/ismej.2017.41
Bouslimani, A., Porto, C., Rath, C. M., Wang, M., Guo, Y., Gonzalez, A., et al. (2015). Molecular cartography of the human skin surface in 3D. Proc. Natl. Acad. Sci. 112, E2120–E2129. doi: 10.1073/pnas.1424409112
Callahan, B. J., Mcmurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., and Holmes, S. P. (2016). Dada2: high-resolution sample inference from illumina amplicon data. Nat. Methods 13, 581–583. doi: 10.1038/nMeth.3869
Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2009). BLAST+: architecture and applications. BMC Bioinformatics 10:421. doi: 10.1186/1471-2105-10-421
Cameron, E. S., Schmidt, P. J., Tremblay, B. J.-M., Emelko, M. B., and Müller, K. M. (2021). Enhancing diversity analysis by repeatedly rarefying next generation sequencing data describing microbial communities. Sci. Rep. 11:22302. doi: 10.1038/s41598-021-01636-1
Carrión, V. J., Cordovez, V., Tyc, O., Etalo, D. W., de Bruijn, I., de Jager, V. C. L., et al. (2018). Involvement of Burkholderiaceae and sulfurous volatiles in disease-suppressive soils. ISME J. 12, 2307–2321. doi: 10.1038/s41396-018-0186-x
Chen, M. Y., Kueneman, J. G., González, A., Humphrey, G., Knight, R., and McKenzie, V. J. (2022). Predicting fungal infection rate and severity with skin-associated microbial communities on amphibians. Mol. Ecol. 31, 2140–2156. doi: 10.1111/mec.16372
Cho, I., and Blaser, M. J. (2012). The human microbiome: at the interface of health and disease. Nat. Rev. Genet. 13, 260–270. doi: 10.1038/nrg3182
Coenye, T. (2014). “The family Burkholderiaceae” in The prokaryotes: Alphaproteobacteria and Betaproteobacteria. eds. E. Rosenberg, E. F. DeLong, S. Lory, E. Stackebrandt, and F. Thompson (Berlin, Heidelberg: Springer), 759–776.
Costello, E. K., Stagaman, K., Dethlefsen, L., Bohannan, B. J. M., and Relman, D. A. (2012). The application of ecological theory toward an understanding of the human microbiome. Science 336, 1255–1262. doi: 10.1126/science.1224203
Culp, C. E., Falkinham, J. O., and Belden, L. K. (2007). Identification of the natural bacterial microflora on the skin of eastern newts, bullfrog tadpoles and redback salamanders. Herpetologica 63, 66–71. doi: 10.1655/0018-0831(2007)63[66:IOTNBM]2.0.CO;2
Dahl, E. M., Neer, E., Bowie, K. R., Leung, E. T., and Karstens, L. (2022). microshades: An R package for improving color accessibility and organization of microbiome data. Microbiol. Resour. Announc, 11:e00795–22. doi: 10.1128/mra.00795-22
Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A., and Callahan, B. J. (2018). Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome 6:226. doi: 10.1186/s40168-018-0605-2
Edwards, C. L., Byrne, P. G., Harlow, P., and Silla, A. J. (2017). Dietary carotenoid supplementation enhances the cutaneous bacterial communities of the critically endangered southern corroboree frog (Pseudophryne corroboree). Microb. Ecol. 73, 435–444. doi: 10.1007/s00248-016-0853-2
Ellison, S., Knapp, R. A., Sparagon, W., Swei, A., and Vredenburg, V. T. (2019). Reduced skin bacterial diversity correlates with increased pathogen infection intensity in an endangered amphibian host. Mol. Ecol. 28, 127–140. doi: 10.1111/mec.14964
Ellison, S., Knapp, R., and Vredenburg, V. (2021). Longitudinal patterns in the skin microbiome of wild, individually marked frogs from the Sierra Nevada, California. ISME Commun. 1:45. doi: 10.1038/s43705-021-00047-7
Fisher, M. C., and Garner, T. W. J. (2020). Chytrid fungi and global amphibian declines. Nat. Rev. Microbiol. 18, 332–343. doi: 10.1038/s41579-020-0335-x
Fisher, M. C., Garner, T. W. J., and Walker, S. F. (2009). Global emergence of Batrachochytrium dendrobatidis and amphibian chytridiomycosis in space, time, and host. Ann. Rev. Microbiol. 63, 291–310. doi: 10.1146/annurev.micro.091208.073435
Fitzpatrick, B. M., and Allison, A. L. (2014). Similarity and differentiation between bacteria associated with skin of salamanders (Plethodon jordani) and free-living assemblages. FEMS Microbiol. Ecol. 88, 482–494. doi: 10.1111/1574-6941.12314
Garner, T. W. J., Schmidt, B. R., Martel, A., Pasmans, F., Muths, E., Cunningham, A. A., et al. (2016). Mitigating amphibian chytridiomycoses in nature. Philos. Trans. R. Soc. B Biol. Sci. 371:20160207. doi: 10.1098/rstb.2016.0207
Ghose, S. L. (2024). Microbial dynamics in amphibian conservation and disease ecology. [dissertation], Davis, CA: University of California, Davis. ProQuest ID: Ghose_ucdavis_0029D_23235. Merritt ID: ark:/13030/m5jn3mpf. Available online at: https://escholarship.org/uc/item/1nf1h072.
Ghose, S. L., and Eisen, J. A. (2025). Skin microbiomes of frogs vary among individuals and body regions, revealing differences that reflect known patterns of chytrid infection. bioRxiv [Preprint]. doi: 10.1101/2025.02.05.636728
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V., and Egozcue, J. J. (2017). Microbiome datasets are compositional: and this is not optional. Front. Microbiol. 8:2224. doi: 10.3389/fmicb.2017.02224
Grice, E. A., Kong, H. H., Conlan, S., Deming, C. B., Davis, J., Young, A. C., et al. (2009). Topographical and temporal diversity of the human skin microbiome. Science 324, 1190–1192. doi: 10.1126/science.1171700
Harris, R. N., Brucker, R. M., Walke, J. B., Becker, M. H., Schwantes, C. R., Flaherty, D. C., et al. (2009a). Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824. doi: 10.1038/ismej.2009.27
Harris, R. N., Lauer, A., Simon, M. A., Banning, J. L., and Alford, R. A. (2009b). Addition of antifungal skin bacteria to salamanders ameliorates the effects of chytridiomycosis. Dis. Aquat. Org. 83, 11–16. doi: 10.3354/dao02004
Hugenholtz, P., Goebel, B. M., and Pace, N. R. (1998). Impact of culture-independent studies on the emerging phylogenetic view of bacterial diversity. J. Bacteriol. 180, 4765–4774. doi: 10.1128/jb.180.18.4765-4774.1998
Jani, A. J., and Briggs, C. J. (2014). The pathogen Batrachochytrium dendrobatidis disturbs the frog skin microbiome during a natural epidemic and experimental infection. Proc. Natl. Acad. Sci. USA 111, E5049–E5058. doi: 10.1073/pnas.1412752111
Jani, A. J., Bushell, J., Arisdakessian, C. G., Belcaid, M., Boiano, D. M., Brown, C., et al. (2021). The amphibian microbiome exhibits poor resilience following pathogen-induced disturbance. ISME J. 15, 1628–1640. doi: 10.1038/s41396-020-00875-w
Jani, A. J., Knapp, R. A., and Briggs, C. J. (2017). Epidemic and endemic pathogen dynamics correspond to distinct host population microbiomes at a landscape scale. Proc. R. Soc. B Biol. Sci. 284:20170944. doi: 10.1098/rspb.2017.0944
Jiménez, R. R., Carfagno, A., Linhoff, L., Gratwicke, B., Woodhams, D. C., Chafran, L. S., et al. (2022). Inhibitory bacterial diversity and mucosome function differentiate susceptibility of Appalachian salamanders to chytrid fungal infection. Appl. Environ. Microbiol. 88:e0181821. doi: 10.1128/aem.01818-21
Jiménez, R. R., and Sommer, S. (2017). The amphibian microbiome: natural range of variation, pathogenic dysbiosis, and role in conservation. Biodivers. Conserv. 26, 763–786. doi: 10.1007/s10531-016-1272-x
Kasai, H., Katsuta, A., Sekiguchi, H., Matsuda, S., Adachi, K., Shindo, K., et al. (2007). Rubritalea squalenifaciens sp. nov., a squalene-producing marine bacterium belonging to subdivision 1 of the phylum ‘Verrucomicrobia.’. Int. J. Syst. Evol. Microbiol. 57, 1630–1634. doi: 10.1099/ijs.0.65010-0
Klindworth, A., Pruesse, E., Schweer, T., Peplies, J., Quast, C., Horn, M., et al. (2013). Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41:e1. doi: 10.1093/nar/gks808
Knapp, R. A., Joseph, M. B., Smith, T. C., Hegeman, E. E., Vredenburg, V. T., Erdman, J. E., et al. (2022). Effectiveness of antifungal treatments during chytridiomycosis epizootics in populations of an endangered frog. PeerJ 10:e12712. doi: 10.7717/peerj.12712
Krog, M. C., Hugerth, L. W., Fransson, E., Bashir, Z., Nyboe Andersen, A., Edfeldt, G., et al. (2022). The healthy female microbiome across body sites: effect of hormonal contraceptives and the menstrual cycle. Hum. Reprod. 37, 1525–1543. doi: 10.1093/humrep/deac094
Kueneman, J. G., Bletz, M. C., Becker, M., Gratwicke, B., Garcés, O. A., Hertz, A., et al. (2022). Effects of captivity and rewilding on amphibian skin microbiomes. Biol. Conserv. 271:109576. doi: 10.1016/j.biocon.2022.109576
Kueneman, J. G., Parfrey, L. W., Woodhams, D. C., Archer, H. M., Knight, R., and McKenzie, V. J. (2014). The amphibian skin-associated microbiome across species, space and life history stages. Mol. Ecol. 23, 1238–1250. doi: 10.1111/mec.12510
Kueneman, J. G., Woodhams, D. C., Harris, R., Archer, H. M., Knight, R., and Mckenzie, V. J. (2016a). Probiotic treatment restores protection against lethal fungal infection lost during amphibian captivity. Proc. Biol. Sci. 283:20161553. doi: 10.1098/rspb.2016.1553
Kueneman, J. G., Woodhams, D. C., Treuren, W. V., and Archer, H. M. (2016b). Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas). Int. Soc. Microb. Ecol. 10, 934–944. doi: 10.1038/ismej.2015.168
Lauer, A., Simon, M. A., Banning, J. L., André, E., Duncan, K., and Harris, R. N. (2007). Common cutaneous bacteria from the eastern red-backed salamander can inhibit pathogenic fungi. Copeia 2007, 630–640. doi: 10.1643/0045-8511(2007)2007[630:CCBFTE]2.0.CO;2
Lee, W.-J., and Hase, K. (2014). Gut microbiota-generated metabolites in animal health and disease. Nat. Chem. Biol. 10, 416–424. doi: 10.1038/nchembio.1535
Lillywhite, H. B., and Licht, P. (1975). A comparative study of integumentary mucous secretions in amphibians. Comp. Biochem. Physiol. A 51, 937–941. doi: 10.1016/0300-9629(75)90077-8
Longo, A. V., Savage, A. E., Hewson, I., and Zamudio, K. R. (2015). Seasonal and ontogenetic variation of skin microbial communities and relationships to natural disease dynamics in declining amphibians. R. Soc. Open Sci. 2:140377. doi: 10.1098/rsos.140377
Loudon, A. H., Venkataraman, A., Van Treuren, W., Woodhams, D. C., Parfrey, L. W., McKenzie, V. J., et al. (2016). Vertebrate hosts as islands: dynamics of selection, immigration, loss, persistence, and potential function of bacteria on salamander skin. Front. Microbiol. 7, 1–11. doi: 10.3389/fmicb.2016.00333
Loudon, A. H., Woodhams, D. C., Parfrey, L. W., Archer, H., Knight, R., McKenzie, V., et al. (2014). Microbial community dynamics and effect of environmental microbial reservoirs on red-backed salamanders (Plethodon cinereus). ISME J. 8, 830–840. doi: 10.1038/ismej.2013.200
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Lozupone, C. A., Hamady, M., Kelley, S. T., and Knight, R. (2007). Quantitative and qualitative B diversity measures lead to different insights into factors that structure microbial communities. Appl. Environ. Microbiol. 73, 1576–1585. doi: 10.1128/AEM.01996-06
Lozupone, C., and Knight, R. (2005). UniFrac: a new phylogenetic method for comparing microbial communities. Appl. Environ. Microbiol. 71, 8228–8235. doi: 10.1128/AEM.71.12.8228-8235.2005
Martin, M. (2011). Cutadapt removes adaptor sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. doi: 10.14806/ej.17.1.200
McKnight, D. T., Huerlimann, R., Bower, D. S., Schwarzkopf, L., Alford, R. A., and Zenger, K. R. (2019). Methods for normalizing microbiome data: an ecological perspective. Methods Ecol. Evol. 10, 389–400. doi: 10.1111/2041-210X.13115
McMurdie, P. J., and Holmes, S. (2013). Phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217. doi: 10.1371/journal.pone.0061217
McMurdie, P. J., and Holmes, S. (2014). Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Comput. Biol. 10:e1003531. doi: 10.1371/journal.pcbi.1003531
Morton, J. T., Toran, L., Edlund, A., Metcalf, J. L., Lauber, C., and Knight, R. (2017). Uncovering the horseshoe effect in microbial analyses. mSystems 2:e00166-16. doi: 10.1128/msystems.00166-16
Muletz Wolz, C. R., Yarwood, S. A., Campbell Grant, E. H., Fleischer, R. C., and Lips, K. R. (2018). Effects of host species and environment on the skin microbiome of plethodontid salamanders. J. Anim. Ecol. 87, 341–353. doi: 10.1111/1365-2656.12726
Myers, J. M., Ramsey, J. P., Blackman, A. L., Nichols, A. E., Minbiole, K. P. C., and Harris, R. N. (2012). Synergistic inhibition of the lethal fungal pathogen Batrachochytrium dendrobatidis: the combined effect of symbiotic bacterial metabolites and antimicrobial peptides of the frog Rana muscosa. J. Chem. Ecol. 38, 958–965. doi: 10.1007/s10886-012-0170-2
North, S., and Alford, R. A. (2008). Infection intensity and sampling locality affect Batrachochytrium dendrobatidis distribution among body regions on green-eyed tree frogs Litoria genimaculata. Dis. Aquat. Org. 81, 177–188. doi: 10.3354/dao01958
Oever, J., and Ten Netea, M. G. (2014). The bacteriome-mycobiome interaction and antifungal host defense. Eur. J. Immunol. 44, 3182–3191. doi: 10.1002/eji.201344405
Oksanen, J., Simpson, G. L., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., et al. (2022). Vegan: Community ecology package.
Parker, J. M., Mikaelian, I., Hahn, N., and Diggs, H. E. (2002). Clinical diagnosis and treatment of epidermal chytridiomycosis in African clawed frogs (Xenopus tropicalis). Comp. Med. 52, 265–268
Parveen, B., Mary, I., Vellet, A., Ravet, V., and Debroas, D. (2013). Temporal dynamics and phylogenetic diversity of free-living and particle-associated Verrucomicrobia communities in relation to environmental variables in a mesotrophic lake. FEMS Microbiol. Ecol. 83, 189–201. doi: 10.1111/j.1574-6941.2012.01469.x
Pessier, A. P., Nichols, D. K., Longcore, J. E., and Fuller, M. S. (1999). Cutaneous chytridiomycosis in poison dart frogs (Dendrobates spp.) and White's tree frogs (Litoria caerulea). J. Vet. Diagn. Invest. 11, 194–199. doi: 10.1177/104063879901100219
Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., et al. (2012). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. doi: 10.1093/nar/gks1219
R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing: Vienna.
Redford, K. H., Segre, J. A., Salafsky, N., del Rio, C. M., and McAloose, D. (2012). Conservation and the microbiome. Conserv. Biol. J. Soc. Conserv. Biol. 26, 195–197. doi: 10.1111/j.1523-1739.2012.01829.x
Robinson, C. J., Bohannan, B. J. M., and Young, V. B. (2010). From structure to function: the ecology of host-associated microbial communities. Microbiol. Mol. Biol. Rev. 74, 453–476. doi: 10.1128/MMBR.00014-10
Rosenberg, E. (2014). “The family Rubritaleaceae” in The prokaryotes: Other major lineages of Bacteria and the Archaea. eds. E. Rosenberg, E. F. DeLong, S. Lory, E. Stackebrandt, and F. Thompson (Berlin, Heidelberg: Springer), 861–862.
Rosenblum, E. B., Voyles, J., Poorten, T. J., and Stajich, J. E. (2010). The deadly chytrid fungus: a story of an emerging pathogen. PLoS Pathog. 6, 4–6. doi: 10.1371/journal.ppat.1000550
Sabino-Pinto, J., Bletz, M. C., Islam, M. M., Shimizu, N., Bhuju, S., Geffers, R., et al. (2016). Composition of the cutaneous bacterial community in Japanese amphibians: effects of captivity, host species, and body region. Microb. Ecol. 72, 460–469. doi: 10.1007/s00248-016-0797-6
Sanchez, E., Bletz, M. C., Duntsch, L., Bhuju, S., Geffers, R., Jarek, M., et al. (2017). Cutaneous bacterial communities of a poisonous salamander: a perspective from life stages, body parts and environmental conditions. Microb. Ecol. 73, 455–465. doi: 10.1007/s00248-016-0863-0
Sanford, J. A., and Gallo, R. L. (2013). Functions of the skin microbiota in health and disease. Semin. Immunol. 25, 370–377. doi: 10.1016/j.smim.2013.09.005
Scheele, B. C., Pasmans, F., Skerratt, L. F., Berger, L., Martel, A., Beukema, W., et al. (2019). Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science 363, 1459–1463. doi: 10.1126/science.aav0379
Scheuermayer, M., Gulder, T. A. M., Bringmann, G., and Hentschel, U. (2006). Rubritalea marina gen. Nov., sp. nov., a marine representative of the phylum ‘Verrucomicrobia’, isolated from a sponge (Porifera). Int. J. Syst. Evol. Microbiol. 56, 2119–2124. doi: 10.1099/ijs.0.64360-0
Schliep, K. P. (2011). Phangorn: phylogenetic analysis in R. Bioinformatics 27, 592–593. doi: 10.1093/bioinformatics/btq706
Schliep, K., Potts, A. J., Morrison, D. A., and Grimm, G. W. (2017). Intertwining phylogenetic trees and networks. Methods Ecol. Evol. 8, 1212–1220. doi: 10.1111/2041-210X.12760
Schloss, P. D. (2023). Rarefaction is currently the best approach to control for uneven sequencing effort in amplicon sequence analyses. mSphere 9:e0035423. doi: 10.1101/2023.06.23.546313
Shibagaki, N., Suda, W., Clavaud, C., Bastien, P., Takayasu, L., Iioka, E., et al. (2017). Aging-related changes in the diversity of women’s skin microbiomes associated with oral bacteria. Sci. Rep. 7:10567. doi: 10.1038/s41598-017-10834-9
Skerratt, L. F., Berger, L., Speare, R., Cashins, S., McDonald, K. R., Phillott, A. D., et al. (2007). Spread of chytridiomycosis has caused the rapid global decline and extinction of frogs. EcoHealth 4, 125–134. doi: 10.1007/s10393-007-0093-5
Sugden, S., St Clair, C. C., and Stein, L. Y. (2021). Individual and site-specific variation in a biogeographical profile of the coyote gastrointestinal microbiota. Microb. Ecol. 81, 240–252. doi: 10.1007/s00248-020-01547-0
Tennessen, J. A., Woodhams, D. C., Chaurand, P., Reinert, L. K., Billheimer, D., Shyr, Y., et al. (2009). Variations in the expressed antimicrobial peptide repertoire of northern leopard frog (Rana pipiens) populations suggest intraspecies differences in resistance to pathogens. Dev. Comp. Immunol. 33, 1247–1257. doi: 10.1016/j.dci.2009.07.004
van Passel, M. W. J., Kant, R., Palva, A., Copeland, A., Lucas, S., Lapidus, A., et al. (2011b). Genome sequence of the Verrucomicrobium Opitutus terrae PB90-1, an abundant inhabitant of rice paddy soil ecosystems. J. Bacteriol. 193, 2367–2368. doi: 10.1128/jb.00228-11
van Passel, M. W. J., Kant, R., Zoetendal, E. G., Plugge, C. M., Derrien, M., Malfatti, S. A., et al. (2011a). The genome of Akkermansia muciniphila, a dedicated intestinal mucin degrader, and its use in exploring intestinal metagenomes. PLoS One 6:e16876. doi: 10.1371/journal.pone.0016876
Varga, J. F. A., Bui-Marinos, M. P., and Katzenback, B. A. (2019). Frog skin innate immune defences: sensing and surviving pathogens. Front. Immunol. 9:3128. doi: 10.3389/fimmu.2018.03128
Voyles, J., Young, S., Berger, L., Campbell, C., Voyles, W. F., Dinudom, A., et al. (2009). Pathogenesis of chytridiomycosis, a cause of catastrophic amphibian declines. Science 326, 582–585. doi: 10.1126/science.1176765
Vredenburg, V. T., Bingham, R., Knapp, R., Morgan, J. A. T., Moritz, C., and Wake, D. (2007). Concordant molecular and phenotypic data delineate new taxonomy and conservation priorities for the endangered mountain yellow-legged frog. J. Zool. 271, 361–374. doi: 10.1111/j.1469-7998.2006.00258.x
Vredenburg, V. T., Knapp, R. A., Tunstall, T. S., and Briggs, C. J. (2010). Dynamics of an emerging disease drive large-scale amphibian population extinctions. Proc. Natl. Acad. Sci. USA 107, 9689–9694. doi: 10.1073/pnas.0914111107
Walke, J. B., Becker, M. H., Loftus, S. C., House, L. L., Cormier, G., Jensen, R. V., et al. (2014). Amphibian skin may select for rare environmental microbes. ISME J. 8, 2207–2217. doi: 10.1038/ismej.2014.77
Walke, J. B., and Belden, L. K. (2016). Harnessing the microbiome to prevent fungal infections: lessons from amphibians. PLoS Pathog. 12:e1005796. doi: 10.1371/journal.ppat.1005796
Wang, Q., Garrity, G. M., Tiedje, J. M., and Cole, J. R. (2007). Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267. doi: 10.1128/AEM.00062-07
Weiss, S., Xu, Z. Z., Peddada, S., Amir, A., Bittinger, K., Gonzalez, A., et al. (2017). Normalization and microbial differential abundance strategies depend upon data characteristics. Microbiome 5:27. doi: 10.1186/s40168-017-0237-y
Woodhams, D. C., Alford, R. A., Antwis, R. E., Archer, H., Becker, M. H., Belden, L. K., et al. (2015). Antifungal isolates database of amphibian skin-associated bacteria and function against emerging fungal pathogens. Ecology 96:595. doi: 10.1890/14-1837.1
Woodhams, D. C., Ardipradja, K., Alford, R. A., Marantelli, G., Reinert, L. K., and Rollins-Smith, L. A. (2007a). Resistance to chytridiomycosis varies among amphibian species and is correlated with skin peptide defenses. Anim. Conserv. 10, 409–417. doi: 10.1111/j.1469-1795.2007.00130.x
Woodhams, D. C., Geiger, C. C., Reinert, L. K., Rollins-Smith, L. A., Lam, B., Harris, R. N., et al. (2012). Treatment of amphibians infected with chytrid fungus: learning from failed trials with itraconazole, antimicrobial peptides, bacteria, and heat therapy. Dis. Aquat. Org. 98, 11–25. doi: 10.3354/dao02429
Woodhams, D. C., Kenyon, N., Bell, S. C., Alford, R. A., Chen, S., Billheimer, D., et al. (2010). Adaptations of skin peptide defences and possible response to the amphibian chytrid fungus in populations of Australian green-eyed treefrogs, Litoria genimaculata. Divers. Distrib. 16, 703–712. doi: 10.1111/j.1472-4642.2010.00666.x
Woodhams, D. C., Rollins-Smith, L. A., Carey, C., Reinert, L., Tyler, M. J., and Alford, R. A. (2006a). Population trends associated with skin peptide defenses against chytridiomycosis in Australian frogs. Oecologia 146, 531–540. doi: 10.1007/s00442-005-0228-8
Woodhams, D. C., Rollins-Smith, L. A., Reinert, L. K., Lam, B. A., Harris, R. N., Briggs, C. J., et al. (2020). Probiotics modulate a novel amphibian skin defense peptide that is antifungal and facilitates growth of antifungal bacteria. Microb. Ecol. 79, 192–202. doi: 10.1007/s00248-019-01385-9
Woodhams, D. C., Voyles, J., Lips, K. R., Carey, C., and Rollins-Smith, L. A. (2006b). Predicted disease susceptibility in a Panamanian amphibian assemblage based on skin peptide defenses. J. Wildl. Dis. 42, 207–218. doi: 10.7589/0090-3558-42.2.207
Woodhams, D. C., Vredenburg, V. T., Simon, M. A., Billheimer, D., Shakhtour, B., Shyr, Y., et al. (2007b). Symbiotic bacteria contribute to innate immune defenses of the threatened mountain yellow-legged frog, Rana muscosa. Biol. Conserv. 138, 390–398. doi: 10.1016/j.biocon.2007.05.004
Wright, E. S. (2016). Using DECIPHER v2.0 to analyze big biological sequence data in R. R J. 8:352. doi: 10.32614/RJ-2016-025
Yilmaz, P., Parfrey, L. W., Yarza, P., Gerken, J., Pruesse, E., Quast, C., et al. (2014). The SILVA and “all-species living tree project (LTP)” taxonomic frameworks. Nucleic Acids Res. 42, D643–D648. doi: 10.1093/nar/gkt1209
Yoon, J., Matsuo, Y., Matsuda, S., Adachi, K., Kasai, H., and Yokota, A. (2007). Rubritalea spongiae sp. nov. and Rubritalea tangerina sp. nov., two carotenoid-and squalene-producing marine bacteria of the family Verrucomicrobiaceae within the phylum ‘Verrucomicrobia’, isolated from marine animals. Int. J. Syst. Evol. Microbiol. 57, 2337–2343. doi: 10.1099/ijs.0.65243-0
Yoon, J., Matsuo, Y., Matsuda, S., Adachi, K., Kasai, H., and Yokota, A. (2008). Rubritalea sabuli sp. nov., a carotenoid-and squalene-producing member of the family Verrucomicrobiaceae, isolated from marine sediment. Int. J. Syst. Evol. Microbiol. 58, 992–997. doi: 10.1099/ijs.0.65540-0
Keywords: microbiome, amphibian, Rana sierrae, skin, captivity, Batrachochytrium dendrobatidis, within-individual heterogeneity
Citation: Ghose SL and Eisen JA (2025) Skin microbiomes of frogs vary among body regions, revealing differences that reflect known patterns of chytrid infection. Front. Microbiol. 16:1579231. doi: 10.3389/fmicb.2025.1579231
Edited by:
Chih-Horng Kuo, Academia Sinica, TaiwanReviewed by:
Liang-chun Wang, National Sun Yat-sen University, TaiwanAdeline Loyau, Institut National Polytechnique de Toulouse, France
Copyright © 2025 Ghose and Eisen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Sonia L. Ghose, c29uaWEuZ2hvc2VAZ21haWwuY29t