Kenen B. Goodwin
Jaren D. Hutchinson2
Zachariah Gompert- 1Department of Watershed Sciences, Utah State University, Logan, UT, United States
- 2Department of Wildland Resources, Utah State University, Logan, UT, United States
- 3Department of Biology, Utah State University, Logan, UT, United States
Host-associated microbiomes play important roles in host health and pathogen defense. In amphibians, the skin-associated microbiota can contribute to innate immunity with potential implications for disease management. Few studies have examined season-long temporal variation in the amphibian skin-associated microbiome, and the interactions between bacteria and fungi on amphibian skin remain poorly understood. We characterize season-long temporal variation in the skin-associated microbiome of the western tiger salamander (Ambystoma mavortium) for both bacteria and fungi between sites and across salamander life stages. Two hundred seven skin-associated microbiome samples were collected from salamanders at two Rocky Mountain lakes throughout the summer and fall of 2018, and 127 additional microbiome samples were collected from lake water and lake substrate. We used 16S rRNA and ITS amplicon sequencing with Bayesian Dirichlet-multinomial regression to estimate the relative abundances of bacterial and fungal taxa, test for differential abundance, examine microbial selection, and derive alpha diversity. We predicted the ability of bacterial communities to inhibit the amphibian chytrid fungus Batrachochytrium dendrobatidis (Bd), a cutaneous fungal pathogen, using stochastic character mapping and a database of Bd-inhibitory bacterial isolates. For both bacteria and fungi, we observed variation in community composition through time, between sites, and with salamander age and life stage. We further found that temporal trends in community composition were specific to each combination of salamander age, life stage, and lake. We found salamander skin to be selective for microbes, with many taxa disproportionately represented relative to the environment. Salamander skin appeared to select for predicted Bd-inhibitory bacteria, and we found a negative relationship between the relative abundances of predicted Bd-inhibitory bacteria and Bd. We hope these findings will assist in the conservation of amphibian species threatened by chytridiomycosis and other emerging diseases.
Introduction
Host-associated microbiomes can interact with their hosts in many ways. Specialized metabolites produced by microbes can influence various aspects of host biology (Sharon et al., 2014), and host production of antimicrobial peptides can in turn influence microbial community structure (McFall-Ngai et al., 2013). Microbial communities are increasingly recognized as providing beneficial and necessary services for their hosts (Dethlefsen et al., 2007; Grice and Segre, 2011), and maintaining and restoring healthy microbiomes can be important for host health (Tosh and McDonald, 2012). Host-associated microbiomes can inhibit pathogens or parasites through competition, the activation of host immune responses, and the production of inhibitory secondary metabolites (Lee and Mazmanian, 2010; Britton and Young, 2014; Grunseich et al., 2019). An imbalance in the host-associated microbiome can permit transient opportunistic pathogens and resident microbes with pathogenic potential to harm the host (Lee and Mazmanian, 2010).
Much attention has been given to the amphibian skin-associated microbiome’s role in innate immunity for its potential in disease management (Walke and Belden, 2016). Chytridiomycosis is a devastating amphibian skin disease caused by the fungal pathogen Batrachochytrium dendrobatidis (hereafter Bd; Longcore et al., 1999; Skerratt et al., 2007). Because numerous amphibian skin-associated bacteria have been found to inhibit the growth of Bd, probiotic bioaugmentation and habitat management have the potential to influence susceptibility to chytridiomycosis (Harris et al., 2009; Kueneman et al., 2016a; Grant et al., 2018). A sound understanding of host-associated microbiomes and their natural range of variation is necessary to select effective probiotics for safe and successful probiotic bioaugmentation strategies (Bletz et al., 2013).
While amphibian skin-associated microbiomes are species-specific, vary with life history stage, and are distinct from environmental microbiomes (i.e., soil, lake substrate, and lake water microbiomes), some variation in the microbiomes is attributable to location and abiotic water quality (McKenzie et al., 2011; Kueneman et al., 2013; Walke et al., 2014; Bletz et al., 2017a,b; Ellison et al., 2019). The composition of skin-associated microbial communities has been found to vary between larval and metamorphosed life stages in both frog and salamander species, with community diversity being higher in the adults of these species than their larvae (Kueneman et al., 2013, 2016b; Sabino-Pinto et al., 2017). Temperature has been found to influence operational taxonomic unit (OTU) richness and the production of antifungal metabolites in amphibian skin-associated microbiomes (Daskin et al., 2014; Muletz-Wolz et al., 2019).
Although many studies have worked to characterize species-specific and spatial variation in the amphibian skin-associated microbiome, fine-scale season-long temporal variation in natural systems remains a major gap in our knowledge of the amphibian skin-associated microbiome with few applicable studies (Sabino-Pinto et al., 2017; Bletz et al., 2017a). Since both Bd infection prevalence and amphibian skin-associated microbiomes show seasonal and year-to-year variation (Savage et al., 2011; Longo et al., 2015; Familiar López et al., 2017; Douglas et al., 2021; Basanta et al., 2022), season-long temporal variation in the amphibian skin-associated microbiome warrants investigation for its implications in disease management.
Using a database of amphibian skin-associated Bd-inhibitory bacterial isolates and their 16S rRNA gene sequences (Woodhams et al., 2015), Bletz et al. (2017a) found that despite significant changes in bacterial community structure on the skin of salamandrid newts, the relative abundances of bacteria with Bd-inhibitory potential did not change significantly during a 12-week sampling period nor across life history stages in two of the three species studied. Sabino-Pinto et al. (2017) found bacterial communities on the skin of two salamandrid newt species to change significantly between months, and also using the database of Woodhams et al. (2015), the study found the relative abundance of putative Bd-inhibitory bacteria to be higher on the skin of larvae compared to adults for one of the species. The database by Woodhams et al. (2015) contains nearly 2,000 bacterial isolates tested for Bd-inhibitory function in vitro assays, with about half of the isolates being from Central-South America. However, the application of this database to predict amphibian skin-associated microbiome Bd-inhibitory function is limited by our knowledge of how these bacterial isolates function on amphibian skin. Observing fungal responses to changes in bacterial abundances could assist in detecting bacterial-fungal relationships.
Despite the focus of many amphibian skin-associated microbiome studies on bacteria, few studies have examined how bacteria interact with non-Bd fungal taxa and other microeukaryotes on amphibian skin (Kueneman et al., 2016b, 2017; Belasen et al., 2021). For example, Kueneman et al. (2016b) found many correlations between bacterial and fungal taxa on the skin of the western toad (Anaxyrus boreas), and the authors proposed that larval stages of amphibians may depend on high relative abundances of antifungal bacteria to confer innate immunity before metamorphosis and the maturation of the host adaptive immune system (Rollins-Smith, 1998). Hence, the interactions between bacteria and fungi on amphibian skin may have substantial implications for host health and disease management.
Broadly, our study aims to investigate temporal variation in the amphibian skin-associated microbiome using the western tiger salamander (Ambystoma mavortium; hereafter salamander) as a model amphibian. In the Rocky Mountains of North America, the western tiger salamander serves as an apex predator in many fishless high alpine lakes. When the snow melts at these lakes, adult salamanders travel from upland to the lakes to breed, and some of these salamanders remain in the lakes throughout the early summer. During the summer months, eggs hatch and larval salamanders may follow several life history strategies, including metamorphosing during the same year as hatching, overwintering as larvae and metamorphosing the following year, and becoming sexually mature in the larval stage as neotenes (Sexton and Bizer, 1978). Due to their local abundance and the presence of at least one life stage throughout the warm months (June to September, hereafter warm season) at fishless high alpine lakes, the western tiger salamander is an ideal amphibian for consistently obtaining skin-associated microbiome samples throughout the warm season.
In this study, we first examine season-long temporal variation of both bacteria and fungi in the salamander skin-associated microbiome between sites and across life history stages, and we consider whether temporal trends are similar between sites and life stages. Based on these data, we identify differentially abundant microbes between salamander skin and the environment and compare the predictive ability of spatiotemporal and water quality covariates on microbial community composition. We then ask (i) whether variation in the salamander skin-associated microbiome influences predicted Bd-inhibitory function, and (ii) whether predicted Bd-inhibitory function is correlated with the relative abundance of Bd.
Materials and methods
Study sites
Salamanders were sampled from the largest of the Gibson Lakes (Franklin County, ID; 447,845 easting, 4,654,056 northing, NAD 83 UTM Zone 12; elevation: 2,579 m) and Ponds Lake (Summit County, UT; 503,020 easting, 4,503,670 northing, NAD 83 UTM Zone 12; elevation: 3,058 m). These lakes were chosen for sampling due to their differences in geology, substrate, and water conditions. We chose to sample lakes with different environmental conditions in order to investigate whether temporal trends in microbiome composition on salamander skin were similar between different lake environments. Both lakes are fishless, have no tributaries or outlets, and are located in different subranges of the Rocky Mountains. Gibson Lakes is a ~ 2.5-ha shallow lake in a limestone basin of the Bear River Mountains. Patches of submerged vegetation cover much of the lake bottom, and the lake substrate is primarily composed of soft sticky mud. Ponds Lake is a ~ 2.3-ha lake in a granitic basin of the Uinta Mountains. The lake substrate is a thick layer of loose vegetative material, and some parts of the shoreline have floating mats of vegetation. The water in Ponds Lake is stained red with dissolved organic carbon.
In 2018, access to Gibson Lakes was blocked due to snow at lower elevations until June 9th, when salamander eggs were observed attached to submerged vegetation. By the next week, when field sampling began, most of the previously observed eggs had hatched. Data from NRCS SNOTEL sites (see Supplementary material) suggest that snow melted at both lakes at about the same time in 2018, possibly within days of each other, and snow typically melts at these lakes about a week apart. Based on these data, it is likely that salamanders laid eggs in 2018 at about the same time at both lakes.
Sampling design
To ensure that the sampled salamanders were distributed throughout the lakes, the lakes were sampled by strata. Gibson Lakes was assigned 4 strata and Ponds Lake was assigned 5 strata (Figure 1). The strata divide the lakes into regions based on natural landmarks for easy recognition in the field. Within a lake, all strata had roughly the same area, and their areas remained roughly the same as each other as water levels dropped throughout the warm season. Three age classes of salamanders could be distinguished based on length and weight measurements, age-0, age-1, and age-2+. These age classes were of distinctly different sizes, with the length and weight of each age class increasing throughout the warm season (Supplementary Figure S1). During each visit to a lake (hereafter sampling event), we collected up to 20 salamanders from each age class with a maximum of five and four salamanders per stratum at Gibson Lakes and Ponds Lake, respectively. Each lake was sampled every other week during the 2018 warm season. Sampling began shortly after snowmelt and continued until the lakes became too cold to safely catch salamanders. Sampling began at Gibson Lakes on June 16th and Ponds Lake on June 23rd. Gibson Lakes was too cold to sample on September 29th, marking the end of the field season.
Figure 1. Strata for Gibson Lakes (left) and Ponds Lake (right). Images of (A) a young-of-year larval salamander and (B) a metamorphosed adult salamander.
Salamanders were considered larvae if they retained any of their larval gill structures, and salamanders were considered metamorphosed individuals once all traces of their gill structures were absorbed. For each age class, larval and metamorphosed individuals were encountered, which we refer to as life stages, and we refer to the six possible combinations of age class and life stage as stage classes. We expect most age-2+ individuals to be sexually mature adults, at which point gilled individuals are considered neotenes.
Data collection
Upon arriving at a lake, environmental microbiome samples and water quality data were collected. During the first visit to each lake, a location was selected just offshore in each stratum to collect these samples and data. These locations were chosen to have relatively homogeneous depths across strata and to minimize the distance which the sampling location would need to move with receding water levels. Water quality data was collected prior to collecting environmental microbiome samples to minimize disturbance to the water. Water temperature, pH, electrical conductivity, and dissolved oxygen (ppm and percent) were measured just below the water surface using handheld meters (Hannah Instruments HI98129 and HI9146). For sampling the lake water microbiome, 500 ml of lake water was collected from the water surface in a laboratory Nalgene bottle. Following collection of a lake water microbiome sample, a lake substrate microbiome sample was collected from the top ~10 cm of pond substrate using a small PVC clam gun. The substrate column was deposited into a 15-ml conical tube, and excess water was decanted. The substrate was thoroughly stirred with a teasing needle, and ~ 1.5 ml of substrate was deposited into a sterile 2-ml microcentrifuge tube. The microcentrifuge tubes containing substrate samples were placed on ice in a cooler while in the field. New latex gloves were worn for each environmental microbiome sample, and the clam gun and teasing needle were rinsed with 95% ethanol between substrate samples. The clam gun and teasing needle were rinsed with 6% bleach solution followed by a thorough rinse with distilled water between sampling events. Nalgene bottles were rinsed thoroughly with distilled water and autoclaved for 20 min at 121°C between holding lake water samples. During each sampling event, the depth of a predefined rock was measured to determine relative lake elevation, the water level of the lake relative to its height at the beginning of the warm season.
After collecting environmental microbiomes and water quality data for all strata, salamanders were captured for each stratum. Salamanders were collected by hand and dip net, and salamanders were stored in 5-gallon buckets filled with lake water. For each stratum, different age classes were stored in separate 5-gallon buckets to reduce the risk of smaller salamanders being harmed from predation or aggression from larger individuals. While storing salamanders from the same age class and stratum together in 5-gallon buckets could have allowed for microbial contamination between individuals, we suspect that potential contamination between individuals was minimal for the following reasons. First, only a few individuals were stored together at a time (an average of 3.04 and maximum of five individuals). Second, the period of time which individuals were stored together was short (typically about 25 min). Finally, individuals tended to disperse themselves relatively evenly within the buckets, and contact between individuals and the ventral surfaces (the body region of interest) of others was rare.
Each salamander was handled with new latex gloves, and snout-vent length (SVL) and weight measurements were taken to verify age classes (Supplementary Figure S1). Sex was determined for age-2+ salamanders. The ventral surface of each salamander was rinsed with 50 ml of distilled water (Bletz et al., 2017a) to remove environmental material and transient microbes (Culp et al., 2007; Lauer et al., 2007), and the salamander’s ventral surface was swabbed with a sterile rayon-tipped swab (MW113 Medical Wire and Equipment, Corsham, United Kingdom). Swabbing was performed by stroking the swab across the ventral surface ten times (one time = an up and back stroke along the full length of the belly; Bletz et al., 2017a). Swabs used to sample salamander skin-associated microbiomes were stored in individual sterile 2-ml microcentrifuge tubes and placed on ice in a cooler while in the field. After processing salamanders for a stratum was complete, the salamanders were released back into the stratum, and salamander collection began at the next stratum. While it is possible that salamanders sampled in one stratum may have been sampled again in another stratum during the same sampling event, few salamanders were observed to have swum far from their point of release, which was away from stratum borders.
For each sampling event, the lake was surveyed for salamanders for a minimum of 5 person-hours divided evenly among the lake’s strata. Salamanders were processed after the stratum minimum sampling time was reached or the maximum number of individuals from all available age classes had been collected, and the search for salamanders then proceeded to the next stratum. After field sampling and while still at the lake, wet and dry negative control swabs were taken. Wet control swabs were sprayed with 50 ml of distilled water, and nothing was done to the dry control swabs. Wet and dry control swabs were placed in individual sterile 2-ml microcentrifuge tubes and stored on ice in a cooler while in the field.
Following field sampling and on the same day, lake water samples were prefiltered through a 5.0-μm prefilter membrane to remove debris followed by filtration with a 0.22-μm filter membrane to catch microbes (Millipore Sigma SVLP02500 and GSWP04700, respectively). Multiple 5.0-μm prefilter membranes were used for each water sample as necessary, whereas samples which experienced clogging on the 0.22-μm filter membrane (three samples) were discarded. Following filtration, 0.22-μm filter membranes were folded and stored in 2-ml microcentrifuge tubes. For autoclavable filtration equipment, the equipment was rinsed thoroughly with distilled water between water samples followed by autoclaving for 20 min at 121°C. Non-autoclavable filtration equipment was rinsed with 6% bleach solution followed by a thorough rinse with distilled water between water samples. Every four or five sampling events, five 500-ml distilled water samples were filtered as negative controls.
All samples were transferred to a −80°C freezer for storage, and the typical time from field collection to freezer storage was about five-and-a-half hours. Salamanders were collected, stored, handled, and released according to an approved Utah State University Institutional Animal Care and Use Committee protocol (#2798), a Utah Division of Wildlife Resources Certificate of Registration (#2COLL10232), and an Idaho Department of Fish and Game Wildlife Collection/Banding/Possession Permit (#180110).
DNA extraction and library preparation
DNA was extracted with the DNeasy PowerSoil Pro Kit (Qiagen, Inc.) following the manufacturer’s protocol, and 12 empty extractions were performed as blank negative controls. Substrate samples were centrifuged for 30 s at 10,000 × g, excess liquid was removed with a pipette, and a scoopula was used to collect 250 mg of substrate from each sample for DNA extraction. Water sample filter membranes were finely diced using scissors and forceps into reagent reservoirs before being transferred to DNA extraction tubes. Swab samples were transferred to DNA extraction tubes using a different pair of forceps than that used for water samples. Pre-DNA extraction sample preparation work was performed under a fume hood, and the scoopula, scissors, and forceps were rinsed with 95% ethanol, flamed, and rinsed thoroughly with distilled water between samples. Reagent troughs were rinsed thoroughly with distilled water and autoclaved for 20 min at 121°C between water samples.
Following DNA extraction, two samples of ZymoBIOMICS Microbial Community DNA Standard (Zymo Research D6305) were added as mock community positive sequencing controls. 6 μl of a control oligo pool was added to 30 μl of full concentration DNA extract. The control oligo pool contained 0.01 pg/μl each of 16S and ITS well-specific cross contamination oligos (hereafter coligos; Hawkins et al., 2018) and 0.03 pg/μl each of 16S and ITS synthetic genes (hereafter synthgenes; Tourlousse et al., 2017). The addition of fixed amounts of 16S and ITS synthgenes to a constant volume of each sample’s DNA extract will be used later in estimating the amount of microbial DNA in each sample. Sample DNA concentrations were measured via absorption and normalized to 10 ng/μl with an automated liquid handler. Combinatorial dual indexing was performed on the samples with two-stage polymerase chain reaction (PCR). First stage PCR amplified the 16S rRNA and ITS genetic barcoding regions, added unique dual index combinations to each sample, and added a portion of the Illumina Nextera adapter. For each sample, two first-stage PCR replicates were performed and subsequently pooled. Second stage PCR completed Illumina adapter addition. The 16S rRNA V4 region was amplified using the primers 515F (forward; Parada et al., 2016) and 806R (reverse; Caporaso et al., 2011). The ITS1 region was amplified using the primers ITS1-F (forward; Gardes and Bruns, 1993) and ITS2 (reverse; White et al., 1990). A modified AxyPrep MagBead PCR Clean-up protocol was used to purify the amplified DNA after each PCR reaction. Library preparation occurred at the University of Wyoming Genome Technologies Laboratory (Laramie, WY). See Supplementary material for library preparation details.
DNA sequencing and processing
Paired-end DNA sequencing of pooled amplicon product was performed on both Illumina MiSeq (v3 600-cycle kit, 2 × 300 base pair [bp] reads) and Illumina NextSeq (v2 300-cycle kit, 2 × 150 bp reads) platforms at the Utah State University Center for Integrated Biosystems (Logan, UT). Both sequencing platforms offer their own advantages for 16S and ITS amplicon sequencing, where Illumina MiSeq produces longer but fewer reads than Illumina NextSeq. The longer MiSeq sequences provide greater taxonomic resolution, and the greater number of NextSeq sequences reduces uncertainty in relative abundance estimates. We leverage the benefits of both sequencing platforms by using the longer-length MiSeq sequences as study-specific 16S and ITS reference libraries to enhance the taxonomic resolution of our shorter but more numerous NextSeq sequences. Illumina MiSeq produced 19 million paired-end reads, and Illumina NextSeq produced 187 million paired-end reads.
MiSeq reads were partitioned into 16S and ITS datasets based on their primer regions using a custom Perl script (version 5.18.1; see Data Availability for script), and index tags were removed. Since variable length index tags were used, MiSeq reads were trimmed to 290 bp using cutadapt (version 2.10; Martin, 2011) to ensure that non-overlapping sequences did not appear different simply due to read length. Using cutadapt, read pairs that contained Ns were removed, and forward primers and reverse complements of reverse primers were trimmed (with a maximum error rate of 0.15, a minimum trimmed length of 1 bp, and discarding untrimmed read pairs), with trimming the reverse primer’s reverse complement being required for 16S read pairs. Variable length index tags reduce amplicon sequencing error on Illumina platforms by increasing heterogeneity in the composition of bases called in each cycle (Fadrosh et al., 2014).
The DADA2 bioinformatics pipeline (version 3.10; Callahan et al., 2016) was used in the R statistical software program (version 4.0.2; R Core Team, 2020) for quality filtering, phiX removal, denoising, merging pairs, chimera removal, and taxonomic assignment of MiSeq reads (see Supplementary material for details). While 16S reads were of appropriate lengths for merging pairs, the variable length of the ITS region resulted in both overlapping and non-overlapping read pairs. DADA2 has the ability to work with both overlapping and non-overlapping read pairs, allowing for the retention of fungal taxa with long ITS genes. Overlapping ITS read pairs were merged, while non-overlapping ITS read pairs were retained in the pipeline as concatenated sequences with 10-N spacers, which DADA2’s implementation of the naïve Bayesian classifier is designed to work with. DADA2’s naïve Bayesian classifier (Wang et al., 2007) was used to classify unique sequences in the MiSeq 16S and ITS datasets using Silva (version 138; Quast et al., 2012) and UNITE (general dynamic FASTA release for fungi; version 8.2; Nilsson et al., 2019) reference libraries, respectively. To create study-specific 16S and ITS reference libraries, NextSeq-length forward and reverse reads were created from the classified MiSeq 16S and ITS sequences, and consensus taxonomies and MiSeq-length sequences (for predicting Bd-inhibitory function) were generated for duplicate reference read pairs (see Supplementary material for details). Integers were appended to reference taxa names to differentiate each amplicon sequence variant (ASV) associated with a taxon.
NextSeq reads were assigned to PCR replicate, barcode region, and sample (i.e., reads were demultiplexed) using Perl while allowing 1 bp mismatches in the index tags (index tags were designed to differ by at least 2 bp). Allowing 1 bp mismatches in the index tags allows reads which experience sequencing error in the index tag regions to be retained if the index tags can still be uniquely identified. We note that allowing sequencing errors in index tag regions is not uncommon during demultiplexing. For example, demultiplexing in cutadapt and QIIME 2 (qiime cutadapt demux-paired command; version 2022.8; Bolyen et al., 2019) allow for 10% mismatches in index tags by default. Following demultiplexing, phiX reads were discarded and index tags were removed using Perl. The following steps were performed sequentially on the NextSeq reads using cutadapt: reads were trimmed to 140 bp to make all reads the same length, read pairs with Ns were removed, forward primers and reverse complements of reverse primers were trimmed (with the same settings as the MiSeq data but without requiring trimming of the reverse primers’ reverse complements).
Using exact matching in R, 21.4 million of 54.3 million NextSeq 16S reads were identified to 15,792 reference sequences, and 60.4 million of 113.9 million NextSeq ITS reads were identified to 3,488 reference sequences. Of the identified NextSeq sequences, 17.0% of 16S sequences were coligos or the synthgene, and 79.5% of ITS sequences were coligos or the synthgene. All samples were checked for between-well cross contamination through use of the coligos. Three salamander samples and one blank control sample were removed from the 16S dataset due to high amounts of between-well contamination (having a ratio of any contaminant coligo to non-contaminant coligo greater than 0.1 after summing coligo read counts across PCR replicates). Two salamander samples were removed from the ITS dataset due to lack of detection of any non-synthgene and non-coligo sequences in both PCR replicates. Coligos were removed from the datasets for all subsequent analyses. In the mock community samples, we observed strong amplification bias in the ITS data (Supplementary Figure S2), and one fungal taxon was split into three substantial ASVs. In an effort to mitigate the potential impact of fungal taxa being split into multiple ASVs, we merged fungal ASVs which were assigned the same taxonomy into the same taxa. We chose to forego rarefaction of our samples as it increases uncertainty in relative abundances (McMurdie and Holmes, 2014).
We performed principal component analyses (PCAs) on the proportional abundances of taxa across PCR replicate and sample type (Supplementary Figures S3–S5). Taxa proportional abundances within samples were similar across PCR replicates, so read counts were summed across PCR replicates for each sample. There were ten salamander samples which grouped closely with wet swab and dry swab negative controls in the 16S PCAs on sample type, so these samples were removed from the 16S data for all subsequent analyses. Following Harrison et al. (2021), we used synthgene read counts to calculate the amount of microbial DNA in each sample relative to the synthgene (i.e., microbial read count divided by synthgene read count), and we compared the amount of microbial DNA in field samples to their associated negative controls (Supplementary Figure S6). Synthgenes were subsequently removed from the datasets. The final datasets for our field samples contained 15,690 bacterial taxa (6,529 for salamander, 8,873 for water, and 14,591 for substrate) and 469 fungal taxa (289 for salamander, 224 for water, and 413 for substrate).
Water quality between sites and through time
To examine how water quality changed throughout the warm season, we fit a linear mixed-effects model for each water quality parameter (i.e., temperature, pH, conductivity, and dissolved oxygen [ppm and %]) using the lmerTest R package (version 3.1.3; Kuznetsova et al., 2017). In these linear mixed-effects models, stratum was treated as a random effect, and site, week, and their interaction were included as fixed effects. Stratum was coded with nine values representing the four strata in Gibson Lakes and the five strata in Ponds Lake. Site was treated as a categorical predictor, and week was treated as a continuous predictor. Week represented the number of weeks since June 9th, 2018. Water quality measurements are displayed in Supplementary Figure S7.
Predicting Bd-inhibitory function
A database of amphibian skin-associated microbiome Bd-inhibitory bacterial isolates (Woodhams et al., 2015) was used to predict which bacteria observed in our datasets exhibit Bd-inhibitory properties (see Supplementary material for details). We trimmed sequences in the database of Woodhams et al. (2015) to the 16S rRNA V4 region using our 16S amplification primers with R, and we aligned the MiSeq 16S sequences of taxa detected in our NextSeq 16S field samples with the Woodhams et al. (2015) sequences using Clustal Omega (version 1.2.4; Sievers et al., 2011). We used FastTree 2 (version 2.1.11; Price et al., 2010) to create a phylogenetic tree, and we used stochastic character mapping with the make.simmap function in the phytools package (version 0.7.70; Revell, 2012) to predict the Bd-inhibition statuses of our observed taxa. Stochastic character mapping extends ancestral state reconstruction to probabilistically predict unobserved traits at the tips of a phylogenetic tree (Bollback, 2006). While existing applications of the Woodhams et al. (2015) database tend to employ local alignment or clustering methods to classify bacterial taxa as “potentially” Bd-inhibitory (e.g., Kueneman et al., 2016b; Bletz et al., 2017a; Kruger, 2020), stochastic character mapping provides the benefit of yielding probabilistic predictions that bacterial taxa are actually Bd-inhibitory. We further note that extended ancestral trait reconstruction is commonly applied in predicting the metabolic function of gut microbiomes (Langille et al., 2013). We visualized our phylogenetic tree with the posterior probabilities of our taxa being Bd-inhibitory using the Interactive Tree of Life (Supplementary Figure S8; version 6.5.4; Letunic and Bork, 2021).
The vast majority of our taxa had low confidence in their Bd-inhibition statuses (99.5% of posterior probabilities were between 47.9 and 52.5%), whereas most posterior probabilities which were < 40% or > 60% were also ≤10% or ≥ 90% (33 of 39). Therefore, we considered our bacterial taxa to be Bd-inhibitory if their posterior probabilities of Bd-inhibition were ≥ 90%, and we considered our bacterial taxa to be non-Bd-inhibitory if their posterior probabilities of Bd-inhibition were ≤ 10%. Otherwise, we considered our bacterial taxa to have an uncertain Bd-inhibition status.
Microbial composition modeling
For both bacterial and fungal communities, we fit Bayesian Dirichlet-multinomial regression models to the salamander, water, and substrate microbiome data to identify differentially abundant microbes and to evaluate differences in overall community composition. Bayesian Dirichlet-multinomial regression estimates the effect of covariates on a set of proportions which sum to one (i.e., a simplex) and uses a set of counts as a multivariate response. In the context of microbiome data, the model uses read counts to estimate the expected proportional abundances of microbial taxa in the community, and the model considers the underlying uncertainty in each sample’s composition, which is dictated by the sample’s total read count (i.e., its sampling effort). Bayesian Dirichlet-multinomial models outperform other analyses of compositional data in detecting differences in community composition, and the model further allows for the identification of the taxa responsible for those differences (i.e., the model allows for differential abundance testing; Harrison et al., 2020).
Our Bayesian Dirichlet-multinomial regression model was adapted from the Bayesian Dirichlet regression model of Sennhenn-Reulen (2018), and we used backwards variable selection by widely applicable information criterion (WAIC) to optimize predictive accuracy. WAIC is a Bayesian analog to AIC and approximates the predictive accuracy of leave-one-out cross-validation (Gelman et al., 2014). In our model, sample read counts are distributed according to the Dirichlet-multinomial distribution. Each taxon receives a linear predictor combination, and the softmax function (a multivariate inverse logit) normalizes linear predictor combinations for all taxa into expected proportions. The last taxon serves as a reference category, and its intercept and regression coefficients are set to zero to allow for model identifiability. A precision parameter controls the degree of overdispersion relative to the multinomial distribution. See Supplementary material for model details.
Our Bayesian Dirichlet-multinomial regression models were computationally intensive to fit, with the number of model parameters and model run time increasing with the number of taxa included. To keep model run-times practical, we opted to select the 100 most proportionally abundant taxa from the salamander samples, plus an “other” category, for inclusion in the composition models. To select these taxa, we calculated the proportion of reads of each taxon in each salamander sample. For each barcode region (i.e., 16S or ITS), we then averaged each taxon’s proportion of reads by combinations of site and life stage. We then averaged across these averages, and we took the 100 taxa with the highest averaged proportion of reads for each barcode region for use in modeling. Other taxa which were not included in the top 100 for each barcode region had their read counts merged into an “other” category. By weighting each combination of site and life stage equally (i.e., by taking averages of averages) in selecting the top 100 taxa, we ensured that the top 100 taxa were not dominated by taxa from one site or life stage simply due to differences in sample size. Including the “other” category in the models ensures that the proportional abundances of the top 100 taxa remain unbiased. We chose to select the top 100 taxa because these comprise the vast majority of reads in salamander samples (93.1% of bacterial and 98.6% of fungal reads). As such, we expect variation in the composition of these taxa to represent most of the variation in community composition. Additionally, following initial testing, we deemed including 101 categories in the composition models (the top 100 taxa plus the “other” category) to be near the upper reasonable limit of our computing capacity on a high-performance computing cluster. Ultimately, our Bayesian Dirichlet-multinomial regression models, coupled with our backwards variable selection approach, took 785 CPU days to run. Datasets used in the modeling of water and substrate microbial communities included the same taxa as used for the salamander modeling, plus their own “other” categories. Proportional abundance estimates from Bayesian Dirichlet-multinomial regression models for water and substrate are later used to identify microbes which are disproportionately abundant on salamander skin relative to the environment. Since not all top microbial taxa in the salamander samples were detected in the water and substrate samples, the water and substrate datasets used in modeling had fewer than 101 taxa.
As water quality was highly correlated with space and time (Supplementary Figure S7), we fit models with two different sets of predictors. One predictor set included spatiotemporal covariates, while the other predictor set substituted spatiotemporal covariates with water quality. The spatiotemporal predictor set included four-way interactions between age, life stage, site, and a second-degree polynomial for week, all lower-level interactions, and the individual predictors. Stratum was also included as a predictor and treated as a hierarchical effect. The water quality predictor set included a five-way interaction between age, life stage, temperature (°C), pH, and dissolved oxygen (ppm), all lower-level interactions, and the individual predictors. Models for water and substrate lacked age and life stage predictors. Site and life stage were treated as categorical predictors, and age and week were treated as continuous predictors. Age took whole integers from zero (age-0) to two (age-2+), and week represented the number of weeks since June 9th, 2018.
Bayesian Dirichlet-multinomial regression models were fit in Stan (version 2.21.0; Carpenter et al., 2017) using the rstan R interface (version 2.21.2; Stan Development Team, 2020) with four Hamiltonian Monte Carlo (HMC) chains, 500 warmup iterations, 500 sampling iterations, and no thinning. Stan was chosen for its efficient HMC algorithm, and HMC chains were run in parallel on a University of Utah high-performance computing cluster. Gelman-Rubin convergence diagnostics () and trace plots of the posteriors were used to assess model convergence. Since interpreting the effect of Dirichlet-multinomial regression coefficients on proportional abundances is not straightforward (see Supplementary material for a discussion), we opted for a graphical interpretation of the best-fit models (i.e., the models selected by backwards variable selection). We generated posterior predictions of proportional abundances for each combination of non-stratum predictors observed in the datasets, where predictions were for the average stratum (see Supplementary material for prediction details). The posterior predictions of taxa proportional abundances were summarized with 95% credible intervals, 50% credible intervals, and their median values.
For salamander microbiomes, Hill’s diversity index with α = 2 (Haegeman et al., 2013) was derived from the posterior predictions of taxa proportional abundances. By treating Hill’s diversity as a derived parameter from the Bayesian Dirichlet-multinomial regression models, we propagated the uncertainty associated with taxa proportional abundances to our diversity index. We chose Hill’s diversity (α = 2) as our alpha diversity index because it is insensitive to the many rare taxa expected in microbial communities, and can therefore be robustly estimated from microbiome data (Haegeman et al., 2013). Because of its insensitivity to rare taxa, we expected our grouping of rare taxa into an “other” category in the Bayesian Dirichlet-multinomial regression models to have a negligible impact on Hill’s diversity. We tested this expectation in the context of our data as follows. By grouping rare taxa into an “other” category, we created a situation of maximum unevenness within the group (i.e., the entire abundance of the “other” category was composed of a single taxon). We evaluated the sensitivity of Hill’s diversity to rare taxa by considering how evenly distributing the abundance of the “other” category across all of its member taxa influenced the index. For each HMC sample, we split the proportional abundance predictions of the “other” category into its individual members with uniform proportional abundances. We then re-derived Hill’s diversity and compared the posterior medians with the original values. After excluding a single combination of stage class and sampling event for bacterial communities, we observed a very strong correlation in Hill’s diversity estimates between the two methods (Pearson correlation coefficient of >0.999 for both bacteria and fungi). The excluded combination of stage class and sampling event was age-0 metamorphosed salamanders at Gibson Lakes on September 15th, which had a high proportional abundance of “other” bacterial taxa (posterior median of 0.196) and whose diversity estimate was found to be sensitive to the grouping of rare taxa into an “other” category (Hill’s diversity increased by 226.6% in the described test). We removed this combination of stage class and sampling event from our Hill’s diversity estimates for bacteria so that all remaining diversity estimates were reliable. The posterior distributions of Hill’s diversity were summarized with 95% credible intervals, 50% credible intervals, and their median values.
To estimate the Bd-inhibitory function of bacterial communities, we summed the read counts of bacteria belonging to each Bd-inhibition category within each sample, and we fit additional Bayesian Dirichlet-multinomial regression models for salamander, water, and substrate samples with three response categories representing the Bd-inhibition statuses (i.e., Bd-inhibitory, non-Bd-inhibitory, and uncertain Bd-inhibition status). These models used the spatiotemporal predictor set, and backwards variable selection by WAIC was again used to optimize predictive accuracy. We again generated posterior predictions of proportional abundances for each combination of non-stratum predictors observed in the datasets, where predictions were for the average stratum. We summarized the proportional abundances of each Bd-inhibition category with 95% credible intervals, 50% credible intervals, and their median values.
To examine which taxa and Bd-inhibition categories were disproportionately abundant on salamander skin relative to the environment, we considered the proportional abundance of microbes that salamanders experience in their environments to be a mixture between water and substrate proportional abundances, with the mixing ratio being a product of salamander behavior. Although we do not know this ratio, we expect that the result of this mixture is between the lower of the 0.025 quantiles (the lower ends of the 95% credible intervals) and the upper of the 0.975 quantiles (the upper ends of the 95% credible intervals) of the proportional abundance posterior predictions from water and substrate Bayesian Dirichlet-multinomial regression models, and we consider this range to represent the proportional abundance of a taxon or Bd-inhibition category in the environment. In determining this range, if a taxon was not detected in the water or substrate samples, and therefore was not included in the Bayesian Dirichlet-multinomial regression modeling for that sample type, it was considered to have 0.025 and 0.975 quantiles of proportional abundance predictions for that sample type of zero. The end result is that, if the proportional abundance of a taxon on salamander skin is higher than this range, then the taxon is disproportionately more abundant on salamander skin compared to both water and substrate. Conversely, if the proportional abundance of a taxon on salamander skin is lower than this range, then the taxon is disproportionately more abundant in both water and substrate than on salamander skin.
Relationship between Bd-inhibitory bacteria and Bd
We detected Bd from ITS amplicon sequencing (i.e., there were fungal microbiome reads which were classified as Bd) on salamander skin at both lakes, and since Bd was absent in all negative control samples, we are confident that this was not the result of contamination. To verify that fungal microbiome reads which the naïve Bayesian classifier assigned to Bd were likely classified correctly, we performed an online nucleotide BLAST search (Zhang et al., 2000) with default settings for each of the 25 Bd ASVs which were previously merged into the Bd taxon. For each Bd ASV, the best-matching BLAST hit (the match with the lowest E-value) was a reference sequence belonging to Bd. 92% of Bd ASVs had 96.7% similarity or greater compared to their best-matching Bd reference sequence, and all ASVs had at least 94.5% similarity compared to their best-matching Bd reference sequence.
We tested for a relationship between the relative abundances of Bd-inhibitory bacteria and Bd on the skin of metamorphosed salamanders by fitting a Bayesian beta-binomial regression model with a logit link. We restricted this analysis to metamorphosed individuals because cutaneous Bd infections do not typically produce disease in larval amphibians (Marantelli et al., 2004). In our samples, Bd was detected on the skin of only 2.2% of larval or neotenic individuals (3 of 139) compared to 57.6% of metamorphosed individuals (38 of 66). The Bayesian beta-binomial regression model can be viewed as a univariate version of our earlier Bayesian Dirichlet-multinomial regression model. The Bayesian beta-binomial regression accounts for uncertainty in response values (i.e., the proportional abundance of Bd in fungal communities) by considering Bd read counts to be beta-binomially distributed, and a precision parameter controls the degree of overdispersion relative to the binomial distribution. Our model additionally accounts for uncertainty in the predictor values (i.e., the proportional abundance of Bd-inhibitory bacteria in bacterial communities) by estimating predictor values within the model from Bd-inhibitory read counts. Within the model, Bd-inhibitory read counts are beta-binomially distributed, and the logit linearizes Bd-inhibitory proportional abundance estimates for use as predictor values. See Supplementary material for model details. The model was fit in JAGS (version 4.3.0; Plummer, 2003) using the rjags R interface (version 4.10; Plummer, 2019) with three Markov chain Monte Carlo (MCMC) chains, 20,000 adaptation iterations, 20,000 warmup iterations, 100,000 sampling iterations, and no thinning. Gelman-Rubin convergence diagnostics and trace plots of the posteriors were used to assess model convergence.
Results
Field sampling
We observed higher amounts of microbial DNA in field samples compared to their associated negative controls (Supplementary Figure S6). Total sample counts are included in Table 1, and a breakdown of salamander skin-associated microbiome samples are displayed in Figure 2. The residency of different salamander age classes varied through time, and age-0 salamanders were too small to sample during the early warm season.
Table 1. Microbiome sample counts.
Figure 2. Counts of salamander skin-associated microbiome samples through time.
Based on observed sizes of males and females, most age-2+ salamanders – and only age-2+ salamanders – are thought to have been adults. Males develop swollen cloacas once sexually mature (Stebbins, 2003), and only one non-male age-2+ salamander (83 mm SVL; assumed to be female) had an SVL less than the smallest male (84 mm), with other small age-2+ salamanders in the range of 85 to 87 mm SVL being a mix of males (3) and females (4). Given the overlap in size between males and females, few subadults are expected to have been included in the age-2+ age class since male salamanders of this size were showing clear signs of sexual maturity. The absence of swollen cloacas from all age-0 and age-1 individuals suggests that only age-2+ individuals were sexually mature. 23 of 55 (41.8%) of sexed age-2+ salamanders were male (36.0% for Gibson Lakes and 46.7% for Ponds Lake).
Parameter estimates, test statistics, and p-values from the linear mixed-effects models for water quality are reported in Supplementary Table S1. These models included stratum as a random effect, and site, week, and their interaction were included as fixed effects. From the linear mixed-effects models, all water quality parameters (i.e., temperature, pH, conductivity, and dissolved oxygen [ppm and %]) changed significantly throughout the warm season (values of p ≤ 0.05 for all regression coefficients for week). We found a significant effect of site for all water quality parameters besides temperature, and we found significant interactions between site and week for conductivity and dissolved oxygen (both ppm and %). These results suggest that we were unable to detect differences in water temperature between Gibson Lakes and Ponds Lake, and temperature at both lakes decreased throughout the warm season (βweek = −0.713, value of p < 0.001). pH was lower at Ponds Lake compared to Gibson Lakes (βsite = −1.048, value of p < 0.001) and increased throughout the warm season at both lakes (βweek = 0.053, value of p < 0.001; non-significant interaction between site and week, value of p = 0.147). For conductivity and dissolved oxygen (both ppm and %), temporal trends were dependent on the lake (p-values ≤0.05 for all site, week, and interaction regression coefficients). Field measurements of water quality are shown in Supplementary Figure S7.
Predictions of Bd-inhibitory function
Only 33 of the 15,690 taxa detected in our NextSeq 16S field samples were classified as Bd-inhibitory or non-Bd-inhibitory (i.e., posterior probabilities ≥90% or ≤ 10%; Supplementary Table S2). Of the 872 Woodhams et al. (2015) sequences which were used in the alignment, there were 361 unique sequences, and 79 of these unique sequences occurred across multiple bacterial isolates. We note that 41 of these 79 sequences (51.9%) had inconsistent Bd-inhibition statuses (i.e., statuses varied across isolates associated with the same sequence). We also note that the aligned Woodhams et al. (2015) sequences provided limited phylogenetic coverage of the bacterial taxa detected in our field samples (Supplementary Figure S8), with only 5,330 (34.0%) of our bacterial taxa belonging to phyla included in the Woodhams et al. (2015) database.
Microbial composition
The proportion of reads in salamander samples belonging to each microbial class for each combination of stage class and sampling event are displayed in Figures 3, 4 for bacterial and fungal communities, respectively. Salamander bacterial communities were dominated by members of the phylum Proteobacteria and class Gammaproteobacteria (Figure 3), which comprised 89.8 and 86.7% of reads, respectively. The proportion of reads belonging to each fungal class (Figure 4) were more balanced than bacteria. Among salamander fungal reads, 41.6% belonged to class Rhizophydiomycetes, 21.4% were unidentified fungi, and 15.2% belonged to class Tremellomycetes.
Figure 3. Proportion of bacterial reads in salamander samples belonging to each class. The ten classes with the highest number of salamander sample reads are displayed along with a category for the other classes. A stacked bar chart is displayed for each combination of stage class and sampling event, and each stacked bar chart represents reads pooled across samples belonging to the combination of stage class and sampling event. A “G” proceeds the dates of sampling events at Gibson Lakes, and a “P” proceeds the dates of sampling events at Ponds Lake.
Figure 4. Proportion of fungal reads in salamander samples belonging to each class. The ten classes with the highest number of salamander sample reads are displayed along with a category for the other classes. A stacked bar chart is displayed for each combination of stage class and sampling event, and each stacked bar chart represents reads pooled across samples belonging to the combination of stage class and sampling event. A “G” proceeds the dates of sampling events at Gibson Lakes, and a “P” proceeds the dates of sampling events at Ponds Lake.
Bayesian Dirichlet-multinomial regression models with spatiotemporal predictors fit better than models with water quality predictors for salamander samples (Table 2), suggesting that our spatiotemporal predictors were better able to predict salamander microbial composition. All but two of the best-fitting Bayesian Dirichlet-multinomial regression models with spatiotemporal predictors included stratum as a predictor, suggesting compositional variation in microbial communities within the lakes, with the models for salamander fungal communities and water Bd-inhibition categories being the exceptions. Except for stratum in the aforementioned models, all best-fitting spatiotemporal Bayesian Dirichlet-multinomial regression models included all individual predictors or their interactions, suggesting that all of our measured variables contributed to our ability to predict the composition of microbial communities.
Table 2. Predictors included in the best-fitting Bayesian Dirichlet-multinomial regression models for microbiome composition.
The ten Bayesian Dirichlet-multinomial regression models from the backwards variable selection process with the lowest WAIC values for each sample type (i.e., salamander, water, or substrate), predictor set (i.e., spatiotemporal or water quality), and microbial community type (i.e., bacterial community, fungal community, or Bd-inhibition categories) are included in Supplementary Tables S3–S7. Four best-fitting Bayesian Dirichlet-multinomial regression models had other models within two WAIC (i.e., the model for substrate bacterial community composition and models for salamander, water, and substrate Bd-inhibition categories). With one exception, all other models within two WAIC contained the same predictors or their interactions as the best-fitting models. For example, a model with age as a predictor and another model with an interaction between age and site both contain age. The exception was one model within two WAIC of the best-fitting model for Bd-inhibition categories in lake water, which included stratum as a predictor while the best-fitting model excluded it.
Throughout our results, we consider non-overlapping 95% credible intervals of posterior predictions to represent differences in taxa proportional abundances, alpha diversity, or the proportional abundances of Bd-inhibition categories, depending on the analysis. These posterior predictions are all from the best-fitting Bayesian Dirichlet-multinomial regression model (the model with the lowest WAIC) for salamander, water, or substrate samples for bacterial communities, fungal communities, or Bd-inhibition categories. The 95% credible intervals of posterior predictions are presented in the referenced figures.
The following results are from the posterior predictions of the best-fitting Bayesian Dirichlet-multinomial regression model for bacterial community composition on salamander skin. We observed temporal, spatial, and ontogenetic variation in the proportional abundances of the top 100 bacterial taxa (Figure 5; Supplementary Figures S9–S13), with the degree of variation depending on the taxon. Proportional abundance trends were often taxon-specific, although patterns were observed across some taxa. Examples of temporal variation include an increase in the proportional abundance of Comamonadaceae 2 (i.e., the second ASV classified as Comamonadaceae) through time in Gibson Lakes age-0 larvae and a decrease in the proportional abundance of Candidatus Methylopumilus 1 through time in Ponds Lake age-1 larvae. Ontogenetic variation is apparent among many of the top 100 bacterial taxa. For example, the proportional abundances of Comamonadaceae 3 and 6 in Ponds Lake were consistently higher for age-2+ metamorphosed salamanders than other stage classes. In Gibson Lakes age-0 individuals, we observed higher proportional abundances of certain bacterial taxa on larvae followed by a sharp decline post-metamorphosis (e.g., Comamonadaceae 2, Limnohabitans 1, Methylotenera 7, Methylotenera 2, and Sphingorhabdus rigui 1).
Figure 5. Proportional abundance predictions from Bayesian Dirichlet-multinomial regression models for the first ten bacterial taxa from the top 100. Points, thick lines, and thin lines represent posterior medians, 50% credible intervals, and 95% credible intervals, respectively. Hatched ranges represent environmental proportional abundances. Taxa without hatched ranges were not detected in either water or substrate. Note the square root scale on the y-axis. Taxa are ordered (left to right, top to bottom) by descending average proportion of reads in salamander samples in which each combination of site and life stage receives equal weight. Proportional abundance prediction plots for the remaining top 100 bacterial taxa can be found in the Supplementary material.
Ninety of the top 100 bacterial taxa detected in salamander samples were also detected in environmental samples. We detected 86 and 82 of the top 100 bacterial taxa in water and substrate samples, respectively. From the posterior predictions of the best-fitting salamander, water, and substrate Bayesian Dirichlet-multinomial regression models for bacterial community composition, 32 of the top 100 bacterial taxa were disproportionately more abundant on salamander skin relative to the environment for the majority of combinations of stage class and sampling event (Figure 5; Supplementary Figures S9–S13), including Sanguibacter 1 and Gracilibacteria 3. Taxa detected exclusively on salamander skin include Pseudochrobactrum kiredjianiae 1 and Roseomonas 2. None of the top 100 bacterial taxa had lower proportional abundances on salamander skin than in the environment for the majority of combinations of stage class and sampling event.
As for bacteria, we observed spatiotemporal and ontogenetic variation in the proportional abundances of the top 100 fungal taxa based on the posterior predictions of the best-fitting Bayesian Dirichlet-multinomial regression model for fungal community composition on salamander skin (Figure 6; Supplementary Figures S14–S18). Naganishia diffluens experienced an increase in proportional abundance through time for age-0 and age-2+ salamanders at both lakes, and the proportional abundance of Vishniacozyma victoriae increased through time for Gibson Lakes age-0 larvae. Notably, the proportional abundance of Bd in Ponds Lake metamorphosed individuals was very high (between 35 and 55%) compared to Ponds Lake larvae and either life stage at Gibson Lakes (all <2.5%; Figure 6). In contrast, metamorphosed individuals at Ponds Lake had lower proportional abundances of Cystobasidium slooffiae compared to Ponds Lake larvae and either life stage at Gibson Lakes, the opposite of the pattern observed for Bd.
Figure 6. Proportional abundance predictions from Bayesian Dirichlet-multinomial regression models for the first ten fungal taxa from the top 100. Points, thick lines, and thin lines represent posterior medians, 50% credible intervals, and 95% credible intervals, respectively. Hatched ranges represent environmental proportional abundances. Note the square root scale on the y-axis. Taxa are ordered (left to right, top to bottom) by descending average proportion of reads in salamander samples in which each combination of site and life stage receives equal weight. Proportional abundance prediction plots for the remaining top 100 fungal taxa can be found in the Supplementary material.
Ninety-three of the top 100 fungal taxa detected in salamander samples were also detected in environmental samples. We detected 67 and 91 of the top 100 fungal taxa in water and substrate samples, respectively. From the posterior predictions of the best-fitting salamander, water, and substrate Bayesian Dirichlet-multinomial regression models for fungal community composition, 27 of the top 100 fungal taxa were disproportionately more abundant on salamander skin relative to the environment for the majority of combinations of stage class and sampling event (Figure 6; Supplementary Figures S14–S18), including Candida sake, Wallemia muriae, and Vishniacozyma. Taxa detected exclusively on salamander skin include Pleosporales, Melanodiplodia tianschanica, and Buckleyzyma aurantiaca. Four of the top 100 fungal taxa had lower proportional abundances on salamander skin than in the environment for the majority of combinations of stage class and sampling event, including Ascomycota, Basidiomycota, and Rozellomycota.
Salamander bacterial and fungal diversity (i.e., Hill’s diversity with α = 2 derived from the posterior predictions of the best-fitting Bayesian Dirichlet-multinomial regression models for salamander bacterial and fungal communities, respectively) are displayed in Figure 7. Salamander bacterial diversity had the highest values in early-season age-1 metamorphosed individuals at Gibson Lakes, age-2+ neotenes at Ponds Lake in July, and age-0 larvae after mid-August at Gibson Lakes (Figure 7). Bacterial diversity for age-1 larvae increased throughout the early warm season at Ponds Lake when most of this stage class was observed. Bacterial diversity in age-2+ metamorphosed individuals tended to decrease through time at both lakes. In late August and early September, bacterial diversity was higher for age-0 larvae at Gibson Lakes than Ponds Lake. Patterns of microbial diversity for fungi differed than those for bacteria. Salamander fungal diversity was highest among age-0 larvae at Gibson Lakes, age-1 larvae at the beginning of the warm season at Ponds Lake, and late-season larvae at Ponds Lake (Figure 7). Fungal diversity increased throughout the early warm season for metamorphosed individuals at both lakes, with metamorphosed individuals at Ponds Lake having lower fungal diversity than at Gibson Lakes (Figure 7). The lower fungal diversity for metamorphosed individuals at Ponds Lake compared to Gibson Lakes may be due to the high proportional abundance of a single fungal taxon – reducing species evenness – on Ponds Lake metamorphosed salamanders (the posterior medians of this taxon, Bd, ranged from 0.390 to 0.508; Figure 6).
Figure 7. Hill’s diversity index (α = 2) for bacterial and fungal communities on salamander skin. Diversity estimates were derived from proportional abundance predictions of bacterial and fungal taxa from Bayesian Dirichlet-multinomial regression models. Points, thick lines, and thin lines represent posterior medians, 50% credible intervals, and 95% credible intervals, respectively. Note that bacterial diversity for age-0 metamorphosed salamanders at Gibson Lakes on September 15th is omitted because, for this combination of stage class and sampling event only, we found Hill’s diversity to be sensitive to our grouping of rare taxa (i.e., those not belonging to the top 100) into an “other” category in the Bayesian Dirichlet-multinomial regression models.
Based on the best-fitting Bayesian Dirichlet-multinomial regression models for the composition of bacterial Bd-inhibition categories (i.e., Bd-inhibitory, non-Bd-inhibitory, and uncertain Bd-inhibition status), Bd-inhibitory taxa were disproportionately more abundant on salamander skin relative to the environment for most combinations of stage class and sampling event (23 of 25; Figure 8). Non-Bd-inhibitory bacterial taxa were disproportionately more abundant on salamander skin relative to the environment for most combinations of stage class and sampling event at Ponds Lake (10 of 15), but we were unable to detect differences in the proportional abundances of non-Bd-inhibitory bacterial taxa between salamander skin and the environment for any combination of stage class and sampling event at Gibson Lakes (Figure 8). Bacterial taxa of uncertain Bd-inhibition status were disproportionately more abundant in the environment compared to salamander skin for most combinations of stage class and sampling event (23 of 25; Figure 8).
Figure 8. Proportional abundance predictions of Batrachochytrium dendrobatidis (Bd) and Bd-inhibitory bacterial taxa from Bayesian Dirichlet-multinomial regression models. Points, thick lines, and thin lines represent posterior medians, 50% credible intervals, and 95% credible intervals, respectively. Hatched ranges represent environmental proportional abundances. Note the square root scale on the y-axis.
Relationship between Bd-inhibitory bacteria and Bd
From our Bayesian beta-binomial regression model, there is a > 99.9% probability that a negative relationship exists between the proportional abundance of Bd-inhibitory bacteria in bacterial communities and the proportional abundance of Bd in fungal communities on the skin of metamorphosed salamanders (i.e., > 99.9% of MCMC samples for the regression coefficient were negative). The posterior median of the regression coefficient was −2.402, and the 95% credible interval was −3.386 to −1.586. The modeled relationship between the relative abundances of Bd-inhibitory bacteria and Bd on the skin of metamorphosed salamanders is shown in Figure 9.
Figure 9. Bayesian beta-binomial regression between the proportional abundances of Bd-inhibitory bacteria and Bd in bacterial and fungal communities, respectively, on the skin of metamorphosed salamanders. Solid and dashed lines represent the medians and 95% credible intervals of posterior predictions, respectively. Points represent observations from metamorphosed salamanders. Note that there is uncertainty associated with both the response and predictor values (i.e., proportional abundances), and the Bayesian beta-binomial regression considers this uncertainty within the model.
Discussion
We observed spatiotemporal and ontogenetic variation in the relative abundances and microbial diversity of bacterial and fungal taxa in the skin-associated microbiome of the western tiger salamander at two high alpine Rocky Mountain lakes. Our best-fitting Bayesian Dirichlet-multinomial regression models for microbial community composition included all predictors or their interactions except for the models of fungal communities on salamander skin and Bd-inhibition categories in lake water, for which the stratum predictor was excluded (Table 2). Because rare taxa (i.e., not members of the top 100) were grouped into an “other” category in our Bayesian Dirichlet-multinomial regression models, variation in the relative abundances of these rare taxa were masked within changes in the relative abundance of the “other” category. Therefore, our results should be considered conservative estimates of the variation in microbial community composition on tiger salamander skin because variation within the “other” category was not considered. Because our top 100 taxa comprise the vast majority of reads in salamander samples (93.1% of bacterial and 98.6% of fungal reads), we expect variation in the composition of these taxa, plus the “other” category, to represent most of the variation in salamander microbiomes. When viewed conservatively, our models already suggest that all of our covariates (except for stratum in the aforementioned cases) contribute to our ability to predict microbial community composition, so we expect this outcome would be largely unaffected by the inclusion of additional taxa outside of the “other” category.
Our findings of variation in microbial community composition between sites, across life stages, and through time is consistent with the results of other studies (Kueneman et al., 2013; Longo et al., 2015; Bletz et al., 2017a,b). The inclusion of stratum as a predictor in the best-fitting composition model of bacteria on salamander skin, as well as for composition models of bacteria and fungi in the environment, suggests that we observed spatial variation in microbial community composition within lakes in addition to between lakes. Furthermore, the inclusion of salamander age as a predictor in the best-fitting models of bacterial and fungal community composition on salamander skin suggests that, within life stages, we also observed variation in microbial community composition with salamander age. For salamander skin-associated bacterial and fungal communities, composition was better explained by spatiotemporal than water quality covariates. In agreement with other amphibian skin-associated microbiome studies, we found that the skin of the western tiger salamander is a selective environment with taxa disproportionately represented compared to their relative abundances in water and substrate (Kueneman et al., 2013; Walke et al., 2014; Bletz et al., 2017a).
Time or its interactions were included as predictors in all best-fitting Bayesian Dirichlet-multinomial regression models for microbial community composition (Table 2), suggesting that microbial communities changed throughout the warm season. Notably, an interaction between time, salamander age, life stage, and lake was included in the best-fitting model for bacterial communities on salamander skin. Additionally, the best-fitting model for fungal communities on salamander skin included both an interaction between time and salamander age and an interaction between time and lake. This suggests that temporal trends in salamander bacterial and fungal community composition varied by lake and salamander stage class. While time was included as a predictor in the best-fitting model for Bd-inhibition categories on salamander skin, we failed to detect changes in the relative abundances of Bd-inhibitory bacteria through time for any stage class (i.e., within a stage class, all 95% credible intervals overlapped; Figure 8). This concurs with the findings of Bletz et al. (2017a), which found changes in bacterial community composition through time but stability in predicted Bd-inhibitory function. We are hesitant to draw the same conclusion, however, both because failing to detect change does not mean that change has not occurred, and because we could only confidently predict the Bd-inhibition status of a small minority of observed bacterial taxa. We did detect changes in the relative abundances of common bacteria (e.g., Comamonadaceae 2 and Candidatus Methylopumilus 1; Figure 5) and fungi (e.g., Naganishia diffluens and Vishniacozyma victoriae; Figure 6) through time on salamander skin. We also detected changes in bacterial and fungal diversity through time for stage classes at both lakes (Figure 7).
Similar to time, life stage or its interactions were also included as predictors in all best-fitting Bayesian Dirichlet-multinomial regression models for microbial community composition (Table 2). Notably, the proportional abundance of the bacterial taxon Comamonadaceae 1 was higher on the skin of non-metamorphosed salamanders (i.e., larval or neotenic individuals) than metamorphosed salamanders (i.e., 95% credible intervals did not overlap) throughout the warm season at both lakes (Figure 5). This taxon was very abundant on the skin of non-metamorphosed salamanders, typically comprising more than 20% of the bacterial community, and sometimes exceeding 40% (Figure 5). Kueneman et al. (2013) also observed a very high relative abundance (>65%) of a single member of Comamonadaceae on a life stage of the Cascades frog (Rana cascadae), but the taxon dominated the skin of metamorphosed frogs instead of tadpoles. In our study, Comamonadaceae 1 was also disproportionately more abundant on the skin of non-metamorphosed salamanders relative to the environment, whereas we were unable to detect differences in the proportional abundance of this taxon between the environment and the skin of metamorphosed salamanders (Figure 5). While we compared the proportional abundances of microbial taxa on salamander skin to environmental proportional abundances in lake water and lake substrate, we note that metamorphosed salamanders – although caught from the water – may have also had access to terrestrial sources of microbiota (e.g., soil) which we did not sample. For fungi, both the proportional abundance of Cryptococcus uniguttulatus and community diversity were higher on the skin of non-metamorphosed salamanders at every time point where both metamorphosed and non-metamorphosed salamanders were observed (Figures 6, 7). This contrasts with the findings of Kueneman et al. (2016b), in which microeukaryote diversity was higher on adult western toads (Anaxyrus boreas) than tadpoles.
We detected Bd on salamander skin at both lakes, with the relative abundance of Bd being highest for age-1 and age-2+ metamorphosed salamanders at Ponds Lake (Figure 6). The higher abundance of Bd on the skin of metamorphosed compared to larval amphibians is supported by other studies and is thought to be the result of increased keratin, a substrate for Bd, in amphibian skin following metamorphosis, during which structural changes to the skin occur (Berger et al., 1998; Marantelli et al., 2004; Frost et al., 2006). We are unsure why differences in the relative abundance of Bd was much less pronounced between larval and metamorphosed individuals at Gibson Lakes. Since Bd was absent in all negative control samples, we are confident that Bd was present at Gibson Lakes and that its detection was not the result of contamination from Ponds Lake samples.
We observed that Bd-inhibitory bacterial taxa were disproportionately more abundant on salamander skin relative to the environment for most combinations of stage class and sampling event (Figure 8). If bacterial taxa for which we have high confidence in their Bd-inhibition statuses can be considered a random sample from both salamander skin and the environment, then this could be taken as evidence that salamander skin selects for Bd-inhibitory bacteria. However, salamander skin also appeared to select for non-Bd-inhibitory bacteria at one lake, and for both lakes, bacteria of uncertain Bd-inhibition status were disproportionately more abundant in the environment than on salamander skin. Since environmental bacteria are not the focus of the Woodhams et al. (2015) database, we suspect that bacteria in this reference database are more likely to be common on amphibian skin than in the environment. This could result in more environmental bacteria having an uncertain Bd-inhibition status, and Bd-inhibitory and non-Bd-inhibitory bacteria would subsequently appear to be disproportionately more abundant on salamander skin than in the environment. Still, the apparent selection for Bd-inhibitory bacteria on salamander skin is stronger than for non-Bd-inhibitory bacteria (Figure 8), suggesting that selection for Bd-inhibitory bacteria may indeed be occurring. Despite harboring Bd, tiger salamanders have been found to tolerate chytridiomycosis (Davidson et al., 2003), and we suggest that selection for Bd-inhibitory bacteria by tiger salamander skin may contribute to this disease tolerance.
When viewed across combinations of stage class and sampling event, we did not observe any noticeable patterns between the relative abundances of Bd-inhibitory bacteria and Bd (Figure 8). That is, across combinations of stage class and sampling event, the relative abundance of Bd was not high or low when the relative abundance of Bd-inhibitory bacteria was high or low. We did, however, observe a negative pattern between the relative abundances of Bd and the fungal taxon Cystobasidium slooffiae (i.e., the relative abundance of Bd was low when the relative abundance of Cystobasidium slooffiae was high; Figure 6). Conversely, we observed positive patterns between the relative abundances of Bd and Comamonadaceae 3 and 6 (i.e., the relative abundance of Bd was high when the relative abundances of these taxa were high; Figures 5, 6). Comamonadaceae has been found to be abundant on the skin of multiple amphibian species, including the tiger salamander (McKenzie et al., 2011), and some members show evidence of Bd-inhibition or negative co-occurrence with fungal taxa (Woodhams et al., 2015; Kueneman et al., 2016b). Despite this, Walke et al. (2015) found a very weak correlation between a member of Comamonadaceae and Bd, and we found positive patterns between the relative abundances of members of Comamonadaceae and Bd. While we were unable to confidently predict the Bd-inhibition statuses of Comamonadaceae 3 and 6, we did predict one member of Comamonadaceae to be Bd-inhibitory (Comamonadaceae 5). Still, we observed no pattern between the relative abundances of this taxon and Bd (Figures 5, 6).
Within metamorphosed salamanders, we found strong evidence (> 99.9% probability from a Bayesian beta-binomial regression) of a negative relationship between the relative abundances of Bd-inhibitory bacteria and Bd in bacterial and fungal communities, respectively (Figure 9). We caution, however, that the mechanism behind this relationship is unclear from our data. We do not know whether Bd-inhibitory bacteria inhibit Bd growth, or if the opposite is true. Infection with Bd can lead to the restructuring of microbial communities on amphibian skin (Jani and Briggs, 2014; Jani and Briggs, 2018), and it is possible that Bd infection may directly or indirectly inhibit the growth of Bd-inhibitory bacteria. Our use of microbiome read counts to test for a relationship between Bd-inhibitory bacteria and Bd produced comparable results to studies which used quantitative PCR to detect and quantify the abundance of Bd. For example, Jiménez et al. (2022) found that Bd infection intensity significantly decreased on the skin of the eastern newt (Notophthalmus viridescens) as the relative abundance of putative Bd-inhibitory bacteria increased. Similarly, Flechas et al. (2019) found lower Bd infection prevalence within post-metamorphic life stages which also had high relative abundances of Bd-inhibitory bacteria in two frog species.
An analysis between the absolute abundances of Bd-inhibitory bacteria and Bd, instead of the relative abundances, would be of greater interest biologically. Following DNA extraction and prior to PCR, fixed amounts of 16S and ITS synthgenes (i.e., synthetic gene spike-ins) were added to a constant volume of each sample’s DNA extract. The synthgene read counts provide a benchmark to compare taxon read counts with, and can serve as the basis for absolute abundance estimation (Harrison et al., 2021). While we used synthgenes to estimate the amount of microbial DNA in our samples relative to negative controls, we were unable to use the synthgenes for estimating the densities (i.e., count per unit area) of microbial taxa on salamander skin because, as we were not aware of synthgenes at the time, we did not measure swabbed area in the field. Furthermore, a length-weight regression suggested that salamanders grow allometrically (i.e., the body does not grow proportionally in all dimensions; see Supplementary material), so an assumption about salamander growth would have to be made in order to derive surrogates of swabbed area from length measurements (i.e., SVL squared could not be used as a surrogate for swabbed area). We also considered limitations in our swabbing protocol. Since our study focused on variation in microbiome composition, we adopted the swabbing protocol of Bletz et al. (2017a), in which a swab is stroked across the ventral surface of the amphibian ten times (one time = an up and back stroke along the full length of the belly). While swabbing, the ten strokes along the length of the belly were distributed across the belly’s width. Due to the fixed size of the swab, this means that the same belly area was swabbed more times for smaller salamanders than for larger salamanders. This implies that even if we had measured swabbed area, we would have to assume that the number of microbes collected asymptotes after a certain swabbing intensity, and we must have further assumed that we reached this threshold of swabbing intensity. We suggest that the need for these assumptions can be avoided by using a different swabbing protocol. For example, instead of stroking a swab across the ventral surface a certain number of times while covering an area of interest, one could swab the full area of interest (e.g., the belly) a certain number of times and measure the swabbed area, a method which is already applied in studies of Bd load (North and Alford, 2008). If such a swabbing protocol were applied, we suggest that taxon density on amphibian skin could be modeled using a negative binomial regression for rates (i.e., taxon read count per unit “time”) – where rate represents taxon density, and synthgene read count and swabbed area serve as measures of “time.”
A key aim of amphibian skin-associated microbiome studies relates to understanding what role microbial communities play in protecting their hosts against cutaneous diseases such as chytridiomycosis. While DNA metabarcoding is commonly employed to characterize the composition of microbial communities, we experienced challenges relating community composition to functional activity. Using 16S rRNA gene sequences, we were unable to predict Bd-inhibition statuses for the vast majority of our bacterial taxa with any reasonable certainty. This is not surprising given that, after trimming to our amplicon region, the majority of sequences in the Woodhams et al. (2015) database which were shared across multiple bacterial isolates had variable Bd-inhibition statuses, and the isolates included in the database provided limited phylogenetic coverage of our bacterial taxa. Similarly, Becker et al. (2015) found bacterial congeners to frequently range from complete inhibition to facilitation of Bd. Another approach to exploring the functional activity of microbial communities involves metatranscriptomics, the sequencing of RNA within a microbiome to investigate gene expression (Nichols and Davenport, 2021). With a metatranscriptomics approach to exploring functional activity, antifungal secondary metabolite production by microbes experiencing real-world biotic and abiotic conditions on salamander skin could be observed, and a precise knowledge of community composition, while still informative, would not be a pre-requisite for inference.
Our study emphasizes two traditionally understudied areas of amphibian skin-associated microbial ecology, temporal variation in community composition and expanding our view of the microbiome to include fungi in addition to bacteria. Temporal variation in community composition could prove challenging for studies examining spatial variation, where temporal and spatial variation may be confounded. We also identified additional sources of variation in community composition which are not typically considered. Within life stages, we identified additional variation with salamander age, and within lakes, we identified additional variation between strata. Furthermore, we observed that the relationships between community composition and spatiotemporal and stage class covariates are interdependent, complex, and best described using interactions.
Through this study, we have gained a greater understanding of microbial ecology on amphibian skin through the examination of season-long temporal variation of bacterial and fungal communities. In addition to identifying further sources of variation in community composition, we have identified differentially abundant taxa, have examined microbial selection by salamander skin, have investigated alpha diversity, and have tested for a relationship between predicted Bd-inhibitory function and Bd. Ultimately, we hope our findings will assist in the conservation of amphibian species threatened by chytridiomycosis.
Data availability statement
Raw sequence reads are deposited in the Sequence Read Archive (BioProject PRJNA843333). Sample metadata are available on DataDryad (doi: 10.5061/dryad.nzs7h44tw). Scripts used in the analyses are available on GitHub (https://github.com/Urodelan/2022_Salamander_Microbiome).
Ethics statement
The animal study was reviewed and approved by the Utah State University Institutional Animal Care and Use Committee.
Author contributions
KG and ZG designed the research. KG and JH conducted field sampling and wrote the manuscript with revisions from all authors. KG performed the lab work and analyses. ZG provided guidance, support, and feedback at all stages of the project. All authors contributed to the article and approved the submitted version.
Funding
This project was supported with funding from the Utah State University (USU) Quinney College of Natural Resources Undergraduate Research Grant, the USU Office of Research Undergraduate Research and Creative Opportunities Grant, the USU Honors Program Honors Research Fund, the Society for the Study of Amphibians and Reptiles Roger Conant Undergraduate Research Grant in Herpetology, and the National Science Foundation (DEB 1638768 to ZG).
Acknowledgments
We would like to thank our volunteers for fieldwork assistance, with special thanks to Sarah Shaw. We would also like to thank Trisha Atwood, Bonnie Waring, and Kevin Landom for their support and feedback. We thank the Uinta-Wasatch-Cache and Caribou-Targhee National Forests for permission to conduct fieldwork on National Forest lands. The support and resources from the Center for High Performance Computing at the University of Utah are gratefully acknowledged. This study has previously appeared online as a preprint on bioRxiv (doi: 10.1101/2022.06.01.494434).
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.
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.2022.1020329/full#supplementary-material
References
Basanta, M. D., Anaya-Morales, S. L., Martínez-Ugalde, E., González Martínez, T. M., Ávila-Akerberg, V. D., Trejo, M. V., et al. (2022). Metamorphosis and seasonality are major determinants of chytrid infection in a paedomorphic salamander. Anim. Conserv. doi: 10.1111/acv.12824
Becker, M. H., Walke, J. B., Murrill, L., Woodhams, D. C., Reinert, L. K., Rollins-Smith, L. A., et al. (2015). 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
Belasen, A. M., Riolo, M. A., Bletz, M. C., Lyra, M. L., Toledo, L. F., and James, T. Y. (2021). Geography, host genetics, and cross-domain microbial networks structure the skin microbiota of fragmented Brazilian Atlantic Forest frog populations. Ecol. Evol. 11, 9293–9307. doi: 10.1002/ece3.7594
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. 95, 9031–9036. doi: 10.1073/pnas.95.15.9031
Bletz, M. C., Loudon, A. H., Becker, M. H., Bell, S. C., Woodhams, D. C., Minbiole, K. P., et al. (2013). Mitigating amphibian chytridiomycosis with bioaugmentation: characteristics of effective probiotics and strategies for their selection and use. Ecol. Lett. 16, 807–820. doi: 10.1111/ele.12099
Bletz, M. C., Perl, R. G., Bobowski, B. T., Japke, L. M., Tebbe, C. C., Dohrmann, A. B., et al. (2017a). 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
Bletz, M. C., Perl, R. G., and Vences, M. (2017b). Skin microbiota differs drastically between co-occurring frogs and newts. R. Soc. Open Sci. 4:170107. doi: 10.1098/rsos.170107
Bollback, J. P. (2006). SIMMAP: stochastic character mapping of discrete traits on phylogenies. BMC Bioinformatics 7, 1–7. doi: 10.1186/1471-2105-7-88
Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, G. A., et al. (2019). Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857. doi: 10.1038/s41587-019-0209-9
Britton, R. A., and Young, V. B. (2014). Role of the intestinal microbiota in resistance to colonization by Clostridium difficile. Gastroenterology 146, 1547–1553. doi: 10.1053/j.gastro.2014.01.059
Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., Holmes, S. P., et al. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583. doi: 10.1038/nmeth.3869
Caporaso, J. G., Lauber, C. L., Walters, W. A., Berg-Lyons, D., Lozupone, C. A., Turnbaugh, P. J., et al. (2011). Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc. Natl. Acad. Sci. 108, 4516–4522. doi: 10.1073/pnas.1000080107
Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., et al. (2017). Stan: a probabilistic programming language. J. Stat. Softw. 76, 1–32. doi: 10.18637/jss.v076.i01
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
Daskin, J. H., Bell, S. C., Schwarzkopf, L., and Alford, R. A. (2014). Cool temperatures reduce antifungal activity of symbiotic bacteria of threatened amphibians – implications for disease management and patterns of decline. PLoS One 9:e100378. doi: 10.1371/journal.pone.0100378
Davidson, E. W., Parris, M., Collins, J. P., Longcore, J. E., Pessier, A. P., and Brunner, J. (2003). Pathogenicity and transmission of chytridiomycosis in tiger salamanders (Ambystoma tigrinum). Copeia 2003, 601–607. doi: 10.1643/CP-02-120R1
Dethlefsen, L., McFall-Ngai, M., and Relman, D. A. (2007). An ecological and evolutionary perspective on human–microbe mutualism and disease. Nature 449, 811–818. doi: 10.1038/nature06245
Douglas, A. J., Hug, L. A., and Katzenback, B. A. (2021). Composition of the North American wood frog (Rana sylvatica) bacterial skin microbiome and seasonal variation in community structure. Microb. Ecol. 81, 78–92. doi: 10.1007/s00248-020-01550-5
Ellison, S., Rovito, S., Parra-Olea, G., Vásquez-Almazán, C., Flechas, S. V., Bi, K., et al. (2019). The influence of habitat and phylogeny on the skin microbiome of amphibians in Guatemala and Mexico. Microb. Ecol. 78, 257–267. doi: 10.1007/s00248-018-1288-8
Fadrosh, D. W., Ma, B., Gajer, P., Sengamalay, N., Ott, S., Brotman, R. M., et al. (2014). An improved dual-indexing approach for multiplexed 16S rRNA gene sequencing on the Illumina MiSeq platform. Microbiome 2, 1–7. doi: 10.1186/2049-2618-2-6
Familiar López, M., Rebollar, E. A., Harris, R. N., Vredenburg, V. T., and Hero, J. M. (2017). Temporal variation of the skin bacterial community and Batrachochytrium dendrobatidis infection in the terrestrial cryptic frog Philoria loveridgei. Front. Microbiol. 8:2535. doi: 10.3389/fmicb.2017.02535
Flechas, S. V., Acosta-González, A., Escobar, L. A., Kueneman, J. G., Sánchez-Quitian, Z. A., Parra-Giraldo, C. M., et al. (2019). Microbiota and skin defense peptides may facilitate coexistence of two sympatric Andean frog species with a lethal pathogen. ISME J. 13, 361–373. doi: 10.1038/s41396-018-0284-9
Frost, D. R., Grant, T., Faivovich, J., Bain, R. H., Haas, A., Haddad, C. F., et al. (2006). The amphibian tree of life. Bull. Am. Mus. Nat. Hist. 2006, 1–291. doi: 10.5531/sd.sp.13
Gardes, M., and Bruns, T. D. (1993). ITS primers with enhanced specificity for basidiomycetes-application to the identification of mycorrhizae and rusts. Mol. Ecol. 2, 113–118. doi: 10.1111/j.1365-294X.1993.tb00005.x
Gelman, A., Hwang, J., and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Stat. Comput. 24, 997–1016. doi: 10.1007/s11222-013-9416-2
Grant, E. H. C., Adams, M. J., Fisher, R. N., Grear, D. A., Halstead, B. J., and Hossack, B. R. (2018). Identifying management-relevant research priorities for responding to disease-associated amphibian declines. Glob. Ecol. Conserv. 16:e00441. doi: 10.1016/j.gecco.2018.e00441
Grice, E. A., and Segre, J. A. (2011). The skin microbiome. Nat. Rev. Microbiol. 9, 244–253. doi: 10.1038/nrmicro2537
Grunseich, J. M., Thompson, M. N., Aguirre, N. M., and Helms, A. M. (2019). The role of plant-associated microbes in mediating host-plant selection by insect herbivores. Plan. Theory 9:6. doi: 10.3390/plants9010006
Haegeman, B., Hamelin, J., Moriarty, J., Neal, P., Dushoff, J., and Weitz, J. S. (2013). Robust estimation of microbial diversity in theory and in practice. ISME J. 7, 1092–1101. doi: 10.1038/ismej.2013.10
Harris, R. N., Brucker, R. M., Walke, J. B., Becker, M. H., Schwantes, C. R., Flaherty, D. C., et al. (2009). 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
Harrison, J. G., Calder, J. W., Shastry, V., and Buerkle, A. C. (2020). Dirichlet-multinomial modelling outperforms alternatives for analysis of microbiome and other ecological count data. Mol. Ecol. Resour. 20, 481–497. doi: 10.1111/1755-0998.13128
Harrison, J. G., Calder, J. W., Shuman, B., and Buerkle, A. C. (2021). The quest for absolute abundance: the use of internal standards for DNA-based community ecology. Mol. Ecol. Resour. 21, 30–43. doi: 10.1111/1755-0998.13247
Hawkins, J. A., Jones, S. K., Finkelstein, I. J., and Press, W. H. (2018). Indel-correcting DNA barcodes for high-throughput sequencing. Proc. Natl. Acad. Sci. 115, E6217–E6226. doi: 10.1073/pnas.1802640115
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. 111, E5049–E5058. doi: 10.1073/pnas.1412752111
Jani, A. J., and Briggs, C. J. (2018). Host and aquatic environment shape the amphibian skin microbiome but effects on downstream resistance to the pathogen Batrachochytrium dendrobatidis are variable. Front. Microbiol. 9:487. doi: 10.3389/fmicb.2018.00487
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
Kruger, A. (2020). Frog skin microbiota vary with host species and environment but not chytrid infection. Front. Microbiol. 11:1330. doi: 10.3389/fmicb.2020.01330
Kueneman, J. G., Parfrey, L. W., Woodhams, D. C., Archer, H. M., Knight, R., and McKenzie, V. J. (2013). 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., Weiss, S., and McKenzie, V. J. (2017). Composition of micro-eukaryotes on the skin of the cascades frog (Rana cascadae) and patterns of correlation between skin microbes and Batrachochytrium dendrobatidis. Front. Microbiol. 8:2350. doi: 10.3389/fmicb.2017.02350
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. R. Soc. B Biol. Sci. 283:20161553. doi: 10.1098/rspb.2016.1553
Kueneman, J. G., Woodhams, D. C., Van Treuren, W., Archer, H. M., Knight, R., and McKenzie, V. J. (2016b). Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas). ISME J. 10, 934–944. doi: 10.1038/ismej.2015.168
Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). lmerTest package: tests in linear mixed effects models. R package version 3.1.3. J. Stat. Softw. 82, 1–26. doi: 10.18637/jss.v082.i13
Langille, M. G., Zaneveld, J., Caporaso, J. G., McDonald, D., Knights, D., Reyes, J. A., et al. (2013). Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat. Biotechnol. 31, 814–821. doi: 10.1038/nbt.2676
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, Y. K., and Mazmanian, S. K. (2010). Has the microbiota played a critical role in the evolution of the adaptive immune system? Science 330, 1768–1773. doi: 10.1126/science.1195568
Letunic, I., and Bork, P. (2021). Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 49, W293–W296. doi: 10.1093/nar/gkab301
Longcore, J. E., Pessier, A. P., and Nichols, D. K. (1999). Batrachochytrium dendrobatidis gen. et sp. nov., a chytrid pathogenic to amphibians. Mycologia 91, 219–227. doi: 10.1080/00275514.1999.12061011
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
Marantelli, G., Berger, L., Speare, R., and Keegan, L. (2004). Distribution of the amphibian chytrid Batrachochytrium dendrobatidis and keratin during tadpole development. Pac. Conserv. Biol. 10, 173–179. doi: 10.1071/PC040173
Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17, 10–12. doi: 10.14806/ej.17.1.200
McFall-Ngai, M., Hadfield, M. G., Bosch, T. C., Carey, H. V., Domazet-Lošo, T., Douglas, A. E., et al. (2013). Animals in a bacterial world, a new imperative for the life sciences. Proc. Natl. Acad. Sci. 110, 3229–3236. doi: 10.1073/pnas.1218525110
McKenzie, V. J., Bowers, R. M., Fierer, N., Knight, R., and Lauber, C. L. (2011). Co-habiting amphibian species harbor unique skin bacterial communities in wild populations. ISME J. 6, 588–596. doi: 10.1038/ismej.2011.129
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
Muletz-Wolz, C. R., Fleischer, R. C., and Lips, K. R. (2019). Fungal disease and temperature alter skin microbiome structure in an experimental salamander system. Mol. Ecol. 28, 2917–2931. doi: 10.1111/mec.15122
Nichols, R. G., and Davenport, E. R. (2021). The relationship between the gut microbiome and host gene expression: a review. Hum. Genet. 140, 747–760. doi: 10.1007/s00439-020-02237-0
Nilsson, R. H., Larsson, K. H., Taylor, A. F. S., Bengtsson-Palme, J., Jeppesen, T. S., Schigel, D., et al. (2019). The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classifications. Nucleic Acids Res. 47, D259–D264. doi: 10.1093/nar/gky1022
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
Parada, A. E., Needham, D. M., and Fuhrman, J. A. (2016). Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ. Microbiol. 18, 1403–1414. doi: 10.1111/1462-2920.13023
Plummer, M. (2003). JAGS: A Program for Analysis of Bayesian Graphical Models Using Gibbs Sampling. Proceedings of the 3rd International Workshop on Distributed Statistical Computing, Vienna. (Vol. 124, No. 125.10, pp. 1–10).
Plummer, M. (2019). rjags: Bayesian graphical models using MCMC. R package version 4-10. Available at: https://CRAN.R-project.org/package=rjags.
Price, M. N., Dehal, P. S., and Arkin, A. P. (2010). FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS One 5:e9490. doi: 10.1371/journal.pone.0009490
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 (2020). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
Revell, L. J. (2012). phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223. doi: 10.1111/j.2041-210X.2011.00169.x
Rollins-Smith, L. A. (1998). Metamorphosis and the amphibian immune system. Immunol. Rev. 166, 221–230. doi: 10.1111/j.1600-065X.1998.tb01265.x
Sabino-Pinto, J., Galán, P., Rodríguez, S., Bletz, M. C., Bhuju, S., Geffers, R., et al. (2017). Temporal changes in cutaneous bacterial communities of terrestrial- and aquatic-phase newts (Amphibia). Environ. Microbiol. 19, 3025–3038. doi: 10.1111/1462-2920.13762
Savage, A. E., Sredl, M. J., and Zamudio, K. R. (2011). Disease dynamics vary spatially and temporally in a North American amphibian. Biol. Conserv. 144, 1910–1915. doi: 10.1016/j.biocon.2011.03.018
Sennhenn-Reulen, H. (2018). Bayesian regression for a Dirichlet distributed response using Stan. arXiv [Preprint]. doi: 10.48550/arXiv.1808.06399
Sexton, O. J., and Bizer, J. R. (1978). Life history patterns of Ambystoma tigrinum in montane Colorado. Am. Midl. Nat. 99:101. doi: 10.2307/2424936
Sharon, G., Garg, N., Debelius, J., Knight, R., Dorrestein, P. C., and Mazmanian, S. K. (2014). Specialized metabolites from the microbiome in health and disease. Cell Metab. 20, 719–730. doi: 10.1016/j.cmet.2014.10.016
Sievers, F., Wilm, A., Dineen, D., Gibson, T. J., Karplus, K., Li, W., et al. (2011). Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol. Syst. Biol. 7:539. doi: 10.1038/msb.2011.75
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. Eco Health 4, 125–134. doi: 10.1007/s10393-007-0093-5
Stan Development Team (2020). RStan: the R interface to Stan. R package version 2.21.2. Available at: http://mc-stan.org/
Stebbins, R. C. (2003). A Field Guide to Western Reptiles and Amphibians, 3rd Edn. New York: Houghton Mifflin Harcourt.
Tosh, P. K., and McDonald, L. C. (2012). Infection control in the multidrug-resistant era: tending the human microbiome. Clin. Infect. Dis. 54, 707–713. doi: 10.1093/cid/cir899
Tourlousse, D. M., Yoshiike, S., Ohashi, A., Matsukura, S., Noda, N., and Sekiguchi, Y. (2017). Synthetic spike-in standards for high-throughput 16S rRNA gene amplicon sequencing. Nucleic Acids Res. 45:e23. doi: 10.1093/nar/gkw984
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., Becker, M. H., Loftus, S. C., House, L. L., Teotonio, T. L., Minbiole, K. P., et al. (2015). Community structure and function of amphibian skin microbes: an experiment with bullfrogs exposed to a chytrid fungus. PLoS One 10:e0139848. doi: 10.1371/journal.pone.0139848
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). Naïve 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
White, T. J., Bruns, T., Lee, S., and Taylor, J. (1990). “Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics” in PCR Protocols: A Guide to Methods and Applications. eds. M. A. Innis, D. H. Gelfand, J. J. Sninsky, and T. J. White (New York, NY: Academic Press, Inc)
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
Keywords: temporal variation, spatial variation, bacteria, fungi, relative abundance, Bayesian analysis
Citation: Goodwin KB, Hutchinson JD and Gompert Z (2022) Spatiotemporal and ontogenetic variation, microbial selection, and predicted Bd-inhibitory function in the skin-associated microbiome of a Rocky Mountain amphibian. Front. Microbiol. 13:1020329. doi: 10.3389/fmicb.2022.1020329
Edited by:
Eria Alaide Rebollar, National Autonomous University of Mexico, MexicoReviewed by:
Maria Delia Basanta, University of Nevada, United StatesObed Hernandez-Gomez, Dominican University of California, United States
Copyright © 2022 Goodwin, Hutchinson and Gompert. 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: Kenen B. Goodwin, dXJvZGVsYW5AZ21haWwuY29t