Exploring Immune-Related Prognostic Signatures in the Tumor Microenvironment of Colon Cancer

Background: Colon cancer is a common malignant tumor with poor prognosis. The aim of this study is to explore the immune-related prognostic signatures and the tumor immune microenvironment of colon cancer.

Methods: The mRNA expression data of TCGA-COAD from the UCSC Xena platform and the list of immune-related genes (IRGs) from the ImmPort database were used to identify immune-related differentially expressed genes (DEGs). Then, we constructed an immune-related risk score prognostic model and validated its predictive performance in the test dataset, the whole dataset, and two independent GEO datasets. In addition, we explored the differences in tumor-infiltrating immune cell types, tumor mutation burden (TMB), microsatellite status, and expression levels of immune checkpoints and their ligands between the high-risk and low-risk score groups. Moreover, the potential value of the identified immune-related signature with respect to immunotherapy was investigated based on an immunotherapeutic cohort (Imvigor210) treated with an anti-PD-L1 agent.

Results: Seven immune-related DEGs were identified as prognostic signatures. The areas under the curves (AUCs) of the constructed risk score model for overall survival (OS) were calculated (training dataset: 0.780 at 3 years, 0.801 at 4 years, and 0.766 at 5 years; test dataset: 0.642 at 3 years, 0.647 at 4 years, and 0.629 at 5 years; and the whole dataset: 0.642 at 3 years, 0.647 at 4 years, and 0.629 at 5 years). In the high-risk score group of the whole dataset, patients had worse OS, higher TMN stages, advanced pathological stages, and a higher TP53 mutation rate (p < 0.05). In addition, a high level of resting NK cells or M0 macrophages, and high TMB were significantly related to poor OS (p < 0.05). Also, we observed that high-risk score patients had a high expression level of PD-L1, PD-1, and CTLA-4 (p < 0.05). The patients with high-risk scores demonstrated worse prognosis than those with low-risk scores in multiple datasets (GSE39582: p = 0.0023; GSE17536: p = 0.0008; immunotherapeutic cohort without platinum treatment: p = 0.0014; immunotherapeutic cohort with platinum treatment: p = 0.0027).

Conclusion: We developed a robust immune-related prognostic signature that performed great in multiple cohorts and explored the characteristics of the tumor immune microenvironment of colon cancer patients, which may give suggestions for the prognosis and immunotherapy in the future.

Introduction

Colon cancer is known as one of the most malignant tumors with a high mortality rate worldwide (Siegel et al., 2018). Despite the recent progress in diagnosis and therapy, the overall prognosis for colon cancer patients remains poor because effective biomarkers for prognosis prediction are lacking (Keum and Giovannucci, 2019). Therefore, it is urgent and essential to explore valuable prognostic signatures and therapeutic targets for colon cancer.

Immunotherapy takes advantage of the body’s own immune system to attack cancer, which has become a powerful and promising clinical strategy for treating various tumors (Riley et al., 2019), including colon cancer(Chalabi et al., 2020; Lichtenstern et al., 2020; Ghonim et al., 2021). Immune checkpoint inhibitors (ICIs), a typical type of immunotherapy, function through inhibiting negative regulatory receptors, such as programmed cell death 1 (PD-1) and cytotoxic T lymphocyte antigen 4 (CTLA4), and thereby activates antitumor immunity (Tolba, 2020). However, only a fraction of patients were benefited from immunotherapy due to the heterogeneity and complexity of the tumor immune microenvironment (Dienstmann et al., 2017; Wang et al., 2019). Although it has been proved that IRGs were associated with the development of colon cancer (Cereda et al., 2016; Yu et al., 2019), these insights have not been applied to clinical practice. Recently, using bioinformatics and machine learning methods, various types of immune-related biomarkers have been found to be associated with the prognosis of colon cancer, such as long non-coding RNAs (Yilin Lin et al., 2020), cell infiltration (Zhou et al., 2019), and IRGs (Chen et al., 2020). However, the molecular characteristics describing the tumor immune microenvironment need to be further investigated due to their potential of prognosis and immunotherapy of colon cancer.

In this study, we constructed and validated a robust immune-related prognostic model based on TCGA-COAD cohorts and two independent GEO datasets. Additionally, we explored the relationship between the constructed prognostic model and colon cancer patients’ clinical and pathological features. In addition, we analyzed the characteristics of the tumor immune microenvironment, including tumor-infiltrating cell composition, TMB, TP53 mutation rates, and the mRNA expression levels of PD-1/PD-L1/CLTA4. Furthermore, the immune-related signature was also significantly associated with OS in patients with anti-PD-L1 treatment, and colon cancer patients with low-risk scores may be more sensitive to ICI therapy. These findings may provide new insights toward novel therapeutic targets for colon cancer.

Materials and Methods

Data Acquiring

TCGA Cohorts and the List of Immune-Related Genes

