Skin bacterial community differences among three species of co-occurring Ranid frogs

Field site

We captured, swabbed, and released amphibians at three different sites on or near the Mianus River Gorge in Bedford, New York, USA, to detect Bd infection and to collect samples of the skin bacterial communities (Table 1; Fig. 1). The sampling took place in June 2017 and 2018. Site 1 was located behind a residential communities, accessible by a dirt road. Site 1 consisted of a stream that varied in forest cover and the density of understory cover and split off into an open wetland area. Additionally, across an open meadow there was a vernal pool with partial forest cover. Site 1 was the only site where skin bacterial communities were examined, and therefore, was divided into five subsites based on differences in habitat (Fig. 1). Amphibians were sampled at site 1 in 2017 and 2018, but skin bacterial samples were only taken from 2017. Site 2 consisted of a roadside pond with no forest cover, on private land, that led to a wetland that in turn fed into another pond, both on MRGP property. The wetland and ponds on MRGP property both had partial forest cover and had a dense understory dominated by Symplocarpus foetidus (skunk cabbage). All three of these water bodies were sampled as a single site (site 2) in 2017. Site 3 was along a stream that started in a forested area near two vernal pools and led into a thick marshland with less forest cover. Amphibians were collected near the marshland and at the vernal pools in 2018.

Amphibian sampling

We collected skin swabs from L. sylvaticus, L. clamitans, L. catesbeianus, and L. palustris. At all three sites, we collected samples from L. sylvaticus and L. clamitans, while L. catesbeianus was only sampled at Site 2, and L. palustris was only found and sampled at sites 1 and 3 (Table 1). We also sampled Lithobates tadpoles, which were not identified to species, at site 1 in 2017 and 2018, at site 2 in 2017, and at site 3 in 2018 (Table 1). The use of animals was approved by Virginia Tech’s Institutional Animal Care and Use Committee (16-193-STAT) and the state of New York (permit #2213).

We sampled only one site and one species per day. Using dip nets, we caught amphibians, which were then placed individually into sterile Whirl-Pak bags until our sample size at the site was reached to avoid re-sampling individuals (adults per site: 20 in 2017 and 25 in 2018; tadpoles per site: 30 in 2017 and 20 in 2018). For tadpoles, we included some water from the sample site in the bags. Once the sample size for that day was reached or it was noon, the amphibians were weighed in the Whirl-Pak bags. Wearing sterile nitrile gloves, we then removed each amphibian from its bag and weighed the Whirl-Pak bag by itself, with amphibian mass as the difference. The amphibian was then rinsed with 50 ml of autoclaved distilled water to remove dirt and environmental bacteria. We then swabbed amphibians using a single sterile rayon swab (MW113, Medical Wire Equipment, Corsham, England). Swabbing of adults consisted of five strokes in one direction on each of the hind feet and thighs, and 20 strokes on the ventral side (total strokes = 40), while for tadpoles, we swabbed around their mouth 25 times, an area commonly infected by Bd in this lifestage (Marantelli et al., 2004). We then placed the swab in a sterile 1.5 ml microcentrifuge tube and placed it on ice while in the field. Lastly, we measured the snout-vent length of the adult frogs, and the tail length and total length of tadpoles, before they were released back in the site. Upon returning from the field (<6 h after sampling), we stored the swabs in a −20 °C freezer. At the end of June, all samples were transferred to Virginia Tech, Blacksburg, VA, and stored in a −80 °C freezer.

Environmental data collection

We took temperature and pH measurements (Oakton Waterproof pHTestr 30) at locations where amphibians were sampled each day. Water temperatures (overall mean = 19.7, sd = 3.51) were similar across the different sampling locations (Fig. S1A; Fig. S2A). The pH also was similar across all sampling sites (all with mean ~7), except the vernal pool subsite at site 1, which had a lower mean pH of 4.75 (Fig. S1B; Fig. S2B).

Sample and data processing

All amphibian samples, across all three sites, were used to determine the prevalence of Bd in MRGP. However, bacterial community data was only obtained from 2017 samples across all five subsites at site 1, and consisted of 66 swabs from three amphibian species; L. sylvaticus, L. clamitans, and L. palustris (Table 1).

Swab DNA extraction

The swabs were transported to Virginia Tech in Blacksburg, VA, where they were all processed by a single individual (ZG). DNA was extracted from the swabs using a Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) with the lysozyme pre-treatment for gram-positive bacteria. The lysozyme pre-treatment for each sample consisted of adding 0.18 ml of lysis buffer and 3.7 g of lysozyme, vortexing, and incubating for one hour at 37 °C. We then stored the extracted DNA in 100 µl of molecular grade water in a −20 °C freezer, and used it as template DNA for Bd PCR and assessing skin bacterial communities using 16S rRNA gene amplicon sequencing.

