Native Bacillus-Based Probiotic Consortia Suppress Vibrio parahaemolyticus and Restructure Hatchery Water Microbiomes in Shrimp Larval Systems

1. Introduction

Shrimp aquaculture has expanded rapidly over the last two decades, yet recurrent disease outbreaks driven by opportunistic bacteria continue to constrain productivity, profitability, and sustainability. Among the most consequential threats, acute hepatopancreatic necrosis disease (AHPND) has been tightly linked to Vibrio parahaemolyticus lineages that became virulent after acquiring a plasmid encoding PirAB-like toxins, which can trigger acute mortality and severe hepatopancreatic lesions in penaeid shrimp. Foundational work demonstrated that pathogenicity can emerge when V. parahaemolyticus acquires a virulence plasmid expressing a lethal toxin [1], and subsequent studies established molecular targets for detection and quantification of Pir-like toxin genes and their association with AHPND outbreaks [2,3,4,5]. Beyond PirAB-mediated pathology, AHPND-causing strains may also maintain antibacterial type VI secretion systems with diverse effector repertoires, enhancing ecological competitiveness within dense aquaculture microbiomes [6]. Whole-genome sequencing has further clarified geographic origins and dissemination of outbreak-associated lineages, reinforcing the dynamic evolution and mobility of virulence determinants in shrimp production systems [7]. Importantly, experimental evidence indicates that PirABVP toxins can aggravate broader vibriosis contexts, emphasizing that management must address both pathogen abundance and virulence potential across the culture cycle [8,9].

Vibrios are naturally abundant in coastal and marine systems and exhibit substantial genetic and ecological diversity, which makes hatchery water and larval reservoirs interfaces where background Vibrio populations can expand under intensive husbandry, feed inputs, and fluctuating physicochemical conditions. Comprehensive syntheses describe the biodiversity of vibrios across ecological niches and the drivers of their emergence in aquatic food webs [10]. For V. parahaemolyticus specifically, pathogenesis is multifactorial and depends on host susceptibility, environmental conditions, and the accessory genome, which carries a range of fitness determinants and virulence factors beyond Pir toxins [11]. Consequently, disease mitigation strategies must target not only the presence of Vibrio spp. but also the ecological context that enables pathogenic dominance, persistence, and expression of virulence determinants.

Sustainable disease control in aquaculture increasingly prioritizes preventive approaches that reinforce microbial and environmental barriers against opportunistic pathogens. Probiotic bacteria are among the most widely investigated tools, acting through competitive exclusion, production of antimicrobials, nutrient competition, and system-level stabilization of rearing conditions. Early frameworks positioned probiotics as biological control agents in aquaculture and highlighted their potential to reduce disease pressure while supporting host performance [12]. At a conceptual level, the effectiveness of probiotics is rarely attributable to a single mechanism; rather, it emerges from shifts in community interactions and resource partitioning that favor beneficial guilds and suppress pathogen success. Comparative sequence-based approaches have long been used to differentiate microbial pathogens across hosts and settings [13], and analogous principles underpin culture-independent profiling in aquaculture: probiotic establishment, pathogen suppression, and functional transitions are best evaluated at the level of whole communities rather than individual isolates alone.

Among candidate probiotics, Bacillus spp. are attractive because spores tolerate stressful rearing conditions and industrial handling, while vegetative cells can produce a spectrum of antimicrobial metabolites and enzymes that influence both pathogens and water quality. Reviews emphasize that reliable probiotics require systematic screening that includes robust identification, safety assessment, and mechanistic evaluation under relevant environmental conditions [14,15,16]. In shrimp systems, Bacillus subtilis supplementation has been associated with improved growth performance, enhanced digestive enzyme activity, modulation of immune gene expression, and increased resistance under disease challenge [17]. Administering Bacillus strains directly in rearing water has also been reported to enhance water quality and increase resistance against Vibrio infection, supporting the idea that probiotics can reshape the rearing environment as well as the host–microbe interface [18]. Broader syntheses further highlight the growing evidence base for Bacillus as probiotics in aquaculture and the mechanistic versatility of this genus [19]. At the molecular level, Bacillus lipopeptides represent a key class of bioactive compounds that can inhibit competitors and shape interbacterial interactions, offering a plausible functional bridge between Bacillus dominance and anti-Vibrio effects in mixed communities [20]. Complementary management strategies, such as biofloc technology, likewise rely on microbial community engineering to improve water quality and reduce pathogen pressure, reinforcing the central role of microbiome structure in disease-resilient production [21].

Despite strong rationale and extensive research, translating probiotics into consistent field performance remains challenging because efficacy depends on strain selection, dosing, formulation, and environmental context. In practice, native consortia derived from the local production environment may offer advantages in ecological compatibility and persistence, but their deployment requires rigorous validation that connects culture-based antagonism phenotypes with community-wide outcomes. Standardized in vitro approaches for antimicrobial activity assessment provide essential first-line evidence of antagonism and help triage candidates prior to in vivo challenge trials or farm-scale implementation [22]. However, culture-based assays alone cannot resolve the non-culturable fraction of hatchery microbiomes, nor can they quantify how probiotics restructure broader bacterial networks that influence pathogen success and functional potential.

High-throughput sequencing and reproducible bioinformatic pipelines now enable community-level evaluation at a resolution suitable for mechanism-oriented inference in aquaculture systems. QIIME 2 provides an extensible framework for microbiome data science with provenance tracking, supporting transparent processing and analysis [23]. DADA2 enables high-resolution inference of amplicon sequence variants (ASVs) without reliance on OTU clustering, improving sensitivity to ecological change across treatments and time [24]. Complementary preprocessing steps, including adapter and quality trimming with Cutadapt and Trimmomatic, support accurate recovery of biological signal from raw reads [25,26], while VSEARCH provides a versatile open-source toolkit for key operations such as dereplication and chimera-related procedures when needed [27]. Taxonomic assignment commonly relies on curated rRNA databases such as SILVA [28], and primer choice is a critical determinant of coverage and bias in bacterial diversity studies [29]. In R, phyloseq supports reproducible microbiome analysis and visualization by integrating feature tables, taxonomy, metadata, and derived diversity metrics in a unified structure [30]. When the analytical aim extends beyond description toward identifying taxa that differ across experimental conditions, DESeq2 offers moderated dispersion estimation and fold-change inference that is widely applied to sequencing count data [31].

