Phylogenetic Analysis of Bacterial Communities in Different Regions of the Gastrointestinal Tract of Agkistrodon piscivorus, the Cottonmouth Snake

Vertebrates are metagenomic organisms in that they are composed not only of their own genes but also those of their associated microbial cells. The majority of these associated microorganisms are found in the gastrointestinal tract (GIT) and presumably assist in processes such as energy and nutrient acquisition. Few studies have investigated the associated gut bacterial communities of non-mammalian vertebrates, and most rely on captive animals and/or fecal samples only. Here we investigate the gut bacterial community composition of a squamate reptile, the cottonmouth snake, Agkistrodon piscivorus through pyrosequencing of the bacterial 16S rRNA gene. We characterize the bacterial communities present in the small intestine, large intestine and cloaca. Many bacterial lineages present have been reported by other vertebrate gut community studies, but we also recovered unexpected bacteria that may be unique to squamate gut communities. Bacterial communities were not phylogenetically clustered according to GIT region, but there were statistically significant differences in community composition between regions. Additionally we demonstrate the utility of using cloacal swabs as a method for sampling snake gut bacterial communities.


Introduction
Vertebrates are metagenomic organisms; they are not only composed of their own genetic material, but also that of their associated microbial communities [1]. The majority of these microorganisms are found in the host intestinal tract, and presumably assist in essential processes of energy and nutrient acquisition [2]. The ecological and evolutionary forces that act on both the host and it's trillions of resident microorganisms sculpt the endogenic microbiome. With the advent of next generation sequencing technologies we are now better able than ever to characterize this observed microbial diversity. However, most studies investigating evolutionary patterns in vertebrate gut microbiomes have focused on mammals [1,2] and even among these studies, many have used captive animals from zoos or farms rather than wild populations. Very few studies have examined the gut microbiome of squamate reptiles (snakes, lizards), despite this being one of the most diverse and successful vertebrate clades.
The cottonmouth (Agkistrodon piscivorus, Serpentes, Viperidae) is a semiaquatic snake widespread throughout southeastern United States. The ecology and demographic history of A. piscivorus has been well studied, and it is often used as a model system in studies of venom evolution [3][4][5][6]. Though A. piscivorus is considered a generalist predator, preying upon reptiles, birds and invertebrates, the diet is dominated by amphibians and fish [6,7]. The paucity of information on the microbiome of wild vertebrate gastrointestinal tract (GIT) regions renders studies of any species meritorious, and exploration of the gut microbiome of A. piscivorus is particularly interesting as almost all other aspects of its ecology and biology are well known. Given the extent of knowledge on this organism's natural history, inferences regarding factors influencing the composition of its GIT microbial community should be possible once that community has been well characterized. In this study we examine the bacterial communities of the small and large intestines, and cloaca of eight individuals of adult A. piscivorus. Our results reveal the presence of distinct bacterial community composition in each GIT region, provide novel insights into the diversity of squamate reptile associated bacterial communities in the wild, and demonstrate the utility of using non-lethal cloacal swabs to sample this diversity.

