Network Toxicology and Molecular Docking to Investigate the non-AChE Mechanisms of Organophosphate-Induced Neurodevelopmental Toxicity

1. Introduction

The widespread use of organophosphate pesticides (OPs) in agriculture still raises concerns about their ecological impact and human health implication [1]. Surface waters and sediments can be contaminated by OPs from rainfall, runoff, improper disposal, or leaching from groundwater. This contamination can affect a wide range of aquatic organisms, such as fish, invertebrates, and algae, as well as non-target terrestrial organisms [2].

There are numerous toxic effects associated with OPs on fish, including hematological disturbances, gill, kidney, and liver alterations [3], oxidative stress, immune disorders, alteration of the intestinal microbiota [4], and behavioral disorders including reduced predator escape response [5]. In humans, exposure to OPs during pregnancy and postnatal periods increases the risk of autism spectrum disorder [6], impaired IQ, and verbal comprehension [7].

Exposure to OPs poses a high risk of neurotoxicity, especially for the developing nervous system. The effects of this type of exposure can include cognitive, motor, and memory impairments, as well as changes in brain morphology, referred to as Neurodevelopmental Toxicity (NDT) [8].

Acute exposure to high doses of OPs, such as self-poisoning, can significantly inhibit the activity of the acetylcholinesterase (AChE) enzyme. As a result, acetylcholine accumulates in the synaptic cleft, and cholinergic receptors in the nervous system and neuromuscular junction are overactivated, leading to a cholinergic crisis. The symptoms of human poisoning include respiratory depression, cardiovascular complications, gastrointestinal problems, and mental confusion, and the mortality rate is 15-30% [9].

When environmental contamination occurs, OPs concentrations are typically low and may not significantly inhibit AChE activity. However, these low concentrations can still cause behavioral alterations and neurological damage, particularly during developmental stages [10]. Therefore, it is important to emphasize that non-AChE mechanisms are also involved in the NDT of OPs [8]. These non-AChE mechanisms may include glutamatergic excitotoxicity, apoptosis [11], disturbances in intracellular Ca2+ homeostasis, neuronal activity [12], modulation of BDNF expression, and effects on the GABAergic and serotonergic systems [13].

Even though AChE inhibition is a common mechanism of action for OPs, there is evidence that the non-AChE effects can vary among different OP types. For instance, chlorpyrifos (CP) and malathion exhibited contrasting effects on zebrafish larvae. While malathion treatment led to morphological abnormalities in the brains of larvae, no such effect was observed with CP. Furthermore, CP-exposed larvae exhibited decreased swimming speed and increased resting time, whereas malathion-exposed larvae showed increased swimming speed and decreased resting time [14].

A comparative analysis of seven OPs (acephate, CP, dichlorvos, diazinon, malathion, parathion, and profenofos) in planarians revealed shared toxicological endpoints among the OPs (except for acephate and parathion), such as abnormal body forms, increased viscosity, writhing behavior, and alterations in swimming speed during exposure to blue light. However, other toxicological endpoints differed between the OPs [15].

This study was motivated by the following question: Do OPs share non-AChE mechanisms that contribute to NDT? This information can be valuable in enhancing toxicity testing protocols, establishing more precise safety margins, and providing targets for more effective treatments in cases of exposure to these compounds.

Computational tools can enhance toxicity studies by offering fast and comprehensive results. A valuable tool among these methods is network toxicology, a branch of network pharmacology. Utilizing techniques like protein-protein interactions, network topology analysis, and virtual prediction, it becomes possible to investigate the toxic effects of compounds based on "toxic characterization - compound - gene - protein" relationship. Integrating network toxicology with molecular docking, it is possible to assesses the interactions between the investigated compounds and their predicted targets, thereby helping to determine their toxicological properties [16, 17].

The aim of this study was to predict non-AChE targets that are shared among OP pesticides that contribute to NDT using network toxicology and molecular docking.

2. Materials and Methods

2.1. Chemical Structures of OPs

The 2D structures of the oxon metabolites of the OP pesticides diazinon (DZ), CP, and parathion were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov). We chose the oxon metabolites over of the original compounds due to their direct association with NDT [18, 19].

