Systems analysis identifies potential target genes to overcome cetuximab resistance in colorectal cancer cells

Cetuximab (CTX), a monoclonal antibody against epidermal growth factor receptor, is being widely used for colorectal cancer (CRC) with wild‐type (WT) KRAS. However, its responsiveness is still very limited and WT KRAS is not enough to indicate such responsiveness. Here, by analyzing the gene expression data of CRC patients treated with CTX monotherapy, we have identified DUSP4, ETV5, GNB5, NT5E, and PHLDA1 as potential targets to overcome CTX resistance. We found that knockdown of any of these five genes can increase CTX sensitivity in KRAS WT cells. Interestingly, we further found that GNB5 knockdown can increase CTX sensitivity even for KRAS mutant cells. We unraveled that GNB5 overexpression contributes to CTX resistance by modulating the Akt signaling pathway from experiments and mathematical simulation. Overall, these results indicate that GNB5 might be a promising target for combination therapy with CTX irrespective of KRAS mutation.


Introduction
Colorectal cancer (CRC) is one of the frequently occurring cancer types in the world. The overall 5-year relative survival rate of people with localized stage in CRC is almost 90%. However, if it begins to progress, the 5-year survival rate decreases to about 20% [1]. The prognosis for patients with metastatic CRC remains poor in spite of therapeutic advancements with a median overall survival of 18-21 months [2]. Due to these, a lot of efforts have been made to identify more effective drug targets for metastatic CRC.
Cetuximab (CTX; Erbitux Ò ) is a monoclonal antibody against epidermal growth factor receptor (EGFR), being used to cure CRC patients by inhibiting EGFR activity. The EGFR signaling pathway plays a key role in multiple cellular functions such as cell growth, differentiation, survival, cell cycle progression and angiogenesis. KRAS is a downstream signaling molecule of the EGFR pathway and is mutated iñ 40% CRC patients [3]. Inhibition of EGFR pathway by CTX or other EGFR inhibitors can lead to tumor regression [4], but constitutively active mutant (MT) KRAS was found to confer the resistance to CTX [5,6]. Thus, the mutation status of KRAS was suggested as an important predictor of the benefit from CTX therapy, and CTX is currently used to those patients with wild-type (WT) KRAS. However, almost half of the patients with WT KRAS failed to respond to treatments with CTX [7], whereas some patients harboring KRAS mutation could still respond to CTX therapy [2]. These results imply that the assessment of KRAS mutation status alone is not enough to predict the CTX responsiveness.
In addition to KRAS, other predictive markers that are correlated with the response rate of CTX treatment have been studied. For instance, the downstream signaling molecules of the EGFR pathway, including BRAF, PI3K, and PTEN, were suggested to predict negative outcomes of anti-EGFR-targeted therapy [8]. However, they were not consistently associated with the response to CTX [8,9]. Recently, it was reported that CTX is more effective against left-sided (distal) CRC than right-sided (proximal) CRC, but understanding of molecular difference depending on tumor sidedness still require additional research [10,11]. Hence, there is a pressing need to identify other predictive marker to improve the responsiveness to CTX treatment.
In this study, from gene expression data analysis and cell line experiments, we have identified five potential target genes (DUSP4, ETV5, GNB5, NT5E, and PHLDA1) for CTX-resistant CRC. In KRAS WT cells, knockdown of any of these five genes increased CTX sensitivity. Intriguingly, GNB5 knockdown was sufficient to increase CTX sensitivity in KRAS MT cells. We have further investigated the role of GNB5 in cancer progression and elucidated that GNB5 contributes to CTX resistance by dominantly modulating Akt signaling pathway. Mathematical simulation of our network model supported such resistance mechanism mediated by GNB5. Taken together, we propose that GNB5 may be a promising target for combination therapy with CTX independent of KRAS mutation.

