Environmental and Host Effects on Skin Bacterial Community Composition in Panamanian Frogs

Research on the amphibian skin microbiota has focused on identifying bacterial taxa that deter a pathogenic chytrid fungus, and on describing patterns of microbiota variation. However, it remains unclear how environmental variation affects amphibian skin bacterial communities, and whether the overall functional diversity of the amphibian skin microbiota is associated to such variation. We sampled skin microbial communities from one dendrobatoid frog species across an environmental gradient along the Panama Canal, and from three dendrobatoid frog species before and after the onset of the wet season in one site. We found frog skin microbial alpha diversity to be highest in frogs from sites with low soil pH, but no clear effect of the onset of the wet season. However, we found frog skin microbial community structure to be affected by soil pH and the onset of the wet season, which also resulted in a decrease in between-sample variation. Across the sampled frog species, bacterial functional groups changed with the onset of the wet season, with certain bacterial functional groups entirely disappearing and others differing in their relative abundances. In particular, we found the proportion of Bd-inhibitory bacteria to correlate with mean soil pH, and to increase in two of the frog species with the onset of the wet season. Taken together, our results suggest that structure and predicted function of amphibian bacterial skin communities may be influenced by environmental variables such as pH and precipitation, site effects, and host effects.

Introduction

The interactions between vertebrate hosts and their microbial communities have recently become a subject of accrued interest (Caporaso et al., 2010b, 2012; Escalona et al., 2016; Rebollar et al., 2016a), because microbiota may play an important role in animal ecology and evolution. For example, microbiota have been found to alter the behavior of their host (Bravo et al., 2011; Ezenwa et al., 2012; Marin et al., 2017), affect digestion (Schnorr et al., 2014) and development (McFall-Ngai et al., 2013), as well as contribute to immune system function (Kueneman et al., 2016a; Woodhams et al., 2016). However, how wildlife-associated microbiota vary according to abiotic factors remains poorly understood (Jiménez and Sommer, 2017; Medina et al., 2017; Pollock et al., 2017).

In amphibians, the skin microbiota are known to vary between host species (Kueneman et al., 2014; Belden et al., 2015; Pollock et al., 2017) and sampling sites (Kueneman et al., 2014; Krynak et al., 2016). It has also been suggested that ontogenetic and seasonal changes contribute to variation in the frog skin microbiota (Longo et al., 2015). Interestingly, microbes abundant on amphibians’ skin are usually rare or present in low abundances in the environment, suggesting that the amphibian skin represents a unique and selective environment (Loudon et al., 2014; Walke et al., 2014; Rebollar et al., 2016b). However, it is also clear that amphibians depend on microbial environmental reservoirs, such as forest soil, to maintain diverse skin microbiota (Fitzpatrick and Allison, 2014; Loudon et al., 2014; Kueneman et al., 2016a).

Both moisture (Hartmann et al., 2014) and pH (Fierer and Jackson, 2006; Lauber et al., 2009; Zhalnina et al., 2015; Schappe et al., 2017) drive patterns of variation in soil microbial communities, but whether these abiotic factors indirectly shape the amphibian skin microbial communities remains unclear. Alternatively, soil pH and moisture could potentially have an effect on microbe–microbe and microbe–host interactions, for example, by affecting amphibian antimicrobial peptide secretions (Rinaldi, 2002; Tennessen et al., 2009). Even though many frog skin microbiota studies were focused on identifying sources of variation, whether changes in skin microbial communities correspond with changes in overall microbial functional groups remains poorly understood.

In amphibians, the skin microbiota appear to constitute a critical component of the immune system (Walke et al., 2015; Woodhams et al., 2015; Kueneman et al., 2016b), whereby certain bacterial taxa are capable of producing anti-fungal metabolites (Brucker et al., 2008; Lauer et al., 2008; Harris et al., 2009b). Some of these anti-fungal metabolites deter the pathogenic chytrid fungus Batrachochytrium dendrobatidis (Bd), which has been associated with certain amphibian declines (Lips et al., 2006; Pounds et al., 2006). One of these agents, the betaproteobacterium Janthinobacterium lividum has been isolated from a number of amphibians (e.g., Plethodon cinereus, Brucker et al., 2008; Becker and Harris, 2010; Gastrotheca excubitor, Burkart et al., 2017), and can reduce Bd-related mortality (Harris et al., 2009a; Becker and Harris, 2010).

To date, only a few studies have empirically assessed variation in the overall microbial functional groups in amphibian hosts (Loudon et al., 2014; Davis et al., 2017). However, ecologically important the anti-Bd function of the skin microbes may be, this constitutes only one of the many functions performed by the amphibian skin microbiota. To have a better understanding of the functional implications of amphibian skin microbiota disruption, for example due to chytridiomycosis or environmental variation, it is necessary to better characterize the composition of these microbiota as a first step in assessing the broad spectrum of functionality that they may represent for the amphibian host. Accordingly, we sought to investigate the possible effects of abiotic environmental variables (i.e., soil pH and precipitation levels) across four sites in one frog species, as well as seasonal and host effects in three species of dendrobatoid frogs within one site in shaping the diversity and predicted function of the skin microbiota. Sites along the Panama Canal inhabited by these frogs differ in soil pH and annual precipitation, with well-defined wet and dry seasons, within a relative small area. Only about 50 km separate the Pacific from the Atlantic termini of the Panama Canal, along which is a steep climatic gradient in precipitation, with average annual precipitation ranging from ca. 1,500 to over 4,000 mm per year (Santiago and Mulkey, 2005; Engelbrecht et al., 2007).

If amphibians require microbial environmental reservoirs, such as forest soil, to maintain diverse skin bacterial communities (Loudon et al., 2014), then environmental variables that drive variation in free-living microbial communities like (1) soil pH, (2) precipitation levels, and (3) the onset of the wet season should correlate with frog skin microbial community structure. Frog skin bacterial diversity should be expected to be highest in frogs sampled from sites with neutral soil pH and high humidity, where soil bacterial diversity is highest (Fierer and Jackson, 2006; Lauber et al., 2009; Hartmann et al., 2014; Zhalnina et al., 2015), in comparison to either pH extreme and dryer sites. We also explored (4) host-specific differences on the skin microbiota. Finally, (5) if soil pH, precipitation levels, the onset of the wet season, and host-related skin microbial community shifts are ecologically relevant to these three dendrobatoid frog species, then we should observe a relationship between them and predicted microbial functional group variation, including the presence or absence of bacteria that may either inhibit or facilitate the growth of Bd.

Materials and Methods

Sampling Sites and Frog Species

In total, we sampled 70 dendrobatoid frogs (Table 1) from four sites along the Panama Canal in April 2016 (Figure 1). Cerro Ancon is a secondary growth urban forest, whereas the other three sites are characterized by old growth forests (Ibáñez et al., 2002). For each site, soil pH was recorded from soil samples (2–10 cm from surface) mixed in water (1:2, soil:water ratio; Table 2). We obtained historical annual mean and daily precipitation records from the Panama Canal Authority weather station closest to each site (20–136 years of records), and we used the daily values to calculate accumulated precipitation 5, 10, and 15 days before sampling. We used these accumulated precipitation values to confirm the beginning of the wet season. We believe that these time intervals could explain microbiota variation, as previous studies have recorded changes in amphibian skin microbiota in as little as 7 days in captivity (Loudon et al., 2014) and 3 days in field experiments (Longo and Zamudio, 2016).

TABLE 1

TABLE 1. Host species habitat and sample sizes of the three dendrobatoid frogs by life stage.

FIGURE 1

FIGURE 1. Sample sites along the Panama Canal in Central Panama. (A) Central Panama showing location of the Panama Canal region (inset). (B) Panama Canal region indicating sample localities.

TABLE 2

TABLE 2. Abiotic factors for the sites sampled along the Panama Canal: soil pH and precipitation levels.