2.2. Acquisition of OPs Targets

The 2D chemical structures of paraoxon (PO), chlorpyrifos oxon (CPO), and diazinon-oxon (DO) were submitted to the Swiss Target Prediction database (http://www.swisstargetprediction.ch) and PharmMapper (http://www.lilab-ecust.cn/pharmmapper/) to predict target genes. We then entered the predicted target genes into the UniProt database (https://www.uniprot.org/) to retrieve the standard gene names, and duplicate entries were eliminated.

2.3. Prediction of NDT Targets

The GeneCards database (https://www.genecards.org/) and Online Mendelian Inheritance in Man (OMIM, http://omim.org/) were used to predict potential NDT targets. In February 2023, the related targets were collected using the keywords "neurodevelopmental abnormalities," "neurotoxicity," and "neurodevelopmental disorders," along with the Homo sapiens species. The UniProt database was utilized to retrieve the standard gene names, and after merging the predicted targets for the three keywords from both databases, duplicate entries were removed.

2.4. Venn Analysis and Construction of Protein-Protein Interaction (PPI) Network

In order to identify the intersection genes as potential targets of each OP related to NDT, a Venn diagram was created by using Draw Venn Diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Protein-protein interaction (PPI) networks were constructed using the STRING database (https://string-db.org/). A minimum score of 0.7 was set, indicating a high level of confidence, and only experimental data were selected as the source of active interaction. The PPI network analyzes direct interactions between the inputted proteins and generates a network of interactions where proteins are represented as nodes and the interactions between them are represented as edges.

2.5. Topological Analysis of PPI Networks

The PPI networks were imported into Cytoscape 3.9.1. The cytoHubba plugin was utilized to compute the top 10 hub nodes using the Maximal Clique Centrality (MCC) method. The Analyze Network tool was used for topological analyses of the network, including degree (k), betweenness centrality (BC), clustering coefficient, closeness centrality (CC), and average shortest pathway (ASPL).

The measure k represents the number of connections a node has with other nodes in the network, i.e., it quantifies the number of edges that connect to it. Nodes with a high k value are important for information transmission in the network as they have more control over the flow of information. On the other hand, BC is a measure that indicates a node's ability to connect other nodes in the network. A high BC value indicates that the node is important for communication between different parts of the network. In turn, the CC measure indicates how many direct neighbors of a node are interconnected, and a high CC value for a node indicates its involvement in a strongly connected cluster. Lastly, the ASPL measure indicates the average shortest path length between a node and all other nodes in the network [20].

2.6. Enrichment Analysis

An enrichment analysis was performed using the Consensus PathDB database (http://cpdb.molgen.mpg.de/) to explore Gene Ontology (GO) and Reactome Pathway. The Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) categories of GO were chosen. GO terms and Reactome Pathways with a p-value <0.01 were considered statistically significant. The top 15 GO terms for BP, MF, and CC, as well as the top 30 Reactome pathways, were selected for each OP. Additionally, terms common to all three OPs that were directly related to the Nervous System superpathway, as classified by Reactome.org, were examined. GO and Reactome pathway graphs were generated with the ggplot2 package in RStudio.

2.7. Molecular Docking:

A molecular docking analysis was conducted to explore the potential interactions between OPs and their hub nodes. The analysis followed the following steps:

a) Macromolecule preparation: The 3D structures of the proteins were obtained from the PDB database (https://www.rcsb.org/) in PDB format. Discovery Studio 2021 was used to remove water molecules and ligands, add polar hydrogen bonds, and obtain the x, y, z coordinates for constructing the grid box. Energy minimization was performed using Swiss-PDB Viewer software.

b) Ligand preparation: The 3D structures of the ligands were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) in SDF format and converted to PDB format using Open Babel.

c) Docking simulation was performed using Autodock Vina [21], and the analysis of the binding and generation of images was carried out using Discovery Studio 2021.

3. Results

3.1. Candidate Targets for OP-Induced NDT

After removing the duplicates, 388 targets were found for DZO, 343 for CPO, and 387 for PO in the Swiss Target Prediction and PharmMapper databases. Additionally, 3060 NDT-related targets were obtained from the Gene Cards and OMIM databases after removing the duplicates. Using a Venn diagram to identify the intersections of targets for each OP and NDT, we found 187, 121, and 181 targets for DZO, CPO, and PO, respectively. These targets were used to construct a PPI network using the STRING database, and the data were imported into Cytoscape for visualization and analysis of the hub nodes.

The analysis of PPI networks revealed 81 nodes and 113 edges for the DZO-NDT network, 37 nodes and 50 edges for the CPO-NDT network, and 71 nodes and 103 edges for the PO-NDT network (Figure 1). The top 10 hub nodes for each network were identified using the MCC algorithm in the cytoHubba plugin of Cytoscape. Topological measurements were obtained using the Network Analyzer tool. HSP90AA1, SRC, MET, and EGFR were hub nodes that were common among the three OPs. All hub nodes are listed in Table 1, along with their topological measurements. The measurements for the entire network are available in the supplementary material (Tables S1–S3).

3.2. GO and Reactome Pathway analysis

3.2.1. GO analysis

To analyze the predicted NDT targets for each OP, we conducted a GO analysis for BP, CC and MF (p-value < 0.01). We found 169 GO terms for DZO, with 102 in BP, 20 in MF, and 47 in CC. For CPO, we identified 143 GO terms, with 89 in BP, 18 in MF, and 36 in CC. PO had 155 GO terms, with 98 in BP, 18 in MF, and 39 in CC. For each OP, we selected the Top 15 processes for each term based on their p-values (Figure 2).

We observed that the top 15 GO terms for BP, MF, and CC were similar for DZO, CPO, and PO. The results indicated that the BP terms were mainly related to cellular response to stress, which could be a result of chemical exposure, as well as cellular processes such as communication, signal transduction, and cell death, along with processes involved in the development of an anatomical structure or organism.

Regarding the MF terms, we found that most targets for the three OPs were associated with binding, which refers to the interaction between one molecule and specific sites of another molecule. The results of the CC terms suggested most targets were found within the intracellular milieu, organelles, plasma membrane, and neurons.

3.2.2. Reactome Pathway Analysis

Based on Reactome Pathway analysis, 419 pathways were identified for DZO, 254 for CPO, and 358 for PO (p-value < 0.01). Figure 3 shows the top 30 Reactome pathways for each OP based on the p-value.

Reactome Pathway is a database that encompasses 2,546 pathways organized into 28 superpathways (e.g., signal transduction, disease, immune system, and developmental biology) and their subdivisions [22]. According to the analyses, signal transduction was identified as the most important superpathway for DZO, CPO, and PO, which included pathways related to signaling by receptor tyrosine kinases, MAPK family signaling cascades, and intracellular signaling by second messengers. Additionally, all three OPs shared the axon guidance pathway, which is directly related to nervous system development, suggesting that OPs interfere with neurogenesis.

We also evaluated the Reactome pathways for DZO, CPO and PO that showed significant values (p-value < 0.01) and were classified under the Neuronal system superpathway to identify the key processes in the nervous system affected by OP exposure. Activation of NMDA receptors and postsynaptic events was the common pathway for all three OPs, indicating the activation of N-methyl-D-aspartate receptors (NMDAR), an ionotropic glutamate receptor.

3.3. Molecular Docking

We used molecular docking to predict the interactions between each OP and their respective hub nodes (presented in Table 1). A lower binding energy between the compound and the target indicates a higher affinity. The binding energies between the OPs and macromolecules ranged from -7.2 to -4.2 kJ/mol (Figure 4). Among the targets shared by the OPs, HSP90AA1 had the lowest binding energy with the OPs, followed by MET, EGFR, and SRC. Various types of bonds were formed between the OPs and the amino acid residues of the targets, including hydrogen bonding and various pi interactions (e.g., Pi-Pi, Pi-sigma, and Pi-alkyl). It is noteworthy thar EGFR and SRC formed hydrogen bonds with all three OPs, suggesting a stronger interaction. Detailed results of the molecular docking analysis can be found in the supplementary material (Table S4 and S5 and Figure S1).

4. Discussion

The classic mechanism of neurotoxicity of organophosphates is the inhibition of the AChE activity, leading to the accumulation of the neurotransmitter acetylcholine in the synaptic cleft and consequently, overstimulation of cholinergic receptors. However, numerous studies have provided evidence that multiple mechanisms can contribute to the neurotoxicity of organophosphates, particularly during neurodevelopment [12].

Reactome pathway and Gene ontology enrichment analyzes were performed to identify potential mechanisms in the NTD of OP. The results indicate that these compounds impair signal transduction and axon guidance, disrupt cellular responses to stress, and activate NMDAR. In the OP-induced NTD PPI network, HSP90AA1, EGFR, MET, and SRC were identified as common hub nodes for DZO, CPO, and PO. For better understanding, the discussion will be divided into several topics according to the mechanism of toxicity.

4.1. Signal Transduction

Signal transduction is the process by which extracellular messengers bind to transmembrane receptors and provide information that triggers cellular responses, such as biochemical and biological transformations, as well as gene expression alterations [23]. According to Reactome pathway enrichment analyses, the major signal transduction pathways affected by organophosphates include the receptor Tyrosine Kinase (RTK) pathway, the mitogen-activated protein kinase (MAPK) pathway and signaling by second messengers. It is important to note that these pathways are interconnected, since RTK activation triggers phosphorylation and recruitment of effector proteins, initiating signaling cascades involving MAPKs, PI3K/Akt, phospholipase C-PKC, and STAT [24].

In vitro studies have explored how OP exposure affects signaling pathways. For instance, exposure of cell cultures to CP and monocrotophos activated the MAPK ERK1/2, JNK, and p38 pathways, thereby inducing apoptosis. Activation of these pathways were associated with reactive oxygen species (ROS) generation and oxidative stress, which may contribute to mitochondrial damage and cell death [25,26,27]. It is plausible that these pathways may also be influenced by the interaction of OPs with RTKs.

Based on our literature review, no experimental data were found indicating that OPs cause neurodevelopmental abnormalities via EGFR, MET, or non-receptor tyrosine kinase protein SRC. However, CP has previously been shown to increase migration and invasion of breast cancer cells through c-SRC pathway activation, amplifying downstream AKT and p-38 signaling [28]. Additionally, CP has promoted the growth of human colorectal adenocarcinoma H508 cells through increased EGFR/ERK1/2 signaling [29]. Although activation of the EGFR/ERK1/2 pathway generally stimulates cell growth, it can also promote neuronal cell death [30, 31].

4.2. Axon Guidance

Axon guidance is the process by which neurons send their axons toward their correct targets for synaptic formation, and it is crucial during neural circuit development. Growing axons have a structure at their tips known as the growth cone, which contains receptors that respond to extracellular messengers that attract or repel the axon, guiding them to their destination [32]. The hub nodes EGFR, MET, and SRC play important roles in axon guidance.

The EGFR (or ErbB1) belongs to the ErbB family, which includes three other members: ErbB2, ErbB3, and ErbB4. Activation of EGFR initiates signaling pathways that are important for neuronal proliferation, differentiation, migration, as well as neural circuit development, including the Ras-Raf-MEK-ERK1/2, STAT3, and PI3K-Akt-mTOR pathways [33]. Additionally, activated EGFR interacts with SRC, forming a complex that enhances EGFR phosphorylation, consequently amplifying downstream signaling [34]. SRC also interacts with other signaling pathways involved in axon guidance, such as the Ephrin receptor (EphA) pathways [35], Sonic hedgehog (Shh) pathway [36], Netrin-1 pathway [37], and MET pathway [38]. Similar to EGFR, the activation of MET triggers downstream signaling pathways of PI3K/Akt, Ras-Raf-MEK-ERK, and STAT3, which are crucial for neurodevelopment and synaptogenesis [38]

Axonal growth is a fundamental process in the formation of neural circuits and the development of neuronal connectivity. Previous studies have reported impairment of axonal growth by OPs through noncholinergic mechanisms. For example, low concentrations of CP and CPO inhibited axonal outgrowth in cultures of rat sympathetic neurons [39]. Moreover, exposure to CPO during neurodevelopment disrupted axonal outgrowth in zebrafish [40]. These results showed that CP and CPO did not cause axon retraction or a decrease in the number of axons, but rather reduced growth rate compared to the untreated control group. Furthermore, no morphological changes in the axonal growth cone were observed. Experimental data suggest that OPs altered axonal growth signaling [39], similar to the findings in our study that OPs impair axon guidance.

4.3. Cellular Response to Stress

Exposure to environmental contaminants during prenatal development can trigger cellular response to stress that is closely related to heat shock and oxidative stress systems [41]. The HSP90AA1 gene encodes the Heat shock protein (Hsps) HSP90, which acts as a chaperone to refold denatured proteins in the presence of stressors such as heat, oxidative stress, hypoxia, and exposure to cytotoxic agents. In the absence of stress, chaperones facilitate protein refolding, stabilization, translocation, and degradation, preventing protein aggregation. In neurodevelopment, Hsps regulate pathways related to cell growth and migration, such as the PI3K/Akt signaling pathway, and they are crucial mediators for axon guidance [42].

Cellular stress increases the expression of HSP90 and diverts its functions to cope with stress responses. This functional diversion or inhibition of HSP90 can alter the activity of various proteins and result in increased mutations, leading to neurodevelopmental disorders [43]. Increased expression of HSP90 has been observed in the liver, muscle, kidneys, and spleen of fish exposed to CP, indicating that OPs exposure can induce cellular response to stress [44, 45]. Although these data were not obtained from nervous system cells, common carp exposed to CP showed increased expression of HSP90 and hypoactivity, which may indicate neurotoxicity [44].

There are several ways cells respond to stress, including activating of survival pathways or apoptosis [46]. There is still much to be clarified about the pathways by which OPs induce apoptosis, but they are known to occur during vulnerable windows of neurodevelopment and is likely not related to AChE inhibition [47]. However, it has been observed that the AKT pathway and oxidative stress are mechanisms by which OPs induce neural apoptosis, but independently of each other [48]. Therefore, the relationship between OPs and HSP90 may be a pathway to understanding how these processes occur.

4.5. Activation of NMDAR

The ionotropic glutamate receptors N-Methyl-D-Aspartate receptors (NMDAR) are calcium-permeable channels. The synaptic activity of NMDAR contributes to axonal and dendritic growth, as well as the maturation of glutamatergic synapses [49]. These processes are crucial for the development of cognitive functions, neuronal plasticity, memory, and learning [50]. Excessive activation of NMDAR can increase calcium influx, leading to cellular damage and neuronal death, likely mediated by the PKC/ERK pathway [51, 52].

During exposure to high concentrations of OPs, ionotropic AMPA and NMDA glutamate receptors are activated due to increased excitatory signaling following AChE inhibition and excessive stimulation of muscarinic receptors, leading to seizures and neuronal death [53]. However, the mechanisms related to low concentrations of these chemicals have not been fully elucidated yet.

Studies indicate that glutamate-mediated excitotoxicity seems to be more related to the toxicity effects of CP than DZ. An NMDAR antagonist, MK-801, attenuated CP toxicity, but not DZ toxicity, in cortical neuronal and glial cultures [11]. Furthermore, in neonatal rats exposed to CP and DZ, both substances upregulated the expression of NMDAR subunits, but CP induced more significant changes than DZ [54]

In our study, we did not assess the direct activation of NMDAR by OPs, nor the genes encoding the subunits of this receptor. However, the activity of NMDAR can be regulated by SRC [55], EGFR [56], and MET [57], which have been identified as hub nodes in the PPI network of OPs. Furthermore, studies with chronic exposure to low doses of OPs is a more realistic way to assess environmental exposure, and further studies are needed to evaluate the relationship between this exposure and NMDAR activity, considering other biomarkers beyond the genes encoding the subunits of this receptor.

5. Conclusions

Using network toxicology and molecular docking analyses, this study identified the hub nodes HSP90AA1, EGFR, MET, and SRC as potential targets for OPs. Additionally, signal transduction, axon guidance, cellular response to stress, and activation of NMDAR were found to be key pathways involved in OP-induced NDT. These findings are important for understanding the mechanisms of neurotoxicity of these substances at environmentally relevant concentrations that do not involve cholinergic pathways.