Prediction of Target Genes and Pathways Associated With Cetuximab Insensitivity in Colorectal Cancer

Introduction

Colorectal cancer (CRC) is among the leading causes of cancer death both in Western Europe and East Asia and is also one of the most intensively studied diseases for kinase inhibitor therapy. ^{1 –3} Until now, the prognostic outcomes of patients with metastatic CRC have been considerably improved due to the introduction of molecular-targeted drugs, such as the angiogenesis inhibitors (bevacizumab and ramucirumab) and chemotherapy agents like oxaliplatin and fluoropyrimidines. ^{4 –6}

Cetuximab, a monoclonal antibody targeting epidermal growth factor receptor (EGFR), is also recommended for metastatic CRC. ² Despite significant improved clinical outcomes from previous studies, ^2,7 the therapeutic responses to cetuximab remain largely varied. The overall therapeutic efficacy of cetuximab is therefore limited by the lack of biomarkers that can spot cetuximab-sensitive patients and maximize the therapeutic benefits. Previously, KRAS/BRAF/PIK3CA mutation statuses, AREG/EREG expression, and Notch/Erbb2 pathways contributed to the molecular mechanisms for effective cetuximab treatment. ^{8 –11} Nonetheless, systematic identification of genes, pathways and protein–protein interaction (PPI) networks underlying the therapeutic insensitivity of cetuximab remain sparse.

With this understanding, a comprehensive in silico analysis strategy was employed for GSE56386 gene expression profile, including 4 clinical samples from responders to cetuximab and 4 from nonresponders. The differentially expressed genes (DEGs) were identified and further annotated by functional gene ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and constructed by PPI networks. The correlations between the hub genes were determined. Among the hub genes, B-cell lymphoma protein 2 (BCL2) and acetyl-coenzyme A carboxylase β (ACACB) were further externally validated by the CRC cohort of The Cancer Genome Atlas (TCGA).

Materials and Methods

Results

GO and KEGG Pathway Enrichment Analysis

The identified DEGs were analyzed with GO and KEGG pathway annotations by DAVID (Figure 1). For BP, upregulated genes were mostly involved in epidermis, epithelium development, and epithelial cell differentiation, while the downregulated genes were most enriched in digestion, lipid metabolic process, and in the single-organism catabolic process (Figure 1A). For CC, upregulated genes were mostly associated with the cornified envelope, extracellular region, and intermediate filament, while the downregulated genes were mostly associated with apical part of cell, cluster of actin-based cell projections, and brush border (Figure 1B). For MF, the upregulated genes were mostly associated with calcium ion binding, structural molecule activity, and structural constituent of cytoskeleton, while the downregulated genes were involved in 3′,5′-cyclic-adenosine monophosphate (AMP) phosphodiesterase activity, ion binding and transcription factor activity, and RNA polymerase II distal enhancer sequence-specific binding (Figure 1C). For KEGG, pathway analysis, amoebiasis, and hippo signaling pathway were enriched in upregulated genes, while pancreatic secretion and renin secretion signaling pathway were most enriched in downregulated genes (Figure 1D). In addition, similar studies were compared in terms of DEGs and pathways (Supplementary table 1).

OPEN IN VIEWER

Figure 1.

Gene ontology and KEGG analysis of differentially expressed genes associated with cetuximab insensitivity of CRC. (A) Biological function of gene ontology in upregulated/downregulated groups; (B) cellular component of gene ontology in upregulated/downregulated groups; (C) molecular function of gene ontology in upregulated/downregulated groups; (D) the KEGG pathway analysis results of differentially expressed genes in upregulated/downregulated groups. BP indicates biological function; CC, cellular component; CRC, colorectal cancer; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Protein–Protein Network Construction and Modules Selection

All the DEGs were processed by STRING. The PPI networks were constructed by the interactions results of the nodes with a degree more than 10, including a total of 200 nodes and 1333 edges (Figure 2). In addition, 3 top-scored modules within the PPI networks were identified by the MCODE with annotation by KEGG, respectively. The module 1 was mainly associated with neuroactive ligand–receptor interaction (Figure 3A and B). The module 2 was mostly associated with calcium signaling pathway (Figure 3C and D). The module 3 was mostly associated with the signaling pathways regulating pluripotency of stem cells (Figure 3E and F).

