A network pharmacology-based investigation of emodin against pancreatic adenocarcinoma

Abstract

Pancreatic adenocarcinoma (PAAD) is one of the most common malignancies worldwide with an increasing incidence and poor outcome due to the lack of effective diagnostic and treatment methods. Emerging evidence implicates that emodin displays extensive spectrum anticancer properties. Differential expression genes in PAAD patients were analyzed by Gene Expression Profiling Interactive Analysis (GEPIA) website, and the targets of emodin were obtained via Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform. Subsequently, enrichment analyses were performed using R software. A protein–protein interaction (PPI) network was constructed by STRING database and Cytoscape software was used to identify the hub genes. Prognostic value and immune infiltration landscapes were explored through Kaplan–Meier plotter (KM plotter) website and the Single-Sample Gene Set Enrichment Analysis package of R. Finally, molecular docking was used to computationally verify the interaction of ligand and receptor proteins. A total of 9191 genes were significantly differentially expressed in PAAD patients and 34 potential targets of emodin were obtained. Intersections of the 2 groups were considered as potential targets of emodin against PAAD. Functional enrichment analyses illustrated that these potential targets were linked to numerous pathological processes. Hub genes identified through PPI networks were correlated with poor prognosis and infiltration level of different immune cells in PAAD patients. Perhaps emodin interacted with the key molecules and regulate the activity of them. We revealed the inherent mechanism of emodin against PAAD with the aid of network pharmacology, which provided reliable evidence and a novel guideline for clinical treatment.

1. Introduction

Pancreatic cancer, with 5-year survival <10%, is one of the clinically common malignancies. For the last few years, surgical resection obviously improved the survival of pancreatic cancer patients. However, the annual death toll of this disease is a continuous rise due to the increasing clinical morbidity. Pancreatic adenocarcinoma (PAAD), which is thought to arise from proliferative premalignant pancreatic intraepithelial neoplasia of the ductal epithelium, accounts for more than 90% of all diagnosed pancreatic cancer.^[1] According to the pathological features of PAAD, only 15% to 20% of patients undergo an operation.^[2] It is difficult to detect early on account of nonspecific symptoms of PAAD patients. Early metastasis is a characteristic of PAAD and over 50% of patients have distant metastasis when they seek medical treatment, which seriously affects the prognosis of patients.^[3] With the development of tumor immune microenvironment research in recent years, infiltration level of multiple immune cells is considered closely associated with diagnosis and prognosis of PAAD. It is urgent to find novel diagnosis and treatment targets because many antineoplastic therapies exhibit limited effectiveness.

Traditional Chinese medicine (TCM) is an inseparable part of Chinese traditional culture, with more than 2000 years of inheritance and accumulation. The value of TCM has been discovered and recognized worldwide. A growing body of research suggests that TCM inhibits tumor growth and increases the sensitivity of tumor cells to chemotherapy through different mechanisms.^[4] Emodin, a natural anthraquinone derivative occurs in numerous widely used Chinese herbal medicine, has diverse pharmacological properties like diuretic, vasorelaxant, anti-ulcerogenic, anti-oxidant, anti-inflammatory, antimicrobial, anti-viral, antidiabetic, and hepatoprotective activities.^[5,6] Previous studies have demonstrated that emodin probably suppressed tumor growth and reduced the metastasis of various cancer, such as hepatocellular carcinoma, cervical, breast, colon, and non-small-cell lung cancer.^[7–10] It has previously been observed that the mechanisms can influence an epithelial-mesenchymal transition, inhibit angiogenesis, and induce necroptosis.^[8,10,11] Besides, emodin improved an immunosuppression of tumor environment via regulating the tumor-promoting feedforward interactions between cancer cells and macrophages.^[12] However, the evidence for emodin-targeted molecules is inconclusive.

