Gut bacterial communities in the freshwater snail Planorbella trivolvis and their modification by a non-herbivorous diet

Rearing of snails

The snails were originally purchased from the ornamental fish market, and a breeding colony were formed in the laboratory. When a sufficient number of newborn snails were reared to about 5–7 mm of shell diameter, a total of 120 juvenile snails with uniform size (5.88 ± 0.14 mm of shell diameter) were selected and randomly distributed into two groups by their diets: herbivorous groups (HV) (60 snails) and non-herbivorous groups (NHV) (60 snails). Each group had three replicate aquariums with 20 snails in each repetition. HV group was fed dry alfalfa contain a high content of cellulose (26%), and NHV group was fed market-purchased pellet feed (Takara sakana-ii, for ornamental fish feeding), which contained low cellulose (5%) (Table 1). The aquaculture water temperature was maintained at 24~25 °C, with a pH of 6.7 ± 0.4. The snails were reared for 60 days and no snails died during rearing.

Sample collection and processing

After 60 days of rearing, the body length and weight of snails were determined first. Next, the intestines were sampled. When the intestines were collected, the snails were cleaned with tap water, and the surface of the shell was subsequently wiped with 70% alcohol. The shells were carefully broken and disassembled from the snail body. After that step, the intestines were extracted carefully to avoid rupturing the gut wall. Next, the intestines were sampled and contained their contents. All dissections were conducted under aseptic conditions. For histological examination, the collected intestines were fixed in 4% formaldehyde for histological examination. Because the intestines of individual snails were too small to meet the sampling requirements, the intestines and contents of the 5–8 snails were packed into a centrifuge tube to form a repeat sample for sequencing.

Growth index

The growth performance of snails after 60 days of rearing was characterized by their weight gain rate (WGR) and specific growth rate (SGR). The WGR was calculated by comparing their weight gain (WG) with their initial average weight (IW). The SGRs were obtained by dividing their weight gain by 60 days. The formulas of WGR and SGR are as follows:

WGR (%) = [(TW − IW)/IW] × 100

SGR (%/d) = [(TW − IW)/d] × 100

WGR: Weight gain rate; SGR: Specific growth rate; TW: terminal weight; IW: initial weight; d: rearing days.

Intestinal morphology analysis

To examine the effect of diets on the intestinal morphology of snails, the intestinal samples were washed in saline solution and fixed in 4% paraformaldehyde solution above 48 h. Next, the fixed samples were successively treated as follows: dehydrated by an alcohol gradient (50%, 70%, 90% and 100%), xylene transparency, and paraffin embedding. Three paraffin sections from each sample were cut and stained by the hematoxylin-eosin (HE) method. Last, the stained slice was sealed by neutral gum and subjected to observation and scenery shooting using a panoramic scanner (3DHISTECH; Pannoramic MIDI, Budapest, Hungary).

Determination of intestinal microbiota diversity

The intestine samples obtained as described above contained the contents used for genome extraction. Total DNA was extracted using the FAST DNA™ Spin Kit for soil (MP-BIO, Santa Ana, CA, USA) following the protocol of the manufacturer. The extracted DNA was subjected to 1.5% agarose-gel electrophoresis, and the DNA concentration and OD260 nm/OD280 nm (OD, optical density) value were determined with a Nanodrop spectrophotometer (Thermo Scientific, Wilmington, DE, USA).

To detect intestinal bacteria, the V3~V4 regions of the bacterial 16S rRNA gene were amplified with DNA primers 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Mori et al., 2013). The 12 bp barcodes on both ends of the primers were used to identify the sequences of different samples. Next, PCR amplification was conducted, and the procedure was as follows: initial denaturation at 95 °C for 3 min followed by 27 cycles of denaturing at 95 °C for 30 s, annealing at 55 °C for 30 s and extension at 72 °C for 45 s, single extension at 72 °C for 10 min, and ending at 4 °C. The amplified products were subjected to electrophoresis and extracted from a 2% agarose gel. Next, the products were purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). After Qubit quantification and detection, amplicons from each PCR sample were normalized to equimolar amounts, and high-throughput sequencing was performed using the 2 × 300 bp protocol on an Illumina MiSeq PE300 platform (Illumina, San Diego, CA, USA).

The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: PRJNA640745).

