Probiotic colonization of Xenopus laevis skin causes short-term changes in skin microbiomes and gene expression

Abstract

Probiotic therapies have been suggested for amelioration efforts of wildlife disease such as chytridiomycosis caused by Batrachochytrium spp. in amphibians. However, there is a lack of information on how probiotic application affects resident microbial communities and immune responses. To better understand these interactions, we hypothesized that probiotic application would alter microbial community composition and host immune expression in Xenopus laevis. Accordingly, we applied three amphibian-derived and anti-Batrachochytrium bacteria strains (two Pseudomonas spp. and one Stenotrophomonas sp.) to X. laevis in monoculture and also as a cocktail. We quantified microbial community structure using 16S rRNA gene sequencing. We also quantified genes involved in X. laevis immune responses using quantitative reverse transcriptase polymerase chain reaction (RT-qPCR) and skin transcriptomics over 1 and 3-week periods. All probiotic treatments successfully colonized X. laevis skin for 3 weeks, but with differential amplicon sequence variant (ASV) sequence counts over time. Bacterial community and immune gene effects were most pronounced at week 1 post-probiotic exposure and decreased thereafter. All probiotic treatments caused initial changes to bacterial community alpha and beta diversity, including reduction in diversity from pre-exposure anti-Batrachochytrium bacterial ASV relative abundance. Probiotic colonization by Pseudomonas probiotic strain RSB5.4 reduced expression of regulatory T cell marker (FOXP3, measured with RT-qPCR) and caused the greatest gene expression changes detected by transcriptomics. Single bacterial strains and mixed cultures, therefore, altered amphibian microbiome-immune interactions. This work will help to improve our understanding of the role of the microbiome-immune interface underlying both disease dynamics and emergent eco-evolutionary processes.IMPORTANCEAmphibian skin microbial communities have an important role in determining disease outcomes, in part through complex yet poorly understood interactions with host immune systems. Here we report that probiotic-induced changes to the Xenopus laevis frog skin microbial communities also result in significant alterations to these animals' immune gene expression. These findings underscore the interdependence of amphibian skin immune-microbiome interactions.

Keywords: amphibian; disease ecology; immune response; metagenomics; microbial ecology; microbiome; probiotics; transcriptomics.

Fig 1

Probiotic persistence. (A) Heatmap abundances of ASV sequence counts per sample by week 0 (day 0; prior to probiotic inoculation), week 1 (day 8), week 2 (day 15), and week 3 (day 22). Sample range by color shade is given in the corresponding legend. Scatterplot of probiotic sequence counts for each probiotic treatment and corresponding linear mixed-effects model equation for (B) ASV18 (Pseudomonas RSB5.4 [P1]), (C) ASV9 (Stenotrophomonas THA2.2 [P2]), (D) ASV10 (Pseudomonas tolaasii RSB5.11, P3), and (E) the three-probiotic cocktail.

Fig 2

Line plots of anti-Bd bacterial ASV richness, by treatment and timepoint. All treatment groups at the beginning of the experiment pre-inoculation, were similar at day 0 (D0). All treatments at week 1 post-inoculation (W1), week 2 post-inoculation (W2), and week 3 post-inoculation (W3) also had similar or slightly decreasing anti-Bd bacterial ASV richness over time, with P4 having slightly higher richness (data shown is excluding observed probiotics by respective individuals at each timepoint). Treatments are matched by color and are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas RSB5.11, and P4 = cocktail. Error bars at each point represent plus or minus one standard deviation.

Fig 3

Line plots of corrected ASV relative abundance of anti-Bd bacteria, by treatment and timepoint. All treatment groups exhibited a lower relative abundance of anti-Bd bacteria over time as compared to the control (data shown is excluding observed probiotics by respective individuals at each timepoint). Treatments are matched by color and are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas tolaasii RSB5.11, and P4 = cocktail. Time points are given for each treatment and are coded as follows: D0 = day 0, W1 = week 1, W2 = week 2, and W3 = week 3. Error bars at each point represent plus or minus one standard deviation.

Fig 4

Principal coordinate analyses (axes 2 and 3; shown to visualize treatment effects) of microbial community beta diversity using the Bray-Curtis dissimilarity metric, by time point. Data ellipses show an 80% confidence and assume a multivariate t-distribution. The treatment group is given by color and shape in the legend. (A) All treatment groups at the beginning of the experiment, pre-inoculation. (B) All treatments at week 1 post-inoculation. (C) All treatments at week 2 post-inoculation. (D) All treatments at week 3 post-inoculation. Treatments are matched by color across panels and are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas tolaasii RSB5.11, and P4 = cocktail.

Fig 5

Boxplots of 2^−ΔΔCT GAPDH normalized expression data at week 1 for (A) CSF1, (B) IL10, (C) FOXP3, and (D) TNFA. The expression (y-axis) scale has been log₁₀ normalized for purposes of visualization only. Bars with * indicate a significant difference of P < 0.05. Treatments are matched by color across panels and are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas tolaasii RSB5.11, and P4 = cocktail.

Fig 6

Upset plot showing overlap between each set of differentially expressed immune-related genes. Horizontal bars show the size of each set, and vertical bars show the size of each intersection between sets (denoted by black points within the central grid). The rows of the central grid are colored by test treatment. Treatments are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas tolaasii RSB5.11, and P4 = cocktail.

Fig 7

An association network of proportional abundance of bacterial genera and gene expression modules across all treatments (A). Nodes are colored by data type (bacterial genus/expression module) and their size is scaled by betweenness centrality. Edge color indicates direction of association (red = negative, blue = positive). The heatmap shows correlations (red = negative, blue = positive) between gene expression modules (B) and bacterial genera. Matrix color shows direction (as in panel A) and strength of correlation; non-significant correlations are set to zero (white). Bacteria are ordered by phylum. Vertical lines highlight the modules with the highest proportion of genes annotated with the “immune system process” GO term.