Eighty adult treefrogs (Polypedates megacephalus) were collected from New Taipei City and Taichung City, Taiwan. All animals were housed in 240-l glass tanks at 23°C under a 8:16 h light:dark cycle. Turkestan cockroach nymphs (Finke, 2013) were fed to the treefrogs at a quantity of 10% of treefrog biomass twice a week. The treefrogs were acclimatized for 3 months prior to experiments.

Fecal samples were collected from treefrog guts within 20 min after euthanasia. To avoid cross-contamination, each sample was collected using a fresh pair of sterile tweezers. The contents of each gut were emptied into a sterile vial and immediately stored at –80°C. DNA was subsequently extracted with the QIAamp DNA Stool Mini Kit (Qiagen, Valencia, CA, USA). The V4 region (292 bp) of 16S rRNA gene was PCR-amplified with 515F (5′−GTGCCAGCMGCCGCGGTAA−3′) and 12-base barcoded 806R (5′−GGACTACHVGGGTWTCTAAT−3′) primers (Caporaso et al., 2012). Following PCR, samples were gel extracted with the NucleoSpin Gel Extraction kit (Macherey-Nagel, Germany). The purified samples were pooled in equal concentrations and sequenced using an Illumina MiSeq (Illumina, San Diego, CA, USA) with a V2 PE500 cartridge (500 cycles). All datasets have been deposited in GenBank under the BioProject ID PRJNA341914 and BioSample ID SAMN05730167 and SAMN05730170.

All paired-end sequences were merged by FLASH (Magoč and Salzberg, 2011), and all merged sequences were further analyzed by conducting the Quantitative Insights Into Microbial Ecology (QIIME) pipeline (Caporaso et al., 2010b). In the QIIME analysis pipeline, the low-quality sequences (sequences that were <200 bp in length, had a quality score <25, contained ambiguous characters, had an unreadable barcode, or did not contain the primer sequence) were removed using the USEARCH quality filter (Edgar, 2010). The UCLUST (Edgar, 2010) function in QIIME was used to cluster the remaining sequences, with a minimum coverage of 99% and minimum sequence identity of 97%. The longest sequences from each phylotype were selected to perform sequence identification, and PyNAST (Caporaso et al., 2010a) and UCLUST were selected to perform the sequence alignment and taxonomy assignment.

A concept for improvement of the investigation of the functional roles of gut microbiota is to explore the microbial world from a simple to complex state, from an initial stage to homeostasis. Hibernation is considered as a survival strategy designed to conserve energy when conditions are harsh. During hibernation, animals can adapt to temperature perturbations and extend their lifespan by slowing their heartbeat, reducing metabolic activity and their energy requirements (Book, 1974). At the same time, most microorganisms in the gut of a host decline in numbers due to extreme temperature and low nutrient supply (Gossling et al., 1982a). Hibernation of treefrogs, therefore, provides a natural model for monitoring microbial growth starting from inoculation (i.e., the beginning of a developing microbial interaction network) to homeostasis, and could provide an opportunity to shed light on our understanding of how microorganisms form their interaction networks.

To stimulate AH, we modified a program previously described in leopard frogs (Rana pipiens) (Gossling et al., 1982a). Fasting is known to reduce microbial complexity in amphibians, as well as fish, reptiles, birds, and mammals (Gossling et al., 1982b; Sonoyama et al., 2009; Costello et al., 2010; Kohl et al., 2014). Therefore, prior to stimulation of AH, food was withheld from the treefrogs kept in the same housing for 7 days to reduce overall diversity. All fasting treefrogs were then transferred into an incubator kept in constant darkness. AH was stimulated initially at 21°C with a relative humidity of 90%. During day 1, the temperature was maintained. During days 2–3, the temperature was gradually reduced to 4°C. Thereafter, treefrogs were housed in the same manner for 7 days. During days 11–12 (i.e., after 7 days of AH described above), the temperature was gradually increased to 21°C, and the treefrogs returned to their active status. Each active treefrog was fed with Turkestan cockroach nymphs (~0.25–0.3 g). Post-feeding samples were collected within 2.5 days. All protocols were approved by the Academia Sinica Biosafety Committee and Institutional Animal Care and Utilization Committee.

To collect time series samples for generating the microbial interaction network, we collected the gut contents of the treefrogs over a total of 12 time points spanning 15 days over the AH period, including day 1 (fasting for 7 days before AH), day 11–11.75 (every 6 h), and days 12–15 (every 12 h). For each time point, three to five treefrogs were euthanized. Details of treefrog body mass, sampling size, and time points of fecal collection over the AH period are presented in Supplementary Table 1.

