Big-belly seahorses were maintained and treated in accordance with the guidelines of Animal Ethics Experimentation approved by the Animal Care and Use Committee of Ludong University (document number: LDU-RB20210308NXY-9). Healthy big-belly seahorses were collected from the Wendeng Seahorse Center of Ludong University, Yantai 264025, Shandong Province, China. Male and female seahorses were maintained in ponds connected to a central circulation system with mechanical and biological filtration, ultraviolet sterilization, and a protein skimmer that continuously aerates water (salinity: 31.5 ± 0.5‰, temperature: 19 ± 0.5°C, pH: 8.2 ± 0.1) at Ludong University for 2 weeks before the experiments. Plastic plants were used as holdfasts. The seahorses were fed three times per day (08:00, 12:00, and 16:00) with frozen Mysis, and residual feed and feces were siphoned out 2 h after each feeding session.

To construct the research model, seven groups were set up (20 male and 20 female seahorses per group) to determine the appropriate combination of seahorse size (wet weight, g) and challenge dose of E. piscicida (cfu/mL) via intraperitoneal (IP) injection (14). In detail, Con represents the combination of 2.5–3.0 g and physiological saline solution, and EPs represent the E. piscicida challenge groups with the combinations 3.5–4.0 g and 1 × 10⁵ cfu/mL, 3.5–4.0 g and 1 × 10⁷ cfu/mL, 3.5–4.0 g and 1 × 10⁹ cfu/mL, 4.1–4.5 g and 1 × 10⁵ cfu/mL, and 4.1–4.5 g and 1 × 10⁷ cfu/mL. After the IP injection of E. piscicida, the male and female seahorses were separately maintained in tanks (50 × 40 × 40 cm³) under the same culture conditions mentioned above, except that the circulating water was turned off. After 24 h, the seahorses were fed normally, feces were siphoned out, and seawater was supplemented to maintain normal water levels. The number of deaths was recorded daily to calculate survival rate, and the seahorses were switched to new tanks every week for 21 days. The appropriate combination for research model construction was determined according to the survival rate and pathological characteristics of the seahorses.

Using the optimal combination of seahorse size and challenge dose, 200 big-belly seahorses (male: female = 1:1) were equally separated into two groups (Con and EP groups) with 50 males and 50 females per group, and a new research model was constructed. Both growth-related and pathological parameters of the big-belly seahorses listed in the established evaluation system were recorded on days 0, 1, 5, 9, 15, and 21 as previously reported (25). In addition, the respiratory rate (RR) of randomly chosen seahorses per group was recorded daily. After anesthetization with 0.035% MS-222 (Sigma-Aldrich, Saint Louis, Missouri, USA) for 2 min, intestinal samples were collected on days 0, 1, 5, 9, 15, and 21 for further histological and quantitative real-time polymerase chain reaction (q-PCR) analysis. After 21 d, the pathological characteristics of E. piscicida-induced enteritis in big-belly seahorses were examined.

The established evaluation system was updated by redefining the scoring range and supplementing the scoring system with new parameters according to our previously established principles (25). Disease activity index (DAI) was determined to reveal the pattern and key pathogenic time points of E. piscicida-induced enteritis in big-belly seahorses. To further reveal the pathogenesis of E. piscicida-induced enteritis in big-belly seahorses, 60 seahorses (male:female = 1:1) were divided into two groups (Con: 10 males and 10 females; EP: 20 males and 20 females). After repeating the model construction steps, intestinal samples were collected at key pathogenic time points for integrated metagenomics and metabolomics analyses.

The intestine samples randomly collected from the seahorses in the Con and EP groups (n = 3) were fixed in Bouin’s solution. Within 24 h, the tissues were dehydrated in an alcohol-xylene series and embedded in paraffin. The samples from different intestinal segments were cut into 8-μm-thick sections and stained with hematoxylin and eosin (Beyotime, Tianjin, Hebei, China) (25).The sections were then examined under a light microscope (BX53; Olympus, Tokyo, Japan) at 400× magnification to visualize the pathogenic changes in the intestine.