Bd detection

Following DNA extraction, we screened all the amphibian swabs collected (2017 and 2018, N = 408) for Bd using PCR (Annis et al., 2004). Each 25 µl PCR reaction, one per sample, contained: 0.5 µl of dNTPs, 0.2 µl of Taq DNA Polymerase, 3.9 µl of Taq Buffer with MgCl₂, 2.5 µl of ITS 1–3 primer, 2.5 µl of 5.8S primer, 2 µl of extracted sample DNA, and 13.4 µl of water. The thermocycler conditions were: 93 °C for 10 min to start, then 30 cycles of 93 °C for 45 s, 60 °C for 45 s, and 72 °C for 1 min, and a final 10 min at 72°C. Every thermocycler run had a positive (extracted DNA from a Bd JEL 404 stock) and negative (molecular grade water) control. Each PCR product was run on a 1% TAE agarose gel. The gels were inspected to ensure the amplification in the positive sample and no amplification in the negative sample. We recorded a sample as Bd positive if a band was seen on the gel.

Skin bacterial communities

To assess the skin microbiome, we completed 16S rRNA gene amplicon sequencing of the 66 samples collected in 2017 at site 1 (17 L. sylvaticus, 25 L. palustris, 12 L. clamitans, and 12 Lithobates tadpoles). We amplified the V4 region of the 16S rRNA gene using the 515F and 806R primers. The 806R primer contained a 12bp error-correcting Golay code to mark individual samples. Each 25 µl reaction contained 0.5 µl of both the 515F and 806R primers (at 10 µM concentration), along with 12 µl of ultra-clean PCR grade water, 10 µl of 5Prime hot master mix, and 2 µl of the template DNA. Each 25 µL PCR reaction was run in triplicate, along with a negative control reaction that did not include any template DNA. The thermocycler conditions were: 94 °C for 3 min to start, then 35 cycles of 94 °C for 45 s, followed by 50 °C for 1 min, and 72 °C for 1.5 min, and a final 10 min at 72 °C. At the end of the PCR run, we combined each sample’s triplicate PCR products into one tube. We visualized the combined PCR product on a 1% TAE agarose gel, where the negative samples were checked for contamination and the sample wells were checked for amplification. The PCR products were quantified using a Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA) with a dsDNA High Sensitivity assay kit. We combined 200 ng of DNA from each sample (N = 66) into a final pool. We then purified the pooled sample using the QIAquick PCR purification kit (Qiagen, Hilden, Germany) and sent it to the Genomics Center at the Dana Farber Cancer Institute of Harvard University for 250 bp single-end sequencing on an Illumina Mi-Seq instrument.

Sequence processing

We processed the 250 bp forward (single-end) reads using the QIIME2 pipeline (Bolyen et al., 2018). We imported the raw reads and demultiplexed them. The reads were of consistently high quality across the full-length, so we did not need to trim them. We denoised the data using DADA2 (Callahan et al., 2016), which included filtering out phiX and chimeric reads, as well as correcting amplicon errors. We used the recommended ‘big data’ parameters for DADA2, which included truncating reads with a quality score cut-off of 11. In addition, we only used 10,000 reads to build the error distribution, which we have found is adequate for our high-quality data and significantly decreases the run time. The resulting amplicon sequence variants (ASV) table contained 45,080 ASVs. We then filtered out any ASVs that were present at less than 0.01% of the total read count, to work with ASVs that are likely more impactful in the system. This left 1,103 ASVs in the table. Taxonomy was assigned to these ASVs using the SILVA v132 database (Quast et al., 2012) and scikit-learn classifier (Pedregosa et al., 2011). We then filtered out any ASVs that were assigned as chloroplasts, mitochondria, or were unassigned; this cut the number of ASVs to 1,079. After visualizing the alpha rarefaction curve, the dataset was rarefied at 20,000 reads, which resulted in the loss of four samples with lower read counts (three L. sylvaticus and one L. palustris). The final table contained 62 samples (Table 1) and 1,079 ASVs.

Differences in bacterial communities between amphibian species and life stage and sampling subsites

