Systematic analysis of potential targets of the curcumin analog pentagamavunon-1 (PGV-1) in overcoming resistance of glioblastoma cells to bevacizumab.

BACKGROUND: Glioblastoma is one of the most aggressive and deadliest malignant tumors. Acquired resistance decreases the effectiveness of bevacizumab in glioblastoma treatment and thus increases the mortality rate in patients with glioblastoma. In this study, the potential targets of pentagamavunone-1 (PGV-1), a curcumin analog, were explored as a complementary treatment to bevacizumab in glioblastoma therapy.
METHODS: Target prediction, data collection, and analysis were conducted using the similarity ensemble approach (SEA), SwissTargetPrediction, STRING DB, and Gene Expression Omnibus (GEO) datasets. Gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted using Webgestalt and DAVID, respectively. Hub genes were selected based on the highest degree scores using the CytoHubba. Analysis of genetic alterations and gene expression as well as Kaplan-Meier survival analysis of selected genes were conducted with cBioportal and GEPIA. Immune infiltration correlations between selected genes and immune cells were analyzed with database TIMER 2.0.
RESULTS: We found 374 targets of PGV-1, 1139 differentially expressed genes (DEGs) from bevacizumab-resistant-glioblastoma cells. A Venn diagram analysis using these two sets of data resulted in 21 genes that were identified as potential targets of PGV-1 against bevacizumab resistance (PBR). PBR regulated the metabolism of xenobiotics by cytochrome P450. Seven potential therapeutic PBR, namely GSTM1, AKR1C3, AKR1C4, PTGS2, ADAM10, AKR1B1, and HSD17B110 were found to have genetic alterations in 1.2%-30% of patients with glioblastoma. Analysis using the GEPIA database showed that the mRNA expression of ADAM10, AKR1B1, and HSD17B10 was significantly upregulated in glioblastoma patients. Kaplan-Meier survival analysis showed that only patients with low mRNA expression of AKR1B1 had significantly better overall survival than the patients in the high mRNA group. We also found a correlation between PBR and immune cells and thus revealed the potential of PGV-1 as an immunotherapeutic agent via targeting of PBR.
CONCLUSION: This study highlighted seven PBR, namely, GSTM1, AKR1C3, AKR1C4, PTGS2, ADAM10, AKR1B1, and HSD17B110. This study also emphasized the potential of PBR as a target for immunotherapy with PGV-1. Further validation of the results of this study is required for the development of PGV-1 as an adjunct to immunotherapy for glioblastoma to counteract bevacizumab resistance.

Introduction

Glioblastoma, a malignant tumor which can grow in the brain and spinal cord, is one of the deadliest types of cancer in the world (Friedmann-Morvinski, 2014). The incidence and mortality rates of glioblastoma have been increasing from year to year (Korja et al., 2019). Bevacizumab is a monoclonal antibody against vascular endothelial growth factor (VEGF), which is an angiogenic factor secreted by endothelial and glioblastoma cells (Tamura et al., 2017). Therapy using bevacizumab is an option for preventing metastasis (Huang et al., 2017). However, the development of resistance by glioblastoma cells causes a decrease in the effectiveness of bevacizumab and increases mortality in glioblastoma patients (Ramezani et al., 2019). Therefore, researchers continue to try to find adjuvant for glioblastoma therapy.

Curcumin is a natural compound that has been widely studied for its anticancer properties in various types of cancer, including breast (Wang et al., 2016), colon (Selvam et al., 2019), and brain cancer (Sordillo et al., 2015). Curcumin has also been proven to increase the anticancer activity of chemotherapeutic agents in breast cancer (Wen et al., 2019), colon cancer (Su et al., 2018), and glioblastoma (Dhandapani et al., 2007). However, the use of curcumin is limited due to its low solubility and bioavailability (Peng et al., 2018). Pentagamavunone-1 (PGV-1), shown in Fig. 1A, is a synthetic analog of curcumin that was developed to overcome the problem of low bioavailability (Dai et al., 2012). PGV-1 exhibited more potent cytotoxicty than curcumin on various cancer cells, especially breast cancer cells (Da’i et al., 2007, Dai et al., 2011, Hermawan et al., 2011, Meiyanto et al., 2014) and colon cancer cells (Meiyanto et al., 2018). This compound increases the toxicity of doxorubicin in breast cancer cells (Meiyanto et al., 2014), inhibits breast cancer metastasis (Meiyanto et al., 2019), and shows better anticancer activity than curcumin on K562 cancer cells (Lestari et al., 2019). More importantly, PGV-1 was known to inhibit the expression of the angiogenic factors, VEGF and COX-2 (Meiyanto et al., 2006). Accordingly, PVG-1 is potential to be developed as a complementary treatment to bevacizumab in glioblastoma therapy.