We used the Metagenomic Microbial Interaction Simulator (MetaMIS) with a user friendly interface (Shaw et al., 2016) to infer microbial relations by introducing a single time-series dataset of microbial composition containing 12 time points. MetaMIS is a tool based on the generalized Lotka-Volterra model (Bhargava, 1989), and designed to infer underlying microbial interactions according to metagenomic abundance profiles. Lotka-Volterra equations have been widely used to infer animal interactions in dynamic systems, and recently have been applied to reveal microbial interacting relationships between operational taxonomic units (OTUs). The detailed algorithms and equations were described by Bucci et al. (Stein et al., 2013). Due to the limitation of computing power (Interl® Core™ i7-4770 CPU @3.40 GHz processor and 32 Gb RAM), it is not feasible to infer interaction network at the species level. Therefore, we assumed that in general, all species in a genus identified in our amplicon study share similar functions. According to the compositional profiles among 12 time points at the genus level, MetaMIS can systematically examine interaction patterns, such as mutualism or competition; only the top 25% of interacting relationships were considered in this study. We used the mfinder tool (Milo et al., 2002) to identify significant 3-node directed motifs that contain two directed edges pointing to the same node. This process was repeated 100 times and all relationships that passed the criteria described above at least 50 times (permutation cutoff > 0.5) were considered to be reliable interacting relationships.

There are several genera of bacteria characterized as probiotic in amphibians (Pasteris et al., 2009a,b; Dias et al., 2010; Becker et al., 2011; Mendoza et al., 2012; Bletz et al., 2013). LAB have been considered as the major probiotics for the treatment of raniculture (Dias et al., 2010; Mendoza et al., 2012). In this study, L. garvieae was selected as a probiotic strain because Lactococcus was present in at least 85% of all the samples with a relative abundance of 0.14 ± 0.39%, and L. garvieae has been identified as a biological control agent in bullfrogs (Mendoza et al., 2012). According to the interacting relationships inferred by MetaMIS, Corynebacterium showed beneficial effects on Lactococcus, while Bacillus was found to have a beneficial interacting relationship with Corynebacterium. We decided to use three representative species of the three mentioned genera for further validation. Therefore, L. garvieae, C. variabile, and B. coagulans were chosen to validate the inferred interacting relationships. These three species were selected by their abundance ranking within the genus according to 16S rRNA amplicon data. We reasoned that if the inferred interacting relations were correct, the cooperative combination would yield a higher probiotic concentration than using a single strain.

To prepare for oral administration, L. garvieae and B. coagulans were grown in Tryptic Soy broth at 30 and 55°C respectively. C. variabile was grown in Brain Heart Infusion broth at 30°C. All strains were grown overnight with agitation in a shaking incubator. After incubation, the cells were harvested by centrifugation at 2,500 × g for 20 min at 4°C. The cell pellets were washed twice with 0.9% saline and resuspended using the same buffer. The measured population level of bacteria in the test diet was 10⁷ CFU g⁻¹ (colony-forming unit).

Treefrogs were divided into five groups (N = 8 per group), and acclimated for 1 week before the start of the trial. The trial was conducted for a 2-week period. Each test group of treefrogs was dosed with L. garvieae (G₁), C. variabile (G₂), or B. coagulans (G₃) singly, or with a combination of the three strains (G₄) once per day by direct oral gavage. The control (G₅) was fed with 0.9% saline during the entire trial period. After the trial, four fecal samples were collected from four treefrogs in each group (Supplementary Table 2) for further 16S rRNA amplicon analysis, and the other four treefrogs were used to perform quantitative PCR for interleukin 10 (IL-10) expression as described in the following section.

