Effects of Bellamya Stocking Density on Water Quality and Bacterial Community Responses in Aquaculture Effluent

1. Introduction

Aquaculture is now a central component of global aquatic food production, but production intensification also increases the generation of nutrient-rich discharge water. Residual feed, excretion, suspended particles and algal biomass can introduce nitrogen, phosphorus and organic matter into receiving waters, thereby increasing the risk of eutrophication and ecological deterioration [1,2,3]. Recent reviews emphasize that aquaculture effluents are difficult to manage because they often occur in large volumes but at relatively low pollutant concentrations; consequently, treatment systems must balance removal efficiency, land demand, energy consumption and operational cost [2,3]. Biological and nature-based approaches are therefore increasingly considered as complements to conventional physical and chemical treatment technologies [2,3,4,5].

Microorganisms are central to the fate of nutrients in aquaculture waters. Ammonia oxidation, nitrite oxidation, denitrification and anammox collectively determine the transformation and removal of inorganic nitrogen, while microbial processing of dissolved and particulate organic matter can influence phosphorus availability and algal-bacterial interactions [6,7,8,9,10]. These processes are environmentally sensitive and may become uncoupled under changing oxygen, carbon and nutrient conditions. In particular, NO₂⁻-N can accumulate when ammonia oxidation proceeds faster than nitrite oxidation or subsequent reduction processes [6,7]. Accordingly, assessment of ecological effluent treatment should combine conventional water-quality indicators with analysis of bacterial community structure and taxon-environment relationships.

Freshwater gastropods can modify aquatic habitats through grazing, particle ingestion, biodeposition and bioturbation. Mesocosm studies have shown that snails can suppress phytoplankton or epiphytic algae and can interact with submerged vegetation to improve water clarity and nutrient conditions [11,12]. Within the genus Bellamya, recent studies have demonstrated effects on sediment organic-matter degradation, nitrogen forms, bacterial community composition and cross-domain microbial interactions [13,14,15,16,17]. Other work has shown that the gut microbiota of Bellamya aeruginosa responds to cyanobacterial blooms and that snail-associated bacterial communities differ from surrounding epiphytic communities [18,19]. These findings indicate that Bellamya can influence both material cycling and microbial niches, but most studies have focused on sediments, host-associated microbiota or plant–snail systems rather than directly evaluating aquaculture effluent. This gap supports an integrated assessment of Bellamya treatment that combines water-quality outcomes with bacterial-community responses.

Recent aquaculture studies also show that biologically based treatment and integrated culture systems can couple changes in nutrient conditions with shifts in microbial communities [20,21]. Although these systems are not direct analogues of Bellamya treatment, they support evaluating organism-based interventions using both water-quality and microbial-community endpoints.

High-throughput 16S rRNA gene sequencing provides an efficient approach for characterizing community diversity, composition and treatment-responsive taxa. Ordination and PERMANOVA can test whole-community differences, LEfSe can identify discriminative taxa, and correlation or constrained ordination can reveal associations between bacterial groups and environmental gradients [22,23,24,25,26,27,28,29]. However, taxonomic profiles and predicted functions should be interpreted conservatively because microbial functions are not uniquely determined by taxonomy and 16S-based functional inference is not equivalent to metagenomic or transcriptional measurement [30,31,32].

The present study therefore investigated the effects of Bellamya stocking density on aquaculture effluent using an integrated water-quality and bacterial-community framework. We hypothesized that Bellamya stocking density would alter nutrient and algal-biomass indicators and that these changes would be accompanied by shifts in bacterial community structure. Specifically, we aimed to (i) evaluate stocking-density effects on SS, inorganic nitrogen, reactive phosphate and Chl-a; (ii) characterize bacterial alpha and beta diversity and taxonomic composition; (iii) identify treatment-responsive taxa and representative genera; and (iv) assess associations between selected bacterial genera and water-quality variables using Spearman correlation and redundancy analysis (RDA).

2. Materials and Methods

2.1. Experimental Design and Sample Collection

The experiment comprised five experimental groups: a no-snail control without Bellamya (CON) and four Bellamya treatments designated as LD, MD, MHD and HD. Each group contained three independent experimental units (n = 3). The CON treatment contained no Bellamya, whereas 15, 25, 30 and 65 individuals were introduced into the LD, MD, MHD and HD treatments, respectively. All experimental units received aquaculture effluent from the same source and were maintained under identical conditions except for Bellamya stocking density.