OPEN IN VIEWER

Figure 2.

Protein–protein interaction network of the DEGs. Red nodes stand for upregulated genes while the blue nodes represent downregulated genes, with the lines representing interactions between each gene. DEGs indicates differentially expressed genes.

OPEN IN VIEWER

Figure 3.

The top 3 modules extracted from the protein–protein interaction network: (A) module 1, (B) the KEGG pathway analysis of module 1, (C) module 2, (D) the KEGG pathway analysis of module 2, (E) module 3, and (F) the KEGG pathway analysis of module 3. KEGG indicates Kyoto Encyclopedia of Genes and Genomes.

Correlations of the Hub Genes in Cetuxmab-Treated and Untreated Genomic Profiles

The top 10 nodes in the PPI networks with the highest degrees were screened as hub genes, including epithermal growth factor (EGF), G-protein subunit β 5 (GNB5), G-protein subunit γ 4 (GNG4), fibroblast growth factor 2 (FGF2), BCL2, ACACB, KIT proto-oncogene receptor tyrosine kinase (KIT), adenylate cyclase 4 (ADCY4), neuropeptide Y (NPY), and neurotensin (NTS; Table 1). Previously, Kim et al reported that chemoresistant genotypes were adaptively enhanced by neoadjuvant chemotherapy in triple-negative breast cancer using single-cell sequencing while transcriptional profiles were reprogramed accordingly. ²⁴ This feature highlighted the significance of chemotherapy during the treatment time course and indicated the potential gene expression alterations do exist in response to chemotherapy. Intriguingly, given the distinct therapeutic time features among the 3 cetuximab-associated profiles, GSE56386, GSE5851, and GSE82236 (Table 2), we further evaluated the pairwise correlations of the hub genes in each profile (Figure 4). Although ADCY4 was absent in GSE5851, FGF absent in GSE65021, and 5 hub genes (GNG4, FGF2, BCL2, NPY, and NTS) were filtered due to limited expression level in GSE82236, we found distinct expression patterns between GSE56386 and GSE5851, which could partially be attributed by specimen variances. However, whether the correlation exists between the patterns and the time courses of cetuximab treatment remained further investigations. Notably, in GSE56386, the highest correlation was found between BCL2 and ACACB, which indicated potential functions associated with cetuximab insensitivity in CRC. Meanwhile, the hub genes expressions between each group in 4 profiles (GSE56386, GSE5851, GSE82236, and GSE65021) were also exhibited (Supplementary figure 1). In addition, the pairwise correlations of the hub genes were validated in GSE65021 (Supplementary figure 2).

OPEN IN VIEWER

Table 1.

The 10 Hub Genes in the PPI.

Gene Symbol	Gene Name	Degree	Betweenness Centrality	Adjacent P Value	Log FCa
EGF	Epithermal growth factor	45	0.1312	.03282	−2.72106
GNB5	G-protein subunit β 5	41	0.0264	.03282	2.53597
GNG4	G-protein subunit γ 4	41	0.0264	.03282	3.77775
FGF2	Fibroblast growth factor 2	39	0.061	.03004	−2.59538
BCL2	B-cell lymphoma protein 2	39	0.0552	.03282	−2.52730
ACACB	Acetyl-coenzyme A carboxylase β	35	0.095	.03856	−2.97898
KIT	KIT proto-oncogene receptor tyrosine kinase	34	0.0704	.03004	−3.79229
ADCY4	Adenylate cyclase 4	33	0.0345	.03979	−3.68243
NPY	Neuropeptide Y	32	0.0137	.03536	−3.20788
NTS	Neurotensin	29	0.0456	.01102	6.28008

Abbreviations: FC, fold-change;PPI, protein–protein interaction.

^a Nonresponders versus responders.

OPEN IN VIEWER

Table 2.

Comparisons of the 3 GSE Profiles.

	GSE56386	GSE5851	GSE82236
