Abstract

Objective: As an anti-angiogenetic treatment, bevacizumab (BVZ) therapy for GBM patients targets vascular endothelial growth factor A (VEGF), which can significantly reduce tumor sizes in the early stage based on imaging studies. However, the related mechanism of the VEGF-associated network has not yet been elucidated. Recently, subtypes of GBM BVZ response were classified and detected, so the treatment effects and related mechanisms could be reassessed individually for these different subtypes. In this study, we want to further explore possible mechanisms, effects, and side effects of BVZ treatment in GBM patients based on bioinformatics analysis. Methods: We performed multiple statistical and bioinformatics studies on several mRNA expression profiles, from hundreds of GBM patients, including an mRNA profile from 20 GBM patients before and after BVZ treatment, based on the classification of GBM patients into BVZ-responsive and BVZ- nonresponsive subtypes. Furthermore, the expression levels of various genes in the VEGF-associated and unassociated networks before and after BVZ treatment were also compared in GBM BVZ subtypes using these datasets. Results: After BVZ treatment, mRNA expression patterns and functional pathways were very different between GBM BVZ-responsive and nonresponsive subtypes. Most importantly, Angpt1, SH2D2A, and LEF1 expression levels were significantly suppressed in the tissues of GBM BVZresponsive patients but not in the tissues of GMB BVZ-nonresponsive patients. According to bioinformatics analysis, there is a network connection between these genes and VEGF;therefore, the inhibition of these genes by BVZ treatment may occur through this network of V EGF. Notably, these results suggested that some aging-related side effects,possibly involving inEPHA1 and ESM1 and related networks, may result from the treatment of BVZ-nonresponsive GBM patients. Conclusions: Ang1, SH1D2A, and LEF1were associated with BVZ treatment in GBM BVZ-responsive patients, and the effects of these genes may act through the VEGF network as VEGF-associated or unassociated pathways. The study also revealed that there may be aging-related side effects caused by the treatment for GBM BVZ-nonresponsive patients, suggesting the need for predetection of GBM BVZ- responsive subtypes. However, further research in the field is needed to better understand the roles of these genes following BVZ treatment.

Keywords: GBM; BVZ; Subtypes; Precise Treatment; VEGF; Bioinformatics; Network; Angiogenesis; Aging; Angpt1; SH2D2A; LEF1

Introduction

Glioblastoma multiforme (GBM) is the deadliest brain tumor with a median lifespan of approximately 15 months after diagnosis[1]. As GBM is a heterogeneous and multivariate disease, despite decades of investigation into the pathogenesis of GBM, the various genetic factors involved in the development and recurrence of the disease and their regulatory pathways remain poorly understood[2-5]. At present, there are many treatment methods for GBM patients, including resection, radiation, chemotherapy, immunotherapy, and antiangiogenic therapy, such as bevacizumab (BVZ) treatment (trade name Avastin or BEV). BVZ treatments target vascular endothelial growth factor A (VEGF), which can significantly reduce the early stage tumor diameters of some GBM patients based on imaging studies[6, 7].However, BVZ alone or in combination with treatment had a significant effect on the transcriptional changes and tumor size decreases in its responders, with approximately 30-33% of all GBM patients; but changes in nonresponders were small [8-10].

Angiogenesis, as a VEGF-related process, plays an important role in the development and recurrence of GBM, and despite numerous reports,it remains unclear in these processes[11]. It requires advanced and precise investigation. Recently, the classification and detection of GBM BVZ-responsive subtypes using machine learning approaches and miRNA biomarkers, including miR-21, miR-10b, and miR-197, allowed further and preciseinvestigation of functional pathways related to VEGF and angiogenesis. Compared with patients with the GBM BVZnonresponsive subtype, patients with the GBM BVZ-responsive subtype had shorter survival times, suggesting that the effect of BVZ may be better than the previously calculated; and lower levels of VEGF methylation, indicating its high expression levels and the high rate of microvascular cell proliferation. Furthermore, differentially expressed (DE) mRNAs and different functional pathways were found between these two subtypes [10, 12]. Therefore, it is worth re-evaluating and precisely exploring existing datasets in light of new classifications of regulatory angiogenesis mechanisms.