The experiment was conducted in transparent acrylic tanks measuring 30 cm × 17 cm × 20 cm. Each tank was filled with 10.2 L (0.0102 m³) of aquaculture effluent collected from a culture pond in Songjiang, Shanghai, China. Stocking numbers, initial total wet weights, water volume and biomass stocking densities are summarized in Table 1. Snails were collected from a crab pond, and individuals of similar size were selected. Before the experiment, the shell surfaces were rinsed, and the snails were held for 24 h in tap water followed by 24 h in filtered pond water to reduce residual gut contents and acclimate them to the experimental water. The experiment lasted 60 days, from 1 September to 30 October 2025, and the water temperature ranged from 15 to 25 °C.

At the end of the experiment, water was collected independently from each experimental unit for physicochemical analysis and 16S rRNA gene sequencing. For microbial DNA sampling, 5 L of water from each tank were sequentially filtered through 5.0-μm and 0.22-μm membranes. The membrane samples were collected aseptically and used for subsequent DNA extraction.

2.2. Water-Quality Measurements

Suspended solids (SS) were determined gravimetrically according to GB 11901-1989. Briefly, 50 mL of thoroughly mixed water was filtered through a 0.45-μm membrane, and the retained material was dried at 103–105 °C to constant weight.

Ammonium nitrogen (NH₄⁺-N) and nitrate nitrogen (NO₃⁻-N) were determined according to HJ 535-2009 and HJ/T 346-2007, respectively, using a DT-3900 multiparameter water-quality analyzer (Beijing Huamei Water Analytical Instrument Technology Co., Ltd., Beijing, China).

Nitrite nitrogen (NO₂⁻-N) was measured spectrophotometrically according to GB 7493-1987. Nitrite was diazotized with sulfanilamide and coupled with N-(1-naphthyl)-ethylenediamine dihydrochloride, and absorbance was measured at 540 nm.

Reactive phosphate (PO₄³⁻-P) was determined using the molybdenum blue method according to GB 17378.4-2007. Water samples were filtered through a 0.45-μm membrane, reacted with ammonium molybdate and ascorbic acid, and absorbance was measured at 882 nm.

Chlorophyll a (Chl-a) was determined spectrophotometrically according to HJ 897-2017. Algal biomass was retained on glass-fiber filters, extracted with 90% acetone, and absorbance of the extract was measured at 750, 664, 647 and 630 nm. Each independent experimental unit was treated as one biological replicate.

2.3. DNA Extraction, PCR Amplification and Sequencing

Total microbial DNA was extracted from the membrane samples using a HiPure Water DNA Kit (Magen Biotechnology Co., Ltd., Guangzhou, China; Cat. No. D3145-02) according to the manufacturer’s instructions. DNA concentration and purity were assessed using a NanoDrop 2000 spectrophotometer, and DNA integrity was evaluated by 1% agarose gel electrophoresis at 5 V cm⁻¹ for 20 min.

The V3–V4 region of the bacterial 16S rRNA gene was amplified using barcode-tagged primers 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). PCR amplification was performed in a 20-μL reaction mixture containing 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 Polymerase, 0.2 μL of bovine serum albumin (BSA), approximately 10 ng of template DNA, and nuclease-free water to a final volume of 20 μL.

The thermal cycling program consisted of an initial denaturation at 95 °C for 3 min; 30 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 45 s; followed by a final extension at 72 °C for 10 min and holding at 10 °C.

Three independent PCR reactions were performed for each sample, and the resulting amplicons were pooled. PCR products were examined by 2% agarose gel electrophoresis and purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). Purified amplicons were quantified using a Quantus™ Fluorometer (Promega, Madison, WI, USA) and pooled in appropriate proportions according to the required sequencing depth.

Sequencing libraries were prepared using the Hieff NGS^® Ultima DNA Library Prep Kit for Illumina^® (Yeasen Biotechnology, Shanghai, China) according to the manufacturer’s instructions. Paired-end sequencing (2 × 250 bp) was performed on an Illumina MiSeq platform.

2.4. Sequence Processing and Taxonomic Annotation

