Since the protein structure of the TonB-dependent siderophore receptor of F. columnare was not available in the PDB database, the SWISS-MODEL, Phyre2, Robetta, and ColabFold web servers were introduced to generate protein structural data (Kelley et al., 2015; Waterhouse et al., 2018; Baek et al., 2021; Mirdita et al., 2022). During the modeling process, all parameters were selected as default parameters. The predicted protein structure was modeled after the amino acid sequence of the TonB-dependent siderophore receptor in the reference genome of F. columnare in the NCBI database (WP_065213071.1) (Agarwala et al., 2017). To validate the model quality, the SAVES web server (Available online at https://servicesn.mbi.ucla.edu/SAVES/, accessed July 10, 2022) was used to analyze the stereochemistry of the obtained protein models and to assess the model quality. The model with the highest overall score was selected for molecular docking. The protein structure was processed with the Protein Preparation wizard in Schrodinger suite. Then, the AUTODOCK 4.0 MGL tool was used to convert the protein PDB files into the required PDBQT format.

The 3D structures of the molecular ligands used in this study were downloaded from the Traditional Chinese Medicine Library (Catalog No. L8300, Selleck Chemicals, Houston, TX, USA). Openbabel was used to split the downloaded SDF files and convert them to the PDB format (O'Boyle et al., 2011). Finally, the AUTODOCK 4.0 MGL tool was used to convert the file to the PDBQT format, assuming a neutral condition (pH = 7.4) to satisfy the requirements for bacterial growth.

We uploaded the protein 3D structure files to the CASTp 3.0 server to predict active pocket locations based on geometric and topological properties (Tian et al., 2018), and labeled amino acid residues to visualize pockets in MGLTools 1.5.6 (https://ccsb.scripps.edu/mgltools/). The position of the catalytic site was manually adjusted to completely wrap the pocket. The final parameters for the catalytic sites were set as center-x = 2.472, center-y = 9.250, and center-z = −1.556, with a grid size of 29 × 29 × 29. Autodock Vina (version 1.2.2) and R (version 4.2.3) were utilized to perform high-throughput molecular docking-based virtual screening (Trott and Olson, 2010). Throughout the docking, the ligand was in a flexible form, while the protein was in a rigid form. The maximum number of binding modes to output was fixed at 9, and the exhaustiveness level (controlling the number of independent runs performed) was set at 8. The conformation with the best affinity was selected as the final docking conformation. The molecular dynamics simulations were carried out with Desmond/Maestro non-commercial version 2022.1 as a molecular dynamics software to evaluate the stability and binding interactions of the complexes (Bowers et al., 2006). TIP3P water molecules were added to the systems, which were then neutralized by a 0.15 M NaCl solution. After the minimization and relaxation of the system, the production simulation was performed for 100 ns in an isothermal–isobaric ensemble at 300 K and 1 bar. The trajectory coordinates were recorded every 100 ps. The molecular dynamics analysis was performed using the Simulation Interaction Diagram from Desmond. Pymol 1.8.6 was used for the visualization of the complexes (Lashkov et al., 2021). Furthermore, the top 20 ranked compounds were selected for further property evaluation.

Online platforms were used to predict the drug properties, bioactivity score, ADME properties, and toxicity profiles of the compounds. Specifically, the ADMETlab 2.0 and SWISS ADME were used to predict the drug properties, ADME properties, and toxicity profile, while Molinspiration, a molecular property calculator based on Java (Available online at http://www.molinspiration.com/, accessed July 22, 2022), was used to predict the bioactivity scores (Daina et al., 2017; Xiong et al., 2021). Florfenicol, a commonly used aquatic antibiotic, was selected as a control to eliminate compounds with poor drug-like properties. The screened compounds participated in the subsequent MIC (minimum inhibitory concentration) assay.

The experimental strain of the pathogen F. columnare was isolated and purified from the gill tissue of dying diseased fish via inoculation with Cytophaga agar containing 0.5 g/L Tryptone, 0.2 g/L Beef extract, 0.5 g/L Yeast extract, 0.2 g/L Na-acetate, and 1% Agar. The strain was maintained as a frozen stock at −80°C in 50% (v/v) glycerol at Sichuan Agricultural University. After resuscitating the frozen stock culture, a single colony was transferred from the Cytophaga agar to 10 mL of Flavobacterium columnare growth media (FCGM, containing 8.0 g/L Tryptone, 0.8 g/L Yeast extract, 1.0 g/L MgSO4·7H₂O, 0.74 g/L CaCl₂·2H₂O, 5.0 g/L NaCl, 1.5 g/L Sodium Citrate) broth by picking it up with an inoculating loop and cultured at 28°C with shaking at 200 rpm for 24 h. The strain was identified as F. columnare through DNA extraction, PCR amplification of the 16SrDNA fragment (forward primer: 5′-AGAGTTTGATCCTGGCTCAG-3′; reverse primer: 5′-GGTTACCTGTTACGGACTT-3′), and a BLASTN comparison of the NCBI database (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

A determination of the MIC of different small-molecule compounds against F. columnare was carried out in 96-well U-bottom polystyrene microtiter plate using the broth dilution method. All compounds were purchased from Selleck Chemicals and dissolved in dimethyl sulfoxide (DMSO) to prepare a 10 mM master batch. The 2-fold serial dilution of each compound ranging from 128 to 0.0625 μg mL⁻¹ was used for the assay (the maximum concentration of Uvaol, which is extremely insoluble in water, is 32 μg mL⁻¹). The FCGM broth (pH = 7.4) was added to all wells as the medium and 0.5 McFarland bacteria were cultured for 24 h (Farmer, 2004). The bacterial growth was observed by using 0.015% Resazurin dye to examine the viability of the bacterial cells. FCGM broth supplemented with DMSO and cultures without the compound addition were used as the negative and positive controls, respectively. Meanwhile, the MIC of Florfenicol was tested in the same way against F. columnare.

F. columnare strains were streaked from freezer stocks onto Cytophaga agar and incubated at 28°C for 24 h. Growth from the Cytophaga agar was then used to inoculate 10 mL of FCGM broth, which was incubated at 28°C with rotation to 0.5 McFarland. Enoxolone was added to the medium to a final MIC (8 μg/mL). For some experiments, ferric chloride (FeCl₃) and ferrous sulfate (FeSO₄) were added to 10 mL of FCGM cultures containing Enoxolone, respectively. Ferric chloride and ferrous sulfate were added in increments of 50 μM or 100 μM. Cultures without Enoxolone addition were used as controls. All cultures were set up in three parallels and incubated for 24 h at 28°C. We added 100 μL of the incubated culture to a 96-well plate. The cell viability was measured in triplicate using CCK-8 (Beyotime Biotechnology, Shanghai, China) at a dose of 10 μL per well⁻¹. The absorbance values were measured at wavelength of 450 nm using an ultraviolet–visible spectrophotometer according to the manufacturer's protocol.

Sixty healthy largemouth bass (0.99 ± 0.31 g) were purchased from a farm in Sichuan. They were acclimated in the laboratory for 2 weeks prior to the experiment. All largemouth bass were divided into 6 groups (1 control and 5 experimental groups), each of which contained 10 largemouth bass. The control group received DMSO immersion (1%), while the experimental group was immersed in aeration water containing 4 μg/mL, 8 μg/mL, 16 μg/mL, 32 μg/mL, and 64 μg/mL of Enoxolone, respectively. During the experiment, the largemouth bass were exposed to an uninterrupted oxygen supply, the pH was marinated at 6.5–8.0, and water temperature was about 25°C. The experimental period was set to 96 h, and the number of deaths and clinical symptoms of the largemouth bass were recorded daily.

Data on the cell viability were analyzed and expressed as mean ± standard deviation. A statistical analysis was carried out using an analysis of variance (ANOVA) followed by multiple comparison tests (LSDs). If the variance among groups was not homogenous, Dunnett's T3 test was performed. A value of p < 0.05 was considered statistically significant and labeled with an asterisk (^*). GraphPad Prism 8 was used for the graphical representations.

The TonB-dependent siderophore receptor protein lacks PDB entry; for this reason, the structure of the protein has been depicted through homology modeling. Amongst the 701 amino acids of TonB-dependent siderophore receptors, 5.85%, 49.50%, and 44.65% of these amino acids construct α-helices, β-strands, and transmembrane helices, respectively. SWISS-MODEL, Phyre2, Robetta, and ColabFold predicted the amino acid data of the proteins to generate structural images. Supplementary Figure S1 shows the 3D conformation of the F. columnare TonB-dependent siderophore receptor proteins constructed using different procedures. All protein models have a similar structure, namely a β-barrel structure formed of multiple β-folds and a plug structure contained in its center. The reasonableness of all modeled target protein topologies was analyzed using SAVES web servers. A Ramachandran plot was used to visualize whether the dihedral angles of the amino acid residues in the homologous modeled protein structures were in reasonable regions. A more rational protein structure has more amino acid residues in the best regions. According to the Ramachandran plot, the percentages of amino acid residues in the best regions for the SWISS-MODEL, Phyre2, Robetta, and ColabFold rank1 models were 84.5%, 77.9%, 88.5%, and 86.6%, respectively. In addition, non-bonded interactions between different atom types were analyzed in the ERRAT program (Colovos and Yeates, 1993). As an important indicator for assessing the three-dimensional structure of proteins via crystallography, the high factors obtained in ERRAT represent the high resolution of the structure. The results showed that the ColabFold rank1 model protein used for the analysis received the highest score of 83.847, followed by SWISS-MODEL (81.719), Robetta (81.584), and Phyre2 (49.461) (Supplementary Figures S2–S5). The PAE (prediction alignment error) results indicated that the ColabFold rank1 model protein has more confidence in the expected positions of amino acid residues in the five 3D structures generated by ColabFold, indicating that the ColabFold rank1 model protein has the highest confidence (Supplementary Figure S2D). In combination with the above data, the ColabFold rank1 model protein was used as the receptor protein for the subsequent screening.

The crystal structure of the target protein was not resolved. Meanwhile, the target protein has low identity with the protein in the PDB database and is highly conserved in the genus (Supplementary Figure S6). Therefore, CASTp was used to characterize and measure the geometric and topological properties of protein structures to complete the visualization of TonB-dependent siderophore receptor protein pocket channels. The results of the CASTp calculation are shown in Supplementary Figure S7. After localizing the active pocket, AutoDock Vina was used to dock the unique conformation of a total of 1,335 TCM ligands to the active pocket in the TonB-dependent siderophore receptor. The docking scores of each protein–ligand complex generated via molecular docking were ranked to screen for stably bound ligands. As revealed in Table 1, the top 20 screened hits were listed according to their docking scores. In addition to the highest docking score that Solamarine (CAS: 20311-51-7) showed in terms of its docking with the receptor, all other compounds had docking scores above −10.9 kcal/mol. To more fully cognize the similarities and differences among the top 20 screens, compounds were searched by type, effect, and target of action to create a Sankey diagram. The results showed that the top 20 compounds were mostly pentacyclic triterpenoids with similar structures (linked by six isoprene units) (Figure 1). Most of the compounds have anticancer and anti-inflammatory effects, with the most frequent targets being enzymes (14; 70%), followed by transporters (four; 20%) and membrane receptors (two; 10%).

Excellent drug-likeness and biological activity contribute to reducing the risk of failure of drug candidates in the production application stage. The ADME/T and biological activity of the top 20 ranked compounds obtained from the virtual screening were predicted and further filtered using Florfenicol as a standard. The results show that the 14 compounds conform to Lipinski's rule and have relatively good characteristic parameters (Table 2). In addition, most compounds do not cross the blood–brain barrier and do not inhibit CYP family enzymes, showing a high safety profile (Table 3). The Molinspiration server was used to calculate the bioactivity scores of the compounds on GPCRs (G protein-coupled receptors), nuclear receptors, ion channels, kinases, proteases, and other enzymes. The results showed that all the compounds screened had more than moderate activity (−5 to 0) against the predicted targets, with some showing significant activity (score > 0) (Table 4). The results of the toxicity prediction show that some of the compounds have potential cardiotoxicity, carcinogenicity, mutagenicity, and aquatic toxicity (Table 5). Since the compounds will be used in an aqueous environment, Lipinski's rule, the blood–brain barrier, CYP family enzymes, and toxicity (hERG, AMES Toxicity Carcinogenicity, and LC50FM) were set as screening conditions in this study to make the compounds bioenvironmentally friendly. Ultimately, Solamargine, Asiaticoside, Obacunone, Tenuifolin, Pseudoprotodioscin, Bardoxolone, Protodioscin, Dioscin, and Saikosaponin A were removed due to ineligible drug-likeness, and the remaining compounds were used in subsequent assays for antimicrobial activity.

We tested the antibacterial activity of the primary screened small-molecule compounds against F. columnare using a microdilution method based on Resazurin staining. The results showed that F. columnare showed strong resistance to Hederagenin, α-Boswellic acid, Waltonitone, Madecassic acid, Oleanonic Acid, and Uvaol. In contrast, Enoxolone, Oleanolic acid, Ursonic acid, Corosolic acid, and Pomolic acid showed antibacterial activity. Among all compounds, Enoxolone had the smallest MIC, which was 8 μg/mL (Figure 2). Meanwhile, the MIC of Florfenicol, the control drug, was determined to be 2 μg/mL. Therefore, Enoxolone was selected as a potential lead compound for subsequent experiments.

The dynamic stability and interactions of the TonB-dependent siderophore receptor–Enoxolone complex were further evaluated with molecular dynamics simulations. To study the dynamic stability of the complex, root mean square deviations (RMSDs) from the starting structures were analyzed. The plots showed that the protein and ligand were in a steady state during the simulations. After equilibration, the RMSD of the protein stabilized at 2.6 Å to 2.8 Å, whereas the RMSD of the ligand stabilized at 1.6 Å to 2.0 Å (Figure 3A). Then, the overall calculated root mean square fluctuations (RMSFs) of the protein and ligand were also small (Figures 3B, C). These small fluctuations indicate the high stability of the initial conformation of the complex. In addition, the RMSD and RMSF between apo-protein and protein-ligand complex were also nearly (Supplementary Figure S8). After docking, ARG647, SER675, GLU688, and TRP693 bind to the ligand via hydrogen bonding. During the simulation, the protein and small molecules formed multiple sets of interactions, including hydrogen bonds, hydrophobic interactions, and water bridges (Supplementary Figure S9). ASN 85, ARG 629, ARG 647, GLU 688, ASN 692, and ARG 694 were recognized as key amino acid residues in the binding process. In particular, the hydrogen bond formed by ARG 629 with small molecules was believed to play a crucial role (Figure 3D). Therefore, Enoxolone displayed a high binding affinity for TonB-dependent siderophore receptors in this conformation.

To further explore the potential mechanism of Enoxolone inhibition in F. columnare, six prokaryotic transcriptome sequencing libraries were constructed, containing three DMSO-treated controls and three Enoxolone-treated groups. The six libraries were sequenced, and the data are summarized in Table 6. An average of 25,720,127 raw reads were obtained per sample, with Q20 and Q30 both over 93%. The percentage of rRNA assessed against the reference genome was < 10% for all samples. In addition, the sequenced sequences were evenly distributed across the genes, indicating unbiased sequencing (Supplementary Figure S10). The quality control data showed that the sequencing results met the requirements for prokaryotic transcriptome analysis. Further, all unigenes obtained in the assembly were blasted against six databases (NR, Swiss-Prot, Pfam, COG, GO, and KEGG), among which 2,715 genes were annotated (Supplementary Figure S11). The expression levels of the genes were normalized to FPKM and TPM. Prior to the DEG analysis, a principal component analysis (PCA) was performed based on the expression matrix. The results showed that there was a significant difference in gene expression between the different groups (Supplementary Figure S12A). DEGs were analyzed based on the read counts, and 76 DEGs (up: 52; down: 24) were identified using edgeR (p-value < 0.05 and |log2FC| ≥ 1 in normalized expression values) (Supplementary Figure S12B).

GO and KEGG annotation were performed to further understand the biological functions of the DEGs. Supplementary Figure S13A shows the annotation of the top 10 terms of the DEGs in the GO database. The annotated results showed that the most significant changes in the expression of genes related to integral components of the membrane were observed after Enoxolone treatment. In addition, there were differences in the expression of genes associated with the DNA binding, plasma membrane, and phosphorelay signal transduction systems. KEGG annotation showed that genes related to membrane transport and amino acid metabolism showed the most significant changes in expression after Enoxolone treatment (Supplementary Figure S13B). The GO and KEGG pathway enrichment analyses were conducted by scripting with the R language, and the pathways were considered to be significantly enriched when p < 0.05. The GO analysis indicated that the main differently expressed pathways in the Enoxolone-treated F. columnare included amino acid anabolism, biogenic amine synthesis, indole-containing compound synthesis, and signal transduction (Figure 4A). The bacterial secretion system, phenylalanine, tyrosine and tryptophan biosynthesis, and cationic antimicrobial peptide (CAMP) resistance were the main enriched KEGG pathways (Figure 4B). Due to the limited number of pathways in the database for the enrichment analysis, DEGs were further analyzed on a case-by-case basis. The analysis showed that the TonB-dependent siderophore receptor and the biopolymer transporter ExbD were labeled as up- and down-regulating DEGs, respectively. In addition, the expression of genes related to drug efflux (TolC, ABC transporter, EmrA) and amino acid synthesis (TrpD, TrpC) was significantly up-regulated; in contrast, the expression of genes related to iron ion transport (TonB receptor plug, ExbD, ExbB, FeoB), DNA synthesis (DNA primase), and the type VI secretion system (VgrG) was significantly down-regulated (Figure 4C). Notably, the TonB-dependent siderophore receptor is expressed in the opposite trend to its plug structure. The down-regulation of the expression of the plug structure as a binding site may indicate that the negative feedback regulation of the strain is due to the inhibition of the activity of the TonB-dependent siderophore receptor upon binding to Enoxolone.

The gene response network of F. columnare following Enoxolone treatment was mapped against the screened genes to clarify the mechanism of inhibition. The results showed that the TonB-dependent iron transport pathway was inhibited after Enoxolone treatment, and that Enoxolone entering the bacterial cytosol via the TonB-dependent siderophore receptor activated the drug efflux system and the amino acid synthesis pathway. In addition, the relative decrease in iron ion content in bacteria inhibited DNA synthesis and reduced bacterial virulence (Figure 5A). According to the mapping of this gene regulation network, we hypothesized that the mechanism of Enoxolone resistance in F. columnare may be related to its competitive binding to the TonB-dependent siderophore receptor, resulting in an inadequate uptake of iron ions by bacteria. To demonstrate the accuracy of the proposed hypothesis, we first performed molecular docking with Enoxolone for other proteins in the TonB-dependent iron transport pathway that were significantly down-regulated. The results showed that the binding score of these proteins to Enoxolone was as low as −9.3 kcal/mol, which is less stable than the binding of the TonB-dependent siderophore receptor to Enoxolone (Figure 5B). Subsequently, we added FeCl₃ and FeSO₄ after the Enoxolone treatment of F. columnare and assayed the cell viability with CCK-8. The results showed that the addition of FeSO₄ had no significant effect on cell viability, while the addition of FeCl₃ significantly increased cell viability in a concentration-dependent manner (Figure 5C). Therefore, competitive binding to TonB-dependent siderophore receptors, leading to the inadequate uptake of iron ions, can be considered as a mechanism of Enoxolone against F. columnare. At the same time, the results based on the prokaryotic transcriptome also demonstrate the accuracy of the molecular docking virtual screening.

To evaluate the toxicity, largemouth bass were selected as experimental animals and exposed to Enoxolone (0–64 μg/mL). No mortality was observed at a lower concentration (4 μg/mL). However, this increased in a concentration-dependent manner. Compared to the control group, the group treated with 64 μg/mL of Enoxolone showed 100% mortality at 12 h. Moreover, only one mortality was observed in the minimal inhibitory concentration of Enoxolone tested (8 μg/mL) compared to that in the control group during the observation period (Figure 6). The biological toxicity exhibited by Enoxolone is acceptable, considering the small size of the experimental fish used and the availability of short and spaced treatments for production applications.