Intestine samples (n = 5) randomly collected from the Con and EP groups on days 0, 1, 5, 9, 15, and 21 were used to extract total RNA using Trizol Reagent (9109, Takara, USA). cDNA was synthesized using the PrimeScript RT Reagent cDNA Amplification Kit with a gDNA Eraser (RR047A; Takara, USA). q-PCR analysis was performed to evaluate gene expression using a SYBR^® Premix Ex Taq™ (Takara, Dalian, China) on a Bio-Rad CFX96 Touch machine (Bio-Rad, USA) according to a previously reported method (29). The sequences of the primers used are listed in Table S1 . The relevant levels of target genes were determined using the 2^−ΔΔCt method (30).

Fifteen seahorses from each group were randomly selected, and the frequency of their gill cover movements per minute was recorded to calculate the RR (31).

To eliminate the potential sex-related differences and obtain better matching results of metagenomic and metabolomic sequencing analyses, eight intestinal samples per group, one male and one female seahorses, were randomly collected and thoroughly ground in liquid nitrogen. Half of each sample was used for metagenomic sequencing, and the other half was used for metabolomic sequencing.

To reveal the effects of E. piscicida infection on the composition, diversity, structure, function, VF, and ARO of intestinal microbiota, four samples from each group were randomly selected for metagenomic sequencing, and a paired-end sequencing approach was used to obtain raw data (32). Thereafter, data quality control, removal of host genome reads, metagenomic assembly, encoding gene prediction, and non-redundancy gene set construction were performed. Detailed information regarding these steps can be found in the Supplementary Materials and Methods . The raw sequencing data of metagenome in this study are available in the Sequence Read Archive database of National Center for Biotechnology Information (NCBI) under accession number PRJNA916642.

To determine the effects of E. piscicida infection on host metabolism, liquid chromatography-mass spectrometry (LC-MS) analysis of all eight samples from each group was used for untargeted metabolomic sequencing. Chromatographic separation was performed on a Waters UPLC Acquity I-Class PLUS (Waters Corp., Milford, Connecticut, USA) system with mass spectrometric detection using a Waters UPLC Xevo G2-XS QTOF attachment (Waters Corp.). The samples were analyzed in both positive and negative ion modes. All raw data were collected using MassLynx software (version 4.2, Milford, Connecticut, USA) and entered into Progenesis QI (version 2.4, Munich, Bavaria, Germany) software for further analysis (33).