In addition to taxonomic restructuring, functional interpretation is increasingly central to evaluating whether probiotic-associated community shifts plausibly translate into pathogen suppression. Tools such as PICRUSt2 can infer functional potential from amplicon data [32], and predicted pathways can be interpreted using reference resources such as KEGG [33], eggNOG [34], and the COG database [35]. Although functional inference from 16S data has limitations, integrating predicted functions with phenotypic assays (e.g., inhibition halos) provides a pragmatic framework to test whether community restructuring aligns with bioactive potential relevant to antimicrobial activity, stress response, or competitive fitness.

To move from descriptive microbiomics to mechanism-oriented hypotheses about disease control, integrating multivariate and structural modeling is particularly valuable. Phylogeny-aware beta-diversity metrics such as UniFrac support robust comparison of communities across experimental conditions [36], while ordination approaches summarize complex ecological gradients and highlight treatment-specific regimes. Principal component analysis (PCA) remains a cornerstone method for dimension reduction and visualization of multivariate patterns [38,39], with standardized workflows available through tools such as FactoMineR [40]. Beyond ordination, structural equation modeling (SEM) provides a formal framework to quantify directed relationships among latent constructs and observed indicators when community structure, environmental variation, and functional outcomes interact. Classical reliability metrics such as Cronbach’s alpha support internal consistency assessment of multi-indicator constructs [41], while convergent and discriminant validity can be evaluated using established criteria including Fornell–Larcker and the heterotrait–monotrait ratio (HTMT) [42,43]. For prediction-oriented, multi-construct models in complex biological systems, variance-based SEM using partial least squares (PLS-SEM) is widely used because it can accommodate complex models and does not require strict distributional assumptions [44,45,46]. In parallel, predictive machine-learning approaches can provide complementary validation by testing whether multivariate indices learned from community and environmental data generalize to held-out observations, strengthening confidence that identified drivers are not merely descriptive but also predictive of antagonistic outcomes.

From an ecological standpoint, pathogen control in shrimp hatcheries is best viewed as management of a dynamic microbial meta-community rather than suppression of a single taxon. Virulence plasmids and toxin repertoires explain why specific V. parahaemolyticus lineages can trigger AHPND [1,2,3,4,5], yet disease expression is shaped by background community competition, resource availability, and physicochemical stressors that modulate host susceptibility and microbial growth rates [6,7,8,9,10,11]. Accordingly, interventions that shift the community toward stable, competitive, and functionally protective states may reduce the probability that toxigenic vibrios reach critical abundances or express virulence determinants. Within this logic, Bacillus-based probiotics remain among the most practical tools for community-level engineering because spores withstand storage and delivery constraints, and Bacillus metabolism can contribute both direct antagonism and indirect system stabilization [14,15,16,17,18,19,20,21]. Nevertheless, variability across farms persists, often reflecting mismatches between probiotic strains and local conditions—an argument for prioritizing native consortia that are already adapted to the production environment.

Objective of this study. In this context, the present study aimed to evaluate native probiotic consortia (CN5 and RS3), individually and as a mixed consortium (MIX), for their capacity to suppress V. parahaemolyticus and restructure hatchery water microbiomes across a 30-day time course, relative to a no-probiotic control (CTRL). Specifically, we combined (i) standardized in vitro antagonism screening to quantify anti-Vibrio activity [22], (ii) 16S rRNA gene (V3–V4) amplicon profiling with reproducible pipelines to characterize taxonomic and diversity shifts [23,24,25,26,27,28,29,30,31], (iii) functional inference to estimate bioactive potential relevant to pathogen suppression [32,33,34,35], and (iv) integrative multivariate and structural modeling (PCA and PLS-SEM) to test mechanistic relationships linking water quality, microbial diversity, Bacillus dominance, Vibrio presence, inferred bioactive function, and antagonistic outcomes [36,38,39,40,41,42,43,44,45,46]. By integrating phenotype, community composition, inferred function, and mechanistic modeling, this work seeks to provide a high-resolution, field-relevant evaluation of native probiotic consortia as scalable biocontrol tools for disease-resilient shrimp hatchery management.

2. Materials and Methods

2.1. Experimental Design

An exploratory and analytical experimental design was implemented to (i) characterize bacterial community dynamics in Penaeus vannamei larval culture water and (ii) quantify the antagonistic potential of native probiotic consortia against Vibrio parahaemolyticus, integrating culture-based screening, 16S rRNA gene sequencing (V3–V4), and multivariate/causal modeling. The target pathogen was selected due to its opportunistic behavior in shrimp aquaculture and its association with AHPND through acquisition of virulence plasmids and Pir-like toxins [1,2,3,4,5]. Additional genomic and ecological evidence supporting AHPND/Vibrio epidemiology and strain diversification was considered when defining the analytical scope [6,7,8,9,10,11].

Water was sampled from a larval reservoir located in Guayas Province (Ecuador) across four sampling moments separated by 10-day intervals (Day 0, 10, 20, and 30). At each time point, four experimental conditions were evaluated: two native probiotic consortia (CN5, RS3), a mixed consortium (MIX = CN5+RS3), and a no-probiotic control (CTRL). Each treatment–time combination included eight independent biological replicates (4 treatments × 4 times × 8 replicates), yielding a total of 128 observations for SEM-aligned indices and downstream multivariate analyses (Table 1).

2.2. Evaluated Variables

Six core response domains were evaluated and harmonized as SEM-aligned indices for integrative analysis shown in the following list and Table 2:

• Anti-Vibrio activity: inhibition halo diameter (mm) from a simultaneous inhibition/competitive exclusion assay; summarized as mean halo per observation. Antimicrobial evaluation procedures followed standardized in vitro guidance [22].
• Vibrio presence: index derived from taxonomic profiles (relative abundance of Vibrio spp., including V. parahaemolyticus when detected), consistent with known Vibrio biodiversity and disease relevance [10,11].
• Bacillus dominance: index derived from relative abundance of Bacillus spp., reflecting the expected probiotic-enriched regime based on Bacillus probiotic use in shrimp [17,18,19].
• Microbial diversity: alpha-diversity index computed from ASV/feature tables (e.g., observed genera and Simpson diversity), supporting ecological interpretation of community restructuring.
• Bioactive function: index based on inferred functional potential using PICRUSt2 [32], interpreted with KEGG pathway hierarchies [33] and complementary orthology resources (eggNOG, COG) [34,35].
• Water quality: composite index derived from dissolved oxygen, salinity, temperature, and pH after z-standardization.

