Diet Rather than Genetic Status Shapes the Gut Microbiota in Two Congeneric Snakes

Background: Various external and internal factors affect the gut microbiota of animals. The colonization and proliferation of gut microbes have been studied in a diverse array of animal taxa but remain poorly known in snakes. Here, we used the 16S rRNA gene sequencing technology on the Roach 454 platform to analyze the gut microbiota composition using fecal samples collected from three snake groups [gravid females, newly hatched (preprandial) hatchlings and postprandial hatchlings] of two congeneric colubrid snake species (Elaphe carinata and E. taeniura) that are sympatric across a wide range in mainland China. We tested two hypotheses. First, the gut microbiota should not differ between the two species at hatching if the maternal or genetic contribution has no role in affecting post-hatching gut microbial colonization. Second, differences in the gut microbiota between newly hatched (preprandial) and postprandial hatchlings should not exist in both species if the dietary contribution has no role in affecting post-hatching gut microbial colonization. Results: The top three dominant phyla were Firmicutes, Bacteroidetes, and Proteobacteria in both species. None of the measured alpha diversity indexes differed among the three snake groups or between the two species. The relative abundance of the gut microbiota differed among the three snake groups and between the two species, and so did the relative abundances of the functions associated with the metabolism, cellular processes and environmental information processing. Evidence from gravid females and hatchlings showed that the gut microbiota composition was similar between the two species. The metabolism held the overwhelming predominance of functional categories at the top level in both species. Conclusion: Only the relative abundance of the gut microbiota differed between the two species, and the gut microbiota composition changed rapidly in postprandial hatchlings and differed among the three snakes groups in both species. From these ndings, we may conclude that the dietary rather than the maternal or genetic contribution affects gut microbial colonization in snakes.


