To identify novel and potent nematode chitinase inhibitors, in this study, a total of 8 million molecules from ZINC library were used for virtual screening. Firstly, in order to improve the screening efficiency, the pesticide-likeness rules (Hao et al., 2011) were used to reasonably reduce the number of compounds and this yielded 3176361 molecules. Then, the reported co-crystal structure of CeCht1 with PP7 (PDB ID: 6LE8) (Chen W. et al., 2021) was selected as the template to generate the pharmacophore model in MOE. We focused on the H-bonding and hydrophobic interactions between ligand and residues Glu179, Trp138, and Trp394 at −1 and +1 subsists of CeCht1. As a result, a pharmacophore model was constructed with the features of a hydrogen bond donor at Glu179, aromatic center at Trp394, a hydrogen-bond acceptor and a hydrophobic center at Trp138 (Figure 2). Based on this pharmacophore model, 1376261 compounds were obtained.

To further reduce the number of the hits, these molecules were refined with docking-based virtual screening. The virtual screening was performed using the MOE modeling software, and the crystal structure of CeCht1 complexed with PP7 was selected as the template for molecular docking. The obtained poses were rank-ordered according to their binding scores (a higher negative score indicated a higher stable interaction). As a result, the top 500 compounds with score ranging from −9.75 to −9.00 were selected. These hits were further refined with visual inspection. Considering the interactions between CeCht1 and the corresponding inhibitors revealed from crystal complexes, only hits that mainly interacted with residues Trp138, Asp177, Asp179, Tyr247, Asp248 and Trp394 were considered. Finally, 15 representative compounds were purchased from Topscience Corporation and were used for biological activity assays. Structures, docking scores of the compounds are listed in Supplementary Table S1. The corresponding sequence of these steps is summarized in Figure 3.

After the 15 potential compounds were purchased, the corresponding inhibitory activities against CeCht1 were assayed at a concentration of 100 μM. As shown in Supplementary Table S1, compounds HAU-4, HAU-7, HAU-9, and HAU-11 exhibited good inhibitory potency against CeCht1 with inhibition rates >75% at 100 μM. Two compounds (namely, HAU-1 and HAU-3) showed moderate efficiency against CeCht1 with inhibition rates >50% at 100 μM. The remaining compounds showed relatively weak efficiency (inhibition rate <50% at 100 μM). These results preliminarily indicate that our virtual screening strategy is effective.

Seven compounds with inhibition rates >50% at the concentration of 100 μM were further assayed for their inhibitory activities. As shown in Table 1, compounds HAU-4 (ZINC09610803) and HAU-7 (ZINC12704597) exhibited excellent potency against CeCht1 with IC₅₀ values of 4.2 μM and 10.0 μM, respectively (Supplementary Figure S1). The remaining compounds displayed inhibition rate <50% at 10 μM against CeCht1. Therefore, compounds HAU-4 and HAU-7 can be used as the novel nematode chitinase inhibitors, and the study of their inhibition mechanisms is of great significance for guiding the development of new nematicides.

To investigate the basis for the potency of compounds HAU-4 and HAU-7 toward CeCht1, the molecular docking studies were carried out. As shown in Supplementary Figure S2A,B,C, compounds HAU-4 and HAU-7 were all well-anchored in the substrate-binding pocket of CeCht1 at −1, +1, and +2 subsites, and stabilized by hydrophobic interactions and hydrogen bonds. Specifically, the dipyrido-pyrimidine moiety of HAU-4 was found to bind to the -1 and +1 subsites of CeCht1, and formed H-bonding interactions with residues Trp138 and Asp248 (Supplementary Figure S2B). The p-methoxyphenyl fragment of HAU-4 was located at the +2 subsite of the active pocket entrance consisting of Trp253 (Supplementary Figure S2B). For compound HAU-7, the pyrrole ring was shown to be inserted into the -1 subsite of CeCht1 via π-π stacking interactions with Trp394, and the amide linker can form H-bonding interactions with residues Trp138 and Asp248 (Supplementary Figure S2C). In addition, the naphthalene ring of HAU-7 extends to the +2 subsite (Supplementary Figure S2C).