2.3. Biological Material

The biological material comprised:

• culture water from penaeid shrimp larval production systems (matrix for microbiome profiling and probiotic screening);
• two native probiotic consortia (CN5 and RS3) assembled from Bacillus-enriched isolates obtained from the same production environment, in line with probiotic selection principles in aquaculture [12,14,15,16];
• a mixed consortium treatment (MIX) prepared by combining CN5 and RS3 at equal proportions; and
• a no-probiotic control (CTRL).

Vibrio parahaemolyticus served as the target pathogen for antagonism assays because AHPND-linked virulence is associated with toxin-encoding plasmids (Pir-like toxins) and related genomic features [1,2,3,4,5,8,9].

2.4. Culture-Based Isolation and Morphotypic Screening

Culture-based isolation was performed to recover native candidate probiotic bacteria from larval reservoir water. Samples were serially diluted (10⁻¹–10⁻²) in sterile saline and plated on tryptic soy agar (TSA) and selective Chromagar™ Bacillus to enrich for Bacillus-like morphotypes. Plates were incubated aerobically until visible colony development. Distinct colonies were selected based on morphology (size, texture, pigmentation, margin, elevation), re-streaked to purity, and preserved (TSA slants and/or glycerol stocks at −80 °C).

Two candidate consortia (CN5 and RS3) were defined as distinct pools of Bacillus-enriched isolates that showed consistent morphotypic profiles and inhibitory activity in preliminary screens, consistent with recommended screening steps for probiotic selection (Table 3) [14,15,16]. Mechanistically, Bacillus-associated antagonism was supported conceptually by lipopeptide-mediated inhibition and competitive interactions [20].

2.5. In Vitro Antagonism Assay and Antagonistic Effectiveness

Antagonistic activity against V. parahaemolyticus was quantified using a simultaneous inhibition (competitive exclusion) assay following standardized antimicrobial evaluation approaches [22]. Briefly, a V. parahaemolyticus lawn was prepared on TSA plates, after which consortia (CN5, RS3, MIX) were inoculated at predefined positions. Plates were incubated under aerobic conditions and inhibition halos were measured (mm) along orthogonal axes and averaged per plate.

Antagonistic effectiveness (AE, %) was computed as the proportion of observed inhibition relative to the maximum measurable radius permitted by the assay geometry. Each treatment–time group included eight replicate plates to support robust estimation of mean halo ± SE and subsequent multivariate modeling.

2.6. Data Analysis

2.6.1. DNA Extraction, Library Preparation, and Sequencing

For microbiome profiling, each water sample replicate was homogenized and an aliquot (250 mL) was filtered through a 0.45 µm membrane using vacuum filtration. DNA was extracted from membranes using DNeasy® PowerWater® (QIAGEN) following the manufacturer’s protocol. DNA quality/quantity were evaluated via NanoDrop™ and Qubit® 4.0, and integrity was checked by 1% agarose gel electrophoresis. The V3–V4 region of the 16S rRNA gene was amplified using universal primers 341F/805R, consistent with common primer evaluation for bacterial diversity studies [29]. Amplicons were purified, indexed, pooled, and sequenced on an Illumina MiSeq platform (2 × 300 bp).

2.6.2. Bioinformatic Processing and Taxonomic Assignment

Raw reads were processed in QIIME2 using DADA2 for denoising, paired-end merging, and chimera removal to generate ASVs [23,24]. Adapter removal and quality trimming followed standard preprocessing tools (Cutadapt and Trimmomatic) [25,26]. Where needed, VSEARCH supported dereplication and related operations [27]. Taxonomy was assigned using a Naïve Bayes classifier trained on SILVA (release 138) for the V3–V4 region [28]. Outputs (ASV table, taxonomy, metadata) were imported into R for downstream ecological processing with phyloseq [30]. Differential abundance testing across treatments and time points used DESeq2 with multiple-testing correction (Table 4) [31].

2.6.3. Functional Inference and Annotation

Functional potential was inferred from ASVs using PICRUSt2 [32]. Predicted pathways were summarized against KEGG hierarchies [33] and complemented with orthology resources (eggNOG and COG) to support functional interpretation [34,35]. A “bioactive function” index was constructed by aggregating inferred functions plausibly related to antimicrobial biosynthesis, stress response, and competitive fitness.

2.6.4. Multivariate Statistics (PCA)

Beta-diversity was computed using UniFrac distances [36]. For integrative visualization, PCA was performed on centered and scaled SEM-aligned indices using standard PCA definitions [38,39] and implemented in R via FactoMineR [40]. Cluster structure in PCA space was explored using k-means on leading PCs to summarize multivariate regimes.

2.6.5. Machine Learning (Random Forest) for Out-of-Sample Validation

To quantify predictive generalization beyond in-sample modeling, a Random Forest regressor was trained to predict anti-Vibrio activity (halo mean, mm) from SEM-aligned indices (e.g., Bacillus dominance, Vibrio presence, microbial diversity, bioactive function, and water quality). Data were split into training and test partitions (e.g., 70/30), preserving treatment and time representation. Hyperparameters (number of trees, mtry/max_features, node size/min_samples_leaf) were tuned via cross-validation on the training set. Predictive performance on the test set was summarized using RMSE, MAE, and R². Model interpretation used permutation importance to rank predictors by their contribution to predictive accuracy (Table 5).

2.6.6. PLS-SEM Specification and Quality Assessment

A variance-based SEM (PLS-SEM) was specified to test a mechanistic cascade in which water quality and microbial diversity influence Bacillus dominance and Vibrio presence, which in turn shape inferred bioactive function and anti-Vibrio activity. Latent constructs were modeled reflectively with observed indicators (e.g., microbial diversity: observed genera and Simpson; water quality: dissolved oxygen, salinity, temperature, pH; anti-Vibrio activity: halo measurements; bioactive function: inferred KEGG/functional indicators). Reliability was assessed using Cronbach’s alpha [41]. Convergent validity used AVE and discriminant validity used Fornell–Larcker and HTMT [42,43]. Structural paths were estimated by non-parametric bootstrapping (e.g., 5000 resamples) and reported with R², effect sizes (f²), and predictive relevance measures, following PLS-SEM guidance [44,45,46].

