Evaluation of Candidatus Liberibacter Asiaticus Efflux Pump Inhibition by Antimicrobial Peptides

Citrus greening, also known as Huanglongbing (HLB), is caused by the unculturable bacterium Candidatus Liberibacter spp. (e.g., CLas), and has caused a devastating decline in citrus production in many areas of the world. As of yet, there are no definitive treatments for controlling the disease. Antimicrobial peptides (AMPs) that have the potential to block secretion-dependent effector proteins at the outer-membrane domains were screened in silico. Predictions of drug-receptor interactions were built using multiple in silico techniques, including molecular docking analysis, molecular dynamics, molecular mechanics generalized Born surface area analysis, and principal component analysis. The efflux pump TolC of the Type 1 secretion system interacted with natural bacteriocin plantaricin JLA-9, blocking the β barrel. The trajectory-based principal component analysis revealed the possible binding mechanism of the peptides. Furthermore, in vitro assays using two closely related culturable surrogates of CLas (Liberibacter crescens and Rhizobium spp.) showed that Plantaricin JLA-9 and two other screened AMPs inhibited bacterial growth and caused mortality. The findings contribute to designing effective therapies to manage plant diseases associated with Candidatus Liberibacter spp.


Introduction
HLB is the deadliest disease threatening citrus production worldwide [1]. The disease is partially associated with a fastidious and unculturable bacterium Candidatus Liberibacter asiaticus (CLas), transmitted between trees by the Asian citrus psyllid Diaphorina citri [2,3]. There is no treatment for the disease, causing citrus trees to become unproductive and die within a few years. After its emergence in the 1920s in China, HLB widely spread to other countries in Africa and Asia. In 2005, the disease was first reported in the U.S. in South Florida. It led to an >70 percent reduction in Florida's citrus production, resulting in economic losses of more than $8 billion [4]. HLB has now spread to many citrus-producing areas in the US, including California, Florida, Georgia, Louisiana, Puerto Rico, South Carolina, Texas, and the U.S. Virgin Islands [5].
Current management practices rely mainly on controlling the insect vector that spreads the bacteria from infected to healthy plants [6]. Since there is no single effective way to control HLB fully, many approaches have been studied, developed, and are being applied to manage the disease [6]. Some treatment strategies include a combination of broadspectrum antibiotics, insecticides, nutritional supplements, and thermotherapy [7,8]. The foliar spray of pesticides can control the primary vector and transmitter of HLB, Asian citrus

Protein Structure Modeling and Validation
Multidrug efflux proteins are resistant to many small molecules, including antibiotics. T1SS is one of the known secretion systems found in CLas, which can secrete offensive enzymes and effectors [39]. Finding an effective AMP to bind and inhibit the outer membrane component may be an effective way to reduce the efflux function of the entire transmembrane protein and is the strategy that was evaluated in this study.
Since no experimental structures were available for the CLas TolC protein, the TolC protein model was prepared using the homology modeling approach based on the amino acid sequence of the CLas TolC protein and the 3D structure of a related protein. A total of 210 efflux pump templates were found by BLAST, and the selected template had 30.16% sequence identity and 0.35 sequence similarity, as shown in Table 1 and Figure S1. The trimer homology model was built by the multidrug efflux (OMP) OprN (PDB: 5iuy) [40,41]. The ERRAT value for the homology model was 92.96, and models with an ERRAT score above 50 were considered to be of good overall quality for non-bonded atomic interactions. The passed PROVE parameter showed that no buried outlier protein atoms were found in the homology model. The VERIFY 3D parameter showed that 51.23% of the residues had an average 3D-1D score >= 0.2. This parameter was below 80%, indicating that the homology model could be optimized by further processing. The predicted homology models were first preprocessed and prepared using Schrödinger's Protein Preparation Wizard. The trimeric TolC protein homology model consists of chains G, H, and I. Then, the Ramachandran plot model was generated by PROCHECK [42]. According to