Fig. 1

(A) Chemical structure of PGV-1. (B) Flowchart of the study. (C). A Venn diagram of predicted PGV-1 target genes and mRNA from GSE81645 (bevacizumab-resistant U87 cells). (D) GO enrichment analysis of potential target genes of PGV-1 in overcoming glioblastoma resistance to bevacizumab.

In this study, we performed a systematic analysis using a bioinformatics approach to identify targets and molecular mechanisms of PGV-1 as a complementary treatment to overcome glioblastoma resistance to bevacizumab (Fig. 1B). PGV-1 targets were predicted using several public databases, while differentially expressed genes (DEGs) in glioblastomas that are resistant to bevacizumab were downloaded from the GEO. Several analyses, including Gene Ontology; KEGG pathway enrichment; hub gene selection; analysis of genetic alterations, gene expression, and Kaplan–Meier survival; and immune infiltration correlations were further conducted to explore the importance of target genes which can serve as the basis for the development of PGV-1 as an adjunct to bevacizumab in combination therapy for glioblastoma.

Material and methods

Target prediction and data collection

Direct targets (DT) of PGV-1 were predicted using the similarity ensemble approach (SEA) () (Keiser et al., 2007) and SwissTargetPrediction (), (Daina et al., 2019) using the default settings of the databases. Indirect targets of PGV-1 were collected from each DT using STRING DB with a confidence level of 0.4 and a maximum of 20 interactions. A total of direct and indirect targets of PGV-1 were collected and further named as potential targets of PGV-1. Microarray data on bevacizumab-resistant xenograft mouse U87 glioma cells were downloaded from GEO datasets, namely, GSE45161, entitled Acquired resistance to anti-VEGF therapy in glioblastoma is associated with a mesenchymal transition (Piao et al., 2013). A Venn diagram was generated from potential targets of PGV-1 and DEGs of bevacizumab-resistant xenograft mice and resulted in potential targets of PGV-1 against bevacizumab resistance (PBR).

Functional annotation