Identification of candidate target genes that can overcome cetuximab resistance
To identify the potential targets that can overcome CTX resistance, we first analyzed the gene expression data from the Khambata-Ford dataset (GSE5851) which annotated responses to CTX monotherapy for 80 metastatic CRC patients (n = 42, resistant; n = 26, sensitive; n = 12, intermediate). We identified 41 differentially expressed genes (DEGs; P-value < 0.001 by moderated t-test) from 13 162 genes between resistant and sensitive patients ( Fig. 1A and Table S1). These genes consist of 11 upregulated and 30 downregulated genes in the resistant patients. We further analyzed the gene expression data (GSE59857) which annotated responses to CTX in 155 cell lines (n = 131, resistant; n = 10, sensitive; n = 14, intermediate). Among the 11 upregulated genes in the resistant patients, seven genes showed consistently upregulated expressions (P-value < 0.05 by one-sided Wilcoxon rank-sum test) in the resistant cell lines (Fig. 1B). To validate these results, quantitative real-time PCR (qRT-PCR) for the seven upregulated genes was performed in both CTX-sensitive (NCI-H508) and CTX-resistant (HCT116, HT29, and SW48) cell lines (Fig. 1C). We found that five genes, DUSP4, ETV5, NT5E, GNB5, and PHLDA1, were significantly upregulated in CTX-resistant cell lines. We further analyzed the relationship between sidedness and the expression of the five genes using the cancer genome atlas (TCGA) clinical data of 623 CRC patients, and found that they were also significantly upregulated (P-value < 0.001 by likelihood ratio test) in right-sided CRC (n = 217), which is less responsive to CTX treatment, than left-sided CRC (n = 338; Fig. 1D). Hence, we focused on these five genes in the subsequent analysis to identify the target genes' inhibition of which can be synergistic with CTX treatment.

Recovery of cetuximab responsiveness by inhibiting the potential target genes
The five potential target genes (DUSP4, ETV5, NT5E, GNB5, and PHLDA1) were proposed to be associated with the responsiveness of CTX in some previous studies [12][13][14][15], but there was no experimental study investigating how cancer cells actually respond to CTX with the targeted inhibition of these genes. We therefore examined whether these five genes are effective combinatorial targets that can induce CTX responsiveness in CRC. We have investigated the effect of knockdown of these five genes using an RNA interference by expressing small hairpin RNAs (shRNAs) against each of them. Depending on KRAS mutation, the cellular response to such knockdown was different. In the CTX-resistant cell lines with KRAS WT, SW48, we found that knockdown of any of the five genes increase the CTX sensitivity ( Fig. 2A-E). We also found that knockdown of any of the five genes also increases susceptibility to another well-known EGFR inhibitor, erlotinib ( Fig. S1A-E).
Interestingly, for the CTX-resistant cell lines with KRAS mutation, HCT116, we found that knockdown of GNB5 increases the sensitivity to CTX (Fig. 3C).
On the other hand, there was no difference in CTX responsiveness between scrambled and DUSP4, ETV5, NT5E, or PHLDA1 knockdown in HCT116 cells, respectively (Fig. 3A,B,D,E). In the case of GNB5 knockdown in HCT116 cells, the drug sensitivity to erlotinib was also increased (Fig. S2C). In the case of DUSP4, ETV5, NT5E or PHLDA1 knockdown in HCT116 cells, there was no difference in erlotinib sensitivity (Fig. S2A,B,D,E). These results indicate that targeting GNB5 might be a new way to overcome resistance to CTX as well as other EGFR inhibitors independent of KRAS mutation.