Protein Structure Modeling and Validation
Multidrug efflux proteins are resistant to many small molecules, including antibiotics. T1SS is one of the known secretion systems found in CLas, which can secrete offensive enzymes and effectors [39]. Finding an effective AMP to bind and inhibit the outer membrane component may be an effective way to reduce the efflux function of the entire transmembrane protein and is the strategy that was evaluated in this study.
Since no experimental structures were available for the CLas TolC protein, the TolC protein model was prepared using the homology modeling approach based on the amino acid sequence of the CLas TolC protein and the 3D structure of a related protein. A total of 210 efflux pump templates were found by BLAST, and the selected template had 30.16% sequence identity and 0.35 sequence similarity, as shown in Table 1 and Figure S1. The trimer homology model was built by the multidrug efflux (OMP) OprN (PDB: 5iuy) [40,41]. The ERRAT value for the homology model was 92.96, and models with an ERRAT score above 50 were considered to be of good overall quality for non-bonded atomic interactions. The passed PROVE parameter showed that no buried outlier protein atoms were found in the homology model. The VERIFY 3D parameter showed that 51.23% of the residues had an average 3D-1D score >= 0.2. This parameter was below 80%, indicating that the homology model could be optimized by further processing. The predicted homology models were first preprocessed and prepared using Schrödinger's Protein Preparation Wizard. The trimeric TolC protein homology model consists of chains G, H, and I. Then, the Ramachandran plot model was generated by PROCHECK [42]. According to the PROCHECK report, 89.3% of the residues were in the most favored regions, 10.6% were in the additional allowed regions, 0.1% were in the generously allowed regions, and 0% were in the disallowed regions, as shown in Figure 2. The residue level in the Ramachandran plot-favored score was 99.9% (1091/1092), which met ideal conditions (>98%) [43]. This indicated that the protein backbones were favorable, and the proposed confirmation was stable. The Z score was predicted by comparison with the non-redundant PDB structure set of the proposed model ( Figure S1). The Z score showed that the model was less than one standard deviation from the mean models. The predicted model was considered reliable for further in silico studies based on the overall structural assessment. the PROCHECK report, 89.3% of the residues were in the most favored regions were in the additional allowed regions, 0.1% were in the generously allowed regio 0% were in the disallowed regions, as shown in Figure 2. The residue level in the chandran plot-favored score was 99.9% (1091/1092), which met ideal conditions [43]. This indicated that the protein backbones were favorable, and the proposed mation was stable. The Z score was predicted by comparison with the non-redunda structure set of the proposed model ( Figure S1). The Z score showed that the mod less than one standard deviation from the mean models. The predicted model w sidered reliable for further in silico studies based on the overall structural assessm

