The pathogenic P. plecoglossicida strain (NZBD9) was isolated from the spleen of naturally infected L. crocea with white spot disease and stored at -80°C (Luo et al., 2019). The P. plecoglossicida strain was routinely grown in Luria Bertani (LB) broth at 18 or 28°C with shaking at 220 rpm. Escherichia coli DH5α was obtained from TransGen Biotech (Beijing, China) and grown in LB broth at 37°C with shaking at 220 rpm.

The RNAi strain was constructed based on previous study (Luo et al., 2020). The sequences of five pairs of short hairpin RNA (shRNA) were designed according to the sequences of the fliS gene using RNAi Designer (http://rnaidesigner.thermofisher.com) and are provided in Supplementary Table 1 . The five pairs of shRNAs were annealed and ligated to the linearized pCM130/tac vector. The recombined pCM130/tac vector was transformed into competent E. coli DH5α by heat shock, and competent P. plecoglossicida was eventually electrotransformed. The tetracycline resistance marker (10 μg/mL; final concentration) was used for screening stably silenced clones. The mRNA levels of fliS in the five mutants were determined by quantitative real-time polymerase chain reaction (qRT-PCR).

Overnight culture of P. plecoglossicida in LB broth was adjusted to an optical density at 600 nm (OD_{600 nm}) of 0.1 and diluted 10 000-fold with fresh LB broth. Aliquots (200 μL) of the diluted bacterial culture were added to 10 wells of a microtiter plate, then incubated at 28°C. The OD_{600 nm} of each well was automatically detected each hour for 48 h.

Swarming motility was determined following previously described methods (Zuo et al., 2019). After overnight cultivation in LB medium, OD_{600 nm} was adjusted to 0.2. The bacterial suspension (1 μL) was inoculated on semisolid agar plates. The plates were incubated at 28°C for 10–16 h to determine colony diameter. Three replicates were performed per group.

Biofilm formation was explored following previous research (Jiao et al., 2021). P. plecoglossicida in the exponential growth period was adjusted to OD_600nm = 0.2 with fresh LB broth. Then, 100 μL of diluted bacterial suspension was added into each well of a 96-well plate and incubated at 28°C for 24 h. Subsequently, each well was washed twice with 200 μL of sterile phosphate-buffered saline (PBS), stained with 125 μL of crystal violet (0.1%) for 15 min, washed three times with sterile PBS, and dried under sterile conditions. Finally, 200 μL of 33% acetic acid was added to dissolve the dyed biofilm. The OD_590nm of each well was measured, and each strain was repeated 10 times independently.

Overnight culture of P. plecoglossicida was adjusted with sterile PBS to OD_600nm = 1.0, and the bacterial suspension (0.25 mL) was then aspirated into a sterile syringe (1 mL). The open end of a capillary tube (inner diameter 0.1 mm, sealed at one end) was filled with E. coioides skin mucus (3 μL), horizontally dipped into the bacterial suspension in the syringe, then incubated at 28°C for 1 h. Finally, the number of bacteria in the mucus was determined by plate counting (Mao et al., 2020). Sterile PBS in a capillary tube was used as the control. Three independent biological replicates were performed per group.

In vitro adhesion was determined following previous study (He et al., 2021). First, sterile mucus of E. coioides (20 μL) was evenly coated onto a slide and fixed with 200 μL of 4% methanol for 30 min. Then, 200 μL of bacterial suspension (OD_600nm = 0.3) was evenly coated on the area with mucus. After the slides were incubated at 28°C for 2 h, they were washed three times with PBS to remove adherent cells. The bacterial cells adhered to the glass slide were fixed with 200 μL of 4% methanol for 30 min and stained with 0.1% crystal violet for 3 min, after which the unstained crystal violet was washed with PBS. Five slides were performed per group, and the adhered bacteria in 20 randomly selected fields were counted under a microscope (×1 000) (Leica DM4000 B LED, Leica, Germany).

Healthy size-matched E. coioides (14.7 ± 0.7 cm in length) were purchased from Zhangzhou city (Fujian, China) and adaptively maintained at 18°C for one week in laboratory seawater (water volume 320 L; pH = 7.8–8.2; continuous oxygen). The seawater temperature remained at 18 ± 1°C throughout the experiment.

For survival rate assay, 60 fish were randomly divided into three groups: i.e., fliS-RNAi strain group, wild type strain group, and PBS group. Each fish in the infection groups was intraperitoneally injected with the corresponding strain at a dose of 5 × 10⁴ colony-forming units (CFU)/fish. For the negative control, 20 fish were injected with 200 μL of PBS. Survival was observed and recorded every 12 h after injection to 10 days post infection (dpi).

For spleen sampling, 240 fish were randomly divided into three groups: i.e., wild type strain infection group, fliS-RNAi strain infection group, and PBS injection group. The fish were treated as described above. Six spleens were randomly sampled from each group at 1, 2, and 3 dpi. Two spleens from each group were randomly mixed into one sample and subjected to bacterial load and gene expression assays. Thirty-six spleens were randomly sampled at 4 dpi, with six spleens from the same group randomly mixed into one sample and subjected to qRT-PCR, pathogen load assay, fliS gene expression determination, and metabolomic analysis. All spleen samples were frozen in liquid nitrogen and stored in a refrigerator at -80°C.

qRT-PCR was performed using a QuantStudio 6 Flex Real-Time PCR system (Life Technologies, Carlsbad, CA, USA). Primers were synthesized by Tsingke Biotechnology Co., Ltd., (China) as shown in Supplementary Table 2 . The fliS gene expression levels of P. plecoglossicida were normalized using 16S rDNA. Pathogen load of P. plecoglossicida in the infected spleens was determined based on the copy number of the housekeeping gene gyrB (Xin et al., 2020). Relative gene expression in the different groups was calculated using the 2^−△△CT method (Luo et al., 2020).

An Illumina TruSeq™ RNA Sample Prep Kit was used for the RNA sequencing (RNA-seq) libraries. Total RNA was extracted from spleen samples, with RNA concentration and purity detected using a Nanodrop 2000 spectrophotometer. The RNA integrity number (RIN) was determined using an Agilent 2100 bioanalyzer. Oligo-dT was used to enrich mRNA, fragmentation buffer was added to fragment mRNA, cDNA was reverse synthesized, and adapters were ligated. Sequencing was performed using the Illumina NovaSeq 6000 sequencing platform by Majorbio Biotechnology Co., Ltd. (Shanghai, China).

Raw Illumina reads were trimmed and quality-controlled (QC) using fastp to obtain high-quality clean data to ensure the accuracy of subsequent analyses. The clean reads of each sample were compared with the reference sequences assembled by Trinity to obtain the mapping results of each sample and for subsequent unigene and transcript quantification of each sample.

To explore the related biological processes, BLAST2GO was used for Gene Ontology (GO) annotation (Diao et al., 2019) and molecular annotation of differentially expressed transcripts of samples. Differentially expressed genes (DEGs) of transcripts between the fliS-RNAi strain infection (experimental) and the wild type strain infection groups (control) were analyzed using edgeR (screening criteria: |log₂FC| ≥ 1 and false discovery rate (FDR) < 0.05). GO enrichment analysis was performed on the obtained transcripts using GOATOOLS. Signaling pathway functional enrichment was analyzed by Kyoto Encyclopedia of Genes and Genomes (KEGG) (Zhao et al., 2018). KOBAS (http://kobas.cbi.pku.edu.cn/home.do) were used for KEGG pathway analysis of differential transcripts. Ten genes were randomly selected to verify the reliability of RNA-seq by qRT-PCR.

A certain amount of sample was accurately weighed and added to the extract obtained from the metabolite extraction process in a low-temperature environment. After centrifugation, the supernatant was collected for liquid chromatography-mass spectrometry (LC-MS) detection. Six replicates were performed per group.

LC-MS was performed using ultra-performance liquid chromatography tandem Fourier transform mass spectrometry (UPLC-MS/MS) with the UHPLC-Q Exactive HF-X system (Thermo Fisher). Chromatographic conditions were as follows: Acquity UPLC HSS T3 column (100 mm × 2.1 mm i.d., 1.8 μm); Waters, Milford, USA); mobile phase A 95% water + 5% acetonitrile (containing 0.1% formic acid) and mobile phase B 47.5% acetonitrile + 47.5% isopropanol + 5% water (containing 0.1% formic acid); injection volume 2 μL; column temperature 40°C. For mass spectrometry, the sample was subjected to electrospray ionization and mass spectrometry signals were collected in positive and negative ion scanning modes. Specific parameters are provided in Supplementary Table 3 .

After UPLC-TOF/MS, the raw data were preprocessed, including original data missing value filtering, missing value filling, data normalization, QC verification, and data conversion. In addition, mass spectrometry information was exported and compared in the Human Metabolome database (HMDB) (http://www.hmdb.ca/) and Metlin database (https://metlin.scripps.edu/) to obtain specific information on metabolites.

The preprocessed data were uploaded to the MSG BioCloud platform (https://cloud.majorbio.com) for data analysis. The R package ropls (v1.6.2) was used for principal component analysis (PCA) and orthogonal least squares discriminant analysis (OPLS-DA). In addition, t-test and difference multiple analysis were performed. differentially expressed metabolites between the fliS-RNAi and wild type groups were screened (VIP > 1, FDR < 0. 05, |FC| ≥ 1), and the pathways involved in the differentially expressed metabolites were obtained using KEGG (https://www.kegg.jp/kegg/pathway.html). The Python software package SciPy was used for pathway enrichment analysis and Fisher’s exact test was used to obtain the most relevant biological pathways.

The VennDiagram package in R was used to analyze KEGG pathways involved in differentially expressed genes (DEGs) and differentially expressed metabolites. The top 10 KEGG pathways with the largest number of DEGs and differentially expressed metabolites were counted for visualization. The DEGs obtained by transcriptomics and differentially expressed metabolites obtained by metabolomics were integrated to analyze KEGG pathway data.

All data are expressed as means ± standard deviation (SD) from at least three sets of independent experiments. Data analysis was implemented using IBM SPSS Statistics v22.0 (New York, USA), with one-way analysis of variance (ANOVA) with Dunnett’s test used. P-values < 0.05 were considered statistically significant.

The RNA-seq data were deposited in the GenBank SRA database under accession numbers SRP 374606 (fliS-RNAi strain infection group), SRP 370956 (NZBD9 strain infection group), and SRP 370960 (PBS injection group). The metabolome data were deposited in the China National GeneBank Database under accession number CNP0003231.

The expression level of fliS in five different RNAi strains was verified by qRT-PCR. As shown in Figure 1A , all five shRNAs significantly reduced fliS mRNA expression by 64.2%, 76.7%, 89.0%, 77.7%, and 84.8%, respectively. Thus, five stably silenced strains of P. plecoglossicida were successfully constructed. Among them, the fliS-RNAi-126 strain (hereinafter the fliS-RNAi strain) showed the best RNAi effect and was selected for subsequent experiments ( Figure 1A ). As shown in Figure 1B , there were no significant differences in the growth curves of the fliS-RNAi and wild type strains.

The LB agar plate-cultured colonies of the fliS-RNAi and wild type strains of P. plecoglossicida are shown in Figure 2A . Results indicated that the fliS-RNAi strain colonies were significantly smaller (6.58 ± 0.11 mm, p < 0.01) than those of the wild type strain (8.37 ± 0.08 mm) ( Figure 2B ).

Biofilm formation of the fliS-RNAi and wild type strains of P. plecoglossicida is shown in Figure 2C . Compared to the wild type strain, the fliS-RNAi biofilm was lighter and deeper than that of LB. Results also showed that the OD₅₉₀ of fliS-RNAi was significantly lower (0.39 ± 0.02 nm, p < 0.001) than that of the wild type strain (0.89 ± 0.04 nm) ( Figure 2D ).

Adhesion of the fliS-RNAi and wild type strains of P. plecoglossicida to the mucus surface of E. coioides is shown in Figure 2E . Results indicated that the adhesion ability of the fliS-RNAi strain was markedly weakened, with significantly fewer (286 ± 6.20 cells/vision, p < 0.001) adhered bacteria compared to the wild type strain (814 ± 22.20 cells/vision) ( Figure 2F ).

Chemotaxis capacity of the fliS-RNAi and wild type strains of P. plecoglossicida is shown in Figure 2G . Results indicated that fliS-RNAi strain chemotaxis was marked weakened, with significantly lower (2.6 × 10⁴ CFU/mL, p < 0.001) chemotaxis compared to the wild type strain (4.7 × 10⁴ CFU/mL).

The E. coioides infected with the wild type strain of P. plecoglossicida began to die at 2 dpi, with 100% mortality reached at 6 dpi. In contrast, no E. coioides infected with the fliS-RNAi strain died by 10 dpi, as also found in the PBS-injected controls ( Figure 3A ). Furthermore, the apparent characteristics of E. coioides spleens at 4 dpi ( Figure 3B ) showed many white nodules on the spleen surface of E. coioides infected by the wild type strain, but almost no white spots on the spleen surface of fish infected by the fliS-RNAi strain or PBS. As shown in Figure 3C , although pathogen load in the fliS-RNAi strain-infected E. coioides spleens did not change significantly from 1 dpi to 4 dpi, it was consistently significantly lower than that of the wild type strain, and relative pathogen load (fliS-RNAi strain pathogen load/wild type strain pathogen load) decreased sharply with the duration of infection. The fliS gene of both strains was more highly expressed in the spleen than in vitro, and fliS expression in the fliS-RNAi strain was consistently lower than that in the wild type strain ( Figure 3D ).

RNA-seq was performed on the spleens of E. coioides infected with the fliS-RNAi or wild type strains of P. plecoglossicida. Quality control of the original sequencing data showed that the Q20 (percentage of bases with quality value ≥ 20) of each sample was above 97.41%, Q30 (percentage of bases with quality value ≥ 30) was above 93.94%, and GC content was about 50% ( Supplementary Table 4 ). The sequencing error rate was less than 0.1%, and sequencing accuracy was high. Thus, the quality of the sequencing data met the requirements for subsequent data processing and analysis.

Through the correlation analysis between biological repeat samples, on the one hand to test whether the variation between biological repeats meets the expectations of experimental design, on the other hand to provide basic reference for DEG analysis. Pearson correlation analysis was used to assess the associations between the wild type strain-infected group and fliS-RNAi strain-infected groups. The minimum and maximum Pearson correlation coefficients between groups were 0.098 and 0.164, respectively, and within groups were 0.99 and 1, respectively ( Supplementary Figure 1 ). Thus, biological repeatability was deemed good, with strong intra-group correlation and small inter-group correlation. The experimental data were therefore reliable and could be used for further analysis of DEGs.

EdgeR was used for differential transcript analysis of E. coioides. The statistical standard of mRNA expression level of transcriptome data was FDR < 0.05 and |log₂FC| ≥ 1 as GO and KEGG pathway enrichment analysis of DEGs. Transcriptomic analysis of E. coioides spleens sampled at 4 dpi identified 52 006 DEGs between the fliS-RNAi and wild type strain-infected spleens, including 38 393 up-regulated DEGs and 13 613 down-regulated DEGs ( Figure 4A ). The top 50 most significantly up-regulated and down-regulated DEGs were selected for analysis ( Figure 4B ).

Figure 4C shows the top 20 GO enrichment results. Based on GO functional annotation, the DEGs were annotated into three categories, i.e., biological process (BP), cell composition (CC), and molecular function (MF), with significant enrichment in protein-containing complex and immune system process and response to stimulus.

Further enrichment analysis identified 23 KEGG pathways involved in the immune response ( Figure 4D ), with 23.8% of DEGs enriched in immune system pathways. Based on Fisher’s exact test, 213 DEGs were significantly enriched in the cytokine-cytokine receptor interaction pathway, including 102 significantly up-regulated DEGs and 111 significantly down-regulated DEGs ( Figure 4E ). Most genes involved in this pathway were chemokines, including class I helical cytokines, class II helical cytokines, IL1-like cytokines, IL17-like cytokines, TNF family, and TGF-β family. Notably, chemokines such as CXCL8, CXCR2, CCL11, and CXCL12 were up-regulated, while certain inflammatory genes such as IL1B, IL1R2, IL17, and IL17RA were down-regulated. Furthermore, TNF gene expression was down-regulated, while TGFB1 and TGFB2 gene expression levels were up-regulated. qRT-PCR verification was performed on 10 randomly selected DEGs from the RNA-seq database. Results showed that the expression trends of the selected DEGs were consistent with those of transcriptome sequencing ( Supplementary Figure 2 ).

To determine the cumulative differences in spleen metabolites between the fliS-RNAi and wild type strain-infected groups, PCA was performed to obtain scores of the two principal components. In the total ion PCA model, PC1 = 60.30%, PC2 = 8.52% ( Figure 5A ), R²X = 0.603, and the cumulative difference explanatory value was greater than 0.5, so the established PCA model was stable and could be used for subsequent metabolic difference analysis. To display the metabolic differences between the two groups more accurately, the supervised OPLS-DA model was used for further analysis. The fliS-RNAi strain-infected samples were mainly distributed in the negative half of PC1, while the wild type strain-infected samples were mainly distributed in the positive half of PC1, so the model could effectively distinguish the two groups of samples ( Figure 5B ). The replacement test of OPLS-DA model showed that R²X and R²Y were 0.753 and 0.991, and Q² was 0.986, which were close to 1 and had good reliability ( Supplementary Table 5 ).

The differentially expressed metabolites between the two infection groups were analyzed. Based on screening criteria (FDR < 0. 05, VIP > 1, FC > 1, FC < 1), 326 differentially expressed metabolites were identified, mainly divided into organic acids, lipids, carbohydrates, nucleic acids, peptides, vitamins and coenzyme factors, steroids, hormones and transmission media, and antibiotics according to their biological function. Compared with the wild type strain-infected group, 107 and 219 differentially expressed metabolites in the fliS-RNAi strain-infected group were up-regulated (FC > 1) and down-regulated (FC < 1), respectively ( Figure 5C ). The top 20 up-regulated and down-regulated differentially expressed metabolites were selected for further analysis ( Figure 5D ).

Based on KEGG pathway enrichment analysis, 105 metabolic pathways were annotated, with 18 found to be significantly enriched (P < 0.05), including glycosylphosphatidylinositol (GPI)-anchor biosynthesis, glycerophospholipid metabolism, arginine biosynthesis, arginine and proline metabolism, ABC transporters, linoleic acid metabolism, alpha-linolenic acid metabolism, sphingolipid signaling pathway, and glutathione metabolism. Six metabolic pathways were highly significantly enriched (P < 0.01), including linoleic acid metabolism, arginine and proline metabolism, glycerophospholipid metabolism, choline metabolism in cancer, purine metabolism, and especially arginine biosynthesis pathway ( Figure 5E ).

The transcriptomic and metabolomic results were correlated. Venn analysis indicated that 101 KEGG pathways were shared between the transcriptomes and metabolomes, accounting for 94.2% of the total metabolome pathways (105) and 29.1% of the total transcriptome pathways (347) ( Figure 6A ). Figure 6B shows the top 10 KEGG pathways with the largest number of DEGs and differentially expressed metabolites. We found that arginine biosynthesis pathway was most involved ( Figure 6B ), and the correlation between metabolites and DEGs was the best ( Figure 6C ). So focused on this pathway to explore the relationship between them, further in-depth analysis was conducted.

We identified 11 DEGs and seven differentially expressed metabolites involved in the arginine biosynthesis pathway ( Figure 6D ). Among them, glnA and argH were up-regulated and glsA and rocF were down-regulated, while glutamine and arginine were up-regulated and ornithine and citrulline were down-regulated. We found that glnA gene and argH gene in the transcriptome involved in arginine biosynthesis pathway, as well as glutamine and arginine in the metabolic group were significantly up-regulated. Among them, glnA gene was positively regulated with glutamine, and argH gene was positively regulated with arginine.