Number of samples	8	80	12
Colorectal	8	4a	12b
Liver	0	61	0
Others	0	15	0
Sample sources	Primary tumors tested in ex vivo platform for response to cetuximab	Pretreatment metastatic biopsy	Cetuximab-sensitive/resistant cell lines from HCA7, 3-dimensional culture
Exposure to cetuximab	Yes	No	Yes
Gene expression profiles type	Microarray	Microarray	RNA profiling by high-throughput sequencing
Sample date	2014	2008	2017
Platform	GPL13607	GPL571	GPL11154/GPL18573
Equipment	Agilent DNA microarray scanner	Affymetrix GeneChip scanner 3000	Illumina NextSeq 500 sequencer
Country	India	United States	United States

Abbreviation: GSE, gene set enrichment.

^a One sample was probably rectal tissue.

^b Cell lines from HCA7.

OPEN IN VIEWER

Figure 4.

The illustration of cetuximab treatment during the time courses with 3 profiles (GSE56386, GSE5851, and GSE82236) associated with the insensitivity of cetuximab in CRC and the pairwise correlation of the hub genes expression. The red circle indicated negative correlation, the blue indicated positive correlation. The values of correlation coefficients were represented by the color bar aside. Color intensity and the circle size were proportional to the correlation coefficients. CRC indicates colorectal cancer.

Network of the Correlations of BCL2 and ACACB in TCGA

To further elucidate the potential correlation among the hub genes in general CRC, the gene expressions and correlation values of the top 10 hub genes were investigated in TCGA by GEPIA platform. The significant correlated genes were screened (P value <.05; Figure 5A). Intriguingly, both BCL2 and ACACB remain the top hub genes with the highest degrees in this network, distinct from the top ranked hub genes in the primary network (Table 1). The NTS gene was not significantly correlated with the other 9 hub genes, and only 1 significant negative correlation was identified between BCL2 and GNG4. Notably, the mRNA expressions of BCL2 and ACACB in CRC were significantly reduced compared to normal tissues (Figure 5B and C). Only the mRNA expression of BCL2 was significantly varied in pathological stages (Figure 5D).

OPEN IN VIEWER

Figure 5.

Integrative analysis of BCL2 and ACACB in TCGA. (A) The correlations of hub genes based on gene expressions in TCGA; the degree was in proportion to the red color; the correlation values were illustrated by the thickness of the connecting lines. (B) The mRNA expression of BCL2 in tumor versus normal tissues (red: tumor; blue: normal); (C) the mRNA expression of BCL2 in tumor versus normal tissues (red: tumor; blue: normal); (D) the stage distribution of the mRNA expression of BCL2; and (E) the stage distribution of the mRNA expression of ACACB. BCL2 indicates B-cell lymphoma protein 2; ACACB, acetyl-coenzyme A carboxylase β; TCGA, The Cancer Genome Atlas; mRNA, messenger RNA.

Integrative Analysis of the Prognostic Analysis of BCL2 and ACACB

Both BCL2 and ACACB were not significantly associated with PFS in GSE5851 (Supplementary figure 3). To further confirm whether BCL2 and ACACB were independent prognostic factors associated with general CRC, the clinicopathological data from TCGA, containing gender, age, TNM stages, and mRNA expression of BCL2 and ACACB, were extracted for both univariate and multivariate Cox analysis. The BCL2 was, in fact, identified as prognostic significant in univariate analysis of RFS. However, both BCL2 and ACACB were not defined as independent prognostic factors in TCGA (Tables 3 and 4). Furthermore, 4 gene expression profiles, TCGA, GSE12945, GSE17536, and GSE17537 were selected for integrative analysis of the prognostic values of BCL2 and ACACB2. In fact, BCL2 was significantly associated with the prognosis (hazard ratio [HR] = 0.62, 95% confidence interval [CI], 0.30-0.95, P = .024, I ² = 68.3%) whereas ACACB was not (HR = 0.90, 95% CI, 0.49-1.31, P = .052, I ² = 61.3%; Figure 6).

OPEN IN VIEWER

Table 3.

Univariate and Multivariate Cox Analysis of Overall Survival in TCGA CRC Cohort.

OS
	Univariate Cox			Multivariate Cox
Characteristics	Hazard Ratio	95% CI	P Value	Hazard Ratio	95% CI	P Value
Gender	1.038	0.641-1.682	.879	–	–	–
Age	2.075	1.151-3.741	.015	2.798	1.523-5.141	.001
T	3.100	1.834-5.240	<.0001	1.959	1.068-3.592	.03
N	1.892	1.419-2.522	<.0001	1.377	0.833-2.277	.212
Metastasis	3.659	2.159-6.201	<.0001	1.875	0.620-5.667	.265
Stage	2.092	1.585-2.761	<.0001	1.304	0.595-2.855	.507
BCL2	0.898	0.730-1.105	.31	–	–	–
ACACB	1.018	0.804-1.289	.882	–	–	–