Virtual Screening of AMPs
Potential AMPs were screened from the APD3 database, including 2619 AM were primarily natural antimicrobial peptides. The number of AMPs in the APD3 was narrowed down by adding conditions compatible with the membrane-active p Short peptide macromolecules of less than fifteen amino acids in length were sele one of the screening criteria due to their excellent bioavailability and stability [46,4 AMPs were further selected based on previous evidence of inhibitory activity again terial efflux pumps. A preliminary screening of 15 antimicrobial peptides associate bacterial membranes was carried out from the filter search. Eight AMPs with mem activity against Gram-negative bacteria were further selected based on the experi data listed in Table 2.
The AMPs were screened based on their affinities to the β-barrel structure of th protein. The β-barrel structure consists of residues Ala 94 to Leu 128 and Tyr 300 338. All the homology models of the peptide candidates had α-helical conformatio positive net charges. These screened peptides came from various sources, such as a insects, plants, and cultured bacteria, including temporins, urechistachykinins, co plantaricin, and darobactin. Temporins are a family of antimicrobial peptides fr European red frog Rana temporaria [48]. These peptides were originally extracte frog skin secretion. The tachykinin-related peptides urechistachykinin I and II ar the echiuroid worm, Urechis unicinctus [49]. They have known effects on neurons

Virtual Screening of AMPs
Potential AMPs were screened from the APD3 database, including 2619 AMPs that were primarily natural antimicrobial peptides. The number of AMPs in the APD3 library was narrowed down by adding conditions compatible with the membrane-active peptide. Short peptide macromolecules of less than fifteen amino acids in length were selected as one of the screening criteria due to their excellent bioavailability and stability [46,47]. The AMPs were further selected based on previous evidence of inhibitory activity against bacterial efflux pumps. A preliminary screening of 15 antimicrobial peptides associated with bacterial membranes was carried out from the filter search. Eight AMPs with membrane activity against Gram-negative bacteria were further selected based on the experimental data listed in Table 2. The AMPs were screened based on their affinities to the β-barrel structure of the TolC protein. The β-barrel structure consists of residues Ala 94 to Leu 128 and Tyr 300 to Gly 338. All the homology models of the peptide candidates had α-helical conformations and positive net charges. These screened peptides came from various sources, such as animals, insects, plants, and cultured bacteria, including temporins, urechistachykinins, colistins, plantaricin, and darobactin. Temporins are a family of antimicrobial peptides from the European red frog Rana temporaria [48]. These peptides were originally extracted from frog skin secretion. The tachykinin-related peptides urechistachykinin I and II are from the echiuroid worm, Urechis unicinctus [49]. They have known effects on neurons and G protein-coupled receptors in animals. Colistin (polymyxin E) has been used for more than 50 years to treat infections caused by Gram-positive bacteria, such as Bacillus subtilis [50,51]. The bacteriocin plantaricin JLA-9 is from Suan-Tsai, a traditional Chinese fermented cabbage [52]. JLA-9 has broad-spectrum antibacterial activity against Gram-positive and Gram-negative bacteria under acidic conditions. The Trp-rich peptide TetraF2W-RK is a synthetic peptide [53]. The peptide exhibits good activity targeting bacterial membranes and inhibits the antibiotic-resistant bacteria Staphylococcus aureus [54]. Darobactin is a newly identified antibiotic isolated from Photorhabdus spp. in 2019 [55]. Darobactin is effective against antibiotic-resistant Gram-negative pathogens and has low toxicity and good pharmacokinetics in vivo.

Molecular Docking Analysis (Glide)
The conformations with the eight AMP-protein complex docking scores were analyzed for intermolecular interactions with Standard Precision (SP) peptide docking. Table 3 shows the binding affinity of AMPs ranging from −2.787 to −9.605 kcal/mol and −7.678 kcal/mol for the efflux pump inhibitor MRL-494 as the positive control [56]. MRL-494 was found as the antibiotic targeting the outer membrane protein of bacteria. All the molecules showed multiple hydrogen bonds with the protein complex. Darobactin, with the highest binding affinity of −9.605 kcal/mol, simultaneously formed five hydrogen bonds with residues ASP314 and ASN315 on chain G, residues ASN315 and SER316 on chain H, and the ASN312 residue on chain I. In addition, eleven polar and four hydrophobic residues were proximal to darobactin at a cut-off distance of 3 Å. There were only two negatively charged residues, ASP314 on chain G and ASP314 on chain H, and no positive residues were involved in the interaction. Urechistachykinin II had a binding score of −9.332 kcal/mol and formed six hydrogen bonds. One π-π stack interaction, 13 polar residues, seven hydrophobic bonds, two negatively charged residues, and two positively charged residues were involved in the interaction of urechistachykinin II with the protein. Plantaricin JLA-9 had a binding affinity of −9.002 kcal/mol and formed six hydrogen bonds. One π-π stack interaction, one π-cation interaction, nine polar residues, four hydrophobic bonds, and two negatively charged residues were involved in the interaction of plantaricin JLA-9 with the protein.
All the screened AMPs and MRL-494 interacted well with all three sub-chains of TolC simultaneously. Similar to MRL-494, darobactin, urechistachykinin, and planataricin JLA-9 formed H bonds with one or more of the residues among ASP314 on chain G and ASP314 and ASP315 on chain H, suggesting that these three peptides bind to the same domain as MRL-494. Table 3. AMPs screened by APD3 and the Standard Precision (SP)-peptide mode docking score and interactions with the receptor. The selection of residue groups located within 3 Å of the AMP ligand.

Definition
Docking Score Glide Model H-Bond π-π Stack π-Cation Since the three AMPs, darobactin, urechistachykinin II, and plantaricin JLA-9, had a good binding affinity and close interaction with TolC, these three AMPs were selected for further kinetic simulation analyses. In Figures 3 and S2-S5, the interaction diagrams of the three AMPs mentioned above showed that the docking position was at the outlet of the efflux pump. This interaction reduced the surface area at the efflux pump's outlet to inhibit the protein's efflux function.

Molecular Dynamics Simulations
To understand the structural and conformational changes during peptide binding, 100 ns molecular dynamics simulations were performed for five complex systems (apoprotein, a known inhibitor, and three selected AMPs). The root-mean-square deviation of atomic positions (RMSD) and root mean square fluctuation (RMSF) were monitored to evaluate the stability of the protein and peptide complexes using the simulation trajectory. In Figure 4, the protein Cα RMSDs were compared, and the peptide-bound complexes showed a lower deviation. The average RMSD value for the overall simulation and its standard deviation of the apo-protein was 4.99 ± 0.487 Å. The average RMSD value for the overall simulation and its standard deviation of the positive control MRL-494 was 4.38 ± 0.338 Å. The average RMSD value for the overall simulations and their standard deviation of darobactin, plantaricin JLA-9, and urechistachykinin II were 4.83 ± 0.510 Å, 4.72 ± 0.407 Å, and 4.11 ± 0.258 Å, respectively. Based on the RMSD values, all these simulations were considered stable after 50 ns.

Molecular Dynamics Simulations
To understand the structural and conformational changes during peptide binding, 100 ns molecular dynamics simulations were performed for five complex systems (apoprotein, a known inhibitor, and three selected AMPs). The root-mean-square deviation of atomic positions (RMSD) and root mean square fluctuation (RMSF) were monitored to evaluate the stability of the protein and peptide complexes using the simulation trajectory. In Figure 4, the protein Cα RMSDs were compared, and the peptide-bound complexes showed a lower deviation. The average RMSD value for the overall simulation and its standard deviation of the apo-protein was 4.99 ± 0.487 Å. The average RMSD value for the overall simulation and its standard deviation of the positive control MRL-494 was 4.38 ± 0.338 Å. The average RMSD value for the overall simulations and their standard deviation of darobactin, plantaricin JLA-9, and urechistachykinin II were 4.83 ± 0.510 Å, 4.72 ± 0.407 Å, and 4.11 ± 0.258 Å, respectively. Based on the RMSD values, all these simulations were considered stable after 50 ns.  In Figure 5, Cα RMSF was used to measure the flexibility and stability of the protein residues. RMSF maps were prepared for each of the three sub-chains of the trimeric efflux pump protein to view the contribution of the ligand to the residues. The residues with high RMSF values were always considered highly flexible regions. The average RMSF at the active binding region of the apo-protein was 10.4 ± 1.53 Å . The average RMSF of MRL-494 was 3.56 ± 0.513 Å . The average RMSF of darobactin, plantaricin JLA-9, and urechistachykinin II were 5.64 ± 0.810 Å , 5.51 ± 0.860 Å , and 4.62 ± 0.775 Å , respectively. All the AMP-protein complexes had lower RMSF values compared to the apo-protein, suggesting that AMP binding can form stable protein-AMP complexes.  In Figure 5, Cα RMSF was used to measure the flexibility and stability of the protein residues. RMSF maps were prepared for each of the three sub-chains of the trimeric efflux pump protein to view the contribution of the ligand to the residues. The residues with high RMSF values were always considered highly flexible regions. The average RMSF at the active binding region of the apo-protein was 10.4 ± 1.53 Å. The average RMSF of MRL-494 was 3.56 ± 0.513 Å. The average RMSF of darobactin, plantaricin JLA-9, and urechistachykinin II were 5.64 ± 0.810 Å, 5.51 ± 0.860 Å, and 4.62 ± 0.775 Å, respectively. All the AMP-protein complexes had lower RMSF values compared to the apo-protein, suggesting that AMP binding can form stable protein-AMP complexes.  In Figure 5, Cα RMSF was used to measure the flexibility and stability of the protein residues. RMSF maps were prepared for each of the three sub-chains of the trimeric efflux pump protein to view the contribution of the ligand to the residues. The residues with high RMSF values were always considered highly flexible regions. The average RMSF at the active binding region of the apo-protein was 10.4 ± 1.53 Å. The average RMSF of MRL-494 was 3.56 ± 0.513 Å. The average RMSF of darobactin, plantaricin JLA-9, and urechistachykinin II were 5.64 ± 0.810 Å, 5.51 ± 0.860 Å, and 4.62 ± 0.775 Å, respectively. All the AMP-protein complexes had lower RMSF values compared to the apo-protein, suggesting that AMP binding can form stable protein-AMP complexes.

Prime MM-GBSA Analysis
The molecular mechanics-generalized Born surface area (MM-GBSA) binding energy method was performed to calculate the binding energies between the AMPs and the TolC protein quantitatively. Table 4 presents the overall average binding free energy (∆G bind ) and the specific binding free energy components. All the listed energies were obtained and calculated from the last 50 ns of the coordinate sampling of the simulated trajectory. The ∆G bind of MRL-494 was −63.55 ± 8.51 kcal/mol. The ∆G bind of darobactin, plantaricin JLA-9, and urechistachykinin II were −42.09 ± 11.20 kcal/mol, −51.84 ± 10.90 kcal/mol, and −46.15 ± 12.09 kcal/mol, respectively. The ligand strain energy represents the energy cost of getting the compound into the active pocket. This is the energy difference between the optimized drug-protein complex and the drug outside the binding pocket. The ligand strain energies for the screened peptides, in increasing order, were plantaricin JLA-9, darobactin, and urechistachykinin II. Lower ligand strain energies are associated with a ligand more easily accommodating and fitting to its associated active site.

Principal Component Analysis (PCA)
To understand the motion of the simulation, a PCA was performed using the pairwise distance method. Five hundred frames were sampled to represent the changes across the simulation ( Figure 6A-E). The first two eigenvectors, PC1 and PC2, in all five PCA tests contributed to approximately 50% of the variance to the motion. During the MD simulation, the conformational changes of the docked peptide ligands and protein receptors were represented by a color gradient from purple to yellow. Figure 6F shows that the first few eigenvectors explained most of the variance, and when the number of eigenvectors reached 100, it explained practically 90% of the total variance. The PC1 vs. PC2 scatter plots explained the structural changes by the AMP interactions. The graphs showed consecutive point changes, and no outliers were found, indicating a smooth simulation. The PC1 of apo-protein exhibited a higher range (~300) than the interacting complexes (~200). This phenomenon also occurred on the apo-protein's PC2. The PCA's clustering of the data points suggested structural stabilization, whereas scattering suggested structural changes. According to the first and second principal components, the interacting complexes had low structural variation. According to the last 50 ns of the simulation, all the AMP groups underwent a point of cluster distribution, which was thought to reduce protein flexibility due to binding, thus limiting the variability of protein conformation. These AMP-protein structures also hindered the cylindrical outlet, reducing or even blocking the enforced excretion of intracellular drugs.

Biological Efficacy Assays
Because CLas is unculturable, the biological efficacy of the three in silico predicted peptides was assessed using two closely related surrogates, a culturable Rhizobium galegae as well as a moderately culturable Liberibacter crescens. Such surrogate models are related to CLas and have been used for the potential drug screening and testing of uncultured Liberibacter [9,13]. For example, a surrogate Sinorhizobium meliloti was used to study the effects of CLas transcription regulator proteins LdtR, RpoH/RpoH1, and VisNR [57][58][59]. CLas SecA protein was studied in silico and in vitro for potential small molecule inhibitors using L. crescens. L. crescens was used to study the mode of action of the transcriptional accessory protein PrbP in interactions with tolfenamic acid [60,61]. In this study, the surrogate bacteria R. galegae and L. crescens were closely related to CLas, and in the same order as Rhizobiales.
According to the Tblastn sequencing analysis, the CLas efflux pump protein aligned well with corresponding homologs in the culturable surrogates [62]. In R. galegae, the aligned segments were in nucleotides 4,025,395 to 4,026,762 and 2,135,053 to 2,135,241. The first segment covered 98% of the CLas TolC protein, and the second covered 16% of the CLas TolC protein. The first segment aligned with CLas. The identity rate was 42%. The max score was 382, with an E value of 9 × 10 −120 In L. crescens, the aligned segment was in nucleotides 1,077,242 to 1,078,594. The aligned sequences showed 98% coverage of CLas proteins and 48% identity. The max score was 416, with an E value of 9 × 10 −132 . The results indicated that the surrogates had similar functional proteins associated with the CLas TolC efflux pump protein.
The efficacy of the peptides on R. galegae was examined by quantifying the growth of bacteria, i.e., colony forming units (CFU) upon treatment with the AMPs. In these assays, Plantaricin JLA-9 significantly inhibited R. galegae (p ≤ 0.05) in all the tested dosages (10, 25, 50, and 100 µg/mL) when compared to the mock control. Darobactin significantly (p ≤ 0.05) inhibited R. galegae at 25, 50, and 100 µg/mL, whereas Urechistachykinin II showed efficacy only at 100 µg/mL ( Figure 7A). In addition, a viability/mortality assay was performed using the viability/cytotoxicity assay kit (Biotium) that uses two fluorescent dyes to detect dead bacteria as well as all bacteria. R. galegae were treated with the three peptides at a dosage of 50 µg/mL and cell mortality rates were estimated by counting the

Biological Efficacy Assays
Because CLas is unculturable, the biological efficacy of the three in silico predicted peptides was assessed using two closely related surrogates, a culturable Rhizobium galegae as well as a moderately culturable Liberibacter crescens. Such surrogate models are related to CLas and have been used for the potential drug screening and testing of uncultured Liberibacter [9,13]. For example, a surrogate Sinorhizobium meliloti was used to study the effects of CLas transcription regulator proteins LdtR, RpoH/RpoH1, and VisNR [57][58][59]. CLas SecA protein was studied in silico and in vitro for potential small molecule inhibitors using L. crescens. L. crescens was used to study the mode of action of the transcriptional accessory protein PrbP in interactions with tolfenamic acid [60,61]. In this study, the surrogate bacteria R. galegae and L. crescens were closely related to CLas, and in the same order as Rhizobiales.
According to the Tblastn sequencing analysis, the CLas efflux pump protein aligned well with corresponding homologs in the culturable surrogates [62]. In R. galegae, the aligned segments were in nucleotides 4,025,395 to 4,026,762 and 2,135,053 to 2,135,241. The first segment covered 98% of the CLas TolC protein, and the second covered 16% of the CLas TolC protein. The first segment aligned with CLas. The identity rate was 42%. The max score was 382, with an E value of 9 × 10 −120 In L. crescens, the aligned segment was in nucleotides 1,077,242 to 1,078,594. The aligned sequences showed 98% coverage of CLas proteins and 48% identity. The max score was 416, with an E value of 9 × 10 −132 . The results indicated that the surrogates had similar functional proteins associated with the CLas TolC efflux pump protein.
The efficacy of the peptides on R. galegae was examined by quantifying the growth of bacteria, i.e., colony forming units (CFU) upon treatment with the AMPs. In these assays, Plantaricin JLA-9 significantly inhibited R. galegae (p ≤ 0.05) in all the tested dosages (10,25,50, and 100 µg/mL) when compared to the mock control. Darobactin significantly (p ≤ 0.05) inhibited R. galegae at 25, 50, and 100 µg/mL, whereas Urechistachykinin II showed efficacy only at 100 µg/mL ( Figure 7A). In addition, a viability/mortality assay was performed using the viability/cytotoxicity assay kit (Biotium) that uses two fluorescent dyes to detect dead bacteria as well as all bacteria. R. galegae were treated with the three peptides at a dosage of 50 µg/mL and cell mortality rates were estimated by counting the number of dead cells (red stained) among all cells (green stained) ( Figure 7B). All three peptides, Plantaricin JLA-9, Darobactin, and Urechistachykinin II caused significant mortality of 55.50, 52.20, and 53.60%, respectively, thus confirming the antimicrobial activity against R. galegae ( Figure 7B). Similarly, efficacy of the three peptides was assayed on L. crescens using the viability/cytotoxicity assay. Plantaricin JLA-9 and Urechistachykinin II showed significant cell mortality (p ≤ 0.01); however, Darobactin did not impact L. crescens viability at the concentrations tested, suggesting bacteria-specific differences (Supplemental Figure S6). number of dead cells (red stained) among all cells (green stained) ( Figure 7B). All three peptides, Plantaricin JLA-9, Darobactin, and Urechistachykinin II caused significant mortality of 55.50, 52.20, and 53.60%, respectively, thus confirming the antimicrobial activity against R. galegae ( Figure 7B). Similarly, efficacy of the three peptides was assayed on L. crescens using the viability/cytotoxicity assay. Plantaricin JLA-9 and Urechistachykinin II showed significant cell mortality (p ≤ 0.01); however, Darobactin did not impact L. crescens viability at the concentrations tested, suggesting bacteria-specific differences (Supplemental Figure S6). 4% v/v) used for dissolving the peptides as a mock control. Asterisks indicate statistically significant differences between control and treatment by Student's t-test (*, 0.01 and **, 0.001). (B) R. galegae cell viability/mortality was estimated after treatment with AMPs (50 µg/mL) and along with a mock (DMSO, 2.4% v/v). Boiled cells were used as a positive control. Percent (%) mortality was estimated by counting the number of dead cells (red stained) among all cells (green stained). Error bars represent ± standard error of the mean (n = 5).
p-values were calculated by a two-sample t-test (one-tailed) relative to mock control. Scale bar = 10 µm.
In summary, the in vitro bioassays indicated that Plantaricin JLA-9, Darobactin, and Urechistachykinin II showed inhibitory activities in one or both surrogate bacteria. Plantaricin JLA-9 showed significant inhibition of R. galegae as low as 10, 25, 50, and 100 µg/mL and L. crescens at 50 and 100 µg/mL. Plantaricin JLA-9, initially isolated from Suan-Tsai, showed broad-spectrum antibacterial activity against Gram-positive and Gram-negative bacteria [52]. Darobactin showed significant inhibition of R. galegae at 25, 50, and 100 µg/mL. Darobactin is a bicyclic peptide discovered in Photorhabdus bacteria in 2019 [55]. In this study, darobactin was synthesized as a linear peptide; thus, the bioactivities could be different. Urechistachykinin II showed inhibition of R. galegae at 50 and 100 µg/mL, whereas L. crescens showed mortality at 100 µg/mL. Urechistachykinin II is a tachykinin-related neuropeptide from the echiuroid worm with bioactivity on neurons and G protein-coupled receptors in animals [49]. In other studies, Blacutt identified three active natural compounds, cladosporols A, C, and D, using the L. crescens diffusion assay [63]. These three compounds were found to have dose-dependent antagonisms to L. crescens. Zhang screened two small bioactive molecules (P684-2850, P684-3808) with strong antimicrobial activities to L. crescens in vitro [60]. The IC50 values for P684-2850 and P684-3808 were 11 µg/mL and 15 µg/mL, respectively. Although our surrogate bioassays verified the bioactivities of the in silico screened molecules, further long-term structure-function relationships and in planta efficacy trials are necessary to understand the inhibitory mechanism on CLas.

