Introduction
The phylum Apicomplexa comprises a large group of protozoan parasites. Many of them, such as Toxoplasma gondii, Plasmodium sp. and Cryptosporidium sp., cause diseases in humans, while others ail domestic animals, as is the case with Neospora caninum (Morrison, Reference Morrison2009). These organisms are characterized by the presence of the apical complex, which is a set of structures formed by rhoptries and micronemes responsible for the parasite's active invasion of the host cell (Frénal et al., Reference Frénal, Dubremetz, Lebrun and Soldati-Favre2017). Some members of this phylum have another unique organelle, the apicoplast, which has a secondary endosymbiotic origin and is involved in several cellular processes. This organelle is essential to the survival of these parasites and is an excellent target for drug development. Another peculiarity of the phylum is the highly complex life cycles of its representatives, which alternate between different hosts. This life cycle is robustly regulated by gene expression that controls the adaptation and alteration among the parasites' life forms (Frénal et al., Reference Frénal, Dubremetz, Lebrun and Soldati-Favre2017).
During this transition, the parasites undergo morphological, biochemical and molecular changes. It is currently known that about 18% of genes are regulated in a stage-specific manner by unknown transcriptional regulators in Toxoplasma sp. (Radke et al., Reference Radke, Behnke, Mackey, Radke, Roos and White2005). Many apicomplexans lack canonical transcription factors. However, a family known as AP2-like has been linked to the Apetala-2 transcription factor in plants (Balaji, Reference Balaji2005). These factors have been studied extensively in order to investigate whether they play a role in transcriptional regulation (Painter et al., Reference Painter, Campbell and Llinás2011; Iwanaga et al., Reference Iwanaga, Kaneko, Kato and Yuda2012; Oberstaller et al., Reference Oberstaller, Pumpalova, Schieler, Llinás and Kissinger2014; Jeninga et al., Reference Jeninga, Quinn and Petter2019; Srivastava et al., Reference Srivastava, White and Sullivan2020). That said, epigenetic regulation has already proven to be efficient in transcriptional regulation. Among all histone modifications, acetylation seems to be the most abundant. In this context, both canonical and histone variants are hyperacetylated, especially at the N-terminal tails. Also, these histones suffer several other modifications (Nardelli et al., Reference Nardelli, Che, de Monerri, Xiao, Nieves, Madrid-Aliste, Angel, Sullivan, Angeletti, Kim, and Weiss and M2013). Several modifiers were identified besides deacetylase and acetyltransferases (ToxoDB.org), and acetylation is one of the most abundant.
Histone deacetylases (HDACs) and histone acetyltransferases (HATs) are mainly responsible for histone acetylation/deacetylation complexes (Inoue and Fujimoto, Reference Inoue and Fujimoto1969; Seto and Yoshida, Reference Seto and Yoshida2014; Sepehri et al., Reference Sepehri, Beacon, Osman, Jahan and Davie2019). However, they can have other non-histone targets (Seto and Yoshida, Reference Seto and Yoshida2014; Verdin and Ott, Reference Verdin and Ott2015). One of the most well-known functions of these enzymes is their role in transcription. Histone acetylation, promoted by HATs, makes chromatin more open and accessible to transcription machinery, favouring gene transcription. On the other hand, HDACs remove acetylations and render chromatin inaccessible, consequently preventing transcription and silencing genes. Although this is the most well-known function of HDACs, they also play a role in other processes such as DNA repair and replication (Seto and Yoshida, Reference Seto and Yoshida2014).
HDACs are traditionally divided into four classes according to their domain, their similarity to yeast sequences, cofactors and inhibition of trichostatin A. Classical zinc-dependent HDACs belong to classes I, II and IV. Contrastingly, non-classical domains are categorized into class III; they are known as sirtuins and use nicotinamide adenine dinucleotide (NAD)+ as a cofactor.
HDACs have been studied as a target for drugs supporting the treatment of several illnesses, including cancer and degenerative diseases such as Alzheimer's disease; many HDAC inhibitors have been approved by the U.S. Food and Drug Administration thus far (Benedetti et al., Reference Benedetti, Conte and Altucci2015). In 1996, Darkin-Rattray and colleagues reported the effectiveness of HDAC inhibitors against Plasmodium sp. Since their publication, other HDAC inhibitors have been tested (Bougdour et al., Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009), and interesting effects have been observed, not only with regards to apicomplexans but for other parasites as well, like the Euglenozoa Trypanosoma sp. and Leishmania sp., cestodes and the trematode Schistosoma sp. (Andrews et al., Reference Andrews, Gupta, Tran, Fairlie, Gobert and Bozdech2012; Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017; Guidi et al., Reference Guidi, Saccoccia, Gennari, Gimmelli, Nizi, Lalli, Paonessa, Papoff, Bresciani and Ruberti2018; Vaca et al., Reference Vaca, Celentano, Toscanini, Heimburg, Ghazy, Zeyen, Hauser, Oliveira, Elissondo, Jung, Sippl, Camicia and Rosenzvit2021). Even so, little is known about the HDACs of apicomplexan species. Thus far, there are no studies focused on thoroughly addressing their structure by employing X-ray crystallography. Moreover, the binding mechanism of these inhibitors in the HDACs' structure is unknown.
Homology modelling is a powerful tool for understanding protein structures that have no solved structure. This approach is based on the principle that proteins with similar amino acid sequences have similar tertiary structures (Chahal et al., Reference Chahal, Nirwan, Kakkar and Ramasami2020). Based on the three-dimensional structure model, it is possible to perform several analyses such as docking, screening and drug discovery; it is a useful technique since many ligands can be used in this field of study.
Therefore, we modelled the HDACs from T. gondii, Plasmodium falciparum, Cryptosporidium parvum, Cyclospora cayetanensis, N. caninum and Hammondia hammondi to explore the specificity of inhibitors that are analogous to FR235222 and apicidin on those HDACs. We chose FR235222 and apicidin compounds because they have already been tested on parasites (Darkin-Rattray et al., Reference Darkin-Rattray, Gurnett, Myers, Dulski, Crumley, Allocco, Cannova, Meinke, Colletti, Bednarek, Singh, Goetz, Dombrowski, Polishook and Schmatz1996; Mori et al., Reference Mori, Urano, Abe, Furukawa, Furukawa, Tsurumi, Sakamoto, Hashimoto, Takase, Hino and Fujii2003; Bougdour et al., Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009) and their analogues were also tested for HDAC inhibition (Colletti et al., Reference Colletti, Myers, Darkin-Rattray, Gurnett, Dulski, Galuska, Allocco, Ayer, Li, Lim, Crumley, Cannova, Schmatz, Wyvratt, Fisher and Meinke2001; Miller et al., Reference Miller, Witter and Belvedere2003; Di Micco et al., Reference Di Micco, Terracciano, Bruno, Rodriquez, Riccio, Taddei and Bifulco2008; Terui et al., Reference Terui, Chu, Li, Ando, Fukunaga, Aoki and Toda2008; Ballante et al., Reference Ballante, Reddy, Zhou and Marshall2017). To this end, we performed docking followed by structural analysis to give information about the probable mode of binding and to provide relevant information for drug development.
Materials and methods
In silico modelling
Table 1 shows 14 Apicomplexan HDACs whose amino acid sequences were obtained from the following databases: Toxodb.org, Plasmodb.org and Cryptodb.org [TgHDAC2 (TGME49_249620); TgHDAC3 (TGME49_227290); TgHDAC4 (TGME49_227790); TGME49_227790; CpHDAC1 (cgd6_80); CpHDAC2 (cgd6_1380); CpHDAC3 (cgd8_480); PfHDAC1 (PF3D7_0925700); CcHDAC1 (cyc_00497); CcHDAC3 (cyc_03294); HhHDAC2 (HHA_249620); HhHDAC3 (HHA_227290); HhHDAC4 (HHA_257790); NcHDAC4 (CEL67040.1)]. These protein sequences were submitted to the HHPRED server to detect homologies and predict structures (Söding et al., Reference Söding, Biegert and Lupas2005). A sequence alignment of the target and protein using atomic structure to present the highest homology was used as an input in the Modeller program (Webb and Sali, Reference Webb and Sali2016). All alignments obtained in this study are available at https://drive.google.com/file/d/135QLmSBR4gIu0cxyhQ6nSIyXBrTuMOF/view?usp=sharing. A total of 50 models was generated for each target, and the best models were selected according to their discrete optimized protein energy (DOPE) score. For molecular modelling, ions bound to template structures were removed. The model's quality was also analysed in the ProSA-web (Wiederstein and Sippl, Reference Wiederstein and Sippl2007) and Ramachandran plot (https://zlab.umassmed.edu/bu/rama/) servers.
Proteins selected for analysis
Proteins used in this study, their gene ID and the corresponding templates used to generate homology models.
It is also worth mentioning that TgHDAC4 was misannotated in the database, and our group noticed that the TgHDAC4 mRNA was longer than expected. Thereupon, Dr Arnab Pain's group (King Abdullah University of Science and Technology, Saudi Arabia) reanalysed the RNA-seq data and helped us to properly reassemble the gene by including two more exons that had been missing in the first annotation (Ramaprasad et al., Reference Ramaprasad, Mourier, Naeem, Malas, Moussa, Panigrahi, Vermont, Otto, Wastling and Pain2015).
Molecular docking
Molecular docking was performed using Chimera (Pettersen et al., Reference Pettersen, Goddard, Huang, Couch, Greenblatt, Meng and Ferrin2004) and AutoDocK Vina (Trott and Olson, Reference Trott and Olson2010) programs. The structure of human HDAC2 linked with the HDAC inhibitor suberoylanilide hydroxamic acid (SAHA) was used as the reference (file code: 4LXZ) to set-up the grid box (Supplementary Fig. S1).
Docking simulations were performed with 12 compounds (1: FR235222; 2: CHEMBL 4099599; 3: CHEMBL 2370998; 4: CHEMBL 358276; 5: CHEMBL 472574; 6: CHEMBL 470325; 7: CHEMBL 492600; 8: CHEMBL 488096; 9: CHEMBL451169; 10: CHEMBL 1793810; 11: CHEMBL 1793988; 12: apicidin, 13: panobinostat, 14: chloroquine) (Fig. 1). These compounds were named ligands 1 to 14. The PDB structures of these compounds were downloaded from the ChEMBL database. They were chosen for being analogues of FR235222 (compounds 4, 5, 6, 7 and 8) and apicidin (compounds 2, 3, 9, 10 and 11). In addition, compounds 13 and 14 are hydroxamate-based inhibitors and were chosen based on previous in vitro assays indicating their ability to inhibit Plasmodium HDACs (Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017).
Ligand structures used in this study. Ligand 1 (FR235222), ligand 2 (CHEMBL4099599), ligand 3 (CHEMBL2370998), ligand 4 (CHEMBL358276), ligand 5 (CHEMBL472574), ligand 6 (CHEMBL470325), ligand 7 (CHEMBL492600), ligand 8 (CHEMBL488096), ligand 9 (CHEMBL451169), ligand 10 (CHEMBL179381), ligand 11 (CHEMBL1793988) and ligand 12 (apicidin); ligand 13 (panobinostat) and ligand 14 (chloroquine). Ligand 7 is the trans-form of ligand 5, and ligand 6 is the cis-form of ligand 8.
The docking simulations were performed as follows: (1) each HDAC model was checked for missing atoms, bonds and contacts; subsequently, energy minimization was assessed using the molecular modelling toolkit (MMTK) package on Chimera version 1.14 employing the following parameters: Amber ff14SB force field; 100 steepest descent steps, with a steepest descend step size of 0.02 Å; 10 conjugate gradient steps, with a conjugate gradient step size of 0.02 Å. The minimized structure was used as the receptor for docking analysis. (2) AutoDock Tools module was used to generate pdbqt input files. (3) AutoDock Vina algorithm was used to perform a docking with the selected ligands. The docking simulation was then run at an exhaustiveness = 8. The docking results were evaluated using the lowest binding affinity score (kcal mol−1) predicted by the built-in scoring function of the AutoDock Vina module.
The feasibility and robustness of our interactions are strengthened by comparing the docking results of human HDAC2 interaction with SAHA and the corresponding crystallographic structure was deposited under the code 4LXZ. Both ligand structures presented a root-mean-square deviation (RMSD) of 1.72 Å, indicating both docking and experimental structures are very similar.
Structural analysis
All structural analyses and creation of images were carried out with PyMOL (the PyMOL Molecular Graphics System, version 2.4 Schrödinger, LLC) and Chimera version 1.14 (Pettersen et al., Reference Pettersen, Goddard, Huang, Couch, Greenblatt, Meng and Ferrin2004). Interactions between protein and ligands were analysed using PDBePISA (proteins, interfaces, structures and assemblies PISA service) at the European Bioinformatics Institute (http://www.ebi.ac.uk/pdbe/prot_int/pistart.html) and/or using PDBsum (Laskowski et al., Reference Laskowski, Jabłońska, Pravda, Vařeková and Thornton2018).
Molecular dynamics
Molecular dynamics simulations were performed for a period of 30 ns for CpHDAC1, CpHDAC2, PfHDAC1 and TgHDAC3 models. GROMACS 2020.4 (Abraham et al., Reference Abraham, Murtola, Schulz, Páll, Smith, Hess and Lindahl2015) and CHARMM36m (version jul-2021) (Huang et al., Reference Huang, Rauscher, Nawrocki, Ran, Feig, de Groot, Grubmüller and Grubmüller2017) force field were used to prepare topologies and coordination files. Protonation states were assigned using H++ server version 4.0 (Anandakrishnan et al., Reference Anandakrishnan, Aguilar and Onufriev2012) considering the pH 7.4. The models were immersed in a dodecahedron simulative box, 15 Å away from the protein surface, filled with transferable intermolecular potential with 3 points (TIP3P) (Jorgensen et al., Reference Jorgensen, Chandrasekhar, Madura, Impey and Klein1983) explicit solvent molecules and rendered electroneutral by the introduction of sodium counterions.
A first round of minimization using the steepest descent algorithm (Haug et al., Reference Haug, Arora and Matsui1976) was performed using several steps where maximum force <1000.0 kJ mol−1 nm−1 in order to regularize the structures. Optimization and relaxation of solvent and ions were performed in two steps (isothermal–isochoric ensemble and isothermal–isobaric ensemble) by simulating the system at a temperature of 310 K and pressure of 1 bar. The simulation was carried out under periodic boundary conditions with a time step of 2.0 fs at a constant temperature of 310 K using the modified Berendsen thermostat (V-rescale algorithm) (Bussi et al., Reference Bussi, Donadio and Parrinello2007). Electrostatic interactions were treated using the particle-mesh Ewald algorithm (Essmann et al., Reference Essmann, Perera, Berkowitz, Darden, Lee and Pedersen1995). The resultant trajectory, RMSDs and root-mean-square fluctuations were processed using the GROMACS 2020.4 package.
Results
Homology models
We generated structural models of 14 apicomplexan HDACs by homology modelling (Table 1, Fig. 2) since all proteins showed identity with their respective templates above 30% (Sánchez and Sali, Reference Sánchez and Sali1997). The low sequence homology impaired the modelling of other HDACs such as TgHDAC1, PfHDAC2 and PfHDAC3. These proteins have a long unique sequence of amino acids which could generate unreliable structural data, so we decided to remove them from the analysis. Initially, we observed that almost all the proteins had insertions (mainly C-terminal) with no fold similarity to other known structures. Therefore, the function of those insertions is unknown and should be further investigated through experimental structural approaches such as X-ray crystallography and nuclear magnetic resonance (NMR) for more accuracy in structure determination. In this study, we focused on HDAC domains showing structural homology. Table 2 displays the non-modelled C-terminal regions of all of the proteins (insertions with no previous structure). The quality assessment of the best models (selected based on DOPE score) analysed by ProSA-web and Ramachandran plot servers showed that at least 91% of the residues fell within the most favourable regions (Table 3) (Anderson et al., Reference Anderson, Weng, Campbell and Jiang2005). The CpHDAC2 and TgHDAC3 were the best models, with more than 99% of the residues in favourable regions of the Ramachandran plot. CpHDAC1 has shown the best Z score value (−9.31) (Wiederstein and Sippl, Reference Wiederstein and Sippl2007). The Z score indicates a more favourable model according to X-ray crystallography and NMR experimental data present in the PDB (Wiederstein and Sippl, Reference Wiederstein and Sippl2007). This approach is used to validate a model and identify an erroneous structure. Protein misfolding models have a high energy and Z score; the lower the Z score, the higher the probability the model is correct (Sippl, Reference Sippl1993). In summary, these analyses show that all the models in this study are of good quality.
Tridimensional models of apicomplexan HDACs.
Deleted regions
Many apicomplexan HDACs have insertions that cannot be modelled because they are not homologous to the template protein; therefore, they were not considered in our study.
Analysis of the quality of the models
Protein models were selected based on the best DOPE scores and were analysed by Ramachandran plot and Z score to validate the model.
Docking
We analysed the binding potential of 12 compounds against 14 proteins. The molecules used for molecular docking were chosen due to their known activity against parasites such as apicidin (Darkin-Rattray et al., Reference Darkin-Rattray, Gurnett, Myers, Dulski, Crumley, Allocco, Cannova, Meinke, Colletti, Bednarek, Singh, Goetz, Dombrowski, Polishook and Schmatz1996) and FR235222 (Mori et al., Reference Mori, Urano, Abe, Furukawa, Furukawa, Tsurumi, Sakamoto, Hashimoto, Takase, Hino and Fujii2003; Bougdour et al., Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009). Other compounds analogue to these were also selected to evaluate the affinity and chemical properties that can affect affinity
Table 4 describes the binding energy (in kcal mol−1) regarding the best placement of 12 compounds against each protein. TgHDAC5 and HhHDAC2 were the proteins containing more ligands with potential binding affinity with values above −7.5 kcal mol−1. Of all of the HDACs, TgHDAC5 and HhHDAC2 were the proteins that had more interactions below −7.5. Our analysis demonstrated that ligands 4, 5 and 9 have the potential to inhibit most apicomplexan HDACs.
Ligand binding energy against apicomplexan HDACs
This table shows the scores for molecular docking binding energy. The molecular docking was performed with 14 ligands against 14 apicomplexan HDACs.
Among the three compounds, ligand 9 was the one that interacted with the largest number of species. Also, ligand 9 was unique to 4 of the 6 apicomplexans analysed. Moreover, no studied compounds displayed a strong affinity to C. cayetanensis and N. caninum HDACs as observed through their binding affinity values above −7.0 kcal mol−1.
Ligand 4 is analogous to FR235222 and, according to the ChEMBL database, it exhibits inhibitory activity against HT-29 and Jurkat human tumour cells at concentrations of 0.025 and 0.019 nm, respectively (Miller et al., Reference Miller, Witter and Belvedere2003). We emphasize the strongest interaction of this ligand with TgHDAC3. Although this protein is known to be inhibited by the FR235222 compound (Bougdour et al., Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009), the interaction with ligand 4 appears to be stronger than what occurs with the original FR235222 compound.
Ligand 5 is also an analogue of FR235222, synthesized to be a simplification of the original compound. Di Micco et al. (Reference Di Micco, Terracciano, Bruno, Rodriquez, Riccio, Taddei and Bifulco2008) observed that this compound presented a better pharmacological profile, better results in inhibition in vitro tests, and therefore a better specificity compared to FR235222, which is in line with our findings. Furthermore, this compound showed strong interactions with several HDACs addressed in our study. The strongest interaction was with TgHDAC5, indicating that ligands 4 and 5 have the broadest inhibiting capacity among the apicomplexan HDACs.
Compound 9 is analogous to apicidin, with about 87% similarity. This compound was described in 2008 by Terui and collaborators through the fermentation of Nonomuraea sp. TA-0426 (Terui et al., Reference Terui, Chu, Li, Ando, Fukunaga, Aoki and Toda2008). Besides an HDAC inhibitor, compound 9 is also a glycine transporter type I (GlyT1) inhibitor, a promising target for an antipsychotic treatment (Terui et al., Reference Terui, Chu, Li, Ando, Fukunaga, Aoki and Toda2008).
We also wish to highlight the interaction of ligand 3 with CcHDAC1, which appears to be specific. The half-maximal inhibitory concentration (IC50) of this compound has already been determined for P. falciparum (0.034 nm) and C. parvum (4 nm) (Miller et al., Reference Miller, Witter and Belvedere2003). This compound was also tested against HeLa cells and effectively inhibited HDACs in this tumour lineage (Miller et al., Reference Miller, Witter and Belvedere2003). Ligand 6 was synthesized as a variation of FR235222, and we observed that for CcHDAC1 and CpHDAC3, ligand 6 seems to be stronger and poses a more specific interaction than the original compound. Interestingly, ligand 8, which is the trans-form of ligand 6, was more specific in the interaction with HhHDAC4, indicating an important difference between the conformation of the ligand and the strength of interaction with the protein. In addition to HhHDAC4, TgHDAC4 and CcHDAC3 are also potential targets of compound 8.
Compound 11 is among those that least interacted with apicomplexan HDACs. This apicidin analogue had already been tested with P. falciparum and Eimeria tenella, presenting negligible inhibitory potential. The IC50 values were 782 718.90 and 711 562.64 nm, respectively, much higher than what was the case for the original compound (Colletti et al., Reference Colletti, Myers, Darkin-Rattray, Gurnett, Dulski, Galuska, Allocco, Ayer, Li, Lim, Crumley, Cannova, Schmatz, Wyvratt, Fisher and Meinke2001). The authors indicated that indole modification of the original compound by 5-bromoindole reduces the biological activity and specificity of the compound in relation to the original, which is supported by our study, where apicidin also showed better results.
Regarding hydroxamate-based inhibitors, panobinostat seemed to have more affinity to tested HDACs with greatest scores for TgHDAC2, TgHDAC4 CpHDAC2, NcHDAC4 and CcHDAC1. Chloroquine theoretical affinity was observed mainly for NcHDAC4.
Molecular dynamics analysis
After 30 ns of simulation, the low RMSD variation indicates the convergence of simulation (Supplementary Fig. S2A). The ligand binding sites (Supplementary Fig. 2B–G) do not present high flexibility (RMSD values of about 0.3 nm). Most flexibility (RMSD values of 0.4 nm) might reside upon loops comprised of residues 16–27 (of CPHDAC1) and 198–205 (of PfHDAC1). The low values of RMSD indicate that conformational flexibility of the binding sites should not substantially alter the results of molecular docking.
Structural analysis of interactions
HhHDAC2 is a class I HDAC that interacts with the FR235222 compound through van der Waals interactions (Fig. 3A). The interaction of HhHDAC2 with ligands 4 and 5 is greater due to the hydrogen bonds (Fig. 3B and C, Table 4). Both compounds interact with D363 and H327 (HhHDAC2 numbering). The difference between the interaction of compounds with this HDAC relies on ligand 4 interactions with H365 and ligand 5 interactions with G488 (Fig. 3B and C). The contribution of interactions between ligand 4 and H327 results in more energy binding compared with ligand 5 (Fig. 3B and C, Table 4).
Interaction map with HhHDAC2. Bidimensional map of the interaction between HhHDAC2 and ligands (A) FR235222, (B) 4 and (C) 5, showing van der Walls interactions (
) and hydrogen bonds (
).
Regarding ThHDAC5, its interaction with FR235222 is potentially more robust than with HhHDAC2. While HhHDAC2 only formed van der Waals interactions, TgHDAC5 formed hydrogen bond interactions with H638, D363, D898, R409 and G943 (TgHDAC5 numbering). An exchange in the interaction of G943 by R574 was observed with ligand 4. The greater number of van der Walls interactions with ligand 5 might explain its slightly higher potential binding affinity for TgHDAC5.
The ThHDAC5–apicidin interaction includes hydrogen bonds with R574 and G943 (ThHDAC5 numbering). We observed that the substitution of R574 by R409 and the addition of G943 to the ligand–ThHDAC5 interface contribute to a greater potential affinity with ligand 9 than with apicidin (Fig. 4). Therefore, we conclude that the potential residues for inhibition are: D363 (HhHDAC2) and R409, D898, H638 (TgHDAC5). R409 seems to be a critical residue in these potential interactions.
Interaction map with TgHDAC5. Bidimensional map of interactions between TgHDAC5 and ligands (A) FR235222, (B) 4, (C) 5, (D) apicidin and (E) 9 showing van der Walls interactions (
) and hydrogen bonds (
).
In CcHDAC1 interaction with FR235222 and its analogues, the binding pattern seemed completely different. FR235222 interacts through hydrogen bonds with F268, H266 and G265 (Fig. 5A, CcHDAC1 numbering). On the other hand, ligand 6 interacts with E248, K207 and N276, which appears to result in a stronger interaction (Fig. 5B, Table 4). Notably, the improvement in the binding of CcHDAC1 to ligand 3 arises from the hydrogen bonds with H241 and G14 when compared to the original apicidin compound, which forms only van der Waals interactions (Fig. 5C and D).
Interaction map with CcHDAC1. Bidimensional map of interactions between CcHDAC1 and ligands (A) FR235222, (B) 6, (C) apicidin and (D) 3 showing van der Walls interactions (
) and hydrogen bonds (
).
An interaction that caught our attention was the one between HhHDAC4 and ligand 8, as this is the trans-form of ligand 6 and proved to be much more relevant for inhibiting this protein as well as for other class IV HDACs. In this protein, ligand 6 interacts with residues Y13, G11 and K119 (Fig. 6A). Of these residues, only K119 remained in the interaction map with ligand 8, which also forms three hydrogen bonds with this residue. The interactions with I70, K69 and L72 also appear to be important (Fig. 6B).
Interaction map with HhHDAC4. Bidimensional map of the interactions between HhHDAC4 and ligands (A) 6 and (B) 8 showing van der Walls interactions (
) and hydrogen bonds (
).
Structural comparison among human and apicomplexan HDACs
Superposing of 14 modelled apicomplexan HDACs and human HDAC crystallographic structures (HDAC 2 PDB 4LXZ, HDAC 1 PDB 4BKX, HDAC 3 PDB 4A69, HDAC 8 PDB 5FCW) showed that ligand binding site is very conserved among HDACs. Human HDAC regions 89–97 (HDAC2 numbering) are not conserved in apicomplexan HDACs (Fig. 7A). Among apicomplexans, the ligand binding site of NcHDAC4 diverge more from human HDACs (Fig. 7B).
Similarity among apicomplexan and human HDACs. (A) Superposing of human HDAC2 (PDB 4LX2, blue) and HhHDAC4 (green) representing the conserved fold among HDACs and highlighting the divergent helix in human HDACs (red). (B) Superposing of human HDAC2 (PDB 4LX2, blue) and NcHDAC4 (magenta). Ligand binding site is shown as blue dots.
Discussion
Apicomplexa is a large phylum comprising more than 5000 species, with many non-pathogenic organisms. However, several species pose a risk to public health, such as T. gondii, Plasmodium sp. and Cryptosporidium sp., which cause toxoplasmosis, malaria and cryptosporidiosis, respectively. Although well known, these diseases are still neglected due to a lack of effective treatment. In the case of toxoplasmosis, there has been no development of new drugs for over 60 years. Moreover, the commonly used drugs are ineffective in the chronic phase of the disease and cause many side-effects, leading many patients to abandon their treatment (Dunay et al., Reference Dunay, Gajurel, Dhakal, Liesenfeld and Montoya2018). On the other hand, malaria presents a high rate of drug resistance, which should prompt more efforts to develop a more effective treatment (WHO Guidelines Approved by the Guidelines Review Committee, 2015).
Inhibitors such as trichostatin A, apicidin and FR235222 have already been tested on some parasites with promising results (Darkin-Rattray et al., Reference Darkin-Rattray, Gurnett, Myers, Dulski, Crumley, Allocco, Cannova, Meinke, Colletti, Bednarek, Singh, Goetz, Dombrowski, Polishook and Schmatz1996; Bougdour et al., Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009; Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017). However, little is known regarding apicomplexan HDACs; thus far, a few articles have been published concerning their structure, and most of them focus on P. falciparum HDAC1, which we also address in this study (Melesina et al., Reference Melesina, Robaa, Pierce, Romier and Sippl2015; Kumar et al., Reference Kumar, Dhar and Subbarao2018). PfHDAC1 is expressed in the nucleus of gametocytes in mature blood stages, and its inhibition leads to a global transcriptional dysregulation. Consequently, PfHDAC1 is considered to be a good target for inhibitors (Joshi et al., Reference Joshi, Lin, Chiang, Goldman, Fujioka, Aikawa and Syin1999; Chaal et al., Reference Chaal, Gupta, Wastuwidyaningtyas, Luah and Bozdech2010).
This study brings new and important structural data regarding several HDACs. We performed homology modelling of 14 HDACs from 6 apicomplexan parasites. However, this method presents limitations since it depends on homology similarity greater than 30%. Therefore, we were unable to model some of these proteins' regions since they have many insertions with unknown function and structure. The Saccharomyces cerevisiae HDACs have a C-terminal coiled-coil domain that functions as a scaffold for the HDAC complex, an evidence that the insertion regions may play a similar role in the Apicomplexa (Lee et al., Reference Lee, Bollschweiler, Schäfer and Huber2021). We suggest more in vitro studies on these proteins should be conducted.
In addition to generating models of HDAC proteins, we also carried out molecular docking with 12 HDAC inhibitors, having FR235222 and apicidin as the positive controls along with 10 other analogous compounds obtained through the ChEMBL database. The validation and analysis of the dockings was carried out, and satisfactory models were obtained as a result.
Some inhibitors tested here have not yet been tested on parasites before, and the docking results form a basis as a guide for future studies. Bougdour et al. (Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009) demonstrated the activity of FR235222 on T. gondii HDAC3, with specific inhibition and increased expression of bradyzoite genes. The authors also detected hyperacetylation of histone H4 in promoter regions.
Ligands 4 and 5, analogous to FR235222, were synthesized to be simpler and more efficient compounds compared to the original compound (Di Micco et al., Reference Di Micco, Terracciano, Bruno, Rodriquez, Riccio, Taddei and Bifulco2008). These compounds showed good results in tumour cells such as HT-29 and Jurkat at a nanomolar scale. Based on the docking scores of the present study, we recommend that their effectiveness against parasites should be tested. This is particularly valid for T. gondii, for which Di Micco et al. (Reference Di Micco, Terracciano, Bruno, Rodriquez, Riccio, Taddei and Bifulco2008) achieved promising results in all docked HDACs, especially TgHDAC5.
Bougdour et al. (Reference Bougdour, Maubon, Baldacci, Ortet, Bastien, Bouillon, Barale, Pelloux, Ménard and Hakimi2009) demonstrated the crucial role of the conserved and Apicomplexa HDAC-specific residue TgHDAC3 T99 in the inhibitory activity of FR235222. We observed that TgHDAC3 might be able to bind not only to this catalytic site with a theoretical affinity of −7.1 kcal mol−1, but also to the region that comprised residues 64–67, 165–168 and 187–192 with similar affinity.
Our interaction analysis showed that in some apicomplexans, such as HhHDAC2, the compound FR235222 interacted only through weak bonds, such as van der Waals. Ligands 4 and 5 interact through stronger hydrogen bonds. Furthermore, we verify that targeting R574 (as observed in TgHDAC5–ligand 4) instead of G943 (as observed in HhHDAC2–ligand 4) is important for increasing affinity. Additionally, the ligand proximity between TgHDAC5 and R409 seems to be important for improving the compound affinity.
Apicidin was the first HDAC inhibitor tested on P. falciparum and has been proven to be an effective drug both in vitro and in vivo when tested on mice (Darkin-Rattray et al., Reference Darkin-Rattray, Gurnett, Myers, Dulski, Crumley, Allocco, Cannova, Meinke, Colletti, Bednarek, Singh, Goetz, Dombrowski, Polishook and Schmatz1996). Since then, many analogues have been developed. In this study, among apicidin analogues, ligand 9 showed the best results. This ligand is an inhibitor of GlyT1 and a target for an antipsychotic treatment (Terui et al., Reference Terui, Chu, Li, Ando, Fukunaga, Aoki and Toda2008). We observed that the difference between the binding of ligand 9 and apicidin resides in the replacement of the interaction with R574 by R409; in addition to this, both interactions are closer in ligand 9.
Ligand 3, an apicidin analogue, enhances potential binding by including several hydrogen bonds with CcHDAC1 that were not observed in apicidin–CcHDAC1 interaction and therefore contributed for better docking scores. Ligand 3 had already been tested in P. falciparum and C. parvum, and its IC50 was 0.022 ng mL−1 and 4 nm, respectively (Miller et al., Reference Miller, Witter and Belvedere2003). However, this compound appears to be even more specific to C. cayetanensis as well as ligand 6, which also proved to be highly specific to CcHDAC1.
The binding of hydroxamate-based inhibitors to PfHDACs was previously tested in vitro, with panobinostat and chloroquine presenting IC50 of 4 and 10 nm, respectively, the best values against the tested inhibitors (Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017). Our docking results corroborate with a stronger affinity to panobinostat than chloroquine. Also, the affinity of panobinostat to TgHDAC2, TgHDAC4 CpHDAC2, NcHDAC4 and CcHDAC1 might be stronger than that observed for Plasmodium HDAC (Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017). Also, chloroquine affinity to NcHDAC4 might be higher than that for Plasmodium HDAC (Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017).
Many studies have reported the effect of HDAC inhibitors on P. falciparum (Andrews et al., Reference Andrews, Gupta, Tran, Fairlie, Gobert and Bozdech2012; Engel et al., Reference Engel, Jones, Avery, Sumanadasa, Ng, Fairlie, Skinner-Adams and Andrews2015; Chua et al., Reference Chua, Arnold, Xu, Lancelot, Lamotte, Späth, Prina, Pierce, Fairlie, Skinner-Adams and Andrews2017). However, very little is known about other apicomplexan parasites of medical interest, such as T. gondii and C. parvum. We have hereby modelled and provided insights regarding these parasites' HDACs, discussing promising drugs for treating these diseases as well as potential residues to be targeted when considering drug development. We emphasize the importance of further studies, especially in vivo, to confirm the inhibitory potential of these compounds and promote a better understanding of the structure of these proteins.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182022000427
Author contributions
C. M. S., S. C. N. and T. A. C. B. S. conceived and designed the study. C. M. S., M. S. I. F. and V. R. S. performed homology modelling, molecular docking and structural analysis of complex with coordination of T. A. C. B. S. I. V. B. conducted and analysed molecular docking simulations. C. M. S., S. C. N., M. S. I. F., I. V. B. and T. A. C. B. S. wrote the article. All authors reviewed the submitted version of the article.
Financial support
This research received funding from FIOCRUZ, Inova-Geração de conhecimento (Fiotec VPPCB-007-FIO-18-2-92), PROEP 442339/2019-4, PRONEX (Programa de Apoio a Núcleos de Excelência: Edital: PI 01/2018 – Acordo FA/CNPqFA 41.2018) and CAPES.
Conflict of interest
The authors declare there are no conflicts of interest.
Ethical Standards
Not applicable.