Gene ontology analysis of the PBR was conducted using ORA from Webgestalt (WEB-based GEne SeT AnaLysis Toolkit (), using default settings (Wang et al., 2017). KEGG pathway enrichment analysis was performed using DAVID () version 6.8 (Huang et al., 2009). Levels of P < 0.05 was selected as the threshold for significance.

Protein-protein interaction network and selection of hub genes

Protein-protein interactions among the PBR was constructed using STRING DB () version 11.0 and visualized with Cytoscape version 3.7.1 (Shannon et al., 2003). Hub genes, which are genes with a high correlation within the network, were selected on the basis of the highest degree score using the CytoHubba plug-in of Cytoscape (Chin et al., 2014).

Analysis of genetic alterations within the glioblastoma study

An analysis of genetic alterations in selected genes was conducted with cBioportal (). Briefly, a gene list, encoded as a set of gene symbols, was submitted as a query in cBioportal and searched among glioblastoma studies. Selected studies with the highest alteration number were selected and analyzed for Oncoprint analysis to observe the alterations among genes. Additional mutual exclusivity was performed among genes using a cutoff value of p < 0.05 as the selection criterion.

Gene expression and Kaplan–Meier survival analysis

Analysis of gene expression and Kaplan–Meier survival curve analysis was conducted using GEPIA () (Tang et al., 2017). Briefly, gene symbols were submitted to GEPIA, glioblastoma tissues and adjacent tumors were analyzed for gene expression, and a Kaplan–Meier survival analysis was performed. Levels of P < 0.05 was selected as the threshold for significance.

Infiltration of immune cells

Immune cell infiltration was analyzed using TIMER 2.0 (). The TIMER database was used to explore the correlation between tumor infiltrates and gene expression (Li et al., 2020). TIMER is a comprehensive resource for the systematic analysis of immune cell infiltration in various cancer cells using samples from the TCGA study. In this study, the correlations between the expression of PBR (GSTM1, AKR1C3, PTGS2, AKR1B1, ADAM10, and HSD17B10) and immune cell infiltration, including B cells, CD8, CD4, neutrophils, dendritic cells, macrophages, and cancer-related fibroblasts (CAF), were determined and compared.

Results

PGV-1 target prediction using the similarity ensemble approach (SEA) and SwissTargetPrediction resulted in 6 (Table 1) and 12 (Table 2) direct targets (DT) of PGV-1, respectively. Indirect targets of PGV-1 were collected from each DT using STRING DB (Supplementary Table 1). A total 374 targets were collected from direct and indirect targets of PGV-1 (Supplementary Table 2). Analysis of microarray data of bevacizumab-resistant xenograft mice U87 glioma cells from GSE45161 using GEO2R resulted in 1139 differentially expressed genes (DEGs), consisting of 350 and 789 upregulated and downregulated genes, respectively (Supplementary Table 3). Further analysis of potential overlapping targets of PGV-1 and DEGs of bevacizumab-resistant xenograft mice using a Venn diagram resulted in 21 genes (Fig. 1C, Supplementary Table 4) that were further named as PBR.

Table 1

Predicted targets of PGV-1, as analyzed by SEA.

Query	Target Key	Target Name	Description	P-Value
TTHY_HUMAN	TTR	Transthyretin	6.66E-16	0.44
LOX5_HUMAN	ALOX5	Arachidonate 5-lipoxygenase	1.06E-11	0.44
ARY1_HUMAN	NAT1	Arylamine N-acetyltransferase 1	1.14E-08	0.33
EF2K_HUMAN	EEF2K	Eukaryotic elongation factor 2 kinase	7.97E-08	0.55
CISD1_HUMAN	CISD1	CDGSH iron-sulfur domain-containing protein 1	4.16E-07	0.3
THRB_HUMAN	F2	Prothrombin	2.71E-06	0.46

Table 2

Predicted targets of PGV-1, as analyzed by SwissTargetPrediction.

Target	Common name	Uniprot ID	Target Class	Probability*	Known actives (3D/2D)
Serine/threonine-protein kinase EEF2K	EEF2K	O00418	Kinase	0.109339753	0/2    Â
Beta-secretase 1	BACE1	P56817	Protease	0.109339753	15/95    Â
Coagulation factor VII/tissue factor	F3	P13726	Surface antigen	0.109339753	2/10    Â
11-beta-hydroxysteroid dehydrogenase 1	HSD11B1	P28845	Enzyme	0.109339753	8/11    Â
Tyrosine-protein kinase SRC	SRC	P12931	Kinase	0.109339753	5/3    Â
Histone acetyltransferase p300	EP300	Q09472	Writer	0.109339753	0/3    Â
Estradiol 17-beta-dehydrogenase 2	HSD17B2	P37059	Enzyme	0.109339753	29/0    Â
Serum albumin	ALB	P02768	Secreted protein	0.109339753	2/0    Â
Glyoxalase I	GLO1	Q04760	Enzyme	0.109339753	0/5    Â
Inhibitor of NF-kappa-B kinase (IKK)	IKBKG IKBKB CHUK	Q9Y6K9 O14920 O15111	Kinase	0.109339753	0/3    Â
Tyrosine-protein kinase JAK3	JAK3	P52333	Kinase	0.109339753	17/0    Â
Tyrosine-protein kinase JAK1	JAK1	P23458	Kinase	0.109339753	17/0    Â

Gene ontology analysis of the PBR was conducted on the basis of three categories consisting of biological process, cellular component, and molecular function (Fig. 1D). PBR were found to regulate the biological processes involved in metabolic pathways, response to stimulus, and cell communication. In terms of cellular localization, PBR were found in the membrane, nucleus, and cytosol. Lastly, PBR played a role in the molecular functions of protein and ion binding and of transferase activity. KEGG pathway enrichment analysis using DAVID showed that PBR regulated xenobiotic metabolism by cytochrome P450 (Table 3).

Table 3

KEGG pathway enrichment analysis of the potential target genes of PGV-1 against bevacizumab resistance in glioblastoma.

Term	P Value	Genes
hsa05010:Alzheimer's disease	0.003974712	HSD17B10, APP, NDUFA2, ADAM10, APOE
hsa00980:Metabolism of xenobiotics by cytochrome P450	0.024899131	GSTM1, AKR1C3, AKR1C4
hsa04610:Complement and coagulation cascades	0.032276953	PLAT, C3, F3

The construction of a protein–protein interaction network using STRING revealed a network that consisted of a node and an edge, with a clustering coefficient (Fig. 2A). Selection of hub genes based on degree score resulted in 10 genes, including APP, APOE, PTGS2, ADAM10, C3, F3, HSD17B10, IRS1, AKR1B1, and PLAT (Fig. 2B, Table 4).

Fig. 2

(A) Protein-protein interaction network of potential target genes of PGV-1 in overcoming glioblastoma resistance to bevacizumab, analyzed by STRING. (B) Top 10 hub genes based on highest degree score, analyzed by CytoHubba.

Table 4

Top 10 hub genes based on highest degree score, analyzed by CytoHubba.

No	Symbol	Degree score
1	APP	8
2	APOE	7
3	PTGS2	6
4	ADAM10	5
5	C3	5
6	F3	4
7	HSD17B10	4
8	IRS1	4
9	AKR1B1	4
10	PLAT	3

Analysis of genetic alterations among glioblastoma studies

Seven genes were selected for analysis of genetic alterations using cBioportal. PTGS2, ADAM10, HSD17B10, and AKR1B1 were selected on the basis of the highest degree scores which represented hub genes in the protein–protein interaction network analysis. GSTM1, AKR1C3, and AKR1C4 were selected from the KEGG pathway enrichment analysis, which revealed xenobiotic metabolism by cytochrome P450 as a pathway regulated by the PBR. Across five glioblastoma studies in the cBioportal database, the GBM study (Mayo PDX) reported the highest number of genetic alterations among others, and therefore was selected for further analysis (Fig. 3A). Genetic alterations in seven genes were found in 1.2%–30% of patients with glioblastoma, including AKR1C3 (1.2%), AKR1C4 (1.2%), PTGS2 (1.2%), ADAM10 (1.2%), HSD17B10 (2.4%), AKR1B1 (6%), and GSTM1 (30%) (Fig. 3B). Additional mutual exclusivity analysis results revealed co-occurrence in six pairs of genes which were AKR1C3-AKR1C4, AKR1C3-PTGS2, AKR1C4-PTGS2, AKR1C3-ADAM10, AKR1C4-ADAM10, and PTGS2-ADAM10 (Table 5).

Fig. 3

(A) Overview of genetic changes in GSTM1, AKR1B1, AKR1C3, AKR1C4, PTGS2, ADAM10, and HSD17B10 across five glioblastoma studies, as analyzed by cBioportal. (B) Summary of alterations in GSTM1, AKR1B1, AKR1C3, AKR1C4, PTGS2, ADAM10, and HSD17B10 across glioblastoma patients using Mayo PDX study.

Table 5

Mutual exclusivity of selected genes among glioblastoma study.

A	B	Log2 Odds Ratio	p-Value	Tendency
AKR1C3	AKR1C4	>3	0.012	Co-occurrence
AKR1C3	PTGS2	>3	0.012	Co-occurrence
AKR1C4	PTGS2	>3	0.012	Co-occurrence
AKR1C3	ADAM10	>3	0.012	Co-occurrence
AKR1C4	ADAM10	>3	0.012	Co-occurrence
PTGS2	ADAM10	>3	0.012	Co-occurrence

Overall, no genetics alterations were found in the protein domain of the seven genes that is responsible for the PGV-1 binding. The genetic alterations in GSTM1 were amplification or overexpression (Supplementary Fig. 2A). In AKR1C3 and AKR1C4, deep deletions were found which led to decreased mRNA expression (Supplementary Fig. 2B-C). PTGS2 had a splice site mutation, namely X57_splice, which is located outside of the peroxidase domain (Supplementary Fig. 2D). A missense mutation in ADAM10, namely V443I, was located in the end of repro lysine domain (Supplementary Fig. 2E). Several missesnse mutations were reported in AKR1B1, namely, G39R, V154M, and T141M which are located in the aldo ketoreductase domain, and E315G which is located outside the aldo ketoreductase domain, and further analysis indicated a few additional copies of AKR1B1 due to those mutations (Supplementary Fig. 2F). Two missense mutations were found in HSD17B10, namely G80R and A97T which is located in the short chain dehydrogenase domain (Supplementary Fig. 2G).

Analysis of PTGS2, ADAM10, HSD17B10, AKR1B1, GSTM1, AKR1C3, and AKR1C4 expression in samples from patients with glioblastoma using the GEPIA database showed that the mRNA expression levels of ADAM10, AKR1B1, and HSD17B10 were significantly upregulated in glioblastoma patients compared to those in normal patients (Fig. 4), whereas no differences in transcript levels were found for GSTM1, AKR1C3, AKR1C4, and PTGS2. Kaplan–Meier survival analysis showed that only patients with low mRNA expression of AKR1B1 had significantly better overall survival than those with high mRNA expression (Fig. 5). In addition, the mRNA expression of GSTM1, AKR1C3, PTGS2, ADAM10, AKR1B1, and HSD17B10 did not affect survival.

Fig. 4

mRNA levels of GSTM1, AKR1B1, AKR1C3, AKR1C4, PTGS2, ADAM10, and HSD17B10 in patients with glioblastoma, as analyzed by GEPIA.

Fig. 5

Overall survival of patients with glioblastoma related to the mRNA levels of GSTM1, AKR1B1, AKR1C3, PTGS2, and ADAM10, as analyzed by GEPIA.

Correlation between PBR and infiltration of immune cells into the glioblastoma

TIMER 2.0 was used to explore the correlation between target gene expression and the level of immune cell infiltration in GBM. PTGS2 expression (p = 7.42 × 10−5) had a negative correlation (r =  − 0.331), whereas HSD17B10 expression (p = 4.07 × 10−5) had a positive correlation (r = 0.342) with the purity of GBM (Table 6, Supplementary Fig. 1). Positive correlations were seen in B cells with both ADAM10 (r = 0.171, p = 0.0462) and AKR1B1 (r = 0.316, p = 1.67 × 10−4), and in CD8+ cells with AKR1B1 (r = 0.193, p = 0.0237). Neutrophils showed a positive correlation with both PTGS2 (r = 0.57, p = 3.52 × 10−13)and ADAM10 (r = 0.2, p = 0.0189). Dendritic cells showed a negative correlation with AKR1C3 (r =  − 0.172, p = 0.045) and positive correlations with both PTGS2 (r = 0.363, p = 1.36 × 10−5) and ADAM10 (r = 0.496, p = 6.85 × 10−10). Macrophages were positively correlated with ADAM10 (r = 0.205, p = 9.11 × 10−5) and AKR1B1 (r = 0.381, p = 4.38 × 10−6). Cancer-associated fibroblasts had a negative correlation with AKR1C3 (r =  − 0.171, p = 0.0453), and a positive correlation with PTGS2 (r = 0.447, p = 4.45 × 10−8) and ADAM10 (r = 0.328, p = 9.11 × 10−5).

Table 6

Correlation analysis of PGV-1 target gene expression with the levels of immune cells.

Gene	Purity		B Cells		CD8		CD4		Neutrophil		Dendritic cells		Macrophage		Cancer Associated Fibroblast
	Cor	P	Cor	P	Cor	P	Cor	P	Cor	P	Cor	P	Cor	P	Cor	P
GSTM1	0.09	0.294	0.018	0.836	0.246	0.0372	0.041	0.637	−0.033	0.699	0.074	0.387	0.186	0.0292	0.154	0.0727
AKR1C3	−0.2	0.184	0.124	0.149	0.176	0.0397	−0.035	0.682	−0.103	0.230	−0.172	0.045	0.052	0.549	−0.171	0.0453
PTGS2	−0.331	7.42e-05	−0.137	0.109	0.021	0.812	−0.032	0.711	0.57	3.52e-13	0.363	1.36e-05	0.105	0.220	0.447	4.45e-08
ADAM10	0.027	0.757	0.171	0.0462	−0.147	0.0873	0.127	0.140	0.2	0.0189	0.496	6.85e-10	0.205	9.11e-05	0.328	9.11e-05
AKR1B1	−0.126	0.140	0.316	1.67e-04	0.193	0.0237	−0.046	0.594	−0.014	0.871	0.055	0.522	0.381	4.38e-06	−0.118	0.16
HSD17B10	0.342	4.07e-05	0.12	0.163	0.045	0.600	0.036	0.674	−0.135	0.116	−0.115	0.181	0.07	0.419	−0.11	0.2

Discussion

This study aimed to explore the targets and molecular mechanisms of PGV-1 in overcoming glioblastoma resistance to bevacizumab. This study produced seven PBR, namely GSTM1, AKR1C3, AKR1C4, PTGS2, ADAM10, AKR1B1, and HSD17B110. Genetic alterations in these genes were found in 1.2%–30% of patients with glioblastoma, with the highest frequency of genetic alterations in GSTM1 at 30% and AKR1B1 at 6%. We did not find any alterations in the protein domain responsible for the PGV-1 binding in these seven PBR. Interestingly, we did not find any annotation related to oncogenic, diagnostic, prognostic and therapeutic levels for these genetic alterations, and therefore the significance of these findings can be further explored in future studies.

Validation of mRNA expression in a sample of patients with GEPIA showed that the expression of ADAM10, AKR1B1, and HSD17B10 was significantly upregulated in glioblastoma patients. This contradicts previous research revealing downregulation of the same three genes. Survival analysis using a Kaplan–Meier plot showed that patients with low AKR1B1 mRNA expression showed significantly better overall survival than patients in the group whose mRNA expression was higher. Additionally, mutual exclusivity analysis showed six pairs of genes with co-occurring mutations. This result indicates the importance of the genes such as AKR1C3, AKR1C4, PTGS2, and ADAM10 in the biological mechanism of PGV-1 in overcoming glioblastoma resistance to bevacizumab.

Glioblastoma resistance to bevacizumab is regulated by several mechanisms, including hypoxia (Mao et al., 2016) and activation of the autophagy pathway via PI3K/Akt signaling (Huang et al., 2018), JAK/STST, NFkB, JNK, and Ras (Itatani et al., 2018). Other mechanisms of glioblastoma resistance to bevacizumab have also been proposed, including upregulation of angiopoietin-2 in endothelial cellswhich leads to reduction of tumor-associated macrophages (Scholz et al., 2016) and activation of the kinase signaling pathway (Ramezani et al., 2019). In addition, resistance mechanism to bevacizumab therapy also involved the increased invasive and metastac properties and recruitment of myeloid cells and stromal cells (Haibe et al., 2020). A previous study showed that in silico analysis revealed the inhibitory effects of curcumin and PGV-1 against receptor tyrosine kinase, including HER2 and EGFR and their downstream targets (Meiyanto et al., 2014). In addition, PGV-1 decreased the expression of the angiogenic factors, VEGF and COX-2, in T47D breast cancer cells (Meiyanto et al., 2006). PGV-1 also enhances sensitivity of colon cancer cells by inhibition of NF-κB activity (Meiyanto et al., 2018). Taken together, PGV-1 could potentially target the mechanisms of glioblastoma resistance to bevacizumab.

Aldo keto reductases (AKRs) are a family of enzymes that are involved in oxidoreduction reactions of endogenous substrates or xenobiotics (Penning, 2017). AKR1B1, AKR1C3, and AKR1C4 encode aldo keto reductase family 1 member B, family 1 member C3, and family 1 member C4, respectively (Chen & Zhang, 2012). Overexpression of AKR1B1 was found to be associated with poor prognosis in patients with acute myelogenous leukemias and multiple myelomas (Laffin & Petrash, 2012). Results of this study also showed upregulation of AKR1B1 in bevacizumab-resistant U87 cells. This is supported by previous studies on doxorubicin-resistant U87 glioblastoma cells which showed that AKR1B1 was upregulated in doxorubicin-resistant U87 cells and that the ectopic expression of AKR1B1 inhibited doxorubicin-induced apoptosis (Han et al., 2016).

In glioma cells, hypoxia leads to upregulation of AKR1C3 (Ragel et al., 2007). AKR1C3 was found to be associated with chemoresistance by inactivation of doxorubicin and oracin (Novotna et al., 2008). AKR1C3 was also found to play a pivotal role in prostate cancer resistance to enzalutamide (Liu et al., 2015). Curcumin has been shown to inhibit AKR1B1 activity in a non-competitive manner (Puppala et al., 2010). A recent study demonstrated that PGV-1 binds to aldo keto reductase family 1 member 1 (AKR1C1) and subsequently induces reactive oxygen species elevation (Lestari et al., 2019).

GSTM1 encodes glutathione S-transferase mu 1, a phase II metabolic enzyme involved in the elimination of xenobiotics or electrophilic compounds via formation of a glutathione conjugate (McCarver & Hines, 2002). Genetic variations in GSTM1 are known to be associated with cancer, including breast (Yu et al., 2017), lung (Yu et al., 2017), and bladder cancer (Yu et al., 2017). A previous study showed that deletion in GSTM1 might be correlated with earlier age of onset of brain tumor (Wiencke et al., 1997). Recently, a meta-analysis showed no association between GSTM1 (null/present) variants and glioma risk (Liu et al., 2019).

ADAM10 encodes a disintegrin and metalloproteinase 10, an enzyme found to be involved in the development of cancer and autoimmune disease (Smith et al., 2020). Overexpression of ADAM10 was observed in patients with glioblastoma (Musumeci et al., 2015). ADAM10 promotes migration in glioblastoma cells and is therefore a potential target for glioblastoma therapy (Siney et al., 2017). ADAM10 is also involved in the cleavage of synaptic membrane proteins, which results in the formation of a soluble proteins that can provide energy for the development of a brain tumor (Endres and Deller, 2017).

Several metalloproteinases have been shown to be involved in glioblastoma resistance to bevacizumab. One study demonstrated that increased metalloproteinase activity, especially of MMP9, is responsible for treatment failure of bevacizumab in glioblastoma (Takano et al., 2010). Another study showed that increased metalloproteinase and Akt activities are associated with glioblastoma resistance to bevacizumab (Ramezani et al., 2017). These findings support the effect of PGV-1 on metalloproteinases reported in a previous study which showed that PGV-1 is able to inhibit metastasis and decreased MMP9 activity and expression in breast cancer cells (Meiyanto et al., 2019, Meiyanto et al., 2021).

PTGS2 encodes prostaglandin-endoperoxide synthase 2, also known as cyclooxygenase-2 (COX-2), an enzyme involved in the production of prostaglandin E2 (PGE2) in the inflammation process (Sakoda et al., 2005). COX-2 is not only known to be overexpressed and to regulate carcinogenesis and chemoresistance development in several type of cancer, but has also been explored as a potential target for cancer therapy (Hashemi Goradel et al., 2019). Moreover, PGE2 produced by COX-2 triggers radiation resistance in glioblastoma cells (Cook et al., 2016).

HSD17B10 encodes human type 10 17β-hydroxysteroid dehydrogenase, which is involved in the metabolism of branched-chain fatty acids and is expressed in the brain (He & Yang, 2006). Moreover, overexpression of HSD17B10 was found in patients with Alzheimer’s disease (He & Yang, 2006). Another study showed that upregulation of HSD17B10 was associated with poor chemotherapy response in patients with osteosarcoma (Salas et al., 2009).

In summary, the roles of these seven genes in cancer biology, disease progression, and anticancer drug response are supported by data from previous studies in glioblastoma and other cancers. However, the precise mechanisms by which PGV-1 targets these genes to overcome bevacizumab resistance in glioblastoma are unclear. Further studies to elucidate these mechanisms are warranted.

Immune cell infiltration is involved in cancer progression in glioblastoma (Diao et al., 2020). The resistance of glioblastoma cells to bevacizumab occurs rapidly, mainly due to modulation of the immune response and myeloid cell infiltration (Soubéran et al., 2019). Glioblastoma resistance to immunotherapy is favored by an immune microenvironment dominated by myeloid cells, including marrow-derived macrophages, microglia, myeloid-derived suppressor cells, dendritic cells, and neutrophils (De Leo et al., 2020). The tumor microenvironment causes glioblastoma cells to survive therapy and evade the T-cell response, leading to disease progression and relapse (Mohme et al., 2020). A recent study showed that changes in the arrangement of the tumor microenvironment of glioblastoma was observed in the development of resistance to antiangiogenic therapy (Ali et al., 2021). Glioblastoma resistance to bevacizumab occurs due to an immunosuppressive microenvironment triggered by VEGF (Tamura et al., 2019). Thus, modulation of the immune response is the strategy of choice for malignant tumors such as glioblastoma (Munhoz et al., 2021).

The results of this study indicated a correlation between PBR expression and immune cell infiltration. A positive correlation of B Cells was seen with ADAM10 and AKR1B1. Signaling Gαi nucleotide was positively correlated with ADAM10 maturation and activity on transitional B cells (Hwang et al., 2018). Another study showed that ADAM10 is involved in murine lupus progression on B cells (Lownik et al., 2019). Moreover, B cells are key predictors of successful treatment in patients with melanoma, in which mature and differentiated B cells are directly correlated with checkpoint inhibitor drugs (Willsmore et al., 2020). In addition, activation of B cells by a B cell based vaccine was shown to induce immunity toward glioblastoma (Lee-Chang et al., 2021). AKR1B1 is overexpressed in cancer cells and its expression is correlated with that of inflammatory mediators (Khayami, Hashemi, & Kerachian, 2020).

In this study, a positive correlation was seen between CD8+ cells and AKR1B1. The infiltration of CD4+ cells and a different CD8+ cell density were observed in recurrent glioblastoma (Anghileri et al., 2021). Moreover, the infiltration of B cells, CD8+ T cells, dendritic cells, macrophages, and neutrophils was positively associated with the prognostic risk score of glioblastomata (Wang et al., 2021). A recent study showed that immunotherapy developed for glioblastoma treatment can boost patients’ immunity by activating autologous CD8+ cells, which have the ability to eradicate glioblastoma cells (Lee-Chang et al., 2021).

A previous study observed that a decrease in the CD8+ count in glioblastoma can decrease chemotherapy response, whereas CD4+ count is associated with angiogenesis and glioblastoma progression (Mu et al., 2017). Another study showed that CD8+ cell infiltration increased due to bevacizumab, while tumor-associated macrophages decreased due to bevacizumab treatment, compared with untreated glioblastoma cells (Tamura et al., 2019). Increased activity of CD8+ is one of the strategies of cancer immunotherapy, as demonstrated in glioblastoma cells treated with silica nanoparticles (Bielecki et al., 2021). Moreover, the Inhibition of CXCL signaling leads to the prevention of myeloid-derived suppressor cell migration and increases the accumulation of CD8+ cells (Hu et al., 2021).

In this study, a positive correlation was seen between neutrophils with PTGS2 and ADAM10. Immune infiltration, including tumor-associated macrophages, neutrophils, and T-lymphocytes, was observed in patients with glioblastoma (González-Tablas Pimenta et al., 2020). In the inflammatory process, both the activation of neutrophils and the production of prostaglandin E2 (PGE2) by PTGS2 (COX-2) occur (Nannoni et al., 2020). In addition, ADAM10 plays a pivotal role in the transendothelial migration of neutrophils by cleaving ICAM-1 (Morsing et al., 2021). Taken together, various immune cell populations affect cancer progression and influence response to treatment. Further investigations to explore the role infiltration by B cells, CD8+ cells, neutrophils, dendritic cells, and tumor-associated macrophages and the effects of PGV-1 on these cells are needed.

In this study, dendritic cells showed a negative correlation with AKR1C3 and a positive correlation with PTGS2 and ADAM10. A previous study demonstrated that prostaglandin E2 induces upregulation of COX-2 or PTGS2 through the MAPK/p38 pathway in human follicular dendritic cell-like cells (Cho & Choe, 2020). The results of this study revealed that macrophages were positively correlated with ADAM10 and AKR1B1. Tumor-associated macrophages are involved in immunosuppressive mechanisms in the tumor microenvironment (Zhu et al., 2020b). Inhibition of AKR1B1 decreases the expression of inflammatory cytokines in murine macrophages (Ramana & Srivastava, 2006).

Cancer-associated fibroblasts are negatively correlated with AKR1C3 and positively correlated with PTGS2 and ADAM10. The tumor microenvironment consists of extracellular matrices, growth factors, cytokines, and CAFs, which are formed by a heterogeneous population of activated fibroblasts and play a pivotal role in cancer progression (Mochizuki et al., 2020, Rai et al., 2019). A previous study demonstrated that CAFs promote tumor development in prostate cancer cells through the upregulation of lipid biosynthesis, and inhibition of AKR1C3 could block those pathways and overcome the resistance of prostate cancer cells to anti-androgen receptor therapy (Neuwirt et al., 2020). ADAM10 is expressed by CAFs of colorectal cancer cells and is involved in cancer progression (Mochizuki et al., 2020). CAFs was shown to increase the expression of COX-2 or PTGS2 in human non-small cell lung cancer cells and thus highlighted the potential role of COX-2 as a target for preventing relapse after chemotherapy (Cho et al., 2020). Another study showed that overexpression of COX-2 in CAFs is not only correlated with poor prognosis but also induces migration and invasion in nasopharyngeal carcinoma cells (Zhu et al., 2020a). In summary, future studies on the axis of CAFs, AKR1C3, ADAM10, and PTGS2 due to PGV-1 in glioblastoma are warranted.

This study has several limitations, one of which is that this study uses a bioinformatics approach to predict the target compound of PGV-1. The prediction is made based on the algorithm built by the database. Therefore, it is necessary to conduct further studies for predictions using other machine learning to get more targets with more controlled parameters. In addition, the results of this study only produced PGV-1 targets, so the results of this study still need to be further validated in vitro, in vivo, and in clinical trials. Nevertheless, this research is very useful in accelerating the discovery of candidate protein targets to develop anticancer drugs in overcoming bevacizumab resistance. Collectively, those findings highlight the potential of PGV-1 to overcome glioblastoma resistance to bevacizumab by targeting several genes.

Conclusion

Using a bioinformatics approach, this study highlighted seven potential therapeutic target genes of PGV-1 against bevacizumab-resistant glioblastoma (PBR) namely, GSTM1, AKR1C3, AKR1C4, PTGS2, ADAM10, AKR1B1, and HSD17B10. Furthermore, this study highlights the potential of PBR as a target for immunotherapy with PGV-1. Further validation of the results of this study is required for the development of PGV-1 as adjuvant and immunotherapy against bevacizumab resistance in glioblastoma.

Funding

This research was funded by the World Class Research (WCR) Program by the Directorate General of Higher Education, Ministry of Education, Culture, Research and Technology, Republic of Indonesia, 2021. Contract Number 4518/UN1/DITLIT/DIT-LIT/PT/2021.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.