- Primary research
- Open Access
Identification of key genes involved in tumor immune cell infiltration and cetuximab resistance in colorectal cancer
Cancer Cell International volume 21, Article number: 135 (2021)
The anti-epidermal growth factor receptor (EGFR) antibody introduces adaptable variations to the transcriptome and triggers tumor immune infiltration, resulting in colorectal cancer (CRC) treatment resistance. We intended to identify genes that play essential roles in cetuximab resistance and tumor immune cell infiltration.
A cetuximab-resistant CACO2 cellular model was established, and its transcriptome variations were detected by microarray. Meanwhile, public data from the Gene Expression Omnibus and The Cancer Genome Atlas (TCGA) database were downloaded. Integrated bioinformatics analysis was applied to detect differentially expressed genes (DEGs) between the cetuximab-resistant and the cetuximab-sensitive groups. Then, we investigated correlations between DEGs and immune cell infiltration. The DEGs from bioinformatics analysis were further validated in vitro and in clinical samples.
We identified 732 upregulated and 1259 downregulated DEGs in the induced cellular model. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses, along with Gene Set Enrichment Analysis and Gene Set Variation Analysis, indicated the functions of the DEGs. Together with GSE59857 and GSE5841, 12 common DEGs (SATB-2, AKR1B10, ADH1A, ADH1C, MYB, ATP10B, CDX-2, FAR2, EPHB2, SLC26A3, ORP-1, VAV3) were identified and their predictive values of cetuximab treatment were validated in GSE56386. In online Genomics of Drug Sensitivity in Cancer (GDSC) database, nine of twelve DEGs were recognized in the protein-protein (PPI) network. Based on the transcriptome profiles of CRC samples in TCGA and using Tumor Immune Estimation Resource Version 2.0, we bioinformatically determined that SATB-2, ORP-1, MYB, and CDX-2 expressions were associated with intensive infiltration of B cell, CD4+ T cell, CD8+ T cell and macrophage, which was then validated the correlation in clinical samples by immunohistochemistry. We found that SATB-2, ORP-1, MYB, and CDX-2 were downregulated in vitro with cetuximab treatment. Clinically, patients with advanced CRC and high ORP-1 expression exhibited a longer progression-free survival time when they were treated with anti-EGFR therapy than those with low ORP-1 expression.
SATB-2, ORP-1, MYB, and CDX-2 were related to cetuximab sensitivity as well as enhanced tumor immune cell infiltration in patients with CRC.
In the United States, colorectal cancer (CRC) is estimated to be the third-most frequently occurring malignant tumor and third leading cause of cancer-related mortality, leading to the deaths of about 53,200 people in 2020 . The incidence and mortality rates in 2015 were 37,000 cases and 19,000 deaths, with a 25- to 28-months overall survival time for Chinese patients suffering advanced CRC . Targeting epidermal growth factor receptor (EGFR), monoclonal antibodies cetuximab and panitumumab are effective therapeutics [3, 4] that shrink tumors for resectability [5,6,7] or relieve symptoms from unresectable masses by repressing tumor growth [8,9,10].
Mutations in the RAS family suggest continuous activation of EGFR downstream signaling [11,12,13,14,15] and RAS genes (KRAS, NRAS) are used as biomarkers to predict poor response to anti-EGFR therapy [16, 17]. However, tumor microenvironment re-plasticity, characterized by tumor immune cell infiltration and upregulation of immune checkpoint-related proteins, has recently been recognized as a novel explanation for drug resistance and treatment failure [18,19,20]. In particular, the interaction between tumor and immune microenvironment triggered by cetuximab or panitumumab in CRC accounts for immunogenic cell death and immune treatment resistance [21, 22].
In this study, we firstly performed high-throughput screening of transcriptomic alterations before and after induction cetuximab administration in a colon cancer cell model (CACO2). Then, we conducted integrated bioinformatics analysis of the transcriptome variations using public datasets from the Gene Expression Omnibus (GEO) database to identify several core differentially expressed genes (DEGs) and investigated whether they were correlated with tumor immune cell infiltration. Finally, we validated the identified core DEGs in vitro and in clinical samples.
Cellular model and culture
Cell lines of human CRC (CACO2, HCT116, HT29, NCIH508, and RKO) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). All cell lines were authenticated by short tandem repeat (STR) profiling (Genetic Testing Biotechnology Corporation, Suzhou, China) and routinely tested for mycoplasma using MycoAlert™ Mycoplasma Detection Kit (Lonza; LT07-218, Rockland, ME, USA). We established an EGFR antagonist-resistant cellular model using cetuximab sensitive cell line (CACO2-CS) and validated the transcriptional changes induced by cetuximab in the other cell lines. Using stepwise induction, we started with a low dose causing 50% cell growth inhibition (IC50). The dose was increased to 10 µg/mL after about two months, 50 µg/mL after another two months, and lastly 300 µg/mL. Cell Counting Kit 8 (CCK8) (Dojindo, Shanghai, China) was used to test the cell viability. For colonies formation assay, cells were cultured in six-well plates with 1000 cells and 2 ml media including cetuximab (300 µg/ml) per well. After 14 days, colonies were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Finally, the established cetuximab resistant cell line (CACO2-CR) was cultured with the maximal dose of cetuximab. CACO2-CS, HCT116, HT29, NCIH508, and RKO in the exponential growth phase were precultured in 12-well tissue culture plates for 24 h. Different concentrations of cetuximab (S20130004, Merck KGaA, Darmstadt, Germany) were added to cells which were incubated for 72 h (0 µg/mL and 25 µg/mL for CACO2, 0 µg/mL and 1 µg/mL for NCIH508, and 0 µg/mL and 50 µg/mL for HT29, HCT116, and RKO). CACO2-CS and CACO2-CR cells were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) (HyClone, Logan, UT, USA) containing 20% fetal bovine serum (FBS) (Gibco, Paisley, UK) 100 U/mL penicillin, and 100 U/mL streptomycin (Gibco, Paisley, UK). HT-29 cells were cultured in McCoy’s 5A (Gibco, Grand Island, NY, USA) supplemented with 10% FBS (Gibco, Paisley, UK), 100 U/mL penicillin, and 100 U/mL streptomycin (Gibco, Paisley, UK). HCT116 and NCIH508 cells were cultured in Roswell Park Memorial Institute Modified Medium (Hyclone) supplemented with 10% FBS (Gibco, Paisley, UK), 100 U/mL penicillin, and 100 U/mL streptomycin (Gibco, Paisley, USA). RKO cells were cultured in high glucose DMEM (Hyclone) supplemented with 10 % FBS (Gibco, Paisley, UK), 100 U/mL penicillin, and 100 U/mL streptomycin (Gibco, Paisley, UK). All cells were cultured in a humidified 5% CO2 incubator at 37 °C. All cell lines experiments were repeated at least three times with three to six replicates.
Microarray screening for transcriptional variation after cetuximab induction
Total RNA was extracted and purified with the miRNeasy Mini Kit (Cat. # 217,004 QIAGEN GmBH, Hilden, Germany) according to the manufacturer’s instructions. The RNA integrity number was determined with the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA), and the RNA was amplified and labeled with the Low Input Quick Amp WT Labeling Kit (Cat. # 5190−2943, Agilent Technologies, Santa Clara, CA, USA). Labeled cRNA was purified with the RNeasy mini kit (Cat. # 74,106). Samples were hybridized with about 1.65 µg Cy3-labeled cRNA by the Gene Expression Hybridization Kit (Cat. # 5188–5242, Agilent Technologies) in a hybridization oven (Cat. # G2545A, Agilent Technologies). After about 17 h, the slides were washed in staining dishes (Cat. # 121, Thermo Shandon, Waltham, MA, USA) with the Gene Expression Wash Buffer Kit (Cat. # 5188–5327, Agilent Technologies) and scanned by the Agilent Microarray Scanner (Cat. # G2565CA, Agilent Technologies). Data were retrieved with Feature Extraction software 10.7 (Agilent Technologies). Quantile normalization of the raw data was achieved by the limma package in R (version 3.4.1) .
Real-time quantitative PCR after cetuximab treatment in vitro
Total cellular RNA of cell lines treated with cetuximab was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) under the manufacturer’s instructions. Recombinant DNase I (Takara Bio, Beijing, China) was used to remove potential genomic DNA contamination. cDNA was generated with the PrimerScript™ RT master mix kit (Takara Bio, Dalian China). Real-time quantitative PCR analysis of the identified DEGs was performed using the TB Green Premix Ex Taq™ kit (Takara Bio, Dalian China). All results were normalized to human GAPDH mRNA expression. The primers were listed in Additional file 1: Table S1. The relative threshold cycle (Ct) method was used to display the results.
Identification of differentially expressed genes (DEGs) in the cellular model
For the microarray screening of CACO2 cells treated with cetuximab, we visualized the principal component analysis (PCA) results using the R package ggbiplot (https://github.com/vqv/ggbiplot) after converting the raw signals into gene expression levels. We set thresholds of p < 0.05 and | log2 (fold change) | > 1 to identify significant DEGs using the limma package and created heat maps using the heatmap package.
Gene expression omnibus (GEO) microarray data and DEGs analysis
To further narrow down candidate DEGs in CRC under EGFR antagonist pressure, we downloaded three gene expression datasets including cell lines and clinical tissue sample data from GEO (http://www.ncbi.nlm.nih.gov/geo). GSE59857 included the transcriptional and pharmacological profiles of 155 CRC cell lines. GSE5851 contained data from 80 clinical advanced CRC samples obtained before cetuximab monotherapy. GSE56386 contained data from eight primary CRC tumor tissue samples, comprising four from responders to cetuximab therapy and four from non-responders. From GSE59857, we downloaded the data for 20 cell lines, including 10 cetuximab-sensitive cases (OXCO2: GSM1448146, NCIH508: GSM1448142, DIFI: GSM1448175, COCM1: GSM1448167, CCK81: GSM1448097, C75: GSM1448201, HCA46: GSM1448177, C99: GSM1448204, HDC82: GSM1448128, and COGA1: GSM1448099) and 10 cetuximab-insensitive cases (SNU1047: GSM1448085, COLO320DM: GSM1448173, HUTU80: GSM1448180, KM12: GSM1448073, HDC143: GSM1448185, KM12SM: GSM1448188, KM12C: GSM1448186, COLO320: GSM1448152, KM12L4: GSM1448187, and C10: GSM1448196), among which HDC143 was deemed invalid and removed after analysis. We then converted the probe IDs into gene symbols using illuminaHumanv4.db. We retrieved the data for eight patients with the wild-type KRAS gene, including four responders to cetuximab treatment (GSM136609, GSM136593, GSM136654, and GSM136626) and four non-responders (GSM136635, GSM136646, GSM136607, and GSM136640) from GSE5851 and converted the probe IDs into gene symbols using the annotation library hgu133a2.db. The procedure for analyzing DEGs was the same as that for analyzing the CACO2 microarray data.
Gene Ontology (GO) term and Kyoto Encyclopedia of genes and genomes (KEGG) pathway enrichment analyses
For the enrichment analyses of the DEGs, a Metascape (https://metascape.org/) online analysis tool was used . The parameters were set as the following: min overlap = 3, p-value cutoff = 0.01, and min enrichment = 1.5. Related terms were selected from the top 20, according to their p-values.
Gene sets Enrichment analysis (GSEA) and Gene set variation analysis (GSVA)
GSEA was used to screen the biological states and processes associated with significantly upregulated or downregulated genes in the resistance group, with H.all.v7.1.symbols.gmt as the reference set . The number of permutations was set to 1000. A p-value < 0.05 was considered indicative of significant enrichment. A normalized enrichment score was established to evaluate the degree of enrichment. GSVA was used to describe the enrichment degree of a biological state or process in a sample . Enrichment scores were compared, and those meeting the standards of |log2 (fold change) | > 1 and p < 0.05 were considered significant results.
Construction of protein‐protein interaction (PPI) network
The PPIs of the central DEGs were generated using the STRING online database (https://string-db.org). To construct a PPI network, we retrieved cetuximab resistance candidates from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org) , confirmed their interactions using the STRING database , and finally visualized the result using Cystoscape .
Prediction of tumor immune cell infiltration
To investigate the relationship between hub DEG expression and tumor immune cell infiltration, the transcriptome landscapes of colonic or rectal adenocarcinoma samples (n = 177) in The Cancer Genome Atlas (TCGA) database were analyzed using the online tool Tumor Immune Estimation Resource Version 2.0 (https://cistrome.shinyapps.io/timer/) [30, 31].
Validation in clinical samples by immunohistochemistry staining
Patients with advanced CRC (n = 102) were diagnosed by pathological examination and administered systemic chemotherapeutics at Zhongshan Hospital of Fudan University. Fifty-two patients received therapy containing cetuximab. The progression-free survival time (PFS) was recorded to evaluate the efficacy of cetuximab. Patients signed informed consent documents, and the ethics committee of Zhongshan Hospital approved the study.
Tumor samples were fixed with 4% paraformaldehyde, embedded in paraffin, cut into sections of about 5 µm, and placed onto glass slides. After the samples were deparaffinized with xylene, hydrophilized, and unmasked, they were blocked with bovine serum albumin, immunostained with primary antibodies against SATB-2 (21307-1-AP, 1:100, Proteintech Group, Inc., Wuhan, China), ORP-1 (bs-17514R, 1:200, Bioss Biological Technology Co., Ltd., Beijing, China), MYB (bs-5978R, 1:200, Bioss Biological Technology Co., Ltd.), CDX-2 (bsm-33063m, 1:200, Bioss Biological Technology Co., Ltd.), CD8 (GB13068, 1:100, Servicebio Technology, Wuhan, China), CD19 (GB11061, 1:500, Servicebio Technology), CD4 (GB13064, 1:100, Servicebio Technology), and CD68 (GB13067-M-2, 1:100, Servicebio Technology) in a humidified environment at about 4 °C overnight and incubated with goat anti-rabbit or anti-mouse secondary antibody (1:200) for about 30 min at about 20 °C. Subsequently, the slides were stained with 3,3’-diaminobenzidine and counterstained with hematoxylin. Antigen–antibody complexes in the whole sample were detected using a panoramic slice scanner (3DHISTECH, Budapest, Hungary), recorded in a file, and viewed using CaseViewer 2.2 (3DHISTECH). To evaluate gene expression in tissues, the following formula was used to calculate the H-score using Quant Center 2.1 (3DHISTECH): H-SCORE = ∑ (PI × I) = (percentage of cells of weak intensity × 1) + (percentage of cells of moderate intensity × 2) + percentage of cells of strong intensity × 3), where PI is the proportion of the positive signal pixel area and I is the coloring intensity.
Most statistical analyses were completed using bioinformatic tools mentioned above in R (version 3.4.1). The Benjamini and Hochberg False Discovery Rate method was utilized to adjust the p- values in screening DEGs from GEO profiles. Fisher’s exact test was employed to identify the significant GO terms and KEGG pathways. The correlation significance was examined by Spearman and Pearson correlation analyses. Differential expression levels of identified DEGs were assessed by a two-tailed Student’s t-test. Kaplan - Meier survival curve and Log - rank test analysis was applied to investigated the predictive valued of identified DEGs for patients with CRC and receiving cetuximab contained therapy. The value of p < 0.050 was considered statistically significant.
Model establishment and screening of DEGs
The study design is shown in Fig. 1. After the cetuximab resistant cellular model was established successfully (Fig. 2a), preprocessing of the raw data revealed a uniform distribution of DEGs between CACO2-CS and CACO2-CR (Fig. 2b). The PCA results had acceptable reproducibility (Fig. 2c). A total of 1991 DEGs were identified, with 732 upregulated genes and 1259 downregulated genes in CACO2-CR. These DEGs are shown in a volcano plot (Fig. 2d) and a heat map (Fig. 2e) and are listed in Table 1. In GO analysis, the most enriched GO terms determined by the cellular component function (CC) were “extracellular matrix,” “apical plasma membrane” and “extracellular matrix component” (Fig. 2h). KEGG signaling analysis suggested that the DEGs were considerably enriched in “protein digestion and absorption,” “cytokine-cytokine receptor interaction,” “retinol metabolism” shown in Fig. 2i and Table 2.
Gene set enrichment
GSEA and GSVA of the DEGs (Tables 3 and 4, respectively) indicated that the downregulated genes were enriched in “HALLMARK_COMPLEMENT” (Fig. 3a), “HALLMARK_BILE_ACID_METABOLISM” (Fig. 3b), “HALLMARK_IL2_STATS_SIGNALING” (Fig. 3c), “HALLMARK_UV_RESPONSE_DN” (Fig. 3d), “HALLMARK_ESTROGEN_RESPONSE_EARLY” (Fig. 3e), “HALLMARK_APOPTOSIS” (Fig. 3f), “HALLMARK_INFLAMMATORY_RESPONSE” (Fig. 3g), and “HALLMARK_ESTROGEN_RESPONSE_LATE” (Fig. 3h). The upregulated genes were enriched in “HALLMARK_PEROXISOME” (Fig. 3i). The GSVA results suggested that the downregulated genes were enriched in “HALLMARK_MYC_TARGETS_V1 and V2” and upregulated genes were enriched in “HALLMARK_MYOGENESIS,” “HALLMARK_ANGIOGE-NESIS,” “HALLMARK_INTERFERON_ALPHA_RESPONSE,” and “HALLMARK_TGF_ BETA_SIGNALING” (Fig. 3j).
Common DEGs and PPI network construction
We identified 708 (Fig. 4a) and 298 (Fig. 4b) DEGs in GSE5851 and GSE59857, respectively. After integrating the DEGs from our own screenings with those from the two GEO datasets, 12 common DEGs were identified (Fig. 4c). These DEGs did not represent any protein-level interactions (Fig. 4d). Then, we combined recognized cetuximab resistance-related genes from the GDSC database (Table 5) with the 12 common DEGs and constructed a PPI network (Fig. 4e). The following six core DEGs were identified: SATB-2, ORP-1, MYB, CDX-2, SLC26A3, and EPHB2. Consequently, we selected the GSE56386 dataset to validate the variations in the expression levels of the 12 DEGs between cetuximab responders and non-responders in clinical terms. SATB-2, MYB, CDX-2, SLC26A3, and FAR2 were downregulated (Fig. 5a, c–e, g) and AKR1B10 was upregulated in the non-responder group (Fig. 5l). Although a trend of downregulation was observed in the other six genes (Fig. 5b, f, h–k), they exhibited no significant differences between two groups, possibly owing to the small sample size. These 12 genes had an optimal prediction accuracy in GSE56386. The receiver operating characteristic curves for the 12 genes are presented in Additional file 2: Fig. S1.
Core DEGs and tumor immune cell infiltration
We observed significant correlations between the expression levels of the core DEGs SATB-2 (Fig. 6a), ORP-1 (Fig. 6b), MYB (Fig. 6c), and CDX-2 (Fig. 6d) and tumor immune cell infiltration represented by the expression of B cell, CD4+ T cell, CD8+ T cell, and macrophage markers in TCGA. Immunohistochemical staining of a tissue microarray containing 102 CRC clinical samples (Fig. 7a) demonstrated that the expression of these four genes was positively associated with tumor-infiltrating immune cell markers such as CD4, CD8, CD19, and CD68 (Fig. 7b; Table 6).
Exploring the expression of identified core DEGs in vitro
We conducted preliminary experiments in several CRC cell lines and detected that ORP-1, MYB, and CDX-2 were downregulated in NCIH508wtRAS/RAF (a cell line sensitive to cetuximab), while SATB-2 expression decreased in CACO2 wtRAS/RAF (a cell line partially sensitive to cetuximab) and HT29 wtRAS/mtRAF (a cell line resistant to cetuximab) at the half-maximal inhibition concentration of cetuximab. A reduction in ORP-1 and CDX-2 expression was observed in most of these cell lines regardless of RAS and RAF status (Fig. 8a).
The roles of identified core DEGs in clinical samples with anti-EGFR therapy
The clinicopathological characteristics of the patients were recorded and are shown in Additional file 3: Table S2. Median follow-up time was 33.1 months (interquartile range, 17.2 to 52.5 months). Compared with patients with low ORP-1 expression, those with high ORP-1 expression in CRC experienced a significantly longer PFS time (median times were 9.0 months and 11.0 months, respectively, hazard rate [HR] = 1.901, p = 0.047) after administration of chemotherapeutics containing cetuximab (Fig. 8b).
Alterations in genes that function in the epidermal growth factor (EGF) signaling pathway, such as increases in EGFR copy number, amplification of ERBB family, overexpression of IGF1 or VEGF, or novel mutations such as point mutations in RAS, BRAF, PI3KCA, or MEK in the EGFR extracellular domain or in the downstream pathway, result in EGFR antagonist resistance. Recently, disturbances in the tumor microenvironment caused by EGFR antibody have also been recognized as factors in treatment failure. Garvey et al. reported that cetuximab causes cancer-associated fibroblasts to secrete more EGF and reactivate mitogen-activated protein kinase signaling in para-CRC cells . Critically, it is important to understand the interaction between tumor cells and the tumor microenvironment in anti-EGFR therapy resistance.
We conducted an integrated bioinformatic analysis to identify nine common DEGs between the cetuximab sensitive and resistant groups combining high-throughput data of the cetuximab resistant cell line model and data of cell lines and clinical samples from GEO profiles. The relationship between DEGs identified in our study and tumor immune cell infiltration was evaluated in TCGA and validated in clinical samples from our hospital. We found that four (SATB-2, ORP-1, MYB, and CDX-2) of nine DEGs were associated with infiltrated T cells, B cells, and macrophages in CRC. The decreasing trend of expression levels of ORP-1, MYB and CDX-2 under cetuximab pressure was observed in vitro (Additional file 4: Fig. S2), which was consistent with the expression changes in the established cetuximab resistant cell line (CACO2-CR). Moreover, patients with high expression levels of these genes, especially ORP-1, exhibited prolonged survival receiving anti-EGFR therapy.
SATB-2 encodes a DNA binding protein and mediates transcription regulation as well as chromatin remodeling. SATB-2 attenuates the activity of MEK5/ERK5 and suppresses tumorigenesis and metastasis . The upregulation of SATB-2 via DNA demethylation of the promoter region and H3K4me3 increases TH1-type chemokine expression and immune cell density in CRC . ORP-1, as a member of the oxysterol-binding protein family involve in human innate immune system, binds to phosphatidylinositol 3-phosphate by interacting with RAB7A and stabilizing GTP-RAB7A and regulates the MHC class I -mediated antigen processing and presentation pathway [35, 36]. In addition, ORP-1 suppresses tumorigenesis via metabolism-associated pathway . MYB is a protein-encoding oncogene that functions as a transcriptional activator. Paradoxically, patients with CRC and high MYB expression exhibit low incidence of distant metastases  and favorable clinical prognosis . Millen et al. demonstrated that a high level of CD8+ tumor immune infiltrating cells and a clinical history of longer relapse-free survival were related to high expression of MYB. An immunomodulatory effect conferred by MYB was also observed in CD8+ TILs in the murine CRC model . CDX-2 encodes a regulator of intestine-specific genes involved in cell growth and differentiation. CDX-2 is downregulated in the invasive part of tumor tissues and is associated with tumor-stroma protein expression as well as inflammatory cytokine release in CRC . Low levels of CDX-2 expression indicate a particularly poor survival prognosis, especially in patients with tumors that have a high stromal content .
Zanella et al. reported that EGFR antagonists inhibit colorectal tumor growth and simultaneously protect the tumor from inhibition by transcriptional regulation . Woolston et al. supported this finding and revealed that cetuximab prompts a transformation from a mutational variation to a mesenchymal transition representative of tumor-associated fibroblast enrichment . It seems that cetuximab evokes immune inflammation via antibody-dependent cell-mediated cytotoxicity but paradoxically weakens this effect via a form of IgG1 antibody-mediated immunogenic cell death  while upregulating immunosuppressive TGF-β expression in CRC . In fact, in a stage Ib/II trial of combined cetuximab and pembrolizumab treatment for patients with advanced CRC, the combination treatment of cetuximab with pembrolizumab significantly increased the density of CD3+, CD8+, and CTLA-4+ lymphocytes and natural killer cells in tumors. In the peripheral blood, the overall density of CD4+ and CD8+ lymphocytes decrease, especially that of the PD1+ memory T cells [45, 46]. In the process of acquiring secondary resistance, the tumor might experience the transition from the immune-inflamed phenotype to the immune-desert phenotype, which features by key target genes involved in tumor immune cell infiltration [47,48,49]. We believed that the core DEGs identified in this study, SATB-2, ORP-1, MYB, and CDX-2, might play critical roles in the transition.
Nevertheless, there were some limitations lying in this study. As models of anti-EGFR antibody resistance are not easy to established, only one cellular model was established and screened in the study. More cetuximab resistant cellular models with wild-type RAS and RAF genotypes are required to be established for validation. The DEGs identified using in silico methods should be further validated in vitro and in vivo. The results drawn from the retrospective clinical cohort with a small sample size also need further verification in a larger population. The specific immune cell types involved in secondary anti-EGFR antibody resistance should be more clearly elucidated.
In summary, we distinguished cetuximab-induced DEGs associated with variations of tumor immune cell infiltration in CRC by establishing a cetuximab resistant cellular model and integrated bioinformatics analysis. The results suggested that transcriptomic alterations and immune landscape remodeling should receive additional scrutiny during anti-EGFR antibody treatment. Furthermore, immunotherapy could be considered in the early stages of cetuximab treatment rather than after resistance has already occurred.
Availability of data and materials
Public datasets of colon adenocarcinoma and rectum adenocarcinoma analyzed during the current study can be retrieved from The Cancer Genome Atlas (TCGA) at https://portal.gdc.cancer.gov/. GSE5851, GSE56368, and GSE59857 datasets can be retrieved from the Gene Expression Omnibus (GEO) at https://www.ncbi.nlm.nih.gov/geo/. Other data that support the findings of this study are available from the corresponding author on reasonable request.
Epidermal growth factor receptor
Gene Expression Omnibus
Differentially expressed genes
Principal component analysis
Kyoto Encyclopedia of Genes and Genomes
Gene Sets Enrichment Analysis
Gene Set Variation Analysis
Genomics of Drug Sensitivity in Cancer
The Cancer Genome Atlas
Tumor immune infiltrating cells
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020;70(1):7–30.
Chen W, Zheng R, Baade PD, Zhang S, Zeng H, Bray F, Jemal A, Yu XQ, He J. Cancer statistics in China, 2015. CA Cancer J Clin. 2016;66(2):115–32.
Heinemann V, von Weikersthal LF, Decker T, Kiani A, Vehling-Kaiser U, Al-Batran SE, Heintges T, Lerchenmuller C, Kahl C, Seipelt G, et al. FOLFIRI plus cetuximab versus FOLFIRI plus bevacizumab as first-line treatment for patients with metastatic colorectal cancer (FIRE-3): a randomised, open-label, phase 3 trial. Lancet Oncol. 2014;15(10):1065–75.
Avallone A, Budillon A. Impact of subsequent therapies on outcome of the FIRE-3/AIO KRK0306 trial. J Clin Oncol. 2016;34(13):1564.
Primrose J, Falk S, Finch-Jones M, Valle J, O’Reilly D, Siriwardena A, Hornbuckle J, Peterson M, Rees M, Iveson T, et al. Systemic chemotherapy with or without cetuximab in patients with resectable colorectal liver metastasis: the New EPOC randomised controlled trial. Lancet Oncol. 2014;15(6):601–11.
Levi FA, Boige V, Hebbar M, Smith D, Lepere C, Focan C, Karaboue A, Guimbaud R, Carvalho C, Tumolo S, et al. Conversion to resection of liver metastases from colorectal cancer with hepatic artery infusion of combined chemotherapy and systemic cetuximab in multicenter trial OPTILIV. Ann Oncol. 2016;27(2):267–74.
Bridgewater JA, Pugh SA, Maishman T, Eminton Z, Mellor J, Whitehead A, Stanton L, Radford M, Corkhill A, Griffiths GO, et al. Systemic chemotherapy with or without cetuximab in patients with resectable colorectal liver metastasis (New EPOC): long-term results of a multicentre, randomised, controlled, phase 3 trial. Lancet Oncol. 2020;21(3):398–411.
Garattini S, Torri V, Floriani I. Cetuximab for metastatic colorectal cancer. N Engl J Med. 2009;361(1):96. author reply 96–97.
Venook AP, Niedzwiecki D, Lenz HJ, Innocenti F, Fruth B, Meyerhardt JA, Schrag D, Greene C, O’Neil BH, Atkins JN, et al. Effect of first-line chemotherapy combined with cetuximab or bevacizumab on overall survival in patients with KRAS wild-type advanced or metastatic colorectal cancer: a randomized clinical trial. JAMA. 2017;317(23):2392–401.
Adam R, Aloia T, Levi F, Wicherts DA, de Haas RJ, Paule B, Bralet MP, Bouchahda M, Machover D, Ducreux M, et al. Hepatic resection after rescue cetuximab treatment for colorectal liver metastases previously refractory to conventional systemic therapy. J Clin Oncol. 2007;25(29):4593–602.
Bregni G, Sciallero S, Sobrero A. HER2 amplification and anti-EGFR sensitivity in advanced colorectal cancer. JAMA Oncol. 2019;5(5):605–6.
Waddell T. Targeting EGFR in colorectal cancer: beyond KRAS exon 2. Lancet Oncol. 2014;15(6):540–1.
Diaz LA Jr, Williams RT, Wu J, Kinde I, Hecht JR, Berlin J, Allen B, Bozic I, Reiter JG, Nowak MA, et al. The molecular evolution of acquired resistance to targeted EGFR blockade in colorectal cancers. Nature. 2012;486(7404):537–40.
Sahin IH, Oktay AA. EGFR signalling and oncogenic pathway signature in colorectal cancer. Gut. 2013;62(11):1670.
Rachiglio AM, Lambiase M, Fenizia F, Roma C, Cardone C, Iannaccone A, De Luca A, Carotenuto M, Frezzetti D, Martinelli E, et al. Genomic profiling of KRAS/NRAS/BRAF/PIK3CA wild-type metastatic colorectal cancer patients reveals novel mutations in genes potentially associated with resistance to anti-EGFR agents. Cancers (Basel). 2019;11(6):859.
Douillard JY, Oliner KS, Siena S, Tabernero J, Burkes R, Barugel M, Humblet Y, Bodoky G, Cunningham D, Jassem J, et al. Panitumumab-FOLFOX4 treatment and RAS mutations in colorectal cancer. N Engl J Med. 2013;369(11):1023–34.
Sorich MJ, Wiese MD, Rowland A, Kichenadasse G, McKinnon RA, Karapetis CS. Extended RAS mutations and anti-EGFR monoclonal antibody survival benefit in metastatic colorectal cancer: a meta-analysis of randomized, controlled trials. Ann Oncol. 2015;26(1):13–21.
Bertotti A, Papp E, Jones S, Adleff V, Anagnostou V, Lupo B, Sausen M, Phallen J, Hruban CA, Tokheim C, et al. The genomic landscape of response to EGFR blockade in colorectal cancer. Nature. 2015;526(7572):263–7.
Misale S, Di Nicolantonio F, Sartore-Bianchi A, Siena S, Bardelli A. Resistance to anti-EGFR therapy in colorectal cancer: from heterogeneity to convergent evolution. Cancer Discov. 2014;4(11):1269–80.
Woolston A, Khan K, Cunningham D, Gerlinger M. Genomic and transcriptomic determinants of therapy resistance and immune landscape evolution during anti-EGFR treatment in colorectal cancer. Cancer Cell. 2019;36(1):35-50 e39.
Pozzi C, Cuomo A, Spadoni I, Magni E, Silvola A, Conte A, Sigismund S, Ravenda PS, Bonaldi T, Zampino MG, et al. The EGFR-specific antibody cetuximab combined with chemotherapy triggers immunogenic cell death. Nat Med. 2016;22(6):624–31.
Giordano GR, Porras A, Pancione AM. Immune resistance and EGFR antagonists in colorectal cancer. Cancers (Basel). 2019;11(8):1089.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523.
Subramanian A, Kuehn H, Gould J, Tamayo P, Mesirov JP. GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics. 2007;23(23):3251–3.
Ferreira MR, Santos GA, Biagi CA, Silva Junior WA, Zambuzzi WF. GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J Biomed Mater Res A 2020.
Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, Bindal N, Beare D, Smith JA, Thompson IR, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41:D955-61 Database issue.
Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, Doncheva NT, Legeay M, Fang T, Bork P, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021;49(D1):D605–12.
Otasek D, Morris JH, Boucas J, Pico AR, Demchak B. Cytoscape automation: empowering workflow-based network analysis. Genome Biol. 2019;20(1):185.
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–14.
Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108–10.
Garvey CM, Lau R, Sanchez A, Sun RX, Fong EJ, Doche ME, Chen O, Jusuf A, Lenz HJ, Larson B, et al. Anti-EGFR therapy induces EGF secretion by cancer-associated fibroblasts to confer colorectal cancer chemoresistance. Cancers (Basel). 2020;12(6):1393.
Mansour MA, Hyodo T, Ito S, Kurita K, Kokuryo T, Uehara K, Nagino M, Takahashi M, Hamaguchi M, Senga T. SATB2 suppresses the progression of colorectal cancer cells via inactivation of MEK5/ERK5 signaling. FEBS J. 2015;282(8):1394–405.
Xu M, Xu X, Pan B, Chen X, Lin K, Zeng K, Liu X, Xu T, Sun L, Qin J, et al. LncRNA SATB2-AS1 inhibits tumor metastasis and affects the tumor immune cell microenvironment in colorectal cancer by regulating SATB2. Mol Cancer. 2019;18(1):135.
Johansson M, Lehto M, Tanhuanpaa K, Cover TL, Olkkonen VM. The oxysterol-binding protein homologue ORP1L interacts with Rab7 and alters functional properties of late endocytic compartments. Mol Biol Cell. 2005;16(12):5480–92.
Suchanek M, Hynynen R, Wohlfahrt G, Lehto M, Johansson M, Saarinen H, Radzikowska A, Thiele C, Olkkonen VM. The mammalian oxysterol-binding protein-related proteins (ORPs) bind 25-hydroxycholesterol in an evolutionarily conserved pocket. Biochem J. 2007;405(3):473–80.
Motazacker MM, Pirhonen J, van Capelleveen JC, Weber-Boyvat M, Kuivenhoven JA, Shah S, Hovingh GK, Metso J, Li S, Ikonen E, et al. A loss-of-function variant in OSBPL1A predisposes to low plasma HDL cholesterol levels and impaired cholesterol efflux capacity. Atherosclerosis. 2016;249:140–7.
Tichy M, Knopfova L, Jarkovsky J, Pekarcikova L, Veverkova L, Vlcek P, Katolicka J, Capov I, Hermanova M, Smarda J, et al. Overexpression of c-Myb is associated with suppression of distant metastases in colorectal carcinoma. Tumour Biol. 2016;37(8):10723–9.
Tichy M, Knopfova L, Jarkovsky J, Vlcek P, Katolicka J, Capov I, Hermanova M, Smarda J, Benes P. High c-Myb expression associates with good prognosis in colorectal carcinoma. J Cancer. 2019;10(6):1393–7.
Millen R, Malaterre J, Cross RS, Carpinteri S, Desai J, Tran B, Darcy P, Gibbs P, Sieber O, Zeps N, et al. Immunomodulation by MYB is associated with tumor relapse in patients with early stage colorectal cancer. Oncoimmunology. 2016;5(7):e1149667.
Olsen J, Espersen ML, Jess P, Kirkeby LT, Troelsen JT. The clinical perspectives of CDX2 expression in colorectal cancer: a qualitative systematic review. Surg Oncol. 2014;23(3):167–76.
Sandberg TP, Sweere I, van Pelt GW, Putter H, Vermeulen L, Kuppen PJ, Tollenaar R, Mesker WE. Prognostic value of low CDX2 expression in colorectal cancers with a high stromal content—a short report. Cell Oncol (Dordr). 2019;42(3):397–403.
Zanella ER, Galimi F, Sassi F, Migliardi G, Cottino F, Leto SM, Lupo B, Erriquez J, Isella C, Comoglio PM, et al. IGF2 is an actionable target that identifies a distinct subpopulation of colorectal cancer patients with marginal response to anti-EGFR therapies. Sci Transl Med. 2015;7(272):272ra212.
Tauriello DVF, Palomo-Ponce S, Stork D, Berenguer-Llergo A, Badia-Ramentol J, Iglesias M, Sevillano M, Ibiza S, Canellas A, Hernando-Momblona X, et al. TGFbeta drives immune evasion in genetically reconstituted colon cancer metastasis. Nature. 2018;554(7693):538–43.
Boland PM, Ma WW. Immunotherapy for colorectal cancer. Cancers (Basel). 2017;9(5):50.
Le DT, Kim TW, Van Cutsem E, Geva R, Jager D, Hara H, Burge M, O’Neil B, Kavan P, Yoshino T, et al. Phase II open-label study of pembrolizumab in treatment-refractory, microsatellite instability-high/mismatch repair-deficient metastatic colorectal cancer: KEYNOTE-164. J Clin Oncol. 2020;38(1):11–9.
Inoue Y, Hazama S, Suzuki N, Tokumitsu Y, Kanekiyo S, Tomochika S, Tsunedomi R, Tokuhisa Y, Iida M, Sakamoto K, et al. Cetuximab strongly enhances immune cell infiltration into liver metastatic sites in colorectal cancer. Cancer Sci. 2017;108(3):455–60.
Wang L, Wei Y, Fang W, Lu C, Chen J, Cui G, Diao H. Cetuximab enhanced the cytotoxic activity of immune cells during treatment of colorectal cancer. Cell Physiol Biochem. 2017;44(3):1038–50.
Sunakawa Y, Yang D, Cao S, Zhang W, Moran M, Astrow SH, Hsiang J, Stephens C, Tsuji A, Takahashi T, et al. Immune-related genes to dominate Neutrophil-lymphocyte Ratio (NLR) associated with survival of cetuximab treatment in metastatic colorectal cancer. Clin Colorectal Cancer. 2018;17(4):e741-9.
This study received funding from the Outstanding Youth Foundation of Zhongshan Hospital (No. 2019ZSYQ21), CSCO-Youth Innovation Research Fund (Y-young2019-057), and National Natural Science Foundation of China (81602038, 81772511 and 81900482). The funding sources had no role in the design of the study; collection, analysis, and interpretation of data; or writing the manuscript.
Ethics approval and consent to participate
All procedures performed in this study involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. The Ethics Committee of Zhongshan Hospital Fudan University approved this study.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Sequences of primers used for real-time quantitative PCR.
ROC curves for the cetuximab sensitivity of the 12 DEGs in GSE56386. DEGs: Differentially expressed genes. ROC: Receiver operating characteristic.
Clinicopathological data of patients with colorectal cancer in the tissue microarray.
SATB-2, ORP-1, MYB, and CDX-2 were downregulated in the resistant CACO2. ##, p < 0.01, ###, p < 0.001.
About this article
Cite this article
Liang, L., Liu, M., Sun, X. et al. Identification of key genes involved in tumor immune cell infiltration and cetuximab resistance in colorectal cancer. Cancer Cell Int 21, 135 (2021). https://doi.org/10.1186/s12935-021-01829-8
- Colorectal cancer
- Drug resistance
- Anti‐epidermal growth factor receptor therapy
- Tumor immune cell infiltration
- Transcriptional alterations