Abbreviations: ACACB, acetyl-coenzyme A carboxylase β; BCL2, B-cell lymphoma protein 2; CI, confidence interval; CRC, colorectal cancer; N, node; T, tumor; TCGA, The Cancer Genome Atlas.

OPEN IN VIEWER

Table 4.

Univariate and Multivariate Cox Analysis of Recurrence-Free Survival in TCGA CRC Cohort.

RFS
	Univariate Cox			Multivariate Cox
Characteristics	Hazard Ratio	95% CI	P Value	Hazard Ratio	95% CI	P Value
Gender	1.237	0.747-2.048	.409	–	–	–
Age	0.813	0.491-1.347	.421	–	–	–
T	3.751	2.075-6.781	<.0001	2.994	1.628-5.506	<.0001
N	1.766	1.299-2.401	<.0001	1.668	0.977-2.847	.061
Metastasis	4.098	2.333-7.198	<.0001	4.930	1.569-15.489	.006
Stage	1.903	1.422-2.547	<.0001	0.561	0.253-1.245	.155
BCL2	0.792	0.642-0.977	.029	0.858	0.689-1.068	.170
ACACB	1.095	0.846-1.417	.492	–	–	–

Abbreviations: ACACB, acetyl-coenzyme A carboxylase β; BCL2, B-cell lymphoma protein 2; CI, confidence interval; CRC, colorectal cancer; N, node; T, tumor; TCGA, The Cancer Genome Atlas.

OPEN IN VIEWER

Figure 6.

The integrative analysis of the prognostic values of BCL2 and ACACB in multiple gene expression profiles. (A) Integrative prognostic analysis of BCL2 and (B) integrative prognostic analysis of ACACB. BCL2 indicates B-cell lymphoma protein 2; ACACB, acetyl-coenzyme A carboxylase β.

Discussion

Colorectal cancer is one of the major cancer-related mortality causes for both Western and Eastern worlds. ^1,3 Despite the fact that a comparably large amount of patients have benefited from cetuximab therapy worldwide, the insensitive subset remain mostly indistinguishable. Therefore, prediction of novel genes and pathways associated with cetuximab insensitivity enables individualized therapeutic management and maximized the clinical outcomes.

From the original study of GSE56386, the Notch and Erbb2 signaling pathways were significantly deregulated in nonresponder tumors comparing to responders. ¹¹ In fact, the similar wtPIK3CA, BRAF, and KRAS gene signatures of the included samples lead to the enrichmented Notch and Erbb2 pathways identified by gene set enrichment analysis (GSEA). ¹¹ However, both pathways were not in the significantly enriched list of our manuscript, possible due to different bioinformatics algorithms. In fact, the distinct features of GSEA and DEG-based methodology had been previous discussed. ²⁵ Given the complex landscape and heterogeneity of cetuximab insensitivity in CRC, we further reevaluated the GSE56386 for additional functional genes and pathways, which complemented the previous findings.

In this study, a total of 1350 DEGs with 298 upregulated genes and 1052 downregulated genes were identified. Noteworthy, the different proportions of DEGs between the upregulated/downregulated clusters may indicate that downregulated genes outweighed the upregulated genes in regulating the therapeutic insensitivity of cetuximab.

