2.1 Alpha diversity and the core microbiota
The ACE, Chao1, Shannon, and Sob indices showed a tendency to first decrease and then increase. The ACE, Chao1, and Shannon indices differed among the control, gentamicin, and recovery groups (Kruskal-Wallis H test; P < 0.05). The ACE, Chao1, and Shannon indices differed between the antibiotic bath group and control group (Wilcoxon rank-sum test; P < 0.05, Fig. 1). However, there were no differences in the ACE, Chao1, Shannon, and Sob indices between the antibiotic bath group and recovery group (Wilcoxon rank-sum test; P < 0.05, Fig. 1). The rarefaction curves tended to plateau, indicating that the amount of sampling is reasonable and that more sampling produces only a small number of new OTUs (Fig. S1).
The Venn diagrams showed that the number of unique OTUs in the gentamicin bath group (21) was lower than that in the recovery group (60) and control group (119) (Fig. S2). As the number of samples increased, the number of core OTUs in the control (C) group decreased to a lesser extent, while that in the gentamicin group decreased to a more significant extent (Fig. S3a). The core OTU numbers in the control, gentamicin, and recovery groups and all 21 frogs were 47, 25, 39, and 16, respectively (Fig. S3b). The 16 core OTUs were from four phyla, including 8 from Bacteroidetes, 3 from Actinobacteria, 3 from Proteobacteria, and 2 from Firmicutes (Table S1). The most abundant core OTUs were OTU198 (Vagococcus), OTU179 (Citrobacter), OTU417 (Bacteroides), OTU38 (Bacteroides), OTU333 (Faecalitalea), and OTU753 (Arthrobacter) (Table S1).
2.2 Beta diversity
The gut bacterial community composition differed between the G and C groups and between the R and C groups based on the Bray-Curtis dissimilarity matrix (Adonis: P < 0.05; ANOSIM: P < 0.05; Table 1 and Fig. 2) and the unweighted UniFrac distances (Adonis: P < 0.05; ANOSIM: P < 0.05; Table 1 and Fig. 2). However, the gut bacterial community composition did not significantly differ between the G and R groups based on the Bray-Curtis dissimilarity matrix (Adonis: P > 0.05; ANOSIM: P > 0.05; Table 1 and Fig. 2) and unweighted UniFrac distances (Adonis: P > 0.05; ANOSIM: P > 0.05; Table 1 and Fig. 2). The gut bacterial community remarkably split into two major groups on the NMDS plot that corresponded to the C group and the G + R group (Fig. 2).
Table 1
Pairwise comparisons showing differences in the gut bacterial community in different groups.
|
Bray-Curtis
|
Unweight UniFrac
|
|
ANOSIM Adonis
|
ANOSIM
|
Adonis
|
R
|
P
|
R2
|
P
|
R
|
P
|
R2
|
P
|
C vs. G
|
0.489
|
0.002
|
0.214
|
0.006
|
0.249
|
0.023
|
0.146
|
0.019
|
G vs. R
|
0.022
|
0.329
|
0.081
|
0.362
|
0.095
|
0.174
|
0.094
|
0.192
|
C vs. R
|
0.338
|
0.010
|
0.164
|
0.018
|
0.270
|
0.037
|
0.154
|
0.033
|
All
|
0.284
|
0.005
|
0.204
|
0.005
|
0.207
|
0.018
|
0.170
|
0.015
|
C indicates the control group, G indicates the gentamicin group, and R indicates the recovery group. |
2.3 Composition of and differences in the gut microbiota
The taxonomic assignment analysis showed that the most abundant phyla in the control, gentamicin, and recovery groups were Firmicutes (C: 38.19%, G: 15.90%, R: 30.09%), Bacteroidetes (C: 35.15%, G: 31.00%, R: 40.94%), Proteobacteria (C: 21.96%, G: 44.74%, R: 23.44%), and Actinobacteria (C: 3.23%, G: 5.37%, R: 3.38%) (Figs. 3a and S4). In total, 10 phyla were shared among all groups, and no bacterial phylum was found to be significantly different in abundance in the control, gentamicin, and recovery groups (Kruskal–Wallis H test, FDR correction, CI: Scheffer, P > 0.05).
At the genus level, the most abundant microbial genera were Bacteroides, Morganella, Vagococcus, Faecalitalea, Parabacteroides, Arthrobacter, Alistipes, Pseudomonas, and Myroides (Fig. 3b and S5). Of all 290 genera, the 5 bacterial genera Crenobacter, Morganella, unclassified_f__Eggerthellaceae, unclassified_f__Veillonellaceae, and Weissella exhibited significant differences among the control, gentamicin, and recovery groups (Kruskal–Wallis H t est, FDR correction, CI: Scheffer, P < 0.05; Fig. S5).
The LEfSe showed that Fusobacteria were significantly enriched in the control group (LDA > 2, P < 0.05, Fig. 4a). This bacterial taxon with significant differences accounted for most of the control group, while the gentamicin and recovery groups had fewer bacterial taxon (Fig. 4a). The LEfSe analysis at the genus level revealed that Morganella, CL500_29_marine_group, Paenarthrobacter, and Plesiomonas were significantly enriched in the gentamicin group and that Butyricicoccus, Corynebacterium_1, Enterococcus, Phascolarctobacterium, Providencia, Vagococcus, and Weissella were significantly enriched in the recovery group (LDA > 2, P < 0.05, Fig. 4a). When LDA > 4, the LEfSe analysis showed that at the genus level, Citrobacter (C group), Morganella (G group), and Vagococcus (R group) were significantly enriched (LDA > 4, P < 0.05, Fig. 4b).
2.4 Indicator taxa of frog gut dysbiosis and potentially pathogenic genera
At the genus level, forty-four indicator species were selected, including 42 species indicating the C group and 2 species (Butyricicoccus and Morganella) indicating the G + R group (Fig. 5). A heat map depicting the normalized abundances of the 44 indicator taxa across the samples was generated and showed their abilities to discriminate among the samples according to the sample grouping process (Fig. 5).
The distribution and comparison of potentially pathogenic genera in the gut of the control, gentamicin, and recovery groups are shown in Table 2 and Fig. S7. Aeromonas, Acinetobacter, and Chryseobacterium differed among the control, gentamicin, and recovery groups (Kruskal-Wallis H test; P < 0.05). The relative abundance of the bacterial genera belonging to Aeromonas, Citrobacter, and Chryseobacterium was significantly decreased in the G group (Wilcoxon rank-sum test; P < 0.01, Table 2). After 7 days of recovery, Aeromonas, Citrobacter, and Chryseobacterium still significantly differed between the C and R groups (Wilcoxon rank-sum test; P < 0.01, Table 2).
Table 2
Changes in the relative abundance of potentially pathogenic genera after antibiotic baths (Wilcoxon rank-sum test).
Genus
|
OTU
|
C group
(mean ± SD)
|
G group
(mean ± SD)
|
R group
(mean ± SD)
|
Acinetobacter
|
OTU666, OTU528, OTU332, OTU41, OTU175, OTU173 OTU160
|
1.51 ± 3.06a
|
0.14 ± 0.18a
|
0.09 ± 0.12a
|
Aeromonas
|
OTU2
|
1.25 ± 1.98a
|
0b
|
0b
|
Citrobacter
|
OTU179, OTU74
|
12.70 ± 11.14a
|
0.85 ± 1.71b
|
3.16 ± 3.82b
|
Chryseobacterium
|
OTU413
|
0.28 ± 0.41a
|
0.01 ± 0.01b
|
0b
|
Proteus
|
OTU405
|
0.01 ± 0.01a
|
0.01 ± 0.02a
|
0a
|
Pseudomonas
|
OTU192, OTU552, OTU185
|
3.06 ± 3.51a
|
1.98 ± 3.81a
|
2.71 ± 5.04a
|
Staphylococcus
|
OTU222, OTU733
|
0.01 ± 0.01a
|
0.02 ± 0.04a
|
0.03 ± 0.03a
|
Streptococcus
|
OTU346, OTU322
|
0.01 ± 0.02a
|
0.01 ± 0.0 a
|
0a
|
2.5 Relationship between the bacterial community structure and function
Three hundred functional pathways were obtained in the C and G + R groups. The principal coordinate analysis did not show significant differences in the functional composition between the C and G + R groups (Fig. 6a). Similarly, the bacterial community similarity test did not show a significant difference in the functional composition between the C and G + R groups (C:G + R, ANOSIM statistic R = -0.041, P = 0.637, Fig. 6a).
The linear regression analysis showed that the gut microbiota composition and functional composition of the C and G + R groups were not significantly and positively correlated (C and G + R: r = 0.125, P = 0.079, Fig. 6b), indicating that changes in the gut microbiota of R. dybowskii did not alter bacteria-mediated physiological functions. Significant differences were observed in six KEGG pathways (cancers: overview, circulatory system, environmental adaptation, excretory system, infectious diseases: bacterial, and substance dependence) among the control, gentamicin, and recovery groups (Kruskal-Wallis H test, P < 0.05; Fig. S8).