Sequence Analysis
Computational predictions of the three-dimensional structures of the proteins were constructed using comparative homology modeling techniques. The CLas TolC efflux pump protein was identified by the NCBI unique identifier CLIBASIA_04145 in CLas [2]. Based on the SWISS-MODEL server [64], the crystal structure of the efflux pump component OprN from Pseudomonas aeruginosa (PDB: 5iuy) [40] was used as a template to construct a structural model of the CLas TolC efflux protein. The Ramachandran plots were visualized by PROCHECK [42] to evaluate the effectiveness of the homology model(s). The overall quality of the homology models was evaluated by ERRAT [65], VERIFY 3D [66,67], and PROVE [68] on the SAVES v6.0 server (https://saves.mbi.ucla.edu accessed on 1 October 2022).

Virtual Screening of AMPs
The initial set of AMPs was screened from the Antimicrobial Peptide Database (APD3) [69]. Several criteria were considered for screening AMPs targeting the membrane protein tripartite efflux pump in CLas. First, the AMPs should be active against Gram-negative bacteria. Second, the sequence length of the AMPs was around 15 amino acids to allow stability and bioavailability. Additionally, AMPs with known membrane protein activities were also included. All these parameters were applied in the search engine in APD3.

Molecular Docking
The binding mode prediction of AMPs with the CLas efflux pump receptor complex was performed on the Schrödinger Glide platform [70]. Both the protein and AMPs were prepared by Schrödinger's Protein Preparation Wizard [71]. The AMPs were docked on the rigid protein receptor to find the potential binding pockets with the highest docking scores. The ligand-receptor interactions were visualized using Schrödinger Maestro [72], and the AMPs with multiple and stable non-covalent interactions were selected for molecular dynamics (MD) simulations. MRL-494, a small-molecule inhibitor of β-barrel assembly machine A, was used as the positive control [56].

Molecular Dynamics Simulations
Promising AMPs from the initial docking studies were selected to elucidate the possible mode of action (MOA) in terms of protein-peptide binding interactions [73]. The MD simulation was performed on the Schrödinger Desmond platform [73]. The system charge was neutralized with Na or Cl ions and a 0.1 M concentration of NaCl. The system was built with the Simple Point-Charge (SPC) solvent model. The boundary conditions were set as an orthorhombic box with a buffer distance of 10 Å. The model membrane of the system was the bilayer 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) membrane. The protein-ligand complex was prepared by Schrödinger's Protein Preparation Wizard [71]. The docked AMPs were exported from Glide. The C and N termini of the protein were capped to stabilize the protein structure. All the missing hydrogens were added, and the hydrogen bonds were optimized. The strained minimization was performed with the OPLS3e force field [74].
Each simulation was performed with Desmond's default relax protocol. The force field for the simulation was also OPLS3e. The system-heavy atoms were first minimized with restraints under 10 K, then increased in temperature to 300 K with restraints, and the final relaxation step under 300 K Normal Pressure and Temperature (NPT) ensemble to obtain the equilibrium status for the system. After relaxation, the simulations were performed under a 300 K NPT ensemble at 1.01325 bar pressure for 100 ns. The recording interval was 200 ps, and 500 frames were saved. Post-simulation trajectory analyses, including complex (RMSD) and ligand/protein (RMSF) complex interactions were performed by Schrödinger Simulation Interactions Diagram.