The mRNA sequencing data, mutation profiling data, and clinical information were downloaded from the UCSC Xena platform (https://xenabrowser.net/datapages/). Subsequently, the samples (n = 471) were divided into normal (n = 39) and tumor groups (n = 432), and the detailed information is shown in Supplementary Table S1. The list of immune-related genes was acquired from the ImmPort database (https://immport.niaid.nih.gov/), with a total of 1509 genes.

GEO Cohort for External Validation

Two independent datasets (GSE39582 and GSE17536) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE39582 included 556 colon cancer samples, and GSE17536 included 177 colon cancer samples, with clinical and survival information. The detailed information is shown in Supplementary Tables S2, S3, respectively.

Immunotherapeutic Cohort

An immunotherapeutic cohort (IMvigor210) was obtained from a published study (Mariathasan et al., 2018), which investigated the clinical activity of the PD-L1 blockade with atezolizumab (anti-PD-L1 McAb) in urothelial cancer. The detailed clinical information and gene expression profile of the cohort were available according to the guideline on http://research-pub.gene.com/IMvigor210CoreBiologies using the IMvigor210CoreBiologies R package. We divided the samples into platinum-treated (N = 105) and non-platinum-treated datasets (N = 237) according to whether they received platinum-based chemotherapy or not, and the detailed information is shown in Supplementary Table S4.

Screening Immune-Related DEGs

DEGs between normal and tumor groups were screened using the limma R package (Ritchie et al., 2015), with the cutoff criteria set as | log2 fold change (FC)| >0.585 and adjusted p-value < 0.05. The immune-related DEGs were obtained by overlapping the IRGs and DEGs. In order to investigate biological pathways correlated with immune response, we performed gene ontoloy (GO) functional annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on immune-related DEGs using the clusterProfiler R package (Yu et al., 2012).

Construction and Validation of the Immune-Related Prognostic Model for Colon Cancer

The whole dataset (n = 432) with all tumor samples was randomly divided into training dataset (n = 216) and test dataset (n = 216) with a 1:1 ratio. As shown in the Supplementary Table S5, there was no significant difference among the whole dataset, the training dataset, and test dataset for most clinical-pathological factors. The training dataset was used to identify the prognostic signature and constructed a prognostic risk model. First, we identified the candidate prognostic signature using the univariable Cox proportional hazards regression model and Survival R package. To avoid over-fitting, all genes with p-value < 0.05 were involved in the subsequent least absolute shrinkage and selection operator (LASSO) analysis using the glmnet R package. The association between the mRNA expression level of the filtered candidate prognostic genes and patients’ OS was further investigated using Kaplan–Meier analysis. Then, the multivariate Cox regression analysis was conducted to determine each independent prognostic indicator. Accordingly, the coefficient of the immune-related indicator was obtained from the multivariate Cox results. A formula for the immune-related risk score model was established to predict patient survival:

$\mathrm{risk}\mathrm{score=\Sigma }Coxcoefficientofgene\mathit{\chi }i\mathrm{\times }scale\mathit{exp}ressionvalueofgene\mathit{\chi }i.$

To evaluate the predictive efficiency of the constructed immune-related risk score model, we plotted a receiver operating characteristic (ROC) curve to quantify the area under the curve (AUC) using the survivalROC R package. Also, we selected the turning point of the ROC curve with the most significant difference between true positive and false positive as the optimal cutoff risk score. Patients above the cutoff value belong to the high-risk group, while patients below it belong to the low-risk group. In addition, Kaplan–Meier curves were plotted to distinguish the two groups using the survminer R package.

Moreover, the test dataset and the whole dataset were used to validate the prognostic capability of the immune-related signature. Similarly, the two datasets were divided into high- and low-risk groups based on the constructed risk score model. Next, the ROC and Kaplan–Meier curves were plotted to validate the predictive accuracy of the risk score model. Then, the nomogram was constructed using the whole dataset based on the risk score model and clinical factors, including the age, gender, microsatellite status, and tumor stage. The constructed nomogram was further assessed by calibration. Additionally, the associations between the immune-related constructed risk score model and the clinical and pathological characteristics, including advanced pathological stages and TNM stages, were explored by using the Wilcoxon test. Additionally, the constructed model was further validated using GEO datasets with accession numbers GSE39582 and GSE17536.

Estimation and Comparison of Tumor-Infiltrating Immune Cell Type Fractions

The whole dataset was divided into high- and low-risk groups according to the constructed risk model, and the CIBERSORT algorithm was conducted to access the proportions of 22 types of tumor-infiltrating immune cells using the normalized gene expression matrixes and running with 1000 permutations (Newman et al., 2015). Subsequently, the comparison of immune landscape between the high- and low-risk groups was evaluated using the unpaired t-test. The significant differential immune cell types (p-value < 0.001) were further assessed for their relationship with OS using Kaplan–Meier curves.

Characteristics of Immunotherapy-Related Predictors for Colon Cancer Patients

We first calculated the TMB value and visualized the mutation profiles of the high- and low-risk groups of the whole dataset using the maftools R package (Mayakonda et al., 2018). The unpaired t-test statistically analyzed the differences of the TMB between the high- and low-risk groups. In addition, the OS between the high- and low-risk groups was calculated using the Kaplan–Meier method. Moreover, the Wilcoxon test was used to compare the mRNA levels of immune checkpoints and their ligands between the high- and low-risk groups.

Exploring the Associations Between the Microsatellite Status and the Constructed Prognostic Model

The whole dataset, after removing the samples without microsatellite status information, was used for further analysis based on the constructed prognostic model. Subsequently, the samples were divided into MSI-H and MSS/MSI-L groups according to the microsatellite status information extracted from the phenotypic data, and the Wilcoxon test was performed to analyze the difference of the level of risk score between the MSI-H and MSS/MSI-L groups. Moreover, the OS between the MSI-H group and MSI-L/MSS group was calculated using the Kaplan–Meier method.

The Role of the Immune-Related Prognostic Signature in Immunotherapy

In order to investigate the potential value of the identified immune-related signature with respect to immunotherapy, we obtained the gene expression profiles and corresponding clinical features from an immunotherapeutic cohort (Imvigor210) treated with anti-PD-L1 agent, which was divided into platinum-treated and non-platinum-treated datasets. We first validated the constructed immune-related prognostic model using the platinum-treated and non-treated datasets, respectively. Then, the complete response (CR) or partial response (PR) patients were categorized as responders and compared with non-responders, who displayed stable (SD) or progressive disease (PD), and the risk score of each patient was calculated based on the constructed risk score model. Subsequently, we statistically analyzed the distribution of the risk score between the responders and non-responders. Moreover, we further compared the tumor mutation load and neoantigen burden between high-and low-risk groups using the Wilcoxon test.

Results

Identification of Immune-Related DEGs

A flow chart of the whole analysis pipeline is shown in Figure 1. A total of 571 DEGs (275 upregulated and 296 downregulated) were screened by comparing between tumor and normal groups (Figure 2A). After the intersection with 1509 IRGs, 102 immune-related DEGs were obtained (Figure 2B), of which 83 genes were downregulated, and 19 genes were upregulated. Detailed information is shown in Supplementary Table S6. Subsequently, functional and pathway enrichment analyses were performed using the clusterProfiler R package. KEGG analysis results indicated that the immune-related DEGs were significantly enriched in terms associated with the cytokine–cytokine receptor interaction, neuroactive ligand–receptor interaction, and IL-17 signaling pathway (Figure 2C), while GO related to humoral immune response was mediated by circulating immunoglobulin, humoral immune response, and immunoglobulin-mediated immune response (Figure 2D).

FIGURE 1

FIGURE 1. The whole flow chart of data analysis.

FIGURE 2

FIGURE 2. Analysis of immune-related DEGs and construction of the immune-related prognostic model. (A) The volcano plot on all DEGs between the tumor and normal samples. The green dots represent downregulated genes, while the red dots represent upregulated genes. (B) The Venn diagram of the intersection between the DEGs and IRGs. (C) The KEGG pathway enrichment analysis of immune-related DEGs. (D) GO analysis of immune-related DEGs. (E) The forest plot showed the results of multivariate Cox analysis. (F) The distribution of the high- and low-risk score groups and their relationship with OS, and the expression pattern of seven prognostic signatures in high- and low-risk score groups. (G) The Kaplan–Meier curve revealed that OS in the low-risk score group was significantly higher than that in the high-risk score group. (H) Time-dependent ROC curve analysis of the immune-related risk score model.

Construction of the Immune-Related Risk Score Model and the Evaluation of its Prognostic Ability

To explore the prognostic value of the immune-related DEGs, we performed the univariate Cox regression analysis. A total of 10 genes were significantly related to OS status, and 7 genes with the maximum prognostic value were further identified using LASSO regression analysis (Supplementary Figures S1A,B). The mRNA expression level of the seven genes was significantly associated with patients’ OS (GUCA2A: p = 0.013; CHGB: p = 0.05; SSTR2: p = 0.017; VIP: p = 0.0074; OXTR: p = 0.001; IL1A: p = 0.0035; and GRP: p = 0.016), and the higher expression level of IL1A and GUCA2A was associated with a better patients’ OS, while the other five genes were opposite (Supplementary Figure S1C–I). Then, we conducted the multivariate Cox regression analysis and established an immune-related risk score model based on the training dataset, and the hazard ratio of each gene is shown in Figure 2E. The colon cancer patients were divided into high- and low-risk groups according to the risk score calculated using the formula described in Materials and Methods. The scatter diagram in Figure 2F revealed that the OS tended to become worse with the increase of risk score, and the proportion of death in the high-risk group (the proportion of red dot and blue dot on the right side) was higher than that in the low-risk group. The heatmap in Figure 2F showed that the expression of IL1A and GUCA2A was low in the low-risk group and high in the high-risk group, while the trend of the other five genes was opposite. The Kaplan–Meier analysis results showed that high-risk score patients had worse OS than low score patients (p < 0.0001, Figure 2G). The prognostic accuracy of the risk score model was investigated as a continuous variable (Figure 2H). The AUC of the prognostic model for OS was 0.780 at 3 years, 0.801 at 4 years, and 0.766 at 5 years, indicating its excellent prediction performance.

Validation and Assessment of the Immune-Related Prognostic Signatures

To determine if the constructed risk core model is consistent in different populations, we performed an identical formula using the test dataset and the whole dataset. Consistent with the findings in the training dataset, patients categorized into the high-risk score group had worse OS than the patients in the low-risk score group (p < 0.05, Supplementary Figures S2A,C). The areas under the curves (AUCs) of the prognostic model were 0.642 for 3-year OS, 0.647 for 4-year OS, and 0.629 for 5-year OS using the test dataset, and 0.626 for 3-year OS, 0.663 for 4-year OS, and 0.661 for 5-year OS using the whole dataset (Supplementary Figures S2B,D). The Wilcoxon test showed that the higher risk score was associated with a higher T stage (p = 0.00009), N stages (p = 0.0018), metastasis (p = 0.0064), and advanced pathological stage (p = 0.0034) based on the whole dataset (Figures 3A–D).

FIGURE 3

FIGURE 3. Exploring the relationship between the risk score of the colon cancer patients and clinical and pathological characteristics, including (A) T stages, (B) M stages, (C) N stages, and (D) advanced pathological stages, based on the whole dataset. Construction and validation of a nomogram. (E) Forest plots showed the associations between patients’ characteristics and OS. (F–H) The calibration plot of the nomogram to predict the probability of OS at 3, 4, and 5 years.

To further validate the robustness of the prognostic signatures and improve the accuracy of the performance of the risk score model, we constructed a nomogram that integrated the immune-related risk score and clinical information, including the age, sex, microsatellite status, and tumor stage to quantitatively predict the prognosis of colon cancer patients in the whole dataset. In the nomogram, the score for each variable can be found on the point scale, so that it is easy to estimate the probability of survival at 3, 4, and 5 years by calculating the total score (Supplementary Figure S2E). The forest plot showed that patient’s characters, including the age (>60), tumor stage (III and IV), and risk score were significantly associated with the OS (p-value < 0.05, Figure 3E). The calibration curves revealed that the predictive curves were close to the ideal curve (Figures 3F–H), indicating good performance. Furthermore, the predictive accuracy of this nomogram (C-index: 0.74) was higher than that of the risk score model (C-index: 0.72).

Exploring the Tumor Immune Microenvironment in Colon Cancer Patients

Based on the CIBERSORT algorithm, we estimated the proportions of 22 types of immune cells in each colon cancer patient. Then, we compared the proportions of immune cells between the low-risk group and high-risk group, and the significant differences were found in resting NK cells, M0 macrophages, M2 macrophages, CD4 memory-activated T cells, plasma cells, resting mast cells, and neutrophils. Among them, the resting NK cells and M0 macrophages were the most significant, with p < 0.001 (Figure 4A). The Kaplan–Meier curve revealed that a high level of resting NK cells or M0 macrophages was significantly related to poor OS (p < 0.05, Figures 4B,C).

FIGURE 4

FIGURE 4. Analysing the immune cell types and mutation profiles in high- and low-risk groups based on the whole dataset. (A) Comparing the difference of the proportions of immune cells between the low-risk group and high-risk group using the Wilcoxon test. The values of P were labeled above each boxplot with asterisks (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). (B, C) The Kaplan–Meier analysis of the relationship between the level of resting NK cells and M0 macrophages with patients’ OS. (D) The mutation profiles of colon cancer patients in high- and low-risk groups. (E) Comparison of the mutation rate between the high-risk group and low-risk group.

The mutation profiles of each colon cancer patient were plotted using the whole dataset. As shown in Figure 4D, the top 20 significantly mutated genes were APC, TP53, TTN, KRAS, SYNE1, MUC16, PIK3CA, FAT4, RYR2, OBSCN, ZFHX4, DNAH5, PCLO, CSMD3, ABCA13, DNAH11, LRP1B, FAT3, USH2A, and CSMD1. Among them, the mutation rate of TP53 was significantly different between the high-risk score group and low-risk group (p = 0.037, Figure 4E). However, the TP53 status was not significantly related to patients’ OS (Figure 5A). Besides, we calculated the TMB of each sample and found that there was no significant difference between the high-risk group and the low-risk group (p = 0.85, Figure 5B). However, we observed that high TMB was significantly related to poor patients’ OS (Figure 5C). Additionally, the Wilcoxon test statistically analyzed the difference in the level of risk scores between the MSI-H and MSI-L/MSS groups, and the result showed the difference was not significant (p = 0.06, Figure 5D). As shown in Figure 5E, the microsatellite status cannot be used as an independent prognostic factor (p = 0.83).

FIGURE 5

FIGURE 5. Exploring the tumor immune microenvironment in colon cancer patients. (A) The association between the TP53 status and patients’ OS. (B) The difference of TMB between the high-risk group and the low-risk group. (C) The association between TMB and patients’ OS. (D) The difference in the level of the risk score between the MSI-H and MSI-L/MSS groups. (E) The association between the microsatellite status and patients’ OS. (F–H) Comparison of the expression levels of the immune checkpoints and their ligands between the high-risk score group and low-risk score group. (F) The expression of PD-L1, (G) the expression of PD-1, and (H) the expression of CTLA-4.

Moreover, the Wilcoxon test was used to compare the expression levels of the immune checkpoints and their ligands between the high-risk score group and low-risk score group. The high-risk score group had a high expression level of PD-L1 (p = 0.001), PD-1(p = 0.00022), and CTLA-4 (p = 0.00021, Figures 5F–H).

Validation of the Prognostic Signature Using GEO Datasets

In order to investigate the applicability in multiple cohorts based on different platforms, we further verified the constructed risk score model using two independent external GEO datasets. We found that patients with high-risk scores demonstrated worse prognosis than those with low-risk scores (GSE39582: p = 0.0023, Supplementary Figure S3; GSE17536: p = 0.0008, Supplementary Figure S4). The AUCs of the prognostic model were 0.577 for 3-year OS, 0.569 for 4-year OS, and 0.568 for 5-year OS using the GSE39582 dataset (Supplementary Figure S3A), and 0.578 for 3-year OS, 0.699 for 4-year OS, and 0.657 for 5-year OS using GSE17536 (Supplementary Figure S4A). The Wilcoxon test showed that a higher risk score was associated with a higher T stage (p = 0.0041), metastasis (p = 0.04), N stages (p = 0.037), and advanced pathological stages (p = 0.033) using GSE39582 (Supplementary Figures S3C–F). It was also found that a higher risk score was associated with a higher advanced pathological stage (p = 0.029) in GSE17536 (Supplementary Figure S4C), without obtaining the TMN stage data.

The Prognostic Signature in the Role of ICI Treatment

In the immunotherapeutic cohort, patients with a low-risk score exhibited a significantly prolonged survival rate (non-platinum-treated dataset: p = 0.0014, Supplementary Figure S5A; platinum-treated dataset: p = 0.033, Supplementary Figure S6A). Patients without platinum treatment indicated marked clinical benefits from immunotherapy in the low-risk score group compared to those with a high-risk score (p = 0.0027, Supplementary Figure S5B), but not significantly in patients with platinum treatment (p = 0.44, Supplementary Figure S6B). Further analysis revealed that a higher tumor mutation load in patients with platinum treatment was significantly associated with a low-risk score (p = 0.039, Supplementary Figure S6C), but not in patients without platinum treatment (p = 0.21, Supplementary Figure S5C). Interestingly, higher neoantigen burden in patients without platinum treatment was significantly associated with a low-risk score (p = 0.025, Supplementary Figure S5D), but not in patients with platinum treatment (p = 0.5, Supplementary Figure S6D).

Discussion

The immune cells within the tumor microenvironment (TME) function play a key role in tumorigenesis (Lei et al., 2020). Growing evidence has revealed the therapeutic potential of ICIs in colon cancer (Kimura et al., 2020; Wang et al., 2020). However, the limited knowledge on the characteristics of the TME, to some extent, hindered the development of the application of immunotherapy. In recent years, many efforts have been made to identify immune-related biomarkers for the diagnosis and prognosis of colon cancer (Zhou et al., 2019; Laghi et al., 2020; Li et al., 2020). However, more reliable biomarkers still need to be explored to maximize the application of immunotherapy.

In this study, we developed a prognostic risk score model based on seven IRGs, named GUCA2A, CHGB, SSTR2, VIP, OXTR, IL1A, and GRP, which has been verified in multiple cohorts across different platforms. Among them, GUCA2A, VIP, and OXTR have been demonstrated to be significantly associated with the prognosis of colon cancer (Zhang et al., 2019; Houxi Xu et al., 2020; Kang Lin et al., 2020; Zhang et al., 2020). A previous study reported that guanylyl cyclase C (GUCY2C) could act as a tumor suppressor and play an important role in orchestrating intestinal homeostatic mechanisms, which could be used as a therapeutic target for colon cancer patients (such as the FDA-approved oral GUCY2C ligand linaclotide, Linzess™) (Pattison et al., 2016). GUCA2A may be considered as a potential marker for the prognosis and therapeutic target in colon cancer by binding and activating GUCY2C. As a precursor of regulatory peptide, the relationship between CHGB and tumor is not clear. However, CHGB was suggested to be an immune-related signature for low-grade glioma (Liu et al., 2021) and head and neck squamous cell carcinoma (Zhang et al., 2021). Previous studies also experimentally demonstrated that an abnormal expression of CHGB was associated with aggressive VHL-associated pancreatic neuroendocrine tumors (validated by immunohistochemistry) (Weisbrod et al., 2013), pancreatic cancer (validated by qPCR) (Jia-Sheng Xu et al., 2020), and small cell lung cancer (validated by immunoblotting and qPCR) (Moss et al., 2009). SSTR2, as a G protein-coupled cell surface receptor, can be activated by extracellular ligands, which leads to the inhibition of cell proliferation (Lechner et al., 2021). Precious studies demonstrated that SSTR2 might serve as a molecular target in the diagnosis and treatment of thyroid cancer (Thakur et al., 2021), small intestinal neuroendocrine tumor (Elf et al., 2021), and neuroendocrine tumors (Si et al., 2021). VIP can provide protection from apoptosis in tumorigenesis (Sastry et al., 2017). OXTR and its ligand oxytocin (OXT) are present in the gastrointestinal system and involved in tumorigenesis (Ma et al., 2019). IL1A was involved in various immune responses, inflammatory processes, and hematopoiesis, which might be associated with colon tumorigenesis (Yoshikawa et al., 2017). To our knowledge, the relationship between GRP and tumorigenesis has not been reported.

Furthermore, we systematically explored the characteristics of the tumor immune microenvironment. The results revealed that the tumor-infiltrating resting NK cells or M0 Macrophages, TP53 mutation rates, and TMB could be independent prognostic signatures for colon cancer. Additionally, we observed that the expression levels of checkpoint genes (PD-L1, PD-1, and CTLA-4) were higher in high-risk score patients, which may suggest that our immune-related risk score model was capable of providing support for immunotherapy. More importantly, the immune-related signature was also significantly associated with OS in patients with anti-PD-L1 treatment. We speculated that patients with a low-risk score might be more sensitive to ICI therapy based on the result of Supplementary Figures S5, S6.

In addition, we compared the performance of our constructed immune-related prognostic model with the published prognostic model of colon cancer based on the cohorts TCGA-COAD, GSE39582, and GSE17536, which is summarized in Supplementary Table S7. Our constructed prognostic model was relatively and effectively validated in more internal and external cohorts, including an immunotherapeutic cohort.

Conclusion

In summary, a robust immune-related prognostic model was constructed, and the characteristics of the tumor immune microenvironment were explored, which may be helpful for the prognosis and immunotherapy of colon cancer patients.

Data Availability Statement

Publicly available datasets were analyzed in this study. These data can be found in the following: the expression profiles and corresponding clinical information of TCGA-COAD were downloaded from https://xenabrowser.net/datapages/. GSE39582 and GSE17536 were downloaded from the NCBI-GEO database (https://www.ncbi.nlm.nih.gov/geo/); and the immunotherapeutic cohort (IMvigor210) was available according to the guideline on http://research-pub.gene.com/IMvigor210CoreBiologies using the IMvigor210CoreBiologies R package.

Author Contributions

LC and EC were involved in the study concept, design, and drafting of the manuscript. HZ and YB put forward some kind suggestions. TL and JY helped in analyzing and interpreting the data.

Funding

This work was supported by a grant from the Key Science and Technology Program of the Shaanxi Province (2019ZDLSF02-05) and the National Natural Scientific Foundation of China (81974378 and 82003115).

Conflict of Interest

HZ and YB were employed by the company Shenzhen Nucleus Gene Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.801484/full#supplementary-material

References

Cereda, M., Gambardella, G., Benedetti, L., Iannelli, F., Patel, D., Basso, G., et al. (2016). Patients with Genetically Heterogeneous Synchronous Colorectal Cancer Carry Rare Damaging Germline Mutations in Immune-Related Genes. Nat. Commun. 7, 12072. doi:10.1038/ncomms12072

PubMed Abstract | CrossRef Full Text | Google Scholar

Chalabi, M., Fanchi, L. F., Dijkstra, K. K., Van den Berg, J. G., Aalbers, A. G., Sikorska, K., et al. (2020). Neoadjuvant Immunotherapy Leads to Pathological Responses in MMR-Proficient and MMR-Deficient Early-Stage colon Cancers. Nat. Med. 26 (4), 566–576. doi:10.1038/s41591-020-0805-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Luo, J., and Guo, J. (2020). Development and Validation of a Five-Immune Gene Prognostic Risk Model in colon Cancer. BMC Cancer 20 (1), 395. doi:10.1186/s12885-020-06799-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Dienstmann, R., Vermeulen, L., Guinney, J., Kopetz, S., Tejpar, S., and Tabernero, J. (2017). Consensus Molecular Subtypes and the Evolution of Precision Medicine in Colorectal Cancer. Nat. Rev. Cancer 17 (2), 79–92. doi:10.1038/nrc.2016.126

PubMed Abstract | CrossRef Full Text | Google Scholar

Elf, A.-K., Johanson, V., Marin, I., Bergström, A., Nilsson, O., Svensson, J., et al. (2021). Evaluation of SSTR2 Expression in SI-NETs and Relation to Overall Survival after PRRT. Cancers 13 (9), 2035. doi:10.3390/cancers13092035

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghonim, M. A., Ibba, S. V., Tarhuni, A. F., Errami, Y., Luu, H. H., Dean, M. J., et al. (2021). Targeting PARP-1 with Metronomic Therapy Modulates MDSC Suppressive Function and Enhances Anti-PD-1 Immunotherapy in colon Cancer. J. Immunother. Cancer 9 (1), e001643. doi:10.1136/jitc-2020-001643

CrossRef Full Text | Google Scholar

Keum, N., and Giovannucci, E. (2019). Global burden of Colorectal Cancer: Emerging Trends, Risk Factors and Prevention Strategies. Nat. Rev. Gastroenterol. Hepatol. 16 (12), 713–732. doi:10.1038/s41575-019-0189-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimura, Y., Sumiyoshi, M., Kiyoi, T., and Baba, K. (2020). Dihydroxystilbenes Prevent Azoxymethane/dextran Sulfate Sodium-Induced colon Cancer by Inhibiting colon Cytokines, a Chemokine, and Programmed Cell Death-1 in C57BL/6J Mice. Eur. J. Pharmacol. 886, 173445. doi:10.1016/j.ejphar.2020.173445

CrossRef Full Text | Google Scholar

Laghi, L., Negri, F., Gaiani, F., Cavalleri, T., Grizzi, F., de’ Angelis, G. L., et al. (2020). Prognostic and Predictive Cross-Roads of Microsatellite Instability and Immune Response to Colon Cancer. Int. J. Mol. Sci. 21 (24), 9680. doi:10.3390/ijms21249680

CrossRef Full Text | Google Scholar

Lechner, M., Schartinger, V. H., Steele, C. D., Nei, W. L., Ooft, M. L., Schreiber, L.-M., et al. (2021). Somatostatin Receptor 2 Expression in Nasopharyngeal Cancer Is Induced by Epstein Barr Virus Infection: Impact on Prognosis, Imaging and Therapy. Nat. Commun. 12 (1), 117. doi:10.1038/s41467-020-20308-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, X., Lei, Y., Li, J.-K., Du, W.-X., Li, R.-G., Yang, J., et al. (2020). Immune Cells within the Tumor Microenvironment: Biological Functions and Roles in Cancer Immunotherapy. Cancer Lett. 470, 126–133. doi:10.1016/j.canlet.2019.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Wen, D., Li, X., Yao, C., Chong, W., and Chen, H. (2020). Identification of an Immune Signature Predicting Prognosis Risk and Lymphocyte Infiltration in Colon Cancer. Front. Immunol. 11, 1678. doi:10.3389/fimmu.2020.01678

PubMed Abstract | CrossRef Full Text | Google Scholar

Lichtenstern, C. R., Ngu, R. K., Shalapour, S., and Karin, M. (2020). Immunotherapy, Inflammation and Colorectal Cancer. Cells 9 (3), 618. doi:10.3390/cells9030618

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, K., Huang, J., Luo, H., Luo, C., Zhu, X., Bu, F., et al. (2020). Development of a Prognostic index and Screening of Potential Biomarkers Based on Immunogenomic Landscape Analysis of Colorectal Cancer. Aging 12 (7), 5832–5857. doi:10.18632/aging.102979

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Y., Pan, X., Chen, Z., Lin, S., and Chen, S. (2020). Identification of an Immune-Related Nine-lncRNA Signature Predictive of Overall Survival in Colon Cancer. Front. Genet. 11, 318. doi:10.3389/fgene.2020.00318

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W., Zou, J., Ren, R., Liu, J., Zhang, G., and Wang, M. (2021). A Novel 10-Gene Signature Predicts Poor Prognosis in Low Grade Glioma. Technol. Cancer Res. Treat. 20, 153303382199208. doi:10.1177/1533033821992084

CrossRef Full Text | Google Scholar

Ma, M., Li, L., Chen, H., and Feng, Y. (2019). Oxytocin Inhibition of Metastatic Colorectal Cancer by Suppressing the Expression of Fibroblast Activation Protein-α. Front. Neurosci. 13, 1317. doi:10.3389/fnins.2019.01317

PubMed Abstract | CrossRef Full Text | Google Scholar

Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayakonda, A., Lin, D.-C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Moss, A. C., Jacobson, G. M., Walker, L. E., Blake, N. W., Marshall, E., and Coulson, J. M. (2009). SCG3 Transcript in Peripheral Blood Is a Prognostic Biomarker for REST-Deficient Small Cell Lung Cancer. Clin. Cancer Res. 15 (1), 274–283. doi:10.1158/1078-0432.CCR-08-1163

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12 (5), 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Pattison, A. M., Merlino, D. J., Blomain, E. S., and Waldman, S. A. (2016). Guanylyl Cyclase C Signaling axis and colon Cancer Prevention. World J. Gastroenterol. 22 (36), 8070–8077. doi:10.3748/wjg.v22.i36.8070

CrossRef Full Text | Google Scholar

Riley, R. S., June, C. H., Langer, R., and Mitchell, M. J. (2019). Delivery Technologies for Cancer Immunotherapy. Nat. Rev. Drug Discov. 18 (3), 175–196. doi:10.1038/s41573-018-0006-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sastry, K. S., Chouchane, A. I., Wang, E., Kulik, G., Marincola, F. M., and Chouchane, L. (2017). Cytoprotective Effect of Neuropeptides on Cancer Stem Cells: Vasoactive Intestinal Peptide-Induced Antiapoptotic Signaling. Cell Death Dis 8 (6), e2844. doi:10.1038/cddis.2017.226

PubMed Abstract | CrossRef Full Text | Google Scholar

Si, Y., Kim, S., Ou, J., Lu, Y., Ernst, P., Chen, K., et al. (2021). Anti-SSTR2 Antibody-Drug Conjugate for Neuroendocrine Tumor Therapy. Cancer Gene Ther. 28 (7-8), 799–812. doi:10.1038/s41417-020-0196-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2018). Cancer Statistics, 2018. CA: A Cancer J. Clinicians 68 (1), 7–30. doi:10.3322/caac.21442

