Identification of differentially expressed genes and pathways in diquat and paraquat poisoning using bioinformatics analysis

Abstract

Objective

In this study, differentially expressed genes (DEGs) and signaling pathways involved in diquat (DQ) and paraquat (PQ) poisoning were identified via bioinformatics analysis, in order to inform the development of novel clinical treatments.

Methods

Raw data from GSE153959 were downloaded from the Gene Expression Omnibus database. DEGs of the DQ vs. control (CON) and PQ vs. CON comparison groups were identified using R, and DEGs shared by the two groups were identified using TBtools. Subsequently, the shared DEGs were searched in the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, using the Database for Annotation, Visualization, and Integrated Discovery. A protein–protein interaction (PPI) network was constructed, and hub genes were identified using the cytoHubba plug-in in Cytoscape software. Finally, circos and contrast plots showing the DEGs shared between mouse and human chromosomes were constructed using TBtools.

Results

Thirty-one DEGs shared by the DQ and PQ groups were identified. Enriched biological process terms included positive regulation of cell proliferation and translation. Enriched cellular component terms included extracellular region, intracellular membrane-bounded organelle and mitochondrion. Enriched molecular function terms included transcription factor activity and sequence-specific double-stranded DNA binding. Enriched KEGG pathways included the interleukin-17 signaling pathway, tumor necrosis factor signaling pathway, and human T-cell leukemia virus 1 infection. The top 10 hub genes in the PPI network were prostaglandin-endoperoxide synthase 2 (Ptgs2), chemokine (C-X-C motif) ligand 2 (Cxcl2), colony-stimulating factor 2 (granulocyte-macrophage) (Csf2), matrix metallopeptidase 13 (Mmp13), amphiregulin (Areg), plasminogen activator, urokinase receptor (Plaur), fos-like antigen 1 (Fosl1), epiregulin (Ereg), activating transcription factor 3 (Atf3), and transferrin receptor (Tfrc). Cxcl2, Csf2, and Atf3 played important roles in the mitogen-activated protein kinase (MAPK) signaling pathway.

Conclusions

These pathways and DEGs may serve as targets for gene therapy.

Keywords:

• Diquat
• Paraquat
• Ptgs2
• Cxcl2
• MAPK

1. Introduction

Diquat (DQ) (Magalhães et al. 2018) and paraquat (PQ) (Mitra et al. 2011) are nonselective and widely used herbicides. DQ (1,1′-ethylene-2,2′-bipyridinium ion) (Magalhães et al. 2018; Jones and Vale 2000) is an analog of PQ (1,1′-Dimethyl-4,4′-bipyridinium dichloride) (Mitra et al. 2011). The 50% of lethal dose (LD₅₀) of PQ in humans is 20–40 mg ion/kg of body weight (Vale et al. 1987); In rats, the LD₅₀ of oral PQ is 150 mg/kg (World Health and International Programme on Chemical 2004), whereas, that of oral DQ is 215-235 mg/kg (Magalhães et al. 2018). Therefore, DQ is considered less toxic than PQ (Fortenberry et al. 2016). PQ and DQ poisoning remain fairly common in developing countries (Xing et al. 2020; Lu et al. 2021).

Inhalation, dermal contact, ingestion, and eye contact with PQ can cause organ failure, pulmonary fibrosis, liver injury, kidney injury, and even death. In survivors, PQ poisoning can lead to Parkinson’s disease (PD) or Alzheimer’s disease (AD) (Singh et al. 2018). PQ poisoning most severely affects the lungs, and most deaths are attributed to pulmonary fibrosis and respiratory failure. DQ is mainly metabolized by the kidneys. Therefore, DQ poisoning can induce renal failure (Jones and Vale 2000). DQ can also damage the liver, heart, gastrointestinal tract, skin, and central nervous system (Magalhães et al. 2018; Jones and Vale 2000), and lead to death (Xing et al. 2020). In contrast to PQ poisoning, pulmonary fibrosis has not been reported for DQ poisoning (McCarthy and Speth 1983; Magalhães et al. 2018; Jones and Vale 2000). Both of them mostly adopt comprehensive treatment, and no specific antidote exists for either type of poisoning (Jones and Vale 2000; Sun and He 2017). In this study, based on the similar structures of PQ and DQ, shared differentially expressed genes (DEGs) and signaling pathways were explored using GSE153959 data, to identify potential treatment targets.