Raw paired-end reads were demultiplexed according to sample-specific barcodes. Quality filtering was performed using the sequencing-provider pipeline with Trimmomatic [33], FLASH and USEARCH. Reads were truncated when the mean quality score within a 10-bp sliding window fell below 20, and reads shorter than 50 bp after trimming were discarded. Paired-end reads were merged using FLASH [34] with a minimum overlap of 10 bp; merged sequences with an overlap mismatch ratio greater than 0.2 were removed. Barcode mismatches were not allowed, and no more than two primer mismatches were permitted. Chimeric sequences were removed using USEARCH in combination with de novo and reference-based detection against the Gold database [35].

High-quality sequences were dereplicated, singleton sequences were removed, and the remaining unique sequences were clustered into operational taxonomic units (OTUs) at 97% sequence similarity using USEARCH v7.0.1090 and the UPARSE workflow [23]. All quality-filtered sequences were subsequently mapped back to representative OTU sequences at ≥97% similarity to generate the OTU abundance table. Taxonomic assignment of representative sequences was performed using RDP Classifier v2.2 implemented in QIIME v1.9.0 against the SILVA Release 132 database [22,24,25], with a confidence threshold of 0.7. All downstream taxonomic analyses were therefore based on 97%-similarity OTUs rather than amplicon sequence variants.

The sequencing-provider workflow retained both the original and normalized OTU tables. The normalized OTU table was obtained by rarefying all samples to the minimum sequencing depth across the dataset and was used for taxonomic composition, beta-diversity, differential-taxon and taxon–environment analyses. Alpha-diversity indices and rarefaction curves were obtained from the provider-generated 97%-similarity OTU output.

2.5. Alpha Diversity, Community Structure and Taxonomic Composition

Alpha diversity was evaluated at 97% OTU similarity using mothur v1.41.0 and included ACE and Chao1 richness estimators, Shannon and Simpson diversity indices, and Good’s coverage [26,36,37,38]. Rarefaction curves based on observed OTU richness were generated to assess sequencing sufficiency. Chao1 and Shannon indices were selected for visualization in the main text, whereas the complete alpha-diversity results are provided in Supplementary Table S1.

Between-sample community differences were quantified using Bray–Curtis dissimilarity [39]. The Bray–Curtis distance matrix generated by the sequencing-provider pipeline was used for non-metric multidimensional scaling (NMDS) in the R package vegan, and differences among the five treatment groups were tested by permutational multivariate analysis of variance (PERMANOVA/adonis) with 999 permutations [27].

Bacterial community composition was summarized from the normalized OTU table at the phylum and genus levels. Dominant phyla were displayed in the main text. The top-20 genus profile was retained as a supplementary overview because a large proportion of low-abundance genera was grouped as “Others”. Five clearly classified genera—Rhodobacter, Ahniella, Novosphingobium, Sphingobacterium and Silanimonas—were selected from significant LEfSe biomarkers for direct sample-level abundance visualization.

2.6. Differential Taxa, Taxon–Environment Associations and Functional Prediction

Differentially enriched bacterial taxa among treatments were identified using linear discriminant analysis effect size (LEfSe) [28]. The final five-group analysis was performed using the Galaxy implementation of LEfSe. The Kruskal–Wallis test was used to screen taxa showing significant differences among groups, followed by linear discriminant analysis to estimate the effect size of each discriminative feature. Taxa with P < 0.05 and an LDA score ≥ 3.0 were considered discriminative biomarkers. The term “taxa” is used because the significant LEfSe output included named genera together with uncultured, unclassified, norank and higher-level lineages.

For focused exploratory taxon–environment analysis, a fixed panel of 15 clearly classified genera was used consistently for both Spearman correlation and RDA. The panel was defined before taxon–environment testing from named genus-level taxa that were detected across multiple samples, showed measurable sample-level abundance, and exhibited treatment-responsive patterns in the abundance or LEfSe results. No genus was added or removed on the basis of the resulting correlation or RDA P values. Because the panel was not defined by a single universal abundance or prevalence cutoff, these analyses were treated as exploratory rather than exhaustive. The selected genera were Rhodobacter, Sphingomonas, Novosphingobium, Bradyrhizobium, Ahniella, Flavobacterium, Acinetobacter, Diaphorobacter, Comamonas, Sphingobacterium, Pelomonas, Tabrizicola, Silanimonas, Labrys and Hirschia. Spearman rank correlations were calculated across all 15 sample-level observations between genus-level relative abundances and SS, NH₄⁺-N, NO₂⁻-N, NO₃⁻-N, PO₄³⁻-P and Chl-a. Raw P values were used to annotate the exploratory correlation heatmap, whereas Benjamini–Hochberg false-discovery-rate (FDR) adjustment was applied jointly across all 90 genus–environment tests; the corresponding adjusted q values are reported in Supplementary Table S3.