Calculation of Binding Free Energy
The binding interactions of peptide-protein complexes can be calculated using (MM-GBSA) binding energy method [75]. The primary MM-GBSA uses the VSGB 2.0 dissolution model with the OPLS3e force field.
The docking poses of the AMPs on the receptor were evaluated by calculating the total binding free energy using Schrödinger Prime. Prime MM-GBSA calculations gave the complex binding energies, which validated the performance of the current binding conformation. The Schrödinger Prime calculates the energy of the AMP-protein system via MM-GBSA using the Desmond simulation trajectory. From the entire 100 ns simulation, the last 50 ns trajectory was chosen for the energy calculation. The free energy of binding ∆G bind was calculated as the energy of the receptor-ligand complex minus the energy of the receptor alone and the ligand alone, as follows:

Principal Component Analysis
To study the overall motion of the protein system throughout the entire simulation trajectory, a principal component analysis (PCA) for a pairwise distance of alpha carbon atoms (Cα) was performed. Visual Molecular Dynamics (VMD) [76] software was used to prepare the topology and trajectory files. The python package MDTraj [77] was used for analysis and visualization. Cα was extracted from 500 snapshots of the 50 ns MD trajectories to construct the structure matrix. Every snapshot was aligned to the starting frame before constructing the eigenvectors. The first two principal components (PC1 and PC2) were projected to represent the protein dynamics in the MD trajectories of the AMP-protein and the apo-protein systems.