2.7. Null and Working Hypotheses

• H0-1 (treatment effect): No differences exist among treatments (CN5, RS3, MIX, CTRL) in anti-Vibrio activity (mean inhibition halo) at any time point (0, 10, 20, 30).H1-1: At least one probiotic treatment (CN5, RS3, MIX) increases anti-Vibrio activity relative to CTRL, potentially varying with time.
• H0-2 (community restructuring): Treatments do not alter Vibrio presence, Bacillus dominance, microbial diversity, bioactive function, or the water-quality index across time.H1-2: Probiotic treatments reduce Vibrio presence and increase Bacillus dominance, with concomitant shifts in diversity, inferred functions, and/or water quality.
• H0-3 (PCA separation): Multivariate profiles (PCA scores) do not differ among treatments.H1-3: Treatments yield reproducible multivariate separation in PCA space consistent with distinct ecological regimes.
• H0-4 (SEM paths): All PLS-SEM structural path coefficients equal zero.H1-4: At least one SEM path is non-zero, consistent with a cascade where water quality/diversity affect Bacillus/Vibrio, which mediate bioactive function and antagonism.
• H0-5 (Random Forest predictability): SEM-aligned indices do not predict anti-Vibrio activity better than chance in held-out data.H1-5: A Random Forest model predicts anti-Vibrio activity with high out-of-sample performance, and predictor importance highlights Bacillus dominance and Vibrio presence as primary drivers.

3. Results

3.1. Isolation, Morphotypic Screening, and In Vitro Antagonism of Native Probiotic Consortia

Cultivable bacteria were recovered from seawater collected from a larval reservoir during four consecutive months of Penaeus vannamei production, allowing the isolation of heterotrophic bacteria and presumptive Bacillus spp. under routine culture conditions. Viable growth was consistently obtained within the 10⁻¹–10⁻² dilution range for total heterotrophs and Bacillus spp., supporting the presence of an active cultivable fraction suitable for downstream functional screening.

Differential growth on TSA and Chromagar™ Bacillus enabled a first discrimination between general heterotrophic morphotypes and spore-forming bacilli, respectively (Figure 1B–D).

From this pool, two native consortia were selected as probiotic candidates (CN5 and RS3), while pathogenic controls for functional assays were defined as Vibrio parahaemolyticus (V.p) and an invasive colony phenotype (C.INV), ensuring traceability across assays.

A preliminary morphotypic evaluation of colonies obtained from CN5 and RS3 showed round, creamy colonies with variable size, compatible with morphotypes of probiotic interest, and these isolates/consortia were preserved for subsequent molecular confirmation (16S rRNA).

Notably, Chromagar™ Bacillus plates revealed visually distinctive colony profiles (e.g., pigmented/blue-toned colonies vs. cream colonies), consistent with a bacilli-enriched community within the selected consortia (Figure 1B).

The antagonistic potential of the selected consortia was then evaluated using a simultaneous inhibition (competitive exclusion) approach (Figure 1A,E), where inhibition zones around inoculation points indicated suppression of pathogen growth in vitro.

In agreement with the study’s overall outcome, inhibition assays showed broad and consistent halos for RS3, CN5, and their combination, reaching ~99% antagonistic effectiveness against V. parahaemolyticus and supporting their candidacy as native biocontrol agents for sustainable aquaculture.

3.2. Treatment Effects and Temporal Dynamics of Antagonism and SEM-Aligned Indices

Across the full dataset (n = 128; 32 observations per treatment), probiotic treatments (CN5, MIX, RS3) maintained consistently high antagonistic activity (mean inhibition halo), whereas the control exhibited minimal inhibition (Table 1). In parallel, the control showed higher Vibrio presence and higher microbial diversity than probiotic treatments, while Bacillus dominance was substantially enriched under CN5/MIX/RS3 (Table 1). Water quality (z-index) remained centered near zero across treatments, suggesting that the main separation between experimental conditions was associated with microbiological rather than physicochemical variation.

Table 1. Treatment-level summary of key variables (mean ± SD; n = 32 per treatment).

Variable	CN5	CTRL	MIX	RS3
Antagonism (mean halo, mm)	20.673 ± 2.208	2.266 ± 1.347	21.068 ± 1.930	20.254 ± 1.990
Vibrio presence (index)	0.0166 ± 0.0043	0.1656 ± 0.0254	0.0142 ± 0.0041	0.0116 ± 0.0025
Bacillus dominance (index)	0.3957 ± 0.0161	0.0622 ± 0.0210	0.4308 ± 0.0060	0.4565 ± 0.0097
Microbial diversity (index)	49.974 ± 4.822	97.182 ± 5.450	47.194 ± 2.157	48.720 ± 9.282
Bioactive function (index)	7.561 ± 1.075	5.519 ± 1.291	4.468 ± 0.439	3.518 ± 0.244
Water quality (z-index)	0.001 ± 0.432	−0.008 ± 0.504	−0.001 ± 0.462	0.008 ± 0.524

Table 1. Treatment-level summary of key variables (mean ± SD; n = 32 per treatment).

Variable	CN5	CTRL	MIX	RS3
Antagonism (mean halo, mm)	20.673 ± 2.208	2.266 ± 1.347	21.068 ± 1.930	20.254 ± 1.990
Vibrio presence (index)	0.0166 ± 0.0043	0.1656 ± 0.0254	0.0142 ± 0.0041	0.0116 ± 0.0025
Bacillus dominance (index)	0.3957 ± 0.0161	0.0622 ± 0.0210	0.4308 ± 0.0060	0.4565 ± 0.0097
Microbial diversity (index)	49.974 ± 4.822	97.182 ± 5.450	47.194 ± 2.157	48.720 ± 9.282
Bioactive function (index)	7.561 ± 1.075	5.519 ± 1.291	4.468 ± 0.439	3.518 ± 0.244
Water quality (z-index)	0.001 ± 0.432	−0.008 ± 0.504	−0.001 ± 0.462	0.008 ± 0.524