To test the effects of soil pH and precipitation levels on skin microbiota, we focused our sampling on juvenile and adult Dendrobates auratus (N = 9–19 per site, four sites), a frog living in the leaf litter where soil pH is likely affecting its microbiota (dataset “Dend.aura”). We included D. auratus life stage as a fixed factor in our models to account for previously reported ontogenetic differences (Longo et al., 2015). To test seasonal effects on skin microbiota in one site (Pipeline Road), we sampled three dendrobatoid frog species before and after the onset of the wet season (D. auratus N_before = 8, N_after = 3; Silverstoneia flotator N_before = 2, N_after = 4; and Allobates talamancae N_before = 2, N_after = 4, dataset “Pipeline”). To test differences between frog species, we used the “Pipeline” dataset and controlled for sampling time (i.e., “before” or “after” the onset of the wet season).

To collect our samples, we caught frogs with gloved hands and rinsed them with 50 ml of distilled autoclaved water to wash off any transient microbes (Loudon et al., 2014; Rebollar et al., 2016b). After rinsing, we swabbed the frogs with two sterile cotton tipped swabs 10 times on the ventral side, 10 times on each leg and each toe once. We used new gloves for each frog to prevent cross-contaminating our samples. We stored the swabs in sterile, autoclaved micro-centrifuge tubes on ice until they could be stored in a -20°C freezer prior to DNA extraction. We used one swab for DNA extraction to test Bd presence and to construct 16S rRNA libraries. We stored the second swab at -80°C as a reference.

Research permits allowing the collection of our samples were granted by the Panamanian authority Ministerio de Ambiente (Permit No. SE/A-38-16) and by the Authority of the Panama Canal (Permit No. 173). All procedures with animals were conducted under approved animal care protocols from the Smithsonian Tropical Research Institute’s Animal Care Committee (2016-0301-2019-2) and by the McGill University’s Animal Care Committee (2000-4569).

DNA Extraction, Bd Analyses, and 16S rRNA Library Preparation

We treated swabs with a lysozyme buffer to break up Gram-positive bacteria’s cell walls, followed by a modified protocol of the QIAGEN DNeasy Blood & Tissue Kit.

To test for Bd presence, we performed the qPCR assay developed by Boyle et al. (2004) with the following modifications. We extracted DNA from swabs using the QIAGEN DNeasy Blood & Tissue Kit (Kosch and Summers, 2013; Bletz et al., 2015). We performed the qPCR reactions on a Roche LightCycler 96 System, and a volume of 20 μl per well (Kriger et al., 2006a). Instead of using Applied Biosystems 2× TaqMan Master Mix, we used Roche FastStart Essential DNA Probes Master and we added 0.25 μl of Roche LightCycler Uracil-DNA Glycosylase per reaction to eliminate potential PCR carryover contamination. We used Applied Biosystems TaqMan Exogenous Internal Positive Control Reagents VIC Probe (Hyatt et al., 2007), and we included five negative and two Bd-positive controls on every plate (Kriger et al., 2006a,b). We performed a duplicate analysis by initially running a 1:10 dilution of the DNA extract from each sample in singlicate, followed by a second run of the undiluted DNA extract in singlicate to avoid false negatives due to a potentially low and undetectable quantity of DNA in the diluted samples. Samples in which Bd was not detected in both runs of the duplicate analysis were considered negative. Despite the fact that the sites we sampled have previously been reported as Bd endemic (Woodhams et al., 2008; Rebollar et al., 2014), and Bd might still be present in other amphibian species, all our samples tested negative for Bd; therefore, we did not include Bd as a factor in our statistical analyses.

To construct our 16S rRNA library, we amplified each sample following Caporaso et al. (2012), using primers F515 and R806. We pretreated the PCR pools with Just-a-Plate^TM PCR 96 Purification and Normalization Kit (Charm Biotech, San Diego, CA, United States). We used equimolar parts of each sample to construct the library and cleaned it with Serapure beads to remove any primer dimers. DNA concentration and quality of the library were assessed with NanoDrop^TM ND-1000 (Thermo Fisher Scientific Inc., 2008, Wilmington, DE, United States), Invitrogen Qubit^® Fluorometer (Thermo Fisher Scientific Inc., 2007, Wilmington, DE, United States), and Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., 2000–2005, Santa Clara, CA, United States). The library was sequenced at the NAOS laboratories of the Smithsonian Tropical Research Institute in Panama using an Illumina MiSeq v3 600 cycles cartridge.

Sequence Reads Processing

We processed all samples in the same manner before separating them into the two data sets. Using Quantitative Insights into Microbial Ecology (QIIME, Caporaso et al., 2010b), we filtered sequences to retain high-quality reads. We clustered sequences according to 97% sequence similarity threshold using the UCLUST method (Edgar, 2010). The resulting operational taxonomic units (OTUs) were classified using the Greengenes database when possible (May 2013 release; McDonald et al., 2012); otherwise they were clustered de novo using a 97% sequence similarity threshold. We assigned taxonomy with the Greengenes database and the RDP classifier (Wang et al., 2007). We used PyNAST (Caporaso et al., 2010a) to align representative sequences to the Greengenes database and constructed a phylogenetic tree with FastTree2 (Price et al., 2010). Lastly, we filtered the resulting OTU table to only include clusters ≥0.001% of the total reads (Bokulich et al., 2013). To standardize sequencing effort, we rarefied the samples to a depth of 6,813 sequences. This initial dataset included 27,120 bacterial OTUs from 70 samples. As demonstrated by plateauing rarefaction curves, the samples had appropriate sequencing coverage, which indicates that our sampling captured most of the bacterial alpha diversity (Supplementary Figure S1). For downstream analyses, we divided the initial dataset into a subset of 22,846 OTUs from 58 samples of D. auratus from the four sites (“Dend.aura”), and a subset of 12,354 OTUs from 23 samples from all three species from Pipeline Road (“Pipeline”). Sequences for each sample were deposited in the NCBI SRA under BioProject ID PRJNA433445.

Comparing Alpha and Beta Diversity

Unless otherwise stated, all statistical and graphical analyses were performed using R 3.3.2 (R Core Team, 2017, Vienna, Austria) and the ggplot2 package (Wickham, 2009). Full test statistics can be found in Supplementary Table S1.

To estimate alpha diversity of skin microbiota per frog species and site, we calculated Shannon’s diversity index, which measures the number of species and their evenness within a sample. Using the lme4 package (Bates et al., 2015) and the “Dend.aura” dataset, we tested for an effect of mean soil pH and precipitation levels on bacterial alpha diversity, including site as a random factor and life stage as a fixed factor. We tested for an effect of the onset of the wet season on alpha diversity using unpaired t-tests (two-tailed) for each of the three host species from Pipeline Road (“Pipeline”) sampled before and after the onset of the wet season. We also tested for host effects on bacterial alpha diversity among the Pipeline Road sample using linear mixed effects model with sampling time as a random factor (i.e., “before” or “after” the onset of the wet season).

We used the data derived from D. auratus (“Dend.aura”) samples to test the effect of mean soil pH and precipitation levels on bacterial beta diversity, and data from Pipeline Road (“Pipeline”) samples of all three host species to test the effect of the onset of the wet season and host species effects on bacterial beta diversity. We calculated the UniFrac Weighted and UnWeighted (Lozupone and Knight, 2005) and Bray–Curtis distances matrices, and obtained NMDS ordinations. To determine differences between groups, we used the adonis function within the vegan package (Oksanen et al., 2016) and 999 permutations. We controlled for site effects by including site as strata, which does not correspond directly to a random factor, but in this case strata controlled for repeated measures within sites.

