Network pharmacology- and molecular docking-based approaches to unveil the pharmacological mechanisms of dihydroartemisinin against esophageal carcinoma

Objective: Dihydroartemisinin (DHA) is an active metabolite of artemisinin and its derivatives, which is a potent drug extensively applied in clinical treatment of malaria. The antitumor properties of DHA have received increasing attention. However, there is no systematic summary on the pharmacological mechanisms of DHA against esophageal carcinoma (ESCA). The present study implemented network pharmacology- and molecular docking-based approaches to unveil the pharmacological mechanisms of DHA against ESCA.

Methods: DHA targets were accessed through integrating the SwissTargetPrediction, HERB, as well as BATMAN-TCM platforms. In TCGA-ESCA dataset, genes with differential expression were screened between 161 ESCA and 11 normal tissue specimens. DHA targets against ESCA were obtained through intersection. Their biological significance was evaluated with functional enrichment analysis. A prognostic signature was established via uni- and multivariate cox regression analyses. DHA-target interactions were predicted via molecular docking. Molecular dynamics simulation was implemented to examine the stability of DHA binding to potential targets. Results: The study predicted 160 DHA targets as well as 821 genes with differential expression in ESCA. Afterwards, 16 DHA targets against ESCA were obtained, which remarkably correlated to cell cycle progression. The ADORA2B- and AURKA-based prognostic signature exhibited the reliability and independency in survival prediction. The stable docking of DHA-ADORA2B and DHA-AURKA was confirmed.

Conclusion: Collectively, this study systematically revealed the basis and mechanism of DHA against ESCA through targeting multi-target and multi-pathway mechanisms, and thus offered theoretical and scientific basis for the clinical application of DHA.

Introduction

Esophageal carcinoma (ESCA) remains one of the most lethal cancers across the globe (Sung et al., 2021), with two primary histological subtypes: squamous cell carcinoma (ESCC) as well as adenocarcinoma (EAC) (Lu et al., 2022). The overall 5-year survival rate is merely 15–25% (Hulshof et al., 2021). The features and etiology of ESCA may vary on the basis of region or ethnicity (Liu et al., 2017). ESCC occupies the predominant subtype of ESCA, notably in the Asian and African regions, which correlates to dietary habits as well as exposure to carcinogens (Zhang et al., 2022). Oppositely, EAC mainly develops from esophageal epithelial intestinal metaplasia caused by chronic gastroesophageal reflux diseases, which is the majority of ESCA in Western countries (Zhang et al., 2022). Currently, endoscopic or surgical resection is appropriate for early-stage patients, while radio- or/and chemotherapies are mainly reserved for those with advanced stage or metastases (Mukherjee et al., 2021; Yang et al., 2022). Chemotherapeutic agents (cisplatin, 5-fluorouracil, doxorubicin, etc.) may cause a range of dose-limiting toxicities (Zhang et al., 2021). Molecular targeted therapy has achieved unprecedented progress in cancer therapy (Lu et al., 2022). Nevertheless, only trastuzumab against HER2, as well as ramucirumab targeting VEGFR-2 have been recommended for the treatment of ESCA (Xu et al., 2021; Safran et al., 2022). Immunotherapy comprising immune checkpoint inhibitors, vaccines, monoclonal antibodies as well as adoptive cellular immunotherapy represents a novel therapeutic option of ESCA (Janjigian et al., 2021; Shitara et al., 2021). Despite the remarkable potential of immunotherapy in cancer therapy, individual patients’ response greatly varies, and only a small percentage of cases respond to immunotherapy (Luo et al., 2021; Sun et al., 2021). Hence, achievement of the goal of effective therapeutic options remains challenging. Rather than exploring a single ESCA-causing gene and therapeutic agents that act only on a single target, the entire drug-disease network should be considered, aiming to determine multi-target agents to lower side effects.

Traditional Chinese medicine artemisinin is extracted from Artemisia annua L, and dihydroartemisinin (DHA) with high water solubility and antimalarial activity is the first-generation derivative of this compound that is an effective and fast-acting antimalarial agent with low toxicity (Zhu et al., 2020). DHA displays the considerable potential for preventing or treating ESCA because of its low toxicity and known safety, and preclinical research has proposed the molecular mechanisms and pharmacological effect underlying DHA against ESCA (Zhu et al., 2020). Zhao et al. reported that DHA exhibits an anti-proliferative effect against ESCC cells, which is capable of down-regulating mTOR cascade pathway partly via binding to AKT1 and p70S6K (Zhu et al., 2020). Pyroptosis is a novel form of pro-inflammatory programmed cell death, and its inducers enable to strengthen antitumor effects (Niu et al., 2022). Song et al. proposed that DHA may trigger pyroptosis of ESCC cells via impacting the activity of PKM2-caspase-8/3-GSDME signaling (Jiang et al., 2021). Administration of chemotherapeutic agents is usually accompanied by resistance. DHA may sensitize ESCC to cisplatin through attenuating the activity of Sonic Hedgehog pathway (Cui et al., 2020). There is an interplay between autophagy and epithelial-mesenchymal transition during tumor progression (Chen et al., 2021; Niu et al., 2021). He at al. found that the migratory capacity and the epithelial-mesenchymal transition process of ESCA cells is mitigated by DHA via triggering autophagy (Chen et al., 2020). PKM2 is a key regulator of glycolysis, and DHA lowers glycolysis of ESCA through down-regulating PKM2 (Li et al., 2019). Additionally, hTERT has been confirmed as a therapeutic target of DHA against ESCC (Li et al., 2021). Altogether, DHA exerts an anti-ESCA effect through targeting multi-target and multi-pathway mechanisms.