Khambata-Ford et al (GSE5851) identified a total of 141 DEGs (FC > 2, P < .05) with extracellular region, serine-type endopeptidase inhibitor activity significantly enriched in GO and PPAR pathway in KEGG results. ^9,26 Consistently, extracellular region, extracellular exosome, and serine-type endopeptidase inhibitor activity were significantly identified in the upregulated GO of GSE56386, whereas PPAR signaling pathway, metabolic pathways, and steroid hormone biosynthesis were significantly enriched in the KEGG of GSE56386 (Supplementary table 1). Schütte et al identified a list of 16 genes classifiers for cetuximab responses. However, the DEGs, KEGG pathways, and GO items associated with cetuximab were not fully disclosed. ²⁷ Noteworthy, clinical criteria that defines responders and progression in diseases status remains controversial. Specifically, response evaluation criteria in solid tumors (RECIST) defines progression by the tumor exceeding 20% of initial volume while Schütte et al. proposes a relative evaluation of tumor volume versus matched untreated control group. ^27,28 These facts reflect the potential inconsistences in cetuximab treatment and public outcomes. Collectively, metabolism-associated items were consistently enriched in GO and KEGG between cetuximab responders and nonresponders (GSE56386 and GSE5851; Supplementary table 1). Noteworthy, metabolic features carry the inflection of patients’ responses and the idiosyncrasies of clinical heterogeneity. Therefore, metabolic landscape with dynamic characteristics could be intriguing targets for cetuximab treatment.

Of note, EGFR signaling had been previously validated as closely associated with development and differentiation, ²⁹ consistent with the results in BP term of upregulated genes, indicating potential association between epithelial development/differentiation and the therapeutic insensitivity of cetuximab. Furthermore, the intermediate filament, responsible for structural molecular activity and a structural constituent of cytoskeleton, significantly enriched in CC and MF terms of upregulated genes, had participated in the regulation of cellular cytoskeleton, offering a possible trait between cytoskeleton regulation and the insensitivity of cetuximab. Interestingly, previous studies highlighted a close association between a tumor suppressor, N-myc downstream-regulated gene 1, and cytoskeleton regulation as well as the ErbB family. ^30,31 Therefore, whether the direct interactions between cetuximab insensitivity and cytoskeleton regulation exist requires further experimental validation.

For the downregulated genes, the digestion/lipid metabolic process/single-organism catabolic process were enriched in BP term, while the apical part of cell/cluster of actin-based cell projections/brush border in CC term, 3′,5′-cyclic-AMP phosphodiesterase activity/ion binding, and transcription factor activity/RNA polymerase II distal enhancer sequence-specific binding functions in MF term and pancreatic section in KEGG pathway analysis.

The BCL2 was a key participator in cellular caspase signaling activation and responsible for the life-or-death switch in cellular function. ³² Previously, cetuximab had been validated as to induce the autophagy of CRC cell lines by reducing the expression of BCL2. ³³ Huang et al provided a comprehensive prognostic analysis of BCL2 in CRC with 40 qualified studies, concluded that BCL2 was significantly associated with favorable prognosis (HR = 0.69, 95% CI, 0.55-0.87, P = .002; ³⁴ Supplementary table 2). In this study, the mRNA expression of BCL2 was reduced in tumor and late pathological stages. Mechanistically, the more advanced stage of CRC with less expression of BCL2 might present less level of cetuximab-induced autophagy, possibly leading to less therapeutic responsiveness of cetuximab, which was also consistent with the findings in the PPI networks that the cetuximab insensitivity was probably associated with the overwhelmingly downregulated genes (Figure 2).

The ACACB had been intensively studied in metabolic syndrome, obesity, and diabetes diseases. Previously, ACACB was one of the hub genes in coexpression network of differentially expressed genes in colon cancer. ³⁵ Meanwhile, it was also one of the hub genes correlated with diabetes and CRC. ³⁶ Nonetheless, ACACB was also targeted by metformin, a commonly diabetes drug recently found may influence patients with CRC. ³⁷ The ACACB was not significantly associated with the prognosis of stage II CRC. ³⁸ However, other studies focused on the prognostic roles of ACACB remained limited. This was the first study to indicate a possible association between ACACB and insensitivity of cetuximab and the prognostic roles of ACACB in CRC.

Although both BCL2 and ACACB were not listed as independent prognostic factors in TCGA CRC cohort based on the inclusion criteria (Tables 3 and 4, Supplementary table 3), the meta-analysis of 4 profiles (TCGA, GSE12945, GSE17536, and GSE17537) indicated that BCL2 was significantly associated with the prognosis of CRC, consistent with the findings from Huang et al.

Remarkably, in breast cancer, both BCL2 and ACACB were featured as lower expression in 41 cases (27 relapsed) with mainly ER₋, HER2⁻, and Ki67^high, comparing to 85 (only 19 relapsed) cases with ER⁺ and Ki67^low. Moreover, both BCL2 (HR = 0.29; 95% CI, 0.17-0.49; P < .0001) and ACACB (HR = 0.32; 95% CI, 0.22-0.48; P < .0001) were significantly associated with prognosis in univariable analysis. ³⁹