To shed further light on the appropriate binding modes of compounds HAU-4 and HAU-7 with CeCht1 and in an effort to assay the stabilities of the potential compounds in the active pocket, 40 ns MD simulations were performed (Figure 4). As shown in Figure 4A, the two simulated systems all achieved dynamic equilibrium at approximately 25 ns and their root-mean-square deviation (RMSD) values ultimately maintained at 1.7 and 1.5 Å for protein-ligand complexes of compounds HAU-4 and HAU-7, respectively. These results indicated that these two systems underwent reasonable conformational changes.

Superimposition of the conformations of HAU-4 and HAU-7 with CeCht1 at 40 ns of MD simulations are shown in Figure 4B and Supplementary Figure S2D. Compound HAU-7 was found to be buried at the entrance of the CeCht1 active pocket, and mainly binds to the +1 and +2 subsites of CeCht1. In contrast, compound HAU-4 was found to enter deeper into the active pocket, and the dipyrido-pyrimidine ring was inserted into the -1 subsite of CeCht1 (Figure 4B and Supplementary Figure S2D).

The specific binding mode of HAU-4 with CeCht1 revealed by MD simulations is shown in Figure 4C. In comparison to the docking conformation (Supplementary Figure S2B), the dipyrido-pyrimidine ring of HAU-4 was demonstrated to move deeper into the CeCht1 pocket via 40 ns MD simulations, and the p-methoxyphenyl moiety was far from the +1 and +2 subsites (Supplementary Figure S2E). In detail, the dipyrido-pyrimidine moiety of HAU-4 stacked well with Trp394, Phe398, Tyr302, and Trp62 at −1 subsite. The oxygen atom on dipyrido-pyrimidine ring formed a hydrogen bond with Trp62. On the other hand, the p-methoxyphenyl fragment of HAU-4 bound in a small hydrophobic cave constructed by residues Tyr101, Phe89, Trp138, and it was further stabilized by forming a hydrogen bond with Tyr101. Furthermore, we found that the ethyl group on the dipyrido-pyrimidine points to the position of the catalytic residues Asp175, Asp177, and Glu179 (Chen Q. et al., 2021). This may suggest the introduction of suitable substituent (such as fragment with hydrogen bond donor) in the ethyl of HAU-4 may lead to increased activity. Thus, this binding pattern provided a clear molecular basis for the potent inhibitory activity of compound HAU-4 (IC₅₀ = 4.2 μM) against CeCht1.

The conformations of compound HAU-7 before and after MD simulations in CeCht1 pocket are superimposed and shown in Supplementary Figure S2F. Compound HAU-7 was found to exhibit a great deal of folding and move to the entrance of the pocket relative to the conformation before MD simulations (Figure 4D, Supplementary Figure S2C,F). The specific binding mode of HAU-7 with CeCht1 revealed by MD simulations showed that the pyrrole moiety anchored in the hydrophobic pocket and formed π-π stacking interactions with Trp138, while the naphthalene ring interacted with Ala252 and Trp253 at +2 subsite of CeCht1. It is worth noting that HAU-7 could only form one hydrogen bond with Trp253, with reduced interactions in its docking mode (Figure 4D and Supplementary Figure S2C). Therefore, the affinity between HAU-7 with CeCht1 was determined to be relatively weak, resulting in a moderate inhibitory potency against CeCht1 (IC₅₀ = 10.0 μM).