Redundancy analysis (RDA) was performed using the interactive redundancy-analysis module of the GenesCloud platform with the sample-level relative abundances of the same 15 selected genera and the six measured water-quality variables as input. The significance of the constrained ordination was evaluated using the permutation test implemented by the platform. Potential bacterial functional profiles were inferred using PICRUSt2 [29]. The default PICRUSt2 reference genome database, derived from 41,926 bacterial and archaeal genomes in the Integrated Microbial Genomes (IMG) database (snapshot dated 8 November 2017), was used for prediction. Predicted gene-family abundances were assigned to KEGG Orthology (KO) and COG functional categories using the annotations integrated with this IMG-derived reference genome set; no separate standalone KEGG or COG release was specified by the sequencing provider. These predictions were treated as functional inference only and not as direct measurements of functional-gene abundance, transcription or pathway activity [30,31,32].

2.7. Statistical Analysis

Univariate statistical analyses were performed using IBM SPSS Statistics version 27.0 (IBM Corp.), and figures were prepared using GraphPad Prism version 9.5.0 (GraphPad Software). Data from independent experimental units are presented as mean ± standard error of the mean (SEM), with the three biological replicates shown individually in the figures. Normality was assessed using the Shapiro–Wilk test and homogeneity of variance using Levene’s test. When parametric assumptions were satisfied, differences among the five groups were analyzed using one-way analysis of variance followed by Tukey’s multiple-comparison test. Otherwise, the Kruskal–Wallis test followed by Dunn’s multiple-comparison test was used. Alpha-diversity indices were compared using the Kruskal–Wallis test. Multivariate community analyses were performed using the sequencing-provider pipeline and the R package vegan where indicated. PERMANOVA, Spearman correlation analysis and RDA were performed as described above. For the correlation analysis, both raw P values and Benjamini–Hochberg FDR-adjusted q values were retained; significance stars in the heatmap represent raw P values, while adjusted q values are provided in Supplementary Table S3. Different lowercase letters in univariate figures indicate significant pairwise differences among groups, whereas groups sharing a letter are not significantly different. Statistical significance was set at P < 0.05.

3. Results

3.1. Bellamya Stocking-Density Treatments Improved Several Water-Quality Indicators

Bellamya stocking density altered the water-quality profile of aquaculture effluent (Figure 1). SS showed an overall downward tendency across treatment groups, although the difference was not statistically significant. NH₄⁺-N decreased significantly in all Bellamya-stocked groups compared with the control, with the lowest level observed in the HD treatment. NO₃⁻-N was also significantly reduced by all treatments, and the LD group showed the lowest concentration. PO₄³⁻-P was lower in the HD and MHD groups than in the control, while the LD group showed greater within-group variation. Chl-a was significantly reduced in the HD and MHD groups. In contrast, NO₂⁻-N showed an upward tendency in several treatment groups, particularly HD, but the overall group difference was not statistically significant. These results indicate that Bellamya stocking-density treatments improved several eutrophication-related variables while potentially allowing transient nitrite accumulation.

3.2. Sequencing Coverage and Alpha Diversity

The 16S rRNA gene sequencing dataset contained 60,006–91,528 valid reads per sample and 1,119–1,334 observed OTUs. Good’s coverage was consistently high (99.64%–99.80%), and rarefaction curves approached saturation, indicating adequate sequencing coverage (Table S1; Figure S3). Chao1 richness differed among treatments (P = 0.023), with lower values in MHD and higher values in LD and CON. Shannon diversity also differed among groups (P = 0.043); LD maintained relatively high diversity, whereas MHD showed the greatest within-group variation (Figure 2).

3.3. Bellamya Stocking Density Significantly Reshaped Bacterial Community Structure

NMDS based on Bray-Curtis dissimilarity showed clear treatment-related separation of bacterial communities (Figure 3). The ordination had a stress value of 0.0822, indicating a reliable two-dimensional representation. PERMANOVA confirmed a significant stocking-density treatment effect on community composition, with the stocking-density factor explaining 62% of the observed variation (R² = 0.62, P = 0.001). CON, HD and LD samples formed relatively distinct clusters, whereas MD and MHD showed greater within-group dispersion.