The effects of DQ and PQ on transforming growth factor-alpha (TGF-alpha) transgenic mouse hepatocytes (TAMHs) (Davis and Stamper 2016) in vitro were evaluated, and mitogen-activated protein kinase (MAPK) pathway microarray data (Stamper et al. 2010) were analyzed to identify shared DEGs and signaling pathways to provide a basis for gene therapy.

2. Materials and methods

2.1. Affymetrix microarray data

Raw data from GSE153959, contributed by Stamper BD and Nguyen NU, were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), using the GPL17400, [MoGene-2_1-st] Affymetrix Mouse Gene 2.1 ST Array [transcript (gene)] version. RNA was isolated from TAMHs cultured with 18 µM DQ, 126 µM PQ, or control culture media at the 9-hour timepoints. Three groups (three samples per group) were established: control 9_h (CON), diquat (18 µM) 9_h (DQ), and paraquat (126 µM) 9_h (PQ).

2.2. Identification of DEGs

Transcriptome Analysis Console (TAC) 4.0.2 (Affymetrix, Santa Clara, CA, USA) was used to construct a list of DEGs. Three comparisons were established: DQ vs. CON, PQ vs. CON, and DQ vs. PQ. A fold change (FC) ≤ −2 or ≥ 2 and a p. value ≤ 0.05 were selected as the threshold for identifying DEGs. The DEGs of the three groups were identified and visualized using R (version 4.1.2, ggplot2 package, Windows environment) after deleting the incomplete genes.

2.3. Analysis of DEGs in the intersection between DQ vs. CON and PQ vs. CON