To further explore the cause of the inhibitory potency against CeCht1, binding free energy analyses of HAU-4 and HAU-7 were conducted by using MM/GBSA calculation methods. The predicted binding affinities of these complexes are summarized in Table 2. The results showed that the binding free energy (ΔG _TOT) were in agreement with the experimental values of HAU-4 and HAU-7. Furthermore, both of the electrostatic (E _ELE) and the van der Waals (E _VDW) contributions are essential for the two compounds when bound to CeCht1. For compound HAU-4, the electrostatic energy (E _ELE = −30.53 kcal/mol) was lower than its van der Waals energy (E _VDM = −29.42 kcal/mol). On the contrary, the van der Waals interactions (E _VDM = −26.63) of HAU-7 acted a more important role than the electrostatic interaction kcal/mol (E _ELE = −13.83 kcal/mol). Moreover, the nonpolar interactions (E _VDM + G _SA) of these two complexes were lower than the polar interactions (E _ELE + G _GB), which suggested that the nonpolar interactions in these two systems were more beneficial for the binding of HAU-4 and HAU-7 to CeCht1.

Free energy decomposition calculations were further performed by using MM/GBSA methods to investigate the detailed contributions of residues in the two systems. Residues near 5 Å of the ligands were selected to compute the decomposition free energy, and the main residue contributions are shown in Figure 5. Compound HAU-4 formed the strongest interaction with Trp62 (−4.32 kcal/mol) and formed moderate interactions with Trp394, Tyr101, Tyr302, and Trp138 (−1.94, −1.92, −1.57, and −1.35 kcal mol⁻¹, respectively). This suggested that the hydrophobic interactions and hydrogen-bonding interactions at −1 subsite were the main reasons for the higher inhibitory potency of compound HAU-4. Compound HAU-7 mainly interacted with Trp253, Val334, and Ala252 at the entrance of the active pocket. This result further demonstrated the shallow binding mode of HAU-7 in the CeCht1 active pocket. It is noteworthy that neither compound bound to the catalytic residue Asp175, Asp177, and Glu179 (Chen Q. et al., 2021), illustrating that the efficiency of compounds HAU-4 and HAU-7 is in the micromolar range.

The ZINC database containing 8 million molecules was used for our virtual screening targeting nematode chitinase CeCht1. Initially, the descriptors of all the molecules in the ZINC library were calculated via the Descriptor module in MOE software (version 2016). Secondly, these molecules were filtered based on “Hao’s pesticide-likeness” rules (Hao et al., 2011), which are MW < 435 Da, ClogP <6, HBA <6, HBD <2, ROB <9, ARB <17. In addition, the synthetic feasibility score in MOE of these compounds was considered. The synthetic feasibility score is the fraction of atoms of each new structure that ultimately appear in a retrosynthetic fragment, and a value of one means the molecule is likely to be synthesizable. Based on these criteria, the unsuitable molecules were gradually filtered out.

The complex crystal structure of CeCht1 (PDB ID: 6LE8) (Chen W. et al., 2021) was selected as the template to generate a pharmacophore model. The protein ligand interaction fingerprints (PLIF) tool in MOE was used to establish a pharmacophore model. The pharmacophore features take the spatial position and interactions between ligands and proteins into account, including hydrogen bond donor and acceptor, and hydrophobic interaction characteristics. Then, virtual screening module in MOE was used to perform the pharmacophore-based screening.

Molecular docking was carried out by using the Dock module in MOE (Dong et al., 2021). The complex crystal structure of CeCht1 with PP7 (PDB ID: 6LE8) (Chen W. et al., 2021) was used as the template for the docking-based screening, and the position of the ligand was defined as the docking site. First, the Structure Preparation module in MOE was used to prepare the protein structure. The typical steps in this process included the geometry and electron-density checks, addition of hydrogen atoms, optimization of their position, and energy minimization. Then, the Triangle Matcher placement method and Rigid Receptor post-placement refinement method were used to generate the poses of ligand. The number of poses returned by each ligand placement were set to a default value of 30. The London dG scoring function was used to estimate their free energy of binding. The other parameters were set default values. The molecules were ranked lowest to highest according to their scores, and the pose with the lowest score was retained.

Nematode chitinase CeCht1 (CeCht1-CAD) was overexpressed in Pichia pastoris GS115 (Invitrogen, Carlsbad, CA) and then purified as described previously (Chen Q. et al., 2021).

