A Robust Noninvasive Approach to Study Gut Microbiota Structure of Amphibian Tadpoles by Feces

Xiaowei SONG, Jinghan SONG, Honghong SONG, Qi ZENG and Keke SHI

1 College of Life Sciences, Xinyang Normal University, Xinyang 464000, China

2 Institute for Conservation and Utilization of Agro-bioresources in Dabie Mountains, Xinyang Normal University,Xinyang 464000, China

1. Introduction

Knowledge about the diversity and composition of gut microbiota has accumulated rapidly by cultureindependent methods, especially the 16S rDNA amplicon high-throughput sequencing technique (Suet al., 2012).However, different experimental protocols concerned with material preparation, DNA extraction and PCR primer selection can generate technical biases in drawing the gut microbiota profiles (Larsenet al., 2015; Lozuponeet al., 2013; Wagner Mackenzieet al., 2015).

Since acquiring high quality genomic DNA is the prerequisite for the downstream analysis of gut microbiota structure, researchers have performed many studies on the effect of material types, sample storage conditions and DNA extraction methods (Chooet al., 2015; Ferrandet al., 2014; Larsenet al., 2015; Songet al., 2016;Wagner Mackenzieet al., 2015). In spite of differences between intestinal and fecal microflora (Eckburget al.,2005; Navaet al., 2011), feces have been utilized as a noninvasive material to measure the intra- and interindividual variation in composition and diversity of intestinal bacterial communities (Arumugamet al., 2011;Dethlefsenet al., 2008; Yatsunenkoet al., 2012). Storage conditions for fecal samples could significantly influence microbiota profiles (Bahlet al., 2012; Chooet al., 2015;Maukonenet al., 2012), though some exceptions exist(Fouhyet al., 2015). DNA extraction protocols can produce technical variations due to different cell lysis and purification methods (Claassenet al., 2013; Maukonenet al., 2012; Wagner Mackenzieet al., 2015; Yuanet al.,2012). However, these technical variations seem not to be large enough to distort the biological variations(Lozuponeet al., 2013; Wagner Mackenzieet al., 2015).

Recent research on gut microbiota has shed new light on the associations between gut microbiota and vertebrate physiology, development and evolution (Leyet al.,2008; McFall-Ngaiet al., 2013; Nicholsonet al., 2012;Yatsunenkoet al., 2012), while the number of intestinal microflora studies in amphibians remains far smaller than that in mammals (Bletzet al., 2016; Colomboet al., 2015;Jiménez and Sommer, 2017; Kohlet al., 2013; Kohlet al.,2014; Mashoofet al., 2013; Venceset al., 2016; Wenget al., 2016; Wenget al., 2017). Intensive studies of the gut microbiota in different taxa including amphibians are essential prerequisites to elucidate host-gut microbiota symbioses, e.g., phylosymbiosis (Brookset al., 2016; Leyet al., 2008; Liet al., 2017a; Liet al., 2017b; Shapira,2016; Venceset al., 2016). In addition, the understanding of gut microbiota in amphibians could help us to take effective measures for amphibian conservation and cultivation (Jiménez and Sommer, 2017). Since many amphibians have been suffering from severe survival conditions (Hofet al., 2011), noninvasive approaches will always be an optimal choice for the study of symbioses between amphibians and gut microbiota. However, it remains unclear whether the feces can be applied as a noninvasive material to study the gut microbiota of amphibians.

In this study, we aimed to test the efficacy of a phenolchloroform method and a commercial fecal reagent kit(TIANGEN Biotech Co., Ltd.) in describing intestinal bacterial communities of Asiatic toad (B. gargarizans)tadpoles by feces. Specifically, DNA extraction quality of different methods and sample types was tested in terms of three parameters, i.e., A260/A280, A260/A230 and DNA yield rate. Furthermore, the structural consistency between bacterial communities was evaluated based on 16S rDNA amplicon high-throughput sequencing.

2. Materials and Methods

2.1. Transplantation of B. gargarizans eggsIn February-March 2016, we sampled seven broods ofB.gargarizanseggs from Xinyang City of Henan Province in China (Table S1). A ～10 cm-length chalaza was taken from each brood of eggs, and then was hatched in labs using plastic cylinders (1 L) filled with about 0.5 L dechlorinated drinking water. Subsequently, about 15 larvae in each brood were reared together with boiled green vegetable leaves rich in cellulose or fish foods rich in protein. We applied semi-natural conditions for rearing all larvae, i.e., a water change and food feeding per three days without controlling the light, humidity and temperature. All procedures used in this study were approved by the Animal Care and Use Committee of Xinyang Normal University.