The DEGs in the intersection between the DQ vs. CON and PQ vs. CON comparison groups and their expression were shown in the venn diagram and the heat map created with TBtools (Chen et al. 2020) (version 1.0986; https://github.com/CJ-Chen/TBtools/releases).

2.4. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis of DEGs in the intersection between groups

The Database for Annotation, Visualization, and Integrated Discovery (DAVID; version 6.8; https://david.abcc.ncifcrf.gov) was used to analyze the DEGs for GO terms and KEGG pathway enrichment. The thresholds for enriched terms were established as ≥ 2 genes and a p. value ≤ 0.1.

2.5. Protein-protein interaction (PPI) network and module analysis

The interactions of the DEGs were predicted by the Search Tool for the Retrieval of Interacting Genes with a link confidence ≥ 0.04 (STRING, version 11.5; https://string-db.org). The PPI network was constructed using Cytoscape (version 3.8.2, Java 11.0.6), and the cytoHubba plug-in (version 0.1) was used to identify hub genes.

2.6. Circos plot of mouse chromosomes and contrast plot of mouse and human chromosomes

A circos plot was constructed to show the chromosomal location of the mouse genes in the PPI. A contrast plot was drawn to reveal the distribution of these genes on human chromosomes. The circos plot was drawn using TBtools (Chen et al. 2020) (version 1.0986; https://github.com/CJ-Chen/TBtools/releases) and MCScanX (Wang et al. 2012). The annotation files and genome sequence files for mice (GRCm39) and humans (GRCh38.p13) were downloaded from https://www.gencodegenes.org.

3. Results

3.1. DEGs of the DQ vs. CON and PQ vs. CON comparison groups

In the DQ vs. CON comparison, 52 DEGs were identified, including 47 upregulated and five downregulated genes, and the volcano plot was shown in Figure 1(A). A total of 43 DEGs were identified in the PQ vs. CON comparison, including 36 upregulated and seven downregulated genes, and the volcano plot was shown in Figure 1(B). No DEGs were identified in the DQ vs. PQ comparison.

Figure 1. Volcano plots of differentially expressed genes (DEGs). (A) DQ vs. CON and (B) PQ vs. CON. Gray, non-differentially expressed genes; red, significantly upregulated genes; blue, significantly downregulated genes (based on a fold change ≤ -2 or ≥ 2 and a p. value ≤ 0.05). FC, fold change.

3.2. Intersection of DEGs between DQ vs. CON and PQ vs. CON

A total of 35 DEGs were shared between the two comparison groups, out of which, four duplicate genes were eliminated. The remaining 31 DEGs included 27 upregulated and four downregulated genes. The Venn diagram was presented in Figure 2(A), and the expression of the 31 DEGs was illustrated in Figure 2(B).

Figure 2. (A) Venn diagram of intersecting DEGs. (B) Heat map of intersecting DEGs.

3.3. Go and KEGG pathway enrichment analysis of intersecting DEGs

28 biological process (BP) terms, six cellular components (CC) terms, and seven molecular function (MF) terms, and, seven KEGG-pathways were enriched. The enriched BP terms included translational elongation, response to glucocorticoid, translation, positive regulation of cell proliferation, response to estradiol, ovulation, negative regulation of calcium ion transport, positive regulation of epidermal growth factor- activated receptor activity, response to vitamin D, and regulation of cell proliferation. The enriched CC terms included extracellular region, intracellular membrane-bounded organelle, extracellular space, mitochondrion, cell surface and brush border membrane. The enriched MF terms included triplet codon-amino acid adaptor activity, growth factor activity, transcriptional activator activity, RNA polymerase II transcription regulatory region sequence-specific binding, sequence-specific double-stranded DNA binding, cytokine activity, transcription factor activity, sequence-specific DNA binding, and epidermal growth factor receptor binding (Figure 3(A)). The enriched KEGG pathways included interleukin (IL)-17 signaling pathway, aminoacyl-tRNA biosynthesis, tumor necrosis factor (TNF) signaling pathway, maturity onset diabetes of the young, Kaposi sarcoma-associated herpesvirus infection, ferroptosis and human T-cell leukemia virus 1 infection (Figure 3(B)).

Figure 3. GO and KEGG enrichment analysis of intersecting DEGs (top 10 hub genes, enriched terms ≥ 2 genes and a p. value ≤ 0.1). (A) GO analysis: Biological Process (BP) (green points), Cellular Component (CC) (red points), and Molecular Function (MF) (purple points). (B) KEGG pathways (green points). Point size indicates the p. value of each point.

3.4. PPI network analysis of intersecting DEGs

A PPI network of 12 nodes and 17 edges was constructed using Cytoscape software with link confidence ≥ 0.04 (Figure 4). The cytoHubba plug-in was used to construct an MCC algorithm for identifying hub genes. The top ten hub genes were prostaglandin-endoperoxide synthase 2 (Ptgs2), chemokine (C-X-C motif) ligand 2 (Cxcl2), colony-stimulating factor 2 (granulocyte-macrophage) (Csf2), matrix metallopeptidase 13 (Mmp13), amphiregulin (Areg), plasminogen activator, urokinase receptor (Plaur), fos-like antigen 1 (Fosl1), epiregulin (Ereg), activating transcription factor 3 (Atf3), and transferrin receptor protein (Tfrc). All of these genes were upregulated.

Figure 4. PPI network of the intersecting of DEGs. Red, upregulated genes; green, downregulated genes. Node size was positively related to the source node attribute, edge color and width were positively related to the combined score.

3.5. Circos and contrast plots of mouse and human chromosomes

The Circos plot revealed that Ptgs2 and Atf3 were located on chromosome (chr) 1; Cxcl2, Ereg and Areg on chr5; Csf2 on chr11; Mmp13 on chr9; Plaur on chr7; Fosl1 on chr19; Tfrc on chr16; and Slc2a2 and Hmgcs2 on chr3 (Figure 5(A)). The contrast diagram illustrated the sequence alignment between the mouse and human genomes, with 11 genes exhibiting collinearity. The aligned sequences were marked in red in Figure 5(B), and the details of the sequence alignment were presented in Table 1.

Figure 5. Distribution of genes on chromosomes. (A) Circos plot: distribution of genes on mouse chromosomes. Gene symbol-number corresponds to different transcripts. Gray lines, collinearity of all mouse chromosomes; colored lines, collinearity of chromosomes of the 12 mouse genes. (B) Contrast plot of mouse and human genes: comparative distribution of the genes on the mouse and human chromosomes. Gray lines, collinearity of all mouse and human genes; red lines, collinearity of 11 mouse and human genes; orange bars, mouse chromosomes; green bars, human chromosomes.

Table 1. Comparative distribution of genes on mouse and human chromosomes.

Transcript name	Position on mouse chromosome	Transcript on mouse	Position on human chromosome	Transcript on human
Csf2-201	chr11	ENSMUST00000019060.6	chr5	ENST00000296871.4
Tfrc-201	chr16	ENSMUST00000023486.15	Chr3	ENST00000360110.9
Slc2a2-201	chr3	ENSMUST00000029240.14	Chr3	ENST00000314251.8
Ereg-210	chr5	ENSMUST00000031324.6	Chr4	ENST00000244869.3
Areg-201	chr5	ENSMUST00000031325.6	Chr4	ENST00000502307.1
Ptgs2-201	chr1	ENSMUST00000035065.9	Chr1	ENST00000367468.10
Atf3-201	chr1	ENSMUST00000027941.14	Chr1	ENST00000366987.6
Hmgcs2-202	chr3	ENSMUST00000120541.8	Chr1	ENST00000369406.8
Fosl1-201	chr19	ENSMUST00000025850.7	Chr11	ENST00000312562.7
Mmp13-201	chr9	ENSMUST00000015394.10	Chr11	ENST00000236826.8
Plaur-201	chr7	ENSMUST00000002284.11	Chr19	ENST00000340093.8

Annotation: gene symbol number corresponds to different transcripts; chr, chromosome

4. Discussion

In this study, we analyzed intersecting DEGs between the DQ vs. CON and PQ vs. CON comparison groups in TAMHs in vitro. Based on |FC| ≥ 2 and a p. value ≤ 0.05, 52 DEGs in the DQ vs. CON group (47 upregulated and five downregulated) and 43 DEGs in the PQ vs. CON group (36 upregulated and seven downregulated) were identified. Finally, 31 common DEGs were screened (27 upregulated and four downregulated).

The enriched BP terms were mainly cell proliferation and translation; the enriched CC terms were the integrity of the plasma, organelle membranes, and mitochondrion; and the enriched MF terms were mainly transcription factor activity, sequence-specific DNA binding, cytokine and growth factor activity. These GO terms were not specific to any particular organ, indicating that the toxicity of DQ and PQ can affect cells of all organs. A previous study (Huang et al. 2019) showed that PQ could reduce the level of DNA methylation, inhibiting the proliferation of human neural stem cells. The enriched KEGG pathways included the IL-17 signaling pathway, TNF signaling pathway, and human T-cell leukemia virus 1 infection (Cao et al. 2015; Zhu et al. 2018).

Ptgs2 participates in multiple biological processes, including inflammation (Jaén et al. 2018), colorectal cancer (Benelli et al. 2018), and cardiovascular disease (Mitchell and Kirkby 2019). Ptgs1 and Ptgs2 are iso forms of the cyclooxygenase enzyme. Ptgs1 is expressed at a low level in vivo and mainly plays a homeostatic role, whereas Ptgs2 is expressed in response to various challenges and stress conditions. Meloxicam, a Ptgs2 inhibitor, can block apoptosis and cause mitochondrial dysfunction in cells treated with DQ (Choi et al. 2018). PQ may be toxic to dopaminergic neurons through Ptgs2 (Yang et al. 2010), and Ptgs2 deletion attenuates PQ-induced motor impairment and dopamine utilization in mice (Litteljohn et al. 2008). PQ poisoning causes lung edema and hemorrhage, liver congestion, glycogen degeneration, and proximal tubular degeneration in the kidneys of rats, and Ptgs2 expression is upregulated in affected organs (Malekinejad et al. 2013).

Cxcl2 plays a crucial role in neutrophil recruitment (Kobayashi 2008), and maintains neutrophil homeostasis (Day and Link 2012; Lou et al. 2017). Csf2 is produced by the oviduct and endometrium (Tríbulo et al. 2018) and regulates embryonic development during the preimplantation phase (Hansen et al. 2014; Shyam et al. 2020). Atf3 is a member of the ATF family, representing a large group of basic-region leucine zipper (bZIP) transcription factors (Thompson et al. 2009). Atf3 is an adaptive-response gene involved in cellular processes that adapts to changes in the body (Hai et al. 2011), and can regulate inflammation, cancer, and myocardial ischemia (Thompson et al. 2009; Yang et al. 2016). Cxcl2, Csf2, and Atf3 are target genes of the MAPK8 (JNK1), MPK14 (p38-alpha) and MAP3K1 (MEKK1) reactions in TAMHs in an in vitro hepatocellular model treated with 3′-hydroxyacetanilide (AMAP) and acetaminophen (APAP) for 2 h (Stamper et al. 2010). Mitogen-activated protein kinases (MAPKs) are members of a three-kinase phosphorylation system expressed in all cells, composed of MAPK, MAPK kinase (MKK), and MAPK kinase kinase (MKKK). MKKKs phosphorylate and activate MKKs, and MKKs phosphorylate and activate MAPKs. The MAPK family mainly contains ERK, JNK, and p38. MAPK pathways monitor cell growth, toxicant metabolism, and the oxidative stress response (Johnson et al. 2005). PQ induces significant phosphorylation of ERK, JNK, and p38 in murine BV₂ microglia cells, producing neuroinflammation and neurotoxicity (Wang et al. 2020). PQ can also cause neuronal cell apoptosis through the JNK pathway in SH-SY5Y cells (Ju et al. 2019) and can damage the lungs through the p38-MAPK pathway (Cui et al. 2019; SreeHarsha 2020). However, few studies on DQ-induced cell death via the MAPK pathway have been conducted. MAPK14 positively regulates Csf2, and while MAPK14 and MAPK13 are known to regulate Cxcl2 and Atf3, respectively, the direction of regulation is unknown. The circos plot illustrated the distribution of DEGs on different transcripts of mouse chromosomes. Based on our comparison of the DEGs on mouse chromosomes with the corresponding genes on human chromosomes, we hypothesize that these intersection DEGs play an important role in DQ or PQ poisoning in humans. Future studies of these DEGs in humans are needed to assess potential therapeutic target genes and pathways. A limitation of this study was the lack of a basic test for DQ or PQ toxicity; tests for toxicity in human- derived liver cells treated with DQ or PQ should be developed, and blood samples from patients with DQ or PQ poisoning should be collected for further studies.

Author contributions

Changqing Miao arranged the datasets and collected the relevant literature; Dandan Fan designed the research, drew the diagrams, and wrote and revised the manuscript. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Ethical approval was obtained from the Biomedical Ethics Committee of the Medical department of Xi’an Jiaotong University (2016-135).

Acknowledgment

The authors acknowledge the institutional support from the First Affiliated Hospital of Xi’an JiaoTong University and the hospital library. We thank everyone who participated in the study, those involved in the study design, and those who provided guidance during the study. We appreciate the public data provided by GEO database.

Disclosure statement

The authors declare no conflict of interest.

Additional information

Funding

The study was supported by the National Natural Science Foundation of China [81601674].