To assess differences in the community structure of ASVs across amphibian species and life stage, the skin bacterial community data was visualized with non-metric multidimensional scaling (NMDS) plots using a Bray-Curtis dissimilarity index (Fig. 2A) and Jaccard dissimilarity index (Fig. 2B). The Bray-Curtis dissimilarity index uses relative abundance data and the Jaccard dissimilarity index uses presence and absence data. Skin bacterial communities were also visualized using a Bray-Curtis dissimilarity index and compared across subsites, were the amphibians were collected from (Fig. 3). We used two permutational multivariate analysis of variance (PERMANOVA; vegan: adonis; Oksanen et al., 2013) to test if the multivariate means (centroids) and variances (dispersion) of the skin bacterial samples, grouped by (1) amphibian species and life stage and (2) subsites, were different. The subsite PERMANOVA included all subsites besides the vernal pool, where the majority of samples were tadpoles, and therefore, development was a confounding variable. Using the pairwise.adonis function, we then performed a post-hoc test on the two PERMANOVAs to determine which species or life stage and subsites were different from each other (Martinez Arbizu, 2020). To check the homogenous dispersion among (1) host species and life stage and (2) subsites, we used the betadisper function which calculates the samples distance from the centroid and uses a permutation approach to tests for statistical differences (vegan: permutest; Oksanen et al., 2013). Again, the vernal pool subsite was not included in the dispersion analysis due to the majority of samples being from tadpoles. Each permutation test was run twice, once for each dissimilarity matrix, with 999 permutations for each test. Lastly, we ran two linear discriminate effects size analyses (LEfSe), using the run_lefse function in the microbiomeMarker package (Cao et al., 2022), to compared bacteria phylum abundances across (1) adults species and tadpoles and (2) sampling subsites.

To determine how bacterial diversity measurements (Shannon diversity and ASV richness) differed across amphibian species and between life stages and sampling sites the five subsites in site (1), we used a series of generalized linear models (GLMs). For each alpha diversity measurement, we fit two models, one for subsites at site 1 and one for amphibian species and life stage, for a total of 4 GLMs. Shannon diversity GLM models were fit assuming a Gamma distribution with an inverse link function as these metrics are continuous and strictly positive. The ASV richness metrics were fit with a negative binomial GLM with a log link function since this metric is discrete data and overdispersed. After fitting the models, we used the emmeans function from the emmeans package (Lenth et al., 2018) to make contrasts, and used the Tukey method to adjust for multiple testing between adult amphibian species and lifestage, and sampling sites for each of the metrics independently.

Differences in species and site infection prevalence

We ran a logistic regression to compare adult species and tadpole’s infection prevalence. Using the emmeans function from the emmeans package (Lenth et al., 2018), we made contrasts to compare species and tadpoles. Additionally, we ran a mixed effect logistic regression, with species as a random effect and site as a fixed effect, using the glmTMB function from the glmTMB (Brooks et al., 2017). Lastly, we again ran contrasts with the eemeans package (Lenth et al., 2018) to compare sites to each other.

Differences between bacterial communities in infected and non-infected Lithobates palustris

To compare bacterial communities between infected and non-infected amphibians, we only used L. palustris since this was the only species with enough samples of infected individuals (2017: N = 12 infected, 13 uninfected) for a meaningful statistical comparison. Differences in bacterial communities between infected and non-infected L. palustris were compared using logistic GLMs with alpha diversity measures set as predictor variables. Each alpha diversity measure addresses a different aspect of diversity (but are correlated) and therefore only one was used in a model at a time. We also examined the bacterial community structure between infected and non-infected individuals using a PERMANOVA and permutation test, as described above.

Differences in bacterial communities between amphibian species and life-stages

Adult amphibians had similar ASV richness (all emmeans contrasts p-values > 0.05). We found L. clamitans had a significantly higher Shannon diversity (Fig. 5B) than L. sylvaticus (coefficient = 0.0894, SE = 0.0310, Z.ratio = −2.882, p-value = 0.018), but both species had similar Shannon diversity compared to L. palustris (all contrasts p-values > 0.05). Lastly, when we compared bacterial communities of adults to tadpoles, we found that tadpoles had a higher ASV richness than the adults from two species (tadpoles-L. palustris; coefficient = 0.3671, SE = 0.128, Z.ratio = 2.873, p-value = 0.0212; tadpoles-L. sylvaticus; coefficient = 0.3889, SE = 0.142, Z.ratio = 2.733, p-value = 0.0319; other p-value > 0.05), and a higher Shannon diversity than the adults of one of these species (tadpoles-L. sylvaticus; coefficient = 0.1187, SE = 0.0294, Z.ratio = 4.039, p-value < 0.001; all other p-value > 0.05).