Temporal profiling (0, 10, 20, and 30 days) confirmed that probiotic-associated antagonism remained stable over time (≈20–22 mm), while the control remained low (≈0–6 mm), with clear separation of the control trajectories from probiotic treatments (Figure 1). Likewise, Vibrio presence remained consistently higher in CTRL, whereas Bacillus dominance remained consistently higher under CN5/MIX/RS3 (Figure 1), supporting treatment-driven stabilization of a Bacillus-dominant regime associated with antagonistic activity.

3.3. Species-Level Community Composition Across Treatments and Time

Species-level metagenomic profiling revealed a pronounced treatment-dependent restructuring of the microbial community (Figure 2 and Figure 3). When visualized at the level of individual observations (n = 128), stacked-bar profiles showed that CTRL samples were characterized by relatively higher contributions of Vibrio spp. (notably V. alginolyticus, V. parahaemolyticus, V. jasicida, and V. xuii), whereas probiotic treatments were consistently dominated by Bacillus spp., with recurrent prominence of B. licheniformis and B. amyloliquefaciens (Figure 2). Importantly, “Unclassified_species” was retained as an explicit component of the composition, ensuring transparency in species-level assignment.

To improve interpretability at the group level, mean composition was summarized for each treatment–time combination (16 groups; n = 8 replicates per group) (Figure 3). This aggregation confirmed that CTRL maintained a Vibrio-enriched profile throughout the time course, while CN5/MIX/RS3 maintained a stable Bacillus-enriched profile. Across probiotic treatments, Bacillus dominance was sustained over time, consistent with the functional phenotype of increased antagonism against Vibrio observed in vitro.

3.4. Multivariate Differentiation by PCA and Unsupervised k-Means Clustering

PCA performed on standardized SEM-aligned indices (antagonism, Vibrio presence, Bacillus dominance, microbial diversity, bioactive function, and water quality) explained 81.8% of the total variance in the first two components (PC1 = 64.7%, PC2 = 17.1%) (Figure 4). PC1 primarily contrasted Vibrio presence and microbial diversity against Bacillus dominance and antagonism, providing a compact multivariate summary that is consistent with the species-level restructuring observed by metagenomics (Figure 2 and Figure 3).

Unsupervised k-means clustering (k = 4) in PCA space supported robust group structure. CTRL samples were assigned to distinct clusters, whereas probiotic treatments concentrated into two major multivariate regimes (Table 2), consistent with probiotic-driven ecological states.

Figure 5. PCA biplot of SEM-aligned indices with k-means clustering (k = 4). Points represent observations; shapes denote treatment groups; convex hulls summarize cluster geometry. Arrows indicate variable loadings (direction and contribution) on PC1 and PC2. Dashed lines denote axes at zero.

Figure 5. PCA biplot of SEM-aligned indices with k-means clustering (k = 4). Points represent observations; shapes denote treatment groups; convex hulls summarize cluster geometry. Arrows indicate variable loadings (direction and contribution) on PC1 and PC2. Dashed lines denote axes at zero.

Table 2. k-Means cluster membership by treatment (counts and within-treatment %; k = 4).

Treatment	Cluster 1	Cluster 2	Cluster 3	Cluster 4
CN5 (n=32)	0 (0.0%)	18 (56.2%)	14 (43.8%)	0 (0.0%)
CTRL (n=32)	24 (75.0%)	0 (0.0%)	0 (0.0%)	8 (25.0%)
MIX (n=32)	0 (0.0%)	10 (31.2%)	22 (68.8%)	0 (0.0%)
RS3 (n=32)	0 (0.0%)	7 (21.9%)	25 (78.1%)	0 (0.0%)

Table 2. k-Means cluster membership by treatment (counts and within-treatment %; k = 4).

Treatment	Cluster 1	Cluster 2	Cluster 3	Cluster 4
CN5 (n=32)	0 (0.0%)	18 (56.2%)	14 (43.8%)	0 (0.0%)
CTRL (n=32)	24 (75.0%)	0 (0.0%)	0 (0.0%)	8 (25.0%)
MIX (n=32)	0 (0.0%)	10 (31.2%)	22 (68.8%)	0 (0.0%)
RS3 (n=32)	0 (0.0%)	7 (21.9%)	25 (78.1%)	0 (0.0%)

3.5. PLS-SEM Structural Model: Explained Variance, Effect Sizes, and Direct Effects

The PLS-SEM exhibited high explanatory power for endogenous constructs (Table 3), with R² values of 0.875 (Bacillus dominance), 0.911 (Vibrio presence), 0.946 (bioactive function), and 0.943 (anti-Vibrio activity). Effect-size assessment (f²) indicated very large contributions of microbial diversity to Bacillus dominance (f² = 6.998) and of Bacillus dominance to bioactive function (f² = 17.431), consistent with a cascade where community structure governs functional outputs and antagonism (Table 3).

Table 3. Structural model quality: explained variance (R², adjusted R²), effect sizes (f²), and predictive relevance (Q²_predict).

Endogenous construct	R²	Adjusted R²	Key f² contributors (toward the endogenous construct)	RMSE	MAE	Q²_predict
Bacillus dominance	0.875	0.873	Microbial diversity → Bacillus (6.998); Water quality → Bacillus (0.300)	0.383	0.288	0.858
Vibrio presence	0.911	0.909	Microbial diversity → Vibrio (0.225); Water quality → Vibrio (0.208)	0.453	0.365	0.801
Anti-Vibrio activity	0.943	0.942	Bacillus → Activity (0.563); Vibrio → Activity (0.026)	0.404	0.329	0.841
Bioactive function	0.946	0.945	Bacillus → Function (17.431); Vibrio → Function (0.947)	0.470	0.394	0.785

Table 3. Structural model quality: explained variance (R², adjusted R²), effect sizes (f²), and predictive relevance (Q²_predict).