To explore the mechanisms of various genes and the VEGFassociated network in GBM, on the basis of the classification of GBM BVZ subtypes and the comparison of before and after BVZ treatment for GBM patients, multiple bioinformatics and statistical methods and mRNA expression datasets were used to analyze the functional pathways and network between these subtypes, as well as DE mRNAs.

Materials and Methods

Experimental datasets

The datasets (GSE79671) of GBM patients with BVZ treatment were downloaded from from gene expression omnibus GEO. A sequencing dataset from GBM tissue before and after BVZ treatment was trimmed with Trimmomatic and mapped to the human genome (hg19) [13]. After processing, 15630 genes were obtained and expressed as fragments per kilobase per million reads (FPKM) [9]. The data of the same patient before and after BVZ treatment were paired for statistical study.

In the Cancer Genome Atlas (TCGA) pilot study, mRNAs and mRNAs expression profiles of GBM were downloaded from TCGA portal (http://cancergenome.nih.gov/dataortal/), containing 426 glioblastoma samples, which were generated using Affymetrix U133A, Affymetrix Exon 1.0 ST, custom Agilent 244 K, and Agilent miRNA array platforms [14].

Bioinformatics and statistical methods: Multiple bioinformatics methods were performed, including Venn diagram analysis [15], STRING (version11) [16], g: Profer (ELIXIR, Tartu, Estonia, https://biit.cs.ut.ee/gprofiler/gost), and gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) [17]. In this study, the software used included MATLAB (R2022a, Math Works Inc., Natick, MA), R-project (version 4.2.2, www.rproject.org), and the statistical methods included Student’s t-test and paired t-test [18, 19].

Patients’ characteristics

This study analyzed BVZ-related network and DE gene expression using 426 GBM patient profiles from which references with complete gene profiles from microarray or high through the sequence. In Table 1, the numbers and percentages of subtypes of patients who are classified in BVZ-responsive or nonresponsive subtypes are shown, along with the classification methods used.

Table 1: Patient Characteristics.

Results

Different gene expression patterns between GBM BVZ subtypes before and after BVZ treatments

Gene expression patterns were significantly different in GBM BVZ-responsive and nonresponsive patients before and after BVZ treatment, as shown in Figure 1. The treatment response was assessed according to the RANO criteria and confirmed on the subsequent follow-up MRI [9], including complete or partial response (CR + PR).After dividing GBM patients into GBM BVZ-responsive and nonresponsive subtypes, paired t-tests before and after BVZ treatment were performed for these two subtypes. For GBM BVZresponsive patients following BVZ treatment, 469 genes were significantly expressed based on apaired t-test (p<0.05), in which 391 gene expression levels were increased, 78 gene expression levels were decreased, and their expression histogramsare shown in Figure 1A and B. In the graph, major genes are located in the nearby regions, indicating specific expression patterns. Alternatively, for GBM BVZ- nonresponsive patients following BVZ treatment, 169 genes were significantly expressed, 60 gene expression levels were increased, and 109 gene expression levels were decreased; their expression histogramsare shown in Figure1B, where there are no specific expression regions like in Figure 1A.

Figure 1: BVZ treatment varied gene expression patterns in GBM BVZ subtypes. A) The histogram shows the gene expression patterns in GBM BVZ- responsive patients, the x-axis shows the fold change before/after treatment, and the y-axis shows the gene counts for the same fold change. B) The histogram shows the gene expression patterns in GBM BVZ-NON-responsive patients, the x-axis and the y-axis are the same as above. C) The graph shows the significant functional pathways in GBM BVZ-responsive patients, according to GO and KEGG analysis using 469 DE genes obtained from paired t-tests before/after treatment, and p=0.01.

For gene ontology (GO) and KEGG pathway analysis of 469 DE mRNAs from GBM BVZ-responsive patients, 84 transcription factors (TF), 10 biological processes (BP), 5 cellular components (CC), and 3 KEGG pathways crossed the threshold (p=0.01) and were involved in BVZ treatment for GBM BVZ- responsive patients, as shown in Figure 1C. The VEGF pathway exists as one of the KEGG pathways involved. The most involved part was TF, which may be an unexpected result. In contrast, using the same analysis procedure for 169 DE mRNAs from GBM BVZ-nonresponsive patients, no functional pathways crossed the threshold (p=0.05) and were involved in BVZ treatment, suggesting that DE mRNAs expression was not specific.