CrossRef Full Text | Google Scholar

Thakur, S., Daley, B., Millo, C., Cochran, C., Jacobson, O., Lu, H., et al. (2021). 177Lu-DOTA-EB-TATE, a Radiolabeled Analogue of Somatostatin Receptor Type 2, for the Imaging and Treatment of Thyroid Cancer. Clin. Cancer Res. 27 (5), 1399–1409. doi:10.1158/1078-0432.CCR-20-3453

PubMed Abstract | CrossRef Full Text | Google Scholar

Tolba, M. F. (2020). Revolutionizing the Landscape of Colorectal Cancer Treatment: The Potential Role of Immune Checkpoint Inhibitors. Int. J. Cancer 147 (11), 2996–3006. doi:10.1002/ijc.33056

CrossRef Full Text | Google Scholar

Wang, H., Wang, X., Xu, L., Zhang, J., and Cao, H. (2019). A Molecular Sub-cluster of colon Cancer Cells with Low VDR Expression Is Sensitive to Chemotherapy, BRAF Inhibitors and PI3K-mTOR Inhibitors Treatment. Aging 11 (19), 8587–8603. doi:10.18632/aging.102349

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Bhave, M. S., Yagita, H., and Cardell, S. L. (2020). Natural Killer T-Cell Agonist α-Galactosylceramide and PD-1 Blockade Synergize to Reduce Tumor Development in a Preclinical Model of Colon Cancer. Front. Immunol. 11, 581301. doi:10.3389/fimmu.2020.581301