Raw fastq files were demultiplexed, quality-filtered by Trimmomatic (Bolger, Lohse & Usadel, 2014) and merged by FLASH. UPARSE (Edgar, 2013) was utilized to conduct operational taxonomic unit (OTU) clustering analysis at 97% identity, and chimaeric sequences were identified and removed using UCHIME (Edgar et al., 2011). Next, the representative OTUs were analyzed on the Qiime platform (Caporaso et al., 2010) against the Silva_132 16S Database (http://www.arb-silva.de) to determine taxonomy.

The heatmap analysis was performed in the heatmap package in R. Venn and PCoA analyses were performed in the vegan package in R. The alpha-diversity of the gut bacterial community, containing indices of Sobs, Shannon, Chao, Simpson, Ace, and coverage, was analyzed by Mothur (Schloss et al., 2009). Biomarker discovery was performed on the linear discriminant analysis (LDA) effect size (LEfSe) to identify the specific organisms whose relative abundances differ between two groups of samples (Segata et al., 2011).

To examine functional differences in the microbiome of snails, the functional attributes of metabolic genes from the snail gut microbiota (KEGG Orthologs, KOs) were predicted by a phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt), which were obtained by macthing the sequencing data against the genomic KEGG database (Langille et al., 2013). As a computational approach, PICRUSt can predict the functional composition of an obtained metagenome using marker gene data and a database of reference genomes (Langille et al., 2013).

Statistical analysis

SPSS 19.0 software was employed for a two-tailed Student’s t-test on the growth data of snails. P values < 0.05 were considered to indicate significant differences, and descriptive statistics were expressed as the mean ± SD.

Snail growth

The initial size and weight did not differ between the two groups (P > 0.05) (Table 2). After 60 days of rearing, the diameter of the final shell of the HV snails was about 10.53 mm, representing an increase of 176.08% compared with the initial shell diameter (Table 2). The final shell diameter of NHV snails was about 11.86 mm, representing an increase of 205.19% compared with the initial shell diameter; thus, and the final shell diameter of snails in the NHV group was significantly larger than that of the HV group (P < 0.05). The endpoint weight of HV group was 0.219 g, and that of NHV group was 0.432 g. The weight gain rate of the NHV group (566.90%) was significantly higher than that of the HV group (321.40%) (P < 0.05). Also, the specific growth rate of snails in the NHV group (3.16%) was significantly higher than that in the HV group (1.94%) (P < 0.05) (Table 2).

Intestinal morphology

To examine the effect of feed on the intestinal morphology of snails, histological examination were conducted (Fig. 1). The differences between treatment groups were clearly shown in the simple columnar epithelium cells (SCEC). For the HV snails, the intestinal structure was intact, and the SCEC arrangement was regulated. The microvilli were tightly attached to the SCEC. For the NHV snail, the intestinal SCECs were filled with fat droplets and thus appear transparent in this image (Figs. 1C and 1D). The fat droplets occupied the cytoplasm of the SCEC, and the nucleus was observed to be indented (Figs. 1C and 1D), indicating that low fiber feeding leads to increased fat deposition in intestinal epithelial cells. The microvilli were detached to a certain extent from the SCEC.

Alpha diversity based on 16S rRNA gene sequencing

A total of 737,284 raw reads were generated using the Illumina MiSeq sequence platform and 485,342 high quality sequences were obtained (following quality control and sequence filtration) with an average length of 459.9 bp. The mean number of sequences per sample was 54,277.60 ± 6,298.01 in HV snails and 39,606.60 ± 6,846.06 in NHV snails (Table 3). The 31,486 rarefied sequences were clustered into 525 OTUs at 97% identity, with 486 and 195 OTUs being identified in HV and NHV, respectively.

The obtained 525 OTUs include members of 419 species, 304 genera, 184 families, 108 orders, 37 classes and 19 phyla. There were 18 phyla, 286 genera, 108 orders, 36 classes, 486 OTUs belonging to HV, and there were 13 phyla, 52 orders, 21 classes, 132 genera, and 195 OTUs belonging to NHV. All of the alpha-diversity indices (Sobs, Shannon, Simpson, Ace, Chao) were significantly different between HV and NHV populations (P < 0.05) (student’s t-test) (Table 3). The Good’s coverage for the observed OTUs was above 99.70%, and the rarefaction curves showed clear asymptotes (Fig. 2A), which, taken together indicated that the given level of sampling effort was sufficient to capture the bacterial communities in snails. Notably, there were clear differences for the rarefaction and rank-abundance curves between the treatment groups (Figs. 2A and 2B). The abundance analysis showed that there were only 5 OTUs in the NHV snails with abundances greater than 1%, while there were as many as 22 OTUs in the HV snails (Fig. 2B). The total abundance of OTUs with abundance greater than 1% reached 79.68% in the HV snails and 83.39% in the NHV snails. The number of core OTUs identified in all samples was 22 (Fig. 2C).

