Prediction of novel target genes and pathways involved in bevacizumab-resistant colorectal cancer

Bevacizumab combined with cytotoxic chemotherapy is the backbone of metastatic colorectal cancer (mCRC) therapy; however, its treatment efficacy is hampered by therapeutic resistance. Therefore, understanding the mechanisms underlying bevacizumab resistance is crucial to increasing the therapeutic efficacy of bevacizumab. The Gene Expression Omnibus (GEO) database (dataset, GSE86525) was used to identify the key genes and pathways involved in bevacizumab-resistant mCRC. The GEO2R web tool was used to identify differentially expressed genes (DEGs). Functional and pathway enrichment analyses of the DEGs were performed using the Database for Annotation, Visualization, and Integrated Discovery(DAVID). Protein–protein interaction (PPI) networks were established using the Search Tool for the Retrieval of Interacting Genes/Proteins database(STRING) and visualized using Cytoscape software. A total of 124 DEGs were obtained, 57 of which upregulated and 67 were downregulated. PPI network analysis showed that seven upregulated genes and nine downregulated genes exhibited high PPI degrees. In the functional enrichment, the DEGs were mainly enriched in negative regulation of phosphate metabolic process and positive regulation of cell cycle process gene ontologies (GOs); the enriched pathways were the phosphoinositide 3-kinase-serine/threonine kinase signaling pathway, bladder cancer, and microRNAs in cancer. Cyclin-dependent kinase inhibitor 1A(CDKN1A), toll-like receptor 4 (TLR4), CD19 molecule (CD19), breast cancer 1, early onset (BRCA1), platelet-derived growth factor subunit A (PDGFA), and matrix metallopeptidase 1 (MMP1) were the DEGs involved in the pathways and the PPIs. The clinical validation of the DEGs in mCRC (TNM clinical stages 3 and 4) revealed that high PDGFA expression levels were associated with poor overall survival, whereas high BRCA1 and MMP1 expression levels were associated with favorable progress free survival(PFS). The identified genes and pathways can be potential targets and predictors of therapeutic resistance and prognosis in bevacizumab-treated patients with mCRC.


Introduction
Colorectal cancer (CRC) is the third most frequently diagnosed cancer and the second leading cause of cancer deaths worldwide, accounting for 10% of the worldwide cancer incidence and mortality [1]. Surgery is the treatment of choice for nonmetastatic CRC; however, approximately 20% of cases present with metastatic disease at the time of diagnosis and half of the patients experience recurrence and metastases even after complete resection of the primary tumor, leading to a poor prognosis and median overall survival (OS) of approximately 24 months [2,3]. The inclusion of cytotoxic agents (irinotecan and oxaliplatin) in fluoropyrimidine (intravenous 5-fluorouracil or oral capecitabine)-based systemic chemotherapy has been reported to improve the associated response rates (RR) from 15%-20% to 30%-40%, time to progression from 5-6 to 8 months, and OS from 10-12 to 20-24 months [3][4][5][6][7]. Furthermore, therapeutic benefits have been demonstrated to increase through the use of targeted drugs, such as angiogenesis inhibitors (bevacizumab, ziv-aflibercept, and ramucirumab) and antiepidermal growth factor receptor antibodies (cetuximab and panitumumab), as the first and second lines of treatment in patients with with K-RAS-wild-type tumors tumors [8][9][10][11][12].
Bevacizumab is the first agent to influence OS in patients with metastatic CRC (mCRC); when combined with irinotecan-based chemotherapy, the median OS improved from 15.6 to 20.3 months, median PFS from 6.2 to 10.6 months and RR from 34.8% to 44.8% [10]. The addition of bevacizumab to oxaliplatin-based chemotherapy improved median PFS from 8.0 to 9.4 months though there was no significant difference in OS(19.9 to 21.3 months) [13], while in previously treated mCRC; oxaliplatin based therapy improved both OS and PFS (10.8 to 12.9 months and 4.7 to 7.3 respectively) [14]. When compared for effectiveness, the irinotecan based chemotherapy has shown to have an edge over oxaliplatin based chemotherapy with the addition of bevacizumab (OS = 31.4 vs 30.1 months, PFS = 12.1 vs 10.7 months) [15]. These results have also been echoed in the MAVERICC trial (OS = 27.5 vs 23.9 months, PFS = 12.6 vs 10.1 months) [16]. Bevacizumab is a humanized monoclonal antibody that binds to vascular endothelial growth factor A (VEGF-A) and thus prevents interaction with its receptors, VEGFR-1 (Flt-1) and VEGFR-2 (Flk-1/KDR), leading to the regression of existing tumor blood vessels, normalization of the remaining blood vessels, and consequently tumor inhibition [17]. However, the therapeutic effects of bevacizumab are strongly affected by the lack of biomarkers that can facilitate selecting a population that might benefit from this medication and can predict therapeutic resistance [18][19][20].
In this study, we investigated the predictive biomarkers and pathways of bevacizumab resistance in mCRC by using microarray data from the Genetic Expression Omnibus (GEO) database. The new biomarkers were assessed for their ability to predict OS and PFS. The identification of predictive and prognostic biomarkers can facilitate improving the therapeutic index of bevacizumab.

