We systematically searched publicly available colon cancer datasets. Studies with no survival or TNM staging information were removed from further assessment. Five cohorts with bulk sequencing (TCGA-COAD, GSE39582, GSE37892, GSE17538, GSE14333) and one single-cell sequencing dataset (GSE178341) were enrolled. The RNA sequencing and clinical data of GSE39582, GSE37892, GSE17538, GSE14333 and GSE178341 were downloaded from the medics Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Gene expression data and corresponding clinical information from The Cancer Genome Atlas (TCGA) colon cancer project were downloaded from the UCSC Xena browser (https://xenabrowser.net/datapages/).

Mutation data of the TCGA-COAD cohort were downloaded from the National Cancer Institute Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/mc3-2017). Only TNM stage III CC patients in each dataset were included in this study. The above datasets were utilized in compliance with the ethical requirements of the GEO and TCGA projects. The study was conducted in accordance with the Declaration of Helsinki.

We performed TME immune scoing for stage III CC patients in the GSE39582 and TCGA-COAD cohorts using R package “TMEscore”. The cut-off value of the TME immune score was selected by X-tile software (version 3.6.1, https://medicine.yale.edu/lab/rimm/research/software/). Based on the cut-off values, patients in both cohorts were divided into high and low TME immune score groups. We compared the difference in relapse-free survival (RFS) between the two groups using the Kaplan-Meier method. The receiver operating characteristic (ROC) curves were then plotted to assess the predictive power of the TME immune score for RFS. Genes with |log2 Fold Change| > 1 and adjusted P value < 0.05 were defined as differentially expressed genes (DEGs). The intersection of the DEGs from the GSE39582 and TCGA-COAD datasets was entered into the STRING database (https://cn.string-db.org/) to construct the protein-protein interaction (PPI) network. The network was then imported into Cytoscape software (version 3.9.0, https://cytoscape.org/), and key gene modules were identified using the Molecular Complex Detection (MCODE) algorithm. Functional and pathway enrichment analysis of key gene modules were performed in the Metascape database (https://metascape.org/gp/index.html#/main/step1).

We then used univariate Cox regression analysis to identify genes significantly associated with RFS in the above key gene modules. The association between each gene and RFS was further assessed by meta-analysis combining the datasets GSE39582 (n = 206), TCGA-COAD (n = 128), GSE14333 (n = 81), GSE17538 (n = 75) and GSE37892 (n = 57). If robust heterogeneity was not observed (I² < 40%, P > 0.05), the fixed-effects model was chosen to pool HRs from different cohorts. Otherwise, the random-effects model was selected.

In the TCGA-COAD dataset, 113 patients with stage III CC had complete somatic mutation data. In contrast, only TP53, KRAS, BRAF and mismatch repair mutations were available in the GSE39582 dataset. We compared the mutations and the tumor mutation burden (TMB) between the high and low EPSTI1 groups in the TCGA-COAD cohort. We also analyzed the mutation status of TP53, KRAS, and BRAF and microsatellite stability in both groups.

To quantify the degree of immune cell infiltration in each stage III CC sample, we applied three widely accepted algorithms for evaluation: TIMER 2.0 (http://timer.cistrome.org/), TIDE (http://tide.dfci.harvard.edu/) and CIBERSORTx (https://cibersortx.stanford.edu). According to the algorithm instructions, we uploaded the prepared gene expression matrix into the web tool to obtain infiltration scores. Only the immune cell types detected in more than 50% of the samples were included in further analysis. In the CIBERSORTx estimation procedure, we ran the web tool with LM22 gene signatures and 1000 permutations.

We performed scRNA-seq analysis using the R package “Seurat” (version 4.1.0). Cells with less than 200 genes and more than 50% mitochondrial counts were excluded from the analysis. The expression matrix was then normalized using the “SCTransform” function, and the top 3000 highly variable genes were subjected to principal component analysis (PCA). We constructed the shared nearest neighbor (SNN) graph and the uniform manifold approximation and projection (UMAP) embedding with the top 20 principal components. The identification of main cell types was consistent with the original literature. Based on the M1/M2 macrophage and pDC/mDC gene signatures summarized in the literature (see Table S3 ), we further subclustered the macrophages and DCs using the R package CelliD (16). The proportion of cells in each subpopulation was then calculated. The expression of M1 versus M2 up- and down-regulated genes (17), and plasmacytoid cell type DC (pDC) versus mDC up- and down-regulated genes (18) were scored using the “AddModuleScore” function. The same function was also used to compute the activity scores of immune-related signaling pathways from the Broad Institute’s Hallmark collection.

Patient groups or cell groups were compared using Welch’s t-test if the continuous variables were normally distributed; otherwise, the Mann-Whitney U test was applied. Categorical variables were compared using the chi-square test. We plotted survival curves using the Kaplan-Meier method and used the log-rank test to compare survival differences. The predictive validity of the model was quantified by the area under the ROC curve (AUC). For correlation analysis, we calculated Pearson or Spearman correlation coefficients as indicated. A P value < 0.05 was considered significant. All statistical analyses were performed using R software (version 4.1.0).

To investigate the relationship between the tumor immune microenvironment and RFS in stage III CC patients, we applied the TMEscore model to perform immune scoring in the GSE39582 (n = 206) and TCGA-COAD (n = 128) cohorts. The cut-off value of the TME immune score in each cohort was determined by X-tile software, and the patients were then divided into high- and low-immune score groups ( Table S1 ). According to the cut-off value, 103 patients in the GSE38582 cohort were assigned to the high immune score group, and the remaining 103 patients were assigned to the low immune score group. In the TCGA-COAD cohort, the numbers were 82 and 46, respectively. The high and low TME score groups had different distribution features on the PCA dimensionality reduction map, reflecting the difference in the expression of immune-related genes between the two groups ( Figures 1 A, D ). In both cohorts, the Kaplan-Meier survival curves showed that RFS was significantly worse in the low-immune score group ( Figures 1 B, E ; log-rank P = 0.006 in GSE39582, log-rank P = 0.04 in TCGA-COAD). To measure the predictive performance of the TMEscore model, we calculated the time-dependent AUC values for both cohorts at 2, 5 and 7 years. The AUCs at these time points were 0.60, 0.58 and 0.61 in the GSE39582 cohort, while they were 0.62, 0.61 and 0.56 in the TCGA-COAD cohort ( Figures 1C, F ). These results suggest that the TME immune score can predict RFS in patients with stage III CC.

We then analyzed the DEGs between the different TME immune scoring groups in the GSE39582 and TCGA-COAD datasets. Compared to the low TME immune score group, patients in the high TME immune score group had 116 genes that were significantly upregulated and 7 genes that were significantly downregulated in GSE39582 ( Figures 2A , Table S2 ), while 182 genes were upregulated and 3 genes were downregulated in the high immune score group in TCGA-COAD ( Figure 2B ). We then intersected the DEGs from the two cohorts to obtain 53 common genes ( Figure 2C ). By entering these genes into the STRING database, we constructed a protein-protein interaction network ( Figure 2D ). The minimum required interaction score was set as medium confidence (0.400). The PPI network was then imported into Cytoscape software, and two key gene modules were identified by the MCODE algorithm: MCODE cluster 1 contained 24 genes ( Figure 2E ), while MCODE cluster 2 contained 16 genes ( Figure 2F ). The enrichment analysis suggests that these genes are mainly related to human immune function, particularly the interferon gamma signaling pathway ( Figure 2G, H ).

Survival analysis was performed via univariate Cox regression to identify a subset of genes closely related to RFS in the MCODE clusters. The results revealed that only EPSTI1 and CXCL11 were significantly associated with RFS in both the GSE39582 and TCGA-COD cohorts ( Figure 3 ). In GSE39582, the hazard ratios (HRs) for EPSTI1 and CXCL11 were 0.79 (95% CI, 0.65-0.96; P < 0.05) and 0.89 (95% CI, 0.80-0.99; P < 0.05), respectively. In TCGA-COAD, the HRs for the two genes were 0.69 (95% CI, 0.50-0.95; P < 0.05) and 0.74 (95% CI, 0.57-0.97; P < 0.05), respectively. Using X-tile software, we determined the cut-off values for EPSTI1 and CXCL11 expression levels, respectively, and divided patients in the GSE39582 and TCGA-COAD cohorts into groups with high and low expression levels of the corresponding genes. Kaplan-Meier survival curves showed that the high EPSTI1 group had prolonged RFS ( Figure 4A, B ; P < 0.05), and the high CXCL11 group also had a better prognosis (Supplementary Figure 1A, B ; P < 0.05). To further confirm the relationship between these two genes and RFS in patients with stage III CC, 547 patients from five datasets, GSE39582 (n = 206), TCGA-COAD (n = 128), GSE14333 (n = 81), GSE17538 (n = 75) and GSE37892 (n = 57), were subjected to meta-analysis. No significant heterogeneity was observed among these datasets (I ² < 40%, P > 0.05), and thus, the fixed effects model was selected for the meta-analysis. The pooled HRs of EPSTI1 and CXCL11 were 0.81 (95% CI, 0.71-0.91) and 0.92 (95% CI, 0.86-0.98), respectively ( Figure 4C ; Supplementary Figure 1C ). These results suggest that high levels of EPSTI1 and CXCL11 expression in stage III CC are significantly associated with prolonged RFS.

To date, a wide range of studies have explored the role of CXCL11 in antitumor immunity in diverse tumors, including CC (19, 20). However, the role of EPSTI1 in antitumor immunity in CC has not yet been reported. Therefore, in the following analysis, we focused on the association of EPSTI1 with the immune microenvironment of CC and its impact on RFS.

The top 20 mutated genes in the high and low EPSTI1 groups from the TCGA-COAD are illustrated in Supplementary Figure 2A . Although the mean TMB was greater in the high EPSTI1 group, there was no significant difference between the two groups ( Supplementary Figure 2B ; P > 0.05). The mutation rate of TP53 exceeded 50% in both groups, but no significant difference was observed ( Supplementary Figures 2C, D ). The proportions of patients with BRAF mutation or microsatellite instability were also not significantly different between the two groups ( Supplementary Figures 2C, D ). In the GSE39582 cohort, the KRAS mutation rate in the high EPSTI1 group was 37%, which was significantly lower than the 61% rate in the low EPSTI group ( Supplementary Figure 2D ). Although the KRAS mutation rate in the high EPSTI1 group was still smaller than that in the low EPSTI group in the TCGA-COAD cohort (33% vs. 41%), the difference was not statistically significant ( Supplementary Figure 2C ).

In this study, we evaluated immune infiltration in tumor tissue based on bulk RNA-seq data using three different approaches, including TIMER 2.0, TIDE and CIBERSORTx. TIMER 2.0 analysis demonstrated that macrophages, myeloid dendritic cells, neutrophils and CD8 T cells were more abundantly infiltrated in the high EPSTI1 group ( Figure 5A ; P < 0.05). TIDE found significantly lower infiltration abundance of M2 tumor-associated macrophages (TAM) in the high EPSTI1 group in both cohorts when assessing immunosuppressive cells ( Figure 5B ; P < 0.05). For MDSCs, the infiltration abundance was significantly lower in the high EPSTI1 group in the TCGA-COAD cohort ( Figure 5B ; P < 0.05). In the GSE39582 cohort, the difference was only marginally significant ( Figure 5B ; P = 0.076), although the high EPSTI1 group still had a smaller infiltration abundance than the low EPSTI1 group. The CIBERSORTx algorithm calculates the infiltration status of 22 types of immune cells. In this study, we excluded a cell type if its infiltrative abundance was calculated to be zero in more than 50% of the samples. Pearson correlation analysis revealed that in the GSE39582 cohort, EPSTI1 expression was positively correlated with the infiltration of memory activated CD4 T cells, follicular helper T cells, and M1 macrophages ( Figure 5C ; R > 0.2, P < 0.05), and negatively correlated with the infiltration of resting memory CD4 T cells, resting NK cells and activated mast cells ( Figure 5C ; R < -0.2, P < 0.05). In the TCGA-COAD cohort, EPSTI1 expression was positively correlated with the infiltration of CD8 T cells, follicular helper T cells and M1 macrophages ( Figure 5D ; R > 0.2, P < 0.05), and negatively correlated with the infiltration of naive B cells and resting NK cells ( Figure 5D ; R < -0.2, P < 0.05). Notably, M1 macrophages displayed a strong correlation with EPSTI1 expression in both cohorts (R = 0.561 in GSE39582; R = 0.394 in TCGA-COAD).

In the bulk RNA-seq datasets above, EPSTI1 was closely associated with immune infiltrates. To further investigate the relationship between EPSTI1 and tumor immune cells, we analyzed 32 stage III CC patients from the single-cell dataset GSE178341. According to the quality control criteria in the original literature, we obtained 114,928 cells from CC tissues. Cell types were identified with reference to the original paper. From the EPSTI1 dimplot and dotplot ( Figure 6A ), we found that EPSTI1 was predominantly expressed in immune cells, both in terms of EPSTI1 expression levels and the proportion of EPST1-positive cells. Further sorting of immune cells revealed that macrophages and DCs were the main cell clusters expressing EPSTI1 ( Figure 6B ), similar to the immune infiltration analysis of the bulk RNA-seq datasets above.

Several studies have indicated that M1macrophages exert tumor-preventing activities, whereas M2 macrophages are associated with immunosuppression (21, 22, 23). At the single-cell level, the ratio of EPSTI1+ macrophages was significantly correlated with that of M1 macrophages (R = 0.405, P = 0.022; Figure 6C , Supplementary Figure 3A ), which was consistent with the CIBERSORTx analysis of bulk sequencing described previously. Moreover, EPSTI1⁺ macrophages scored significantly lower on the M2 but significantly higher on the M1 signature modules ( Figures 6D, E ; P < 0.05). In the scoring of cell signaling pathway activity, EPSTI1⁺ macrophages had a higher mean score for immune-related signaling pathways ( Figure 6F ), showing a different function pattern from EPSTI1^- macrophages.

DCs were another major cell cluster that expressed EPSTI1 in our study. By origin, DCs can be classified into mDCs, which are derived from common myeloid progenitors (CMPs) that also produce macrophages, and pDCs, which are derived from common lymphocyte progenitors (CLPs) that also produce B cells, T cells and NK cells. Although the proportion of EPSTI1+ DCs was positively correlated with that of mDCs without statistical significance (R = 0.302, P = 0.093; Supplementary Figure 3B, C ), EPSTI1⁺ DCs had significantly lower pDC but higher mDC feature scores ( Supplementary Figures 3D, E ; P < 0.05). Similar to EPSTI1⁺ macrophages, the mean scores of immune-related signaling pathways were higher in EPSTI1+ DCs ( Figure 6G ). The association between DCs and macrophages in EPSTI1 expression was also inspected by Pearson correlation analysis. In the 32 patients with stage III CC in GSE178341, macrophages and DCs were significantly correlated (R > 0.6, P < 0.05), both in terms of the average level of EPSTI1 expression and the proportion of EPSTI1⁺ cells ( Figures 6H, I ).