In this study, differentially expressed genes (DEGs) of PAAD issues were identified and targets of emodin were obtained via Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP). The intersection 2 groups were considered as potential targets of emodin against PAAD. Functional enrichment analyses were performed by the clusterProfiler package of R software. The protein–protein interaction (PPI) network was constructed through the STRING database. Plugins of Cytoscape software were used to recognize the key molecules of emodin against PAAD. Subsequently, the correlation between key molecules and prognostic value and immune infiltration level were analyzed via the Single-Sample Gene Set Enrichment Analysis package in R. Molecular docking was performed to verify an interaction of emodin and key molecules, suggesting that emodin played an important role in regulating the activity of key molecules. We disclosed the therapeutic value and potential mechanism of emodin in PAAD patients through network pharmacology, hoping to provide more evidence and guidance for PAAD treatment.

2. Materials and methods

2.1. Differential expression gene analysis

Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) is a widely used interactive website for analyzing the RNA sequencing expression data of 9736 tumors and 8587 normal samples from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and the Genotype-Tissue Expression projects 28407145.^[13] We used GEPIA to compare gene expression in PAAD tissues and normal samples from TCGA database. On account of adjust P value of DEGs were arbitrarily close to zero, we transformed it to −1*Log₁₀(adjust P value) and constructed a volcano map by GraphPad Prism (8.0.1). Adjusted P value <.05 and |log2 fold change (FC)| > 1 were selected as thresholds to screen out DEGs between tumor and normal groups. Lowly expressed genes were shown as green triangles and red dots represented highly expressed genes.

2.2. Acquisition of drug targets

