Gut microbiota is related to host fitness, and influenced by geographical locations and habitats. Pomacea canaliculata is a malignant invasive alien snail that threatens agricultural production and ecosystem functions worldwide. Clarifying the general rules of the gut microbial community structure and function of the snails in different geographical locations and habitats is of great significance for understanding their invasion at different spatial scales. This study used high-throughput sequencing technology to compare and analyze the differences in community structure and function of gut microbiota in P. canaliculata from five geographical locations (Liuzhou, Yulin, Nanning, Wuzhou, and Hezhou) and three different habitats (pond, paddy field, and ditch) in Guangxi Province. The results showed that the intestinal microbial alpha diversity of P. canaliculata was higher in Liuzhou, Yulin, lower in Nanning, Wuzhou, and Hezhou, and higher in ponds compared with paddy fields and ditches. The dominant phyla of gut microbiota in snails were Firmicutes, Cyanobacteria, Proteobacteria, Fusobacteriota, and Bacteroidota, and the dominant genus was Lactococcus. The community structure of gut microbiota in snails varied significantly across different geographical locations and habitats, and the phyla Firmicutes and Cyanobacteria had significantly higher relative abundance in snails collected from Nanning and Yulin, respectively. Moreover, the relative abundance of gut functional microbiota associated with human disease in P. canaliculata was significantly affected by geographical locations and habitats, and with the highest abundance in ponds. However, the relative abundance of functional microbiota related to metabolism, genetic information processing, organizational systems, environmental information processing, and cellular processes was only significantly affected by geographical locations. Collectively, geographical locations and habitats had significantly different effects on the community structure and function of gut microbiota in P. canaliculata, and the greater differences were caused by geographical locations rather than by habitats.

Sample collection

In August to September 2022, we collected P. canaliculata in the field from 5 geographical locations (Nanning, NN (107°77′E, 23°09′N); Liuzhou, LZ (109°31′E,  24°37′N); Yulin, YL (109°99′E, 22°32′N); Wuzhou, WZ (110°30′E, 23°54′N); Hezhou, HZ (111°67′E, 24°35′N)) in Guangxi Province, China (Figure 1A). Three sites containing 3 habitats (pond, paddy field, ditch, Figure 1D, E, F) simultaneously were randomly selected in each geographical location, and 5 quadrats (1m²) were set in each habitat (Figure S3). One adult snail was taken from each quadrat for intestinal sample collection, and the distance between the sample quadrat was about 10 meters. We distinguished between male and female when sampling the snails and 23 female snails and 22 male snails were collected from each geographical location. A total of 225 P. canaliculata were collected in this study (5 geographical locations Í 3 sites Í 3 habitats Í 5 replicates). All the testing snails collected from five geographical locations were preliminarily discerned by shell morphological analysis (Hayes et al., 2012) and using primers LCO1490 or HCO2198 to amplify cytochrome C oxidase subunit I (COI) gene to identify P. canaliculata (Yang et al., 2019) which could be used to sequence for gut microbiota. The body weight, shell height (Table S2), shell width, and shell mouth width of each P. canaliculata were also measured. All sampling individuals were wiped with 75% ethanol three times and followed by rinsing twice in distilled water to sanitize the surface prior to dissection. The entire intestinal contents were extracted carefully to avoid rupturing the gut wall. Each sample was stored in a sterile tube using liquid nitrogen and later stored in a freezer of -80°C.

DNA extraction and sequencing

The total genomic DNA (gDNA) of each sample was extracted using the cetyltrimethylammonium bromide (CTAB) method (Allen et al., 2006). The V3-V4 hypervariable region of the 16S rDNA genes was amplified using specific bacterial primers 341F (CCTAYGGGRBGCASCAG) and 806R (GGACTACNNGGGTATCTAAT) by polymerase chain reactions (PCRs). All PCR mixtures contained 15 µL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 µM of each primer and 10ng target DNA, and cycling conditions consisted of a first denaturation step at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, primer annealing at 50°C for 30 s and extension at 72°C for 30 s, with a final extension step carried out at 72°C for 5 min to ensure complete amplification. The PCR products were purified with a Qiagen Gel Extraction Kit (Qiagen, Germany). Sequencing libraries were generated with NEBNext® Ultra™ IIDNA Library Prep Kit (Cat No. E7645) following the manufacturer’s recommendations, and the library quality was evaluated on the Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. The library was sequenced on an Illumina NovaSeq platform.

Statistical and bioinformatics analyses

Firstly, paired-end reads were assigned to samples based on their unique barcodes and were truncated by cutting off the barcodes and primer, and merged using FLASH (Version 1.2.11). Quality filtering on the raw tags was performed using the fastp (Version 0.20.0) software to obtain high-quality clean tags which were compared with the SILVA 123 database using Vsearch (Version 2.15.0) to detect the chimera sequences, and the chimera sequences were removed to obtain the effective tag. Denoise of the effective tags was performed with DADA2 to obtain initial Amplicon Sequence Variants (ASVs), and then ASVs with an abundance of less than 5 were filtered out. Secondly, species annotation was performed using QIIME2 software (Version QIIME2-202006) based on the SILVA 123 database, and multiple sequence alignment was performed to study the phylogenetic relationship of each ASV and the differences of the dominant species among different samples. Finally, all samples were rarefied to the sequencing depth of the lowest sample (26970 clean reads).