This result indicated that after BVZ treatment, the gene expression patterns and related functional pathways in GBM BVZresponsive patients were specific and beneficial to the patients; in contract, in the treatment for GBM BVZ-nonresponsive patients, the gene expression patterns and related functional pathways were not specific and may not be beneficial to the patients, suggesting the need for prescreened for BVZ therapy.

VEGF-associated and -unassociated angiogenesis pathways

Angiogenesis plays a critical role in the development and recurrence of GBM with VEGF-associated and unassociated pathways. According to Venn diagram analysis, for GBM BVZresponsive patients, angiopoietin1 (Angpt1) and Sh2 domain containing 2A (SH2D2A) existed among angiogenesis (315 genes), VEGF (101 genes), and the differential expression (DE) mRNA profile (469 genes) of BVZ treatment, as shown in Figure 2A, indicating that they are VEGF-associated genes. In the graph, STRING provided 315 angiogenesis genes and 101 ranked with high confidence VEGF-associated genes, and 5 differentially expressed genes were obtained by paired t-test between pre-and post-BVZ treatment for GBM BVZ- responsive subtypes. The angiogenesis pathway existed among these genes, these two genes connected directly with VEGF with thick lines suggested confident not predicted connections, as shown in Figure 2B. Angpt1 has been reported for its complex angiogenesis function [20], but there is almost no report about the angiogenesis function of SH2D2A although it was decreased seven times (p=0.04, n=5) after BVZ treatment.

Lymphoid enhancer-binding factor 1 (LEF1), phosducinlike3 (PDCL3), and a thrombospondin type 1 domain containing 7A (THSD7A) are angiogenesis genes but not in the VEGFassociated gene group, so they belong to the VEGF-unassociated angiogenesis pathway. As shown in Figure 2B, they are not directly connected with VEGF but through other genes, THSD7A was not even found to connect with any other VEGF-related genes in this network. There were no reports about the angiogenesis function of three genes except LEF1, which was decreased approximately 1.5 times (p=0.037, n=5) after BVZ treatment. A recent study in zebrafish showed that the transcription factor LEF1 positively regulated the specification of embryonic endothelial cells during embryonic cell growth through the Wnt pathway [21]. This result suggested that LEF1 may be positively correlated with the growth of vascular endothelial cells, so BVZ treatment suppressed the expression level of LEF1, indicating that it may be beneficial to inhibit the growth of tumor blood vessels. As shown in Figure 2B, LEF1 is not directly associated with VEGF, so its effects may be VEGF indirect.

Figure 2: Angiogenesis pathways related to BVZ treatment in GBM BVZ-responsive patients. A) Venn diagram showing the relationship of angiogenesis, VEGF- associated genes, and DE mRNAs before/after BVZ treatment. B) Angiogenesis network between VEGF- associated or unassociated genes, where thick lines show confident connections but thin lines show predicted connections. C) Angpt1, SH2D2A, and LEF1 expression patterns are shown in GBM BVZ-responsive (n=130) and non-responsive (n=286) patients (p>0.05). D) Angpt1 expression patterns are shown in GBM BVZ-responsive or nonresponsive patients before and after BVZ treatment (p<0.01, n=5).

Comparing the expression levels of Angpt1, SH2D2A, and LEF1 between GBM BVZ- responsive and nonresponsive patients before treatment as shown in Figure 2C, there were no significant differences among these three genes. However, after treatment, their expression levels changed significantly, among which Angpt1 expression levels are shown in Figure 2D, suggesting these genes may serve as therapeutic indicators to monitor treatment.

Figure 3: Aging and VEGF related network in GBM BVZ-nonresponsive patients. A) Venn diagram showing the relationship of aging- associated genes, VEGF-associated genes, and DE mRNAs before/after BVZ treatment in GBM BVZ-nonresponsive patients. B) Network between VEGF, aging-associated and DE mRNAs before/after BVZ treatment in GBM BVZ-nonresponsive patients.