The target of cetuximab, EGFR, belongs to the ErbB family including EGFR (ErbB1), HER2 (ErbB2), HER3 (ErbB3), and HER4 (ErbB4), with each being able to homodimerize or heterodimerize with the rest members. ^40,41 The EGFR is associated with HER2-4 and could initiate dimerization with the phosphorylation of the intracellular tyrosine kinase domain and further activate the downstream signaling cascades, such as the MAPK/ERK/MERK, therefore regulating the cell growth and proliferation. ^41,42 The aberration of EGFR signaling, by gene alteration or signaling malfunction, is among the most common molecular features in several cancer types, including breast cancer and CRC. ^43,44 The inhibition of EGFR by cetuximab could considerably block the aberrant activation of tumor growth. ^2,7 However, current therapeutic management and clinical outcomes of kinase inhibitor for patients with CRC remain challenging. ⁹ It is well perceived now that EGFR insensitivity has been influenced by multiple pathways. A simplified method describing the picture of EGFR signaling does not capture the full complexity of the EGFR-mediated functions in cellular level. Therefore, linear or isolated analysis of biomarkers or signal pathway may not accurately predict the therapeutic response of cetuximab treatment.

Increasing evidence indicates that DEGs enriched by GO and KEGG pathway analysis are important for the knowledge of therapeutic insensitivity of cetuximab. ^{43 –45} Thus, this study provided a systematic exploration of DEGs, PPI network, and hub genes associated with cetuximab insensitivity for patients with CRC, complementing the actionable targets spectrum. However, additional bioinformatics strategy, like GSEA analysis, and other types of cancer associated with cetuximab also contribute to the knowledge of cetuximab insensitivity (Supplementary table 4).

In order to achieve individualized management, molecular subtypes of CRC by gene expression profiling had been explored previously, focusing on 5 subtype classifications (goblet-like, enterocyte, stem-like, inflammatory, transit-amplifying). ⁴⁶ Although the included samples were small in this study, the results highlighted a multidimensional analysis process of hub genes network and subsequent clinical validation with external cohort, complementing the knowledge of insensitivity of cetuximab in CRC treatment.

However, there remain some factors that may require further clarification. First, the diversity of biopsy tissues in GSE5851 may partially confound to the influence induced by cetuximab. Indeed, potential organ-specific heterogeneity is also one of the major limitations for cetuximab treatment in different malignancies. Second, the diversity of the races in different data sets potentially confounded the conclusion. The races of the samples in GSE5851 included white (81.25%), African American (12.5%), and Asian and Others (6.25%), whereas the samples of GSE56386 came from South Asia intrinsic sensitivity or insensitivity may exit in different races. Third, the clinical criteria that defined the response to cetuximab remained inconsistent in the definition of control group. ^27,28

In fact, this manuscript only indicated the varied results between treated and untreated data sets. The impact of the treatment extent of cetuximab, for instance, the therapeutic periods and the working doses on the expression of hub genes, remained largely vacant. ⁹ Generally, a standard cetuximab regimen included 400 mg/m² loading dose and 250 mg/m² working dose per week for the first 3 weeks, dose escalation for next every 3 weeks until more than grade 2 skin rash was recorded. The median therapeutic time courses were 9 weeks. ⁹ Similarly, the ex vivo platform in GSE56386 offered consecutively 72 hours cetuximab treatment prior to the gene expression microarrays. ¹¹ In fact, dynamic alterations of hub genes may exist in tumor in response to the therapeutic management of cetuximab. On the evolutionary basis, the frequency of the cetuximab-resistant or cetuximab-insensitive phenotype may be increased (Marusyk et al, 2014).

In all, this multidimensional in silico analysis provided a novel perspective on potential therapeutic targets, pathways and mechanisms of cetuximab insensitivity in CRC. Further experimental validation is required.

Conclusion

This bioinformatics analysis provided novel insights for systematic exploration of possible target genes and pathways associated with cetuximab insensitivity. The BCL2 was associated with favorable prognosis in CRC.

Supplemental Material