To determine which OTUs were driving the reported beta diversity differences, we performed linear discriminant analyses of effect sizes using LEfSe and normalizing the sum of values to 1M as recommended by Segata et al. (2011). We used 0.05 thresholds for both the Kruskal–Wallis (among class differences) and Wilcoxon (between classes differences) tests, and we ran a linear discriminant analysis (LDA) considering LDA scores >3 as ecologically relevant. LDA scores are an estimate of the effect size of features that are differentially abundant among the categories tested (Segata et al., 2011). For bacterial OTUs beta diversity and functional group richness analyses, we tested whether groups differed in their dispersion patterns by running the betadisper function from the vegan package (Oksanen et al., 2016), which is equivalent to a Levene’s test comparing homogeneity of variances between groups. We only reported significant betadisper results.

Comparing Bacterial Diversity and Predicted Functional Group Variation

To assess whether the alpha and beta diversity differences associated with sites that differed in soil pH and precipitation levels, timing in relation to the onset of the wet season, and host species corresponded with changes in predicted microbial functional groups, we compared our OTU table against the Functional Annotation of Prokaryotic Taxa (FAPROTAX) database (Louca et al., 2016). FAPROTAX extrapolates functions of cultured prokaryotes (identified at the genus or species level) to the rest of the prokaryotic genus to estimate putative function. There are two main limitations of applying this approach to our data: (1) the FAPROTAX database was constructed mainly to analyze biogeochemistry of water bodies and (2) the FAPROTAX database is non-exhaustive; thus, only a small percentage of our OTUs may be assigned to at least one functional group. Even considering these caveats, we believe that predicting microbial functional groups using FAPROTAX is an appropriate alternative, especially when a metatranscriptomics approach is not plausible. Additionally, because FAPROTAX only considers cultured bacterial species with characterized functions, it may be superior to genomic prediction approaches like PICRUSt (Langille et al., 2013) and PAPRICA (Bowman and Ducklow, 2015).

Estimating the Bd-Inhibitory and Bd-Enhancing Functions of the Skin Microbiota

Finally, to identify bacterial isolates that potentially inhibit or enhance Bd, we used the pick_closed_reference_otus.py script in QIIME (Caporaso et al., 2010b) and standard settings to compare our results to the database published by Woodhams et al. (2015). It should be noted that (a) the database only includes culturable bacteria and thus is not exhaustive, (b) the database is likely biased toward bacterial isolates that inhibit Bd, (c) some of the bacteria were isolated from non-Panamanian frogs (Woodhams et al., 2015), and (d) microbiota functions may be context dependent (Medina et al., 2017). To standardize the resulting OTU table, we calculated proportion of Bd-inhibitory and Bd-enhancing bacterial sequences per sampled frog. To do this, we considered each frog and computed the sum of the relative abundances of all Bd-inhibitory and Bd-enhancing bacteria separately. We compared these proportions between host species, sites, and in relation to the onset of the wet season using generalized linear models with the function glm (R Core Team, 2017). Because we were comparing proportions bound between 0 and 1 and to account for overdispersion, we used quasibinomial distributions in the case of Bd-inhibitory bacteria. Due to their low proportions, we used quasipoisson distributions to analyze Bd-enhancing bacteria.

Results

Among-Site Skin Microbiota Variation in D. auratus

In D. auratus, Shannon diversity correlated negatively with mean soil pH and was highest in frogs sampled from the site with lowest soil pH (Figure 2B, P < 0.05). However, precipitation levels did not have an effect on Shannon diversity (P > 0.05; nor on observed OTUs or Phylogenetic Diversity, Supplementary Table S1). In terms of bacterial community structure, we found the sites with the lowest mean soil pH (Fort Sherman) and highest mean soil pH (Barro Colorado Island) to plot on opposite sides of the ordination, with both intermediate mean soil pH sites (Cerro Ancon and Pipeline Road) in the middle. The pattern is less clear with regards to precipitation levels. The two sites with highest precipitation (Fort Sherman and Barro Colorado Island) plotted on opposite sides of the ordination, whereas the two intermediate precipitation level sites (Cerro Ancon and Pipeline Road) plotted in the middle of the ordination. Nevertheless, we found soil pH and precipitation levels as statistically significant factors explaining variation in bacterial community structure (Figure 2C, PERMANOVA UniFrac Weighted, Bray–Curtis, and UniFrac UnWeighted, P < 0.05). Out of 2,084 OTUs, we found 33 most significant OTUs best explaining site differences in D. auratus skin bacterial beta diversity (LDA > 3, Figure 2A). In particular, we found clear overrepresentation of Pseudomonadaceae, Oxalobacteraceae, and Xanthomonadaceae in Fort Sherman; Sphingomonadaceae and Weeksellaceae in Cerro Ancon; Brevibacteriaceae in Pipeline Road; and Moraxellaceae and Sphingobacteriaceae in Barro Colorado.

FIGURE 2

FIGURE 2. Skin microbiota in Dendrobates auratus in relation to sample site and soil pH. (A) Bacterial operational taxonomic units (OTUs) that best define bacterial community structure in D. auratus across sites. Bars represent most significant linear discriminant analysis (LDA) scores computed using LEfSe, i.e., LDA > 3. (B) Box-and-whisker plot comparison of skin bacterial Shannon alpha diversity in D. auratus across sites in relation to average soil pH. Box indicates median, first, and third quartiles; points extending past 1.5 times the inter-quartile length were considered outliers and depicted as dots. (C) Non-metric multi-dimensional scaling (NMDS) plot of Bray–Curtis distances in skin bacterial beta diversity between individuals of D. auratus across sites. In all charts, sites are color-coded as follows: Fort Sherman is red, Cerro Ancon is yellow, Pipeline Road in green, and Barro Colorado is dark blue.

At higher bacterial taxonomic levels, though, the bacterial skin communities among D. auratus populations appeared to be similar. To illustrate, OTUs belonging to the Gammaproteobacteria were abundant in all sites, whereas Betaproteobacteria were abundant in all sites except in Pipeline Road. Anaerolineae bacteria were abundant in Fort Sherman and Barro Colorado Island, the two sites with the lowest and highest mean soil pH, respectively, among the four sites. Bacilli were abundant in Cerro Ancon and in Barro Colorado Island, which are the sites with the most human traffic (i.e., tourists and residents in Cerro Ancon, and researchers and visitors in Barro Colorado Island). The rest of the overrepresentations at higher taxonomic levels were site-specific. In Fort Sherman, Actinobacteria and Alphaproteobacteria were abundant; whereas Deltaproteobacteria, Acidobacteria, Flavobacteriia, and Pseudanabaenales were abundant in Cerro Ancon. At this taxonomic resolution, we did not detect any site-specific overrepresentations in Pipeline Road. Gemmatimonadetes and Thaumarchaeota were abundant in Barro Colorado Island.

Within-Site Skin Microbiota Variation in Relation to the Wet Season and Host Species

We found no effect of the onset of the wet season on skin bacterial alpha diversity of Pipeline Road D. auratus or A. talamancae (both P > 0.05), but there was a barely significant decline in alpha diversity of S. flotator (Figure 3B; t = 2.77, df = 4, P = 0.05). Within Pipeline Road, we found frog skin bacterial beta diversity to be assorted into distinct clusters in relation to the onset of the wet season (Figure 3C; PERMANOVA UniFrac Weighted and UniFrac UnWeighted P < 0.001; Bray–Curtis P < 0.01). After the onset of the wet season, the frog skin microbiota became more homogeneous, which was illustrated by the wet season samples forming a tighter cluster than the dry season samples (Figure 3C; betadisper ANOVA on UniFrac UnWeighted distances F_(1,21) = 6.791, P < 0.05). Out of 1,679 OTUs, we found 20 most significant OTUs explaining the differences in relation to the onset of the wet season (LDA > 3, Figure 3A). Actinobacteria OTUs were abundant both before and after the onset of the wet season. However, the abundant Actinobacteria before the onset of the wet season was a Micrococcaceae, whereas after the onset of the wet season the abundant OTU belonged to the Pseudonocardiaceae. The frogs sampled before the onset of the wet season were characterized by an abundance of OTUs belonging to Pedosphaerae, Synechococcophycidae, Clostridia, Alphaproteobacteria, Betaproteobacteria, and Brachyspirae. On the other hand, Gammaproteobacteria, Thermophilia, Flavobacteriia, and Cyanobacteria were abundant after the onset of the wet season. Overall, some bacterial classes appeared to be more abundant before the onset of the wet season (e.g., Actinobacteria and Bacilli), and decreased after the onset of the wet season when other classes became relatively more abundant (e.g., Flavobacteria and Sphingobacteria, Supplementary Figure S2).