Endogenous construct	R²	Adjusted R²	Key f² contributors (toward the endogenous construct)	RMSE	MAE	Q²_predict
Bacillus dominance	0.875	0.873	Microbial diversity → Bacillus (6.998); Water quality → Bacillus (0.300)	0.383	0.288	0.858
Vibrio presence	0.911	0.909	Microbial diversity → Vibrio (0.225); Water quality → Vibrio (0.208)	0.453	0.365	0.801
Anti-Vibrio activity	0.943	0.942	Bacillus → Activity (0.563); Vibrio → Activity (0.026)	0.404	0.329	0.841
Bioactive function	0.946	0.945	Bacillus → Function (17.431); Vibrio → Function (0.947)	0.470	0.394	0.785

Bootstrapping results (Table 4) indicated a strong positive effect of Bacillus dominance on anti-Vibrio activity (β = 0.803; t = 8.967; p < 0.001) and on bioactive function (β = 0.972; t = 188.864; p < 0.001). Microbial diversity strongly decreased Bacillus dominance (β = −0.947; t = 39.715; p < 0.001) and increased Vibrio presence (β = 0.323; t = 3.830; p < 0.001). Bioactive function showed a strong negative relationship with Vibrio presence (β = −0.653; t = 7.845; p < 0.001). Water-quality paths were not statistically supported at α = 0.05 (p ≥ 0.065), while the direct path from Vibrio presence to anti-Vibrio activity was marginal (β = −0.172; p = 0.057), suggesting that antagonism was primarily driven by the Bacillus-dominant state rather than by water chemistry alone.

Table 4. Bootstrapping results for direct effects in the structural model.

Path	β (Original sample)	Mean (bootstrap)	STDEV	t	p
Bacillus dominance → Anti-Vibrio activity	0.803	0.825	0.090	8.967	<0.001
Bacillus dominance → Bioactive function	0.972	0.972	0.005	188.864	<0.001
Vibrio presence → Anti-Vibrio activity	−0.172	−0.149	0.090	1.901	0.057
Water quality → Bacillus dominance	−0.196	−0.122	0.106	1.842	0.065
Water quality → Vibrio presence	0.146	0.113	0.096	1.521	0.128
Bioactive function → Vibrio presence	−0.653	−0.689	0.083	7.845	<0.001
Microbial diversity → Bacillus dominance	−0.947	−0.934	0.024	39.715	<0.001
Microbial diversity → Vibrio presence	0.323	0.287	0.084	3.830	<0.001

Table 4. Bootstrapping results for direct effects in the structural model.

Path	β (Original sample)	Mean (bootstrap)	STDEV	t	p
Bacillus dominance → Anti-Vibrio activity	0.803	0.825	0.090	8.967	<0.001
Bacillus dominance → Bioactive function	0.972	0.972	0.005	188.864	<0.001
Vibrio presence → Anti-Vibrio activity	−0.172	−0.149	0.090	1.901	0.057
Water quality → Bacillus dominance	−0.196	−0.122	0.106	1.842	0.065
Water quality → Vibrio presence	0.146	0.113	0.096	1.521	0.128
Bioactive function → Vibrio presence	−0.653	−0.689	0.083	7.845	<0.001
Microbial diversity → Bacillus dominance	−0.947	−0.934	0.024	39.715	<0.001
Microbial diversity → Vibrio presence	0.323	0.287	0.084	3.830	<0.001

3.6. Out-of-Sample Validation: Random Forest Prediction of Antagonism

A Random Forest regressor trained on SEM-aligned indices achieved strong predictive performance on the held-out test set (RMSE = 2.194; MAE = 1.87; R² = 0.93) (Figure 6). Permutation importance ranked Bacillus dominance as the most influential predictor of antagonism, followed by Vibrio presence and microbial diversity, whereas water quality contributed minimally (Figure 7). This ranking converges with the SEM results, reinforcing Bacillus dominance as the primary driver of the antagonistic phenotype.

Figure 8. Random Forest permutation importance for antagonism prediction. Higher values indicate a stronger decrease in predictive performance when the corresponding predictor is permuted.

Figure 8. Random Forest permutation importance for antagonism prediction. Higher values indicate a stronger decrease in predictive performance when the corresponding predictor is permuted.

3.7. Model Diagnostics: Global Fit and Collinearity

Global fit indices indicated SRMR = 0.149 for both saturated and estimated models, with discrepancy indices reported in Table 5. Collinearity diagnostics identified elevated VIF values for Temp_C (5.642) and pH (6.168), whereas most indicators exhibited acceptable collinearity levels (VIF ≈ 1–4) (Table 5). These diagnostics suggest redundancy in the water-quality block and should be considered when interpreting water-quality effects and when refining model specification.

Table 5. Global fit and collinearity diagnostics (VIF). (A) Global fit indices

Index		Saturated model	Estimated model
SRMR		0.149	0.149
${\mathit{D}}_{\mathit{U}\mathit{L}\mathit{S}}$		3.381	3.394
${\mathit{D}}_{\mathit{G}}$		2.682	3.049
χ2		1057.135	1157.327
NFI		0.681	0.651
(B) Variance inflation factors (VIF) for indicators
Indicator	VIF	Indicator	VIF
DissolvedO2_mgL	1.394	Salinity_ppt	1.273
Genera_observed	1.103	Simpson	1.103
Halo_CINV_CHROM	3.893	Species_level_%	1.116
Halo_CINV_TSA	4.206	Temp_C	5.642
Halo_Vp_CHROM	3.760	pH	6.168
Halo_Vp_TSA	3.961	rel_B_amyloliquefaciens	1.149
KEGG_AntimicrobialScore_meta	1.116	rel_B_licheniformis	1.349
rel_B_sonorensis	1.319	rel_V_alginolyticus	3.217
rel_V_xuii	3.217

Table 5. Global fit and collinearity diagnostics (VIF). (A) Global fit indices

Index		Saturated model	Estimated model
SRMR		0.149	0.149
${\mathit{D}}_{\mathit{U}\mathit{L}\mathit{S}}$		3.381	3.394
${\mathit{D}}_{\mathit{G}}$		2.682	3.049
χ2		1057.135	1157.327
NFI		0.681	0.651
(B) Variance inflation factors (VIF) for indicators
Indicator	VIF	Indicator	VIF
DissolvedO2_mgL	1.394	Salinity_ppt	1.273
Genera_observed	1.103	Simpson	1.103
Halo_CINV_CHROM	3.893	Species_level_%	1.116
Halo_CINV_TSA	4.206	Temp_C	5.642
Halo_Vp_CHROM	3.760	pH	6.168
Halo_Vp_TSA	3.961	rel_B_amyloliquefaciens	1.149
KEGG_AntimicrobialScore_meta	1.116	rel_B_licheniformis	1.349
rel_B_sonorensis	1.319	rel_V_alginolyticus	3.217
rel_V_xuii	3.217