Bacterial community structure was different between amphibian species and life stages (Fig. 2; PERMANOVA Bray-Curtis: F_(3,58) = 4.8837, R² = 0.202, p-value = 0.001; Jaccard: F_(3,58) = 4.488, R² = 0.188, p-value = 0.001). Specifically, we found that L. sylvaticus bacterial community structure was different than both L. palustris (p-value = 0.006) and L. clamitans (p-value = 0.006). But, L. palustris and L. clamitans had similar bacterial community structures (p-value > 0.05). Additionally, we found that tadpole bacterial community structures were different than all adult amphibians species (all p-values < 0.01). Further, we found that dispersion was significantly different among amphibians species and life stage (permutest Bray-Curtis: F_(3,58) = 7.305, p-value = 0.001; Jaccard: F_(3,58) = 7.305, p-value = 0.001). We found that L. palustris had a greater dispersion, a greater variation in bacterial communities among samples, than L. sylvaticus, and L. clamitans, and L. palustris also showed a greater dispersions than tadpoles. This potentially indicates that L. palustris bacterial community structural differences might be due to differences in dispersion rather than centroids, multivariate means.

Differences in species and site infection prevalence

Infection prevalence was greatest in L. palustris at site 1, with 48% infected in 2017 and 59% infected in 2018 (Table 1). From the logistic regression, we found that infection prevalence was greater in L. palustris compared to L. sylvaticus (coefficient = −2.090, p-value < 0.001), L. clamitans (coefficient = −3.411, p-value < 0.001), and tadpoles (coefficient = −1.629, p-value < 0.001). However, L. palustris infection prevalence was not significantly greater than L. catesbeianus (coefficient = −2.170, p-value = 0.2454), which was only sampled at one site and year. Lastly, when we used species as a random effect in a logistic regression with site as a fixed effect we found that site 2 had a higher infection prevalence than both site 1 (coefficient = 1.64, p-value = 0.0014) and site 3 (coefficient = 3.23, p-value < 0.001). While, site 1 had a higher infection prevalence than site 3 (coefficient = 1.60, p-value = 0.0054).

Differences between bacterial communities in infected and non-infected Lithobates palustris

A comparison between the infected and non-infected L. palustris bacterial communities suggested no differences in any of the alpha diversity measurements (Shannon diversity: glm, b = 0.2032, SE = 0.2854, Z = 0.712, p-value = 0.467; ASV richness: glm, b = −0.0005, SE = 0.0041, Z = −0.122, p-value = 0.903). We also did not find any differences in bacterial community structure (Figs. 2A and 2B) between infected and non-infected L. palustris (PERMANOVA Bray-Curtis: F_(1,22) = 0.717, R² = 0.032, p-value = 0.871; PERMANOVA Jaccard: F_(1,22)= 0.832, R² = 0.0362, p-value = 0.873).

Differences between amphibian bacterial communities across site 1

We found that adult amphibian skin bacterial community structured was different across sampling sties (Fig. 3; PERMANOVA Bray-Curtis: F_(3,43) = 2.2579, R² = 0.1340, p-value = 0.001; Jaccard: F_(3,43) = 1.305, R² = 0.092, p-value = 0.002). Adult amphibians found at the bog site had a difference skin bacteria structure than adult amphibians from other subsites (all p-values < 0.05). Additionally, adults sampled from the tunnel subsite had different skin bacterial community structures than adults sampled form the Ann’s meadow subsite (p-value = 0.012). Lastly, we found no significant difference in dispersions between subsites were adult skin bacterial communities were sampled (Bray-Curtis: F_(3,43) = 2.3994, p-value = 0.089; Jaccard: F_(3,43) = 1.018, p-value = 0.387). The vernal pool subsite was not included in the subsite beta diversity analysis due to the majority of the samples being from tadpoles, and thus, development and site are confounded. However, in the subsite NMDS, Lithobates tadpoles and L. clamitans, collected in or around the vernal pool, clustered together, and separately from the other samples (Fig. 3). Lastly, we found that amphibians in site 1 at different subsites had similar skin bacterial ASV richness and Shannon diversity (all GLM p-values > 0.05; Figs. 5C and 5D).