FIGURE 3

FIGURE 3. Skin microbiota in three dendrobatoid frog species from Pipeline Road in relation to the onset of the wet season. (A) Bacterial OTUs that best define bacterial community structure among species before and after the onset of the wet season. Bars represent most significant linear discriminant analysis (LDA) scores computed using LEfSe, i.e., LDA > 3. (B) Box-and-whisker plot comparison of skin bacterial Shannon alpha diversity in dendrobatoid frog species before and after the onset of the wet season. Box indicates median, first, and third quartiles; points extending past 1.5 times the inter-quartile length were considered outliers and depicted as dots. ^∗ indicates a significant difference in alpha diversity before and after the onset of the wet season. (C) NMDS plot of Bray–Curtis distances in skin bacterial beta diversity between individual frogs of three dendrobatoid species at Pipeline Road. In all charts, brown indicates before the onset of the wet season, whereas light blue indicates after the onset of the set season.

We found Shannon diversity to differ between frog species in Pipeline Road (P < 0.05). However, these three species did not differ in observed OTUs or Phylogenetic Diversity (Supplementary Table S1). Out of the three frog species sampled, D. auratus had the lowest alpha diversity indexes, followed by A. talamancae and then S. flotator. Across frog species, microbial communities were dominated by Gammaproteobacteria, Alphaproteobacteria, and Betaproteobacteria (Supplementary Figure S2). The two stream-associated frogs, A. talamancae and S. flotator, also appeared to have higher relative abundances of Sphingobacteria than D. auratus. Flavobacteria and Actinobacteria seemed to be more abundant in D. auratus. After controlling for sampling time effects, we found no host effect on beta diversity (PERMANOVA Bray–Curtis, UniFrac Weighted and UnWeighted P > 0.05).

Skin Microbiota Predicted Function Variation

We could assign 489 out of 2,184 bacterial OTUs (22.39%) to at least one microbial functional group using the FAPROTAX database, and most functional groups were represented by more than one OTU (Table 3). D. auratus skin microbiota had different microbial functional groups at sites differing in soil pH and precipitation levels (PERMANOVA Bray–Curtis P < 0.001). These site-to-site differences were best explained by 17 functional groups (Table 3), including fermentation in Fort Sherman, aerobic chemoheterotrophy in Cerro Ancon, anoxygenic photoautotrophy in Pipeline Road, and human pathogens in Barro Colorado Island.

TABLE 3

TABLE 3. Bacterial functional groups that best defined the site differences in the microbiota of Dendrobates auratus frogs between sampling sites, the number of bacterial operational taxonomic units (OTUs) identified for each functional group, and the relative abundance of bacteria in each functional group, presented in terms of linear discriminant analysis (LDA) scores computed using LEfSe.

At Pipeline Road and after controlling for frog species, skin microbial functional groups were different before and after the onset of the wet season (PERMANOVA Bray–Curtis P < 0.01), with 16 functional groups associated with this difference (Figure 4). These differences were explained by certain functional groups changing in their relative abundance or entirely disappearing after the onset of the wet season. To elaborate, functional groups involving nitrite, nitrate, and photoautotrophy were more abundant before the onset of the wet season, whereas chemoheterotrophy, cellulolysis, and nitrogenous compounds respiration were more abundant after the onset of the wet season.

FIGURE 4

FIGURE 4. Mean relative abundances of skin microbial functional groups before and after the onset of the wet season in three species of dendrobatoid frogs, Dendrobates auratus, Silverstoneia flotator, and Allobates talamancae, from Pipeline Road, Panama. Relative abundances are depicted in terms of color intensity from white (=0) to darkest blue (=0.433). Microbial functional groups more abundant before the onset of the wet season are preceded by a brown bar, whereas groups more abundant after the onset of the wet season are preceded by a blue bar. Asterisks (^∗) indicate microbial functional groups more abundant in D. auratus.

Variation of Potential Bd-Interacting Microbiota

We could assign 185 out of 2,184 OTUs (8.47%) to one of the Bd-interacting categories using the Woodhams et al. (2015) database. The majority of these OTUs were isolated from Panamanian frogs (118/185 or 63.38%), and 14.05% (26/185) were isolated from the same species we report in this study (3 from A. talamancae, 6 from D. auratus, and 17 from S. flotator). D. auratus sampled from different sites differed in the proportion of Bd-inhibitory (Figure 5A, ${\mathrm{\chi }}_{\mathrm{57,3}}^{\mathrm{2}}$ = 6.144, P < 0.010) and Bd-enhancing microbiota (Figure 5B, ${\mathrm{\chi }}_{\mathrm{57,3}}^{\mathrm{2}}$ = 0.196, P < 0.05); yet, the two Bd-interacting microbial categories followed opposite trends. Bd-inhibitory bacterial OTUs were relatively more abundant in D. auratus from Pipeline Road, followed by Cerro Ancon, Fort Sherman, and Barro Colorado, whereas Bd-enhancing bacterial OTUs were relatively more abundant in Fort Sherman, followed by Cerro Ancon, Pipeline Road, and Barro Colorado Island.

FIGURE 5

FIGURE 5. Proportions of Batrachochytrium dendrobatidis (Bd) inhibitory (A) and enhancing (B) microbiota in D. auratus across sites in central Panama, arranged in order of increasing average soil pH values (A,B), and in relation to the onset of the wet season in three Panamanian dendrobatoid frog species at Pipeline Road (C,D).

In Pipeline Road, the proportion of Bd-inhibitory bacteria differed in relation to the onset of the wet season (${\mathrm{\chi }}_{\mathrm{22,1}}^{\mathrm{2}}$ = 2.420, P < 0.05) and between frog hosts (Figure 5C, ${\mathrm{\chi }}_{\mathrm{20,2}}^{\mathrm{2}}$ = 2.420, P < 0.05). Even though Bd-enhancing bacteria appeared to be relatively more abundant before the onset of the wet season across the three dendrobatoid species, this difference fell short from statistical significance (Figure 5D, P > 0.05), and we found no discernible difference in the proportion of Bd-enhancing bacteria among frog species (P > 0.05).

Discussion

In this study, we assessed the role of two abiotic factors (i.e., soil pH and precipitation levels), seasonal, and host effects in shaping the amphibian skin microbiota diversity and predicted functional groups. To address our hypotheses, we sampled D. auratus across four sites differing in soil pH and precipitation levels, and three dendrobatoid frog species before and after the onset of the wet season within one site.

Our results demonstrate that variation in the skin microbiota of the Panamanian dendrobatoid frogs, D. auratus, S. flotator, and A. talamancae, from the Panama Canal region correlates with environmental variation in soil pH and changes in rainfall associated with the onset of the wet season. This is consistent with the hypothesis that these frogs depend on microbial environmental reservoirs to maintain their diverse skin microbiota (Fitzpatrick and Allison, 2014; Loudon et al., 2014; Kueneman et al., 2016a; Rebollar et al., 2016b). Thus, inhabiting an environment with conditions that promote higher microbial alpha diversity could translate into the frog’s skin being colonized by more microbial species. Typically, soil bacterial alpha diversity is expected to correlate with soil pH following a quadratic model, with highest diversity in nearly neutral soils with pH ≈ 6.0 than soils with pH <4.5 or >8 (Fierer and Jackson, 2006; Lauber et al., 2009; Zhalnina et al., 2015). However, our findings that D. auratus from the lowest soil pH site, Fort Sherman, had the highest alpha diversity and that frogs from Fort Sherman and Barro Colorado, which had the most extremely acidic and basic soil pHs, respectively, among our sample sites had the most differentiated bacterial communities, suggest that the frog host may nevertheless influence the composition of its skin microbiota. This finding is in line with evidence that amphibian skin microbiota can be composed disproportionately of microbes that are otherwise present in relatively low abundances in the environment (Walke et al., 2014; Rebollar et al., 2016b).