2.2. Preparation of intestinal and stool samplesWhen tadpoles developed into the lower limb stage (Gosner 35–40), we collected intestinal and stool samples. We first collected 10 samples (S1–S10) of mixed fecal sediments from cultivation water by using aseptic injectors (Table S1). These mixed stool samples spin-dried in a centrifuge were applied for comparing DNA extraction quality of the phenol-chloroform method and the TIANamp Stool DNA Kit. In addition to the DNA extraction quality, 23 tadpoles (S11–S33) were further selected to compare the microbiota structure between the phenol-chloroform method and the commercial stool kit (Table S1). Tadpoles(S11–S16 and S22–S27) were individually cultivated for less than 24 hours after a water change. Then we got their spin-dried stool samples through the same approach for mixed stool samples. We sacrificed tadpoles using 75%ethanol solution before the extraction of the gut samples into sterile microcentrifuge tubes. All intestinal and stool samples were stored in a –20°C freezer before the DNA extraction.

2.3. Protocols of DNA extraction and DNA quality evaluationPhenol-chloroform method. ―The 500 μL SDS (1%) and each weighed sample were mixed in a sterile 1.5 mL microtube and bathed in the water for 10 min at 60°C. During the water bath, microtubes were overturned and blended for three times. Subsequently,30 μL EDTA (0.5 M) and 20 μL protease K (20 mg/mL)were added into each microtube, and the mixture was bathed in the water for 1 hour at the same temperature of the previous step. And then the samples were centrifuged at 12 000 rpm for 5 min at 4°C. The supernatant fraction was transferred to a new 1.5 mL microtube, and blended with an equal volume of phenol-chloroformisoamyl alcohol (25:24:1 in volume). The samples were centrifuged at 12000 rpm for 5 min at 4°C once again.The supernatant fraction was transferred to a new 1.5 mL microtube, and blended with an equal volume of chloroform-isoamyl alcohol (24:1 in volume). The samples were centrifuged at 12000 rpm for 5 min at 4°C once again. The supernatant fraction was transferred to a new 1.5 mL microtube, and blended with twofold absolute ethanol. After cooled at –20°C for 20 min, the samples were centrifuged at 12 000 rpm for 10 min at 4°C. After the removal of the supernatant fractions, the DNA precipitates were washed in 200 μL ethanol (70%)for 5 min, and then the mixtures were centrifuged at 12 000 rpm for 5 min at 4°C. This step was repeated three times. The dry DNA was dissolved in 40 μL TE buffer solution.

TIANamp stool DNA kit.―TIANamp Stool DNA Kit simplifies the DNA isolation by a fast spincolumn procedure. The protocol recommended by the manufacturer was utilized to extract the metagenomic DNA.

DNA quality evaluation.―The A260/A280, A260/A230 and DNA yield rate were applied for the DNA quality evaluation. The A260/A280, A260/A230 and concentration of DNA products were determined by using NanoVue Plus Spectrophotometer (GE Healthcare Inc.,Germany). The DNA yield rate was given by the ratio of the DNA concentration to the sample weight. The DNA products extracted from feces or intestines of tadpoles(S11–S16 and S22–S27) were subsequently stored in a–20°C freezer until the high-throughput sequencing of 16S rDNA amplicons.