In order to test the immune response after oral administration with probiotic, the expression level of IL-10 was measured after LPS (lipopolysaccharide) stimulation (Qi et al., 2015). To characterize the change of treefrog IL-10 expression after LPS stimulation, four treefrogs in each group were injected intraperitoneally (i.p.) with LPS (150 μg/100 g body weight). Animals were anesthetized and euthanized 24 h after injection. All fecal contents were removed and the remaining tissue samples of gut were collected (weight range 0.03 to 0.08 g) (Supplementary Table 3). All tissue samples were homogenized and total RNA was extracted from homogenized samples using Trizol reagent (Invitrogen, USA), quantified using a Nanodrop–1000 spectrophotometer, and reverse transcribed into cDNA using the Superscript II reverse transcription system (Invitrogen, USA) according to the manufacturer's instructions. Quantitative real time PCR was performed using Power SYBR green PCR Mastermix (Applied Biosystems) on a real-time instrument (ABI mode 7300 Sequence Detector) in 96-well reaction plates. The reaction mixture included 10 ml of Power SYBR green PCR Mastermix, 1 μl of forward and reverse primer (10 μM each) and 1 ml of cDNA, and then brought up to a final total volume of 20 μl with ultra pure water. The sequence of IL-10 in P. megacephalus was described in previous study (Huang et al., 2016) and the forward and reverse primers were designed by NCBI Primer-BLAST (Ye et al., 2012). β-actin was used as a housekeeping control. The primer sequences for amplification of IL-10 and β-actin were as follows: (F, 5′-ACGACCCTGCTCACGTTATG-3; R, 5′-TCCGGGATGGAGTAAGAGGG-3′) and (F, 5′-GGTCGCCCAAGACATCAG-3; R, 5′-GCATACAGGGACAACACA-3′) (Hamdan et al., 2016), respectively. The relative expression of IL-10 in gut tissue samples was normalized to the expression of β−actin. The change of gene expression was expressed as fold change (log base 2) and calculated as described (Qi and Nie, 2008; Qi et al., 2010, 2015). A paired Student's t-test was applied to analyze the significance and fold change (log base 2), with a p-value less than 0.05 considered to be statistically significant.

To estimate the change of microbial complexity throughout the AH period, alpha-diversity was determined using the Shannon index, Simpson index, and the Inversed Simpson index. To determine the differences in bacterial community composition during AH, we used the Bray-Curtis similarity index (a taxonomic metric), which provides a measure of phylogenetic distance between communities from individual samples (Lozupone et al., 2007). To test the differences in richness and phylogenetic indices between time points, Student's t-test was used to determine the significance (p-value < 0.05) between time points. For comparison between each bacterial challenge and controls, we determined the fold-change in relative abundance to demonstrate a response to the stimulus relative to the background. Fold-change significantly higher than 2 or smaller than 0.5 was considered to be relevant. Differences between the two groups were analyzed for significance (p-value < 0.05) by Wilcoxon's test.

To infer microbial interactions, we focused on the level of genera. Overall, 325 genera that were characterized via the QIIME pipeline were used to construct the interacting relationships network using MetaMIS. In each run, pairwise interacting relationships with defined criteria (such as commensalism or amensalism) were generated. We obtained interacting paired relationships for 26,320 to 26,327 pairs in 100 permutations. After analysis, 10,314 inferred interaction pairs (IIPs) passed the permutation criterion (cutoff > 0.5 described in Methods), including 1,568 commensal, 1,737 amensal, 3,777 mutual, and 3,232 competitive relationships (Supplementary Table 4). For reference, we also constructed the interaction network at the species level (Supplementary Table 5). Most inferred relations are consistent in both networks generated in the level of genera and species.

To choose probiotic bacteria for the validation of IIPs, we surveyed the literature and selected three genera, i.e., Lactococcus, Lactobacillus, and Pediococcus that have been commonly used in probiotic treatments in raniculture (Dias et al., 2010; Mendoza et al., 2012). In our dataset, there were 84, 64, and 60 IIPs correlated to Lactococcus, Pediococcus, and Lactobacillus, respectively. Therefore, we selected Lactococcus as the target probiotic. In order to select the commensal partners for Lactococcus, we further surveyed the genera that directly imposed beneficial effects on Lactococcus (inferred from the interaction network). Corynebacterium was selected due to the fact it had the highest beneficial effects in the network. Additionally, Bacillus imposed beneficial effects on Corynebacterium and was selected from 69 IIPs of Corynebacterium. Consequently, in this study, we performed oral administration using three representatives of these genera, i.e., L. garvieae, C. variabile, and B. coagulans to validate the inferred interacting relations in gut microbes of the treefrogs.

To further describe probiotic networks, we identified specific IIPs for Lactococcus, Corynebacterium, and Bacillus from 1,012 IIPs to pinpoint the genera that directly or indirectly interact with our selected bacteria (Supplementary Table 6). There were 84 IIPs interacting with Lactococcus, including one commensal, 41 mutual, one amensal, and 41 competitive. From the 84 interacting partners with Lactococcus, three of them were predominant genera (relative abundance on average was larger than 5%), Bacteroides (20.8 ± 16.2%), Citrobacter (6.6 ± 10.4%), and Shewanella (16.3 ± 24.1%), and were all inferred as competing partners, suggesting that a considerable population that colonized the treefrog intestine may intrinsically inhibit the genus Lactococcus. The number of IIPs with Corynebacterium and Bacillus were 69 and 60 respectively. It is worth noting that compared with Lactococcus, Corynebacterium had more mutual relations with other intestinal bacteria; there were no amensal relationships and only 21 of 69 IIPs were competitive. A small group of microbes that interact with Lactococcus, Corynebacterium, and Bacillus is shown in Figure 2.