4. Discussion

Acute hepatopancreatic necrosis disease (AHPND) remains one of the most disruptive syndromes in shrimp aquaculture because it can trigger rapid mortality while simultaneously reshaping microbial equilibria in culture systems, making prevention strategies strongly dependent on ecological control rather than single-pathogen suppression [47]. Recent syntheses emphasize that effective mitigation requires combining pathogen monitoring with interventions that reduce the probability of toxigenic Vibrio dominance under fluctuating physicochemical conditions typical of hatchery and nursery operations [48]. Under this framework, the present study is best interpreted as an integrated evaluation of (i) the antagonistic phenotype of native probiotic consortia (CN5, RS3, and MIX) against Vibrio parahaemolyticus and (ii) the extent to which those consortia are associated with microbiome restructuring consistent with protective ecological states.

4.1. Relevance of the Findings for AHPND-Risk Management in the Americas

Phylogenomic evidence supports that AHPND-associated V. parahaemolyticus lineages in Latin America may arise via multiple introduction pathways and/or local evolutionary trajectories, reinforcing the need for regionally validated prevention measures rather than assuming uniform strain ecology across continents [49]. Field-based detection work in the Americas further demonstrated that AHPND can become established in production landscapes where environmental and operational drivers (e.g., water exchange, organic loading, and temperature variability) create recurrent opportunities for toxigenic Vibrio expansion [50]. In this context, the high inhibitory performance observed for CN5/RS3/MIX in the competitive exclusion assay (Figure 3.1; inhibition halos clearly visible across plates) is not merely a laboratory phenotype—it represents a practical screening signal for identifying probiotic candidates capable of counteracting V. parahaemolyticus under conditions that often precede AHPND expression in farms.

A complementary concern is that AHPND systems may act as reservoirs for mobile genetic elements that carry not only virulence determinants but also antimicrobial resistance, which can be co-selected under farm pressures and complicate treatment options [51]. This strengthens the rationale for non-antibiotic strategies: your results support a pathway in which native consortia—selected from the same production environment—can generate strong in vitro suppression while also aligning with community patterns consistent with reduced Vibrio presence.

4.2. Water-Quality Context and “System State” Interpretation

Biofloc-based and related microbial-management approaches highlight that disease risk is often mediated by a system-level balance among C/N regime, microbial assimilation of nitrogenous wastes, and community competition, rather than by the pathogen alone [52]. Although your work is not framed as a full biofloc trial, the observed trajectories in the water-quality composite index across sampling moments can be interpreted as part of a broader “system maturation” process in which physicochemical stabilization may create a narrower niche window for opportunistic Vibrio proliferation. Importantly, your SEM structure (water quality → community structure → functional potential → anti-Vibrio activity) is consistent with the idea that environmental conditions shape microbial assembly, which then determines protective capacity.

4.3. Consortia Performance: From Inhibition Halos to Community Restructuring

Evidence from AHPND-related intervention studies shows that probiotic administration can improve survival and reshape bacterial communities when applied after or during pathogen pressure, but outcomes depend on strain ecology and environmental compatibility [53]. Your findings (higher inhibition halos under probiotic treatments vs. control; Figure 3.1) align with the general expectation that multi-strain or community-based interventions can outperform single-mechanism approaches when competitive dynamics in the water column are central.

Consortium-based biocontrol has also been demonstrated in “hybrid” designs that combine microalgae and bacteria, where the protective effect emerges through both direct antagonism and resource/oxygen modulation [54]. Similarly, non-Bacillus probiotics (e.g., Pseudoalteromonas) have been reported to increase resistance of P. vannamei to AHPND-causing Vibrio, indicating that multiple taxonomic routes can reach a protective functional outcome [55]. In your case, however, the sequencing-aligned indices suggest that protection is strongly coupled to Bacillus dominance, which is mechanistically plausible given this genus’ capacity for persistent colonization and metabolite production.

The MIX treatment conceptually resembles “synbiotic-like” logic reported for combined functional additives, where synergy between components can broaden antimicrobial and ecological effects [56]. Interpreted through your PCA separation and SEM paths, MIX can be discussed not only as “CN5 + RS3” but as a potential mechanism-complementation strategy, where distinct isolate pools contribute overlapping but not identical ecological functions.

4.4. AHPND as a Dysbiosis Trigger and the Meaning of Diversity Shifts

AHPND is increasingly recognized as a dysbiosis-associated process: disease states can correlate with destabilized community composition and functional signatures that favor opportunists, particularly vibrios, under stressful or nutrient-rich conditions [57]. This is critical for interpreting your diversity results. If probiotic treatments reduce alpha-diversity while increasing Bacillus dominance, this does not necessarily indicate “worse” ecology; rather, it can represent a protective domination state where competitive exclusion constrains pathogen expansion (especially if water quality remains stable and antagonistic potential increases).

Experimental infection work further indicates that shrimp-associated microbiomes can respond differently to pathogenic vs. non-pathogenic V. parahaemolyticus exposures, supporting the idea that not all Vibrio signals are equal and that community response patterns can be diagnostic of risk states [58]. In larval contexts specifically, microbiome composition can shift markedly during disease, and community-based biomarkers have been proposed to distinguish healthy vs. compromised larval systems [59]. Your design—combining culture-based antagonism with water-microbiome profiling—fits directly within this diagnostic-preventive paradigm.

Beyond Bacillus, actinomycete-based probiotics (e.g., Streptomyces formulations) have shown capacity to modulate shrimp gut microbiota and improve performance metrics, reinforcing that functional outcomes can be achieved through different microbial guilds and metabolite suites [60]. This supports a key interpretive point for your discussion: the central question is not whether the system becomes more diverse, but whether it becomes more functionally defensive and less permissive to Vibrio proliferation.

