Madison JD et al. (2025), Probiotic colonization of Xenopus laevis skin c...

XB-ART-61325

Infect Immun 2025 Apr 13;935:e0056924. doi: 10.1128/iai.00569-24.

Show Gene links Show Anatomy links

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

Madison JD , Osborne OG , Ellison A , Garvey Griffith CN , Gentry L , Gross H , Gratwicke B , Grayfer L , Muletz-Wolz CR .

???displayArticle.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.

???displayArticle.pubmedLink??? 40172536
???displayArticle.pmcLink??? PMC12070741
???displayArticle.link??? Infect Immun
???displayArticle.grants??? [+]

Species referenced: Xenopus laevis
Genes referenced: csf1 foxp3 gapdh il10 il34 tlr2 tlr6

???attribute.lit??? ???displayArticles.show???

	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 log10 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.
	Figure S1. Coverage given (pre-rarefaction) with sequence count plotted against richness. Each line represents a sample. Sample type is given in the legend. C = control, P1 (Pseudomonas RSB5.4), P2 (Stenotrophomonas THA2.2), P3 (Pseudomonas RSB5.11), and P4 (cocktail) are the corresponding probiotic treatments.
	Figure S2. Bar plots of ASV richness (rarefied) by treatment and timepoint. A slight increase in ASV richness was observed over time across treatments and control. Treatments are matched by color/shape and are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas tolaasii RSB5.11, P4 = cocktail. Timepoints 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.
	Figure S3. Principal coordinate analyses (axes 1 and 2) of microbial community beta diversity using the BrayCurtis dissimilarity metric. Treatment group is given by color and time is given by shape, as indicated in the legend. Data ellipses show an 80% confidence and assume a multivariate t-distribution. Treatments are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas RSB5.11, P4 = cocktail.
	Figure S4. Principal coordinate analyses (axes 2 and 3) of microbial community beta diversity using the Jaccard dissimilarity metric. Treatment group is given by color and time is given by shape, as indicated in the legend. Data ellipses show an 80% confidence and assume a multivariate t-distribution. Treatments are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas RSB5.11, P4 = cocktail.
	Figure S5. Upset plot showing overlap between each set of all differentially expressed 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 as in figures 1–7. Treatments are coded as follows: C = no-probiotic control, P1 = Pseudomonas RSB5.4, P2 = Stenotrophomonas THA2.2, P3 = Pseudomonas RSB5.11, P4 = cocktail.
	Figure S6. A heatmap showing the proportion of genes in each expression module annotated with each of the level one biological process GO terms. Horizontal lines highlight the “immune systems process” GO term and vertical lines highlight the modules with the highest proportion of genes annotated to it.
	Figure S7. A barplot showing a subset of differentially expressed genes of interest (as compared to control) at weeks 1 and 3 by respective treatment group (P1–P4) and control (C). This subset of genes includes IgJ in (A) and (B), IL12B in (C) and (D), IL17C in (E) and (F), and LCK in (G) and (H). All expression data is normalized by library size. The left column shows week 1 and the right column shows week 3 for the corresponding gene. 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 RSB5.11, P4 = cocktail.
	Figure S7. A barplot showing a subset of differentially expressed genes of interest (as compared to control) at weeks 1 and 3 by respective treatment group (P1–P4) and control (C). This subset of genes includes IgJ in (A) and (B), IL12B in (C) and (D), IL17C in (E) and (F), and LCK in (G) and (H). All expression data is normalized by library size. The left column shows week 1 and the right column shows week 3 for the corresponding gene. 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 RSB5.11, P4 = cocktail.
	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 log10 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.