In the current study, we developed a new integrated strategy to probe out the key targets and mechanisms of DHA in the treatment of ESCA on the basis of network pharmacology and molecular docking, which offered the theoretical and scientific basis for the clinical application of DHA against ESCA.

Materials and methods

Dihydroartemisinin target prediction

The chemical structure of DHA was acquired from the PubChem chemical information resource (https://pubchem.ncbi.nlm.nih.gov) (Kim et al., 2021). Potential targets of DHA were searched from the SwissTargetPrediction (http://www.swisstargetprediction.ch) (Daina et al., 2019), as well as HERB (http://herb.ac.cn/) (Fang et al., 2021) web servers. Additionally, DHA targets were screened through the BATMAN-TCM (http://bionet.ncpsb.org/batman-tcm) platform in accordance with the screening conditions of score cutoff = 5 and adjusted p-value < 0.05 (Liu et al., 2016). DHA targets were matched and corrected on the basis of the universal Protein Resource (UniProt, (2021); http://sparql.uniprot.org/). Afterwards, DHA targets obtained from three platforms were merged, and deduplicated.

The cancer genome atlas (TCGA)-esophageal carcinoma data acquisition

RNA transcriptome profiles as well as clinicopathological and survival information of ESCA cases were downloaded from TCGA-ESCA dataset (https://cancergenome.nih.gov/). In total, 161 ESCA and 11 normal tissue specimens were enrolled for subsequent analysis.

Differential expression analysis

In TCGA-ESCA dataset, limma package was adopted for screening differentially expressed genes (DEGs) between ESCA and normal tissue specimens. Adjusted p-value < 0.01 together with |log2 fold-change|>1 were set as the screening criteria. Shared genes between DHA targets and DEGs were intersected through the Venny 2.1 online tool, named DHA targets against ESCA.

Protein-protein interaction (PPI) analysis

DHA targets against ESCA were uploaded to the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) (Szklarczyk et al., 2015). In accordance with hide unconnected targets along with highest confidence of 0.900, PPI network of DHA targets against ESCA was obtained. Hub genes among DHA targets against ESCA were further evaluated with cytoHubba plugin of Cytoscape tool 3.7.2 (Chin et al., 2014). Degree of each hub gene was also computed.

Establishment and verification of a prognostic signature

In TCGA-ESCA dataset, univariate cox regression analysis was implemented to screen which DHA targets against ESCA significantly correlated to ESCA cases’ overall survival (OS). Prognostic DHA targets against ESCA with p-value <0.05 were enrolled in a multivariate cox regression model. The prognostic signature was conducted on the basis of the linear combination of regression coefficients derived from the multivariate Cox regression model multiplied with expression values. The risk score of each case was computed. TCGA-ESCA cases were equally separated into discovery set as well as validation set. In each set, ESCA cases were further assigned to high- and low-risk subpopulations. Kaplan-Meier curves of OS between two subpopulations were plotted, and OS difference was estimated with log-rank test. Alive and dead status of cases in each set was displayed. Receiver operating characteristic (ROC) curve was drawn for evaluation of the prediction accuracy of the prognostic signature. Area under the curve (AUC) was also computed. TNMplot.com web tool was employed for comparison of survival difference between high and low expression subpopulations of DHA targets against ESCA in TCGA-ESCA dataset (Bartha and Győrffy, 2021). Expression values of DHA targets against ESCA across distinct pathological stages were displayed with ggpubr package, and the difference was evaluated using one-way analysis of variance.

Uni- and multivariate cox regression analysis

Clinicopathological factors comprising gender, pathological stage, TNM staging as well as prognostic signature were enrolled for estimating their associations with TCGA-ESCA cases’ OS by use of uni- and multivariate cox regression analysis.

Functional enrichment analysis

Biological features of DHA targets against ESCA were probed out utilizing clusterProfiler package (Yu et al., 2012). Gene Ontology (GO) covering biological process, cellular component and molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were separately implemented.

Gene set enrichment analysis (GSEA)

The potentially altered biological processes between high and low expression subpopulations of DHA targets against ESCA were determined with GSEA tool (Subramanian et al., 2005). The “c5. go.v7.5.1. symbols” and “c2. cp.kegg.v7.5.1. symbols” gene sets acquired from the Molecular Signatures Database Molecular Signatures Database (Liberzon et al., 2015) were utilized as reference sets.

Molecular docking

The 3D structures of DHA targets against ESCA were accessed from the PDB database and saved in pdb format. PyMOL software was used for protein pretreatment. In addition, the 3D structure of DHA was created using Chem3D software and saved in mol*2 format. The pre-processed targets and DHA were imported into AutoDockTools 1.5.6 software for molecular docking, and the results were saved in pdbqt format. Vina scripts were run for computing molecular binding energies as well as visualizing molecular docking results. Additionally, Discovery Studio 2019 software was implemented for searching for docking sites as well as calculating LibDockScore for flexible binding. The molecular docking results were imported into PyMOL software for visualizing the molecular docking conformation.

Molecular dynamics simulation

Molecular dynamics simulation was implemented by use of Discovery Studio 2019 software. The molecular structure of DHA was added to the CHARMM force field by “simulation” module, and DHA targets against ESCA were solvated by “solvation” module. The molecular dynamics simulation parameters were set by “standard dynamics cascade module for the targets added with the solvent system. The system temperature was elevated from 50 K to 300 K under 100 ns of analog sampling; the time step was set to 1 ns; and other parameters were set to default values. The Root Mean Square Deviation (RMSD) as well as Root Mean Square Float (RMSF) values were analyzed.

Results

Network pharmacology determines dihydroartemisinin targets against esophageal carcinoma

From the PubChem chemical information resource, we accessed the 2D chemical structure of DHA, as illustrated in Figure 1A. Through integrating the SwissTargetPrediction, HERB, as well as BATMAN-TCM platforms, potential therapeutic targets of DHA were predicted. In accordance with the SwissTargetPrediction platform, the target class distribution for the first 25 DHA targets was visualized (Figure 1B). Enzyme was the predominant target class, occupying 28.0%. In total, 111, 3, and 51 DHA targets were separately predicted by the SwissTargetPrediction, HERB, as well as BATMAN-TCM platforms (Figure 1C). Under merging and deduplication, 160 DHA targets were finally acquired (Supplementary Table S1).

FIGURE 1

FIGURE 1. Network pharmacology determines DHA targets against ESCA. (A) The 2D chemical structure of DHA from the PubChem chemical information resource. (B) The target class distribution for the first 25 DHA targets predicted by the SwissTargetPrediction platform. (C) The distribution of DHA targets predicted by the SwissTargetPrediction, HERB, as well as BATMAN-TCM platforms. (D) Volcano diagram of the DEGs between 161 ESCA and 11 normal tissue specimens in TCGA-ESCA dataset. Red, up-regulation; green, down-regulation. (E) Heatmap of the expression values of the first 100 most significantly DEGs between 161 ESCA and 11 normal tissue specimens in TCGA-ESCA dataset. (F) Venn diagram of the shared genes between DHA targets and DEGs. (G,H) The PPI network of DHA targets against ESCA through the (G) STRING online tool or (H) cytoHubba plugin of Cytoscape tool.

In TCGA-ESCA dataset, 821 DEGs (comprising 429 genes with up-regulation as well as 392 genes with down-regulation) between 161 ESCA and 11 normal tissue specimens were screened in line with the screening condition of adjusted p-value < 0.01 together with |log2 fold-change|>1 (Figure 1D). The specific information of DEGs was listed in Supplementary Table S2. The first 100 most significantly DEGs were depicted in Figure 1E. Venn diagram illustrated 16 DHA targets against ESCA shared by DHA targets and DEGs (Figure 1F), comprising MMP9, CA2, FBP1, CCNA1, CCNA2, CDK2, AOC3, MMP1, TOP2A, AURKA, CDK1, CCNB1, CCNB2, KCNQ1, ADORA2B, and GABRE. Their closely interactions were investigated in the PPI network (Figures 1G,H and Table 1).

TABLE 1

TABLE 1. DHA targets against ESC ranked by degree.

Establishment and verification of a reliable prognostic signature for esophageal carcinoma

The prognostic implication of DHA targets against ESCA was investigated in TCGA-ESCA dataset. Univariate cox regression analysis demonstrated that AURKA, and ADORA2B significantly correlated to ESCA cases’ OS. AURKA, and ADORA2B were utilized for establishing a multivariate cox regression model. The risk score was computed on the basis of the coefficients and expression values of AURKA, and ADORA2B. TCGA-ESCA cases were equally divided into discovery set and validation set. In discovery set, we assigned cases into high- and low-risk subpopulations. Survival difference between subpopulations was assessed. High-risk subpopulation exhibited worse OS outcome in comparison to low-risk subpopulation (Figure 2A). Survival status was also compared between subpopulations. More dead cases were investigated in high-risk subpopulation (Figure 2B). For evaluating whether the prognostic signature enabled to estimate ESCA cases’ OS, ROC curve was plotted. In Figure 2C, the AUC value was 0.637, indicative of the excellent performance of the prognostic signature in estimating ESCA cases’ OS. The prognostic signature was further verified in validation set and entire set. As expected, it exhibited the remarkable advantage in ESCA cases’ prognosis prediction (Figures 2D–I), indicating the clinical generalizability of this signature.

FIGURE 2

FIGURE 2. Establishment and verification of a reliable prognostic signature for ESCA cases in TCGA-ESCA dataset. (A) Kaplan-Meier OS curves and log-rank test between high- and low-risk subpopulations in discovery set. (B) Distribution of alive and dead status in high- and low-risk subpopulations in discovery set. (C) ROC curve on the basis of the prognostic signature in discovery set. (D) Kaplan-Meier OS curves and log-rank test between high- and low-risk subpopulations in validation set. (E) Distribution of alive and dead status in high- and low-risk subpopulations in validation set. (F) ROC curve of the prognostic signature in validation set. (G) Kaplan-Meier OS curves and log-rank test between high- and low-risk subpopulations in entire set. (H) Distribution of alive and dead status in high- and low-risk subpopulations in entire set. (I) ROC curve based on the prognostic signature in entire set. (J) Forest diagram of univariate cox regression analysis for the associations of clinicopathological factors and the prognostic signature with ESCA cases’ OS. (K) Forest diagram of multivariate cox regression analysis for the associations of clinicopathological factors and the prognostic signature with ESCA cases’ OS.

The prognostic signature independently predicts esophageal carcinoma patients’ prognosis

Utilizing univariate cox regression analysis, we investigated the associations of clinicopathological factors and the prognostic signature with ESCA cases’ OS. As illustrated in Figure 2J, the prognostic signature along with pathological stage, N stage and M stage were linked to worse ESCA cases’ OS. Multivariate cox regression analysis was conducted for estimating which factors independently predicted ESCA cases’ OS. As a result, the prognostic signature acted as an independent risk factor (Figure 2K).

Dihydroartemisinin targets against esophageal carcinoma correlate to cell cycle progression

Biological features of DHA targets against ESCA were probed out utilizing GO and KEGG enrichment approach. For biological process, histone phosphorylation, mitotic DNA integrity checkpoint, G2/M transition of mitotic cell cycle, cell cycle G2/M phase transition, mitotic nuclear envelope disassembly, DNA damage response, signal transduction by p53 class mediator resulting in cell cycle arrest, cell cycle G1/S phase transition, signal transduction involved in mitotic G1 DNA damage checkpoint, intracellular signal transduction involved in G1 DNA damage checkpoint and signaling transduction involved in mitotic cell cycle checkpoint were significantly linked to DHA targets against ESCA (Figures 3A,B). Cellular components of cyclin-dependent protein kinase holoenzyme complex, serine/threonine protein kinase complex, protein kinase complex, transferase complex, transferring phosphorus-containing groups, condensed chromosome, pronucleus, chromosomal region, condensed nuclear chromosome, centromeric region, spindle microtubule, and microvillus were remarkably enriched by DHA targets against ESCA (Figures 3C,D). Also, they correlated to molecular functions of histone kinase activity, cyclin-dependent serine/threonine kinase regulator activity, protein kinase regulator activity, kinase regulator activity, cyclin-dependent protein serine/threonine kinase activity, cyclin-dependent protein kinase activity, cyclin binding, protein heterodimerization activity, metalloendopeptidase activity as well as protein serine/threonine kinase activity (Figures 3E,F). As illustrated in KEGG pathway enrichment results, progesterone-mediated oocyte maturation, cell cycle, cellular senescence, oocyte meiosis, p53 signaling pathway, hepatitis B, viral carcinogenesis, human T-cell leukemia virus one infection, AMPK signaling pathway as well as FoxO signaling pathway (Figures 3G,H). Altogether, DHA targets against ESCA might participate in mediating cell cycle progression.

FIGURE 3

FIGURE 3. Biological features of DHA targets against ESCA. (A) The first ten biological process terms enriched by DHA targets against ESCA. (B) Interactions between the first ten biological process terms. (C) The first ten cellular component terms enriched by DHA targets against ESCA. (D) Interactions between the first ten cellular component terms. (E) The first ten molecular function terms enriched by DHA targets against ESCA. (F) Interactions between the first ten molecular function terms. (G) The first ten KEGG pathways enriched by DHA targets against ESCA. (H) Interactions between the first ten KEGG pathways.

Dihydroartemisinin targets against esophageal carcinoma (AURKA and ADORA2B) correlates to prognosis and pathological stage of esophageal carcinoma

Among DHA targets against ESCA, DORA2B, and AURKA significantly correlated to ESCA prognosis. Up-regulated ADORA2B was linked to better OS, while up-regulated AURKA was correlated to worse OS in TCGA-ESCA dataset (Figures 4A,B). Further analysis for the associations of ADORA2B and AURKA with pathological stage of ESCA was carried out. In Figure 4C, stage II cases exhibited higher ADORA2B expression. Additionally, AURKA expression displayed positive associations to pathological stage (Figure 4D), indicating that AURKA was linked to ESCA progression.

FIGURE 4

FIGURE 4. DHA targets against ESCA (ADORA2B and AURKA) correlate to prognosis and pathological stage of ESCA as well as tumorigenic pathways. (A) Kaplan-Meier OS curves between high and low ADORA2B expression subpopulations in TCGA-ESCA dataset. (B) Kaplan-Meier OS curves between high and low AURKA expression subpopulations in TCGA-ESCA dataset. (C,D) ADORA2B and AURKA expression across distinct pathological stages. (E,F) GO enrichment results of ADORA2B and AURKA. (G,H) KEGG pathway enrichment results of ADORA2B and AURKA.

Mechanisms underlying ADORA2B and AURKA

Mechanisms underlying ADORA2B and AURKA were further probed out utilizing GSEA. For GO enrichment results, negative regulation of metallopeptidase activity, sensory perception of smell, intermediate filament, intermediate filament cytoskeleton and olfactory receptor activity exhibited positive interactions to ADORA2B (Figure 4E). Meanwhile, alternative mRNA splicing via spliceosome, negative regulation of wound healing, sensory perception of smell, odorant binding and olfactory receptor activity displayed positive correlations to AURKA (Figure 4F). For KEGG pathway enrichment results, folate biosynthesis, maturity onset diabetes of the young, olfactory transduction, sulfur metabolism as well as taurine and hypotaurine metabolism were positively linked to ADORA2B (Figure 4G). Also, cytosolic DNA sensing pathway, graft versus host disease, maturity onset diabetes of the young, olfactory transduction, and regulation of autophagy exhibited positive correlations to AURKA (Figure 4H).

Molecular docking between dihydroartemisinin and potential targets

Utilizing Chem3D, we obtained the 3D structure of DHA with mol*2 format, and downloaded the 3D structure of potential targets from the PDB database with pdb format. Through AutoDockTools 1.5.6 software, the 3D structure of DHA and potential targets were converted to pdbqt format, and thus searched for the active pocket. Afterwards, we run the Vina script to compute the binding energies of DHA and potential targets. As listed in Table 2, the binding energies of the docking bodies formed by DHA targets against ESCA (MMP9, CA2, FBP1, CCNA1, CCNA2, CDK2, AOC3, MMP1, TOP2A, AURKA, CDK1, CCNB1, CCNB2, KCNQ1, ADORA2B, and GABRE) and DHA were all <−5.0 kcal mol⁻¹, indicative of the stable docking. Additionally, we computed LibDockScore to dock DHA and its potential targets by use of Discovery Studio 2019 software. As a result, in each potential target, docking site of DHA can be found, and the LibDockScore of the docking model formed by DHA and each target was >50. The 3D and 2D molecular docking models between DHA and potential targets (MMP9, CA2, FBP1, CCNA1, CCNA2, CDK2, AOC3, MMP1, TOP2A, AURKA, CDK1, CCNB1, CCNB2, KCNQ1, ADORA2B, and GABRE) were separately displayed as Figures 5A–P and Figures 6A–P.

TABLE 2

TABLE 2. Molecular docking between DHA and potential targets.

FIGURE 5

FIGURE 5. The 3D molecular docking models between DHA and potential targets. (A) MMP9; (B) CA2; (C) FBP1; (D) CCNA1; (E) CCNA2; (F) CDK2; (G) AOC3; (H) MMP1; (I) TOP2A; (J) AURKA; (K) CDK1; (L) CCNB1; (M) CCNB2; (N) KCNQ1; (O) ADORA2B; and (P) GABRE.

FIGURE 6

FIGURE 6. The 2D molecular docking models between DHA and potential targets. (A) MMP9; (B) CA2; (C) FBP1; (D) CCNA1; (E) CCNA2; (F) CDK2; (G) AOC3; (H) MMP1; (I) TOP2A; (J) AURKA; (K) CDK1; (L) CCNB1; (M) CCNB2; (N) KCNQ1; (O) ADORA2B; and (P) GABRE.

Molecular dynamics simulation

The conformations of ADORA2B-DHA and AURKADHA molecular docking were used for subsequent molecular dynamics simulation analysis. The ADORA2B-DHA complex was added by 16,738 water molecules, 44 sodium as well as 53 chloride, while the AURKA-DHA complex was added by 5,688 water molecules, 15 sodium as well as 20 chloride. We further investigated the structural stability of ADORA2B-DHA and AURKADHA complexes during molecular dynamics simulation. Meanwhile, the RMSD values of the two complexes were computed during 100 ns molecular dynamics simulation. The ADORA2B-DHA and AURKADHA complexes were stable following 100 ns molecular dynamics simulation, as illustrated in Figure 7A. The RMSD value of ADORA2B-DHA complex mainly fluctuated from 1.65612 to 2.15719, and the mean RMSD value was 1.89383; while the RMSD value of AURKA-DHA complex mainly fluctuated from 1.05842 to 1.44505, and its mean RMSD value was 1.26881 (Figure 7A). The RMSD fluctuation values of the two complexes were all within a reasonable range, demonstrating that the ADORA2B-DHA and AURKA-DHA complexes were in a stable state during the process of molecular dynamics simulation. For analyzing the volatility of distinct amino acids in the two complexes during the process of molecular dynamics simulation, we also computed the RMSF values of all amino acids during the simulation. ADORA2B-DHA complex fluctuated greatly around amino acids Ser6, Ala7, Phe8, Leu315, and Phe316 (Figure 7B); AURKA-DHA complex fluctuated greatly around amino acids Gln127, Trp128, Ala129, Phe329, Glu330, Ala331, and Asn332 (Figure 7C). The heatmaps of hydrogen bonds during the ADORA2B-DHA and AURKADHA molecular docking were separately exhibited in Figures 7D,E. Hydrogen bonding interactions were present in most conformations of ADORA2B-DHA and AURKADHA molecular docking, demonstrating that these hydrogen bonds were persistent and stable. The molecular dynamics simulation results of ADORA2B-DHA and AURKA-DHA complexes were separately displayed in Figures 8A,B.

FIGURE 7

FIGURE 7. Molecular dynamics simulation of ADORA2B-DHA and AURKA-DHA complexes. (A) Alterations in RMSD values during the process of molecular dynamics simulation of ADORA2B-DHA (red) and AURKA-DHA (green) complexes. (B,C) Alterations in RMSF values during the process of molecular dynamics simulation of (B) ADORA2B-DHA and (C) AURKA-DHA complexes. (D,E) Heatmaps of hydrogen bonding interactions of (D) ADORA2B-DHA and (E) AURKA-DHA complexes.

FIGURE 8

FIGURE 8. Molecular dynamics simulation of ADORA2B-DHA and AURKA-DHA complexes. (A) ADORA2B-DHA complex; (B) AURKA-DHA complex.

Discussion

At present, the paradigms of development of single-targeted therapeutic agents are challenging, primarily because of lack of effectiveness as well as resistance. Therefore, natural compounds selectively acting on two or more targets exhibit higher efficacy in comparison to single-target drugs. The current study employed network pharmacology as well as molecular docking approaches to predict multi-targets and multi-pathways of DHA against ESCA.

Previous studies have proposed that DHA exhibits the anti-ESCA property through targeting distinct molecules as well as pathways (Li et al., 2019; Ma et al., 2020). Thus, we presented a systematic analysis for unveiling DHA targets against ESCA. Through integrating the SwissTargetPrediction, HERB and BATMAN-TCM platforms, 160 DHA targets were finally acquired. Meanwhile, 821 DEGs between 161 ESCA and 11 normal tissue specimens in TCGA-ESCA dataset were screened, which might participate in ESCA initiation as well as progression. Following the intersection of DHA targets and DEGs, we eventually determined 16 DHA targets against ESCA, comprising MMP9, CA2, FBP1, CCNA1, CCNA2, CDK2, AOC3, MMP1, TOP2A, AURKA, CDK1, CCNB1, CCNB2, KCNQ1, ADORA2B, and GABRE. These DHA targets against ESCA remarkably correlated to cell cycle progression. Evidence demonstrates that DHA triggers cell cycle arrest during tumor progression (Lin et al., 2016; Saidi et al., 2021).

Among them, ADORA2B as well as AURKA correlated to ESCA cases’ OS. AURKA is a membrane of serine/threonine kinase family, and its activation is essential for cell division process by modulating mitosis (Du et al., 2021). AURKA accelerates ESCC progression via improving the activity of distinct pathways such as EGFR-PI3K-Akt (Du et al., 2020; Shi et al., 2021). Nevertheless, whether ADORA2B participates in ESCA remains uncharted. On the basis of ADORA2B as well as AURKA, a prognostic signature was developed for survival prediction. Despite the reliability and independency of the prognostic signature in ESCA cases’ OS, large prospective ESCA cohorts are required for verifying it.

Molecular docking uses mapping software to place small molecule compounds on the binding region of macromolecular targets, and then calculates parameters to predict the binding ability and binding mode of the two (Tao et al., 2020). Through the strength of the binding ability, the possible mechanisms of action of the small molecule compound can be preliminarily inferred, and the interaction mode between the small molecule compound and the potential target can be quickly and accurately described, which provides a scientific basis for the preparation of derivatives, and thus shortens the drug development cycle and reduces research and development costs. In this study, through use of molecular docking, the stable binding of DHA to 16 potential targets were verified, and thus speculated the main anti-ESCA targets and interactions. Molecular dynamics simulation can simulate the combination of the two in the natural state, and can restore some reaction processes that cannot be detected by current technical means. We employed molecular dynamics simulations for confirming the stability of the binding of DHA to protein targets (ADORA2B and AURKA). Despite this, experiments are essential for confirming ADORA2B-DHA and AURKA-DHA interactions in ESCA.

Conclusion

Collectively, the current research unveiled that DHA exerted an anti-ESCA effect through targeting multiple targets (especially ADORA2B and AURKA) as well as multiple pathways (especially cell cycle progression) on the basis of network pharmacology and molecular docking approaches, and thus provided the theoretical basis for the pharmacological research of DHA against ESCA.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

HW conceived and designed the study and conducted all the experiments and data analysis, and wrote the manuscript. All authors reviewed and approved the manuscript.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.1017520/full#supplementary-material

Supplementary Table S1 | DHA targets predicted by the SwissTargetPrediction, HERB, as well as BATMAN-TCM platforms.

Supplementary Table S2 | The specific information of 821 DEGs between 161 ESCA and 11 normal tissue specimens in TCGA-ESCA dataset.

Abbreviations

ESCA, esophageal carcinoma; ESCC, squamous cell carcinoma; EAC, adenocarcinoma; DHA, dihydroartemisinin; TCGA, The Cancer Genome Atlas; PPI, Protein-protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes; OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene set enrichment analysis; RMSD, Root Mean Square Deviation; RMSF, Root Mean Square Float.

References

Bartha, Á., and Győrffy, B. (2021). TNMplot.com: A web tool for the comparison of gene expression in normal, tumor and metastatic tissues. Ijms 22 (5), 2622. doi:10.3390/ijms22052622

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Niu, X., Qiao, X., Liu, S., Ma, H., Shi, X., et al. (2021). Characterization of interplay between autophagy and ferroptosis and their synergistical roles on manipulating immunological tumor microenvironment in squamous cell carcinomas. Front. Immunol. 12, 739039. doi:10.3389/fimmu.2021.739039

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., He, L. Y., Lai, S., and He, Y. (2020). Dihydroartemisinin inhibits the migration of esophageal cancer cells by inducing autophagy. Oncol. Lett. 20 (4), 1. doi:10.3892/ol.2020.11955

PubMed Abstract | CrossRef Full Text | Google Scholar

Chin, C. H., Chen, S. H., Wu, H. H., Ho, C. W., Ko, M. T., and Lin, C. Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol. 8, S11. doi:10.1186/1752-0509-8-s4-s11

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, W., Fang, T., Duan, Z., Xiang, D., Wang, Y., Zhang, M., et al. (2020). Dihydroartemisinin sensitizes esophageal squamous cell carcinoma to cisplatin by inhibiting sonic Hedgehog signaling. Front. Cell Dev. Biol. 8, 596788. doi:10.3389/fcell.2020.596788

PubMed Abstract | CrossRef Full Text | Google Scholar

Daina, A., Michielin, O., and Zoete, V. (2019). SwissTargetPrediction: Updated data and new features for efficient prediction of protein targets of small molecules. Nucleic Acids Res. 47 (W1), W357–w364. doi:10.1093/nar/gkz382

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, R., Huang, C., Chen, H., Liu, K., Xiang, P., Yao, N., et al. (2020). SDCBP/MDA-9/syntenin phosphorylation by AURKA promotes esophageal squamous cell carcinoma progression through the EGFR-PI3K-Akt signaling pathway. Oncogene 39 (31), 5405–5419. doi:10.1038/s41388-020-1369-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, R., Huang, C., Liu, K., Li, X., and Dong, Z. (2021). Targeting AURKA in cancer: Molecular mechanisms and opportunities for cancer therapy. Mol. Cancer 20 (1), 15. doi:10.1186/s12943-020-01305-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, S., Dong, L., Liu, L., Guo, J., Zhao, L., Zhang, J., et al. (2021). Herb: A high-throughput experiment- and reference-guided database of traditional Chinese medicine. Nucleic Acids Res. 49 (D1), D1197–d1206. doi:10.1093/nar/gkaa1063

PubMed Abstract | CrossRef Full Text | Google Scholar

Hulshof, M., Geijsen, E. D., Rozema, T., Oppedijk, V., Buijsen, J., Neelis, K. J., et al. (2021). Randomized study on dose escalation in definitive chemoradiation for patients with locally advanced esophageal cancer (ARTDECO study). Jco 39 (25), 2816–2824. doi:10.1200/jco.20.03697

PubMed Abstract | CrossRef Full Text | Google Scholar

Janjigian, Y. Y., Shitara, K., Moehler, M., Garrido, M., Salman, P., Shen, L., et al. (2021). First-line nivolumab plus chemotherapy versus chemotherapy alone for advanced gastric, gastro-oesophageal junction, and oesophageal adenocarcinoma (CheckMate 649): A randomised, open-label, phase 3 trial. Lancet 398 (10294), 27–40. doi:10.1016/s0140-6736(21)00797-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, M., Wu, Y., Qi, L., Li, L., Song, D., Gan, J., et al. (2021). Dihydroartemisinin mediating PKM2-caspase-8/3-GSDME axis for pyroptosis in esophageal squamous cell carcinoma. Chemico-Biological Interact. 350, 109704. doi:10.1016/j.cbi.2021.109704

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., et al. (2021). PubChem in 2021: New data content and improved web interfaces. Nucleic Acids Res. 49 (D1), D1388–d1395. doi:10.1093/nar/gkaa971

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Ma, Q., Xu, L., Gao, C., Yao, L., Wen, J., et al. (2021). Human telomerase reverse transcriptase as a therapeutic target of dihydroartemisinin for esophageal squamous cancer. Front. Pharmacol. 12, 769787. doi:10.3389/fphar.2021.769787

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Huang, P., Gan, J., Ling, X., Du, X., Liao, Y., et al. (2019). Dihydroartemisinin represses esophageal cancer glycolysis by down-regulating pyruvate kinase M2. Eur. J. Pharmacol. 854, 232–239. doi:10.1016/j.ejphar.2019.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The molecular signatures database hallmark gene set collection. Cell Syst. 1 (6), 417–425. doi:10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, R., Zhang, Z., Chen, L., Zhou, Y., Zou, P., Feng, C., et al. (2016). Dihydroartemisinin (DHA) induces ferroptosis and causes cell cycle arrest in head and neck carcinoma cells. Cancer Lett. 381 (1), 165–175. doi:10.1016/j.canlet.2016.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Zhang, M., Ying, S., Zhang, C., Lin, R., Zheng, J., et al. (2017). Genetic Alterations in esophageal tissues from squamous dysplasia to carcinoma. Gastroenterology 153 (1), 166–177. doi:10.1053/j.gastro.2017.03.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Guo, F., Wang, Y., Li, C., Zhang, X., Li, H., et al. (2016). BATMAN-TCM: A bioinformatics analysis tool for molecular mechANism of traditional Chinese medicine. Sci. Rep. 6, 21146. doi:10.1038/srep21146

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Z., Wang, J., Shu, Y., Liu, L., Kong, L., Yang, L., et al. (2022). Sintilimab versus placebo in combination with chemotherapy as first line treatment for locally advanced or metastatic oesophageal squamous cell carcinoma (ORIENT-15): Multicentre, randomised, double blind, phase 3 trial. Bmj 377, e068714. doi:10.1136/bmj-2021-068714

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, H., Lu, J., Bai, Y., Mao, T., Wang, J., Fan, Q., et al. (2021). Effect of camrelizumab vs placebo added to chemotherapy on survival and progression-free survival in patients with advanced or metastatic esophageal squamous cell carcinoma. Jama 326 (10), 916–925. doi:10.1001/jama.2021.12836

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Q., Liao, H., Xu, L., Li, Q., Zou, J., Sun, R., et al. (2020). Autophagy-dependent cell cycle arrest in esophageal cancer cells exposed to dihydroartemisinin. Chin. Med. 15, 37. doi:10.1186/s13020-020-00318-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukherjee, S., Hurt, C., Radhakrishna, G., Gwynne, S., Bateman, A., Gollins, S., et al. (2021). Oxaliplatin/capecitabine or carboplatin/paclitaxel-based preoperative chemoradiation for resectable oesophageal adenocarcinoma (NeoSCOPE): Long-term results of a randomised controlled trial. Eur. J. Cancer 153, 153–161. doi:10.1016/j.ejca.2021.05.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, X., Chen, L., Li, Y., Hu, Z., and He, F. (2022). Ferroptosis, necroptosis, and pyroptosis in the tumor microenvironment: Perspectives for immunotherapy of SCLC. Seminars Cancer Biol. doi:10.1016/j.semcancer.2022.03.009

CrossRef Full Text | Google Scholar

Niu, X., Huang, B., Qiao, X., Liu, J., Chen, L., and Zhong, M. (2021). MicroRNA-1-3p suppresses malignant phenotypes of ameloblastoma through down-regulating lysosomal associated membrane protein 2-mediated autophagy. Front. Med. 8, 670188. doi:10.3389/fmed.2021.670188

PubMed Abstract | CrossRef Full Text | Google Scholar

Safran, H. P., Winter, K., Ilson, D. H., Wigle, D., DiPetrillo, T., Haddock, M. G., et al. (2022). Trastuzumab with trimodality treatment for oesophageal adenocarcinoma with HER2 overexpression (NRG oncology/RTOG 1010): A multicentre, randomised, phase 3 trial. Lancet Oncol. 23 (2), 259–269. doi:10.1016/s1470-2045(21)00718-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Saidi, H., Murtaza, B., Khan, A. S., Koceir, E. A., Hichami, A., and Khan, N. A. (2021). DHA induces Jurkat T-cell arrest in G2/M phase of cell cycle and modulates the plasma membrane expression of TRPC3/6 channels. Biochimie 181, 169–175. doi:10.1016/j.biochi.2020.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, K., Zhang, J. Z., Yang, L., Li, N. N., Yue, Y., Du, X. H., et al. (2021). Protein deubiquitylase USP3 stabilizes Aurora A to promote proliferation and metastasis of esophageal squamous cell carcinoma. BMC Cancer 21 (1), 1196. doi:10.1186/s12885-021-08934-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shitara, K., Özgüroğlu, M., Bang, Y. J., Di Bartolomeo, M., Mandalà, M., Ryu, M. H., et al. (2021). Molecular determinants of clinical outcomes with pembrolizumab versus paclitaxel in a randomized, open-label, phase III trial in patients with gastroesophageal adenocarcinoma. Ann. Oncol. 32 (9), 1127–1136. doi:10.1016/j.annonc.2021.05.803

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, J. M., Shen, L., Shah, M. A., Enzinger, P., Adenis, A., Doi, T., et al. (2021). Pembrolizumab plus chemotherapy versus chemotherapy alone for first-line treatment of advanced oesophageal cancer (KEYNOTE-590): A randomised, placebo-controlled, phase 3 study. Lancet 398 (10302), 759–771. doi:10.1016/s0140-6736(21)01234-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA A Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi:10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, Q., Du, J., Li, X., Zeng, J., Tan, B., Xu, J., et al. (2020). Network pharmacology and molecular docking analysis on molecular targets and mechanisms of Huashi Baidu formula in the treatment of COVID-19. Drug Dev. Industrial Pharm. 46 (8), 1345–1353. doi:10.1080/03639045.2020.1788070

PubMed Abstract | CrossRef Full Text | Google Scholar

UniProt (2021). The universal protein knowledgebase in 2021. Nucleic Acids Res. 49 (D1), D480–d489. doi:10.1093/nar/gkaa1100()

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, R. H., Zhang, Y., Pan, H., Feng, J., Zhang, T., Liu, T., et al. (2021). Efficacy and safety of weekly paclitaxel with or without ramucirumab as second-line therapy for the treatment of advanced gastric or gastroesophageal junction adenocarcinoma (RAINBOW-Asia): A randomised, multicentre, double-blind, phase 3 trial. Lancet Gastroenterology Hepatology 6 (12), 1015–1024. doi:10.1016/s2468-1253(21)00313-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Li, B., Yi, J., Hua, R., Chen, H., Tan, L., et al. (2022). Robot-assisted versus conventional minimally invasive esophagectomy for resectable esophageal squamous cell carcinoma. Ann. Surg. 275 (4), 646–653. doi:10.1097/sla.0000000000005023

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Liang, H., Li, Z., Xue, Y., Wang, Y., Zhou, Z., et al. (2021). Perioperative or postoperative adjuvant oxaliplatin with S-1 versus adjuvant oxaliplatin with capecitabine in patients with locally advanced gastric or gastro-oesophageal junction adenocarcinoma undergoing D2 gastrectomy (RESOLVE): An open-label, superiority and non-inferiority, phase 3 randomised controlled trial. Lancet Oncol. 22 (8), 1081–1092. doi:10.1016/s1470-2045(21)00297-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Wang, Y., and Meng, L. (2022). Comparative genomic analysis of esophageal squamous cell carcinoma and adenocarcinoma: New opportunities towards molecularly targeted therapy. Acta Pharm. Sin. B 12 (3), 1054–1067. doi:10.1016/j.apsb.2021.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, L., Chen, X., Zhu, Y., Qin, J., Niu, T., Ding, Y., et al. (2020). Dihydroartemisinin inhibits the proliferation of esophageal squamous cell carcinoma partially by targeting AKT1 and p70S6K. Front. Pharmacol. 11, 587470. doi:10.3389/fphar.2020.587470

PubMed Abstract | CrossRef Full Text | Google Scholar