Characteristics of GNB5 in CRC from public databases
Our results showed that GNB5 could regulate the response to CTX regardless of KRAS status. We therefore analyze public databases to expound the various features of GNB5. To examine the relationship between GNB5 expression and KRAS mutation status, we first analyzed the gene expression data (GSE59857; Fig. 4A). In CTX-resistant cell lines, GNB5 expressions was not significantly different regardless of KRAS status (n = 57, KRAS WT; n = 69, KRAS mutation), whereas CTX-sensitive cell lines (n = 10) showed lower expressions of GNB5 than CTX-resistant cell lines. Therefore, we thought that GNB5 is a predictor of CTX responsiveness in CRC independent of KRAS.
Since biomarkers are often overexpressed or mutated in cancer cells [16], we searched for genetic alterations in CRC from cBioPortal (http://cbioportal. org) [17]. As shown in Fig. 4B, there was very less amplification, mutation, or deletion of the GNB5 gene in CRC. To examine whether GNB5 is highly expressed in CRC than normal tissue, we investigated the gene expression data of 98 colorectal tumor samples to their matched adjacent normal mucosa (GSE44076) [18]. However, there was no significant difference in GNB5 expresson between them (Fig. 4C).
To examine the association between GNB5 expression and clinical outcome, we analyzed the overall survival of 623 CRC patients using TCGA clinical data. Intriguingly, the patient with high GNB5 expression showed worse overall survival among the CRC patients (Fig. 4D), indicating that elevated GNB5 expression significantly associated with a poor prognosis. With regard to the clinical relevance about tumor, node, and metastasis staging, the size of the primary tumor (T) was slightly correlated with GNB5 expression, albeit the others were not (Fig. S3). CRC are known that can be classified into four distinct consensus molecular subtypes (CMS) [19], and we found that GNB5 expression is associated with CMS1 (characterized as immune activation and worse survival after relapse) and CMS4 (characterized as mesenchymal transition and worse overall survival) (Fig. 4E). Additionally, GNB5 expression was clearly associated with microsatellite instable (MSI) tumors than microsatellite stable (MSS) tumors (Fig. 4F). To identify the signaling pathways that are regulated by GNB5, we performed gene set enrichment analysis (GSEA) using the RNA-seq data from TCGA. The genes associated with GNB5 upregulation were significantly enriched in the genes upregulated by epithelial cell proliferation, MEK, and Akt pathways (Fig. 4G). The genes associated with GNB5 upregulation were also significantly enriched in the gene sets related to hypoxia and epithelial-mesenchymal transition (Fig. 4H).

Biological role of GNB5 in CRC
GNB5 is a downstream molecule of G protein-coupled receptors (GPCRs) which are the largest receptor family that mediate physiologically important signaling [20]. One of the Gbc subunits, GNB5, transduces signals from GPCRs to various downstream effectors including Ras-MEK-ERK and PI3K-Akt pathways [21]. To investigate the signaling effect of GNB5, we examined the phosphorylation of ERK and Akt by GNB5 overexpression in NCI-H508 cells (Fig. 5A). The phosphorylation of Akt was significantly increased while the phosphorylation of ERK was modest by GNB5 overexpression. By colony-forming assay, GNB5 overexpressing cells showed a significant increase in colony number compared to control cells (Fig. 5B). We investigated the role of GNB5 in response to CTX treatment. For the CTX treatment to Knockdown efficiency of target genes was confirmed by RT-PCR (insets). GAPDH or b-actin was used as a loading control. To analyze cell death, PI was added to cells at a final concentration of 1 lgÁmL À1 . Images were taken 72 h after cetuximab treatment (C, right). **P < 0.01 by Student's t-test. GNB5-overexpressing NCI-H508 cells, we found that the cellular viability was considerably increased than control (Fig. 5C). From this result, we thought that GNB5 contributes to CTX resistance and proliferation by dominantly affecting the Akt signaling than the ERK signaling.
A mathematical model for the cetuximab resistance by GNB5 To investigate the GNB5 mechanism associated with CTX resistance using a systems biology approach [22][23][24][25][26][27][28][29][30][31][32][33][34], we have reconstructed the EGFR-GNB5 signaling network (see Materials and methods). In our regulatory network model, canonical Ras-MEK-ERK and PI3K-Akt pathways were integrated by including the interaction of GNB5 with KRAS and PI3K [35,36], as well as the cross-regulation between ERK and PI3K [37]. All the regulations in the network model were mathematically described using ordinary differential equations (ODEs). To describe the extent of cell survivability, the survival variable was defined as a linear combination of ERK and Akt activities. In the case of KRAS WT, the cell survivability of the network model was decreased by CTX treatment (Fig. 6A, left). As experimentally observed in Fig. 4A, GNB5 overexpression mainly increased the Akt activity (Fig. 6A, middle), which subsequently led to CTX resistance (Fig. 6A, right). In the case of KRAS mutation, the model network initially showed CTX resistance (Fig. 6B, left). Although GNB5 knockdown alone was not effective (Fig. 6B, middle), it showed a synergistic effect in combination with CTX by exhibiting CTX sensitivity (Fig. 6B, right). Our network model was verified by comparing its simulation results with the experimental measurements of ERK and Akt phosphorylation levels of KRAS MT SW620 cells in a timedependent manner (Fig. 6C). The experimental data were qualitatively in accord with the corresponding simulated time courses obtained from our network model (Fig. 6D,E). The simulation results for the KRAS MT model showed a comparable tendency with the experimental results of SW620 cells (Fig. 6F). We note that these results are robust to random parameter perturbations (Fig. S4) and to variations in the model structure (Fig. 6G). Overall, we suggest that GNB5 is not only a resistance marker for CTX by upregulating the Akt signaling (Fig. 6H) but also a promising target for overcoming CTX resistance independent of KRAS mutation.

Discussion
In this study, we demonstrated that the proposed predictive markers for CTX not only indicate the responsiveness of drug treatment but also provide combinatorial targets to increase the responsiveness. Inhibition of any of the five potential targets, DUSP4, ETV5, GNB5, NT5E, or PHLDA1, exhibited synergistic effects with CTX. Knockdown of DUSP4, ETV5, NT5E, or PHLDA1 could recover the CTX responsiveness in KRAS WT CRC cells, but targeting GNB5 made CRC cells become sensitive to CTX both in KRAS WT and KRAS MT CRC cells.
Since the expressions of DUSP4, ETV5, and PHLDA1 were considered as signatures of activating the RAS pathway [38,39], overexpression of these genes may sustain the RAS pathway signaling and  thereby result in resistance to CTX. There was a previous study showing that knockdown of NT5E in glioblastoma recovers the sensitivity to chemotherapy with vincristine [40]. Thus, NT5E may mediate resistance to multiple drugs including CTX.
From our cell experiments combined with mathematical simulation, we found that GNB5 determines CTX responsiveness in both KRAS WT and MT cells by modulating Akt activity. PI3K-Akt signaling pathway in KRAS MT CRC was revealed to be dominantly controlled not by KRAS but by other receptors such as IGFR [41], providing a possible rationale about why GNB5 is independent of KRAS mutation. The drug resistance mechanism called 'bypass track' means that, despite a drug inhibiting a receptor, the downstream pathway of the receptor can still be activated by alternative routes from other receptors [42]. Such a resistance mechanism has been observed in many cases with the treatment against EGFR inhibitors [43]. In case of CRC, overexpression of various receptors, such as MET, HER2, FGFR2, or IGFR, was revealed to contribute to the resistance to the EGFR inhibitor [44][45][46][47]. We contemplated that GPCRs might contribute to bypass track resistance to EGFR inhibitors, including CTX and elrotinib, via GNB5.
In addition to the signaling effect of inhibiting EGFR signaling, CTX can exert the immune effect of antibody-dependent cellular cytotoxicity (ADCC) [48]. CTX-mediated ADCC was known to induce cell apoptosis by death receptor ligands such as tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) and Fas ligand (FasL) [49,50]. KRAS mutation was found to have a resistance mechanism to CTX-mediated ADCC by blocking FasL-induced cell death in CRC [50]. Intriguingly, GNB5 was found to decrease TRAIL-induced cell death and knockdown of GNB5 restored the apoptotic response to TRAIL in CRC [51]. These previous studies show the possibility that GNB5 may inhibit CTX-mediated ADCC.
We suggest that CTX treatment will not be effective if one of the DUSP4, ETV5, GNB5, NT5E, or PHLDA1 genes is highly expressed in KRAS WT patients. Our findings indicate that inhibition of the upregulated gene among the five genes might be a promising strategy to treat CRC for combination therapy with CTX. We also suggest that CTX treatment may be effective by inhibiting GNB5 in KRAS mutated patients. Thus, GNB5 might be a novel promising target that can overcome CTX resistance independent of KRAS mutation.