Peptide Synthesis
The three in silico predicted peptides (Plantaricin JLA-9, Darobactin, and Urechistachykinin II) were synthesized using PeptideSyn TM Peptide Synthesis Technology (LifeTein, LLC, Somerset, NJ, USA). The purity of the peptides was >95%, as determined by High-Performance Liquid Chromatography (HPLC) and Mass Spectrometry (MS) (Supplementary Dataset). All the peptides were dissolved in 100% DMSO and diluted further for in vitro efficacy testing.

Biological Efficacy Assays
Because CLas is unculturable, the biological efficacy of the three peptides was determined in vitro using closely related surrogates, Rhizobium galegae and Liberibacter crescens. The sequence analysis and comparison of the CLas TolC efflux protein CLIBASIA_04145 and surrogate hosts R. galegae and L. crescens were performed using Tblastn [78]. For R. galegae assays, cells were grown at 28 • C in YM broth (Yeast extract 3 g/L, Malt extract 3 g/L, Dextrose 10 g/L, Peptone 5 g/L, pH 6.2) to a cell density (OD 600 ) of 1.0. Next, the R. galegae cells were diluted into a fresh YM media to a cell density of 1 × 10 6 . Multiple concentrations of the three peptides (10, 25, 50, and 100 µg/mL) were mixed in a 96-well plate (Nunc F96 microtiter plate) with a bacterial suspension (at a final concentration of 1 × 10 6 colony forming units [CFU] per ml for bacteria) to a total volume of 100 µL. Three biological replicates for each concentration and peptide were used. Because pure DMSO was used to dissolve the peptides, a mock control containing equivalent concentrations of final DMSO concentration (0.2, 0.6, 1.2, and 2.4% v/v), corresponding to the different peptide dosage (10, 25, 50, and 100 µg/mL), was also included. The cells were treated by incubation at 28 • C for 2 h, and then 50 µL was plated on YM agar plates. Subsequently, the plates were incubated at 28 • C for 48 h and CFUs were estimated.
L. crescens is still a moderately fastidious bacterium and generally forms a lawn instead of colonies. Hence, we conducted a cell viability/mortality assay using dyes that distinguishes between live and dead cells (Biotium, Fremont, CA, USA). Selected dosages of AMPs (50 and 100 µg/mL) were tested on L. crescens as well as R. galegae. Briefly,~3 mL of L. crescens culture was grown in a BM7 (α-Ketoglutarate 2 g/L, ACES buffer 10 g/L, KOH 3.75 g/L, Fetal Bovine Serum 15%, TMN-FH Medium 30%, pH 6.8) medium for 4 days. The cells were harvested by centrifuging 1 mL of cultures at 7000× g. The cells were washed three times with 0.85% NaCl and diluted to OD600 of 0.1 with 0.85% NaCl. The R. galegae cells were grown overnight at 28 • C, 200 rpm in a 3 mL TY medium. All cells were treated further with antimicrobial peptides at 50 or 100 µg/mL for 3 h at 28 • C in a plate shaker at 280 rpm. Cells were collected from the microplate into a 1.5 mL tube and centrifuged at 10,000× g for 10 min, the supernatant was discarded, and cells were resuspended with 100 µL of 0.85% NaCl. As a control treatment, cells were collected and boiled for 5 and 10 min for R. galegae and L.crescens, respectively. After treatment, 1 µL of dye mix was added to each treatment and incubated at room temperature for 15 min in the dark. Cells were imaged using an Olympus BX51 (Olympus LifeScience, Waltham, MA, USA) fluorescence microscope under a 60× (Immersion oil) objective using FITC and TRITC filter sets. Representative images were taken, and the percent mortality was calculated from an average of 6 field of views from each treatment.