Although precipitation levels have no apparent effect on skin microbiota alpha diversity in D. auratus, they do appear to correlate with skin bacterial community dispersion patterns, as frogs sampled from sites with wetter climates showed lower between sample variation than frogs sampled from sites with drier climates. This is in line with the thought that moisture promotes more diverse environmental microbial communities (Lauber et al., 2009; Hartmann et al., 2014). The site differences may be explained both in terms of the relative abundances and the presence (or absence) of particular OTUs. Three of the bacterial classes that we found to be overrepresented, Actinobacteria, Gammaproteobacteria, and Betaproteobacteria, were previously reported to be more abundant in frogs from Bd-positive sites, whereas Sphingobacteria were abundant in frogs from a Bd-negative site (Rebollar et al., 2016b).

A dependence on microbial environmental reservoirs to maintain frog skin bacterial communities is further indicated by the significant change in skin bacterial alpha diversity, particularly in S. flotator at Pipeline Road, coincident with the onset of the wet season rains. The streamside habitat of these frogs tends to dry almost completely during the dry season but floods during the wet season, which may result in an increase in stress hormones, which has been related to a reduction in alpha diversity in the oral microbiota of certain species (Stothart et al., 2016). Our results also indicate that skin microbial communities of different species of frogs can become more similar with the onset of the wet season, which may be due to the more homogeneously damp environment after the start of the rains, allowing the frogs to forage more widely. Such seasonal changes in microbial community structure have also been reported in the skin microbiota of temperate ranid frogs (Lithobates yavapaiensis, Longo et al., 2015), and in the gut microbiota of other vertebrate groups (Keenan et al., 2013). Moreover, whether these seasonal patterns result in changes in the overall functional groups of the frog skin microbiota remained largely unexplored.

The redundancy that we found in microbial functional groups, whereby similar functions may be attributed to more than one bacterial OTU, may corroborate with evidence of an amphibian skin functional core (Hamady and Knight, 2009; Human Microbiome Project Consortium, 2012; Davis et al., 2017). From an evolutionary ecology perspective, such redundancy of microbial functional groups may be expected to be beneficial to the host in the face of potentially rapid microbial turnover (Loudon et al., 2014; Hughey et al., 2016), as an essential metabolic function, if carried out by only a single OTU could easily be lost, to the detriment of the host’s fitness. Such redundancy of ecologically relevant microbial functional groups has also been observed in soil (Wellington et al., 2003). In our study, the most abundant microbial functional groups were chloroplasts in Fort Sherman, aerobic chemoheterotrophs in Cerro Ancon, anoxygenic photoautotrophs in Pipeline Road, and human pathogens in Barro Colorado Island. The various bacteria representing these functional groups can exploit different carbon and energy sources (Balkwill et al., 1989; Johnston et al., 2009), which might indicate the position they occupy in the microbial biofilm and whether they interact predominantly with the frog host or with other microbes. The overabundance of human pathogens in frogs from Barro Colorado Island may also suggest that amphibians can play a role in the transmission of zoonotic diseases (Daszak et al., 2000) or, reciprocally, that human activities can affect the microbial communities of wildlife (Hacioglu and Tosunoglu, 2014). The shift in the most abundant functional groups from before to after the onset of the wet season from bacteria associated with photoautotrophy and nitrification to those associated to chemoheterotrophy and aerobic chemoheterotrophy suggests that when the environmental context of the host becomes more homogeneous, the skin microbiota follow a similar trend. That the microbial functional group differences we observe are ecologically relevant to the frog host remain, however, an assumption. How these functional groups actually interact with the host remains largely unknown.

The skin microbiota of these frogs may also be shaped by the properties of the frog’s skin, as well as the environment. All amphibians produce and sequester a rich variety of nitrogenous compounds, including alkaloids, in skin granular glands that could represent a strong microbial selective force. Indeed, alkaloids purified from anuran skin extracts have been found to act as antimicrobials in in vitro trials (Mina et al., 2015). Like other species of dart-poison frogs, D. auratus sequesters dietary-derived alkaloids in its skin to a far greater extent than either S. flotator or A. talamancae. The relative high abundances of aromatic compound degrading bacteria, such as Rhodococcus, Nocardioides and Acinetobacter, that we found in the skin microbiota of D. auratus compared to those of S. flotator and A. talamancae suggest that certain microbes might interact with the frog’s toxic alkaloids. This is similar to the case of the Midwife Toad, Alytes obstetricans, which secretes potent skin peptides that are also thought to play a role in shaping its microbiota (Davis et al., 2017). The potential for interactions between biogenically active amphibian skin compounds and amphibian skin microbiota regarding the frog’s toxicity and its immunity from disease may have considerable ecological and evolutionary implications.

Even though all the frogs we sampled tested negative for Bd, this chytrid amphibian pathogen is known to be present on both sides of the Panama Canal (Woodhams et al., 2008; Rebollar et al., 2014) and our results suggest that environmental factors may mediate the abundance of Bd-inhibitory and Bd-enhancing microbiota. Variation in the proportion of Bd-inhibitory bacteria in the skin of D. auratus correlated with mean soil pH, and decreasing proportions of Bd-enhancing bacteria in the skin of all three frog species in concert with the onset of the wet season, as we have found, may help explain seasonal variation observed in Bd prevalence (Lips et al., 2006; Longo et al., 2010; Whitfield et al., 2012). The lower proportions of Bd-inhibitory bacteria and higher proportions of Bd-enhancing bacteria we observed in the microbiota of A. talamancae compared to D. auratus, may signify the presence of differential effects on population persistence in these frogs (Lam et al., 2010). These findings may offer an explanation for previous suggestions that A. talamancae is more susceptible to Bd than D. auratus (Lips et al., 2006), and that susceptible frog species have distinct skin microbiota (Rebollar et al., 2016b).

Both the FAPROTAX database (Louca et al., 2016) and the Bd-interacting bacteria database (Woodhams et al., 2015) that we used to identify frog skin microbial functions have their limitations. As they are not exhaustive, they may be able to assign only a small percentage of discovered bacterial OTUs to one or more functional groups, enabling us to predict the function of less that 25% of the frog skin bacterial community. Furthermore, a particular bacterial OTU might carry out distinct functions depending on environmental context (Medina et al., 2017) and thus our predictions could be underestimating bacterial functional diversity. However, it is unlikely that we reported a bacterial functional group that was actually absent as most of these were largely redundant. Nevertheless, because both of these databases only consider cultured bacterial species with characterized functions, they may be considered superior to genomic approaches, like PICRUSt (Langille et al., 2013) or PAPRICA (Bowman and Ducklow, 2015), in minimizing Type I errors in predicting microbial functional groups.

We note that as both soil pH and precipitation levels are associated to the sampling sites, to be cautious we could interpret the differences discussed here as site differences. Similarly, given our reduced sample sizes and lack of seasonal replicates, we should be cautious when generalizing the results associated with the onset of the wet season. Despite these limitations, we present evidence that abiotic factors, site effects, and host effects shape bacterial skin community structure and predicted function in three tropical frog species.

The variability we observe in the skin microbiota of these three species of Panamanian frogs appears clearly to be correlated with abiotic factors, even though there may not yet be sufficient evidence to assign cause and effect. Frog skin microbial communities may plausibly be affected by interactions of constituent microbes with the environment, the host, and/or other co-occurring microbial species. Thus, there remains much more to be understood about such wildlife-associated microbiota and the complexities of their interactions with their amphibian hosts.

