I sampled the cutaneous microbiota of spring peepers (Pseudacris crucifer), wood frogs (Lithobates sylvaticus), green frogs (Lithobates clamitans), bullfrogs (Lithobates catesbeianus), carpenter frogs (Lithobates virgatipes), and Pine Barrens treefrogs (Hyla andersonii) across five sites: Success Lake in Colliers Mills Wildlife Management Area (WMA), Jackson, NJ, Albertson Bog in Wharton State Forest, Hammonton, NJ, Morin Pond and Kai Pond in Somerset, NJ, and Imlaystown Bog in Assunpink WMA, Allentown, NJ (Supplementary Figure S1). Frog species were selected based on preliminary visual surveys, such that there was a sufficiently large population at each site to make sampling feasible. Sites were further selected based on permit regulations and pond accessibility. Two species – green frogs and bullfrogs – were sampled at more than one site, and only green frogs were sampled at more than two sites. Frogs were sampled in the spring and summer of 2016 and 2017, and I aimed to sample 10 individuals of target species at each site (Table 1). Sampling methods were as previously described (Kruger, 2020). Briefly, each frog was handled with new nitrile gloves and rinsed twice with sterile deionized water to exclude transient matter that is not part of the skin-associated microbial community (McKenzie et al., 2012; Kueneman et al., 2014). I swabbed each frog with a sterile cotton swab (Medline MDS202000) 20 times (five streaks each on the dorsal side, ventral side, and each hind limb). All amphibians were released immediately after sampling. Environmental bacteria were sampled by swirling a swab in pond water at approximately 10 cm depth for 5 s (n = 9 across sites, Table 1). Swabs were placed in sterile centrifuge tubes on ice and subsequently preserved at −80°C until DNA extraction. The New Jersey DEP approved this protocol for sampling amphibian skin bacteria (NJDEP Scientific Collecting Permit Nos. SC 2016093 and 2017053), and amphibian sampling methods were approved by Rutgers University’s Institutional Animal Care and Use Committee (Protocol #14-080).

I extracted DNA from swabs using the Qiagen DNeasy PowerSoil kit following kit instructions with modifications according to the Earth Microbiome Project for low DNA amounts (Berg-Lyons et al., 2018). After extraction, an initial round of PCR of the full 16S gene with 27F and 1492R primers was performed to deal with low DNA concentrations. Each 20 μL reaction contained 0.02 U/uL Phusion^® DNA polymerase, 200 μM dNTPs, 0.5 μM forward primer, 0.5 μM reverse primer, and 1 μL DNA template. The mixtures were amplified at 98°C for 30 s, followed by 35 cycles of 98°C for 10 s, 62°C for 30 s, 72°C for 45 s and a final elongation for 7 min. I checked products for successful amplification using gel electrophoresis. PCR products were cleaned using the AxyPrep Mag PCR cleanup kit per kit instructions and sent to the Integrated Microbiome Resource lab at Dalhousie University (Halifax, Canada) for library preparation and sequencing. A second round of PCR of the V4-V5 region using the 515F and 926R primers was performed, and pooled PCR amplicons were sequenced by an Illumina MiSeq using 2x300 bp paired-end v3 chemistry. Sequencing procedural details were as described by Comeau et al. (2017).

Sequence assembly was performed by the Integrated Microbiome Resource lab at Dalhousie University according to the Microbiome Helper standard operating procedures workflow (Comeau et al., 2017) using QIIME v. 1.9.1 (Caporaso et al., 2010). Briefly, sequences were clustered into operational taxonomic units (OTUs) based on 97% sequence similarity using open-reference OTU picking with QIIME wrapper scripts (Caporaso et al., 2010), and taxonomy was assigned using the Greengenes database (DeSantis et al., 2006). Sequences without taxonomic matches were clustered de novo at the 97% sequence similarity level. Sequences were quality filtered such that low-confidence OTUs making up <0.1% of reads and chimeric reads were removed. Sequencing depth per sample ranged from 530 to 77,737 reads. I rarefied all samples to a depth of 2500 reads to standardize sampling effort while maximizing sample inclusion, and samples that did not reach this threshold (n = 5) were removed from subsequent analyses. After filtering, 76 frog skin samples and 8 environmental samples remained (Table 1). Sequences have been deposited in the SRA database (Bioproject accession number: PRJNA601697).

I sent extracted DNA to Pennsylvania State University-Altoona College for Bd-testing using qPCR, where primers developed by Boyle et al. (2004) were used. All samples were run in triplicate with positive and negative controls using qPCR conditions described by Julian et al. (2016). A synthetic gBlock^® fragment of the ITS-1 and 5.8 ribosomal genes of Bd (GenBank accession number AY598034) was used to generate a standard curve ranging from 9 to 9000 zoospore equivalents to allow estimation of Bd concentrations. Individuals were considered Bd-positive if at least two out of three qPCR wells were positive (Kerby et al., 2013; Richards-Hrdlicka et al., 2013; Becker et al., 2015a) at a detection threshold of >1.0 zoospore equivalents (Supplementary Table S1). If one of the three wells was positive, the reaction was run in triplicate again. If at least one of the wells was positive on the second run of PCR, this was considered a positive detection of Bd. Bd infection intensity for each sample was calculated as log-transformed mean zoospore equivalents across technical (positive) replicates.

Shannon index, Chao1, Faith’s phylogenetic distance (Faith, 1992), and observed OTU richness alpha diversity metrics were calculated with QIIME (Caporaso et al., 2010). Beta diversity was calculated using weighted and unweighted Unifrac distances (Lozupone et al., 2011), Bray–Curtis dissimilarity (Bray and Curtis, 1957), and Jaccard indices in the phyloseq package (McMurdie and Holmes, 2013) in R version 3.5.1 (R Core Team, 2016). Unifrac metrics encompass information on phylogenetic relatedness (Lozupone et al., 2011), while Jaccard and Bray–Curtis do not. Jaccard and unweighted Unifrac use only a presence/absence matrix, while Bray–Curtis and weighted Unifrac use a matrix containing relative abundances of OTUs sampled on each individual. I visualized results using non-metric multi-dimensional scaling (NMDS).

I compared the bacterial OTUs identified in this study to OTUs in a published database containing bacterial isolates known to inhibit or enhance Bd activity (Woodhams et al., 2015). I used a custom blast search in Geneious version 11.1.3 (Kearse et al., 2012) to find which potentially inhibitory Bd OTUs were present in my dataset based on >97% sequence similarity. Sequence similarity does not necessarily imply Bd-inhibitory function, but rather suggests this function is likely (Muletz-Wolz et al., 2017a). I used the megablast tool and had Geneious return only the top hit from the antifungal isolate database (Muletz-Wolz et al., 2017b).

I conducted indicator species analyses using the “multipatt” function in the indicspecies package (De Cáceres and Legendre, 2009) in R (R Core Team, 2016). This analysis identifies species that characterize a group of sites, or in this case, bacterial OTUs that characterize different frog species, based on their relative abundance and frequency of occurrence (Dufrene and Legendre, 1997). OTU tables were filtered before analyses to exclude environmental bacterial samples, thereby only including amphibian-associated bacteria. I identified OTUs that had high association values (IndVal > 0.7, Becker et al., 2017) and were significantly (p < 0.05) associated with host species. I also performed an indicator analysis between Bd positive and negative individuals to determine which OTUs were strongly associated with individuals based on Bd status. I compared the indicator OTUs to the antifungal isolate database (Woodhams et al., 2015) to determine how many indicator OTUs were putative anti-Bd bacteria.

All statistical analyses were performed using R (R Core Team, 2016). I analyzed the effect of host species identity on each alpha diversity metric with a mixed effects model using the lme4 package (Bates et al., 2015). Host species was assessed as a fixed effect, site was included as a random effect, and significance of effects was determined using the “Anova” function in the car package (Fox et al., 2016). I used Shapiro–Wilk tests to assess normality and evaluated model fit using residual plots. Chao1 index data were square-root transformed and analyzed with a linear mixed effects model, and OTU richness data were analyzed using a generalized linear mixed effects model with a negative binomial error distribution. I used Kruskal–Wallis tests followed by FDR-corrected Wilcoxon tests to determine if bacterial community alpha diversity varied based on Bd status (positive vs. negative individuals) and sample type (frog vs. environment). Because green frogs were the only species sampled at more than two sites, I analyzed patterns across these populations to determine the effect of sampling site on alpha diversity using Kruskal–Wallis tests.

I determined the effect of host species on cutaneous bacterial beta diversity using permutational multivariate analysis of variance (PERMANOVA) with the “adonis” function (999 permutations) in vegan (Oksanen et al., 2013). When testing the importance of frog species in influencing community beta diversity, the “strata” argument accounted for sampling site. A separate PERMANOVA was used to compare beta diversity between Bd-positive and Bd-negative individuals. As described above, I analyzed patterns in the three green frog populations to determine the effect of sampling site on bacterial community structure. I also used a PERMANOVA to compare beta diversity between environment and frog-associated bacterial communities. Pairwise PERMANOVAs compared significant effects. I used the false discovery rate (FDR) procedure to correct for multiple comparisons (Benjamini and Hochberg, 1995). The “betadisper” function in vegan (Oksanen et al., 2013) was used to run tests for homogeneity of dispersions for frog species and sites (using green frog populations), and Tukey’s HSD test was used for post hoc analyses. Only significant betadisper results are reported, since heterogeneous dispersion among groups can influence PERMANOVA results.

To analyze the diversity (i.e., number and total relative abundance) of potentially anti-Bd bacteria, I used mixed effects models with Bd status and host species as fixed effects and site as a random effect. Among the green frogs, I used Kruskal–Wallis tests to determine if the number and total relative abundance of putative anti-Bd OTUs varied across sites. I compared the proportion of putative anti-Bd OTUs among all frog-associated OTUs to the proportion of these OTUs among indicator taxa using a Chi-square goodness of fit test to determine if more indicator OTUs were inhibitory than expected based on their prevalence among all frog-associated OTUs.

All alpha diversity metrics yielded similar results. Hence, here I present the diversity estimates of Chao1, a metric commonly used for measuring microbial diversity, and OTU richness. OTU richness was variable across individuals, ranging from 4 to 312 OTUs per frog. OTU richness (Figure 2; GLMM: χ² = 22.68, df = 5, p < 0.001) and Chao1 index (LMM: χ² = 11.54, df = 5, p = 0.042) varied significantly among host species. Among green frogs, the only frog species sampled at more than two sites, there were significant differences in alpha diversity (KW – OTU richness: χ² = 10.02, p = 0.007; Chao1: χ² = 10.29, p = 0.006) among populations. These results indicate that both host species and sampling location influenced individual-level bacterial community diversity.

Bray–Curtis, Jaccard, unweighted Unifrac, and weighted Unifrac yielded consistent results, therefore only Bray–Curtis results are presented for simplicity. Host species (Pseudo-F = 5.34, df = 5, R² = 0.26, p = 0.001) was significantly associated with bacterial beta diversity differences among individuals (Figure 3). Compositional differences among frog host species’ microbiota were significant between bullfrogs and green frogs, bullfrogs and spring peepers, and wood frogs and green frogs (pairwise PERMANOVA: p < 0.05). There were also site-level differences in beta diversity among the three green frog populations (Figure 3; PERMANOVA – Pseudo-F = 4.28, df = 2, R² = 0.25, p = 0.001).

I found significant differences in beta dispersion (distance to centroid) among host species (F_5,70 = 14.03, p = 0.001). These differences were driven by the relatively small distance to centroid observed among Pine Barrens treefrogs (Figure 4), indicating that Pine Barrens treefrogs’ skin bacterial communities were more compositionally similar to one another than individuals of other species. This pattern is also evident in the NMDS plots of Bray–Curtis distances (Figure 3), where Pine Barrens treefrog individuals cluster closely in ordination space and have a relatively narrow confidence ellipse.

No clinical signs of chytridiomycosis were observed during frog swabbing. Nine individuals tested positive for Bd (Table 1 and Supplementary Table S1), making Bd prevalence 11.8% across all individuals sampled (n = 76). One individual, a green frog at Morin Pond, yielded an inconclusive result with only one of three PCR wells testing positive for Bd. Because there was not enough DNA to re-test this individual, it was conservatively considered negative for Bd.

Bacterial alpha diversity did not differ between Bd-positive and Bd-negative individuals (KW – OTU richness: χ² = 1.89, p = 0.17; Chao1: χ² = 2.02, p = 0.15). Among spring peepers, the only species with roughly equal numbers of Bd positive (n = 4) and negative (n = 3) individuals, there were also no differences in alpha diversity based on Bd status (KW – OTU richness: χ² = 0.8, p = 0.37; Chao1: χ² = 1.13, p = 0.29), suggesting that this trend is not simply an artifact of unbalanced sample sizes. Using Pearson correlation, I did not find that either alpha diversity metric correlated with Bd infection intensity (Supplementary Table S1; OTU richness: R = 0.11, p = 0.34; Chao1: R = 0.11, p = 0.33).

Bd status explained very little variation in community structure and composition and was not significant in PERMANOVA results based on Bray–Curtis dissimilarities (Pseudo-F = 1.17, df = 1, R² = 0.02, p = 0.75). Furthermore, Bd-positive individuals were nested within overall community structure in ordination space (Figure 3). Taken together, these results indicate that differences among skin bacterial communities were not related to Bd status.

Of the 2194 frog-associated OTUs, 11.9% (260) matched inhibitory isolates in the antifungal isolate database. These putative anti-Bd OTUs made up 32.6% of bacterial abundance within host-associated bacterial communities. The mean number of anti-Bd OTUs across individuals was 21 (SD ±14), and the majority (80%) of anti-Bd OTUs were classified as Proteobacteria, although members of Bacteroidetes, Actinobacteria, and Firmicutes were also present.

There were significant differences in the number of anti-Bd OTUs across host species and between Bd-positive and Bd-negative individuals (Figure 5; GLMM – Species: χ² = 39.6, df = 5, p < 0.001; Bd status: χ² = 27.5, df = 1, p < 0.001). However, these differences were no longer significant when considering the relative abundance of anti-Bd OTUs among host species and between Bd-positive and Bd-negative individuals (Figure 5; GLMM – Species: χ² = 6.40, df = 5, p = 0.27; Bd status: χ² = 0.37, df = 1, p = 0.54), suggesting that OTU richness is not always related to abundance. There were significant differences in both the number (KW χ² = 12.81, p = 0.002) and relative abundance (KW χ² = 10.67, p = 0.005) of anti-Bd OTUs across the three green frog populations (Figure 6).

PERMANOVA results indicated that host species identity (Pseudo-F = 4.25, df = 5, R² = 0.21, p = 0.001) significantly influenced anti-Bd community beta diversity (Supplementary Figure S2). Additionally, there were significant site-level differences in putative anti-Bd OTUs among the three populations of green frogs (Figure 6; PERMANOVA – Pseudo-F = 4.79, df = 2, R² = 0.27, p = 0.001). However, Bd status had no significant effect on anti-Bd community beta diversity (Pseudo-F = 1.12, df = 1, R² = 0.015, p = 0.8). Similar to the overall beta diversity results, there were significant differences in beta dispersion for anti-Bd OTU communities among host species (F_5,70 = 10.4, p = 0.001). Once again, these differences were driven by the relatively small dispersion among Pine Barrens treefrogs. These trends suggest that patterns in anti-Bd community diversity reflect patterns in overall community diversity.

There were 71 OTUs strongly associated (IndVal > 0.7) with at least one frog species (Figure 7). Of these indicator taxa, 30 (42.3%) were putative anti-Bd OTUs, which is significantly more than expected based on the proportion of putative anti-Bd OTUs among the 2194 OTUs in the larger dataset (χ² = 62.94, df = 2, p < 0.001). The majority (84.5%) of indicator OTUs belonged to the phylum Proteobacteria. Only 13 of all indicators were associated with more than one frog species (Supplementary Table S2).

In the indicator analysis based on host Bd status, six OTUs were highly associated with Bd-positive individuals. All six OTUs belonged to the family Comamonadaceae (Phylum: Proteobacteria, Class: Betaproteobacteria, Order: Burkholderiales). Three out of six of these OTUs were putative anti-Bd isolates based on consensus matches in the antifungal isolate database, indicating that several potential anti-Bd OTUs were strongly associated with Bd-positive individuals.

Within the environmental samples, there were no significant differences in OTU richness (KW χ² = 1.83, p = 0.77) or Chao1 index (KW χ² = 2.75, p = 0.6) among sites. Mean OTU richness was significantly higher in environmental samples compared to frog samples (KW χ² = 4.68, p = 0.03), but there were no differences in Chao1 index (KW χ² = 3.8, p = 0.051). Beta diversity was significantly different between environmental and frog samples (PERMANOVA – Pseudo-F = 1.63, df = 1, R² = 0.02, p = 0.02) when samples were pooled across species and sites. Taken together, these results suggest differences in community structure between frog and environmental bacterial communities.