Microarray data acquisition and analysis
The gene expression data of metastatic CRC patients treated with CTX monotherapy (GSE5851) were obtained from gene expression omnibus (GEO; http://www.ncbi. nlm.nih.gov/geo). DEGs between CTX-sensitive and -resistant patients were identified by moderated t-test using 'limma' package in R version 3.4.4 [52]. The gene expression data of human CRC cell lines with CTX treatment (GSE59857) were obtained from GEO and analyzed by one-sided Wilcoxon rank-sum test. The paired gene expression data of CRC and normal adjacent mucosa (GSE44076) were obtained from GEO and analyzed by moderated t-test.

TCGA expression and clinical data analysis
Clinical features associated with CRC patients were obtained from TCGA. In the Kaplan-Meier survival plot, we divided patients into two groups according to the median value of GNB5 expressions and performed Kaplan-Meier survival analysis using 'survival' package in R version 3.4.4 [53]. The log-rank test was used to assess the statistical significance. The RNA-seq data of CRC patients from TCGA was analyzed using 'edgeR' package in R version 3.4.4 [54]. The likelihood ratio test was used to assess the statistical significance. According to the information of primary tumor site for each patient, cancers derived from cecum, ascending colon, and hepatic flexure are denoted as right-sided, and from splenic flexure, sigmoid colon, descending colon, rectosigmoid junction, and rectum are denoted as left-sided.

Gene set enrichment analysis
The input datasets for GSEA analysis were prepared as follows: (a) The RNA-seq data of CRC patients were obtained from TCGA; (b) The phenotype label is determined based on the median of the GNB5 expression; (c) Gene ontology (GO) gene sets, oncogenic gene sets, and hallmarks gene sets were downloaded from MSigDB (http://www.broadinstitute. org/gsea/msigdb). Other parameters were set by their default values. GSEA analysis was performed using javaGSEA Desktop Application (http://www.broadinstitute.org/gsea/ index.jsp) with 3000 permutations.

Mathematical modeling
The ODEs of the network model were constructed mainly based on Michaelis-Menten kinetics [55] as follows: where <molecule>a or <molecule>i denote the concentration of active or inactive form of <molecule>, respectively. The total concentration of both forms of a molecule was set by <molecule>tot. Inhibition of EGFR by CTX was modeled by the competitive inhibition kinetics [56]. Survival variable is considered as a linear combination of ERKa and AKTa, resulting in a quantity with an arbitrary unit. KRAS mutation in this model was employed by decreased degradation of KRAS. Here, we note that our network model was constructed to explain the underlying regulatory mechanism in a qualitative way, and not to rigorously reproduce the quantitative amounts of actual biochemical reactions. Thus, the nominal or conditioned values of kinetic parameters and concentrations of our simplified model were simply set to represent the causal relationship between signaling molecules.
To investigate the effects of network structural variations, we included the direct link of EGFR to PI3K [57] and revised the model equation of PI3K as follows: and excluded the link between ERK-PI3K by adjusting the corresponding parameter value. All the information of parameters used in the analysis is provided in Tables S2 and S3.

Total RNA extraction and qRT-PCR
Total RNA was extracted from cells by using PureLink Ò RNA Mini Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA), according to the manufacturer's standard protocol, and treated with RNase-free DNase I (Thermo Fisher Scientific Inc.) to remove contaminating genomic DNA. cDNA was then synthesized from total RNA by reverse transcription (RT) using a DiaStar RT kit (Solgent, Daejeon, Korea). RT-PCR analysis was performed using the PCR system (Veriti 96well Thermal Cycler; Applied Biosystems, Waltham, MA, USA) with the following primers for human (Table S4). qRT-PCR analysis was performed using the QuantStudio 5 real-time PCR system (Applied Biosystems) with the corresponding primers.

Plasmid construction and virus production
To generate lentiviral particles, HEK293T cells were cotransfected with relevant lentiviral plasmid and packaging mix (pLP1, pLP2, and pLP/VSVG) using Lipofectamine (Invitrogene

Colony-forming assay
The NCI-H508 cells (200 cells) were seeded into six-well plate in triplicate and incubated at 37°C for 14 days. Cells were fixed with 3.7% PFA for 15 min, stained with a 0.5% crystal violet solution for 1 h, rinsed with distilled water, and air-dried. Cell-bound crystal violet was dissolved with 1% SDS and the absorbance measured at 590 nm using xMark TM microplate spectrophotometer (Bio-Rad, Hercules, CA, USA).