From the inferred interactions, 84, 69, and 60 IIPs were correlated to Lactococcus, Corynebacterium, and Bacillus, respectively. After oral administration of the three representative probiotics of these genera (L. garvieae, C. variabile, and B. coagulans), 65, 57, and 52 IIPs of Lactococcus, Corynebacterium, and Bacillus were identified in the trial by 16S rRNA amplicon analysis. The IIPs with abundance data of zero were excluded from the following analysis. Oral administration of three distinct bacterial species representative of each genus led to successful gut colonization, and led to a reasonable change in IIPs. In the G₁ treatment group (oral administration of L. garvieae), Lactococcus significantly increased in relative abundance 9.88 ± 2.1-fold compared with control (G₅ treatment group) (Supplementary Table 6). In the G₂ treatment group (oral administration of C. variabile), Corynebacterium significantly increased in relative abundance 160.51 ± 58.54-fold (data insufficient to test the change of Lactococcus). On the other hand, in the G₃ treatment group (oral administration of B. coagulans), the relative abundance of Bacillus decreased 0.59 ± 0.49-fold, with no change in Corynebacterium, and a 1.91 ± 0.49-fold increase in Lactococcus. In the G₄ treatment group (oral administration of a combination of L. garvieae, C. variabile, and B. coagulans), as expected, Lactococcus significantly increased in relative abundance 26.17 ± 6.96-fold (Wilcoxon's test, p-value < 0.05), 2.64 times greater increase than see for L. garvieae alone, reflecting the positive effects of the commensal partners. Also in the G₄ treatment group, Corynebacterium increased 95.99 ± 9.57-fold in relative abundance (Wilcoxon's test, p-value = 0.0625), while Bacillus decreased 0.09 ± 0.04-fold in relative abundance in comparison with controls, reflecting the inhibitory effect of Lactococcus. Figure 3 and Supplementary Table 7 illustrated the relative abundances of microbes that interact with Lactococcus, Corynebacterium, and Bacillus after the two-week oral trials.

The inferred network describes the potential interactions among genera, in addition to the interactions mentioned above, more information remained to be discussed (Figure 3). For example, a recent in vitro experiment showed an inhibition of L. garvieae K2 against Klebsiella pneumoniae U11468 (Olaoye, 2016), and our network inference also suggests a similar interaction between these two species. To test the whole immune response after oral probiotic administration, IL-10, an immunoregulatory cytokine involved in immune response in amphibians, was used as an index. The frog IL-10 contains conserved amino acid residues and motifs that are essential for bioactivity. The same residues have been proved to be necessary for immunostimulatory function of human IL-10. Studies have been done using IL-10 expression to examine the host immune response to bacterial infection (Qi et al., 2015). In our study, we found that IL-10 expression variation was consistent with the level of Lactococcus in G₁ to G₄ treatment groups (Table 2). For example, IL-10 expression was increased in the G₁ treatment group, while the corresponding level of Lactococcus in G₁ was upregulated by 9.88, and the change was significant compared with controls (Wilcoxon's test, p-value < 0.05). In addition, IL-10 level in the G₄ treatment group was significantly higher than in controls, also reflecting the high Lactococcus level in G₄. These results also support the idea that a cooperative combination of Lactococcus triggered significantly higher expression of IL-10 than observed when using a single strain.

To validate the inferred interaction network, the corresponding changes of all IIPs were evaluated by calculating the fold changes between the test and control groups. In the G₁ group, after oral administration of L. garvieae for 2 weeks, 52.3% (34 out of 65) of IIPs correlated with Lactococcus responded consistently with our inferred relations, including one commensal, 10 competitive, and 23 mutual interactions (Table 3 and Supplementary Table 6). There were 47.4% (27 out of 57) and 46.2% (24 out of 52) of IIPs correlated with Corynebacterium and Bacillus, respectively, which also responded consistently with our inferred relations. However, there were 23.1% (15 out of 65), 31.6% (18 out of 57), and 34.6% (18 out of 52) IIPs that correlated with Lactococcus, Corynebacterium, and Bacillus, presenting a conflicting response with our inferred relations, respectively. Overall, in the G₁ treatment group, the ratio of IIPs that responded according to the inferred relationship was 48.9% (85 out of 174), the ratio of IIPs that showed a conflicting response to the inferred relationships was 29.3% (51 out of 174) (Table 3). In the evaluations for the G₂, G₃, and G₄ treatment groups, we found 43.1–54.4% of corresponding changes of IIPs were consistent with the inferred interaction (Table 3).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.