2.4. 16S rDNA amplicon high-throughput sequencingThe library construction of 16S rDNA amplicons and high-throughput sequencing on MiSeq (Illumina Inc.,USA) were achieved in a commercial company (Genergy Inc., China). Specifically, the hypervariable regions of V3–V4 in the bacterial 16S rDNA were amplified from the microbiota DNA products using the universal primer pair 341F–CCTACGGGNGGCWGCAG and 785R–GACTACHVGGGTATCTAATCC (Klindworthet al.,2013). The amplicons were generated by a two-step,tailed PCR on the DNA products in terms of the 16S Metagenomic Sequencing Library Preparation protocol with some modifications (http://support.illumina.com/downloads/16s_metagenomic_sequencing_library_preparation.html). The volume of each PCR solution was 25 μL, which consisted of 10 ng DNA, 2.5 μL 10 ×Takara Ex Taq Buffer, dNTP (2.5 mmol/L each), 5 pmol/L forward primer, 5 pmol/L reverse primer, 0.1 μL Takara Ex Taq, and ddH2O. The 1stPCR condition was 94°C for 3 min, 20 cycles (94°C for 10 s, 55°C for 15 s, 72°C for 30 s) and 72°C for 7 min. The 2ndPCR condition was 94°C for 3 min, 5 cycles (94°C for 10 s, 55°C for 15 s,72°C for 30 s) and 72°C for 7 min. Finally, the 16S rDNA amplicons quantified by Quant-iT? PicoGreen? dsDNA Assay Kit (Life Technologies, USA) were paired-end sequenced on the Illumina MiSeq platform.

2.5. Bioinformatic analyses of 16S rDNA ampliconsThe raw paired-end reads, after the removal of barcodes,were filtered by Trimmomatic software with three minimum thresholds (Bolgeret al., 2014), i.e., terminal base quality score equal to 25, average quality score equal to 25 in sliding windows of 50 bp with a step of 1 bp, fragment length equal to 100 bp. The paired-end sequences were merged with the flash software for a minimum overlap of 10 bp and a maximum mismatch proportion of 0.2. The merged sequences including ambiguous bases were excluded. Subsequently, the 16S rDNA reference database downloaded from NCBI was utilized for the validation of V3–V4 regions. In this study,16S rDNA amplicons were successfully sequenced in 32 samples excluding S14_S, S17_G and S20_G (Tables 1,S1).

Table 1 Sequencing results for successfully sequenced 32 samples.

Operational taxonomic units (OTUs) were assigned to the 16S rDNA amplicons using UPARSE software with the identity threshold of 97% (Edgar, 2013). The annotation of representative OTU sequences was executed through sequence alignments on the RDP database using the RDP classifier (confidence threshold = 0.8) (Wanget al., 2007).

The mothur program was used to calculate α diversity indexes (Schlosset al., 2009), i.e., richness (OTU observed, ACE and Chao1), diversity and evenness(Shannon, Simpson and Shannoneven), and Good’s coverage. The rarefaction and Shannon curves were drawn to measure whether the sequencing depth was appropriate for the richness and diversity calculation. To statistically analyze the differences in α and β diversity indexes between samples, 16 881 (i.e., minimum number of sequences in all samples) sequences in each sample were subsampled. In addition, we drew stackbars and Venn charts to show the taxonomic composition and abundance of samples. Furthermore, we used Past software (v3.14) to execute principal coordinates analysis(PCoA) based on Bray-Curtis dissimilarity calculated from taxonomic abundances to analyze the β diversity of microbiota structure (Hammeret al., 2001). The transformation exponent was set to the default value (i.e.,2).

2.6. Statistical analyses of the DNA extraction quality and the microbiota structureTo compare the DNA extraction quality of phenol-chloroform method and the TIANamp Stool DNA Kit, we utilized SPSS software(v20.0, IBM Corporation) to perform Mann-Whitney Rank Sum Test and two-way ANOVA (two factors:experimental material and DNA extraction method) in terms of A260/A280, A260/A230 and DNA yield rate.Pairwise multiple comparisons were performed in twoway ANOVA with Bonferroni correction when overall significance level was less than 0.05. In addition, fitting curves, which sketched the relationship of the DNA yield rate and sample weight, were compared between the phenol-chloroform method and the commercial stool kit.

To test the effect of experimental material and DNA extraction method on the microbiota structures, we executed two-way ANOVA to compare α diversity indexes of intestinal and fecal microbiota in S11–S33 tadpoles. We hypothesized that the factors, i.e., host genetic background and diet, could bias the comparative analysis of microbiota structure. Therefore, we used pairedt-Test or Wilcoxon Signed Rank Test to compare α diversity indexes of gut microbiota between 10 pairs of littermate tadpoles (S11–S16, S18–S21 and S24–S33),i.e., S11_G–S16_G, S18_G–S21_G versus S24_G–S33_G(Gut_ph vs Gut_kit). Subsequently, the comparisons of α diversity indexes between intestinal and fecal microbiota were performed on 12 tadpoles (S11–S16 and S22–S27),i.e., S11_S–S16_S versus S11_G–S16_G (Stool_ph vs Gut_ph) and S22_S–S27_S versus S22_G–S27_G (Stool_kit vs Gut_kit). The intestinal/fecal metagenomic DNA of S11–S16 and S22–S27 was extracted using the phenolchloroform method and TIANamp Stool DNA Kit,respectively.

In addition to permutational multivariate ANOVA(PERMANOVA), we also performed a Mantel test on the Bray-Curtis dissimilarities calculated from taxonomic(i.e., OTU, genus and phylum) abundances in the above three cases, i.e.,Gut_ph vs Gut_kit,Stool_ph vs Gut_phandStool_kit vs Gut_kit. The number of permutations was set to 9999 for PERMANOVA and Mantel tests. To test which taxonomies significantly affected the structural divergence of microbiota in these three cases, Lda Effective Size (LEfSe) were executed with recommended options (http://huttenhower.sph.harvard.edu/galaxy/).

3. Results

3.1. Efficacy comparison in terms of DNA purity and yield rateAmong the tests of the between-method difference in A260/A280, A260/A230 and DNA yield rate, a significant difference was detected in A260/A230(Mann-Whitney Rank Sum Test:U= 205,n1= 23,n2=27,P= 0.04) and DNA yield rate (U= 66,n1= 15,n2=17,P= 0.02) but not A260/A280 (U= 300.5,n1= 23,n2= 27,P= 0.85) (Table S1; Figures 1, 2). Whereas when the two-way ANOVA was applied (Table 2), a significant between-method difference was detected in A260/A230 and DNA yield rate but not in A260/A280. A significant between-material difference was detected in A260/A280 and A260/A230 but not in DNA yield rate. No significant interaction was detected between these two factors (i.e.,experimental material and DNA extraction method). In the pairwise multiple comparisons, A260/A230 (avg.± std. err.) of stool samples (1.05 ± 0.11,n= 12) was significantly different from that of gut samples (1.51 ±0.08,n= 23) and mixed stool samples (1.51 ± 0.10,n=15). Furthermore, the fitting curves for DNA yield rate and sample weight show that the DNA yield rate in the phenol-chloroform method and the commercial stool kit decreases along with saturated DNA extraction materials in solutions (Figure S1). The phenol-chloroform method has a greater DNA yield rate than the commercial stool kit does at a given sample weight.

Figure 1 Comparisons of A260/A280 and A260/A230 between TIANamp Stool DNA Kit (n = 23) and phenol-chloroform method (n = 27).*P ＜ 0.05.

Table 2 Two-way ANOVA on A260/A280, A260/A230 and DNA yield rate.

3.2 Efficacy comparisons in terms of microbiota structuresInGut_ph vs Gut_kit, total number of OTUs was 121, among which 4 belonged toGut_ph(3.3%), 4 belonged toGut_kit(3.3%) and 113 belonged to between groups (93.4%) (Figure 3). InStool_ph vs Gut_ph, total number of OTUs was 100, among which 6 belonged toStool_ph(6%), 8 belonged toGut_ph(8%) and 86 belonged to between groups (86%). InStool_kit vs Gut_kit, total number of OTUs was 114, among which 9 belonged toStool_ kit(7.9%), 7 belonged toGut_ kit(6.1%) and 98 belonged to between groups (86%). The Venn chart indicates that the majority of OTUs could be extracted by using the phenol-chloroform method or the commercial stool kit from intestines or stools.

Figure 2 Comparison of DNA yield rate between TIANamp Stool DNA Kit (n = 17) and phenol-chloroform method (n = 15). *P ＜ 0.05.

No significant differences in α diversity indexes were detected in between- material groups and between-method groups by using two-way ANOVA (Table S2). When host genetic background and diet factors were taken into account, we detected no significant differences in almost all α diversity indices inGut_ph vs Gut_kit,Stool_ph vs Gut_phandStool_kitandGut_kit(Table 3; Figure 4). One exception was that the richness values (i.e., number of observed OTUs) showed a significant difference betweenStool_kitandGut_kit(P= 0.02,n= 6) (Table 3; Figure 4).

None of the Bray-Curtis dissimilarities calculated from OTU, genus and phylum abundances showed a significant between-group difference inGut_ph vs Gut_kit(Table 4; Figures 5, S2, S3). However, only the Bray-Curtis dissimilarity based on genus abundances gave rise to a weak between-group correlation (R= 0.32,P= 0.04). In the PERMANOVA forStool_ph vs Gut_ph(Table 4; Figures 5, S2, S3), the OTU, genus and phylum abundances also showed an insignificant between-group variation in terms of Bray-Curtis dissimilarity. In addition,both Bray-Curtis dissimilarities calculated from OTU and genus abundances showed a significant between-group correlation, i.e.,R= 0.64 (P= 0.04) andR= 0.58 (P=0.02). As for theStool_kit vs Gut_kit(Table 4; Figures 5,S2, S3), both Bray-Curtis dissimilarities calculated from OTU and genus abundances also showed a significant between-group correlation, i.e.,R= 0.65 (P= 0.01) andR= 0.64 (P= 0.01). Nevertheless, significant betweengroup variations was detected in both Bray-Curtis dissimilarities based on OTU and genus abundances(PERMANOVA:F1,10= 4.33,P= 0.02) andF1,10= 4.37,P= 0.03).

Table 3 Alpha diversity indexes (avg. ± std. dev.) calculated from OTU tables rarefied to 16 881 sequences per sample in three cases, i.e.,Gut_ph vs Gut_kit (n = 9), Stool_ph vs Gut_ph (n = 5) and Stool_kit vs Gut_kit (n = 6).

Figure 3 The distribution of OTUs in fecal and intestinal microbiota based on two DNA extraction methods.

Table 4 PERMANOVA and Mantel test on the Bray-Curtis dissimilarities calculated from taxonomic (i.e., OTU, genus and phylum)abundances in three cases, i.e., Gut_ph vs Gut_kit (n = 9), Stool_ph vs Gut_ph (n = 5) and Stool_kit vs Gut_kit (n = 6).

Figure 4 Comparisons of α diversity indices in three cases, i.e.,Gut_ph vs Gut_kit, Stool_ph vs Gut_ph and Stool_kit vs Gut_kit. The bars and error bars represent mean values and standard deviations, respectively. *P ＜ 0.05.

InGut_ph vs Gut_kit, we identified by the LEfSe analysis that three genera of phylum Proteobacteria showed significant between-group divergence, i.e.,Sphingorhabdus biased toGut_ph, and Coxiella biased toGut_kitalong with an unnamed genus in order Rhizobiales (Figure 6). As forStool_ph vs Gut_ph, five genera (i.e., Hydrogenophaga, Rhizobium,Brevundimonas, an unnamed genus in order Rhizobiales,an unnamed genus in family Sphingomonadaceae) in phylum Proteobacteria biased toStool_ph, and genus Clostridium_XlVa in phylum Firmicutes biased toGut_ph(Figure 6). In addition, significant betweengroup divergences of four phyla (i.e., Proteobacteria,Bacteroidetes, Verrucomicrobia, Firmicutes) composed of 15 genera were detected inStool_kit vs Gut_kit, i.e.,12 genera in the phyla of Proteobacteria, Bacteroidetes and Verrucomicrobia biased toStool_kit, and 3 genera in the phyla of Proteobacteria and Firmicutes biased toGut_kit(Figure S4). The shared between-material bias of genera inStool_ph vs Gut_phandStool_kit vs Gut_kitincluded four genera in phylum Proteobacteria biased toStool(i.e., Hydrogenophaga, Rhizobium, Brevundimonas,an unnamed genus in family Sphingomonadaceae) and one genus in phylum Firmicutes biased toGut(i.e.,Clostridium_XlVa).

4. Discussion

Figure 5 Principal coordinates analysis (PCoA) based on Bray-Curtis dissimilarity calculated from OTU abundances in three cases, i.e.,Gut_ph vs Gut_kit, Stool_ph vs Gut_ph and Stool_kit vs Gut_kit..

Figure 6 LEfSe results in the case of Gut_ph vs Gut_kit and Stool_ph vs Gut_ph.

The phenol-chloroform method is a classic and costeffective approach to extract eukaryotic or prokaryotic genomic DNA (Ausubel, 2002). The DNA yield rate of the phenol-chloroform method is significantly higher than those of several commercial methods (Wagner Mackenzieet al., 2015). Similarly, this study demonstrated the DNA extraction efficiency of the phenol-chloroform method outweighed that of TIANamp Stool DNA Kit in terms of DNA yield rate. In addition, the phenol-chloroform method produced a better A260/A230 than TIANamp Stool DNA Kit did, which probably resulted from more effective DNA washing in the last step of phenolchloroform method. Due to the poor operationality and technical overlook of mucosal bacteria in separating gut contents from intestines, we extracted the mixture of gut microbiota DNA and host genomic DNA from intestinal samples. Nevertheless, we applied the universal primer pairs to specifically amplify V3–V4 regions of bacterial 16S rDNA (Klindworthet al., 2013). It has been reported that the mechanical treatments (e.g., bead beating) of samples can produce more efficient cell lysis of Grampositive bacteria and yield high amounts of DNA (Ferrandet al., 2014; Guo and Zhang, 2013). However, no method has become a gold standard for 16S rDNA highthroughput sequencing (Yamagishiet al., 2016). Methods with a mechanical treatment face a tradeoff between cell lysis efficiency and DNA disruption level. Due to the easy digestion of tadpole intestines and feces we did not take the mechanical treatments into account.

The composition and diversity of gut microbiota are undoubtedly affected by multiple factors, e.g., host genetic background, dietary profile and environmental situation(Dabrowska and Witkiewicz, 2016; Davenport, 2016; Jinet al., 2017; Voreadeset al., 2014). The annotation on the gut microbiota structure can be susceptibly biased by the experimental protocols (Chooet al., 2015; Claassenet al., 2013; Yuanet al., 2012). However, many studies have demonstrated that biological variations outweigh technical variations generated by DNA extraction methods (Salonenet al., 2010; Wagner Mackenzieet al., 2015; Wesolowska-Andersenet al., 2014) and sample storage conditions (Blekhmanet al., 2016; Fouhyet al., 2015). Here we homogenized and minimized biological variations in between groups as far as possible.Nevertheless, it is impossible to eliminate each interfering factor. For instance, the littermate tadpole pairs used inGut_ph vs Gut_kitpossess a similar but not an identical genetic background, thereby the between-method heterogeneity is possibly enhanced. From the analyses ofGut_ph vs Gut_kit, the technical variation generated by DNA extraction methods was outweighed by the inter-subject variation. However, the phenol-chloroform method and the commercial stool kit probably resulted in a significant inconsistency in the structural composition of microbiota, e.g., OTU and phylum abundances. To ensure biological differences outweigh systematic biases,we had better use the identical standardized protocols for the comparative analysis of gut microbial consortia.

Although the microbiota structure shows a significant variation between feces and intestinal contents in the case of the sophisticated intestines (Guet al., 2013), feces has been applied as an effective noninvasive material for the study of gut microbiota in mammals. Larsenet al.(2015) used ribosomal intergenic spacer analysis to reveal that the microbiota structure of fishes was similar but significant different in feces and intestines. However, their preparation procedure of fecal and intestinal samples,i.e., the fecal samples were squeezed from the intestinal samples, probably enhance the between-material variations. In this study, the fecal microbiota based on the phenol-chloroform method can more efficiently reflect the gut microbiota in terms of composition and diversity. On the contrary, the feces and intestines possess more inconsistent microbiota structures deduced from the commercial stool kit. When we applied the PCoA to compare the samples fed with different food types, the experimental material factor seemed to be dominant rather than genetic background and food type factors in the case of commercial stool kit. We will explore and discuss the effects of multi-factors (e.g., genetic background and food type) on the microbiota structure of Asiatic toad tadpoles in a further study. Even though no significant difference between fecal and intestinal microbiota was detected in the phenol-chloroform method, betweenmaterial variations do exist, e.g., four genera in phylum Proteobacteria were more abundant in feces and one genus in phylum Firmicutes in intestines. We argue that the inconsistency possibly resulted from the moderate interference of microorganisms in water and tadpole skin to feces.

5. Conclusion

According to the DNA extraction quality and structural comparisons between fecal and intestinal microbiota,the phenol-chloroform method is probably more robust than a commercial fecal reagent kit (TIANGEN Biotech Co., Ltd.) in evaluating the gut microbiota structure of amphibian tadpoles with feces. To the best of our knowledge, this study provides the first evidence that feces of amphibian tadpoles can be applied as an effective noninvasive material for the study of gut microbiota.

AcknowledgementsWe thank Weizhao YANG(Lund University), Yunfei WANG (MD Anderson Cancer Center), Bruce WALDMAN (Seoul National University) and anonymous referees for their comments and suggestions. This study was supported by the National Natural Science Foundation of China (NSFC 31600104), Key Scientific Research Project of Higher Education in Henan Province (No. 17B180004), National Undergraduate Training Program for Innovation and Entrepreneurship (No. 201610477013), Ph.D. Research Startup Foundation of Xinyang Normal University (No.0201424), and Nanhu Scholars Program for Young Scholars of Xinyang Normal University.

Arumugam M., Raes J., Pelletier E., Le Paslier D., Yamada T.,Mende D. R., Fernandes G.R., Tap J., Bruls T., Batto J.M.,Bertalan M., Borruel N., Casellas F., Fernandez L., Gautier L., Hansen T., Hattori M., Hayashi T., Kleerebezem M.,Kurokawa K., Leclerc M., Levenez F., Manichanh C., Nielsen H. B., Nielsen T., Pons N., Poulain J., Qin J., Sicheritz-Ponten T., Tims S., Torrents D., Ugarte E., Zoetendal E. G.,Wang J., Guarner F., Pedersen O., de Vos W. M., Brunak S.,Dore J., Antolin M., Artiguenave F., Blottiere H. M., Almeida M., Brechot C., Cara C., Chervaux C., Cultrone A., Delorme C., Denariaz G., Dervyn R., Foerstner K. U., Friss C., van de Guchte M., Guedon E., Haimet F., Huber W., van Hylckama-Vlieg J., Jamet A., Juste C., Kaci G., Knol J., Lakhdari O.,Layec S., Le Roux K., Maguin E., Merieux A., Melo Minardi R., M’Rini C., Muller J., Oozeer R., Parkhill J., Renault P.,Rescigno M., Sanchez N., Sunagawa S., Torrejon A., Turner K., Vandemeulebrouck G., Varela E., Winogradsky Y., Zeller G., Weissenbach J., Ehrlich S. D., Bork P.2011. Enterotypes of the human gut microbiome. Nature, 473: 174?180

Ausubel F. M.2002. Short protocols in molecular biology: Wiley

Bahl M. I., Bergstrom A., Licht T. R.2012. Freezing fecal samples prior to DNA extraction affects the Firmicutes to Bacteroidetes ratio determined by downstream quantitative PCR analysis.FEMS Microbiol Lett, 329: 193?197

Blekhman R., Tang K., Archie E. A., Barreiro L. B., Johnson Z. P., Wilson M. E., Kohn J., Yuan M. L., Gesquiere L.,Grieneisen L. E., Tung J.2016. Common methods for fecal sample storage in field studies yield consistent signatures of individual identity in microbiome sequencing data. Sci Rep, 6:31519

Bletz M. C., Goedbloed D. J., Sanchez E., Reinhardt T., Tebbe C. C.2016. Amphibian gut microbiota shifts differentially in community structure but converges on habitat-specific predicted functions. Nat Commun, 7: 13699

Bolger A.M., Lohse M., Usadel B.2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 30:2114?2120

B

rooks A. W., Kohl K. D., Brucker R. M., van Opstal E.J., Bordenstein S. R.2016. Phylosymbiosis: relationships and functional effects of microbial communities across host evolutionary history. PLoS Biol, 14: e2000225

Choo J. M., Leong L. E., Rogers G. B.2015. Sample storage conditions significantly influence faecal microbiome profiles. Sci Rep, 5: 16350

Claassen S., du Toit E., Kaba M., Moodley C., Zar H. J., Nicol M. P.2013. A comparison of the efficiency of five different commercial DNA extraction kits for extraction of DNA from faecal samples. J Microbiol Methods, 94: 103?110

Colombo B. M., Scalvenzi T., Benlamara S., Pollet N.2015.Microbiota and mucosal immunity in amphibians. Front Immunol, 6: 111

Dabrowska K., Witkiewicz W.2016. Correlations of Host Genetics and Gut Microbiome Composition. Front Microbiol, 7: 1357

Davenport E. R.2016. Elucidating the role of the host genome in shaping microbiome composition. Gut Microbes, 7: 178?184

Dethlefsen L., Huse S., Sogin M. L., Relman D. A.2008. The pervasive effects of an antibiotic on the human gut microbiota,as revealed by deep 16S rRNA sequencing. PLoS Biol, 6: e280

Eckburg P. B., Bik E. M., Bernstein C. N., Purdom E.,Dethlefsen L., Sargent M., Gill S. R., Nelson K. E., Relman D. A.2005. Diversity of the human intestinal microbial flora.Science, 308: 1635?1638

Edgar R. C.2013. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods, 10: 996?998

Ferrand J., Patron K., Legrand-Frossi C., Frippiat J. P., Merlin C., Alauzet C., Lozniewski A.2014. Comparison of seven methods for extraction of bacterial DNA from fecal and cecal samples of mice. J Microbiol Methods, 105: 180?185

Fouhy F., Deane J., Rea M. C., O’Sullivan O., Ross R. P.,O’Callaghan G., Plant B. J., Stanton C.2015. The effects of freezing on faecal microbiota as determined using MiSeq sequencing and culture-based investigations. PLoS One, 10:e0119355

Gu S., Chen D., Zhang J.N., Lv X., Wang K., Duan L. P., Nie Y., Wu X. L.2013. Bacterial community mapping of the mouse gastrointestinal tract. PLoS One, 8: e74957

Guo F., Zhang T.2013. Biases during DNA extraction of activated sludge samples revealed by high throughput sequencing. Appl Microbiol Biotechnol, 97: 4607?4616

Hammer ?., Harper D. A., Ryan P. D.2001. PAST:Paleontological Statistics Software Package for Education and Data Analysis. Palaeontologia Electronica, 4: 1?9

Hof C., Araujo M. B., Jetz W., Rahbek C.2011. Additive threats from pathogens, climate and land-use change for global amphibian diversity. Nature, 480: 516?519

Jiménez R. R., Sommer S.2017. The amphibian microbiome:natural range of variation, pathogenic dysbiosis, and role in conservation. Biodivers Conserv, 26: 763?786

Jin Y., Wu S., Zeng Z., Fu Z.2017. Effects of environmental pollutants on gut microbiota. Environ Pollut, 222: 1?9

Klindworth A., Pruesse E., Schweer T., Peplies J., Quast C.,Horn M., Glockner F.O.2013. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and nextgeneration sequencing-based diversity studies. Nucleic Acids Res, 41: e1

Kohl K. D., Cary T. L., Karasov W. H., Dearing M. D.2013.Restructuring of the amphibian gut microbiota through metamorphosis. Environ Microbiol Rep, 5: 899?903

Kohl K. D., Amaya J., Passement C. A., Dearing M. D., McCue M. D.2014. Unique and shared responses of the gut microbiota to prolonged fasting: a comparative study across five classes of vertebrate hosts. FEMS Microbiol Ecol, 90: 883?894

Larsen A. M., Mohammed H. H., Arias C. R.2015. Comparison of DNA extraction protocols for the analysis of gut microbiota in fishes. FEMS Microbiol Lett, 362

Ley R. E., Lozupone C.A., Hamady M., Knight R., Gordon J.I.2008. Worlds within worlds: evolution of the vertebrate gut microbiota. Nat Rev Microbiol, 6: 776?788

Li H., Qu J., Li T., Yao M., Li J., Li X.2017a. Gut microbiota may predict host divergence time during Glires evolution. FEMS Microbiol Ecol, 93: fix009-fix009

Li T., Long M., Li H., Gatesoupe F. J., Zhang X., Zhang Q.,Feng D., Li A.2017b. Multi-omics analysis reveals a correlation between the host phylogeny, gut microbiota and metabolite profiles in cyprinid fishes. Front Microbiol, 8

Lozupone C. A., Stombaugh J., Gonzalez A., Ackermann G.,Wendel D., Vazquez-Baeza Y., Jansson J. K., Gordon J.I., Knight R.2013. Meta-analyses of studies of the human microbiota. Genome Res, 23: 1704?1714

Mashoof S., Goodroe A., Du C. C., Eubanks J. O., Jacobs N.,Steiner J. M., Tizard I., Suchodolski J. S., Criscitiello M. F.2013. Ancient T-independence of mucosal IgX/A: gut microbiota unaffected by larval thymectomy in Xenopus laevis. Mucosal Immunol, 6: 358?368

Maukonen J., Simoes C., Saarela M.2012. The currently used commercial DNA-extraction methods give different results of clostridial and actinobacterial populations derived from human fecal samples. FEMS Microbiol Ecol, 79: 697?708

McFall-Ngai M., Hadfield M. G., Bosch T. C., Carey H. V.,Domazet-Loso T., Douglas A. E., Dubilier N., Eberl G.,Fukami T., Gilbert S. F., Hentschel U., King N., Kjelleberg S., Knoll A. H., Kremer N., Mazmanian S. K., Metcalf J.L., Nealson K., Pierce N. E., Rawls J. F., Reid A., Ruby E.G., Rumpho M., Sanders J. G., Tautz D., Wernegreen J. J.2013. Animals in a bacterial world, a new imperative for the life sciences. Proc Natl Acad Sci U S A, 110: 3229?3236

Nava G. M., Friedrichsen H. J., Stappenbeck T. S.2011. Spatial organization of intestinal microbiota in the mouse ascending colon. ISME J, 5: 627?638

Nicholson J. K., Holmes E., Kinross J., Burcelin R., Gibson G., Jia W., Pettersson S.2012. Host-gut microbiota metabolic interactions. Science, 336: 1262?1267

Salonen A., Nikkila J., Jalanka-Tuovinen J., Immonen O.,Rajilic-Stojanovic M., Kekkonen R. A., Palva A., de Vos W.M.2010. Comparative analysis of fecal DNA extraction methods with phylogenetic microarray: effective recovery of bacterial and archaeal DNA using mechanical cell lysis. J Microbiol Methods,81: 127?134

Schloss P. D., Westcott S. L., Ryabin T., Hall J. R., Hartmann M., Hollister E. B., Lesniewski R. A., Oakley B. B., Parks D. H., Robinson C. J., Sahl J. W., Stres B., Thallinger G. G.,Van Horn D. J., Weber C. F.2009. Introducing mothur: opensource, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol, 75: 7537?7541

Shapira M.2016. Gut microbiotas and host evolution: scaling up symbiosis. Trends Ecol Evol, 31: 539?549

Song S.J., Amir A., Metcalf J. L., Amato K. R., Xu Z. Z.,Humphrey G., Knight R.2016. Preservation methods differ in fecal microbiome stability, affecting suitability for field studies.mSystems, 1

Su C., Lei L., Duan Y., Zhang K.Q., Yang J.2012. Cultureindependent methods for studying environmental microorganisms: methods, application, and perspective. Appl Microbiol Biotechnol, 93: 993?1003

Vences M., Lyra M. L., Kueneman J. G., Bletz M. C., Archer H. M., Canitz J., Handreck S., Randrianiaina R. D., Struck U., Bhuju S., Jarek M., Geffers R., McKenzie V. J., Tebbe C. C., Haddad C. F., Glos J.2016. Gut bacterial communities across tadpole ecomorphs in two diverse tropical anuran faunas.Naturwissenschaften, 103: 25

Voreades N., Kozil A., Weir T. L.2014. Diet and the development of the human intestinal microbiome. Front Microbiol, 5: 494

Wagner Mackenzie B., Waite D. W., Taylor M. W.2015.Evaluating variation in human gut microbiota profiles due to DNA extraction method and inter-subject differences. Front Microbiol, 6: 130

Wang Q., Garrity G. M., Tiedje J. M., Cole J. R.2007. Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol, 73: 5261-5267

Weng F. C., Yang Y. J., Wang D.2016. Functional analysis for gut microbes of the brown tree frog (Polypedates megacephalus) in artificial hibernation. BMC Genomics, 17: 1024

Weng F. C. H., Shaw G. T. W., Weng C. Y., Yang Y. J., Wang D.2017. Inferring microbial interactions in the gut of the Hong Kong whipping frog (Polypedates megacephalus) and a validation using probiotics. Front Microbiol, 8

Wesolowska-Andersen A., Bahl M. I., Carvalho V., Kristiansen K., Sicheritz-Ponten T., Gupta R., Licht T. R.2014. Choice of bacterial DNA extraction method from fecal material influences community structure as evaluated by metagenomic analysis.Microbiome, 2: 19

Yamagishi J., Sato Y., Shinozaki N., Ye B., Tsuboi A., Nagasaki M., Yamashita R.2016. Comparison of boiling and robotics automation method in DNA extraction for metagenomic sequencing of human oral microbes. PLoS One, 11: e0154389

Yatsunenko T., Rey F. E., Manary M. J., Trehan I., Dominguez-Bello M. G., Contreras M., Magris M., Hidalgo G.,Baldassano R. N., Anokhin A. P., Heath A. C., Warner B.,Reeder J., Kuczynski J., Caporaso J. G., Lozupone C. A.,Lauber C., Clemente J. C., Knights D., Knight R., Gordon J. I.2012. Human gut microbiome viewed across age and geography. Nature, 486: 222?227

Yuan S., Cohen D. B., Ravel J., Abdo Z., Forney L. J.2012.Evaluation of methods for the extraction and purification of DNA from the human microbiome. PLoS One, 7: e33865