Author Contributions

BV designed the study, collected the data, performed the analyses, and wrote the first version of the manuscript. DL, RI, and DG revised the manuscript. BV and DG prepared the figures. All authors read and approved the final manuscript.

Funding

This project was supported by the National Science and Engineering Research Council of Canada Scholarships: CGS-M-481960-2015 and CSM-MSFSS-488268 to BV, and Discovery Grants 312126-2012 to DL and 106154-2013 to DG. RI was supported by the Panama Amphibian Rescue and Conservation Project and the Sistema Nacional de Investigación of Panama.

Conflict of Interest Statement

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.

Acknowledgments

We thank Marta Vargas for the technical guidance for constructing 16S rRNA libraries, Estefany Illueca for the technical laboratory assistance for conducting the qPCR assays, and Dr. Ben Turner for kindly providing the soil pH values. We also thank Dr. James Barnett, Catherine Fauvelle, and Colleen Bobbie for the comments on the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2018.00298/full#supplementary-material

References

Balkwill, D. L., Fredrickson, J. K., and Thomas, J. M. (1989). Vertical and horizontal variations in the physiological diversity of the aerobic chemoheterotrophic bacterial microflora in deep southeast coastal plain subsurface sediments. Appl. Environ. Microbiol. 55, 1058–1065.

PubMed Abstract | Google Scholar

Bates, D., Maechler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48. doi: 10.18637/jss.v067.i01

CrossRef Full Text | Google Scholar

Becker, M., and Harris, R. N. (2010). Cutaneous bacteria of the redback salamander prevent morbidity associated with a lethal disease. PLoS One 5:e10957. doi: 10.1371/journal.pone.0010957

PubMed Abstract | CrossRef Full Text | Google Scholar

Belden, L. K., Hughey, M. C., Rebollar, E. A., Umile, T. P., Loftus, S. C., Burzynski, E. A., et al. (2015). Panamanian frog species host unique skin bacterial communities. Front. Microbiol. 6:1171. doi: 10.3389/fmicb.2015.01171

PubMed Abstract | CrossRef Full Text | Google Scholar

Bletz, M. C., Rebollar, E. A., and Harris, R. N. (2015). Differential efficiency among DNA extraction methods influences detection of the amphibian pathogen Batrachochytrium dendrobatidis. Dis. Aquat. Organ. 113, 1–8. doi: 10.3354/dao02822

PubMed Abstract | CrossRef Full Text | Google Scholar

Bokulich, N. A., Subramanian, S., Faith, J. J., Gevers, D., Gordon, J. I., Knight, R., et al. (2013). Quality-filtering vastly improves diversity estimates from illumina amplicon sequencing. Nat. Methods 10, 57–59. doi: 10.1038/nmeth.2276

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowman, J. S., and Ducklow, H. W. (2015). Microbial communities can be described by metabolic structure: a general framework and application to a seasonally variable, depth-stratified microbial community from the coastal west Antarctic Peninsula. PLoS One 10:e0135868. doi: 10.1371/journal.pone.0135868

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, D. G., Boyle, D. B., Olsen, V., Morgan, J. A., and Hyatt, A. D. (2004). Rapid quantitative detection of chytridiomycosis (Batrachochytrium dendrobatidis) in amphibian samples using real-time Taqman PCR assay. Dis. Aquat. Organ. 60, 141–148. doi: 10.3354/dao060141

PubMed Abstract | CrossRef Full Text | Google Scholar

Bravo, J. A., Forsythe, P., Chew, M. V., Escaravage, E., Savignac, H. M., Dinan, T. G., et al. (2011). Ingestion of Lactobacillus strain regulates emotional behavior and central GABA receptor expression in a mouse via the vagus nerve. Proc. Natl. Acad. Sci. U.S.A. 108, 16050–16055. doi: 10.1073/pnas.1102999108

PubMed Abstract | CrossRef Full Text | Google Scholar

Brucker, R. M., Harris, R. N., Schwantes, C. R., Gallaher, T. N., Flaherty, D. C., Lam, B. A., et al. (2008). Amphibian chemical defense: Antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus. J. Chem. Ecol. 34, 1422–1429. doi: 10.1007/s10886-008-9555-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Burkart, D., Flechas, S. V., Vredenburg, V. T., and Catenazzi, A. (2017). Cutaneous bacteria, but not peptodes, are associated with chytridiomycosis resistance in Peruvian marsupial frogs. Anim. Conserv. 20, 483–491. doi: 10.1111/acv.12352

CrossRef Full Text | Google Scholar

Caporaso, J. G., Bittinger, K., Bushman, F. D., DeSantis, T. Z., Andersen, G. L., and Knight, R. (2010a). PyNAST: a flexible tool for aligning sequences to a template alignment. Bioinfomatics 26, 266–267. doi: 10.1093/bioinformatics/btp636

PubMed Abstract | CrossRef Full Text | Google Scholar

Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010b). Qiime allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336. doi: 10.1038/nmeth.f.303

PubMed Abstract | CrossRef Full Text | Google Scholar

Caporaso, J. G., Lauber, C. L., Walters, W. A., Berg-Lyons, D., Huntley, J., Fierer, N., et al. (2012). Ultra-high-throughput microbial community analysis on the illumina hiseq and miseq platforms. ISME J. 6, 1621–1624. doi: 10.1038/ismej.2012.8

PubMed Abstract | CrossRef Full Text | Google Scholar

Daszak, P., Cunningham, A. A., and Hyatt, A. D. (2000). Emerging infectious diseases of wildlife–threats to biodiversity and human health. Science 287, 443–449. doi: 10.1126/science.287.5452.443

CrossRef Full Text | Google Scholar

Davis, L. R., Bigler, L., and Woodhams, D. C. (2017). Developmental trajectories of amphibian microbiota: response to bacterial therapy depends on initial community structure. Environ. Microbiol. 19, 1502–1517. doi: 10.1111/1462-2920.13707

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R. C. (2010). Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461. doi: 10.1093/bioinformatics/btq461

PubMed Abstract | CrossRef Full Text | Google Scholar

Engelbrecht, B. M. J., Comita, L. S., Condit, R., Kursar, T. A., Tyree, M. T., Turner, B. L., et al. (2007). Drought sensitivity shapes species distribution patterns in tropical forests. Nature 447, 80–82. doi: 10.1038/nature05747

PubMed Abstract | CrossRef Full Text | Google Scholar

Escalona, M., Rocha, S., and Posada, D. (2016). A comparison of tools for the simulation of genomic next-generation sequencing data. Nat. Rev. Genet. 17, 459–469. doi: 10.1038/nrg.2016.57

PubMed Abstract | CrossRef Full Text | Google Scholar

Ezenwa, V. O., Gerardo, N. M., Inouye, D. W., Medina, M., and Xavier, J. B. (2012). Microbiology. Animal behavior and the microbiome. Science 338, 198–199. doi: 10.1126/science.1227412

PubMed Abstract | CrossRef Full Text | Google Scholar

Fierer, N., and Jackson, R. B. (2006). The diversity and biogeography of soil bacterial communities. Proc. Natl. Acad. Sci. U.S.A. 103, 626–631. doi: 10.1073/pnas.0507535103

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitzpatrick, B. M., and Allison, A. L. (2014). Similarity and differentiation between bacteria associated with skin of salamanders (Plethodon jordani) and free-living assemblages. FEMS Microbiol. Ecol. 88, 482–494. doi: 10.1111/1574-6941.12314

PubMed Abstract | CrossRef Full Text | Google Scholar