PubMed Abstract | CrossRef Full Text | Google Scholar

Weisbrod, A. B., Zhang, L., Jain, M., Barak, S., Quezado, M. M., and Kebebew, E. (2013). Altered PTEN, ATRX, CHGA, CHGB, and TP53 Expression Are Associated with Aggressive VHL-Associated Pancreatic Neuroendocrine Tumors. Horm. Canc 4 (3), 165–175. doi:10.1007/s12672-013-0134-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, H., Ma, Y., Zhang, J., Gu, J., Jing, X., Lu, S., et al. (2020). Identification and Verification of Core Genes in Colorectal Cancer. Biomed. Res. Int. 2020, 1–13. doi:10.1155/2020/8082697

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J. S., Liao, K. L., Wang, X., He, J., and Wang, X. Z. (2020). Combining Bioinformatics Techniques to Explore the Molecular Mechanisms Involved in Pancreatic Cancer Metastasis and Prognosis. J. Cel. Mol. Med. 24 (24), 14128–14138. doi:10.1111/jcmm.16023

CrossRef Full Text | Google Scholar

Yoshikawa, T., Wu, J., Otsuka, M., Kishikawa, T., Suzuki, N., Takata, A., et al. (2017). Repression of MicroRNA Function Mediates Inflammation-Associated Colon Tumorigenesis. Gastroenterology 152 (3), 631–643. doi:10.1053/j.gastro.2016.10.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