4.5. Mechanistic Plausibility: Bacillus-Driven Protection, Quorum Quenching, and Bioactive Metabolites

A closely comparable line of evidence shows that dietary supplementation with Bacillus velezensis can modulate shrimp microbiota and enhance resistance under AHPND-relevant contexts, supporting the plausibility of your Bacillus-dominance → protection pathway [61]. Mechanistically, quorum quenching and interference with Vibrio signaling has emerged as a credible route by which Bacillus strains reduce virulence expression and colonization efficiency, in addition to direct growth inhibition [62]. Inhibition of quorum-sensing–regulated behaviors (e.g., biofilm formation, motility, and coordinated virulence) is especially relevant because AHPND risk is linked not only to abundance but also to expression dynamics in the pathogen.

Recent evidence of quorum-sensing inhibition by B. velezensis against shrimp-pathogenic Vibrio spp. supports interpreting your inhibition halos as a composite of (i) antibacterial metabolites and (ii) anti-virulence interference that reduces competitive success [63]. This interpretation is consistent with the broader strategy of disrupting quorum sensing as a disease-control approach in aquaculture systems [64], and with the policy-level direction toward alternatives to antibiotics that reduce selection pressure for resistance while maintaining productivity [65]. Reviews focused on vibrios specifically emphasize that quorum-sensing interference can target key behaviors central to pathogenic success, making it an attractive complement to community engineering [66].

At the metabolite level, Bacillus lipopeptides and biosurfactants are well-established as multifunctional compounds that can inhibit competitors and influence microbial surface interactions—mechanisms that align strongly with the visible inhibition halos and with a SEM pathway linking Bacillus dominance to inferred bioactive functional potential [67]. In your discussion, this supports framing Bacillus dominance not as a taxonomic endpoint, but as a functional driver capable of reshaping both competitive outcomes and predicted pathway signatures.

4.6. Interpreting Inferred Functions Responsibly

Because your “bioactive function index” is inferred from amplicon profiles, it should be discussed with appropriate caution. Predictive functional profiling (e.g., PICRUSt-style inference) can generate useful hypotheses about pathway shifts when shotgun metagenomics is not available, and it has been widely used to bridge composition-to-function narratives [68]. However, functional redundancy and taxonomic-function decoupling can limit the resolution of inference: distinct communities may encode overlapping functional capacities, and function can vary at strain level within a genus [69]. Therefore, your discussion should position inferred functions as directional evidence consistent with antagonism phenotypes, not as definitive proof of metabolite production.

A second interpretive safeguard is the compositional nature of relative abundance data: changes in one taxon necessarily affect the proportional representation of others, which can amplify apparent shifts unless interpreted using appropriate normalization and modeling logic [70]. This is relevant for explaining why Bacillus dominance may coincide with decreased alpha-diversity and reduced relative abundance of other groups, without implying absolute elimination.

4.7. Integrated Analytics: PCA + Random Forest + PLS-SEM as Convergent Evidence

Your multivariate results can be framed as convergent evidence that treatments correspond to distinct ecological regimes. Random forests provide a robust, nonlinear framework for classification/regression that can capture interactions among predictors—important in microbiome datasets where effects are rarely purely additive [71]. In microbial ecology, random forest approaches have been successfully used to identify discriminant taxa or features and to support predictive interpretation of community shifts [72]. At the same time, best-practice frameworks emphasize that microbiome machine learning must be interpreted carefully (avoiding leakage, ensuring cross-validation discipline, and treating feature importance as suggestive rather than causal) [73]. In your case, Random Forest can be positioned as a confirmatory ranking tool that highlights which SEM-aligned indices (e.g., Bacillus dominance, Vibrio presence, water quality, diversity) most strongly discriminate treatments/time points.

Because implementation details can influence reproducibility and performance, noting the use of a well-established random forest implementation for high-dimensional data is defensible this context [74]. The central narrative then becomes:

• PCA: demonstrates separation of treatment regimes in reduced dimensional space (ecological-state visualization).
• Random Forest: ranks which indices best predict regime membership (predictive corroboration).
• PLS-SEM: tests an explicit directed mechanism (pathway-based explanation).

To strengthen the SEM interpretation, you can cite the importance of consistent estimation logic and reliability in variance-based SEM contexts [75], while also acknowledging that predictive performance and out-of-sample relevance are increasingly recognized as core evaluation dimensions for PLS-oriented models [76]. This allows you to argue that your model is not only statistically coherent but also aligned with prediction-oriented objectives typical of applied aquaculture microbiome management.

4.8. Limitations and Implications for Application

Two limitations should be stated clearly and framed constructively. First, functional inference should be treated as hypothesis-generating and ideally validated in follow-up work using targeted metabolomics (e.g., lipopeptides), qPCR of functional genes, or shotgun metagenomics, consistent with recognized limits of inference and redundancy [68,69]. Second, the inhibition-halo assay is a controlled approximation of competitive outcomes; translating effects to operational settings requires validation under farm-like complexity and dosing regimes, consistent with the broader experience that probiotic efficacy is context-dependent [48,53].

5. Conclusions

This study demonstrates that native Bacillus-based probiotic consortia (CN5, RS3, and their combination) are consistently associated with strong in vitro antagonistic activity against Vibrio parahaemolyticus and with pronounced restructuring of shrimp hatchery water microbiomes toward Bacillus-dominated states. By integrating culture-based inhibition assays with amplicon sequencing, functional inference, and multivariate and structural modeling, the results support a coherent ecological pattern in which Bacillus dominance and reduced Vibrio presence co-occur with increased antagonistic potential. Importantly, the findings do not establish direct causality between specific metabolites, gene products, or in vivo disease outcomes, nor do they demonstrate protection against AHPND under farm or challenge conditions. Functional predictions were inferred from 16S rRNA gene data and therefore represent potential rather than experimentally confirmed metabolic activity. Water-quality effects were limited within the modeled system and should be interpreted cautiously. Consequently, while the data support the suitability of native probiotic consortia as promising biocontrol candidates and as drivers of protective microbial community states, further validation under controlled infection trials and commercial-scale conditions is required. Overall, this work contributes mechanistic and ecological evidence that community-level microbial management, rather than single-pathogen targeting, is a viable framework for improving disease resilience in shrimp hatchery environments.