Computational insights into organophosphate interactions with oligochaeta hemoglobin: Implications for environmental safety and molecular design
Article information
Abstract
Pesticides, including organophosphates, have been reported to cause important environmental impact through effects over different taxonomic groups. Oligochaeta are often used as bioindicators, but little is known regarding molecular level pesticide interactions in this group. In our study we present a comprehensive in silico analysis of the interactions between 52 organophosphates and Lumbricus erythrocruorin hemoglobin. We performed a molecular docking analysis with GOLD software, to assess the two organophosphates most likely to interact with the studied protein, being Glufosinate ammonium and Fonofos. Next, we used Desmond software for Molecular Dynamics Simulation (MDS), to elucidate the potential mechanistic effects of these widely used pesticides. MDS of both ligands showed the potential of interacting with heme groups, and with residues important for chains interface, which may affect hemoglobin functioning. Our findings advocate for the application of computational methodologies in environmental toxicology, aiming to guide the development of agrochemicals that minimize ecological damage. It underscores the critical need for environmentally conscious chemical design and calls for further research into the subtle molecular interactions affecting non-target species within agricultural ecosystems.
Introduction
For many years, more than 500 pesticides have been extensively used worldwide (Pesticide Action Network International's 2022). Despite its potential impacts to the environmental and health, in an important food producer like Brazil, only 133 of the more than 500 pesticides were banned. Many of these molecules have a high potential to be absorbed and harm the life cycle of organisms that should not be the target [1]. Considering this fact, the indiscriminate use of pesticides can have a significant socio-environmental impact on fauna, flora, and microbiota.
To assess environmental impacts, some organisms are used as bioindicators, such as Oligochaeta (Annelida, Clitellata) [2]. Numerous studies have been conducted to investigate the impact of pesticides on the life cycle of these animals, which play a crucial role in soil ecology [3–5]. These organisms may be affected by pesticides from different ways. In the case of hemoglobin, the protein responsible for oxygen transport, some allosteric modulators, including organic phosphates, can cause changes in oxygen transport [6]. However, there is a lack of research focused on directly analyzing the influence these substances have on hemoglobin.
The hemoglobin of the species Lumbricus terrestris Linnaeus, 1758, previously known as Lumbricus erythrocruorin, is a dodecamer, meaning it has 12 subunits, which are influenced by the binding of substances to this macromolecule [7]. As an allosteric protein, it is very important to understand the influence that certain molecules can have when binding with the hemoglobin. Various studies have shown that, in addition to the heme groups, some molecules can bind to other sites of hemoglobin, affecting its oxygen transport activity [8,9].
In this sense, we investigated the interaction between 52 organophosphates and the hemoglobin of L. terrestris, employing in silico approaches such as molecular docking and molecular dynamics (MD).
Materials and Methods
Molecular docking
We retrieved the crystal structure of carboxyhemoglobin of L. terrestris from PDB (PDB ID 1X9F) at 2.60 Å resolution [7]. As the protein is a dodecamer, subunits A and D were chosen for further studies. Due to the lack of a canonical UniProt sequence for this oligomeric assembly, we adopted the residue numbering defined in the crystallographic model. This ensures consistency with both our docking and MD simulations and enables accurate reproduction of the results.
The ligand binding site was predicted using PrankWeb [10]. Fifty-two organophosphates were used in this study [11]. The ligands structures were built in silico with Discovery Studio Visualizer (v.2020.1), with the protonation states at 7.4, predicted using MarvinSketch software (“MarvinSketch, Version 12.2.27.0. ChemAxon, San Diego, CA, USA.”, [s.d.]). The energy of the structures were minimized through MOPAC PM7 semi-empirical method [12]. We performed the molecular docking with GOLD suite software (V 5.4) [13], using ChemPLP as scoring function. Protein structure was prepared by adding hydrogen atoms. The 10 Å search radius was defined as the center of PrankWeb predicted binding site, between both subunits.
Molecular Dynamics (MD)
We selected for the molecular dynamics, the two ligands with the highest scoring in molecular docking. Thus, three experiments were prepared: Apo-Hem, which consists of hemoglobin without the presence of the studied ligands; Glufo-Hem, hemoglobin in the presence of ammonium glufosinate; Fono-Hem, hemoglobin in the presence of fonofos. All steps related to MD were carried out using the Desmond Engine software [14], with the OPLS-2005 force field. The system was similarly prepared for Apo-Hem, Glufo-Hem, and Fono-Hem. Each of the models was solvated using the TIP3 solvent model, with boundary conditions in an orthorhombic box with periodic boundary conditions, maintaining a distance of 10 Å from the protein. Na+ ions were added to neutralize the system, and the system was constructed at a concentration of 0.15 M NaCl. The system was minimized in a 100 ps simulation. The MD simulation was carried out for 100 ns using NPT at a temperature of 300 K, and the system was relaxed using the default protocol of the software.
We calculated analysis parameters for the results such as Root mean square deviation (RMSD), Root mean square fluctuation (RMSF), for protein and ligands, as well as other parameters, using the Desmond Engine software. The analysis of interactions was done with the assistance of the Discovery Studio Visualizer software (v.2020.1), and figures related to interactions were made in the PyMol software (v. 2.5.2) (Schrödinger LCC, New York, NY, USA).
Results and Discussion
Our study predicted an allosteric binding site for hemoglobin, as well as the possibility of its interaction with 52 organophosphates, which are still used worldwide. MD also indicates a high potential of organophosphates interfering in oxygen transport in Oligochaeta.
The presence of allosteric binding sites in hemoglobins has been well documented [15–17], and in our study, the predicted binding cavity for Lumbricus terrestris hemoglobin, identified via PrankWeb, presented a conservation score of 2.265 (maximum 4). In fact, Anashkina et al. [18] demonstrated through molecular docking and experimental data that human hemoglobin contains at least four internal binding pockets for small organic molecules such as glutathione. These sites, located within subunit interfaces, are involved in allosteric regulation through noncovalent interactions. Although L. terrestris hemoglobin has a dodecameric arrangement, the predicted site in our model shares topological similarities with these conserved cavities. Therefore, our proposed binding site is not solely based on computational prediction, but is supported by structural and functional analogies with experimentally validated sites in other hemoglobins [18,19]. Despite the predicted binding site in L. terrestris (PDBID 1X9F) is structurally distinct from the 2,3-BPG site (PDBID 1B86), both lie within internal protein cavities and involve noncovalent interactions interfering in protein stabilization and modulation.
The study of molecular docking highlighted the potential interactions of all organophosphates, and the scoring table observed for all ligands is available as Supplementary Table 1. Table 1 lists the top 5 organophosphates with the highest docking scores. It is important to emphasize that docking scores do not represent absolute binding affinities; rather, they serve as a comparative metric to prioritize compounds within large libraries. In this context, we selected ligands for MD simulations not only based on score ranking, but also considering the number and types of predicted interactions with key hemoglobin residues, particularly those associated with heme stabilization and inter-subunit interfaces.
Best ranked docked ligands using ChemPLP scoring function.
Following the results of molecular docking, we carried out the three molecular dynamics simulations. The Apo-Hem simulation was conducted to evaluate whether ligand binding induces structural perturbations in hemoglobin. Due to the large size and complexity of the native dodecameric assembly, our MD simulations were restricted to two representative subunits (chains A and D), which is related to the predicted binding site and relevant interfacial residues. While this reduction limits the assessment of global allosteric effects and long-range cooperativity, it provides sufficient resolution to explore local interactions within the binding pocket. Future studies should consider simulations of the full dodecamer or employ coarse-grained approaches to better understand potential allosteric propagation mechanisms across the multimeric complex.
Regarding glufosinate ammonium, the molecule forms various hydrogen bonds with both chains of hemoglobin, especially with chain A. Additionally, being a molecule with 4 rotatable bonds and three charges, two negative and one positive, it demonstrates a high potential for interactions with various residues. Water molecules appear to play an essential role in establishing these hydrogen bonds, with 5 out of 11 being formed through water bridges (Figure 1).
MD predicted interactions between Lumbricus terrestris hemoglobin (PDB ID 1X9F) and glufosinate ammonium. Atom colors are shown as follows: nitrogen in blue, oxygen in red, iron in red, sulfur in yellow, phosphorus in orange and carbon in gray (amino acids) or white (ligand). Non-polar hydrogen atoms were omitted for clarity. Interactions are shown as dotted lines: hydrogen bonds in green. Chains A and D are indicated with residue identification.
The Fonofos molecule has a total of 5 rotatable bonds, but as it is not a charged molecule, it seems to have a lower potential for interaction with amino acid residues compared to Glufo-Hem. Although it makes only one interaction with amino acid residues, a π-stacking with HIS89 of chain D, fonofos forms alkyl-alkyl, π-alkyl interactions, and a hydrogen bond, through a water bridge, with the heme group of chain D, as well as an alkyl-alkyl interaction with the heme group of chain A (Figure 2). According to the experiment, the aromatic ring of fonofos plays an essential role in the positioning of the ligand in the binding site, with its N-ethyl groups being responsible for alkyl-alkyl interactions with both heme groups.
MD predicted interactions between Lumbricus terrestris hemoglobin (PDB ID 1X9F) and fonofos. Atom colors are shown as follows: nitrogen in blue, oxygen in red, iron in red, sulfur in yellow, phosphorus in orange and carbon in gray (amino acids) or white (ligand). Non-polar hydrogen atoms were omitted for clarity. Interactions are shown as dotted lines: hydrogen bonds in green, alkyl-alkyl in red, π-alkyl in black and π-π stacking in blue. CO(D) is shown as not bond to HEM(D), but it is due to the representation of PyMol, the reported distance occurred in all simulations, including Apo-Hem. Chains A and D are indicated with residue identification.
During our simulations, both the RMSD and RMSF of the global alpha carbon chains did not vary more than 3.5 Å, suggesting that the interaction with the ligands does not induce major conformational changes in the protein. Despite this, the regions with higher RMSF are those that do not form α-helices, particularly in the C-terminal region of chain A and the N-terminal of chain D, with the RMSF being slightly higher in simulations in the presence of ligands (Figure 1). The Protein Secondary Structure Elements (SSE), composed of α-helices, remained close to 60% in all simulations, being 62.11% for Apo-Hem, 62.83% for Glufo-Hem, and 59.56% for Fono-Hem.
The residues involved in interactions with the ligands, which are also related to the stabilization of the heme group in the protein subunits, underwent changes that may be significant (Figures 1 and 2). Both ligands interact with the F3 position of chain D, the H89 residue, and Glufo-Hem also interacts with F3 of chain A, H96, and the R65 residue (D), corresponding to the E10 position of the E helix. Both the F3 positions of the F helix and E10 of the E helix are associated with the stabilization of the heme group in hemoglobins [7,20–22]. Table 2 shows the distances found in the different structures, original and simulations, for the residues of the A/D interface, F3, and E10, and heme groups.
Distance (Å) between the residues and heme groups involved in the interface between chains A and D, related to the interactions with the studied ligands.
We observed that the major difference in distances from the crystallographic structure to Apo-Hem is related to the behavior of the carboxyl groups of the heme, which face the binding site cavity. This difference is due to the fact that in the crystallographic structure, the mentioned group is stretched, while in the Apo-Hem simulation, it is folded (Figure 3). In the case of Glufo-Hem, the interaction between the NH3+ group of the ligand and the carboxyl group of the heme in chain A stabilized the carboxyl, keeping it close to the ligand and away from H89(D) (Figure 1). Although there are no interactions between glufosinate ammonium and HEM(D), the presence of the ligand prevented further torsion of the carboxyl, keeping heme (D) distant from H96(A) (Figure 1). On the other hand, in the case of Fono-Hem, the ligand, through a π-stacking interaction with H89(D), stabilized this amino acid residue, keeping it more distant from the interaction with R72(A). Additionally, the presence of Fonofos in the binding site favored the stretching of the carboxyl of heme (D) through a water bridge (Figure 2).
Superposition of Lumbricus terrestris hemoglobin crystal structure (PDB ID 1X9F) and MD simulation of Glufo-Hem and Fono-Hem. Heme groups, F3 and E10 residues are shown. Atom colors are shown as follows: nitrogen in blue, oxygen in red, iron in red, and carbon in white (Crystal structure), green (Glufo-Hem), or pink (Fono-Hem). Hydrogens were omitted for clarity. Chains A and D are indicated with residue identification.
Our results agree with other studies on the relation of pesticides and hemoglobin. For instance, cartap hydrochloride, chlorpyrifos, cypermethrin, phenoxyherbicides and paraquat are known to interact with hemoglobin, affecting oxygen transport [23–27]. Regarding Glufosinate ammonium and fonofos, the first was banned from 29 countries and the latter from 41, but they are both still authorized in Brazil, for example [28]. Doroudian et al. (2022), studied tetraethyl Pyrophosphate (TEPP), an organophosphate, interactions with human hemoglobin, both in silico and in vitro. The referred work pointed out that TEPP induces hemoglobin aggregation and may cause structure changes, which may affect protein function. Although there is some information regarding the effects of pesticides on hemoglobin, there is a lack of studies with other taxonomic groups rather than vertebrates at the molecular level. In fact, bioindicators groups, such as Oligochaeta, could be used to assess pesticide risk in laboratory tests, since they are bioindicators intimately related to soil [29]. The approach we present here, based on molecular docking and molecular dynamics, may be a cost-effective pre-analysis on pesticide assessing risk that can be applied simultaneously on several bioindicator taxa, likely improving the detection of potential harming molecules.
This research not only establishes a basis but also advocates for the expanded application of computational methodologies to explore the potential harm which widespread substances may inflict on the natural environment. Moreover, it underscores the value of harnessing these computational tools for the intelligent design of chemicals that precisely aim at their intended targets, thereby reducing environmental degradation. This strategy exemplifies a forward-thinking approach to mitigating the ecological footprint of chemical agents, paving the way for more sustainable practices in substance design and application.
Conclusions
This study used molecular docking and dynamics simulations to evaluate the interaction of 52 organophosphate pesticides with Lumbricus terrestris hemoglobin. A conserved internal binding site was predicted, with structural features compatible with previously described allosteric cavities in other hemoglobins. Two compounds, glufosinate ammonium and fonofos, were selected for MD simulations based on their binding mode and interactions with residues close to the heme and interfacial regions.
Although the simulations were carried out using only two subunits of the dodecameric complex, they provided important insights into the local effects of ligand binding. Both ligands induced structural changes in residues involved in heme stabilization and at subunit interfaces, suggesting possible modulation of hemoglobin conformation. These results support the hypothesis that some organophosphates may interfere with the structural integrity of hemoglobin through noncovalent interactions.
Further studies should explore the effect of these compounds on oxygen affinity and structural stability using experimental approaches. Simulations involving the full hemoglobin complex or coarse-grained models may help to understand large-scale conformational changes and allosteric propagation. Given the ecological importance of L. terrestris, these findings may contribute to a better understanding of the molecular impacts of pesticides in soil invertebrates.
Our findings highlight the need for strategic design of molecules used in agricultural systems, leveraging established methodologies from Computational Chemistry and Biology. Moreover, Oligochaeta can serve not only as bioindicator models but also as in silico models for assessing the environmental impact of organophosphates and additional pesticides in terrestrial and aquatic ecosystems.
Notes
Acknowledgement
This work was supported by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES, Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP (process number 2018/12069-9) and Fundação Nacional de Desenvolvimento do Ensino Superior Particular - FUNADESP.
Conflict of interest
The authors declare that there is not any conflict of interest within this paper.
CRediT author statement
AMG: Conceptualization, Data Curation, Methodology, Writing-Review & Editing, Supervision; AVN: Data Curation, Validation, Writing-Original Preparation; MLO: Data Curation, Validation, Writing-Review & Editing; GRG: Supervision, Data Validation, Writing–Review & Editing.
Supplementary Material
This material is available online at www.eaht.org.