Ethics Statement
Adult snakes were collected and sacrificed in accordance with IACUC protocols approved by the Committee for Animal Care and Use at the University of Mississippi (#13-02 & #13-04). Euthanasia was performed using an overdose of the anesthetic lidocaine, injected into the brain (UM IACUC SOP #13-02). Collecting permits were obtained from the Mississippi Department of Wildlife, Fisheries, and Parks (permit #'s 0827101 & 1009112).

Microbial Sampling
Snakes were sampled from two sites in Winston County (one individual, N32.98463 X W088.9980) and Lafayette County (seven individuals, N34.427238 X W089.38631), Mississippi, in spring of 2011 and spring of 2012 (Table 1). All snakes sampled were encountered during nighttime surveys along small streams leading to larger bodies of water. Snakes were collected by hand and safely restrained with clear plastic tubing placed over the head during sample collection. Once restrained, snakes were palpated to evaluate whether prey items were present in the GIT then the exterior cloaca of the snake was cleaned using a sterile alcohol pad. This sterilization step was to ensure that the cloacal sample primarily included cloaca associated microbes rather than environmental or transient microbes. Following cleaning, cloacal swabs were collected by inserting a sterile polyester-tipped applicator (Fisher Cat# 23-400-122) into the cloaca, taking care not to insert beyond the coprodeum and into the large intestine, then turning the swab several times before withdrawing. Once withdrawn the applicator was immediately placed into a sterile 2 ml tube and placed on ice before being transferred to a -20°C freezer prior to DNA extraction. Three individuals (samples 103, 110, and 111) were sampled in more detail to determine the bacterial communities of their small and large intestines. These snakes were transported to the Department of Biology at the University of Mississippi where they were humanely euthanized. Immediately following death, a mid-ventral incision was made to expose the GIT, which was then removed. None of the individuals had identifiable prey items present in the GIT. Incisions were made in proximal and distal ends of both the small and large intestines, which were then swabbed with sterile polyester-tipped applicators that were immediately placed in sterile 2 ml collection tubes and frozen (-20°C) until DNA extraction. The remainder of the snake was preserved in 10% buffered formalin and whole specimens were deposited at the Sam Noble Oklahoma Museum of Natural History (Table 1).

DNA Extraction
Microbial DNA was extracted by a bead beating procedure using MoBio Power Soil Extraction kits (MoBio Laboratories, Carlsbad, CA, USA). We followed the manufacturer's standard DNA extraction protocol with minor adjustments. Thawed applicator tips were placed into bead tubes containing lysis buffer, and 50-100 μl of the lysis buffer was used to rinse any remaining particles that may have become dislodged from tips out of the 2ml collection tube and into the bead tube. Additional adjustments included incubating samples at 65°C for 15 min after the addition of Solution C1, and vortexing bead tubes horizontally for 25 min.

PCR Amplification and Analysis
We used a nested PCR approach and bacterial specific primers to amplify a variable region of the 16S rRNA gene for initial analysis of the gut bacterial community by denaturing gradient gel electrophoresis (DGGE). This nested approach was necessary as DGGE can only be performed on fragments <500bp, and our template DNA was of low quantity rendering direct amplification of small fragments difficult. We first amplified near full-length fragments of the bacterial 16S rRNA gene using primer sets Bac8f and Univ1492r, and used this as our template for amplifying a shorter region for use in DGGE with the primer sets Bac1070f and Univ 1392GCr. Primers, PCR protocols, and cycle conditions have been previously described [8]. Amplification products from the second round of amplifications (bases 1070-1392) were analyzed in DGGE gels using a 40% to 70% denaturant gradient, and electrophoresis for 20 h at 80 V. Following electrophoresis, gels were stained with SYBR Green I and visualized by UV transillumination using a Kodak Gel Logic 200 system running Molecular Imaging Software 4.0 (Eastman Kodak, Rochester, NY, USA). Banding patterns were converted to binary data based on presence or absence of specific bands in each sample. Binary data was then used to create a distance matrix, showing similarity between samples (Yue & Clayton theta index [9]), and these relationships were visualized by ordination of samples through non-metric multidimensional scaling (NMDS). All analyses were performed using the bioinformatics software Mothur [10].
Pyrosequencing DGGE analyses revealed that the proximal and distal ends of the large and small intestines possessed near identical banding patterns, therefore we combined the DNA extraction from these regions such that each dissected individual had a single small and single large intestinal sample analyzed by pyrosequencing, along with the cloacal samples from both dissected and released individuals. The initial (bp 8-1492) 16S rRNA amplicons from each sample were sequenced using bacterial tag-encoded FLX amplicon 454 pyrosequencing (bTEFAP) [11] at the Research and Testing Laboratory sequencing facility (Lubbock, TX). Library amplification was performed with the bacterial primers 939f and 1392r [8,12] under the following conditions: 95°C for 5 minutes, followed by 35 cycles at 95°C for 30s, 54°C for 40s, an extension at 72°C for 1 minute, followed by a final elongation of 10 minutes at 72°C. Sequencing was performed on a Roche 454 FLX titanium instrument using standard reagents and following manufacturer's guidelines. Sequences were deposited in the NIH NCBI Sequence Reads Archive (SRA Accession numbers SAMN03287547-SAMN03287560).

Sequence Analysis
Pyrosequence data was accessed and processed in the program Mothur [10] following general procedures recommended by Schloss et al. [13]. Following denoising and barcode removal, sequences were aligned using the greengenes reference database (http://greengenes. secondgenome.com/, May 2013 version) and sequences differing by only a single nucleotide were grouped together. Sequences were checked for chimeras using the chimera.slayer command in Mothur and potential chimeras were removed. The resulting sequences were classified according to the greengenes database and any non-bacterial 16S sequences were removed. The remaining sequences were clustered into operational taxonomic units (OTUs) based on 97% similarity.
The distribution of OTUs in each sample was used for analyses of diversity and comparisons of community structure. Because samples varied in the number of final valid sequences obtained, all analyses were performed on a randomly sampled subset of the total dataset for each sample, which corresponded to the number of sequences in the smallest sample (i.e. all samples were standardized to be equivalent to the sample with the lowest number of valid reads). This random subsampling was performed 1,000 times for each analysis, with the composite outcome reported. Alpha diversity within each sample was determined by Schao and inverse Simpsons indices, and rarefaction and collection curves were used to visualize whether our procedures included enough sampling (enough reads) to assess this diversity. Beta diversity (comparisons of bacterial community structure between samples) was examined using the Yue & Clayton theta similarity index [9] which accounts for proportional abundance of OTUs in a sample. Similarity between samples was visualized by NMDS, Venn diagrams, as well as dendrogram construction. We tested the spatial separation of samples observed in NMDS through analysis of molecular variance (AMOVA) and analysis of similarities (ANOSIM). Dendrograms were constructed based on Yue & Clayton theta (thetayc) distances and a 95% majority rule consensus tree was generated from the distribution of 1000 trees without burn in using the program TreeAnnotator [14].

DGGE Analyses
While based on a limited number of samples from just two locations, NMDS of community similarity based on DGGE binary data showed a clear pattern of distinct small intestine, large intestine, and cloacal communities, with only slight overlapping of multidimensional space between the large and small intestine samples (Fig 1A). There was no apparent association in GIT bacterial community structure among individuals or localities, although our ability to test this was limited by the number of individuals sampled. All regions of the GIT show similar levels of diversity and richness as inferred from DGGE banding patterns. Band number, used as a proxy for richness [15], ranged from 8-29, with a mean of 17.

Sequence Analyses
After alignment, gap removal, and potential chimera removal we recovered 40,317 valid sequence reads, representing 7,137 unique sequences with a median length of 267 bp from our pooled dataset. The classified sequences were binned into OTUs and rarefaction suggested that we recovered >98% of the diversity in all but one of our samples, and with that sample (individual 103l) we recovered 95% of the diversity ( Table 2).  Bacterial diversity recovered in our 7,137 unique sequences was binned into 503 distinct OTUs spanning 14 bacterial phyla, with <0.002% (just 92 sequences) of the 40,317 sequence dataset designated as "unclassified Bacteria". Among individuals for which all three regions were sampled the large intestine harbored the most diverse bacterial communities in two of the three individuals (Table 2). In terms of community composition, the large intestine samples were dominated by sequences affiliated with the Bacteroidetes, followed by Firmicutes, Proteobacteria, and Lentisphaerae, while members of the Proteobacteria were the dominant group in both the small intestine and cloaca samples, followed by sequences classified as Firmicutes and Bacteroidetes (Fig 2). Gammaproteobacteria were the dominant subphylum of Proteobacteria in all but two samples; one small intestine (110) sample and one cloacal sample (111) had Deltaproteobacteria and Betaproteobacteria, respectively, as their dominant Proteobacteria. Importantly, the dominant bacterial phyla found in both the small and large intestines were also found in cloacal samples (Fig 2).
NMDS plots based on pyrosequence data did not show as tight of a grouping of samples by GIT region as was produced by DGGE analyses (Fig 1A and 1B). ANOSIM and AMOVA tests of spatial separation did suggest a significant difference in community composition among all three regions as a group (AMOVA: df = 2 F s = 1.902 p = 0.005, ANOSIM: R 2 = 0.606 p = 0.001), but none of the individual pairwise comparisons between different GIT regions were significant. Next generation sequencing is likely to detect ultra-rare species that are typically not detected by DGGE [16], and these may not have an important role in the GIT community, or could represent transitory cells that are not permanent members of the gut community. Therefore we also analyzed the pyrosequence data after sequentially removing the rarest 1%, 5% and 10% of sequences. NMDS of communities with the rarest 1% and 5% of sequences removed resembled the structuring seen in our DGGE data, although that signal was lost once the rarest 10% of sequences are removed (Fig 1C (1%), D (10%), 5% removed not shown in figure for clarity). AMOVA following the removal of the rarest 5% of sequences still suggested overall community differences among GIT regions (p = 0.007), without any individual pairwise comparisons between GIT regions being significant. However, ANOSIM did detect significant differences for all pairwise comparisons of GIT regions once any level (1%, 5%, or 10%) of the rarest sequences was removed (Table 3).
Because we had cloacal samples from both dissected and released specimens, we focused on the cloaca as our primary GIT region of interest, as this would allow future sampling to avoid sacrificing individuals. Furthermore, of the 22 OTUs representing greater than 1% of the sequence reads, all but four were found in the cloaca as well as other regions (Table 4). Based on proportions of sequence reads obtained, the most abundant bacterial lineage in the cloaca varied among individuals: Members of the Bacteroidetes were dominant in two samples, Proteobacteria (primarily subphyla Gamma and Beta) were dominant in two, similar levels of abundance of Bacteroidetes and Firmicutes were found in three, while one sample had similar proportions of Bacteroidetes, Proteobacteria and Firmicutes (Fig 3). At a finer taxonomic scale, the most abundant sequences affiliated with Bacteroidetes were classified as members of the genus Bacteroides, with two different species dominating in the large intestine and cloaca (Table 4). While the most abundant Proteobacteria (specifically Gammaproteobacteria) sequences obtained from the cloaca were classified as either Salmonella enterica, or as unclassified Enterobacteriaceae, those of the small intestine were identified as Pseudomonas veronii and Aeromanadaceae (Table 4). Both Enterobacter and Acinetobacter were found in all regions of the GIT. The most numerous Firmicutes sequences were in the family Clostridiaceae and those we were able to classify to genus were Clostridium (Table 4).
Venn diagrams revealed that some OTUs in each region of the GIT were shared by all three individuals sampled, but the majority of OTUs in each region were unique to a given individual  4). Dendrogram analysis did not reveal a definitive pattern of clustering by GIT region sampled or individual, although cloacal samples tended to group together (Fig 4).

Discussion
The aim of this study was to characterize the gut microbiome of a non-captive squamate reptile species. Only one other study has examined bacterial community structure in the GIT of wild snakes [17]. That study included two individuals of A. piscivorus, the species studied here, but the authors focused on the most abundant bacterial populations in the gut community as determined from DGGE followed by traditional Sanger sequencing [17], rather than incorporating the high throughput sequencing approach that we used. Despite the difference in methodology, that study also found the Firmicutes and Bacteroidetes to be the most dominant phyla in the GIT of A. piscivorus, although no samples of the cloaca were taken and, presumably because of methodological constraints, the study did not attempt to assess diversity or abundance. [17].
Our study demonstrates the utility of both standardized sampling via cloacal swabs and the greater information in community composition that can be obtained by using next generation sequencing rather than traditional methods of bacterial community analyses. Both DGGE and next generation sequencing recovered distinct bacterial communities corresponding to discrete regions of the GIT (small intestine, large intestine, cloaca), although this was only apparent when the more common sequence types were considered. As might be expected given the cloaca's terminal location, cloacal samples captured the bacterial diversity found in earlier regions of the GIT, but not necessarily the proportional abundance of specific taxa. While the dominance of members of the Bacteroidetes in the large intestine is consistent with other findings in snakes [17,18], the dominance of members of the Proteobacteria in the small intestine and cloaca was unexpected. Previous work has suggested that the small intestine of captive snakes is dominated by members of the Firmicutes and Bacteroidetes, regardless of physiological processes (e.g. active digestion) [18]. Our findings suggest this may not be the case in wild snakes. It is well known that the mammalian gut is dominated by members of the Bacteroidetes and Firmicutes [19] and the limited amount of previous work on wild reptiles (e.g. marine iguanas) and captive snakes suggests that this is the same for reptiles [18,20]. However, a recent meta-analysis compared the GIT bacteria of insects, birds and mammals and found that the bacterial community in insects and birds is more likely to be dominated by Proteobacteria and Firmicutes rather than Bacteroidetes [21]. The dominance of Proteobacteria in the GIT of snakes in this study suggests that predatory snakes may show more similarities to birds in terms of their GIT bacterial communities than to other vertebrate organisms.
In one of the few previous studies to investigate bacterial community structure in wild squamates, Hong et al. found that sequences affiliated with the Firmicutes dominated fecal samples from both marine and land iguanas in the Galapagos, with the classes Bacteroidia and Clostridia being dominant in both host species [20]. They suggested that the similarity between these herbivorous reptiles to mammals in the general composition of their GIT bacterial assemblages is because of the prevalence of herbivory in mammals, and that abundance of members of the Bacteroidetes and Firmicutes in the gut aids in the breakdown and nutrient acquisition of complex polysaccharides [2,20]. Wood eating fish species have also been shown to contain similar GIT diversity to herbivorous mammals [22] and these similarities (at the level of phylum) have led some authors to suggest that fish served as the original vertebrate hosts to these bacterial communities [23,24]. Although A. piscivorus is a predatory species (no snakes are herbivorous), Clostridia were the dominant Firmicutes recovered in our analysis, possibly reflecting a broader function of this subset of Firmicutes in the GIT that is not restricted to herbivory. It may be that this class of obligate fermentative bacteria is so well adapted to the anoxic conditions of the intestines that they are likely to be dominant in the GIT of any larger organism. Cottonmouth snakes are generalist predators and are known to feed on carrion; and some authors suggest that certain populations of A. piscivorous my acquire food primarily by scavenging [25]. Many Clostridia are associated with the fermentation of animal material and this may convey benefits to host species that prey on or scavenge other vertebrates.
All of the individuals sampled in our study were adults, and no detectable prey items were found in dissected individuals or through palpitation of non-dissected individuals. That said, the broad variation in dominant bacterial OTUs we recovered from cloacal samples leads us to suspect that abundance may be influenced by active digestion of prey items that were undetected. Although the bacterial species present in the GIT has been shown to remain constant throughout digestion and subsistence in captive snakes, the relative abundance of these species can be correlated with active physiological processes such as digestion [18]. An additional factor relating to diet variation in the gut bacterial community that has yet to be explored is ontogenetic change. Diet has been shown to change through ontogeny in A. piscivorus [26], and it's not unreasonable to suspect that this could lead to changes in bacterial community structure. Acquisition of a new diet has been hypothesized to be a fundamental driver for species' diversification and concomitant gut microbiota evolution [19] and dietary changes during ontogeny should also drive changes in gut microbiota. Lastly, it has been shown that in mammalian hosts which undergo torpor during hibernation both bacterial community composition and richness vary dramatically during hibernation [27]. It would be interesting to explore whether this holds true for ectothermic species that undergo hibernation, which A. piscivorus does in certain areas of its range.
Our study lays the groundwork for further investigation of gut bacterial community structure in squamate reptiles. We identified discrete bacterial communities that correspond to regions along the GIT and show that cloacal samples encompass the breadth of bacterial diversity found in the snake gut. Additionally we show the utility of our standardized sampling method using cloacal swabs. Interestingly, and in contrast to previous investigations of snakes, we find that at the phylum level the snake gut microbiome shows more similarities to that of birds than to other vertebrates. Future studies should include broader sampling of host species for more detailed comparative analyses, and test whether gut bacterial communities function as either evolutionary or ecological traits.
with permits and the Mississippi Department of Wildlife, Fisheries and Parks for granting collecting permits.