Microarray data
The gene expression profile of GSE86525 was obtained from the GEO (http://www.ncbi.nlm. nih.gov/geo/) database [21], which was sequenced on the GPL16699 platform of Agilent-039494 SurePrint G3 Human GE v2 8 × 60K Microarray 039381 (Agilent Technologies, Santa Clara, CA, USA). The GSE86525 dataset includes microarray gene expression data derived from three bevacizumab-resistant HT29 xenograft tumors and three untreated HT29 xenograft tumors as controls. In brief, HT29 cells (1 × 10 7 ) suspended in phosphate-buffered saline were subcutaneously injected into the flanks of BALB/c nude mice, and the tumor-bearing mice were treated with bevacizumab (5 mg/kg, twice a week) for 3 weeks to obtain bevacizumab-resistant tumors. MTT colorimetric assays were used to determine the 50% inhibitory concentration for bevacizumab-resistant and untreated xenograft tumors; the tumor sizes were compared between the two groups. The sample tissues were immediately frozen under liquid nitrogen after isolation. Total RNAs were extracted from the samples, evaluated, labeled and hybridized, using a SurePrint G3 Human GE 8 × 60K microarray (Agilent Technologies). Array images were captured using a DNA microarray scanner (Agilent Technologies), and the data were analyzed using Feature Extraction Software (Agilent Technologies) to obtain background-corrected signal intensities. The expression data were further analyzed using Gene-Spring GX software (version 11.0, Agilent Technologies), and the differentially expressed genes (DEGs) between the bevacizumab-resistant HT29 tumors vs untreated control were compared using the Fisher exact test, followed by multiple corrections using the Benjamini and Hochberg false discovery rate (FDR) method [22]. Gene sets with an FDR q-value of <0.05 were considered statistically significant, and all experiments were performed in triplicate.

Data preprocessing and DEGs screening
The data were recalculated using the GEO2R analytical tool to identify the DEGs associated with acquired bevacizumab-resistant CRC [23,24]. The t test and Benjamini and Hochberg method were used to calculate the P values and FDR, respectively [22]. The genes were considered to be differentially expressed for an FDR value of <0.05 and fold change (FC) of >2 or <-2 (log2FC > 1 or < -1). The DEG expression data were extracted, and a bidirectional hierarchical clustering plot was constructed using MultiExperiment Viewer (MeV; version 4.8) software [25].

Construction of PPI networks
Protein-protein interaction (PPI) networks were plotted using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; version 10.0; http://www.string-db.org/), an online database comprising comprehensive known and predicted interactions, to determine the interactive relationships among the DEG-encoded proteins. A combined score of >0.7 (high confidence) was used as the cutoff criterion [26]. PPI pairs were visualized using Cytoscape software (version 3.4.0; http://www.cytoscape.org/), and the CytoNCA tool was used to subcluster the plotted PPI networks [27][28][29][30]. Highly connected proteins with important biological functions were identified by calculating the degree (number of line connections between proteins) and the betweenness value (fraction of the number of shortest paths that pass through each node; A measure of how often nodes occur on the shortest paths between other nodes) of each node with a degree cutoff criterion of !2.

Enrichment analysis of DEGs
The Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david. abcc.ncifcrf.gov/) was used to classify the DEGs involved in the PPI networks according to their biological processes, molecular functions, or cellular components by using the Gene Ontology (GO) Consortium Reference (http://www.geneontology.org/) [31,32]. Gene sets with a P value of <0.05 and FDR value of <0.05 were considered statistically significant. In addition, the DAVID tool was used for pathway enrichment analysis, and the reference pathways were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG; http:// www.genome.jp/kegg/) database website to perform KEGG pathway enrichment analysis for the DEGs involved in the PPI networks, with a P value of <0.05 and FDR value of <0.05 being considered statistically significant [33,34].

Clinical validation of the DEGs
The clinical assessment of DEGs associated with bevacizumab resistance was performed using the SurvExpress tool [ with a total of 808 cases, was used in this study. Survival profiles were compared on the basis of a high or low mRNA expression level of a particular gene, and they were censored independently for OS and PFS in months and stratified further according to TNM clinical stages 3 and 4. A log-rank P value of <0.05 was considered statistically significant, and the data were analyzed using SPSS for Macintosh (version 21, IBM Corp Armonk, NY, USA; www-01.ibm.com) for plotting Kaplan-Meier survival curves.

Gene co-expression in colorectal cancer data
The Cancer Genome Atlas (TCGA; https://cancergenome.nih.gov/) was used to obtain CRC data containing gene expression profiles. Level 3 RNASeq data containing gene expression profiles of 635 CRC cases (colon adenocarcinoma, N = 463; and rectal adenocarcinoma, N = 172) were obtained. The standard Pearson correlation coefficients (-1 to 1) and the coefficient of variation (the ratio of standard deviation to mean) of the desired gene pairs were calculated using SPSS for Macintosh (version 21, IBM Corp., Armonk, NY, USA; https://www-01.ibm.com). A P value of <0.05 was considered statistically significant and was used as the cutoff criterion.

DEGs screening and heat map clustering analysis
The GEO2R tool was used to identify DEGs from the data derived from the GPL16699 oligonucleotide microarray platform, comprising 62,976 probe sets. A total of 124 DEGs were determined to be associated with bevacizumab resistance, with 57 being upregulated and 67 being downregulated, as determined according to their log2FC and FDR values (S1 and S2 Tables). MeV software was used to construct a heat map to obtain the bidirectional hierarchical clustering of the DEGs and summarize the upregulated and downregulated DEGs (Fig 1).

Functional enrichment analysis
The DAVID tool was used to classify the DEGs involved in the PPI networks according to their common biological processes, molecular functions, or cellular components. Of the 1,454 GO gene sets included from the reference database, 111 were significantly enriched (P < 0.05; FDR < 0.05). Table 3 lists the top five gene sets, with those involved in the negative regulation of phosphate metabolic process and positive regulation of cell cycle process being the most significant and they include DUSP5, CDKN1A (p21 and Cip1), KIRREL, PDGFA, TLR4, PSMB11, BRCA1, and PBX1.

Fig 2. Protein-protein interaction (PPI) network of differentially expressed genes(A) up-regulated genes and (B)
down-regulated genes. The PPI pairs were imported into Cytoscape software as described in methods and materials. Pink nodes represent up-regulated genes while green nodes represent down-regulated genes. The lines represent interaction relationship between nodes. The highlighted DEGs represents degree = >2.

Survival analysis of the enriched DEGs
The SurvExpress tool was used to assess the enriched DEGs for their ability to predict OS and PFS in mCRC. High PDGFA expression levels were associated with poor OS, whereas high BRCA1 and MMP1 expression levels were associated with favorable PFS. However, the expression levels of CD19, CDKN1A, and TLR4 were neither associated with OS nor PFS (Figs 4 and 5).

Mechanism of gene correlation in tumor tissues
To elucidate the mechanism underlying the gene-gene correlation of the DEGs, TCGA RNA-Seq level 3 CRC data were used. BRCA1 was negatively correlated with PDGFA, CDKNA1, CD19, and TLR4 and positively correlated with MMP1. Moreover, PDGFA was negatively correlated with CDKNA1, BRCA1, MMP1, and TLR4. TLR4 was positively correlated with CDKNA1 and MMP1 and negatively correlated with CD19 and BRCA1. Furthermore, CD19 was positively correlated with CDKNA1 and negatively correlated with BRCA1, MMP1, and TLR4. However, PDGFA and CD19 were not significantly correlated (Fig 6).

Discussion
The overall mortality of CRC has remained unchanged over the past decades, despite advances in surgical and medical therapy [40,41]. This is due to the difficulties associated with early detection of the disease and the development of acquired therapeutic resistance, leading to ineffective treatment in patients with metastatic diseases [42][43][44]. Therefore, the etiological factors and mechanisms of acquired therapeutic resistance must be explored to improve   pathway, bladder cancer, and microRNAs in cancer; however, only high PDGFA expression levels were associated with poor OS, whereas high BRCA1 and MMP1 expression levels were associated with favorable PFS. These discrepancies may be because the study cohort was not specifically on bevacizumab treatment, thus suggesting that biomarkers that predict OS do not Prediction of the novel target genes and pathways involved in bevacizumab-resistant colorectal cancer specifically predict PFS. Therefore, to confidently interpret the study results, these biomarkers require further assessment in patients specifically treated with bevacizumab. The results of this study reveal PDGFA overexpression to be associated with bevacizumab resistance and the prognosis of patients with mCRC. These results are consistent with those of a previous study, which identified PDGFA as a potential predictor of therapeutic resistance and an individual prognostic marker for bevacizumab treatment, because PDGFA expression was observed to be decreased after single-dose bevacizumab treatment in responders but remained unchanged in nonresponders [47]. PDGFA targeting with the PDGF receptor has been reported to increase chemotherapeutic sensitivity in different cancers [47][48][49][50]. Therefore, our study supports the current understanding that PDGFA acts not only as a predictor of treatment response but also as a prognostic factor, because PDGFA upregulation not only limited the response to bevacizumab but also affected the prognosis of patients with mCRC in this study. Notably, PDGF overexpression has been implicated in bevacizumab resistance and poor prognosis in bevacizumab-treated patients because the PDGF pathway is considered an alternative pathway in the development of bevacizumab resistance [51, 52]. The expression levels of MMP1 and BRCA1 were associated with PFS in patients with mCRC. Although this study is the first to demonstrate the aforementioned relationship in mCRC, MMPs have received attention in terms of their role in the mechanism underlying resistance to antiangiogenic therapy, because increased MMP2 and MMP9 expression levels have been associated with resistance to the anti-VEGF and antiplacental growth factor drug aflibercept and with poor OS [53, 54]. Furthermore, MMP1 expression has been strongly associated with tumor metastasis and adverse outcomes in mCRC and has been suggested as a potential prognostic and therapeutic target [55-59]. A previous study reported that BRCA1 is associated with early onset CRC and functions as a DNA repair gene to cytotoxic drugs [60]. BRCA1 has been considered as a predictor of treatment response and prognosis in breast, ovarian, and lung cancers [61][62][63][64][65][66]; however, its role in mCRC and bevacizumab resistance is yet to be explored. The present results suggest that BRCA1 may exert protective effects in mCRC; therefore, BRCA1 should be thoroughly studied because BRCA1 targeting might not only increase the prognostic and therapeutic effects of bevacizumab but also affect the expression levels of its associated genes, namely PDGFA, CDKN1A, TLR4, and MMP1.
CD19, CDKN1A, and TLR4 have also been reported to influence therapeutic resistance or overall prognosis in cancer. CD19 has been associated with chemotherapy and multidrug resistance in many hematological tumors, and plays a central role in targeted therapeutics against B-cell malignancies (because of its expression patterns throughout the B-cell lineage), and against most B-cell malignancies with successful preclinical experiments and first-generation clinical trials [67][68][69][70][71][72]. CDKN1A has been implicated in cell cycle regulation, cell death, DNA repair, and cell motility [73]. Studies have demonstrated CDKN1A overexpression to be associated with poor prognosis in gastric and esophageal carcinomas [74,75]. Furthermore, studies have reported that TLR4 plays a role in CRC; polymorphisms increasing TLR4 signaling led to a highly aggressive CRC, whereas those reducing TLR4 signaling exerted protective effects [76,77]. In addition, high TLR4 expression levels have been associated with highly advanced grades of colonic neoplasia and with lower OS, a high probability of CRC relapse, and the presence of liver metastases in humans [78-81]. Studies have also suggested TLR4 to promote angiogenesis in different cancers by activating the PI3K-AKT signaling pathway to induce VEGF expression. In addition, TLR4 inhibition is associated with VEGF inhibition [82-84]. This finding can explain TLR4 downregulation in the bevacizumab-resistant tumors in this study; however, in vitro validation of this finding is required.
Notably, five of the six genes that were commonly enriched as well as associated with bevacizumab resistance belonged to the PI3K-AKT signaling pathway. Therefore, we suggest that the PI3K-AKT signaling pathway is responsible for restraining the therapeutic efficacy of bevacizumab in mCRC. This observation is in accordance with the results of previous studies, which have suggested that modifications in the PI3K-AKT signaling pathway increase bevacizumab resistance as an alternative pathway to VEGF inhibition [85-87]. Moreover, the occurrence of mutations in the PI3K-AKT signaling pathway remains the main challenge for mCRC treatment with new biological agents [86, 88, 89].
The present findings provide novel data that could predict bevacizumab treatment response and the emergence of resistance. Furthermore, this approach can predict patient prognosis; however, additional studies are required to validate the study findings and determine their clinical applicability.
Supporting information S1