3.4. Community Composition and Representative Bacterial Genera

At the phylum level, the bacterial community was dominated by Proteobacteria, Actinobacteriota, Bacteroidota, Planctomycetota, Verrucomicrobiota, Chloroflexi and Cyanobacteria, with treatment-related differences in relative abundance (Figure 4A). The top-20 genus profile contained a large “Others” fraction, demonstrating substantial contribution from numerous low-abundance genera (Figure S1).

Five LEfSe-supported genera were visualized to clarify treatment responses at genus level (Figure 4B–F). Rhodobacter was more abundant in CON than in all Bellamya treatments. Ahniella was relatively enriched in CON and LD, whereas HD and MHD showed lower abundance. Novosphingobium increased in the treatment groups and reached its highest abundance in MHD, while CON had the lowest abundance. Sphingobacterium was selectively enriched in HD. Silanimonas was most abundant in CON and was lower in HD and MHD, with MD and LD showing intermediate values.

3.5. LEfSe Identified Treatment-Specific Discriminative Taxa

LEfSe identified 25 discriminative taxa among the five groups using P < 0.05 and LDA score ≥ 3.0 (Figure 5; Table S2). CON was characterized by taxa including Rhodobacter, Ahniella and Silanimonas. MHD was characterized by Novosphingobium and several additional lineages, whereas HD was characterized by Sphingobacterium and other treatment-specific taxa. Because some biomarkers were annotated as norank, unclassified or uncultured lineages, the LEfSe results were interpreted at the level of discriminative taxa rather than exclusively as named genera.

3.6. Selected Genera Were Associated with Water-Quality Gradients

Spearman analysis revealed pronounced associations between the selected genera and water-quality variables (Figure 6). The strongest and most consistent patterns involved NH₄⁺-N and Chl-a. Genera including Sphingomonas, Novosphingobium, Bradyrhizobium, Acinetobacter, Diaphorobacter, Comamonas, Sphingobacterium and Pelomonas generally showed negative associations with NH₄⁺-N and Chl-a, whereas Rhodobacter, Ahniella, Silanimonas, Tabrizicola, Labrys and Hirschia showed the opposite pattern. Several treatment-responsive genera were positively associated with NO₂⁻-N, consistent with the contrasting behavior of nitrite relative to the other water-quality indicators. Several strong associations with NH₄⁺-N and Chl-a remained significant after Benjamini–Hochberg correction (Table S3).

RDA further demonstrated a significant multivariate association between the selected genera and the measured environmental variables (permutation test, P = 0.043; Figure 7). RDA1 and RDA2 explained 55.73% and 6.36% of the constrained variation, respectively. NH₄⁺-N, PO₄³⁻-P and Chl-a pointed toward the region occupied mainly by CON samples, whereas NO₂⁻-N showed an opposing orientation and was more closely aligned with several Bellamya-treated samples.

3.7. PICRUSt2-Predicted Functional Profiles

PICRUSt2-based prediction suggested that the bacterial community was dominated by broad functional categories related to metabolism, membrane transport, genetic-information processing, energy metabolism, carbohydrate metabolism and amino-acid metabolism (Figure S2). These results provide a broad functional context for the observed taxonomic changes but were not used as direct evidence of functional-gene abundance or pathway activity.

4. Discussion

4.1. Bellamya Stocking-Density Treatments Improved Nutrient- and Algal-Related Water-Quality Conditions

The present study demonstrates that Bellamya stocking-density treatments can improve several water-quality variables that are directly relevant to aquaculture effluent management. The consistent reduction in NH₄⁺-N and NO₃⁻-N, together with lower PO₄³⁻-P and Chl-a in the HD and MHD treatments, indicates that the presence of Bellamya was associated with a shift toward lower nutrient and algal-biomass pressure. This finding is consistent with the broader rationale for biological treatment of aquaculture discharge water, in which animals, plants and microbial consortia are used to transform, retain or assimilate waste nutrients with lower energy demand than intensive physicochemical systems [2,3,4,5].