Hacioglu, N., and Tosunoglu, M. (2014). Determination of antimicrobial and heavy metal resistance profiles of some bacteria isolated from aquatic amphibian and reptile species. Environ. Monit. Assess. 186, 407–413. doi: 10.1007/s10661-013-3385-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamady, M., and Knight, R. (2009). Microbial community profiling for human microbiome projects: tools, techniques, and challenges. Genome Res. 19, 1141–1152. doi: 10.1101/gr.085464.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, R. N., Brucker, R. M., Walke, J. B., Becker, M. H., Schwantes, C. R., Flaherty, D. C., et al. (2009a). Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824. doi: 10.1038/ismej.2009.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, R. N., Lauer, A., Simon, M. A., Banning, J. L., and Alford, R. A. (2009b). Addition of antifungal skin bacteria to salamanders ameliorates the effects of chytridiomycosis. Dis. Aquat. Organ. 83, 11–16. doi: 10.3354/dao02004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartmann, M., Niklaus, P. A., Zimmermann, S., Schmutz, S., Kremer, J., Abarenkov, K., et al. (2014). Resistance and resilience of the forest soil microbiome to logging-associated compaction. ISME J. 8, 226–244. doi: 10.1038/ismej.2013.141

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughey, M. C., Walke, J. B., Becker, M. H., Umile, T. P., Burzynski, E. A., Minbiole, K. P., et al. (2016). Short-term exposure to coal combustion waste has little impact on the skin microbiome of adult spring peepers (Pseudacris crucifer). Appl. Environ. Microbiol. 82, 3493–3502. doi: 10.1128/AEM.00045-16

PubMed Abstract | CrossRef Full Text | Google Scholar

Human Microbiome Project Consortium (2012). Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214. doi: 10.1038/nature11234

PubMed Abstract | CrossRef Full Text | Google Scholar

Hyatt, A. D., Boyle, A. H. D., Olsen, V., Boyle, D. B., Berger, L., Obendorf, D., et al. (2007). Diagnostic assays and sampling protocols for the detection of Batrachochytrium dendrobatidis. Dis. Aquat. Organ. 73, 175–192. doi: 10.3354/dao073175

PubMed Abstract | CrossRef Full Text | Google Scholar

Ibáñez, R., Condit, R., Angehr, G., Aguilar, S., García, T., Martínez, R., et al. (2002). An ecosystem report on the panama canal: monitoring the status of the forest communities and the watershed. Environ. Monit. Assess. 80, 65–95. doi: 10.1023/A:1020378926399

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiménez, R. R., and Sommer, S. (2017). The amphibian microbiome: natural range of variation, pathogenic dysbiosis, and role in conservation. Biodivers. Conserv. 26, 763–786. doi: 10.1007/s10531-016-1272-x

CrossRef Full Text | Google Scholar

Johnston, D. T., Wolfe-Simon, F., Pearson, A., and Knoll, A. H. (2009). Anoxygenic photosynthesis modulated proterozoic oxygen and sustained earth’s middle age. Proc. Natl. Acad. Sci. U.S.A. 106, 16925–16929. doi: 10.1073/pnas.0909248106

PubMed Abstract | CrossRef Full Text | Google Scholar

Keenan, S. W., Engel, A. S., and Elsey, R. M. (2013). The alligator gut microbiome and implications for archosaur symbioses. Sci. Rep. 3:2877. doi: 10.1038/srep02877

PubMed Abstract | CrossRef Full Text | Google Scholar

Kosch, T. A., and Summers, K. (2013). Techniques for minimizing the effects of PCR inhibitors in the chytridiomycosis assay. Mol. Ecol. Resour. 13, 230–236. doi: 10.1111/1755-0998.12041

PubMed Abstract | CrossRef Full Text | Google Scholar

Kriger, K. M., Hero, J. M., and Ashton, K. J. (2006a). Cost efficiency in the detection of chytridiomycosis using PCR assay. Dis. Aquat. Organ. 71, 149–154.

PubMed Abstract | Google Scholar

Kriger, K. M., Hines, H. B., Hyatt, A. D., Boyle, D. G., and Hero, J. M. (2006b). Techniques for detecting chytridiomycosis in wild frogs: comparing histology with real-time Taqman PCR. Dis. Aquat. Organ. 71, 141–148.

PubMed Abstract | Google Scholar

Krynak, K. L., Burke, D. J., and Benard, M. F. (2016). Landscape and water characteristics correlate with immune defense traits across blanchard’s cricket frog (Acris blanchardi) populations. Biol. Conserv. 193, 153–167. doi: 10.1016/j.biocon.2015.11.019

CrossRef Full Text | Google Scholar

Kueneman, J. G., Parfrey, L. W., Woodhams, D. C., Archer, H. M., Knight, R., and Mckenzie, V. J. (2014). The amphibian skin-associated microbiome across species, space and life history stages. Mol. Ecol. 23, 1238–1250. doi: 10.1111/mec.12510

PubMed Abstract | CrossRef Full Text | Google Scholar

Kueneman, J. G., Woodhams, D. C., Harris, R., Archer, H. M., Knight, R., and Mckenzie, V. J. (2016a). Probiotic treatment restores protection against lethal fungal infection lost during amphibian captivity. Proc. Biol. Sci. 283:20161553. doi: 10.1098/rspb.2016.1553

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, B. A., Walke, J. B., Vredenburg, V. T., and Harris, R. N. (2010). Proportion of individuals with anti-Batrachochytrium dendrobatidis skin bacteria is associated with population persistence in Rana muscosa. Biol. Conserv. 143, 529–531. doi: 10.1016/j.biocon.2009.11.015

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Lauber, C. L., Hamady, M., Knight, R., and Fierer, N. (2009). Pyrosequencing-based assessment of soil ph as a predictor of soil bacterial community structure at the continental scale. Appl. Environ. Microbiol. 75, 5111–5120. doi: 10.1128/AEM.00335-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Lauer, A., Simon, M. A., Banning, J. L., Lam, B. A., and Harris, R. N. (2008). Diversity of cutaneous bacteria with antifungal activity isolated from female four-toed salamanders. ISME J. 2, 145–157. doi: 10.1038/ismej.2007.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Lips, K. R., Brem, F., Brenes, R., Reeve, J. D., Alford, R. A., Voyles, J., et al. (2006). Emerging infectious disease and the loss of biodiversity in a neotropical amphibian community. Proc. Natl. Acad. Sci. U.S.A. 103, 3165–3170. doi: 10.1073/pnas.0506889103

PubMed Abstract | CrossRef Full Text | Google Scholar

Longo, A. V., Burrowes, P. A., and Joglar, R. L. (2010). Seasonality of Batrachochytrium dendrobatidis infection in direct-developing frogs suggests a mechanism for persistence. Dis. Aquat. Organ. 92, 253–260. doi: 10.3354/dao02054

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Longo, A. V., and Zamudio, K. R. (2016). Environmental fluctuations and host skin bacteria shift survival advantage between frogs and their fungal pathogen. ISME 11, 349–361. doi: 10.1038/isme.2016.138

PubMed Abstract | CrossRef Full Text | Google Scholar

Louca, S., Parfrey, L. W., and Doebeli, M. (2016). Decoupling function and taxonomy in the global ocean microbiome. Science 353, 1272–1277. doi: 10.1126/science.aaf4507

PubMed Abstract | CrossRef Full Text | Google Scholar

Loudon, A. H., Woodhams, D. C., Parfrey, L. W., Archer, H., Knight, R., Mckenzie, V., et al. (2014). Microbial community dynamics and effect of environmental microbial reservoirs on red-backed salamanders (Plethodon cinereus). ISME J. 8, 830–840. doi: 10.1038/ismej.2013.200

PubMed Abstract | CrossRef Full Text | Google Scholar

Lozupone, C., and Knight, R. (2005). Unifrac: a new phylogenetic method for comparing microbial communities. Appl. Environ. Microbiol. 71, 8228–8235. doi: 10.1128/AEM.71.12.8228-8235.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Marin, I. A., Goertz, J. E., Ren, T., Rich, S. S., Onengut-Gumuscu, S., Farber, E., et al. (2017). Microbiota alteration is associated with the development of stress-induced despair behavior. Sci. Rep. 7:43859. doi: 10.1038/srep43859

