Three species of freshwater snails, namely 40 RP (Radix plicatula), 40 BU (Boreoelona ussuriensis) and 36 SQ (Sinotaia quadrata), were collected from the West Liao River (43° 36’ 51” N; 122° 13’ 28” E) at Tongliao, Inner Mongolia, China. All samples were taken from shallow water along the waterline, with the main food source of these snails was benthic algae, such as diatoms, as well as humus from emergent plants, aquatic grasses and decaying plant matter. As Table 1 indicates, the average height of SQ, BU and RP shells was 24.78 ± 1.9 mm, 12.2 ± 0.6 mm and 21.3 ± 1.8 mm, respectively. After collection, the snails were quickly transported to the laboratory where their surfaces were scrubbed with 70% alcohol before extracting their intestines under sterile conditions. In this case, based on the size of the intestinal samples, they were combined in groups of five for RA and BU and groups of three for BP. Therefore, a total of 8 RP samples, 8 BU samples and 12 SQ samples were obtained. Species identification was confirmed by the mitochondrial cytochrome c oxidase subunit I (COI) gene, which was amplified with primers LCO1490 and HCO2198 (Folmer et al., 1994). The gene sequences were registered at NCBI GenBank (Boreoelona ussuriensis: OR990564, Sinotaia quadrata: OR994503 and Radix plicatula: OR994310) ( Supplementary Figure S1 ). This study was approved by the ethics committee of Inner Mongolia Minzu University (IACUC-202211-073).

Microbial genomic DNA was extracted from snail gut samples using the E.Z.N.A.^® soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer’s instructions. Then, the extracted DNA was subjected to 1% agarose gel separation before determining its concentration, and purity was assessed with a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA). For bacterial community analysis, the 16S rRNA genes were amplified using the universal bacterial primers 27F (5’-AGRGTTYGATYMTGGCTCAG-3’) and 1492R (5’-RGYTACCTTGTTACGACTT-3’) (Weisburg et al., 1991) to which PacBio barcode sequences were attached to distinguish between samples. In this case, the PCR reactions were performed in triplicate, and each reaction volume was 20-μL comprising 4 μL of 5 ×FastPfu buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu DNA Polymerase, 10 ng of template DNA and DNase-free water. The amplification program was then set up on a 9700 PCR thermocyclerm (ABI GeneAmp^®, CA, USA) as follows: initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturation at 95°C for 30 s, annealing at 60°C for 30 s and extension at 72°C for 45 s. The reaction ended with a final extension at 72°C for 10 min. Subsequently, the resulting PCR products were checked via electrophoresis prior to purification using AMPure^® PB beads (Pacific Biosciences, CA, USA) and quantification with a Quantus™ Fluorometer (Promega, WI, USA).

After pooling purified DNA products in equimolar amounts, the DNA libraries were constructed using the SMRTbell^® Express Template Prep Kit 2.0 (Pacifc Biosciences, CA, USA), as instructed by the manufacturer. These libraries were subsequently sequenced on the Pacbio Sequel II System (Pacific Biosciences, CA, USA) by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).

PacBio raw reads were processed using the SMRTLink analysis software (version 8.0) to obtain demultiplexed circular consensus sequence (CCS) reads with a minimum of three full passes and 99% sequence accuracy. Then, based on the barcodes, these CCS reads were identified and filtered according to their length. Reads with sequences of <1,000 bp or >1,800 bp removed. Subsequently, the optimized-CCS reads were de-noised using the recommended parameters of the DADA2 (Callahan et al., 2016) plugin in the Qiime2 (Bolyen et al., 2019) (version 2020.2) pipeline to generate amplicon sequence variants (ASVs) at single nucleotide resolution. Furthermore, to minimize the effects of sequencing depth on the diversity metrics, the number of 16S rRNA gene sequences in each sample was rarefied to 13994, with this value yielding an average good coverage of 96.74%. Moreover, the taxonomic assignment of ASVs was performed using the Naive Bayes (Vsearch, or BLAST) consensus taxonomy classifier implemented in Qiime2 as well as the SILVA 16S rRNA database (v138).

