Introduction

The vector-borne disease leishmaniasis is caused by various species of Leishmania, a eukaryotic protozoan parasite. The disease is transmitted by female phlebotamine sandflies, which feed on the blood of humans and other mammals. A World Health Organization (WHO) report from 2021 indicates that 79 countries were affected by visceral leishmaniasis (VL). Around 87% of global cases of VL were reported from eight countries: Brazil, Eritrea, Ethiopia, India, Kenya, Somalia, South Sudan, and Sudan (https://www.who.int/data/gho/data/themes/topics/topic-details/GHO/leishmaniasis). A population of more than 1 billion people lives in areas prone to leishmaniasis. VL and cutaneous leishmaniasis (CL) become increasingly common each year (Tzani et al. 2021). It is estimated that 30,000 new cases of VL and 1 million new cases of CL occur every year (https://www.who.int/news-room/fact-sheets/detail/leishmaniasis). International travel and climate change contribute to the increase in leishmaniasis cases in the past (Purse et al. 2017). However, despite a reduction in international travel due to COVID-19, leishmaniasis cases continue to increase (Le Rutte et al. 2021). Therefore, the reduction of leishmaniasis cases is of great concern. The most common drugs used for treating leishmaniasis are pentavalent antimonial, amphotericin B, miltefosine, paromomycin, and pentamidine (Chakravarty and Sundar 2019; Singh and Kumar 2014). However, these drug treatments suffer from one or more limitations, such as treatment failure (Ponte-Sucre et al. 2017), resistance (Horacio et al. 2021), longer half-life (Reimao et al. 2020), relapses after treatment (Araújo et al. 2021), toxicity (Arenas et al. 2017), or higher cost (Azim et al. 2021). Considering the absence of a vaccine, a number of alternative strategies are being researched from a range of significant pathways to create a new drug for the disease (Ranjan et al. 2020; Vijayakumar et al. 2019; Vijayakumar and Das 2018). Some common pathways being targeted include, but are not limited to, purine salvage pathway (Brindha et al. 2021), ergosterol biosynthesis pathway (McCall et al. 2015), polyamine biosynthesis pathway (Da Silva et al. 2012), and thiol metabolism pathway (Singh et al. 2017). Among many important pathways, trypanothione metabolism pathway is considered to play a crucial role in maintaining intracellular redox balance and protecting cells from oxidative stress (Battista et al. 2020). However, the synthesis of glutathione and trypanothione in the pathway is primarily dependent on cysteine availability. In Leishmania, there are two routes by which cysteine can be generated (Williams et al. 2009). The first is the de novo biosynthesis pathway, or sulfhydrylase pathway, and the second is the reverse trans-sulfuration pathway, which converts homocysteine to cysteine in two steps (Fig. 1). Due to the absence of the de novo biosynthesis pathway in mammals, enzymes involved in these pathways are viewed as potential drug targets. De novo cysteine biosynthesis pathway involves two enzymes: serine sulfhydrase and cysteine synthase, which requires sulfide. While plants and bacteria utilize sulfide via the sulfide reduction pathway (Gotor et al. 2015; Pellerin et al. 2019), Leishmania lacks the gene encoding the enzymes of this pathway (Williams et al. 2009), indicating that the organism uses a different source of sulfide. A study (Williams et al. 2003) demonstrated that the major source of sulfide for Leishmania is provided by the action of 3-mercaptopyruvaate sulfurtransferase (3MST) on 3-mercapto pyruvate (3MP). 3MST belongs to the sulfurtransferase family of enzymes, which catalyzes the transfer of sulfur from acceptor molecules to form product. Sulfurtansferases are widely present in plants, animals, and bacteria and have a wide range of roles (Marciano et al. 2010; Nakamura et al. 2000). Among their functions are detoxification, the formation of an iron-sulfur cluster, and the synthesis of 4-thiouridine and 4-thiamine (Mueller et al. 2001; Palenchar et al. 2000). However, in Leishmania evidence suggest that 3MST plays some role in antioxidant and sulfur amino acid metabolism (Singh et al. 2017). It has also been observed that the rate of utility of sulfide arising from the activity of 3MST and when sulfide itself is provided to form cysteine were similar (Williams et al. 2009), indicating the 3MST as a major source of sulfide for the parasite. Further, the enzyme 3MST has been reported to be expressed in all stages of Leishmania life cycle and its expression is upregulated during peroxide or sulfur stress (Williams et al. 2003). This indicates that 3MST plays a role in parasite survival. As well, the structure of Leishmania 3MST (PDB ID. 1OKG) differs significantly from that of the human equivalent (PDB ID 4JGT). As compared to mammalian 3MST, Leishmania 3MST contains 70 additional amino acids in the C-terminal domain (Alphey et al. 2003). Moreover, these additional residues are conserved among trypanosomatids parasite, which signifies their importance. The presence of these additional residues has been reported to be necessary for the folding of the 3MST (Alphey et al. 2003). The above facts make 3MST an appealing drug target. Despite this, as far as the authors are aware, no studies have been performed that attempted to screen inhibitors against Leishmania 3MST (especially L. donovani) in the literature. Hence, in this study, we explored a wide range of natural compounds against the L. donovani 3MST (Ld3MST) through modeling and in silico docking, with the optimal lead-protein complex further validated via a long-term MD simulation. In addition, the in silico validated compound was evaluated for in vitro toxicity and drug efficacy. As a result of our findings, we believe that Rutin, with its non-toxic properties at high concentrations as well as its leishmaniacidal properties at low concentrations, will be a promising candidate for further investigation as a potential drug.

Fig. 1
figure 1

Cysteine biosynthesis pathway in Leishmania sp

Full size image

Materials and methods

Homology modeling

Due to the lack of the Ld3MST X-ray/NMR structure, a homology model of the same is built. The sequence of the enzyme Ld3MST has been obtained from the National Center for Biological Information (accession ID XP_003858284). By using the BLASTp tool, a local alignment search (BLAST) was conducted with the Ld3MST sequence in the protein data bank (PDB) to identify a suitable template. Those templates with maximum identity and maximum query coverage were considered. To model the Ld3MST, the open source modeling tool Modeler V9.0 was used (Webb and Sali 2016). Five models were generated and the structural quality of each model was evaluated using ERRAT (Colovos and Yeates 1993), VERYFY3D (Eisenberg et al. 1992), and the Ramachandran Plot (Laskowski et al. 1993) using SAVES V6.0 (https://saves.mbi.ucla.edu/). Validation measures equaling or better than the template model were deemed the best models to proceed with further analysis.

Molecular docking

Molecular docking was performed on the validated 3D structure of Ld3MST. The study utilized TIPdb (Lin et al. 2013), a natural compound library that contains 5284 drug-like chemicals. The structure of all drug-like chemicals was downloaded from TIPdb (https://cwtung.kmu.edu.tw/tipdb/) in SDF format. For molecular docking, the Genetic Optimization for Ligand Docking (GOLD) software was used (Verdonk et al. 2003). Both ligands and proteins were prepared using the GOLD “preparation wizard”. Active site residues for docking were obtained from the previous study (Alphey et al. 2003). Between the N-terminal domain and central domain, the active site contains residues from seven segments of the polypeptide viz., 35–39, 72–75, 104–109, 180–185, 193–196, 209–221, and 253–258 (Fig. 2). However, the active site floor is formed by a loop, with a catalytic side chain of cysteine 253 (CYS253) in the center (Fig. 2). CYS253 C-alpha coordinates were thus chosen for study, and a radius of 10 Angstroms (Å) was taken as a distance around it to describe the cavity for docking, which by definition includes residues from all seven segments of polypeptide (Alphey et al. 2003). For each ligand, ten dock poses were generated. 3MP is a natural substrate for the enzyme 3MST (Williams et al. 2003) and was used as a positive control. Dock score of 3MP was used as the cutoff score, and only ligands scoring above the cutoff were considered. LIGPLOT+ (Laskowski and Swindells 2011) was used to study the interactions between the top scored ligand complex and substrate complex. Molecular dynamics simulation (MD) was performed to examine the stability of the complexes.

Fig. 2
figure 2

Structure of Ld3MST. Pink: 3-MP; blue: active site helix; green: active site sheet; brown: active site loops; and yellow: active site center residue cysteine 253

Full size image

Molecular dynamic simulation (MD)

The MD simulations and energy minimization steps were carried out using GROningen MAchine for Chemical Simulations (GROMACS) tool. GROMOS96 43a1 force field (Schuler et al. 2001) was utilized and in order to produce protein topography files, GROMACS built-in utilities were used. Topology files of ligands were generated by the PRODRG server (Van Aalten et al. 1996). Initially, the protein–ligand complexes were placed in the periodic boundary conditions, using the 10 Å distance between the protein and box edges, and each system was solvated independently using the single point charging solvent model. A required number of sodium ions (Na +) or chlorine ions (Cl −) were added, in order to balance the total charge in the system. The energy minimization for each system was achieved by using the steepest descent algorithm for 50,000 steps. To constrain all types of bond angles, the linear constraint solver (LINCS) algorithm (Hess et al. 1997) was used, and to treat long-range electrostatics, the Particle Mesh Ewald method (Essmann et al. 1995) was used. In addition, equilibration was conducted in two phases. During the first phase, the number of particles, volume, and temperature (NVT) were kept constant and one nanosecond (ns) of NVT ensemble was run using the V-rescale temperature coupling method to attain 300 K. In the second phase, the number of the particles, pressure, and temperature (NPT) were kept constant and 1 ns of NPT ensemble was run using Parrinello–Rahman pressure coupling method (Martoňák et al. 2003) to attain 1.0 bar pressure. Production MD was run for 100 ns after equilibration of the system at defined temperature and pressure. The simulation was conducted with three systems, i.e., Ld3MST alone (unbound), Ld3MST with 3MP (substrate), and Ld3MST with the best docked natural compound (Rutin). The corresponding MD trajectories were analyzed to study the root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), radius of gyration (Rg), hydrogen bond interactions, and binding energy. For binding energy calculation, PRODIGY-LIGAND (Vangone et al. 2019) was utilized.

Macrophage cell line and parasite cultures

In this study, the macrophage cell lines THP-1 and RAW-264.1 were used. The promastigote form of L. donovani (MHOM/IN/1983/AG83) parasite was used. Cells were cultured in RPMI-1640 supplemented with 10% fetal bovine serum (FBS), containing (0.1 mg/mL streptomycin and 100 IU/mL penicillin G) 1% antibiotic. In a CO2 incubator, RAW-264.1 cells were maintained at 37 °C in 5% CO2 (pH 7.4). Sub-culturing of the RAW-264.1 cell line was performed every 3 days. In the BOD incubator, promastigotes were cultured at 25 °C using M-199 medium with 10% FBS (pH 7.4).

In vitro assay on promastigotes and amastigotes

Leishmania donovani (MHOM/IN/1983/AG83) promastigotes were cultured as described above. Rutin was obtained from Sigma-ALDRICH® and as Rutin has high solubility in methanol (Kraujalis et al. 2013), 10 mg of Rutin was dissolved in 1 mL of methanol, to make a stock solution. Log phase cells (1 × 106/ml) were treated with Rutin (8.19, 16.38, 40.95, 81.90, 122.85, 163.80, and 327.60 μM) for 3 days and cell viability is measured by Trypan blue exclusion and methyl thiazole tetrazolium (MTT) reduction method after every 24 h using standard procedures (Tabrez et al. 2021). IC50 for the drug was defined as a concentration that killed 50% of the parasite after 48 h. The control studies were performed with Amphotercin B (AmB) at concentrations of 5, 10, 20, 30, and 50 nM and phosphate-buffered saline (PBS). The morphology of the cells was viewed using Nikon TS100-inverted microscope at the magnification of × 40.

Macrophage cells were maintained in RPMI-1640 medium and infected with promastigotes as described in the previous section. Parasite-infected macrophages, grown for 24 h, were treated with different concentrations of Rutin (16.38–131.04 µM), for 48 h. AmB (0.5 µM) was taken as positive control. Chamber slides were washed and supplemented with fresh medium and kept in CO2 incubator for additional 12 h. Untreated parasite-infected macrophages were used as control. The amastigotes from 100 macrophage nuclei per well were counted under the oil immersion objective of a light microscope (Nikon Eclipse TS100, Japan) after methanol fixation and Giemsa staining.

Cytotoxicity assays

In order to determine the cytotoxic effects, the THP-1 human macrophage cells (5 × 104 cells) were seeded in the presence of Rutin concentrations (40.95 to 819.00 μM) and AmB concentration (5 to 50 nM) in 96-well culture plates in 5% CO2 for 48 h at 37 °C. Cell viabilities were measured using colorimetric methyl thiazole tetrazolium (MTT) assay as it was described earlier (Tabrez et al. 2021). Absorbance ratios were quantified at O.D. 570 nm. PBS and AmB were utilized as controls. These results were indicated as the mean percentage (%) reduction of macrophages in comparison with untreated control samples × 100. Finally, the concentration producing 50% cytotoxicity (μM) was determined.

Cloning, characterization, and protein purification of recombinant Ld3MST

Leishmania donovani genomic DNA was isolated using a DNA isolation kit (Qiagen) according to the manufacturer’s instructions. Genomic DNA was used as a template to amplify 3MST gene using forward GTGTGGATCCATGTCTGCTCCCGCTGCTGC and reverse primer GTGTAAGCTTTTAGTTGGGCACTGCCGGCG with the help of polymerase chain reaction (PCR). Vector-pET28a was used to prepare the clone. The MST-PCR product and vector-pET28a were digested with BamHI and HindIII. In a 1:7 molar ratio, purified vector and insert were ligated by T4 DNA ligase. The ligated products were transformed into Escherichia coli DH5α and transformed colonies were selected on luria broth (LB) agar plates with antibiotics (30 mg/ml kanamycin). PCR was used to confirm the clone and the release of the insert was achieved by a double digestion using the same set of enzymes used for cloning. In addition, positive colonies were incubated overnight at 37 °C at 150 rpm in 10 ml LB-kanamycin (30 μg/ml) media. Sub-cultured cells (O.D. reached to 0.4–0.5) were treated with isopropyl β-d-1-thiogalactopyranoside (IPTG) (0.5 mM) incubated for 4 h; later, pellets were collected for SDS-PAGE analysis. In 200 μl of 2 × SDS gel loading buffer, the pellets were resuspended, boiled for 10 min at 100 °C, and centrifuged at 12,000 g for 1 min. An electrophoresis at 70 V constant current was performed with 10 μl of supernatant loaded on 10% gel. Five microliters of unstained protein marker was loaded on the same gel. Gel was stained with Coomassie brilliant blue R-250 staining solution for 1 h. The gel was destained with destaining solution and was observed under white light in the gel documentation system. The protein was purified using the Ni + NTA agarose bead system which binds specifically to histidine through the imidazole competitive binding mechanism.

Sulfurtransferase activity of recombinant 3-MST

The optimized reaction mixture (1 ml) contained 0.1 M Tris–HCl, pH 7.3, 10 mM Ferrocyanate, and recombinant enzyme 100 ng. The reaction was initiated by adding sodium thiosulfate at 5 mM and incubated for 10 min at 37 °C. The reaction was stopped by the addition of 500 μl of formaldehyde and 1.5 ml of ferric nitrate-50 g/L. Further, 65% HNO3 was added for FeSCN− complex formation, with sulphur coming from sodium thiosulphate and cyanide coming from ferro cyanide. As a result, this complex will have a red color and absorb 460 nm wavelength light (Williams et al. 2003). This assay system was then used to investigate Rutin enzyme activity inhibition. Different concentrations (1.64–40.95 μM) of Rutin were used to evaluate the enzyme activity.

Statistical methods

The data collected from in vitro, enzymatic, and cytotoxicity assays are analyzed with Prism 5.0 software (GraphPad Prism, San Diego, California, USA). The means of experimental groups were compared with the means of controls for statistical significance using the one-way analysis of variance and Dunnett’s multiple comparison test. In terms of p-value, statistical significance is denoted by “*” for p < 0.05, “**:” for p < 0.01, “***” for p < 0.001, and NS for non-significance. For this study, p < 0.05 is considered as significance. The other basic statistics and graphs on MD trajectory data were calculated using R Version 4.

Results and discussion

Structure of L. donovani 3MST

Due to the lack of X-ray or NMR structures, the structure of Ld3MST was built using homology modeling. To obtain a reliable model through homology modeling, a structural template having high identity with the query is a must. Hence for Ld3MST, template search was performed using BLASTp and opting PDB as database. The results had structure of L. major 3MST (PDB ID. 1OKG) as a best hit having a query coverage of 99%, identity of 95%, and with no gaps (Supplementary Fig. 1). Hence, homology modeling was performed with 1OKG as a template using Modeler. Five models were generated, and upon performing structural quality analysis using ERRAT, verify 3d, and Ramachandran plot (Supplementary Table 1), it has been observed that the modeled one had best quality score on all measures. The model and template structural quality measures were all compared, and the results (Supplementary Fig. 2−4) suggested that the model have on par better scores than that of template on every measures. More noticeably, the Ramachandran plot of the model had all the residues in favored/allowed region, except one residue in generously allowed region (Supplementary Fig. 2B). In addition, the ERRAT score (89.326) and verify 3d score (96.22) also found favoring the model (Supplementary Fig. 3b and 4b). Moreover, upon superimposing the model with template (Supplementary Fig. 5), it has been found that the RMSD was 0.56 Å. In general, a RMSD of less than 1 Å is considered good when the model and the template were superimposed (Modi et al. 2016). The above results suggest that the build model has above average quality and hence the built Ld3MST model was used for further analysis.

Molecular docking

The validated Ld3MST model was utilized for docking. Since this study involves screening of large number of natural compounds, a cutoff value for filtering best ligands become mandatory. In order to achieve this, initially, 3MP (natural substrate of Ld3MST) (Williams et al. 2003) was docked to the active site (as described in “Materials and methods” section) of Ld3MST. It was observed that the substrate 3MP had a GOLD score of 40.58 and prodigy binding score of − 6.0 kcal/mol (Table 1). In addition, 3MP docked pose was found interacting with key residues of the active site (Fig. 3A). Hence, in order to screen the library of natural compounds, a cutoff score of 40.58 was set, i.e., natural compounds that had a docking score of greater than or equal to 41.0 were alone considered. The 5284 drug-like natural chemicals from TIPdb were subjected to docking with Ld3MST. A total of 275 natural compounds were found having a dock score greater than the assigned cutoff. However, only five compounds had the dock score of greater than 50 and prodigy score of less than − 8.5 (Table 1). The top scoring compound (TIP000059) was found to be quercetin 3-rutinoside (Rutin), which is a flavanol glycoside reported to be found in more than 10 plants (Supplementary Table 2). Among the top five compounds, Rutin was found to have more H-bonds with the active site residues than others (Table 1). Also Rutin was found to have either hydrophobic interaction or hydrogen bonds with the majority of Ld3MST active site cavity residues (Fig. 3B). Moreover, the dock score (55.30) and binding energy score (− 10.5 kcal/mol) of Rutin were significantly greater than the dock score and binding energy score of the natural substrate 3MP. The observation from the docking result suggested that the screened compound Rutin has the potential to interact with the key active site residues of Ld3MST as well as possess better affinity towards the active site cavity in comparison to that of Ld3MST substrate.

Table 1 Docking solutions of screened ligands with L. donovani 3-MST
Full size table
Fig. 3
figure 3

Comparison of residues of Ld3MST interacting with A 3MP and B Rutin. Dotted lines: H-bond interactions. Spiked half circle: Hydrophobic interactions

Full size image

MD simulation of docked complexes

In order to observe the behavior of unbound Ld3MST (unbound), Ld3MST with 3MP (substrate), and Ld3MST with Rutin, a 100-ns MD simulation was performed for each of the system. Each system trajectory was analyzed and compared.

Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) analysis

The backbone RMSDs of every system were computed (Fig. 4) using 100-ns trajectory inputs. As a result, Ld3MST bound to substrate and Rutin shows greater stability than 3MST that is unbound. As can be seen from Fig. 4A1−A3, the Ld3MST-Rutin complex (Fig. 4A3) had one sharp peak. Thus, the Ld3MST-Rutin complex appears to be relatively stable and possess less deviations than the other two systems. Furthermore, the Ld3MST-Rutin complex had a lower average RMSD (~ 0.35 nm) than other systems during the entire simulation period. To observe whether the ligands are stable during simulation, we plotted the RMSD for all atoms of the substrate and Rutin (Fig. 4B). Results indicate that Rutin and substrate have lower overall deviations (< 0.25 nm). When observing how often the data points vary, it appears that Rutin maintains a relatively stable pose throughout the simulation in comparison to the substrate. Also, the residue-wise fluctuations of all three systems were computed (Fig. 4C). In contrast to the unbound and substrate-bound forms of Ld3MST, the residues of Rutin-bound complex fluctuated relatively less. While the RMSF differs between the systems, it appears that the patterns of fluctuation are similar between the Ld3MST-substrate and Ld3MST-Rutin complexes. This indicates that Rutin and natural substrates both impact the Ld3MST in a comparable way. Hence, based on the RMSD and RMSF analysis, it is obvious that the Ld3MST-Rutin complex is more stable compared to that of the unbound and substrate complex.

Fig. 4
figure 4

RMSD analysis of 100 ns trajectory. A Backbone RMSD of all systems. A1 RMSD density plot of unbound 3MST. A2 RMSD density plot of 3MST-3MP complex. A3 RMSD density plot of 3MST–Rutin complex. B All atoms RMSD of ligands. C RMSF of all residues of all the systems

Full size image

Radius of gyration (Rg) and solvent accessible surface area (SASA) analysis

Radius of gyration (Rg), a measure of the protein’s compactness (Hong and Lei 2009), was computed for all systems over the entire simulation period (Fig. 5). It has been observed that both Rutin and Substrate Ld3MST complexes Rg have converged over time and maintained stability better than the unbound Ld3MST. From Fig. 5A1−A3, it is evident that unbound Ld3MST was loosely packed compared to the other two complexes. The Ld3MST-substrate complex, however, was found to be marginally more compact than the Ld3MST-Rutin complex. Perhaps the marginal difference can be attributed to the smaller size of the substrate (11 atoms) in comparison to Rutin (73 atoms). For all systems, entire protein SASA was calculated to determine if loose packing influenced solvent accessibility (Fig. 6A). There were no significant differences in SASA observed among the simulated systems, since the average SASA value of the three systems was around ~ 181 nm2 (Fig. 6A1−A3, black broken lines). SASA calculations were also performed to determine if there had been any changes to the active site (Fig. 6B). As expected, the average SASA of unbound Ld3MST (~ 48 nm2) was relatively higher than that of substrate-bound (~ 44 nm2) and Rutin-bound Ld3MST (~ 44 nm2). In accordance with the Rg (Fig. 5A), SASA (Fig. 6A), and active site SASA (Fig. 6B), the Ld3MST-Rutin complex exhibits compactness and a more stable SASA throughout the simulation period than unbound and Ld3MST-substrate complex.

Fig. 5
figure 5

Radius of gyration (Rg) analysis of 100 ns trajectory. A Rg of all protein atoms. A1 Rg density plot of unbound 3MST. A2 Rg density plot of 3MST–3MP complex. A3 Rg density plot of 3MST–Rutin complex

Full size image
Fig. 6
figure 6

Solvent accessible surface area (SASA) analysis of 100 ns trajectory. A SASA of all protein atoms; A1 SASA density plot of unbound 3MST; A2 SASA density plot of 3MST–3MP complex. A3 SASA density plot of 3MST–Rutin complex. B SASA of active site. B1 SASA density plot of active site of unbound 3MST. B2 SASA density plot of active site of 3MST–3MP complex. B3 SASA density plot of active site of 3MST–Rutin complex

Full size image

Hydrogen bond and binding energy analysis

Using the simulated trajectory, hydrogen bond analysis and binding energy analysis were performed to observe the interaction between Rutin and Ld3MST. A control for the analysis was the interaction of the Ld3MST with its natural substrate. Before studying the H-bond from the trajectory, we want to examine how far the ligand is from the active site throughout the entire simulation time period. For this purpose, we computed and plotted the distance between the center of the active site cavity of the Ld3MST and the center of the corresponding ligands (Fig. 7A). It is observed that the substrate and Rutin maintained a distance less than 2 Å throughout the simulation. By comparison, we found that substrate and Rutin are not far apart in terms of distance to the active site and the difference is observed to be meager. Possibly, the mild difference observed is due to the difference in sizes of the substrate and Rutin, making their centers different. However, the results indicate that Rutin was in close proximity to the active site cavity to establish stable interactions, suggesting that Rutin has the same or better capacity to interact with Ld3MST as substrate. This was also supported by results of H-bond analysis (Fig. 7B,C), where the average number of H-bonds formed between the Ld3MST and Rutin (Fig. 7C) was higher than the average number formed between the Ld3MST and substrate. In addition, the H-bond occupancy was calculated for both complexes from the entire simulation trajectory (Table 2). As can be seen from the Table 2, Rutin has the ability to form H-bonds consistently with the Ld3MST active site residues on par with that of the substrate. Moreover, Rutin is able to interact with different segments of the active sites, as opposed to the substrate, which primarily interact with the center segment of the active site. However, substrate did not achieve H-bond occupancy of more than 10% with the central active site residue CYS253, even though it formed H-bond interactions with it in the dock pose (Fig. 3A). In contrast, Rutin did not display H-bond interaction with the CYS253 in dock pose (Fig. 3B), but had 12% H-bond occupancy in the simulation trajectory. In light of the H-bond analysis, it is evident that Rutin can interact stably in all segments of the active site as an inhibitor, blocking the active site cavity from entering the substrate binding site. As the H-bond is not the only interaction a ligand–protein complex can have, we computed the binding energy of the complexes extracted from MD trajectories every 10 ns (Table 3). The results show that substrate had better binding energy (− 6.6 kcal / mol) during the 0 and 90 ns frames and Rutin had better binding energy (− 10.6 kcal / mol) during the 20 and 40 ns frames. Upon comparing the binding energies of Rutin and substrate, the former has a better average binding energy (− 9.9 ± 0.5 kcal/mol) than the latter (− 6.3 ± 0.17 kcal/mol). Overall, based on the H-bond and binding energy analysis discussed above, it is evident that Rutin has the ability to interact with the active site and could act as a potential inhibitor of Ld3MST.

Fig. 7
figure 7

A Distance of active site of Ld3MST with 3MP and Rutin. Number of hydrogen bonds of Ld3MST with B 3MP and C Rutin

Full size image
Table 2 Hydrogen bond occupancy of 3-MST-Rutin and 3-MST-3MP complex computed from 100 ns MD trajectory
Full size table
Table 3 Binding energies of L. donovani 3MST with 3MP and Rutin from MD trajectory
Full size table

Leishmaniacidal effect and cytotoxicity study of Rutin

In addition to examining the interaction of Rutin with Ld3MST computationally, this study also assessed Rutin’s leishmaniacidal effect. The effects of different concentrations of Rutin on survival of L. donovani promastigotes and amastigotes were studied, in comparison with a negative control group without Rutin treatment and a positive control group with AmB treatment. After 48 h of treatment, effect based dose finding analysis revealed a concentration of 40.95 μM of Rutin is necessary to kill 50% of L. donovani promastigotes (Fig. 8A). A high concentration of 327.6 μM of Rutin resulted in less than 20% survival of promastigotes. A highly statistically significant difference (p < 0.001) was observed between Rutin’s leishmaniacidal effect and the untreated. AmB, however, exhibited an IC50 value of 30 nM. Similarly, after 48 h of treatment, a concentration of ~ 90.09 µM of Rutin (p < 0.001) was found necessary to kill 50% of L. donovani amastigotes (Fig. 8C). An appreciable increase in the IC50 value of Rutin for amastigotes when compared to promastigotes may be the result of reduced uptake of Rutin through macrophage membranes. However, the cytotoxicity study (Fig. 8B) revealed that the cytotoxicity of Rutin (tested at different concentrations from lower to higher) was found to have negligible toxicity to THP-1 than AmB at all concentrations. The selectivity index (ratio of cytotoxicity to the IC50) was computed and compared between AmB and Rutin. The maximum concentration tested for cytotoxicity of Rutin in this study is 500 μg/ml (819 μM), which resulted in only ~ 10% cytotoxicity. Since 50% cytotoxic concentration is required to calculate the selectivity index, we considered ~ 900 μM concentration of Rutin to be 50% cytotoxic. In addition, a single concentration of AmB (0.5 μM) was used as positive control, which resulted in ~ 80% cell death (Fig. 8C), and therefore cannot be regarded as the IC50 for AmB. Hence, in accordance with the previously reported study (Vermeersch et al. 2009), AmB’s IC50 was considered as 0.1 μM. The selectivity index of AmB on amastigotes (1 μM / 0.1 μM = 10) was observed not to be more different to that of Rutin on amastigotes (> 900 μM / 90 μM ≥ 10). In addition, Rutin concentrations as high as ~ 819 μM were not significantly different (p > 0.05) from the control, indicating that it does not affect the normal cells, even at several times of its IC50. On the other hand, AmB was toxic to normal cells at concentrations as low as 1 μM. The morphological differences between the control cells and those treated with Rutin were measured under the microscope (Fig. 8E). While the control cells are normal (oblong and flagellated) (Ambit et al. 2011) (Fig. 8E1), Rutin-treated cells are malformed (round and small) (Fig. 8E2 and E3). The change in cells might be caused by a disruption of membrane fluidity (Fernandes et al. 2017) and apoptosis phenomena (Verma et al. 2007), both of which could be triggered by Rutin. Nevertheless, the leishmaniacidal and cytotoxic effects of Rutin demonstrated in this study are more favorable and promising for its potential use as a lead against leishmaniasis.

Fig. 8
figure 8

Antileishmanial activity and cytotoxicity of Rutin. A Cell viability MTT assay on promastigotes. B Cytotoxicity assay on THP-1. C Cell viability MTT assay on amastigotes. D Activity assay of recombinant Ld3MST on substrate, sodium thiosulfate. * p < 0.05; ** p < 0.01; *** p < 0.001; NS—non-significant; AmB—AmBisome as positive control; and NC—phosphate buffer solution as negative control; C—observation of cell morphology under microscope (× 40); C1—control untreated promastigote cells with PBS; C2—control untreated THP-1 cells; C3—control untreated amastigote cells with PBS; C4—enzyme activity of Ld3MST in the absence of Rutin. E1—control untreated; E2—treated with 40.95 μM Rutin; E3—treated with 327.6 μM Rutin

Full size image

Inhibition of 3-MST activity by Rutin

Initially the Ld3MST was cloned and purified (Supplementary Fig. 6). PCR was used to amplify the Ld3MST gene using L. donovani genomic DNA as a template. In a 1% agarose gel, a distinct most intense band of about 1.1 kb was seen at 64 °C (Supplementary Fig. 6A). Purified PCR product was double digested with BamH1 and HindIII restriction enzyme. Further, digested product ran on 1% agarose gel and specific band was purified by gel extraction kit. Correspondingly, pET28a plasmid was digested by same restriction enzymes and product (~ 5.3 kb) was purified in a similar fashion. Selected colonies from transformed colony plate were subjected to colony PCR using forward and reverse primer. Further, it was confirmed by double digestion of cloned plasmid by BamH1 and HindIII and ran a gel. ~ 5.3 kb and ~ 1.1 kb band was confirmed (Supplementary Fig. 6B). The Ld3MST gene was inserted in pET28a. Protein expression was performed in BL-21 cells by IPTG induction. 0.5 mM IPTG-induced cells showed the highly expressed protein band at ~ 41 kDa as expected in SDS-page (Supplementary Fig. 6C). Purification was accomplished utilizing Ni–NTA agarose beads, protein-tagged histidine, and imidazole-based chemistry to get a particular enriched pure recombinant 3MST product. The presence of a ~ 41-kDa band on the elution indicated that the required protein had been purified specifically. It was compared to samples of uninduced, induced, and induced supernatant (Supplementary Fig. 6D).

In spite of Rutin’s leishmaniacidal activity (promastigotes and amastigotes) at various concentrations, experiments were done to see if it has the ability to inhibit the 3-MST. This was accomplished by measuring the enzyme activity (in terms of absorbance, as described in the method section) in the presence of 3-MST and its substrate (sodium thiosulfate). This was regarded as 100% enzyme activity. Furthermore, the enzyme activity was measured using different concentrations of Rutin added to the reactions. Reduction in the enzyme activity was observed to be ~ 20, ~ 35, ~ 50, and ~ 87% for Rutin concentration 1.64, 8.19, 16.38, and 40.95 μM, respectively (Fig. 8D). Except at concentration 1.64 μM, Rutin inhibited the enzyme activity significantly (p < 0.01) at other concentrations. This result suggests that Rutin holds the potency to interact with Ld3MST and interfere in the enzyme activity. In summary, both in silico results and enzyme activity results of this study demonstrate that Ld3MST is a potent target, as well as Rutin an efficient inhibitor for Ld3MST.

Conclusion

Three elements that make Ld3MST an appealing therapeutic target include its expression in all stages of Leishmania life cycle, upregulation during peroxide or sulfur stress, and structural differentiation from its human equivalent. In light of previous studies showing 3MST is important for parasite survival, especially during sulfide stress, this study was conducted to explore Ld3MST as a drug target. To the best of our knowledge, this was the first attempt to screen a small molecule against the Ld3MST. Accurate modeling of Ld3MST, screening of potent natural compounds, and in vitro testing of potential compounds were undertaken. Long-time MD simulations were performed to verify the modeling and docking results. Based on the MD simulation, the modeled protein appears to be structurally stable, and the screened compound (Rutin) is able to bind more effectively to the Ld3MST than its substrate (3MP). Specifically, the H-bond and occupancy analysis of the MD trajectories revealed that Rutin interacted consistently with various active site segments across the simulation period, suggesting that Rutin could actively block the active site cavity of Ld3MST. The binding energy results of Rutin also indicated the presence of other strong non-H-bond (in addition to H-bond) interactions with Ld3MST, in comparison to the substrate. Moreover, an evaluation of Rutin’s leishmaniacidal effect and cytotoxicity on cell lines revealed that Rutin had an IC50 of 40.95 μM (promastigote), 90.09 μM (amastigote), and exhibited no cytotoxicity even at high concentration 819 μM. Further, the study revealed Rutin’s ability to interfere with Ld3MST and reduce its enzyme activity. However, these results do not mean that Rutin could be an alternative to AmB. Rather, Rutin may be a useful compound to be used in synergy with other antileishmanial drugs or it can be redesigned to have better activity. The fact that the screened compound is non-toxic at high concentrations and leishmaniacidal at low concentrations makes it a potent small molecule of choice for further investigations.