Background
The past few years have witnessed a great advance in high-throughput sequencing technologies and novel analytical methods that have been widely used to interpret gut microbial structure and abundance and to explore the symbiotic relationship between gut microbes and their host animals. The gut microbiota affects many physiological functions of the host, including growth [1], behavior [2], metabolism [3], reproduction [4] and in ammatory/immune responses [5]. The host can enhance its tness by providing and maintaining an in vivo environment suitable for gut microbes [6]. Gut microbial dysbiosis may trigger chronic in ammatory diseases and certain types of cancer [7][8][9].
Geophagy (soil-eating) and coprophagy (feces-eating) both are important ways through which animals acquire gut microbes. Soil is a dynamic reservoir of microorganisms, playing a key role in maintaining and regulating the gut microbiota in animals. Geophagy explains why gut microbes are similar to the soil microorganisms in the giant armadillo Priodontes maximus [24,25]. Another example for the microbe colonization through geophagy is that newly hatched iguanas (Iguana iguana) ingest soil within the nest before they digest food, and they even harvest gut microbes by ingesting feces of their seniors after the second or third week of their life [26]. Animals ingest soil in manners that make them better able to digest certain substances such as plant material [26] and chitin [25].
Numerous studies have suggested that complex mechanisms and factors are involved in the colonization and maintenance of gut microbes, and that the gut microbiota changes rapidly soon after the external environment changes. The gastrointestinal tracts in vertebrates are considered to be essentially sterile at the time of birth, and their microbiomes colonize and develop from very early days of life [6, 27,28]. The original gut microbiota plays a vital role in the human gastrointestinal ecosystem and has a long-term effect on its dynamics [27]. The diversity and composition of gut microbes increase with age in ostriches (Struthio camelus), and this change is accompanied by the multiple colonization and extinction events of microbial species [6]. Evidence from juvenile lizards shows maternally derived gut microorganisms, as revealed by the fact that they share ~35% of their gut microbes with their mothers [16]. Clearly, the original gut colonizers can be produced and rapidly proliferate according to the host pro les [6,27]. However, further studies on more species are needed to better understand the colonization and maintenance mechanisms of the gut microbiota in the host.
Ontogenetic shifts in the structure and composition of the gut microbiota have been documented in a diverse array of animal taxa [6, 16, [21][22][23]. For example, the gut microbial diversity changes rapidly in human newborns, taking only six days to change Lactococcus lactis in the gut into Enterobacter coli [29]. Another example is that, relative to other gut microbes, the abundances of bacteria of the classes Verrucomicrobiae and Erysipelotrichia are high in one-week-old ostriches but reduce successively when they grow to two-week-old [6]. Changes in gut microbes are closely related to the environmental and dietary conditions, and feeding mode can shape the acquisition and structure of the original gut microbiota [30]. Taken together, previous studies have shed some light on the ecology of gut microbes and their symbiotic relationships with the host, yet more relevant experiments should be done to test if there are some patterns generalisable to animal hosts.
Here, we describe a study using 16S rRNA gene sequencing on the Roach 454 platform to examine whether the composition and abundance of gut microorganisms differ between the king ratsnake Elaphe carinata and the stripe-tailed ratsnake Elaphe taeniura and display ontogenetic shifts. Speci cally, we used fecal samples collected from gravid females, newly hatched (preprandial) hatchlings and postprandial hatchlings of each species to test two hypotheses. First, the gut microbiota should not differ between the two species at hatching if the maternal or genetic contribution has no role in affecting posthatching gut microbial colonization. Second, differences in the gut microbiota between newly hatched (preprandial) and postprandial hatchlings should not exist in both species if the dietary contribution has no role in affecting post-hatching gut microbial colonization. We used these two congeneric colubrid snakes as the model animals for two reasons. First, they are sympatric across a wide range in mainland China [37,38] and therefore ideally suited to the studies examining interspeci c differences in the gut microbiota with minimized bias from phylogenetic and/or ecological differences. Second, both E. carinata [39][40][41][42] and E. taeniura [43][44][45][46] have been well studied, and individuals of different ages can be easily maintained in the laboratory for data collection.
Egg incubation and collection of hatchling fecal sample Freshly laid eggs were always collected less than 12 h post-laying and then incubated at 27 ℃ in a substrate of moist vermiculite at a water potential of −12 kPa [42]. Hatchlings were collected, weighed, and measured (for SVL and tail length) less that 6 h post-hatching.
Seven hatchling E. carinata (347-386 mm SVL, and 15.9-22.3 g) and seven hatchling E. taeniura (334-377 mm SVL, and 21.1-28.8 g) were randomly chosen to collect fecal samples. These hatchlings were individually housed in 160 × 120 × 100 mm ethanol-disinfected plastic containers until the end of the experiment. The containers contained the ethanol-disinfected ceramic refuges, and water was provided every day. Of the 14 hatchlings, three of each species were used to collect fecal samples at hatching (numbered ECH1 to ECH3 for E. carinata, and ETH1 to ETH3 for E. taeniura). The remaining four of each species were fed ad libitum with newborn mice 3 wk post-hatching, and their fecal samples were individually collected two days later (numbered ECP1 to ECP4 for E. carinata, and ETP1 to ETP4 for E. taeniura). Fecal samples from hatchlings were also stored at −20 ℃ until later extraction of genomic DNA.
We ended the experiment in late September 2013 and released all snakes (females and their hatchlings) to the sites where gravid females were collected. Our experimental procedures were approved by the Institutional Animal Care and Use Committee of Nanjing Normal University (IACUC-200422).

DNA extraction, PCR ampli cation and sequencing of fecal samples
The fecal genomic DNA was extracted using the QIAamp® PowerFecal® Pro DNA Kit (QIAGEN, Germany) according to the manufacturer protocols. The concentration and purity of the extracted DNA were measured using the ultraviolet spectrophotometer. The integrity of DNA samples was detected by 0.8% agarose gel electrophoresis at 120V for 20 min.
Polymerase chain reaction (PCR) was performed in the 25 µl reaction system consisting of 10 µM of each primer, 5 mM dNTP, 10 ⋅ PCR reaction buffer, 0.625 U Pyrobest DNA polymerase (Takara) and 2 µL (40 ng) genomic DNA. The PCR thermal cycling conditions were as follows: initial denaturation at 94 °C for 4 min, followed by 22 cycles of denaturation at 94 °C for 30 s, at 55 °C for 30 s and at 72 °C for 45 s, and a nal extension at 72 °C for 8 min. Afterwards, the PCR products were visualized using the 1.0% agarose gel and cleaned up by the Agencourt AMPure xp Beads (Beckman Coulter, USA) according to the manufacturer's instructions. Additionally, the amplicon DNA library was quanti ed using the Quant-iT PicoGreen dsDNA reagent and kit (Life Tech, Carlsbad, USA) in accordance with manufacturer's protocols.
Pyrosequencing was then performed on a 454 Life Sciences Genome Sequencer FLX instrument (Roche) (Personalbio, Shanghai, China) according to titanium chemistry following the manufacturer's instructions.
Control standards were used to optimize the quality of pyrosequencing data. Only the raw reads that were almost perfectly matched with the primers were used, and the barcode sequences with quality average > Q25 were retained. Typically, the valid reads were selected in line with the following criteria: 1. length of at least 260 bp nucleotides (barcodes and primers excluded); 2. no undetermined bases; and 3. no over six consecutively identical bases. A total of 284,056 high-quality reads were obtained. These sequences were then submitted to the National Center for Biotechnology Information (NCBI) Bioproject database (SRA accession number PRJNA642086).

Data analysis
OTU cluster and data standardization Using Usearch 7.0, we clustered the available sequences into the operational taxonomic units (OTUs) at the similarity threshold of 97% [69]. We then used RDP Classi er 2.2 [70] to analyze the representative sequence for each OUT against the Silva 132 Database at the con dence threshold of 70%. To avoid large partial sample deviations, we retained OTUs with the number of OTUs greater than 5 in at least six samples for further analysis.
The sequencing depth index was calculated by mothur 1.30.2 [71] and visualized using R 3.6 [72], thereby evaluating the adequacy of sequencing data. Afterwards, the OTUs abundance information was standardized using the sample with the least sequence number for further analysis. Later, the species accumulation box plot was used to assess whether the microbiological diversity in the 23 fecal samples re ected the diversity of the whole gut microorganisms in snakes.

Alpha and beta diversity
To analyze alpha diversity, mothur 1.30.2 was rst used to calculate the community richness (Chao1 index), community diversity (Shannon diversity index) and community evenness (Shannon even index) for each fecal sample. We then used the R software platform to visualize these indexes. We used ANOVA to examine the differences in alpha diversity among snake groups. Prior to ANOVA, data were tested for homogeneity of variances using Bartlett test and for normality using Shapiro test.
To measure beta diversity, we performed principal co-ordinates analysis (PCoA) and and the partial least squares discriminant analysis (PLS-DA) to compare the differences in structure (the relative abundance of OTUs) of gut microbiota among different groups based on OTUs. We further conducted permutational multivariate analysis of variance (PERMANOVA) with Adonis based on Bray_Curtis distance matrices with 999 permutations to analyze the bacterial similarity between different groups.
We used the linear discriminant analysis effect size (LEfSe) [73] to compare the microbial abundances from phylum to family levels among different groups, thus determining whether the gut microbiota differed between species and among snake groups. Hereafter, we used a series of pairwise tests to investigate the biological consistency among diverse subgroups. We used linear discriminatory analysis (LDA) to evaluate the effect size for each selected classi cation. In this study, only the bacterial taxa with a log LDA score > 4 (over 4 orders of magnitude) were used. We used Kruskal-Wallis test to perform multiple-sample analyses at different taxonomic levels based on the results of LEfSe.
PICRUSt was applied in searching protein sequences and predicting gene functions based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [74]. These genes were then classi ed and allocated to the corresponding KEGG pathways [75]. The numbers of functional genes in each pathway were counted to assess their relative abundances in different groups. Additionally, the LEfSe was utilized to compare the abundances of KEGG gene functions from KOs gene level 1 to level 3 among different groups, so as to determine the differences in functional genes. In addition, LDA analysis was performed to assess the effect size for each selected level. In this study, only gene classi cation with a log LDA score > 3 was used. All values were presented as mean ± SE, and all statistical analyses were conducted at the signi cance level of α = 0.05.
Differences in the gut microbiota among snake groups Within each species the three snake groups differed in the gut microbial structure and similarity (Adonis: r 2 = 0.37, F 5,17 = 2.03, p = 0.002; Fig. 2) but not in any measured diversity index (all p > 0.05; Fig.S4).

Discussion
Gut microorganisms affect many physiological functions of their hosts, and their roles were determined by the community composition and relative abundance that are related to genetic pro les and microhabitats of the host [10,11,47]. Gut microbiota is mainly composed of bacteria of the phyla Firmicutes and Bacteroidetes in amphibians [48], reptiles [49] and mammals [11,50], and bacteria of the phyla Proteobacteria and Firmicutes in insects [47], sh [10,51] and birds [52]. Bacteria of the phylum Bacteroidetes are of importance in the fermentation and degradation of sugars and plant materials [53], bacteria of the phylum Firmicutes are related to the production of diverse digestive enzymes and the potential of vitamin B biosynthesis [54,55], and bacteria of the phylum Proteobacteria are associated with cellulose activity and degradation of various aromatic compounds [56].
In this study, the top three dominant phyla were Firmicutes, Bacteroidetes, and Proteobacteria, which accounted for 90.3% of all phyla identi ed in E. carinata and 88.1% of all phyla identi ed in E. taeniura (Fig. 1A). Similar results have been reported for the cottonmouth Agkistrodon piscivorus [57], the timber rattlesnake Crotalus horridus [58], the red-necked keelback snake Rhabdophis subminiatus [34] and four species of farm-raised snakes (the ve-paced pit viper Deinagkistrodon acutus, E. carinata, the Chinese cobra Naja atra, and the common ratsnake Ptyas mucosa; [59]). In both species studied herein, the relative abundance of gut microbes exceeded 3% in seven families and seven genera ( Fig. 1B and C). The top three dominant families were Bacteroidaceae, Enterobacteriaceae and Fusobacteriaceae in E. taeniura (the total relative abundance > 61%), and Bacteroidaceae, Enterobacteriaceae and Lachnospiraceae in E. carinata (the total relative abundance > 59%). Bacteria of the genus Bacteroides were most abundant in both species, accounting for 33% in E. carinata and 24% in E. taeniura (Fig. 1C).
Interspeci c differences in the gut microbiota are caused by multiple factors [16,32,49,59]. Genetic and dietary correlates of the gut microbiota have been detected in farm-raised snakes feeding on the same food [59], captive bank voles fed with different items of food [16], and wild toad-headed lizards (Phrynocephalus vlangalii) from different altitudes [13]. In this study, we found that the gut microbiota was similar between the two species but differed among gravid females, newly hatched hatchlings and postprandial hatchlings in both species (Fig. 2). These ndings con rm that the dietary rather than the genetic or maternal contribution has a more important role in shaping the gut microbiota.
The ratio of Firmicutes to Bacteroidetes (F/B ratio) is related to energy consumption and body-mass change and is affected by multiple factors [13,14,50,60,61]. Host animals with a greater F/B ratio are better able to acquire energy rapidly. For instance, Rhesus macaques (Macaca mulatta) from Tibet have a greater F/B ratio than those from low-altitude (and thus more temperate and less anoxic) localities as an adaptation to cold and anoxic environments in this region [14]. In this study, data from hatchlings showed that the ratio was greater in E. taeniura (0.88) than in E. carinata (0.62), suggesting an enhanced ability of hatchlings to acquire energy in the former species. Interestingly, functional genes related to metabolism were also more abundant in the gut microbiota of hatchling E. taeniura (Fig. 4), suggesting that the difference in the F/B ratio between the two species is related to the functional genes involved in digestion and absorption of energy.
The gut microbial pro les at early stages of life are complicated by the natural dynamics and complex interactions between intrinsic and extrinsic factors [62]. Factors such as early growth, household location and antibiotic experience during pregnancy determine the early gut microbial composition in humans [62]. Similar results have also been reported for lizards [26] and birds [6]. For instance, the gut microbiota of hatchling Iguana ihuana colonizes and proliferates through ingesting the surrounding soil and feces of their relatives [26]. In newborn lambs, feeding modes have been found to shape the acquisition and structure of the initial gut microbiota [30].
It has been widely recognized that diet affects the gut microbiota at early stages of life. However, none of the measured diversity indexes differed among the three snake groups and between the two species. Instead, the relative abundances of microbial members across different taxonomic levels showed signi cant differences (Fig. 3). In both species, bacteria of the phylum Fusobacteria had a higher relative abundance in gravid females. bacteria of the phylum Fusobacteria are usually found in the gut of infected animals or scavengers such as the American alligator Alligator mississippiensis [35], The intraspeci c differences in bacterial abundance observed in gravid females, newly hatched hatchling and postprandial hatchlings suggest that the feeding process may result in changes in the relative abundances of gut microorganisms in both species (Figs. 2 and 3). Similar results have also been reported for humans where signi cant changes in the membership and relative abundance of the gut microbiota occur at their early stages of life [29]. In ostriches, bacteria of the family Bacteroidaceae promote juvenile growth, whereas bacteria of the families Enterobacteriaceae, Enterococcaceae and Lactobacillaceae retard juvenile growth [6]. In this study, the relative abundance of the gut microbiota differed signi cantly between the two snake species (Fig. 3). Postprandial hatchlings differed from newly hatched hatchlings and gravid snakes in the bacterial abundance within each species. Postprandial hatchlings showed similar variation trends in both species, indicating that these two closely related snake species display a similar pattern of the gut microbial colonization and proliferation.
The gut microbial transmission between individuals is a common phenomenon [16, 53,66]. The composition and structure of the gut microbiota in juveniles tend to mature with age, and the gut microbiota in their early life is more similar to the microbial structure in the environment [6, 16, 66]. Coprophagy and food consumption are the key ways of microbial transmission, which contribute to the maturation of gut microbiota in juveniles. Coprophagy is a common phenomenon in reptiles, leading to gut microbial transmission in the gopher tortoise Gopherus polyphemus [67] and three species of green iguanas, I. iguana [26], P. palluma [68] and P. williamsi [16]. In this study where eggs were individually incubated to avoid the microbial transmission from their mothers and between hatchlings, we only detected differences among gravid females, newly hatched hatchling and postprandial hatchlings within two species. This nding further suggests that the dietary rather than the genetic or maternal contribution has a more important role in affecting the gut microbiota.

Conclusions
By comparing data collected from two closely related snake species, we found that the relative abundance of the gut microbiota changed rapidly in the postprandial hatchlings. This nding suggests diet affects the colonization and proliferation of the gut microbiota in snakes. Interspeci c differences in the relative abundance of the gut microbiota were slight in gravid females, newly hatched hatchlings, and postprandial hatchlings. This suggests that the genetic or maternal contribution has no important role in shaping the gut microbiota in the two closely related snakes.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Availability of data and materials
All 16S rRNA gene sequences obtained in this study have been deposited in the NCBI Sequence Read Archive under the BioProject accession number PRJNA642086.
Authors' contributions YFQ and XJ conceived and designed the experiment. XJ supervised the study. YFQ, YQW and YJJ performed the experiment and analyzed the data. YFQ and XJ wrote the manuscript. All authors have read and approved the manuscript. Figure 1 Relative abundances of the gut microbiota in each fecal sample at the phylum (A), family (B), and genus (C) levels. One color indicates one taxa in each plot. The color for others indicates all other taxa not listed in in Plot B or C. ECA1-5: gravid E. carinata numbered 1 to 5; ECH1-3: newly hatched E. carinata numbered 1 to 3; ECP1-4: postprandial hatchling E. carinata numbered 1 to 4; ETA1-4: gravid E. taeniura numbered 1 to 4; ETH1-3: newly hatched E. taeniura numbered 1 to 3; ETP1-4: postprandial hatchling E. taeniura numbered 1 to 4 Figure 2 Results of principal coordinates analysis of Bray-Curtis distance matrix (A) and the partial least squares discriminant analysis (B) for bacterial relative abundance. ECA: gravid E. carinata; ECH: newly hatched E. carinata; ECP: postprandial hatchling E. carinata; ETA: gravid E. taeniura; ETH: newly hatched E. taeniura; ETP: postprandial hatchling E. taeniura  Differences in bacterial taxa among the six groups are determined by LEfSe (A). LDA scores show the differences in relative abundance among six snake groups (B). Letters p, c, o, f and g indicate phylum, class, order, family and genus, respectively. See Figure 2 for de nitions for ECA, ECH, ECP, ETA, ETH and ETP Figure 4 Relative abundances of gene functions in the gut microbiota for each sample at the top (A), second (B) and third (C) levels. One color indicates one functional classi cation in each plot. Sample grouping information is showed in