3.1 Overview of alpha diversity and the core microbiota
Rana amurensis and R. dybowskii were sampled for sequencing, resulting in 2,496,447 high-quality sequences with a mean of 43,042 sequences per sample. A total of 2,607 OTUs were obtained with an average length of 442.74 bp per read. The average number of OTUs per sample was 325.78, with a range of 72 (sample dyL02) to 1,140 (sample dyS02). Rarefaction analysis showed that the sequenced samples mostly reached the plateau phase, particularly the amL and amA samples (Fig. S1).
There were no significant differences in the Sob index between groups (Wilcoxon rank sum test, P > 0.05, Fig. 1a). The Shannon index differed significantly between the amL and dyL groups and between the amA and amS groups (Wilcoxon rank sum test, P < 0.05; Fig. 1b). The PD index differed significantly between the dyL and dyS groups (Wilcoxon rank sum test, P < 0.05) (Fig. 1c).
With respect to the core microbiota of all frogs, there were 8 core OTUs (found in > 90% of all frogs) (Fig. 1d and 2). Twenty and 7 species-specific core OTUs were identified in R. amurensis and R. dybowskii, respectively (Figs. 2). Most notable were Arthrobacter, Chryseobacterium, Erysipelatoclostridium, Hafnia-Obesumbacterium and Pseudomonas, as they were widely distributed between species, habitats and fasting times and were dominant in abundance (Fig. 1d and 2).
3.2 Effects of captive and wild environments and host species on the gut microbiota
The effects of host and habitat on the gut microbiota were assessed using two frog species that hibernated in different habitats, and the results showed that the overall composition and structure of the gut microbiota were strongly associated with the effects of host species (two-way PERMANOVA: Bray-Curtis, P = 0.001; unweighted UniFrac, P = 0.001) and habitat (two-way PERMANOVA: Bray-Curtis, P = 0.003; unweighted UniFrac, P = 0.003), yet no significant effects were detected for their interactions (two-way PERMANOVA: Bray-Curtis, P = 0.079; unweighted UniFrac, P = 0.130) (Table 2; Fig. 3a and b).
Table 2
Results of PERMANOVA examining the effects of host, habitat, and host × habitat interaction on two beta diversity metrics (Bray-Curtis and unweighted UniFrac distance matrices)
| Bray-Curtis | Unweighted UniFrac unweighted UniFrac |
| F | R2 | P | F | R2 | P |
Host | 3.933 | 0.075 | 0.001 | 4.253 | 0.080 | 0.001 |
Habitat | 3.055 | 0.058 | 0.003 | 3.388 | 0.064 | 0.003 |
Host × habitat | 1.546 | 0.029 | 0.079 | 1.431 | 0.027 | 0.130 |
The gut microbiota differed between host species. The composition and structure of the gut microbiota of the same host (R. amurensis and R. dybowskii) differed significantly between the captive and wild environments (R. amurensis, Adonis: Bray-Curtis, R2 = 0.091, P = 0.022; unweighted UniFrac, R2 = 0.088, P = 0.037; R. dybowskii, Adonis: Bray-Curtis, R2 = 0.094, P = 0.014; unweighted UniFrac, R2 = 0.109, P = 0.015) (Fig. 3a and b; Table S1).
The gut microbiota differed between captive and wild environments. In natural habitat, the composition and structure of the gut microbiota differed significantly between the two species (Adonis: Bray-Curtis, R2 = 0.121, P = 0.001, unweighted UniFrac, R2 = 0.091, P = 0.011; Table S1). In captive habitat, there were significant differences in the composition and structure of the gut microbiota based on unweighted UniFrac distances (Adonis: unweighted UniFrac, R2 = 0.142, P = 0.004), but there were no significant differences based on Bray-Curtis distances (Adonis: Bray-Curtis, R2 = 0.081, P = 0.058) (Fig. 3a and b; Table S1).
The dominant phyla in the gut microbiota of both species in both habitats ranked by abundance (> 5%) were Bacteroidetes, Proteobacteria, Firmicutes, Actinobacteria and Spirochaetes, accounting for 96.34 ± 0.65% (amL), 96.59 ± 0.52% (amS), 92.81 ± 0.52% (dyL) and 98.06 ± 0.67% (dyS) of the total phyla, respectively (Figs. 4a and S2). The Wilcoxon rank sum test showed that host had a significant effect on Actinobacteria (amL vs. dyL, P < 0.05; Table S2) and Spirochaetes (amS vs. dyS, P < 0.05) and that habitat had a significant effect on Actinobacteria (dyL vs. dyS, P < 0.05; Table S2). Pseudomonas, Bacteroides, Chryseobacterium, Erysipelatoclostridium and unclassified_f__Erysipelotrichaceae were the most abundant genera among the amL, dyL, amS, and dyS groups (Figs. 4b and S3). Of the 745 genera, 116 (e.g., Bacteroides, Chryseobacterium and Deefgea) were significantly different (Kruskal-Wallis H test, P < 0.05; Table S3).
The LEfSe results revealed significant differences in the gut bacterial taxa, mainly in the phyla Actinobacteria (dyL group), Bacteroidetes (amS group) and Spirochaetes (amL group) (LDA > 4, p < 0.05; Fig. 5a). The LEfSe results at the genus level indicated that 5 of the 745 genera present in the dataset were differentially abundant in the two species between different habitats. Flavobacterium, unclassified_f__Comamonadaceae and unclassified_f__Flavobacteriaceae were most abundant in the dyL group, Acinetobacter was most abundant in the amS group, and Deefgea was most abundant in the dyS group (LDA > 4, P < 0.05; Fig. S4).
3.3 The gut microbiota of frog entering captivity environments occurred more significant changes
The gut microbiota differed with fasting time in wild environments. In R. amurensis, the composition and structure of the gut microbiota differed significantly between the early and late stages of hibernation (Adonis: Bray-Curtis, R2 = 0.179, P = 0.001; unweighted UniFrac, R2 = 0.145, P = 0.014; Fig. 3a and b). In R. dybowskii, the composition of the gut microbiota differed significantly between the early and late stages of hibernation according to Bray-Curtis distances (Adonis: Bray-Curtis, R2 = 0.101, P = 0.032; Fig. 3a) but not unweighted UniFrac distances (Adonis: unweighted UniFrac, R2 = 0.090, P = 0.084; Fig. 3b).
The gut microbiota differed with fasting time in captive environments. In R. amurensis, the composition and structure of the gut microbiota differed significantly between the early and late stages of hibernation (Adonis: Bray-Curtis, R2 = 0.198, P = 0.005; unweighted UniFrac, R2 = 0.147, P = 0.022; Fig. 3a and b). In R. dybowskii, the composition of the gut microbiota differed significantly between the early and late stages of hibernation (Adonis: Bray-Curtis, R2 = 0.273, P = 0.001; unweighted UniFrac, R2 = 0.267, P = 0.016) (Table S4; Fig. 3a and b).
The Bray-Curtis dissimilarities and unweighted UniFrac distances between the early and late stages of hibernation were greater in captive environments than in wild environments (Fig. 3a and b). Based on Bray-Curtis and unweighted UniFrac distance and comparison of the pairwise distance of the gut microbiota between the early and late stages of hibernation of captive and wild, the pairwise distances of the early and late stages of hibernation between captive and wild were significant (Wilcoxon rank-sum test, P < 0.05; Fig. 5a).
We used Venn analysis to pairwise analysed the shared OTUs in each sample in the early of hibernation (amA and dyA) and each sample in the late stages of hibernation (amL, dyL, amS, and dyS groups). The average of the shared OTUs in amA and amL groups was 78.64, and that of the shared OTUs in amA and amS groups was 116.90. There were significant differences between the two groups in the average of the shared OTUs (ANOVA, F = 119.425, df = 1, P < 0.01; Fig. 5b). The average of the shared OTUs in dyA and dyL groups was 90.16, and that of the shared OTUs in dyA and dyS groups was 115.10. There were significant differences between the two groups in the average of the shared OTUs (ANOVA, F = 10.818, df = 1, P < 0.01; Fig. 5b).
3.4 The best variables for predicting microbiota structure as selected by the RDA model
Constrained ordination (RDA) was used to explore the potential relationships between microbial communities and the environment, fasting time, and host species. Three impact factors with VIF values below 3 were included in the RDA model, including host (VIF: 1.001), fasting (VIF: 2.5385), and habitat (VIF: 2.386). The variation partitioning analysis (VPA) results showed that 5.98% of the variance in gut microbiota structure could be explained by the host (3.94%), fasting (1.43%), and habitat (1.33%) (Fig. 6a). However, the AIC-based model results suggested that host and habitat were the variables that needed to be included in the RDA to explain the variance in taxa (AIC: -25.778, P < 0.05; Fig. 6b).