TCMSP (https://old.tcmsp-e.com/tcmsp.php) captures the relationships between drugs, targets, and diseases.^[14] We obtained targets of emodin through TCMSP and compared them with high-expressed genes in PAAD tissues. Genes of the intersection were considered as potential targets of emodin in the treatment of PAAD. The results were visualized via a Venn diagram which was constructed by the ggplot2 package in R software.

2.3. Functional enrichment analyses

Gene Ontology (GO) is a comprehensive source of digital data relating to the functions of genes in 3 independent categories: molecular function, biological process, and cellular component. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database that stores information about genomes, biological pathways, diseases, and drugs. The clusterProfiler package in R was used to perform GO and KEGG analysis of potential targets. A false discovery rate <0.05 was considered as a statistically significant standard. Enrichment results were visualized using the R ggplot2 package.

2.4. PPI network and hub gene analyses

The STRING database (https://string-db.org/), which covers 24’584’628 proteins from 5’090 organisms, is used to predict and construct PPI networks of input genes.^[15] We constructed a PPI network of potential targets via the STRING database, and the minimum required interaction score was 0.4. Subsequently, the PPI network was reconstructed by Cytoscape software (3.8.0).^[16] PPI network contains a few nodes and edges which represent their interactions. Among the nodes that has the most interaction is considered as a hub gene. CytoHubba (http://apps.cytoscape.org/apps/cytohubba), a plugin of Cytoscape software (3.8.0), was used to point out hub genes for predicted PPI network of potential targets.^[17] Molecular Complex Detection (http://apps.cytoscape.org/apps/mcode) plugin of Cytoscape software is used to detect the most profound modules from the PPI network.^[18]

2.5. Prognosis and immune infiltration analyses

Kaplan–Meier plotter (KM plotter, http://kmplot.com/analysis/index.php?p=service) is an online survival analysis website which can carry out research with 54,675 genes and 18,674 cancer samples, involving breast cancer, lung cancer, etc. In order to further explore the impact of potential target expression on PAAD patients, we used KM plotter website to analyze prognostic value of key molecules which were included in the most profound module. Patients were divided into high- and low-expression groups according to the median expression level of mRNA. Log rank test was used to evaluate significant differences of survival time between the 2 groups. Subsequently, the RNAseq FPKM data of PAAD patients was downloaded from TCGA database and the marker genes of 24 immune cell types were obtained from a paper of immunity.^[19] Immune filtration landscapes of key molecules were quantified by Single-Sample Gene Set Enrichment Analysis method, which was carried out using GSVA package from R software (http://www.bioconductor.org).^[20] The correlation between immune cells and key molecules was analyzed by Spearman rank‐correlation coefficient.^[21]

2.6. Molecular docking

Molecular docking is a technique used to mimic the interaction between small ligand molecules and receptor protein macromolecules. The affinity can be represented by binding energy which lower than 0 indicates combining spontaneously and stable conformations. The three-dimensional structures of receptor protein were downloaded from the RCSB PDB database (https://www.rcsb.org/). The three-dimensional structure of molecule ligand was obtained from TCMSP platform as MOL2 format files. Pymol software (version 2.4.0) was used to preprocess the receptor proteins, including removing ligands and water. AutodockTools software (version 1.5.6) was performed to process receptor proteins and small molecule ligand as follows: add hydrogens, compute Gasteiger, assign AD4 type, set torsion tree, and adjust the grid box to include all protein structures. Subsequently, Autogrid4 and Autodock4 programs were operated to get the docking results. Converting File Formats were fulfilled by OpenBabel software (version 2.4.1). The optimum conformations were visualized in Pymol software.

3. Results

3.1. Potential targets of emodin in PAAD

We compared gene expression in PAAD tissues with normal tissues based on GEPIA database. A total of 9191 genes were differentially expressed in PAAD, among which 8716 genes were upregulated and 475 genes were downregulated (Fig. 1A). Simultaneously, 34 targets of emodin were retrieved from TCMSP website. Highly-expressed genes and emodin targets were compared by R programming language and there were 26 common molecules: ACTA2, BTK, CALM1, CASP3, CSF2, F10, FLT1, FLT4, HSP90AA1, IGHG1, IL1B, KDR, MAOB, MMP1, MMP9, NCOA2, PIK3CG, PPARG, PRKACA, PRKCD, PTGS1, PTGS2, SLC2A1, TGFB1, TOP2A, and TP53 (Fig. 1B). These 26 molecules were considered as potential targets of emodin against PAAD for further study.

Figure 1.: Potential targets of emodin against pancreatic adenocarcinoma (PAAD) and functional enrichment analyses. (A) Differentially expression genes (DEGs) of PAAD; (B) intersection of DEGs in PAAD patients and emodin targets; (C) gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of potential targets.

3.2. Functional enrichment analysis of potential targets

To further explore the function of these 26 potential targets and identify key candidate pathways, GO and KEGG pathway enrichment analyses were performed. The top 20 enrichment terms were listed in Table 1. The major enriched terms were positive regulation of nucleocytoplasmic transport, positive regulation of reactive oxygen species metabolic process, membrane raft, membrane microdomain, protein tyrosine kinase activity, IL-17 signaling pathway, and MAPK signaling pathway (Fig. 1C). Supplemental Table 1, https://links.lww.com/MD/I791 depicted the entire enrichment items. Results of enrichment analyses indicated that these potential targets were associated with pathological progress such as cell metabolism, intracellular and extracellular signaling transduction, and immune response.

Table 1 - Functional enrichment analyses of potential targets.

Ontology	ID	Description	GeneRatio	BgRatio	P value	p.adjust	qvalue
BP	GO:0046824	Positive regulation of nucleocytoplasmic transport	6/26	62/18670	2.29e-10	5.13e-07	2.05e-07
BP	GO:2000379	Positive regulation of reactive oxygen species metabolic process	6/26	102/18670	4.83e-09	4.06e-06	1.62e-06
BP	GO:0046822	Regulation of nucleocytoplasmic transport	6/26	104/18670	5.43e-09	4.06e-06	1.62e-06
BP	GO:0018108	Peptidyl-tyrosine phosphorylation	8/26	363/18670	2.18e-08	1.04e-05	4.16e-06
BP	GO:0018212	Peptidyl-tyrosine modification	8/26	366/18670	2.32e-08	1.04e-05	4.16e-06
CC	GO:0045121	Membrane raft	6/26	315/19717	2.79e-06	1.49e-04	1.11e-04
CC	GO:0098857	Membrane microdomain	6/26	316/19717	2.84e-06	1.49e-04	1.11e-04
CC	GO:0098589	Membrane region	6/26	328/19717	3.52e-06	1.49e-04	1.11e-04
CC	GO:0042629	Mast cell granule	2/26	22/19717	3.80e-04	0.010	0.007
CC	GO:0044853	Plasma membrane raft	3/26	109/19717	3.89e-04	0.010	0.007
MF	GO:0004713	Protein tyrosine kinase activity	5/26	134/17697	1.34e-06	2.41e-04	1.10e-04
MF	GO:0097718	Disordered domain specific binding	3/26	33/17697	1.49e-05	0.001	6.12e-04
MF	GO:0043621	Protein self-association	3/26	56/17697	7.41e-05	0.004	0.002
MF	GO:0004714	Transmembrane receptor protein tyrosine kinase activity	3/26	62/17697	1.01e-04	0.005	0.002
MF	GO:0019199	Transmembrane receptor protein kinase activity	3/26	79/17697	2.07e-04	0.007	0.003
KEGG	hsa04657	IL-17 signaling pathway	7/25	94/8076	9.36e-09	1.73e-06	8.86e-07
KEGG	hsa04010	MAPK signaling pathway	8/25	294/8076	1.77e-06	1.46e-04	7.48e-05
KEGG	hsa04915	Estrogen signaling pathway	6/25	138/8076	3.03e-06	1.46e-04	7.48e-05
KEGG	hsa05418	Fluid shear stress and atherosclerosis	6/25	139/8076	3.16e-06	1.46e-04	7.48e-05
KEGG	hsa05323	Rheumatoid arthritis	5/25	93/8076	8.05e-06	2.98e-04	1.52e-04

3.3. PPI between the potential targets

To illustrate the functional association of potential targets, we constructed a PPI network within potential targets through STRING website. There were 24 nodes and 109 edges involved in the network (Fig. 2A). The top 10 genes identified by Cytohubba were TP53, IL1B, CASP3, KDR, MMP9, PTGS2, TGFB1, CSF2, FLT1, and MMP1. Hub genes were shown in different colors from high to low according to the interaction score (Fig. 2B). In addition, we extracted a highly dense module (11 nodes, 53 edges) using the Molecular Complex Detection plug-in for subsequent study (Fig. 2C).

Figure 2.: Protein-protein interaction (PPI) networks. (A) PPI network of potential targets; (B) hub genes identified by Cytohubba; (C) a highly dense module extracted by molecular complex detection.

3.4. Prognostic value of hub genes

To reveal the prognostic value of hub genes which were included in the highly dense module, we analyzed the relationship between gene expression and overall survival (OS) and relapse free survival in PAAD patients via KM plotter website. Results showed that high expression of CASP3 (logrank P = .017), CSF2 (logrank P = .0042), IL1B (logrank P = .037), MMP1 (logrank P = .0011), MMP9 (logrank P = .046), PPARG (logrank P = .014), and PTGS2 (logrank P = .034) were significantly correlated with poor OS. Low expression of KDR (logrank P = .03) was associated with poor OS in PAAD patients (Fig. 3). Simultaneously, high expression of CASP3 (logrank P = .042), CSF2 (logrank P = .031), FLT1 (logrank P = .0077), MMP1 (logrank P = .0091), MMP9 (logrank P = .034), PPARG (logrank P = .013), PTGS2 (logrank P = .018), and TGFB1 (logrank P = .022) were significantly related to poor relapse free survival (Fig. 4).

Figure 3.: The relationship of hub gene expression and overall survival.

Figure 4.: The relationship of hub gene expression and relapse free survival.

3.5. Immune infiltration level analysis of key molecules

We ulteriorly explored whether these key molecules with prognostic value played more important roles in PAAD. The correlation was analyzed between expression of key molecules and infiltration of 24 immune cell types. Results showed that these key molecules were significantly correlated with the infiltration levels of various immune cell subpopulations, among which the infiltration level of Th1 and Th2 cells had a widely positive correlation with gene expression (Fig. 5, Supplemental Table 2–11, https://links.lww.com/MD/I792; https://links.lww.com/MD/I793; https://links.lww.com/MD/I794; https://links.lww.com/MD/I795; https://links.lww.com/MD/I796; https://links.lww.com/MD/I797; https://links.lww.com/MD/I798; https://links.lww.com/MD/I799; https://links.lww.com/MD/I800; https://links.lww.com/MD/I801). The infiltration level of Th17 and NK CD56bright cells was negatively associated with the expression of key molecules.

Figure 5.: the correlation of immune infiltration and hub gene expression.

3.6. Results of molecular docking

Molecular docking was used to reveal whether emodin played an important role in regulation of these 10 key molecules. Protein structures of candidate molecules were found from PDB database as follows: 3KJF,^[22] 2GMF,^[23] 3HNG, 3LTQ,^[24] 2OH4,^[25] 2CLT,^[26] 1L6J,^[27] 6C5T,^[28] 5F19,^[29] and 6P7J.^[30] Results of molecular docking showed that emodin (MOL000472) had a strong connection with each candidate receptor protein. The docking energies were listed in Table 2 and the docking details were demonstrated in Figure 6.

Table 2 - Energy of molecular docking.

Ligand	Target name	PDB ID	Binding energy
Emodin	CASP3	3KJF	−4.41
	CSF2	2GMF	−5.62
	FLT1	3HNG	−5.18
	IL1B	3LTQ	−5.29
	KDR	2OH4	−5.93
	MMP1	2CLT	−4.84
	MMP9	1L6J	−5.87
	PPARG	6C5T	−5.92
	PTGS2	5F19	−4.03
	TGFB1	6P7J	−4.66

Figure 6.: molecular docking of emodin and key molecules.

4. Discussion

PAAD is a type of highly invasive tumor with significant morbidity, mortality, and poor prognosis.^[31] The curative effect of common treatment methods, such as surgery, chemotherapy, and targeted therapies, seems limited due to recurrence and drug resistance.^[32] The medical value and mechanism of TCM have been studied extensively. The antitumor effect of emodin, which is a natural anthraquinone derivative widely found in various Chinese medicinal herbs, has been proved by a growing body of research.^[5] In our study, the potential targets of emodin against PAAD were identified by network pharmacology. Functional enrichment analyses were performed to reveal the biology progress of potential targets. We explored the correlation between key molecules and prognostic value and immune infiltration landscape. Molecular docking was used to exhibit the interaction of emodin and key molecules, suggesting that emodin played an important role in regulating their activity.

Prognosis has always been an important aspect in the area of cancer. In this research, high-expression of 9 targets, namely CASP3, CSF2, FLT1, IL1B, MMP1, MMP9, PPARG, PTGS2, and TGFB1, were related to poor prognosis. Low-expression of KDR suggested poor outcome. Caspase family has been well known for their apoptotic roles. It was also found that CASP3 promoted tumor proliferation, metastasis, and angiogenesis.^[33,34] CSF2 acted as an oncogene in colon cancer and its high-expression indicating poor prognosis.^[35] Activation of FLT1 promoted the metastasis of breast cancer.^[36] IL1B drove the bone metastasis of breast cancer.^[37] MMP1 and MMP9 were generally upregulated in all cancers and they played essential roles in angiogenesis in turn to accelerate cancer cell proliferation and migration.^[38] PTGS2 played a nonnegligible role in tumorigenesis and resistance of cancer cells to radiotherapy and chemotherapy.^[39] PPARG promoted cancer growth in bladder, liver, and prostate cancers.^[40] The above conclusions were consistent with results of our research. However, as a key regulator of lipid metabolism in adipocytes, PPARG had opposite function in cancer tissues: it inhibited tumor proliferation in breast and lung cancers.^[40] High-expression of KDR was correlated with poor prognosis in lung cancer patients.^[41] These findings were different from results in our study, suggesting that pathologic progresses of different tumors were disparate.

Immune filtration of various cells in tumor microenvironment has different effects on tumor development. In our study, infiltration levels of 4 immune cells, involving Th1, Th2, Th17, and NK CD56 bright cell, were significantly associated with expression of key molecules. Th1 cells produced IFN-γ to increase antigen presentation and promote anticancer immunity.^[42] Although IFN-γ also promoted tumor growth by inducing PD-L1 expression, the infiltration level of Th1 was positively correlated with prognosis in CRC patients.^[43] IL-4 and IL-13 produced by Th2 cells activated the mutated KRAS and JAK1/-STAT6/-MYC pathway and further regulated glycolysis to promote tumorigenesis.^[44] The role of Th17 cells playing in tumors was complex and depended on the type of cancer.^[45] NK CD56 bright cells produced high level of immunoregulatory cytokines and had poorly cytotoxic so they were considered as silencers in antitumor response and devoted in immunomodulation.^[46,47] These results indicated that key molecules of emodin against PAAD maybe affect the prognosis of patients via regulating infiltration level of different immune cells and thereby influencing immune response.

Functional enrichment analyses exhibited that potential targets were associated with oxygen metabolism, tyrosine kinase related biological progresses, IL-17, and MAPK signaling pathway. Cancer cells mainly produce energy through aerobic glycolysis, also named the Warburg effect, due to the histopathological characteristics of abnormal vasculature and a large quantity of energy for proliferation.^[48,49] Hypoxia is also a momentous cause of high aggressiveness and treatment-resistant.^[50] Previous studies revealed that disorder of receptor tyrosine kinases was associated with progression and metastasis of tumors.^[51] Tyrosine phosphorylation was also essential for metabolic reprogramming of cancer cells and tumor progression.^[52] As a major proinflammatory cytokine, IL-17 also plays an important role in metabolism-modulating including metabolic regulation within immune cells and effects of immune response on systemic metabolism. Activity of IL-17 is critical for the development of many diseases. Research showed that IL-17 targeted epithelial cells and mesenchymal cells and regulated glycolysis to promote tumorigenesis.^[53] Perhaps IL-17 also induced reprogramming of embryonic stem cells to stimulate development and invasion of PAAD.^[54] MAPK signaling pathway was vital for cell proliferation and located downstream of many growth-factor receptors.^[55] Potential targets significantly correlated with these biological progresses indicated that emodin possibly regulated metabolism and immune response in PAAD patients to impact growth and invasion of tumor.

In this study, we explored the anticancer effect and mechanism of emodin against PAAD via network pharmacology. Data involved in this research was originated from TCGA database and TCMSP platform, so that all results were based on bioinformatics and reasonable in theory and need consist on experimental verification. We will intensively study the authenticity and reliability of underlying mechanisms revealed in this research. In vitro and in vivo experiments are supposed to be conducted to ascertain the influence of emodin on the development of PAAD.

5. Conclusion

Above all, we investigated the potential targets of emodin against PAAD using network pharmacology. A total of 26 potential targets were identified, which suggested that emodin probably inhibited tumor development by affecting tumor metabolism and immune infiltration. Besides, the expression of key molecules was also closely associated with poor prognosis in PAAD patients, which emphasized the nonnegligible function of emodin targets.

Acknowledgments

The authors sincerely acknowledge the Helixlife for some bioinformatics approaches to explore databases.

Author contributions

Conceptualization: Xueying Shi, Akio Mizushima, Shilin Xia.

Data curation: Jingyi Zhu.

Formal analysis: Jingyi Zhu.

Investigation: Bingqian Huang.

Methodology: Bingqian Huang, Shilin Xia.

Project administration: Akio Mizushima, Shilin Xia.

Writing – original draft: Xueying Shi.

Writing – review & editing: Xueying Shi, Bingqian Huang, Jingyi Zhu, Takuji Yamaguchi, Ailing Hu, Masahiro Tabuchi, Daisuke Watanabe, Seiichiro Yoshikawa, Shinobu Mizushima, Akio Mizushima, Shilin Xia.