Conclusions
In this study, we screened AMPs as potential inhibitors for CLas by blocking TolC, a key protein in the efflux pump of the T1SS system. Multiple in silico approaches, such as homology modeling, molecular docking, MD simulations, MM-GBSA calculations, and PCA were used to uncover potential AMPs that target the outer membrane protein TolC. Initial docking studies revealed three AMPs with good binding affinities to TolC. MD simulations for 100 ns confirmed the ability of the AMPs to bind tightly and interact with TolC receptors. Based on further PCA and Gibbs free energy calculations, plantaricin JLA-9 showed the ability to block the β barrel entrance of the TolC protein. The in vitro bioassay indicated that Plantaricin JLA-9, Darobactin, and Urechistachykinin II showed inhibitory activities in one or both surrogate bacteria. Plantaricin JLA-9 showed significant inhibition even at a lower concentration to R. galegae at 10 µg/mL; darobactin showed significant inhibition of R. galegae at 25, 50, and 100 µg/mL, and urechistachykinin II showed inhibition at 50 and 100 µg/mL. In the case of L. crescens, the Plantaricin JLA-9 and Urechistachykinin II showed significant mortality at 50 and 100 µg/mL. In conclusion, the strategy and results from this study could help design effective therapies to manage plant diseases caused by Candidatus Liberibacter spp.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27248729/s1, Figure S1. Homology modeling results. Figure S2. Interactions between MRL-494 (positive control) and receptors under SP-peptide docking mode. Figure S3. Interactions between darobactin and receptors under SP-peptide docking mode. Figure S4. Interactions between plantaricin JLA-9 and receptors under SP-peptide docking mode. Figure S5. Interactions between urechistachykinin II and receptors under SP-peptide docking mode. Figure S6. Efficacy of antimicrobial peptides against L. crescens viability and its mortality. Table S1. The detailed information of the screened AMPs. And AMP purity by HPLC in supplementary.pdf.