For the bioinformatic analysis, α-diversity analyses were conducted to evaluate the diversity, evenness, and richness of the intestinal microbiota, especially for the bacteria and their function, VFs and AROs, and metabolites at different levels by using GraphPad Prism (version 8.0, San Diego, California, USA), and β-diversity analyses were conducted to evaluate the structure of the intestinal microbiota, especially for the bacteria and their function, VFs and AROs, and metabolites at different levels by BMK Cloud platform (http://www.biocloud.net/). Linear discriminant analysis effect size (LEfSe) was used to screen biomarkers of intestinal microbiota using the criteria linear discriminant analysis (LDA) > 3 and P < 0.05 (34).

Network analysis was employed to visualize the potential connections between bacterial genera and VFs (35). After annotation using the Nr, KEGG, GO, eggNOG, VFDB, and CARD databases (36, 37), correlations between bacterial genera and KEGG functional pathways of intestinal microbiota, as well as VFs and KEGG functional pathways of intestinal microbiota, were assessed using Pearson correlation analysis (https://cloud.majorbio.com/page/tools/). The affiliations of VFs and AROs with bacterial species were determined using data from the non-redundant gene set construction. Potential metabolic biomarkers (PMBs) and KMBs were identified according to previously reported criteria (variable importance in projection [VIP] > 1.0, P < 0.05, and VIP > 1.5, P < 0.05) (33). The fold change (FC) of KMBs was calculated using the formula E9D/Con.

Experimental data are presented as the mean ± standard deviation and analyzed using an independent samples t-test in SPSS software (version 23.0, Chicago, Illinois, USA). Structural differences in the intestinal microbiota, function, VF levels, and AROs were determined using the PERMANOVA algorithm (http://www.biocloud.net/). Results with P < 0.05 were considered significant and those with P < 0.01 were considered highly significant.

The survival rate of big-belly seahorses showed a weight- and dose-dependent relationship ( Figure 1A ). The survival rate of the seahorses (4.1–4.5 g) injected with 1 × 10⁵ cfu/mL of E. piscicida was the same as that of the Con group (100%), which is suitable for model construction. Whereas, in the other EP groups, the time to death and survival rate appeared to positively and negatively correlate with body weight and challenge concentration, respectively. After the reconstruction of the research model, typical pathological changes associated with bacterial enteritis, such as intestinal epithelial dissolution, focal bleeding, villus atrophy, separation between the lamina propria and submucosa, thickened lamina propria and muscularis mucosae, vascular distorted congestion, and a large amount of inflammatory necrosis in the muscle layer and serosal surface, were observed ( Figure S1 ). Compared with those of the Con group, the body weight and length of seahorses in the EP group (4.1–4.5 g and 1 × 10⁵cfu/mL) were reduced, especially on day 9 (P < 0.05) ( Figures 1B, C ). The gene expression of intestinal proinflammatory cytokines (IL-1;interferon1, INF1; tumor necrosis factor-α, TNF-α, interleukin-1β, IL-1β; IL-1β receptors), anti-inflammatory factors (IL-10 and IL-2), antimicrobial peptides (hepcidin; liver-expressed antimicrobial peptide, LEAP; piscidin; lysozyme), and Toll-like receptor 5 (TLR5) significantly increased on day 9 and decreased to different extents thereafter ( Figures S2A, B ). At the same time, the gene expression patterns of ZO-1 were similar to those of the immune genes ( Figure S2C ), whereas the other two tight junction genes claudin 5 and occludin showed the opposite expression patterns ( Figures 1D, E ). The RR of the EP group declined from days 2 to 9, especially from days 3 to 9 (P < 0.05), recovered thereafter, and remained at the level of the Con group ( Figure 1F ). As the RR and gene expression of both claudin 5 and occludin varied consistently along with the pathogenic process and other parameters, we updated our previously published evaluation system by adding them to the DAI scoring system ( Table S2 ). We found that the DAI of the EP group increased from day 1, became significantly higher than that of the Con group from day 5 (P < 0.05), peaked on day 9, and decreased there after until day 21 ( Figure 1G ).

A total of 256,229 genes were identified via sequencing, belonging to 6 kingdoms, 95 phyla, 119 classes, 250 orders, 552 families, 1329 genera, and 3171 species ( Figure 2A ; Table S4 ). The number of sequences could reflect the structure and diversity of the microbial community because the corresponding rarefaction curves had already reached the saturation plateau ( Figure S3A ). The α-diversity and β-diversity of the intestinal microorganisms of the E9D group were significantly different from those of the Con and E21D groups at the species level (P < 0.01), whereas the difference between the Con and E21D groups was not significant (P > 0.05) ( Figures S3B, C ).

The intestinal microbiota composition at the kingdom level was numerically dominated by bacteria in each group ( Figure 2A ). As the strains used for stress were bacterial pathogens, all subsequent analyses were based on bacterial levels. There were 2066 shared bacterial species in the Con, E9D, and E21D groups, with 231, 227, and 380 endemic species, respectively ( Figure S3D ). As shown in Figure 2B , the β-diversity of the intestinal microbiota in the E9D group was significantly different from that of the Con and E21D groups at the species level (P < 0.01) but not between the Con and E21D groups (P > 0.05). The Shannon and Simpson indices of α-diversity in the E9D group were significantly lower than those of the Con and E21D groups (P < 0.01), but the ACE and Chao richness indices were not significantly different among the three groups (P > 0.05) ( Figure 2C ). This result indicates that pathogen infection altered the structure, and significantly reduced the evenness and diversity of intestinal microbiota in seahorses.

The composition of the intestinal microbiota at different levels changed after E. piscicida infection ( Figure S4 ). The LEfSe analysis revealed significant differences in the phylogenetic distribution and biomarkers of the microbiota among the groups (P < 0.05) ( Figure 3A ). Notably, Edwardsiella and Edwardsiella piscicida accounted for 70.10% and 32.74% of the relative abundance in the identified bacteria, respectively ( Table S5 ). Edwardsiella abundance was significantly positively correlated (P < 0.05, |r| > 0.7) with Chlamydia, Enterobacter, Yersinia, Pantoea, Salmonella, Xenorhabdus, Klebsiella, and Arthrobacter but negatively correlated (P < 0.05, |r| > 0.7) with Lactobacillus, Microbacterium, Enterococcus, Acinetobacter, Mycobacteroides, Flavobacterium, and Thalassococcus ( Figure 3B ).

The structure and activity of the functional pathways of the intestinal microbiota in the E9D group were significantly different from those of the Con and E21D groups at all three annotated KEGG, GO, and eggNOG levels (P < 0.01) ( Figure 4A ; Figure S5 ). As shown in Figure 4 , compared with that of the Con group, the relative abundance of Edwardsiella, Chlamydia, Enterobacter, and Arthrobacter (P < 0.05) and the activities of 28 positively correlated functional pathways (P < 0.01, |r| > 0.7) were significantly increased (P < 0.01), whereas the relative abundance of Lactobacillus, Enterococcus, Microbacterium, Acinetobacter, Mycobacteroides, Aeromonas, and Burkholderia and the activities of 11 positively correlated functional pathways (P < 0.05, |r| > 0.7) were significantly decreased (P < 0.05) in the E9D group. Notably, the activities of eight functional pathways, including bacterial chemotaxis, flagellar assembly, lipopolysaccharide biosynthesis, phosphotransferase system, bacterial secretion system, QS, ABC transporters, and TCS, increased significantly (P < 0.01) in the E9D group ( Figure 4 ), indicating the crucial role of these pathways during E. piscicida infection.

Metagenomic analysis identified 10 categories, 181 VFs, and 478 virulence genes ( Table S6 ; Figure S6A ). The structure of VFs in the E9D group was significantly different from that of the Con and E21D groups (P < 0.01) ( Figure S6B ), which was similar to the results of the intestinal microbiota. Significant differences in α-diversity were found between the Con and E9D groups (P < 0.01) ( Figure S6C ; Table S6 ), consistent with the results of the intestinal microbiota and their functions.

Compared with that of the Con group, the relative abundance of 128 VFs in the E9D group significantly increased (P < 0.05), including 126 VFs that were extremely significantly increased (P < 0.01) ( Figure 5A ; Table S6 ). Edwardsiella could express 123 VFs covering all 10 categories ( Table S6 ) and all six pathogenic processes. Significant positive correlations were observed among 50 VFs (P < 0.05, |r| > 0.7) ( Figure 5B ). The relative abundance of 15 of them, Flagella, Type IV pili, Lap, Bsa T3SS, Dot/Icm, FarAB, Capsule1, Hp-NAP, Legiobactin, IraAB, Hpt, LPS3, PhoP, BfmR, and BopD, was in the top 50 and complex nodes (P < 0.05) ( Figure 5B ; Table S6 ), indicating their pivotal role in regulating other VFs. As shown in Figure 5C , the top 50 most abundant VFs were significantly positively (P < 0.01, |r| > 0.7) and negatively (P < 0.05, |r| > 0.7) correlated with 28 significantly increased and 11 significantly decreased (P < 0.05) intestinal microbiota functions ( Figure 4 ), respectively, indicating the importance of these VFs in regulating intestinal microbiota function. In addition, 48 AROs were annotated ( Figure S7A ; Table S7 ). The structure of AROs in the E9D group was significantly different from that in the Con and E21D groups (P < 0.05) ( Figure S7B ), which is similar to the results of the intestinal microbiota and VFs. The Chao and Shannon indices of α-diversity in the E9D group were significantly higher than those in the control group (P < 0.01) ( Figure S7C ). Six AROs with significantly increased abundance (P < 0.05) in the E9D group were identified ( Figure 5D ), five of which were expressed by Edwardsiella, adeF, msbA, Ecol_EFTu_PLV, Ecol_GlpT_FOF, and Hinf_PBP3_BCA ( Table S7 ).

Metabolomic data revealed that the structure of metabolites of the E9D groups was different from that of the Con group ( Figures S8A–C ). A total of 17 categories, 4789 metabolites, and 1768 PMBs (VIP > 1.0 and P < 0.05; 623 upregulated and 1145 downregulated) were identified ( Figures 6A, B , S8D ; Table S8 ). At the same time, 491 KMBs (VIP > 1.5 and P < 0.05; 189 upregulated and 302 downregulated) were identified ( Figure S8E ), indicating that E. piscicida infection significantly affected metabolic processes in big-belly seahorses.

Forty-eight out of the 491 KMBs ( Table S9 ) were enriched in six major functional categories and 110 functional pathways ( Figure 6C ). The activities of 27 functional pathways from 5 major functional categories were also significantly different (P < 0.05) ( Figure 6D left panel). Eleven KMBs (abundance TOP 15 and |log2FC | > 0.55) were identified and enriched in 10 significantly different functional pathways (P < 0.05) ( Figure 6D ). Compared with those of the Con group, KMBs associated with histidine metabolism, citrate cycle (TCA cycle), and taurine and hypotaurine metabolism pathways were mostly upregulated. In contrast, KMBs associated with ABC transporters, TCS, ferroptosis, glutathione metabolism, and pyrimidine metabolism pathways were all downregulated, whereas KMBs related to glyoxylate and dicarboxylate metabolism and butanoate metabolism pathway were mostly downregulated in the E9D group ( Figure 6D ). Notably, eight KMBs, alpha-Ketoglutarate, Taurine, Glutathione, Uridine, L-Glutamate, L-Malic acid, Glutathione, and Oxidized glutathione, were involved in regulating multiple functional pathways ( Figure 6D right panel), suggesting essential roles during E. piscicida infection.

Centered on the eight closely related key functions of the intestinal microbiota, we found 34 KVFs significantly positively correlated with six of them (P < 0.01, |r| > 0.8), including15 were core VFs (CVFs) and the other 19 of the top 50 VFs in abundance. As shown in Figure 7A , 34 KVFs among the top 50 VFs with relatively complex nodes correlated with six key intestinal microbiota functions and seven KMBs after E. piscicida infection. For instance, in the E9D group, the relative abundance of VFs associated with motility, adhesion, and invasion of pathogens significantly and simultaneously increased with the functional activity of flagella assembly and bacterial chemotaxis (P < 0.01). The relative abundance of the effector delivery system-related VFs significantly and simultaneously increased with the functional activity of the bacterial secretion system (P < 0.01). The relative abundance of the iron uptake system-related VFs increased with the functional activity of ABC transporters of the intestinal microbiota (P < 0.01), whereas the levels of L-Proline, Taurine, Glutathione, Uridine, Xanthosine, and L-Glutamate significantly decreased (P < 0.05). Similarly, the relative abundance of toxin-related VFs significantly and simultaneously increased with the functional activity of LPS biosynthesis. The relative abundance of the regulation-related VFs significantly and simultaneously increased with the functional activity of TCS of the intestinal microbiota (P < 0.01), whereas the levels of KMBs L-Malic acid and L-Glutamate significantly decreased (P < 0.05) ( Figures 4D left panel; 5C ; 7A ). In addition, seven of the eight KMBs with central roles were enriched and significantly negatively correlated with two key intestinal microbiota functions, TCS and ABC transporters (P < 0.05). The content of L-Malic acid and L-Glutamate related to TCS significantly decreased after E. piscicida infection (P < 0.05), as did the levels of L-Proline, Taurine, Glutathione, Uridine, Xanthosine, and L-Glutamate in relation to ABC transporters (P < 0.05) ( Figure 7A ).