Wild Bd infection dynamics
Bd load was lower in all enzootic populations relative to the epizootic (ANOVA F(3, 52) = 33.65, p < 0.001, SI Table 2, Fig. 1a). Bd prevalence was 100% in all populations except Ibon Acherito (79%).
Host skin bacterial and fungal taxonomic composition is associated with disease dynamic and population
Bacterial Shannon diversity differed based on population (ANOVA F(3, 52) = 8.58, p < 0.001), with Acherito being significantly higher than Arlet (p < 0.001) and Lhurs (p < 0.001). Fungal Shannon diversity did not differ by population (p > 0.05). Disease dynamic and population varied in beta diversity of bacteria and fungi (Bacteria: disease dynamic PERMANOVA Pseudo-F(1, 54) = 7.08, R
2 = 0.12, p = 0.001, population PERMANOVA Pseudo-F(3, 52) = 6.04, R
2 = 0.26, p = 0.001; Fungi: disease dynamic PERMANOVA Pseudo-F(1, 33) = 2.89, R
2 = 0.08, p = 0.002, population PERAMANOVA Pseudo-F(3, 31) = 2.79, R
2 = 0.21, p = 0.001, Fig. 1b, c).
Twenty-nine bacterial ASVs were associated with epizootic dynamics and a single ASV was associated with enzootic populations from sPLS-DA (Fig. 1d, Supplementary Data 1). Nine fungal OTUs were associated with epizootic dynamics on sPLS-DA component 1 (Fig. 1e), whilst component two had 22 and 8 OTUs linked to enzootic and epizootic dynamics respectively (Supplementary Data 2). Seven out of nine fungal OTUs on component 1 were classified as Batrachochytrium (Fig. 1e), whilst OTUs on component 2 belonged to eight classified fungal classes (Supplementary Data 2). Differential abundance analysis using ALDEx2 identified 21 bacterial biomarkers for disease state, of which 1 and 20 were associated with enzootic and epizootic disease dynamics respectively (Supplementary Data 3). Discriminatory bacterial taxa spanned four phyla and seven classes. For the enzootic populations, the single discriminatory ASV (also identified using sPLS-DA) was assigned to the Massilia genus (order: Burkholderiales), whilst ASVs in the epizootic population belonged to the Sphingobacteriales, Actinomycetales, Burkholderiales, Chlamydiales, Flavobacteriales, Rhizobiales and Xanthomonadales. ALDEx2 identified one fungal OTU (Batrachochytrium sp.) significantly associated with epizootic dynamics (Supplementary Data 4).
Host bacterial community function associates with disease dynamic and population
Procrustes analysis showed significant correlation between shotgun bacterial taxonomic data and 16S data (Procrustes m12 squared = 0.43, correlation in symmetric Procrustes rotation = 0.76, p = 0.001). Shannon diversity of KO hits differed based on population (ANOVA F(3, 43) = 9.56, p < 0.001) with Acherito being significantly higher than Arlet (p < 0.001) and Lhurs (p < 0.001), whilst Puits was also higher than Arlet (p = 0.03). KO beta diversity differed based on disease dynamic and population (disease dynamic PERMANOVA Pseudo-F(1, 45) = 6.80, R
2 = 0.13, p = 0.001, population PERMANOVA Pseudo-F(3, 43) = 4.50, R
2 = 0.24, p = 0.001, Fig. 2a–d). Analysis of bacterial KOs using sPLS-DA found the main sources of variation was in amino acid metabolism, carbohydrate metabolism and energy metabolism (Fig. 2e, Supplementary Data 5). ALDEx2 identified 72 KOs, of which 17 and 55 were associated with epizootic and enzootic populations respectively (Supplementary Data 6).
Host skin metabolome maps to disease dynamic and population
2014 metabolite features were measured in negative ion mode. Disease dynamic and population were significant factors driving variability in metabolome (PERMANOVA dynamic Pseudo-F(1, 54) = 10.93, R
2 = 0.17 p = 0.001, population Pseudo-F(3, 52) = 13.14, R
2 = 0.43, p = 0.001, Fig. 3a). Univariate analysis identified 870 differentially abundant metabolite features, of which 183 had a log2 fold change > 1.5 (Fig. 3b, Supplementary Data 7). PLS-DA based on disease dynamic yielded 121 metabolite features with a VIP score > 2, of which 111 were considered most informative (VIP > 2, q < 0.05, log2 fold change > 1.5, Supplementary Data 8). Following putative annotation of the most informative metabolite features, one of the features was suspected to be the anti-Bd metabolite indole-3-carboxaldehyde. More detailed investigation of the MS/MS data for this metabolite (m/z 144.04606 RT 6.96 min) kindled our interest in a second, following the peak of the same m/z in the spectra at 8.73 min, present in the XCMS matrix but originally filtered out as thought to be due to some higher blank signal. However, upon closer inspection, the intensities of this second peak behaved similarly to m/z 144.04606 with higher intensity in the epizootic population (Wilcoxon’s test W = 17, p < 0.001, SI Fig. 2). Whilst both ions were putatively annotated as the anti-Bd metabolite indole-3-carboxaldehyde, the MS/MS data for the second signal resembled those in spectral databases (https://mona.fiehnlab.ucdavis.edu/spectra/display/PR100508; https://mona.fiehnlab.ucdavis.edu/spectra/display/PB000507) and from an authentic standard fragmented in-house, thus providing evidence for the presence (and increased abundance in the epizootic population) of indole-3-carboxaldehyde on wild A. obstetricans skin.
Microbiome/metabolome interactions in the wild
Procrustes analysis showed a significant association between bacterial/fungal taxonomy and metabolic profile (bacteria: m12 square = 0.45, correlation in symmetric Procrustes rotation = 0.74, p = 0.001; fungi: m12 square = 0.60, correlation in symmetric Procrustes rotation = 0.63, p = 0.004). DIABLO confirmed enrichment of eight bacterial ASVs in the epizootic population that were also identified from a single omics analysis. Of the five informative metabolites recovered from our analysis, three were indicators of epizootic disease dynamics and one was linked to enzootic dynamics (Wilcoxon’s test q < 0.05, log2 fold change > 1.5, Supplementary Data 9). Network analysis identified a single cluster representative of epizootic dynamics based on the presence of epizootic discriminatory taxa and their positive correlations with the epizootic enriched metabolites (Fig. 3c). Meanwhile, negative correlations were found between epizootic taxa and the single metabolite feature that was enriched in the enzootic populations.
Unsupervised sPLS regression supports findings from supervised analysis with 31 of the 40 bacteria-metabolite pairs from DIABLO also recovered in sPLS. We identified 450 correlated pairs of bacterial ASVs/metabolite features (39 negative, 411 positive) involving 20 ASVs and 157 metabolite features forming two sub-networks (SI Fig 3a,b, Supplementary Data 10). Bacterial sub-network A was reflective of epizootic dynamics based on a dominance of epizootic discriminative ASVs that were positively correlated with epizootic enriched metabolite features. Sphingobacterium had the greatest betweenness centrality (a measure of how often a node occurs on all shortest paths between two other nodes and therefore an indicator of functional importance). In bacterial sub-network B, an ASV belonging to the Sinobacteraceae family and Haliscomenobacter exhibited the greatest betweenness centrality.
Fungal-metabolome analysis also yielded two sub-networks comprising 461 fungal-metabolite pairs (152 negative, 309 positive) involving 11 fungal OTUs and 109 metabolite features (SI Fig 4a,b, Supplementary Data 11). Fungal sub-network A comprised 8 OTUs, all of which were discriminant for epizootic dynamics from sPLS-DA, with ALDEx2 also identifying one OTU (sp633) as discriminant for epizootic dynamics. Fungal sub-network A was dominated by chytrid with 7 OTUs classified as Batrachochytrium sp. Twenty-seven metabolite features in fungal sub-network A were positively correlated with fungal taxa that were discriminative for epizootic dynamics, whilst twenty-three metabolite features were negatively correlated with fungal taxa in subnetwork A that were discriminative for enzootic dynamics. Fungal sub-network B comprised three OTUs—two belonging to the Massarinaceae and an unknown Dothidiomycete.
Experimental infection alters skin microbiome state
Nineteen out of twenty Bd-exposed animals became infected with Bd by day 30, with Bd infection intensity ranging from 0.88 to 427.08 genomic equivalents (GE) (SI Fig 5). On day 60, a single Bd-exposed animal remained infected, with all other animals clearing infection. There was no mortality.
We show no effect of Bd exposure on bacterial Shannon diversity at any time point measured (p > 0.05). Fungal Shannon diversity was reduced in Bd-exposed animals compared to the control group on day 30 of the experiment (W = 109, p = 0.005) but not on day 0 or 60 (p>0.05). We confirm no significant difference in beta diversity for bacteria and fungi at the start of the experiment (PERMANOVA p > 0.05). Beta diversity differed by day, Bd exposure and their interaction for bacteria (PERMANOVA day F(1,116) = 34.24, R
2 = 0.22, p = 0.001, Bd F(1,116) = 4.92, R
2 = 0.03, p = 0.001, day*Bd F(1,116) = 3.15, R
2 = 0.02, p = 0.005) and fungi (PERMANOVA day F(1,80) = 7.52, R
2 = 0.08, p = 0.001, Bd F(1,80) = 2.05, R
2 = 0.02, p = 0.034, Day*Bd F(1,80) = 2.30, R
2 = 0.03, p = 0.015, Fig. 4a, b). Twenty bacterial ASVs and two fungal OTUs were found to drive differences in beta diversity on day 30 according to sPLS-DA (Supplementary Data 12, 13). The bacterial ASVs spanned three phyla and eight classes, whilst fungal OTUs all belonged to the Basidiomycete phylum and the Class Tremellomycetes. ALDEx2 found no differentially abundant bacterial or fungal taxa at day 0. At day 30, three bacteria and one fungal taxa differed between treatments (Supplementary Data 14, 15). Sphingobacterium and Acinetobacter were associated with the control group, whilst Comomonas was increased in the Bd exposure group. The fungus Trichosporon had significantly higher abundance in the control group compared to Bd-exposed animals on day 30. For day 60, we identified three discriminatory taxa (Comamonas, Brevundimonas and Gemmobacter, Supplementary Data 16), all of which were associated with Bd exposure. No fungal taxa were differentially abundant at day 60.
KO Shannon diversity at day 30 was lower in the Bd-exposed group than the control (t-test df = 17.43, t = 2.96, p = 0.009). KO beta diversity differed by Bd exposure (PERMANOVA F(1, 18) = 5.69, R
2 = 0.24, p = 0.001, Fig. 5a–c). Forty discriminatory KOs were identified from sPLS-DA (Fig. 5d, Supplementary Data 17). ALDEx2 identified 303 discriminatory KOs, with 114 and 189 increased in the control and Bd-exposed groups respectively (Supplementary Data 18). Differentially abundant KO hits were mostly associated with amino acid metabolism, carbohydrate metabolism, membrane transport, signal transduction and metabolism of cofactors and vitamins.
UHPLC-MS metabolomics resulted in a peak list with 2329 features. Metabolome differed based on Bd exposure (PERMANOVA F(1, 38) = 4.50, R
2 = 0.11, p = 0.003). Wilcoxon’s tests identified 729 features of which 37 had a log2 fold change > 1.5 and q < 0.05 (Supplementary Data 19). PLS-DA produced separation based on Bd exposure with 111 metabolite features having a VIP score > 2, of which 24 had a log2 fold change > 1.5 and q < 0.05 (Supplementary Data 20). Screening for anti-Bd metabolites found an increased abundance in the Bd-exposed group of a metabolite feature putatively annotated as indole-3-carboxaldehyde (q < 0.05).
Microbiome/metabolome integration using DIABLO [48] identified bacterial and metabolite features discriminating control and exposure treatments (Fig. 6a,b, Supplementary Data 21). Bd exposure was associated with enrichment of bacterial ASVs and metabolite features (Fig. 6c, d) indicating that infection drives changes in microbe-metabolite interactions. All three discriminatory bacterial ASVs (Acinetobacter, Sphingobacterium, Comamonas) identified using ALDEx2 were also identified by DIABLO.
Biomarkers in the laboratory and field
We identified common biomarkers in the laboratory and field for bacterial 16S, shotgun and metabolome, but not fungal ITS2 data (Supplementary Data 22). Three bacterial ASVs (Sphingobacterium, Stenotrophomonas, Comamonadaceae) were discriminatory in the laboratory and field, with increased abundance in the epizootic population but reduced abundance in the Bd-exposed experimental group (in which all animals survived), indicating an association with negative health outcomes. Of these taxa, Stenotrophomonas abundance was significantly positively correlated with Bd infection intensity on day 30 of the experiment (SI Fig 6). Eighteen bacterial KOs were discriminatory in both the laboratory and field. Six metabolite annotations from the laboratory and field were identical including indole-3-carboxaldehyde, a bacterial-derived anti-Bd metabolite [12] that was increased in the epizootic population and in the Bd-exposed laboratory treatment.