Yu, G., Wu, Y., Wang, W., Xu, J., Lv, X., Cao, X., et al. (2019). Low-dose Decitabine Enhances the Effect of PD-1 Blockade in Colorectal Cancer with Microsatellite Stability by Re-modulating the Tumor Microenvironment. Cell Mol Immunol 16 (4), 401–409. doi:10.1038/s41423-018-0026-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Du, Y., Wang, Z., Lou, R., Wu, J., and Feng, J. (2019). Integrated Analysis of Oncogenic Networks in Colorectal Cancer Identifies GUCA2A as a Molecular Marker. Biochem. Res. Int. 2019, 1–13. doi:10.1155/2019/6469420

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Qin, C., Gan, H., Guo, X., and Zhang, L. (2020). Construction of an Immunogenomic Risk Score for Prognostication in Colon Cancer. Front. Genet. 11, 499. doi:10.3389/fgene.2020.00499

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Chen, P., Zhou, Q., Wang, H., Hua, Q., Wang, J., et al. (2021). A Novel Immune-Related Prognostic Signature in Head and Neck Squamous Cell Carcinoma. Front. Genet. 12, 570336. doi:10.3389/fgene.2021.570336

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, R., Zhang, J., Zeng, D., Sun, H., Rong, X., Shi, M., et al. (2019). Immune Cell Infiltration as a Biomarker for the Diagnosis and Prognosis of Stage I-III colon Cancer. Cancer Immunol. Immunother. 68 (3), 433–442. doi:10.1007/s00262-018-2289-7

PubMed Abstract | CrossRef Full Text | Google Scholar