Optimized-CCS reads were also clustered into operational taxonomic units (OTUs) at a 97% sequence similarity level using UPARSE 7.1 (Edgar, 2013). After selecting each OUT’s most abundant sequence as a representative sequence, those assigned to chloroplast and mitochondrial sequences were removed from all samples. Then, the taxonomy of the remaining OTU representative sequences was determined using the RDP Classifier version 2.2 (Wang et al., 2007) against the 16S rRNA gene database (Silva v138) at a confidence threshold of 0.7. Additionally, based on OTU representative sequences, the metagenomic function was predicted using PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) (Douglas et al., 2020) based on OTU representative sequences. In this case, BugBase (https://bugbase.cs.umn.edu/index.html) was also used to predict the microbiome phenotype of snails (Ward et al., 2017).

Bioinformatic analysis of the gut microbiota was performed using the Majorbio Cloud platform (https://cloud.majorbio.com). Rarefaction curves and alpha diversity indices, including Simpson index, Shannon index and Sobs index, were then calculated with Mothur v1.30.1 (Schloss et al., 2009) based on the ASVs’ information. Furthermore, based on Bray-Curtis dissimilarity, the similarity between the microbial communities of different samples was determined by principal coordinate analysis (PCoA) using the Vegan v2.5-3 package. For hierarchical cluster analysis, the beta diversity distance matrix was computed in Qiime before generating a corresponding tree using R language (version 3.3.1). The latter was also used for statistics, mapping Venn diagrams, generating bar charts of community composition based on taxonomy summaries as well as for heatmap and bubble map analysis with the Vegan package (version 2.4.3). Additionally, ternary analysis was performed by GGTERN (http://www.ggtern.com/). Lastly, the linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011) (http://huttenhower.sph.harvard.edu/LEfSe) was carried out to identify significantly abundant bacterial taxa (phylum to genera) among the different groups (LDA score > 4, P < 0.05).

Quality-filtering of reads generated a total of 1,030,244 sequences, with an average length of 1472 bases, for 28 samples of freshwater snails (8 samples for BU and RP, respectively, and 12 samples for SQ). Additional filtering was performed to remove sequences related to Chloroplast and Mitochondria as well as those from OTUs having a sequence count less than 5. This step yielded 721,211 sequences which were then clustered into 2,898 operational taxonomic units (OTUs) at a 97% sequence-identity threshold. Subsequent taxonomics classified these OTUs into 26 phyla, 50 classes, 124 orders, 203 families, 434 genera and 776 species.

The gut microbiota of the three species of freshwater snails were then analyzed based on the Sobs, Shannon, and Simpson indexes at the OTU level. The Shannon diversity index for SQ, BU and RP were 6.08, 6.24, and 5.72 respectively, while the observed species (Sobs) count was 1545, 1654, and 1500 respectively. Whereas the Simpson index, the values were 0.008, 0.007, and 0.016 respectively. The different diversity indices suggested that the alpha-diversity of BU was higher than RP (P < 0.05), while no significant differences were noted when comparing the Shannon and Simpson diversity metrics of SQ and BU or SQ and RP (P > 0.05) ( Figure 1 ). The rarefaction and Shannon curves further showed that the sequencing depth was adequate for microbial diversity analysis ( Supplementary Figure S2 ).

PCoA analysis showed significant differences between the intestinal bacterial communities of the three snail species ( Figure 2A : unweighted unifrac analysis of variance, R = 0.99; P=0.001; Figure 2B : weighted unifrac analysis of variance, R = 0.93, P=0.001). Furthermore, the hierarchical clustering analysis divided the samples into three distinct groups that were significantly different based on host species ( Figure 2C ).

The VENN analysis revealed that 566 bacteria, accounting for 72.94% of the total bacteria, were common to the three snails. Of these, 145 were identified up to the species level (25.61%) and included Verrucomicrobia_bacterium_ADurb.Bin474, Limnoraphis_robusta_CS-951, Exiguobacterium_sp._MH3, Chroococcidiopsis_thermalis_PCC_7203 and Exiguobacterium_sp. ZWU0009 ( Supplementary Figure S3 , Supplementary Table S1 ). In addition, 165 bacteria were common to only two snail species; of these species, a level of 47 (28.48%) was identified. More specifically, the species that were common to SQ and BU included Lactococcus_garvieae_ATCC_419156, Lactobacillus_amylovorus, Lactobacillus_pontis, Acinetobacter_sp._CIP_64.7, Burkholderia_stabilis and Streptococcus_alactolyticus, while those common to SQ and RP consisted of cyanobacterium_scsio_7.2, Synechococcus_sp._WH__5701 and Exiguobacterium_sp. Similarly, Scytonematopsiscontorta__HA4292-MV4, Chroococcidiopsis sp. PCC 6712 and Parasegetibacter luojiensis were among the species that were common to BU and RP. Furthermore, 45 bacteria were only present in one snail species ( Figure 2D ). Of these, with only 4 identified up to the species level. For instance, Brevinema_andersonii and Synechocystis_aquatilis_SAG_90.79 were specific to the SQ snail, while Melainabacteria_bacterium_MEL.A1 and Pantoea_vagans_C9-1 were only identified in the BU snail ( Supplementary Figure S4 , Supplementary Table S1 ).

At the phylum level, 10 dominant phyla, with > 1% abundance, were identified, including Cyanobacteria, Proteobacteria, Firmicutes, and Planctomycetota, accounting for 93.47, 86.22, and 94.34% of the total reads of SQ, BU, and RP, respectively ( Figure 3A ). The distribution of microbes at the OTUs level between the three snail species was also determined based on ternary analysis. In this case, the OTUs of Proteobacteria were mostly present in SQ, while the OTUs of Planctomycetota were mainly found in SQ. Similarly, the OTUs of Firmicutes and Cyanobacteria were widely identified in all three snail species ( Figure 2E ).

At the genus level, 33 genera with an abundance of more than 3% were identified, with some examples being Phomidium_IAM_M-71, Leptolyngbya_ANT.L52.2, Exiguobacterium, Romboutsia and Clostridium_sensu_strito_1, Bacillus ( Figure 3B ).

The abundance of bacterial communities in the different snail species was analyzed and visualized using a heatmap and bubble map ( Figure 4 ), with the results showing differences in the composition and distribution of bacteria across the samples.

In particular, the bubble map analysis showed that the top five abundant species of Cyanobacteria mainly occurred in SQ and RP snails, while the top three abundant species from Proteobacteria distributed in RP ( Figure 4A ). Verrucomicrobia_bacterium_ADurb_Bin474, the largest abundant species in Planctomycetota, was only distributed in SQ. In the case of the heatmap which displayed the 50 most abundant OTUs, the results not only reflected the differences in the bacterial composition of the three groups of samples but also highlighted a close similarity between the intestinal microbiota of SQ and RP ( Figure 4B ). In addition, a number of lactic acid bacteria, such as Lactobacillus_rhamnosus_GG, Lactobacillus_amylovorus, Lactobacillus_mucosae_LM1 and Lactobacillus_pontis were found to be common to SQ and BU ( Supplementary Table S1 ).

The intestinal microbiota of SQ, BU and RP were also compared to identify significant differences between them. The results showed that most of the microbiota were significantly different between the three snail species from phylum to the species level ( Supplementary Figures S5 , S6 ).

In particular, LEfSe analysis (LDA threshold > 4) revealed significant enrichment of specific bacterial species in each snail species (Figure 5). These included phyla Planctomycetota, class Planctomycetes, orders Leptolyngbyales, Isosphaeralesin and Alicyclobacillales in SQ samples, orders Rhizobiales and genus Phormidium_IAM_M-71 in RP samples, while class Bacilli, orders Chitinophagales, diplorickettsiaceae and Burkholderiales in BU (Figure 5).

Functional prediction of the intestinal bacterial communities of the three snail species was performed using PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) with the results shown in Figure 6 . Overall, the abundance of the KEGG metabolic pathways was consistent across the three species, with 18 pathways having an abundance of more than 1%. Of these, Global and overview maps were the most abundant, followed by Carbohydrate metabolism, Amino acid metabolism and Energy metabolism, just to name a few. A number of these metabolic pathways, including Global and overview maps, Amino acid metabolism, Energy metabolism and Metabolism of cofactors and vitamin amongst others, were significantly different between the three snail species (P < 0.05).

BugBase predicted a number of microbiome phenotypes for the snails, including biofilm formation, pathogenicity, stress tolerance, aerobic, anaerobic or facultative anaerobic respiration, presence of mobile elements as well as gram-positive or gram-negative features which were all different between the three species (P < 0.05) ( Figure 7A ). More specifically, RP showed significant differences in the abundance of bacteria that contributed to stress tolerance, pathogenicity, biofilm formation, the presence of mobile elements and facultative anaerobic phenotypes (P < 0.05) ( Figure 7A ). Whereas bacteria that contributed to the aerobic phenotypes were significantly more abundant in BU than in RP, although no differences were noted when compared with SQ (P < 0.05) ( Figure 7A ). Finally, anaerobic phenotypes were significantly more abundant in SQ compared with BU and RP (P < 0.05) ( Figure 7A ).

Specific bacterial species were markedly correlated with the different phenotypes. For instance, in RP, unclassified_Enterobacteriaceae and Lelliottia were linked to stress tolerance, biofilm formation, potential pathogenicity, mobile elements, and facultatively anaerobic phenotypes, with a contribution of 28.01% and 23.06%, respectively ( Figures 7B–G ). Similarly, Romboutsia and Clostridium_sensu_stricto_1 respectively contributed 17.01% and 3.70% to the anaerobic phenotype predicted for SQ, while in the case of RP, the degree of contribution was 7.52% and 1.95%, respectively ( Figure 7F ). In addition, norank_f_Isosphaeraceae and Pirellula were correlated with the aerobic phenotype in SQ, contributing 9.43% and 9.01% respectively. In BU, the most significant contribution came from Pirellula, Dinghuibacter, Acinetobacter, Brevundimonas and Stenotrophomonas ( Figure 7H ).

A total of 68 GH (glycoside hydrolase) genes were predicted by PICRUSt, with their abundance accounting for 1.26%, 1.55%, and 1.50% of the total genes in SQ, BU and RP, respectively. Table 2 shows 31 highly abundant genes, of which 21 were significantly different between the three groups, while 10 indicated no significant differences.

Of these, in BU, the relatively abundant genes included 6-phospho-beta-glucosidase (EC 3.2.1.86), Beta-glucosidase (EC 3.2.1.21), Beta-N-acetylhexosaminidase (EC 3.2.1.52), Alpha-glucosidase (EC 3.2.1. 20), Beta-galactosidase (EC 3.2.1. 23), Cellulase (EC 3.2.1.4), Alpha-galactosidase (EC 3.2.1. 22) and Alpha-L-fucosidase (EC 3.2.1. 51) ( Table 2 ). Whereas, their abundance was relatively low in SQ, and intermediate in RP.

In this study, a number of genes encoding carbohydrate-digesting enzymes were found in the three species of snails. These enzymes included monosaccharide-digesting ones, such as Beta-glucosidase, Alpha-glucosidase, Beta-galactosidase, Alpha-galactosidase and Alpha-L-rhamnosidase (3.2.1.40), disaccharide-digesting ones, such as Alpha-mannosidase (3.2.1.24) and Beta-mannosidase (3.2.1.25) as well as those that act on polysaccharides, namely Cellulase (3.2.1.4), Cellulose 1,4-beta-Cellobiosidase (3.2.1.91) and Endo-1,4-beta-xylanase (3.2.1.8). In addition, genes of enzymes that digest starch, including Alpha-amylase (3.2.1.1) and Pullulanase (3.2.1.41), as well as for those that are specific to aquatic plants, such as alpha-trehalase (3.2.1.28) and Chitinase (3.2.1.14), were also present ( Table 2 ).

The identification of pathways related to Xenobiotics biodegradation and metabolism indicated the presence of Drug metabolism, Benzoate degradation, Chloroalkane and chloroalkene degradation, Aminobenzoate degradation, Metabolism of xenobiotics by cytochrome P450 and Naphthalene degradation pathways in all three snail species ( Supplementary Figure S7 ). Furthermore, the presence of pathways related to Infectious diseases suggested that Tuberculosis, Salmonella infection, Legionellosis, Pathogenic Escherichia coli infection and Pertussis could also be abundant in those snails ( Supplementary Figure S8 ).

In this study, SMRT sequencing of full-length 16S rRNA genes was used to compare the gut microbiota of three species of wild herbivorous snails living in the same freshwater environment. Although three snail species, namely Sinotaia quadrata (SQ), Boreoelona ussuriensis (BU) and Radix plicatula (RP), are common in China (Wu et al., 2017; Ye et al., 2021); however, their intestinal microbial communities have not yet been reported. Through the successful identification of a number of bacteria up to the species level, this work yielded distinct gut microbial community profiles for the three snail species.

The sequencing results revealed a high diversity of intestinal microbes in the three snails, as evidenced by the abundance of operational taxonomic units (OTUs) and the alpha diversity metrics. Among the indices, the Shannon diversity index indicated a particular diversity of 6.08, 6.24, and 5.72 for SQ, BU, and RP, respectively. After clustering sequences at a 97% similarity-level threshold, a total of 2898 OTUs (SQ: 2532, BU: 2341, RP: 2443) were obtained for the intestinal microbiota of the three snail species. Furthermore, the alpha-diversity metrics results were higher than that of Pomacea canaliculate with the Shannon and Chao indices of 3.35 and 1106, respectively (Chen et al., 2021).

The complex food sources of herbivorous freshwater snails are rich in indigestible lignocelluloses; therefore, a large diversity of microbial species is required in their digestive tracts to efficiently digest and utilize these substances (Harris, 1993). For example, in cattle, a complex microbiota in the rumen helps to decompose plant fibers, with over 5000 OTUs in the rumen microbiota in second-aged cattle (Jami et al., 2013). The results of this study also demonstrated that the diversity of the gut microbiota was significantly higher in SQ compared with RP. This difference could, in fact, be related to the body size and intestinal morphology of freshwater snails, as SQ species have larger and longer intestines, which allow them to retain a higher diversity of microorganisms.

At the phylum level, Cyanobacteria, Proteobacteria, Firmicutes and Planctomycetota were common to all three snail species. In contrast, a study on the gut microbiota of R. auricularia (Shannon 5.12) reported Proteobacteria, Cyanobacteria, Chloroflexi, Firmicutes, and Actinobacteria as dominant phyla (Hu et al., 2018), while in another study on Planorbella trivolvis, Proteobacteria, Bacteroidetes and Actinobacteria were predominant (Hu et al., 2021). Nevertheless, the results of this study, along with the above reports, suggest that Proteobacteria tend to be a major bacterial group in the intestinal microbiota of freshwater snails (Walters et al., 2022; Zhou et al., 2022). However, this may not always be the case. Indeed, the intestinal microbiota of P. canaliculate, collected from culture ponds, had a high abundance of Tenericutes (45.7%) and Firmicutes (27.85%), while Proteobacteria (11.86%) ranked third in abundance (Chen et al., 2021). These differences could be attributed to variations in the food composition or water environment. Although Planctomycetota was found to be highly abundant in all three snail species and also found in another study of Planorbidae snails (Van Horn et al., 2012), it has rarely been reported in other freshwater snails (Hu et al., 2018; Hu et al., 2021; Liu and Li, 2022; Zhou et al., 2022). In fact, it was only reported in 1 study on the intestinal microbiota of terrestrial Hydrobiidae (Walters et al., 2022).

As far as Cyanobacteria was concerned, it accounted for 784 of the 2,898 total OTUs, and its abundance reached 42, 32, and 42% in SQ, BU, and RP, respectively. Cyanobacteria is a naturally prevalent microbe that is usually attached to fallen leaves, aquatic grasses, humus, and rocky soil. During ingestion, freshwater snails take in a large number of these Cyanobacteria, resulting in a high abundance of this phylum among the intestinal microorganisms of freshwater snails. However, this abundance is dependent on the environmental abundance of Cyanobacteria, with significant changes likely affecting the physiological functions of the host (Lyu et al., 2022).

Identifying the intestinal microbiota of freshwater snails is helpful to accurately understand their role. However, most previous studies were based on the high-throughput sequencing of partial 16S rRNA genes (e.g., the V3-V4 region), resulting in sequence lengths that did not allow accurate identification up to the species level. Here, full-length 16S rRNA genes were sequenced to determine the intestinal microorganisms of three freshwater snails. Yet, limitations in the amount of information available in the database restricted the identification of OTUs to the species level to only 25.62%. Therefore, the addition of comprehensive information in the database is expected to improve accuracy. In addition, the identification rate of bacterial species was also impacted by the availability of research data on freshwater snail intestinal microbiota, hence emphasizing the need for further research in this area.

A total of 196 species were identified in this study, and this provided an understanding of their roles in intestinal microbiota. For example, Tumebacillus_flagellatus, abundant in SQ snails (1.93%), is an alpha-amylase/pullulanase-producing bacterium that plays an important role in starch digestion (Wang et al., 2013). Similarly, Limnoraphis_robusta, a nitrogen-fixing and carotenoid-rich cyanobacterial species, is a common bloom species that does not produce cyanotoxins (Komarek et al., 2013), while Clostridium_estertheticum, common to all three snail species, can secrete a bacteriocin with potential antibacterial activity against Bacillus cereus, Staphylococcus aureus (methicillin-resistant), Escherichia coli and Pseudomonas aeruginosa (Wambui et al., 2022). Therefore, full-length 16S rRNA gene sequencing is crucial for identifying the role of strains in the gut.

Surprisingly, some Lactobacillus species were present in SQ and BU species but not in RP, such as Lactobacillus_rhamnosus_GG, Lactobacillus_amylovorus, Lactobacillus_mucosae_LM1 and Lactobacillus_pontis. Lactobacillus are considered as probiotics that are beneficial for animal digestion and health (Borchers et al., 2002). Lactobacillus rhamnosus, in particular, has been widely studied as a prebiotic for its positive effects on gut health and immune enhancement (Segers & Lebeer., 2014). Similarly, by producing lactate and acetate, Lactobacillus_amylovorus is also an important intestinal probiotic that contributes to gut health (Jing et al., 2022), especially in conferring health benefits, such as anti-hyperammonemia and antiviral properties (Singh et al., 2018; Sunmola et al., 2019). Furthermore, Lactobacillus_mucosae has also been studied for its immunomodulatory or nutritional effects (Ryan et al., 2019; Wang et al., 2022), while Lactobacillus_pontis is commonly found in fermentation processes, such as sourdough or silage (Wu et al., 2014; Wang et al., 2020). Therefore, the intestinal microbiota of freshwater snails acts as a reservoir for probiotic species.

Most of the studies on the intestinal microbiota of mollusks were mainly focused on host species of commercial value or were performed in artificial environments. In contrast, only few involved wild freshwater snails. Therefore, differences in the intestinal microbiota of freshwater snail hosts with similar feeding habits and from the same natural environment are yet to be understood.

The gut microbiome of animals is influenced by many factors, such as the host species, the water environment or food, and therefore, in the case of freshwater snails, much of the experimental data tend to come from artificially controlled environments (Chalifour and Li, 2021; Hu et al., 2021; Li et al., 2022; Zhou et al., 2022). This is because in natural water bodies, it is difficult to model the effects of host species on the microbiota due to a number of naturally occurring uncontrollable factors. Hence, in this work, freshwater snails with similar feeding habits were selected from the same water environment to reduce the influence of environmental and food factors. The results subsequently confirmed that, under similar natural conditions, the intestinal microbiota of freshwater snails was dependent on the host species.

The study of Walters et al. (2022) on the intestinal microbes of two snail species from the same family (Hydrobiidae) showed that both the surrounding water environment and host species influenced the gut microbiota. However, in addition to hosts’ influence, differences in food composition may also affect the gut microbiota, and in this context, studies have shown that the main food source of herbivorous freshwater snails was detritus, followed by decaying macrophytes, diatoms and filamentous algae (Madsen, 1992). For instance, Zhou et al. (2022) fed the same vegetable (Brassica rapa) to two freshwater snails, Pomacea canaliculata and Cipangopaludina chinensis (under laboratory rearing conditions) and studied their intestinal microbiota, which revealed differences between the intestinal microbiota of the two species. Here, although a number of bacteria were found to be common to the three snails, the findings of this study highlighted that their abundance was different. In addition, the dominant species and genera were also different. Thus, altogether, the results suggest that some microorganisms are shared by two or more hosts, while others are specific to a single host.

Bacterial phenotypes determine the functions of gut microorganisms. For instance, biofilms can improve the host’s resistance against pathogen invasion, while the presence of mobile elements can improve the host’s tolerance to drug resistance (Palmer et al., 2010; Stokes and Gillings, 2011). Stress tolerance can further improve adaptability to environmental changes, such as temperature, salinity and contaminants (Fahy et al., 2008; Li et al., 2022; Liu et al., 2023). In this study, significant differences were found between the intestinal microbiome phenotypes of the three snail species, as evidenced from the results of BugBase. In the case of RP, Lelliottia and the genus unclassified_Enterobacteriaceae contributed to phenotypes, such as biofilm formation, potential pathogenicity, stress tolerance, the presence of mobile elements as well as facultative anaerobic characteristics. Lelliottia is a widespread strain in natural environments such as water (Kämpfer et al., 2018) but can also be found in Giant pandas or as an endophyte of desert cactus (Eke et al., 2019; Zhang et al., 2022). Similarly, in addition to its wide distribution in soil and water, Enterobacteriaceae are also largely found in humans, animals and plants (Cohen et al., 2020). Pirellula, assigned to the phylum Planctomycetes, is a chemoheterotrophic aerobe that is partly responsible for the aerobic phenotype in snails. In fact, many Planctomycetes are aerobes and, as such, they may play a central role in the degradation and digestion of heteropolysaccharides from plants and algae (Fuerst, 2017). In contrast, Romboutsia contributed to the anaerobic phenotype of the three snails, especially in SQ and RP. The members of this genus are gram-positive and spore-forming obligate anaerobes that are present not only in human and animal guts but also in natural environments such as coastal estuarine mud and alkaline-saline lake sediment (Gerritsen et al., 2014; Wang et al., 2015; Ricaboni et al., 2016; Maheux et al., 2017).

The intestinal microbiota plays an important role in host metabolism, immunity, digestion, growth and development (Brestoff and Artis, 2013; Lin et al., 2017). In the current study, it was found that the intestinal microorganisms of freshwater snails involved in various important metabolic pathways, such as carbon metabolism, amino acid synthesis, energy metabolism, and metabolism of cofactors and vitamin, just to name a few. These metabolic functions enable microorganisms to participate in food digestion and absorption, protein synthesis, as well as in the supply of energy, vitamins and cofactors in the host’s gut. Consistent with this study’s findings, the microbial community associated with A. fulica was also highly enriched in Carbohydrates Metabolism, Amino Acids Metabolism, Cofactors and Vitamins Metabolism, Protein Metabolism and Membrane Transport (Cardoso et al., 2012).

The food of plant-eating gastropods is rich in cellulose that is difficult to digest. Hence, the gut microbiota secretes large amounts of enzymes to assist in the digestion process. In this study, a large number of genes encoding carbohydrate-digesting enzymes were found to be present in the three species of snails. These included monosaccharide-, disaccharide-, polysaccharide- and starch-digesting enzymes, along with those that are specific to aquatic plants. In addition, cellulase (β-1, 4-glucan-4-glucan hydrolase) (EC. 3.2.1.4) is important for converting cellulose into glucose, and interestingly, the gene encoding this enzyme was highly abundant in the intestinal microbiota of all three snail species, thereby highlighting the importance of these bacteria in helping the host to digest lignocellulose. It has been reported that the intestinal microbiota of the herbivorous snail Achatina fulica is rich in cellulase- and hemicellulose-coding sequences (Cardoso et al., 2012). These cellulose-degrading genes can actually be of value in the development of lignocellulose-derived biofuel.

The enrichment of pathways related to Xenobiotics biodegradation and metabolism indicated that the gut microbiota could play an important role in helping freshwater snails decompose foreign harmful substances such as drugs and pollutants that are discharged through industrial, agricultural or human activities (Liu et al., 2023). In addition, the enrichment of pathways related to Infectious diseases further suggested that some gut microbes could be pathogenic to humans and animals, thereby highlighting their significance for public health.

It should be pointed out that PICRUSt2 has certain limitations in the prediction of microbial function. PICRUSt2 predicts the approximate function of a community based on community membership identified by the 16S rRNA gene; however, it does not provide direct information on microbial genes and their functional composition. Furthermore, the prediction cannot fully recapitulate metagenomics-based functional metabolic profiles, distinguish strain-specific functionality, or identify rare environment-specific functions (Douglas et al., 2020; Gehrig et al., 2022). Due to constraints inherent with 16S sequencing, BugBase has the same limitations as the PICRUSt2, such as lower prediction accuracy, unlike to identify OTUs at the strain level and to distinguish horizontal gene transfer between OTUs (Ward et al., 2017).