Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas)

Abstract

Increasingly, host-associated microbiota are recognized to mediate pathogen establishment, providing new ecological perspectives on health and disease. Amphibian skin-associated microbiota interact with the fungal pathogen, Batrachochytrium dendrobatidis (Bd), but little is known about microbial turnover during host development and associations with host immune function. We surveyed skin microbiota of Colorado's endangered boreal toads (Anaxyrus boreas), sampling 181 toads across four life stages (tadpoles, metamorphs, subadults and adults). Our goals were to (1) understand variation in microbial community structure among individuals and sites, (2) characterize shifts in communities during development and (3) examine the prevalence and abundance of known Bd-inhibitory bacteria. We used high-throughput 16S and 18S rRNA gene sequencing (Illumina MiSeq) to characterize bacteria and microeukaryotes, respectively. Life stage had the largest effect on the toad skin microbial community, and site and Bd presence also contributed. Proteobacteria dominated tadpole microbial communities, but were later replaced by Actinobacteria. Microeukaryotes on tadpoles were dominated by the classes Alveolata and Stramenopiles, while fungal groups replaced these groups after metamorphosis. Using a novel database of Bd-inhibitory bacteria, we found fewer Bd-inhibitory bacteria in post-metamorphic stages correlated with increased skin fungi, suggesting that bacteria have a strong role in early developmental stages and reduce skin-associated fungi.

Figure 1

Major skin microbial taxa on boreal toads from Trout Creek at each life stage, collected on the same day. (a) The proportional abundance of bacterial OTU sequences per major taxon across lifestages: tadpoles (n=8), metamorphs (n=8), subadults (n=8) adults (n=10). (b) The proportional abundance of microeukaryote OTU sequences per major taxon across lifestages: tadpoles (n=5), metamorphs (n=4), subadults (n=5), adults (n=5). Data in a is based on 10 750 16S sequences per sample and data in b is based on 2400 18S sequences per sample. Taxonomy is designated by p=phylum and c=class.

Figure 2

Shannon diversity of microbial communities on boreal toads, by lifestage. (a) Alpha diversity (Shannon) of bacteria by life-stage from Trout Creek (n=32). (b) Alpha diversity (Shannon) of microeukaryotes by life stage from Trout Creek (n=19).

Figure 3

Analysis of co-occurring bacterial and fungal OTUs found on boreal toads at one site. This analysis included 19 toad individuals from Trout Creek for which we had successful sequencing yields for both 16S and 18S: n=5 tadpoles, n=4 metamorphs, n=5 subadults and n=5 adults. (a) Network analysis depicting only significantly co-occurring bacterial and fungal OTUs. All square nodes represent OTUs (either bacteria or fungi) with significant co-occurrence with another OTU in the network. Red lines indicate negative co-occurrences between two OTUs. Turquoise lines indicate positive co-occurrences between two OTUs. Yellow=Betaproteobacteria, purple=Actinobacteria, blue=Fungi, white=other bacterial OTUs. Bacteria are represented by taxa that were observed to be more abundant on tadpoles (Betaproteobacteria) or on adult stages (Actinobacteria, see Figure 1a). (b) Heatmap depicting the number of negatively co-occurring OTUs between bacterial and fungal taxa. (All positively correlated interactions are excluded from the heatmap, but are available in Supplementary Data via DRYAD). The numbers indicate the frequency of negative co-occurrences between bacterial and fungal groups calculated as the sum of negative OTU co-occurrences per bacterial order and fungal class, which are also depicted as red lines in a.

Figure 4

Beta diversity of skin bacterial communities of A. boreas from all sites, differentiated by life-stage. Dark blue=adults (n=77), red=subadults (n=28), green=metamorphs (n=47), black=tadpoles (n=16). Diversity patterns are visualized using principle coordinate plots of weighted UniFrac distances. Each point represents the skin bacterial community of an individual boreal toad; symbol color and shape indicate life-stage.

Figure 5

Bd-inhibitory OTUs found on boreal toads and in environmental samples. (a) Mean percent of bacterial sequencing reads with Bd-inhibitory function, from boreal toad life-stages or environmental samples (across all sites, sample sizes in parentheses). Bd-inhibitory function was determined by matching OTUs to the database described in the methods and Woodhams et al., 2015. Four different bioinformatic algorithms were used to identify Bd-inhibitory OTUs as described in the methods section: exact match sequences, trimmed exact matches, exclusively inhibitory matches, and expected inhibitory matches. All four methods showed similar patterns and significant differences among sample type (Kruskal–Wallis tests, P's<0.01) with tadpoles, and to a lesser extent metamorphs, having a significantly higher percentage of their skin bacteria comprised of Bd-inhibitory OTUs. (b) A heatmap showing the 14 bacterial taxa that contain the ‘exclusively Bd-inhibitory matches' (479 OTUs, shown in aqua in Figure 5a) and how they are distributed across boreal toad lifestages and environmental samples. The lifestage and environmental sample x axis labels apply to both a and b.