The water-quality response is also consistent with previous ecological studies of freshwater snails. Snail grazing can reduce epiphytic or planktonic algae and can indirectly improve light conditions and macrophyte performance [11,12]. Studies specifically involving Bellamya purificata have reported altered nitrogen pools, organic-matter degradation and microbial community structure in sediments [13,14,15,16,17]. Although the present study focused on the water column rather than sediment, the direction of the response supports the view that Bellamya can modify aquatic nutrient and microbial environments through multiple, interacting processes. These may include direct ingestion of suspended particles and algae, biodeposition, excretion and indirect modification of microbial niches. The relative contribution of each process cannot be separated by the current design and should therefore not be presented as a proven mechanism.

The absence of a significant SS response despite a downward numerical trend further illustrates why an ecological treatment should not be judged by a single endpoint. Large within-group variation, especially in some treatment groups, can obscure treatment effects when replication is limited. In contrast, the clearer responses of dissolved inorganic nitrogen, reactive phosphate and Chl-a suggest that Bellamya primarily influenced nutrient and biological processes rather than producing a uniform response across all water-quality variables. Similar endpoint-specific responses are common in biological aquaculture treatment systems, where pollutant removal depends on substrate form, microbial activity and food-web interactions [2,3,4,5].

4.2. Nitrite Accumulation Suggests Asynchronous Nitrogen Transformation

NO₂⁻-N behaved differently from NH₄⁺-N and NO₃⁻-N. Although the treatment effect was not statistically significant, several Bellamya groups showed higher nitrite concentrations than the control. Nitrite is an intermediate shared by nitrification, denitrification and other nitrogen-transformation pathways, and its accumulation can occur when the rates of production and consumption are temporarily unbalanced [6,7,8]. The present pattern therefore does not negate the reductions in ammonia and nitrate; instead, it suggests that the nitrogen network may have shifted without reaching a new steady state by the sampling endpoint.

Recent work on Bellamya purificata also indicates that snail density and bioturbation can alter nitrogen forms and bacterial communities in ways that differ among nitrogen pools [14,15]. In those systems, changes in dissolved oxygen, nitrate and nitrite were closely related to shifts in bacterial community composition [15]. The current result is compatible with such asynchronous responses, but direct attribution to nitrification or denitrification would require functional-gene or process-rate measurements. A conservative interpretation is therefore that Bellamya treatment modified the conditions under which nitrogen transformations occurred, producing lower ammonia and nitrate together with a possible transient nitrite bottleneck.

4.3. Bacterial Community Restructuring Was a Major Response to Bellamya Stocking Density

The bacterial-community response was one of the strongest findings of the study. NMDS and PERMANOVA showed pronounced separation among treatments, indicating that the ecological treatment affected the community as a whole rather than only a small number of taxa. Aquaculture microbial communities are highly responsive to nutrient availability, organic-matter quality, algal production and animal-associated habitat modification [4,5,8,9,10]. The significant community shift observed here is therefore ecologically plausible and provides an important bridge between the measured water-quality responses and taxonomic changes.

The dominant phyla in the present effluent samples are commonly reported in freshwater and aquaculture environments. Proteobacteria and Bacteroidota often respond to labile organic substrates, whereas Actinobacteriota, Planctomycetota, Verrucomicrobiota and Chloroflexi occupy diverse roles in freshwater carbon and nutrient cycling. Importantly, the large “Others” fraction at genus level indicates that many low-abundance taxa collectively contributed to community structure. This supports the decision to use the top-20 genus plot as supplementary context rather than as the main ecological interpretation. Microbial-ecology studies increasingly recognize that low-abundance taxa can be disproportionately important in interaction networks and biogeochemical responses, and that community function cannot be inferred simply from the most abundant organisms [9,32,40].

The observed community response is also in line with studies of Bellamya-associated systems. Bellamya purificata cultivation has been shown to alter sediment bacterial composition and function, while the presence of the snail can reshape interactions among bacteria, fungi and protists [15,17]. The gut microbiota of Bellamya aeruginosa also changes in response to cyanobacterial bloom conditions, highlighting the ecological sensitivity of snail-associated microbial systems [18]. Together with the present water-column results, these studies suggest that Bellamya can influence microbial communities across multiple habitat compartments, although the direction and magnitude of change depend on the system, stocking level and environmental context.

4.4. Representative Genera Link Treatment Responses to Environmental Gradients

The genus-level analysis provided a more focused view of treatment responses. Rhodobacter, Ahniella and Silanimonas were associated mainly with the control condition, whereas Novosphingobium and Sphingobacterium were enriched in specific Bellamya treatments. These patterns were supported independently by sample-level abundance plots and LEfSe. The concordance between abundance trends and biomarker analysis strengthens their value as treatment-responsive indicators, but it does not establish direct functional roles in effluent purification.