PubMed Abstract | CrossRef Full Text | Google Scholar

McDonald, D., Price, M. N., Goodrich, J., Nawrocki, E. P., DeSantis, T. Z., Probst, A., et al. (2012). An improved greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea. ISME J. 6, 610–618. doi: 10.1038/ismej.2011.139

PubMed Abstract | CrossRef Full Text | Google Scholar

McFall-Ngai, M., Hadfield, M. G., Bosch, T. C., Carey, H. V., Domazet-Loso, T., Douglas, A. E., et al. (2013). Animals in a bacterial world, a new imperative for the life sciences. Proc. Natl. Acad. Sci. U.S.A. 110, 3229–3236. doi: 10.1073/pnas.1218525110

PubMed Abstract | CrossRef Full Text | Google Scholar

Medina, D., Hughey, M. C., Becker, M. H., Walke, J. B., Umile, T. P., Burzynski, E. A., et al. (2017). Variation in metabolite profiles of amphibian skin bacterial communities across elevations in the Neotropics. Microb. Ecol. 74, 227–238. doi: 10.1007/s00248-017-0933-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Mina, A. E., Ponti, A. K., Woodcraft, N. L., Johnson, E. E., and Saporito, R. A. (2015). Variation in alkaloid-based microbial defenses of the dendrobatid poison frog Oophaga pumilio. Chemoecology 25, 169–178. doi: 10.1007/s10886-018-0930-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Oksanen, J., Blancher, F. G., Friednly, M., Kindt, R., Legendre, P., Mcglinn, D., et al. (2016). Vegan: Community Ecology Package. R package Version 2.4-1. Available at: https://CRAN.R-project.org.package=vegan

Google Scholar

Pollock, F. J., Wada, N., Torda, G., Willis, B. L., and Bourne, D. G. (2017). White syndrome-affected corals have a distinct microbiome at disease lesion fronts. Appl. Environ. Microbiol. 83:e02799-16. doi: 10.1128/AEM.02799-16

PubMed Abstract | CrossRef Full Text | Google Scholar

Pounds, J. A., Bustamante, M. R., Coloma, L. A., Consuegra, J. A., Fogden, M. P., Foster, P. N., et al. (2006). Widespread amphibian extinctions from epidemic disease driven by global warming. Nature 439, 161–167. doi: 10.1038/nature04246

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, M. N., Dehal, P. S., and Arkin, A. P. (2010). FastTree2 – approximately maximum-likelihood trees for large alignments. PLoS One 5:e9490. doi: 10.1371/journal.pone.0009490

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2017). R: a Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Google Scholar

Rebollar, E. A., Antwis, R. E., Becker, M. H., Belden, L. K., Bletz, M. C., Brucker, R. M., et al. (2016a). Using “omics” and integrated multi-omics approaches to guide probiotic selection to mitigate chytridiomycosis and other emerging infectious diseases. Front. Microbiol. 7:68. doi: 10.3389/fmicb.2016.00068

CrossRef Full Text | Google Scholar

Rebollar, E. A., Hughey, M. C., Harris, R. N., Domangue, R. J., Medina, D., Ibáñez, R., et al. (2014). The lethal fungus Batrachochytrium dendrobatidis is present in lowland tropical forests of far eastern Panamá. PLoS One 9:e95484. doi: 10.1371/journal.pone.0095484

PubMed Abstract | CrossRef Full Text | Google Scholar

Rebollar, E. A., Hughey, M. C., Medina, D., Harris, R. N., Ibáñez, R., and Belden, L. K. (2016b). Skin bacterial diversity of panamanian frogs is associated with host susceptibility and presence of Batrachochytrium dendrobatidis. ISME J. 10, 1682–1695. doi: 10.1038/ismej.2015.234

PubMed Abstract | CrossRef Full Text | Google Scholar

Rinaldi, A. C. (2002). Antimicrobial peptides from amphibian skin: an expanding scenario. Curr. Opin. Chem. Biol. 6, 799–804. doi: 10.1016/S1367-5931(02)00401-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Santiago, L. S., and Mulkey, S. S. (2005). Leaf productivity along a precipitation gradient in lowland Panama: patterns from leaf to ecosystem. Trees 19, 349–356. doi: 10.1007/s00468-004-0389-9

CrossRef Full Text | Google Scholar

Schappe, T., Albornoz, F. E., Tuner, B. L., Neat, A., Condit, R., and Jones, F. A. (2017). The role of soil chemistry and plant neighbourhoods in structuring fungal communities in three Panamanian forests. J. Ecol. 105, 569–579. doi: 10.1111/1365-2745.12752

CrossRef Full Text | Google Scholar

Schnorr, S. L., Candela, M., Rampelli, S., Centanni, M., Consolandi, C., Basaglia, G., et al. (2014). Gut microbiome of the hadza hunter-gatherers. Nat. Commun. 5:3654. doi: 10.1038/ncomms4654

PubMed Abstract | CrossRef Full Text | Google Scholar

Segata, N., Izard, J., Waldron, L., Gevers, D., Miropolsky, L., Garrett, W. S., et al. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12:R60. doi: 10.1186/gb-2011-12-6-r60

PubMed Abstract | CrossRef Full Text | Google Scholar

Stothart, M. R., Bobbie, C. B., Schulte-Hostedde, A. I., Boonstra, R., Palme, R., Mykytczuk, N. C., et al. (2016). Stress and the microbiome: linking glucocorticoids to bacterial community dynamics in wild red squirrels. Biol. Lett. 12:20150875. doi: 10.1098/rsbl.2015.0875

PubMed Abstract | CrossRef Full Text | Google Scholar

Tennessen, J. A., Woodhams, D. C., Chaurand, P., Reinert, L. K., Billheimer, D., Shyr, Y., et al. (2009). Variations in the expressed antimicrobial peptide repertoire of northern leopard frog (Rana pipiens) populations suggest intraspecies differences in resistance to pathogens. Dev. Comp. Immunol. 33, 1247–1257. doi: 10.1016/j.dci.2009.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Wellington, E. M., Berry, A., and Krsek, M. (2003). Resolving functional diversity in relation to microbial community structure in soil: exploiting genomics and stable isotope probing. Curr. Opin. Microbiol. 6, 295–301. doi: 10.1016/S1369-5274(03)00066-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitfield, S. M., Kerby, J., Gentry, L. R., and Donnelly, M. A. (2012). Temporal variation in infection prevalence by the amphibian chytrid fungus in three species of frogs at La Selva, Costa Rica. Biotropica 44, 779–784. doi: 10.1111/j.1744-7429.2012.00872.x

CrossRef Full Text | Google Scholar

Wickham, H. (2009). Ggplot2: Elegant Graphics for Data Analysis. New York, NY: Springer-Verglag. doi: 10.1007/978-0-387-98141-3

CrossRef Full Text | Google Scholar

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. Ecol. Arch. 96:595. doi: 10.1890/14-1837.1

CrossRef Full Text | Google Scholar

Woodhams, D. C., Bletz, M., Kueneman, J., and Mckenzie, V. (2016). Managing amphibian disease with skin microbiota. Trends Microbiol. 24, 161–164. doi: 10.1016/j.tim.2015.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Woodhams, D. C., Kilburn, V. L., Reinert, L. K., Voyles, J., Medina, D., Ibáñez, R., et al. (2008). Chytridiomycosis and amphibian population declines continue to spread eastward in Panama. Ecohealth 5, 268–274. doi: 10.1007/s10393-008-0190-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhalnina, K., Dias, R., De Quadros, P. D., Davis-Richardson, A., Camargo, F. A., Clark, I. M., et al. (2015). Soil pH determines microbial diversity and composition in the park grass experiment. Microb. Ecol. 69, 395–406. doi: 10.1007/s00248-014-0530-2

PubMed Abstract | CrossRef Full Text | Google Scholar