Venn analysis showed that the bacteria belonging to 12 phyla, 77 families, and 114 genera were shared by two groups: 6 phyla, 97 families, and 172 genera were unique in HV snails, and 1 phylum, 10 families, and 18 genera were unique in NHV snails (Figs. 2D and 2F).

Composition of snail intestinal microbiota

The phyla Proteobacteria (52.97%) and Bacteroidetes (28.75%) were predominant among the bacteria in the HV group (Fig. 3A). The phyla with abundances above 1% were Actinobacteria, Verrucomicrobia, Chloroflexi, Cyanobacteria, and Chlamydiae. However, only Proteobacteria (95.23%) was predominant among the bacteria in NHV snails. Bacteroidetes was the other phylum with an abundance above 1%.

At the family level, there were 20 families and 8 families with abundances above 1% in the HV and NHV snails, respectively. Weeksellaceae (24.6%) and Rhodobacteraceae (12.38%) showed high abundance in HV snails. Aeromonadaceae (85.40%) dominated the microbiota in NHV snails (Fig. 3B).

At the genus level, the composition of intestinal microbiota differed between the HV group and the NHV group (Fig. 3C). In the HV group, there were 22 genera with abundances above 1%. Cloacibacterium (24.60%) showed the highest abundance in the HV groups. Others, such as Pseudomonas (4.47%), OM6ON (6.12%), and Rhodobacter (3.78%) were also abundant. In the NHV snail, there were 7 bacteria with abundances higher than 1%. Aeromonas, exhibited the highest abundance, accounting for 85.40%.

The heat map analysis constructed from the top 35 abundant genera, reflected that the intestinal microbiota in the two groups was different (Fig. 4). The Wilcoxon rank-sum test showed that the abundance of 26 genera was significantly higher in HV snails than in NHV snails (P < 0.05). Only two genera (Aeromonas and Comamonas) with significantly higher abundances were observed in NHV snails.

Clustering of the snail gut bacterial community

Clustering analysis of intestinal microbiota associated with diet were performed. Principal coordinates analysis (PCoA) of sequencing data using pairwise weighted and unweighted UniFrac distances showed that the bacterial community structure of HV was different from that of NHV (Fig. 5). UniFrac clearly separated different microbiota efficiently by diet, demonstrating the importance of food as a driver of community composition in this freshwater snail.

LEfSe analysis (the LDA threshold was 4) was used to screen out microorganisms with significant differences in the snail intestines of the two groups fed different diets. There were 11 bacterial genera that were significantly enriched in the intestinal samples of the HV group, and there was one genus (Aeomonas) that was enriched in the NHV group (Fig. S1). When the LDA threshold was set to 2, 114 bacterial genera were significantly distributed in the intestinal samples of the HV group, while only 4 bacterial genera were significantly distributed in the NHV group (Table S1).

Although the intestinal microbiota of snails was influenced by food type, we found that 27 genera existed as core microbiota of all samples in both groups. Among these genera, 19 were affiliated with Proteobacteria, 4 were affiliated with Bacteroidetes, 3 were affiliated with Actinobacteria, 2 were affiliated with Verrucomicrobia and 1 was affiliated with Chloroflexi (Table 4).

Predictive metagenomic profiling

Phylogenetic investigation of communities by reconstruction of unobserved states analysis was performed to generate a predictive functional profile to gain insight into the metabolic capacity of the enteric microbiome (Fig. 6). Based on this analysis, the relative abundance of several COG functions was significantly different between feeding regimens, such as the enrichment of genes involved in energy production and conversion, lipid transport and metabolism in HV, and the enrichment of genes involved in amino acid transport and metabolism, signal transduction mechanisms in NHV. No significant differences were observed for the relative abundance of genes involved in carbohydrate transport and metabolism, inorganic ion transport and metabolism, and transcription.

Because plant fiber degradation requires a diverse suite of enzymes for complete hydrolysis, the abundance of genes for several lignocellulose-active enzymes was also examined (Fig. 7). Eighteen enzymes determined to be more enriched in HV snails than in NHV snails, acting upon cellulose, hemicellulose, lignin, and cello-oligosaccharides, respectively. Six enzymes were more enriched in NHV than HV snails, acting upon hemicellulose, and cello-oligosaccharides, respectively (Fig. 7). Enzymes related to protein degradation were also analyzed. The result showed that a large number of genes related to aminopeptidase, dipeptidase, carboxypeptidase and endopeptidase were more enriched in HV snails than in NHV snails (Fig. 8).