The correlation structure nevertheless provides useful ecological context. Several genera that increased under treatment were negatively related to NH₄⁺-N and Chl-a and positively related to NO₂⁻-N. Members of Comamonas are widespread in aquatic and engineered environments, and some strains have documented denitrification capacity, while community-level studies have linked the genus with nitrogen-removal performance [41,42]. Free-living Bradyrhizobium lineages have also been isolated from freshwater sediments and can possess nitrogen-related metabolic potential [43,44]. However, genus-level 16S assignments encompass substantial functional diversity, and functional properties reported for one species or strain cannot be automatically transferred to all sequences assigned to the same genus. The present associations should therefore be interpreted as hypotheses about ecological response groups rather than proof of specific pathways.

The RDA result reinforces this interpretation at the multivariate level. NH₄⁺-N, PO₄³⁻-P and Chl-a were oriented toward the control-associated community, whereas NO₂⁻-N showed an opposing gradient. This pattern closely mirrors the univariate water-quality results and indicates that the selected genera responded coherently to the environmental changes generated by treatment. Similar use of RDA and taxon-environment correlations has revealed strong environmental control of Bellamya-associated microbiomes and aquaculture-water bacterial communities [10,15,18]. Nevertheless, shared responses to treatment can generate correlations even when no direct causal relationship exists. For this reason, the present study uses terms such as “associated with”, “responsive to” and “potentially involved” rather than “drove” or “mediated”.

4.5. Functional Inference, Study Limitations and Application Perspective

PICRUSt2 suggested that metabolism, membrane transport, energy metabolism and carbohydrate- and amino-acid-related functions were prominent in the bacterial community. These categories are broadly compatible with an aquatic microbial community responding to changes in nutrient and organic-matter conditions. However, 16S-based functional inference is inherently indirect [29]. Taxonomic turnover does not necessarily imply equivalent functional turnover, because different taxa can share functions and closely related taxa can differ substantially in gene content [30,31,32]. Accordingly, the predicted functional profile was treated as supplementary evidence and was not used to claim activation of specific nitrogen- or phosphorus-cycling genes.

The present study has several limitations. First, each treatment contained three independent replicates, which is acceptable for a controlled microcosm experiment but limits statistical power and makes visualization of individual data points essential. Second, the experiment was evaluated at a single endpoint, so temporal succession of water quality and microbial communities could not be resolved. Third, the study did not directly measure nitrogen-cycling genes, enzyme activities or process rates; therefore, microbiome data support ecological associations rather than a complete mechanistic pathway. Fourth, the relative contributions of grazing, biodeposition, excretion and microbial habitat modification were not separated. Similar limitations are common in controlled mollusk studies and environmental bioremediation experiments, where laboratory designs provide control but do not fully reproduce the complexity of field systems [11,12,13,14,15,16,17].

Despite these limitations, the study provides a coherent application-oriented dataset: multiple water-quality indicators improved, the bacterial community shifted significantly, treatment-responsive taxa were identified, and the direction of taxon-environment relationships was consistent with the measured environmental gradients. This level of evidence is appropriate for evaluating Bellamya as a potential ecological component of aquaculture effluent management. Future work can test long-term stability, field-scale performance and functional mechanisms, but these questions extend rather than invalidate the current findings.

5. Conclusions

Bellamya stocking-density treatments improved several eutrophication-related water-quality indicators in aquaculture effluent, with clear reductions in NH₄⁺-N and NO₃⁻-N and lower PO₄³⁻-P and Chl-a at selected stocking densities. NO₂⁻-N showed a contrasting upward tendency, suggesting that nitrogen transformation may have been asynchronous or incomplete at the sampling endpoint. Bacterial alpha diversity changed among treatments, and whole-community structure was significantly reshaped. Representative genera and LEfSe biomarkers showed treatment-specific abundance patterns, while Spearman correlation and RDA linked the selected genera to the measured nutrient and algal-biomass gradients. Collectively, these results support Bellamya stocking-density management as a promising ecological approach for improving aquaculture effluent quality and reorganizing bacterial communities. The microbial results indicate community-level associations and potential ecological involvement rather than direct proof of functional mechanisms.