The target compounds were assayed for their inhibitory activities against CeCht1 in end-point experiments using 4-methylumbelliferyl -N, N′-diacetyl-β-D-chitobioside (4-MU-(GlcNAc)₂, Sigma, St. Louis, MO) as the substrate. The CeCht1 was assayed in 20 mM sodium phosphate buffer (pH 6.0). In a final volume of 100 μl, the reaction mixture containing buffer, DMSO (1% (v/v)), CeCht1-CAD (10 nM), 4-MU- (GlcNAc)₂ (4 μM) and inhibitor was incubated at 25°C for 20 min. Then, the reaction was stopped by the addition of sodium carbonate solution (0.5 M, 100 μl). The fluorescence of the liberated 4-methylumbelliferone was quantitated at an excitation of 366 nm and emission of 440 nm. Experiments were performed in triplicate. For determination of the IC₅₀ values, the inhibitory activity of the compounds against CeCht1 was monitored by changing the concentrations of inhibitors.

Molecular Dynamics (MD) simulations were carried out to assess the validity of the binding interactions of the selected compounds and CeCht1 by using AMBER14 package. The ligands and protein were selected for GAFF force field (Duan et al., 2003) and AMBER03 force field (Wang et al., 2004), respectively. The protein-inhibitor system was immersed in a radius truncated octahedral box with TIP3P water molecules extended 10 Å from the complex, and the system was neutralized by adding counterions (Na⁺). Then, two-step energy minimization was performed using the first 2,500 steps steepest-descent and last 2,500 steps conjugated gradient algorithms in Sander module. Next, the system was gradually heated from 0 to 300 K in the NVT ensemble and the temperature was maintained at 300 K during 100 ps with a restrain force constant of 5 kcal/mol/Ų. The SHAKE algorithm was used to constrain the hydrogen atoms (Ryckaert et al., 1977), and the particle mesh Ewald algorithm was applied to calculate the long-range electrostatic interactions with default cutoff of 8.0 Å under periodic boundary conditions (Darden et al., 1993). Finally, 40 ns MD simulations were carried out at the constant temperature and pressure (300 K and 105 Pa) using the PMEMD module in the AMBER14.

The molecular mechanics/generalized Born surface area (MM-GBSA) molecular method was performed to calculate the binding free energy (ΔG _bind) in the post-processing trajectory analysis in Amber14. The last 5 ns simulations were selected as the binding equilibrium conformation to calculate binding free energy. The ΔG _bind of CeCht1-inhibitor complexes was evaluated by energy minimization of the ligand, receptor, and the complex structure as follows: ∆ G b i n d = G c o m p l e x – ( G r e c e p t o r + G l i g a n d ) , ( 1 ) ∆ G b i n d = ∆ E M M + ∆ G s o l v a t i o n – T ∆ S , ( 2 ) ∆ E M M = ∆ E i n t e r n a l + ∆ E e l e c t r o s t a t i c + ∆ E V D W , ( 3 ) ∆ G s o l v a t i o n = ∆ G p o l a r   s o l + ∆ G nonpolarsol , ( 4 ) where G _complex, G_receptor, and G_ligand represent the free energy of complex, CeCht1, and ligand, respectively. ΔG _bind was evaluated by gasphase interaction energy (ΔE _MM), solvation energy term (ΔG _solvation) and entropy term (TΔS) between the candidate molecules and CeCht1. ΔE_MM is the gasphase interaction energy and it contains internal (ΔE _internal), electrostatic (ΔE _{electrostatic}), and van der Waals energies (ΔE _VDW). ΔG _solvation is the solvation free energy, which is the sum of the polar solvation contribution and nonpolar solvation contributions. TΔS shows the change of conformational entropy upon ligand binding. The energy decomposition was carried out to evaluate the contribution of each residue of the CeCht1 via the MM-GBSA method in AMBER14.