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).

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.

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).

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.

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.

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).

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.

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.