Identification of the molecular characteristics associated with microsatellite status of colorectal cancer patients for the clinical application of immunotherapy

Background: Mismatch repair-proficient (pMMR) microsatellite stability (MSS) in colorectal cancer (CRC) indicates an unfavorable therapeutic response to immunotherapy with immune checkpoint inhibitors (ICIs). However, the molecular characteristics of CRC patients with pMMR MSS remain largely unknown.

Methods: Heterogeneities between mismatch repair-deficient (dMMR) microsatellite instability (MSI) and pMMR MSS CRC patients were investigated at the single-cell level. Next, an MSS-related risk score was constructed by single-sample gene set enrichment analysis (ssGSEA). The differences in immune and functional characteristics between the high- and low-score groups were systematically analyzed.

Results: Based on the single-cell RNA (scRNA) atlas, an MSS-specific cancer cell subpopulation was identified. By taking the intersection of the significant differentially expressed genes (DEGs) between different cancer cell subtypes of the single-cell training and validation cohorts, 29 MSS-specific cancer cell marker genes were screened out for the construction of the MSS-related risk score. This risk score signature could efficiently separate pMMR MSS CRC patients into two subtypes with significantly different immune characteristics. The interactions among the different cell types were stronger in the MSS group than in the MSI group, especially for the outgoing signals of the cancer cells. In addition, functional differences between the high- and low-score groups were preliminarily investigated.

Conclusion: In this study, we constructed an effective risk model to classify pMMR MSS CRC patients into two completely different groups based on the specific genes identified by single-cell analysis to identify potential CRC patients sensitive to immunotherapy and screen effective synergistic targets.

Introduction

Colorectal cancer (CRC) is one of the leading causes of cancer-related death worldwide(Sung et al., 2021; Bray and Parkin, 2022). The treatments of CRC have made great progress in recent years, such as targeted therapy based on gene detection and tumor heterogeneity(Adamopoulos et al., 2021; Kotani et al., 2022). However, the prognosis of CRC patients with advance-stage disease, recurrence and metastasis is still unsatisfactory(André et al., 2022; Ciardiello et al., 2022; Topham et al., 2022).

Immunotherapy represented by immune checkpoint inhibitors (ICIs) has shown good curative effects in many malignant tumors, but its application in CRC is still limited to patients with mismatch repair-deficient (dMMR) microsatellite instability (MSI) status(Jin and Sinicrope, 2022; Weng et al., 2022). The lack of immune infiltration and low tumor mutation burden (TMB) are the main reasons for the insensitivity of CRC patients with mismatch repair-proficient (pMMR) microsatellite stability (MSS) status to immunotherapy(Ganesh et al., 2019). However, the genomic features of approximately 85% of CRC patients are pMMR MSS(Boland and Goel, 2010). Thus, it is urgent to clarify the molecular characteristics of pMMR MSS CRC patients to find effective synergistic targets for immunotherapy.

The differences in immune landscapes between dMMR MSI and pMMR MSS CRC patients have been preliminarily explored at the single-cell RNA (scRNA) level(Bao et al., 2020). However, the overall differences between the two groups have not been fully elucidated, and whether the specific molecular features of pMMR MSS patients are associated with insensitivity to immunotherapy remains largely unknown. Several studies have divided pMMR MSS CRC patients into “hot” and “cold” subtypes according to the abundances of immune cell infiltration based on single-cell RNA (scRNA) and bulk RNA sequencing data, and potential synergistic targets of immunotherapy for “hot” tumors with high abundances of immune cells have also been identified(Picco et al., 2021; Ren et al., 2022; Wang et al., 2022). However, there is still a lack of a single-cell level risk score signature to classify pMMR MSS patients into different subgroups and thus judge their sensitivity to immunotherapy.

In this study, the heterogeneities between dMMR MSI and pMMR MSS CRC patients were systematically elucidated by the scRNA atlas. A specific cancer cell cluster was identified in MSS CRC patients. Furthermore, a specific MSS-related gene list was assembled by integrated analysis of the training and validation scRNA cohorts, which were subsequently used for constructing the MSS-related risk score. This MSS-related risk score could efficiently separate CRC patients into two subgroups with completely different immune characteristics in the external transcriptome cohorts. CRC patients in the low-score group were identified as “hot” tumors with high abundances of immune cells and an inhibitory immune microenvironment, which means they are potentially sensitive to immunotherapy.

Materials and methods

Data acquisition

Five CRC patients (SMC03-T, SMC06-T, SMC10-T, SMC24-T and KUL01-T) with MSI-H status (MSI group) and seven CRC patients (SMC23-T, SMC25-T, KUL19-T, KUL21-T, KUL28-T, KUL30-T and KUL31-T) with MSS status (MSS group) were used as a single-cell training cohort (GSE132465 and GSE144735)(Lee et al., 2020), while two CRC patients (T_ cac7 and T_ cac12) with MSI-H status and four CRC patients (T_cac1, T_cac2, T_cac4 and T_cac10) with MSS status were used as a single-cell validation cohort in this study (GSE200997)(Khaliq et al., 2022). In detail, the training cohort comprised the MSI group with 9366 high-quality single cells, and the MSS group with 9603 high-quality single cells. The validation cohort comprised the MSI group with 5555 cells, and the MSS group with 5665 cells. All of the scRNA data were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), and the clinical information of these CRC samples is shown in Supplementary Table S1. In addition, transcriptome data in transcripts per million (TPM) of 601 CRC patients (185 patients with MSI status and 416 patients with MSS status) were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) as a transcriptome training cohort (TCGA-CRC cohort), while the microarray data of the other 527 CRC patients (73 patients with dMMR status and 454 patients with pMMR status) were downloaded from the GEO database (GSE39582) as a transcriptome validation cohort (GEO-CRC cohort)(Marisa et al., 2013).

ScRNA data analysis

All of the scRNA data were analyzed by the “Seurat” R package in this study(Satija et al., 2015). Specifically, data normalization, screening of highly variable genes (HVGs), removal of interbatch differences and data integration were conducted by the “SCTransform” function. Subsequently, dimension reduction and cluster analyses were conducted by the uniform manifold approximation and projection (UMAP) and “FindClusters” methods respectively. Significantly highly expressed genes in each cell cluster were identified by the “FindAllMarkers” method. After that, all of the cell clusters were annotated to specific cell types by the specific markers as previously reported (Supplementary Table S2)(Lee et al., 2020; Liu et al., 2022). Moreover, significant differentially expressed genes (DEGs) between different cell clusters were identified by the “FindMarkers” method with a min.pct of 0.5 and logfc.threshold of 2. Finally, cell-cell communication analysis among different cell types was conducted by the “CellChat” R package(Jin et al., 2021).

Calculation of the risk scores of CRC samples

The risk score of each CRC sample was calculated by the single-sample gene set enrichment analysis (ssGSEA) algorithm based on transcriptome data according to the specific single-cell marker genes. CRC patients in the transcriptome training and validation cohorts were divided into high- and low-score groups according to the median value of the risk score. Survival analysis between the two groups was conducted by the “survMisc” and “survminer” R packages.

Immune analysis and gene set variation analysis (GSVA)

Abundances of 28 immune cells in each CRC sample were calculated by the ssGSEA method according to the specific markers as previously reported (Supplementary Table S3). Correlations between the abundances of immune cells and the expression levels of single-cell marker genes were calculated by the Spearman correlation analysis method. In addition, functional differences between the high- and low-score groups in the transcriptome training and validation cohorts were identified by GSVA(Hänzelmann et al., 2013).

Drug prediction

CellMiner is a web tool that predicts the sensitivity of target genes to drugs based on the NCI-60 cell line set(Reinhold et al., 2019). In this study, potential drugs that have been approved by the Food and Drug Administration (FDA) targeting these single-cell marker genes were screened by CellMiner.

Statistical analysis

All of the data were analyzed by R 4.1.0 in this study. Continuous quantitative data between two groups were compared by Student’s t-test or the Wilcoxon test. Correlation analyses were conducted by the Spearman method. The survival data were analyzed by the log-rank test. A p-value <0.05 was considered statistically significant.

Results

Single-cell atlas of CRC patients with MSI and MSS status

A total of 9366 single cells from CRC patients with MSI status and 9603 cells with MSS status were used to construct the single-cell atlas as a training cohort. The expression profiles of the CRC patients with MSI and MSS status, as well as the correlation between the numbers of features and counts of the single cell, are visualized in Supplementary Figures S1A, B. After data normalization, data integration, dimensionality reduction and clustering analyses, the CRC cells were stratified into 21 clusters (Figure 1A and Supplementary Figure S1C), and the top three significantly highly expressed genes in each cluster are shown in Figure 1B. The general differences between the MSI and MSS groups are visualized in Supplementary Figure S1D, while there were no significant batch effects among the different GEO sources, patients and cell cycles (Supplementary Figures S1E–G). The absolute amounts and percentages of these 21 cell clusters from CRC patients with MSI or MSS status are visualized in Figures 1C, D. We note that Clusters 4, 15 and 19 account for a significantly higher proportion in the MSS group (exceeding 75%), while Clusters 3, 13, 17 and 18 account for a significantly higher proportion in the MSI group. To address the functional characteristics of these cell clusters, the cell types were annotated. These 21 cell clusters could be annotated to 9 cell types (B cells, cancer cells, CD4 T cells, CD8 T cells, endothelial cells, fibroblasts, mast cells, myeloid cells and other T cells) according to specific cell markers (Figures 1E, F), and the marker genes for each cell type were visualized by violin plots (Figure 1G). It is noteworthy that Clusters 3, 4, 15, 17 and 18 are all cancer cells, while Cluster 13 is CD8 T cells and Cluster 19 is mast cells. In addition, the proportion of no known cell types was more than 75% in the MSS group, except for mast cells (Supplementary Figures S2A, C).

FIGURE 1

FIGURE 1. Single-cell atlas of MSI and MSS CRC patients in the scRNA training cohort. (A) UMAP and cluster analyses (Clusters 0–20) of 9366 single cells from CRC patients with MSI status and 9603 cells with MSS status. (B) The heatmap visualizes the top three significantly highly expressed genes in each cell cluster. (C) Cell numbers of MSI and MSS CRC-derived cells in these 21 cell clusters. (D) Cell proportions of MSI and MSS CRC-derived cells in these 21 cell clusters. (E) Dot plots visualize the specific cell marker genes in each cell cluster. (F) UMAP plot of cell annotation results. (G) Violin plots of the specific cell marker genes.

It has been reported that the enrichment of CD8 T cells in CRC patients with MSI status is an important reason for their relatively higher sensitivity to immunotherapy with ICIs(Sui et al., 2021; Zou et al., 2021). Our analysis results further confirmed this knowledge at the single-cell level. Considering that mast cells only account for a very small proportion of total cells (80/18,969), we mainly focused on cancer cells in this study for subsequent analysis. To further investigate the specific characteristics of the CRC patients with MSS status, DEGs between MSS-specific cancer cell Clusters 4 and 15 and other cancer cell clusters were identified and 153 significant DEGs were screened out (Supplementary Table S4).

Cell-cell communication analysis

To investigate the differences in the integrated role among different cell types between the MSI and MSS groups, cell-cell communication analysis was conducted separately. The number and strength of the cell-cell interactions in the MSI group are shown in Figures 2A, B, and the specific intracellular signal information of the MSI group is shown in Figure 2C. We note that the interactions among different cell types were generally weak in the MSI group, and a relatively strong interaction only existed between myeloid cells and endothelial cells (Figures 2B, C). However, the interactions among different cell types were significantly strengthened in the MSS group (Figures 2D, E), and the outgoing signals of various cell types were significantly enhanced, especially for cancer cells (Figure 2F). Moreover, the detailed intercellular signal information of the MSI and MSS groups is visualized in Supplementary Figures S4, S5. These results, combined with the results of the single-cell atlas, revealed the heterogeneities between the MSS and MSI groups, and the potential core roles of specific cancer cell subtypes in MSS CRC patients were also preliminarily elucidated.

FIGURE 2

FIGURE 2. Visualization of cell-cell communication results in the scRNA training cohort. (A) Number of interactions among the nine cell types in the MSI group. (B) Interaction strength among the 9 cell types in the MSI group. (C) The potential incoming and outgoing signaling pathways among the 9 cell types in the MSI group. (D) Number of interactions among the 9 cell types in the MSS group. (E) Interaction strength among the 9 cell types in the MSS group. (F) The potential incoming and outgoing signaling pathways among the 9 cell types in the MSS group.

External validation of the specific cancer cell subtypes in CRC patients with MSS status

To validate the specific cancer cell subtypes in the MSS group, a single-cell validation cohort including 11,220 cells from CRC patients (5555 cells in the MSI group and 5665 cells in the MSS group) was analyzed. The expression profiles of the CRC cells, as well as the correlation between the numbers of features and counts, are visualized in Supplementary Figures S3A, B. After processing with the same method as in the training cohort, the CRC cells were divided into 20 clusters (Figure 3A and Supplementary Figure S3C), and the top three significantly highly expressed genes in each cluster are shown in Figure 3B. The general differences between the MSI and MSS groups are visualized in Supplementary Figure S3D, while there were no significant batch effects among different patients and cell cycles (Supplementary Figures S3E, F). The absolute amounts and percentages of these 20 cell clusters from CRC patients with MSI and MSS status are visualized in Figures 3C, D. It can be observed that Clusters 12, 13, 14 and 16 account for a significantly higher proportion in the MSS group (exceeding 75%), while Cluster 18 accounts for a significantly higher proportion in the MSI group. After cell annotation, 9 cell types were identified according to their specific cell markers (Figures 3E, F), and the marker genes for each cell type were visualized by violin plots (Figure 3G). For these cell clusters with significant differences, Cluster 12 was annotated as endothelial cells, Cluster 13 was annotated as fibroblasts, Cluster 14 was annotated as cancer cells, Cluster 16 was annotated as myeloid cells and Cluster 18 was annotated as CD4 T cells. These results at least further confirmed the existence of a specific cancer cell subpopulation in the MSS group. Moreover, the proportions of the 9 cell types in the MSI and MSS groups are visualized in Supplementary Figures S2D–F.

FIGURE 3

FIGURE 3. Single-cell atlas of MSI and MSS CRC patients in the scRNA validation cohort. (A) UMAP and cluster analyses (Clusters 0–19) of 5555 single cells from CRC patients with MSI status and 5665 cells with MSS status. (B) The heatmap visualizes the top three significantly highly expressed genes in each cell cluster. (C) Cell numbers of MSI and MSS CRC-derived cells in these 20 cell clusters. (D) Cell proportions of MSI and MSS CRC-derived cells in these 20 cell clusters. (E) Dot plots visualize the specific cell marker genes in each cell cluster. (F) UMAP plot of cell annotation results. (G) Violin plots of the specific cell marker genes.

Next, DEGs between specific cancer cell Cluster 14 and other cancer cell clusters were identified (Supplementary Table S5). By taking the intersection of the DEGs between different cancer cell subtypes of the single-cell training and validation cohorts, 29 single-cell marker genes were identified in the MSS-specific cancer cell subtype (Supplementary Table S6).

Calculation of the risk scores of CRC samples and immune analysis

First, the expression status of the 29 single-cell marker genes was analyzed in the TCGA-CRC cohort. The results showed that most of these genes were significantly differentially expressed between the MSS and MSI groups (Supplementary Figure S6A). Next, the risk score of each CRC sample with MSS status was calculated by the ssGSEA algorithm based on transcriptome data of these 29 specific single-cell marker genes. Then, the CRC patients were divided into high- and low-score groups according to the median risk score, and most of the single-cell marker genes were also significantly differentially expressed between the two risk groups (Supplementary Figure S6B). The results of survival analysis showed that there was no significant difference in prognosis between the high- and low-score groups (Supplementary Figure S6C). To further investigate the differences in immune characteristics between the two groups, the immune cell abundances of each CRC sample were also calculated by the ssGSEA algorithm based on the expression levels of the corresponding cell marker genes (Supplementary Figure S6D). The differences in immune and clinicopathological characteristics between the high- and low-score groups are visualized in Figure 4A and Supplementary Figure S7. The results showed that most of the significantly differentially expressed immune cells were highly expressed in the low-score group (Figure 4A). Moreover, it is noteworthy that several classical inhibitory immune cells were significantly enriched in the low-score group, such as regulatory T cells (Tregs) and type 2 T helper cells (Th2 cells), while classical inflammatory immune cells were enriched in the high-score group, such as type 17 T helper cells (Th17 cells). These results suggest that CRC patients with MSS status are not sensitive to immunotherapy, which may be caused by the inhibitory immune microenvironment (low-score group) or lack of immune cells (high-score group). To address this hypothesis further, the expression status of the classical chemokines, immune checkpoint molecules and cytotoxic molecules between the two risk groups was analyzed(Park et al., 2022; Propper and Balkwill, 2022; Josan and Neilan, 2023). The results showed that most of these molecules were significantly highly expressed in the low-score group (Figures 4B–D). Taken together, these results indicated that CRC patients with MSS status could be separated into two subtypes according to the MSS-related risk score, and the CRC patients in the low-score group were potentially sensitive to immunotherapy because of their characteristics of “hot” tumors but inhibitory immune microenvironments.

FIGURE 4

FIGURE 4. Differential analyses between the high- and low-score groups separated by the MSS-related risk score in the TCGA-CRC cohort. (A) The expression levels of the 28 immune cell types in the two groups were visualized by boxplot. (B–D) The expression status of classical chemokines (B), immune checkpoint molecules (C) and cytotoxic molecules (D) between the two risk groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, not significant.

Under this condition, how to reverse the inhibitory immune microenvironments of the CRC patients in the low-score groups becomes the most pressing question. The correlations between the immune cells and single-cell marker genes were then analyzed. The results showed that these genes have a strong correlation with immune cells, especially CD55 with Th2 cells and LCN2 with Th17 cells (Figure 5A). These results indicated that these single-cell marker genes might be potential target genes to regulate the immune microenvironments of CRC patients with MSS status to become synergistic targets to enhance sensitivity to immunotherapy. To search for potential therapeutic drugs targeting these genes, drug sensitivity analyses were conducted by CellMiner. As a result, a total of 956 potential effective drug-gene pairs were identified (Figure 6 and Supplementary Table S7). In addition, to investigate functional differences between the high- and low-risk groups, GSVA was conducted. The results preliminarily revealed that the high-score group was enriched in metabolism-related processes, such as oxidative phosphorylation and tyrosine metabolism, while the low-score group was enriched in genomic stability, such as mismatch repair and non-homologous end joining (Figure 5B).

FIGURE 5

FIGURE 5. Correlation analysis and GSVA in the TCGA-CRC cohort. (A) Correlation heatmap between immune cells and single-cell marker genes. (B) GSVA heatmap between the high- and low-score groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, not significant.

FIGURE 6

FIGURE 6. Drug prediction for the 29 single-cell marker genes. Sixteen representative scatter plots of the relationship between drug sensitivity and the expression levels of the 29 genes.

External validation of the characteristic differences between the high- and low-score groups in the GEO-CRC cohort

In the GEO-CRC cohort, these 29 single-cell marker genes were also significantly differentially expressed between the dMMR and pMMR groups (Supplementary Figure S8A). The risk scores of the CRC samples with pMMR status were calculated and grouped as in the TCGA-CRC cohort. Most of the single-cell marker genes were significantly differentially expressed between the high- and low-score groups (Supplementary Figure S8B). The results of survival analysis showed that the CRC patients in the low-score group had a significantly worse prognosis (Supplementary Figure S8C). Next, the immune cell abundances of CRC samples were also calculated by the ssGSEA algorithm (Supplementary Figure S8D). The differences in immune and clinicopathological characteristics between the high- and low-score groups are visualized in Figure 7A and Supplementary Figure S9. These results further validated that several classical inhibitory immune cells represented by Tregs, Th2 cells and myeloid-derived suppressor cells (MDSCs) were significantly enriched in the low-score group, while the classical inflammatory immune Th17 cells were enriched in the high-score group. In addition, most of the classical chemokines, immune checkpoint molecules and cytotoxic molecules were also significantly highly expressed in the low-score group (Figures 7B–D). Furthermore, correlation analysis results between the immune cells and single-cell marker genes revealed that these genes have a strong correlation with immune cells (Figure 8A). The GSVA results also confirmed that the high-score group was enriched in metabolism-related processes and that the low-score group was enriched in genomic stability (Figure 8B). Taken together, these results further strongly confirmed the analysis results of the TCGA-CRC cohort that pMMR MSS CRC patients could be separated into two subtypes, and the patients in the low-score group had “hot” tumors with inhibitory immune characteristics, which are potentially more sensitive to immunotherapy. These single-cell marker genes are potential synergistic targets to enhance the sensitivity of CRC patients to immunotherapy.

FIGURE 7

FIGURE 7. Differential analyses between the high- and low-score groups separated by the MSS-related risk score in the GEO-CRC cohort. (A) The expression levels of the 28 immune cell types in the two groups were visualized by boxplot. (B–D) The expression status of classical chemokines (B), immune checkpoint molecules (C) and cytotoxic molecules (D) between the two risk groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, not significant.

FIGURE 8

FIGURE 8. Correlation analysis and GSVA in the GEO-CRC cohort. (A) Correlation heatmap between immune cells and single-cell marker genes. (B) GSVA heatmap between the high- and low-score groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, not significant.

Discussion

Immunotherapy represented by ICIs has been approved as the second-line treatment for metastatic CRC with dMMR MSI status, but only a few CRC patients with pMMR MSS status can benefit from immunotherapy due to the lack of immune infiltration and low TMB status(Ganesh et al., 2019; Marcus et al., 2019; Johnson et al., 2022; Weng et al., 2022). Heterogeneities between dMMR MSI and pMMR MSS CRC patients have been preliminarily investigated at the single-cell level, and they have mainly focused on immune cell status(Bao et al., 2020). In this study, we comprehensively dissected the differences between MSI and MSS CRC patients by constructing a single-cell atlas. It is noteworthy that the proportion composition of cancer cell subpopulations was significantly different between the MSI and MSS groups. Moreover, the interactions among different cell types were significantly strengthened in the MSS group compared with the MSI group, and the outgoing signals of the cancer cells were significantly enhanced. These results emphasized the potential core roles of specific cancer cell subtypes in CRC patients with MSS status. Under this condition, the cancer cell subpopulation with a significantly higher proportion in the MSS group was defined as an MSS-specific cancer cell cluster. To elucidate the characteristics of these MSS-specific cancer cells, DEGs between MSS-specific cancer cells and other cancer cell subpopulations were identified. As a result, 29 MSS-specific single-cell marker genes were screened out.

Some studies have aimed to divide pMMR MSS CRC patients into different subgroups to screen potential patients who are sensitive to immunotherapy or a combination therapy scheme(Liang et al., 2021; Gao et al., 2022; Wang et al., 2022). However, these studies have not systematically elucidated the differences between pMMR MSS and dMMR MSI CRC patients at the single-cell level. Moreover, there is no corresponding risk score signature based on differential molecules at the single-cell level to classify CRC patients, so as to determine which patients may be potentially effective for immunotherapy. Under this condition, an MSS-related risk score system was constructed in this study based on the expression levels of the 29 specific single-cell marker genes. This risk score could efficiently divide pMMR MSS CRC patients into two subgroups with significantly different immune characteristics. The high-score group could be identified as the “cold” or “immune-desert” subtype, which is characterized by less immune cell infiltration and low expression of chemokines, immune checkpoint molecules and cytotoxic molecules. In contrast, the low-score group could be identified as the “hot” or “immune-exhausted” subtype, which is characterized by a high abundance of immune cell infiltration and high expression levels of chemokines, cytotoxic molecules and CD8 T cell exhaustion markers, such as CTLA4 and CD274. Moreover, classical inhibitory immune cells (Tregs and Th2 cells) were also enriched in the low-score group. High expression of immune checkpoint molecules is an important predictor of sensitivity to immunotherapy, and it has been reported that reversing the inhibitory immune microenvironment is helpful for enhancing the efficacy of sensitization immunotherapy(Gurjao et al., 2019; Sen et al., 2019; Hosein et al., 2022; Long et al., 2022). Based on this knowledge, these results indicated that the high-score group of pMMR MSS CRC patients could be absolutely insensitive to immunotherapy, while the low-score group CRC patients were the potentially sensitive to immunotherapy if synergistic targets could be identified to reverse their inhibitory immune microenvironment.

In addition, it is noteworthy that some of these 29 specific single-cell marker genes have been studied in tumor immunotherapy(Macor et al., 2015; Strati et al., 2017; Shen et al., 2021; Ahmed et al., 2022; Min et al., 2022), but their role in the response of CRC to immunotherapy is still unclear. Combined with the results of CellMiner, it is suggested that these genes may be potential synergistic targets for enhancing the sensitivity of immunotherapy in CRC patients.

There are still some limitations of the present study. First, whether this MSS-related risk score can efficiently divide pMMR MSS CRC patients into “cold” or “hot” subtypes needs to be verified in a large number of samples from multiple centers. Second, whether the targeted intervention of these 29 single-cell marker genes can reverse the inhibitory immune microenvironment of “hot” CRC patients and thus enhance their sensitivity to immunotherapy needs to be confirmed by basic and clinical studies.

Conclusion

Taken together, the heterogeneities between the dMMR MSI and pMMR MSS CRC patients were systematically elucidated at the single-cell level for the first time. As a result, a specific MSS cancer cell cluster and a specific MSS cancer cell gene list was identified. Next, an MSS-related risk score was constructed based on the 29 specific single-cell marker genes. This risk score model is useful to categorize CRC samples into “cold” or “hot” tumors, which will help us to judge the potential sensitivity of pMMR MSS CRC patients to immunotherapy. Moreover, this model will also provide fundamental knowledge to screen synergistic targets to sensitize pMMR MSS CRC patients to immunotherapy and improve the overall prognosis of CRC patients.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics statement

The studies involving human participants were reviewed and approved by Clinical Research Ethics Committee of the Second Xiangya Hospital, Central South University.

Author contributions

JF and SD conceived the study and analyzed the data. XJ, WC, ZC, PW, WX, and YL collected the data and clinical information. All of the authors drafted the article and approved the final manuscript.

Acknowledgments

The authors want to thank the TCGA and GEO databases for free use.

Conflict of interest

The 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/fphar.2023.1083449/full#supplementary-material

References

Adamopoulos, C., Ahmed, T. A., Tucker, M. R., Ung, P. M. U., Xiao, M., Karoulia, Z., et al. (2021). Exploiting allosteric properties of RAF and MEK inhibitors to target therapy-resistant tumors driven by oncogenic BRAF signaling. Cancer Discov. 11 (7), 1716–1735. doi:10.1158/2159-8290.Cd-20-1351

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahmed, K. M., Veeramachaneni, R., Deng, D., Putluri, N., Putluri, V., Cardenas, M. F., et al. (2022). Glutathione peroxidase 2 is a metabolic driver of the tumor immune microenvironment and immune checkpoint inhibitor response. J. Immunother. Cancer 10 (8), e004752. doi:10.1136/jitc-2022-004752

PubMed Abstract | CrossRef Full Text | Google Scholar

André, T., Lonardi, S., Wong, K. Y. M., Lenz, H. J., Gelsomino, F., Aglietta, M., et al. (2022). Nivolumab plus low-dose ipilimumab in previously treated patients with microsatellite instability-high/mismatch repair-deficient metastatic colorectal cancer: 4-year follow-up from CheckMate 142. Ann. Oncol. 33, 1052–1060. doi:10.1016/j.annonc.2022.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, X., Zhang, H., Wu, W., Cheng, S., Dai, X., Zhu, X., et al. (2020). Analysis of the molecular nature associated with microsatellite status in colon cancer identifies clinical implications for immunotherapy. J. Immunother. Cancer 8 (2), e001437. doi:10.1136/jitc-2020-001437

PubMed Abstract | CrossRef Full Text | Google Scholar

Boland, C. R., and Goel, A. (2010). Microsatellite instability in colorectal cancer. Gastroenterology 138 (6), 2073–2087. doi:10.1053/j.gastro.2009.12.064

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, F., and Parkin, D. M.African Cancer Registry Network (2022). Cancer in sub-saharan africa in 2020: A review of current estimates of the national burden, data gaps, and future needs. Lancet Oncol. 23 (6), 719–728. doi:10.1016/s1470-2045(22)00270-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciardiello, F., Ciardiello, D., Martini, G., Napolitano, S., Tabernero, J., and Cervantes, A. (2022). Clinical management of metastatic colorectal cancer in the era of precision medicine. CA Cancer J. Clin. 72 (4), 372–401. doi:10.3322/caac.21728

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganesh, K., Stadler, Z. K., Cercek, A., Mendelsohn, R. B., Shia, J., Segal, N. H., et al. (2019). Immunotherapy in colorectal cancer: Rationale, challenges and potential. Nat. Rev. Gastroenterol. Hepatol. 16 (6), 361–375. doi:10.1038/s41575-019-0126-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Wu, Z., Zhao, M., Zhang, R., Li, M., Sun, D., et al. (2022). Allosteric inhibition reveals SHP2-mediated tumor immunosuppression in colon cancer by single-cell transcriptomics. Acta Pharm. Sin. B 12 (1), 149–166. doi:10.1016/j.apsb.2021.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurjao, C., Liu, D., Hofree, M., AlDubayan, S. H., Wakiro, I., Su, M. J., et al. (2019). Intrinsic resistance to immune checkpoint blockade in a mismatch repair-deficient colorectal cancer. Cancer Immunol. Res. 7 (8), 1230–1236. doi:10.1158/2326-6066.Cir-18-0683

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

Hosein, A. N., Dangol, G., Okumura, T., Roszik, J., Rajapakshe, K., Siemann, M., et al. (2022). Loss of Rnf43 accelerates kras-mediated neoplasia and remodels the tumor immune microenvironment in pancreatic adenocarcinoma. Gastroenterology 162 (4), 1303–1318.e18. doi:10.1053/j.gastro.2021.12.273

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, S., Guerrero-Juarez, C. F., Zhang, L., Chang, I., Ramos, R., Kuan, C. H., et al. (2021). Inference and analysis of cell-cell communication using CellChat. Nat. Commun. 12 (1), 1088. doi:10.1038/s41467-021-21246-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Z., and Sinicrope, F. A. (2022). Mismatch repair-deficient colorectal cancer: Building on checkpoint blockade. J. Clin. Oncol. 40 (24), 2735–2750. doi:10.1200/jco.21.02691

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, B., Haymaker, C. L., Parra, E. R., Soto, L. M. S., Wang, X., Thomas, J. V., et al. (2022). Phase II study of durvalumab (anti-PD-L1) and trametinib (MEKi) in microsatellite stable (MSS) metastatic colorectal cancer (mCRC). J. Immunother. Cancer 10 (8), e005332. doi:10.1136/jitc-2022-005332

PubMed Abstract | CrossRef Full Text | Google Scholar

Josan, K., and Neilan, T. G. (2023). Immune checkpoint inhibitors: Acute and chronic cardiovascular complications. Nat. Rev. Cardiol. 20 (2), 73–74. doi:10.1038/s41569-022-00827-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaliq, A. M., Erdogan, C., Kurt, Z., Turgut, S. S., Grunvald, M. W., Rand, T., et al. (2022). Refining colorectal cancer classification and clinical stratification through a single-cell atlas. Genome Biol. 23 (1), 113. doi:10.1186/s13059-022-02677-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Kotani, D., Nakamura, Y., Fujisawa, T., Bando, H., Sakamoto, N., Johns, A. L., et al. (2022). ICGC-ARGO precision medicine: Targeted therapy according to longitudinal assessment of tumour heterogeneity in colorectal cancer. Lancet Oncol. 23 (4), 463–464. doi:10.1016/s1470-2045(22)00142-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H. O., Hong, Y., Etlioglu, H. E., Cho, Y. B., Pomella, V., Van den Bosch, B., et al. (2020). Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nat. Genet. 52 (6), 594–603. doi:10.1038/s41588-020-0636-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, Y. H., Chen, K. H., Tsai, J. H., Cheng, Y. M., Lee, C. C., Kao, C. H., et al. (2021). Proteasome inhibitors restore the STAT1 pathway and enhance the expression of MHC class I on human colon cancer cells. J. Biomed. Sci. 28 (1), 75. doi:10.1186/s12929-021-00769-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Zhang, Q., Xing, B., Luo, N., Gao, R., Yu, K., et al. (2022). Immune phenotypic linkage between colorectal cancer and liver metastasis. Cancer Cell 40 (4), 424–437.e5. doi:10.1016/j.ccell.2022.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, J., Wang, D., Wang, A., Chen, P., Lin, Y., Bian, J., et al. (2022). A mutation-based gene set predicts survival benefit after immunotherapy across multiple cancers and reveals the immune response landscape. Genome Med. 14 (1), 20. doi:10.1186/s13073-022-01024-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Macor, P., Secco, E., Mezzaroba, N., Zorzet, S., Durigutto, P., Gaiotto, T., et al. (2015). Bispecific antibodies targeting tumor-associated antigens and neutralizing complement regulators increase the efficacy of antibody-based immunotherapy in mice. Leukemia 29 (2), 406–414. doi:10.1038/leu.2014.185

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcus, L., Lemery, S. J., Keegan, P., and Pazdur, R. (2019). FDA approval summary: Pembrolizumab for the treatment of microsatellite instability-high solid tumors. Clin. Cancer Res. 25 (13), 3753–3758. doi:10.1158/1078-0432.Ccr-18-4070

PubMed Abstract | CrossRef Full Text | Google Scholar

Marisa, L., de Reyniès, A., Duval, A., Selves, J., Gaub, M. P., Vescovo, L., et al. (2013). Gene expression classification of colon cancer into molecular subtypes: Characterization, validation, and prognostic value. PLoS Med. 10 (5), e1001453. doi:10.1371/journal.pmed.1001453

PubMed Abstract | CrossRef Full Text | Google Scholar

Min, H. Y., Cho, J., Sim, J. Y., Boo, H. J., Lee, J. S., Lee, S. B., et al. (2022). S100A14: A novel negative regulator of cancer stemness and immune evasion by inhibiting STAT3-mediated programmed death-ligand 1 expression in colorectal cancer. Clin. Transl. Med. 12 (7), e986. doi:10.1002/ctm2.986

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, H. J., Lee, S. W., Park, Y. H., Kim, T. C., Van Kaer, L., and Hong, S. (2022). CD1d-independent NK1.1(+) Treg cells are IL2-inducible Foxp3(+) T cells co-expressing immunosuppressive and cytotoxic molecules. Front. Immunol. 13, 951592. doi:10.3389/fimmu.2022.951592

PubMed Abstract | CrossRef Full Text | Google Scholar

Picco, G., Cattaneo, C. M., van Vliet, E. J., Crisafulli, G., Rospo, G., Consonni, S., et al. (2021). Werner helicase is a synthetic-lethal vulnerability in mismatch repair-deficient colorectal cancer refractory to targeted therapies, chemotherapy, and immunotherapy. Cancer Discov. 11 (8), 1923–1937. doi:10.1158/2159-8290.Cd-20-1508

PubMed Abstract | CrossRef Full Text | Google Scholar

Propper, D. J., and Balkwill, F. R. (2022). Harnessing cytokines and chemokines for cancer therapy. Nat. Rev. Clin. Oncol. 19 (4), 237–253. doi:10.1038/s41571-021-00588-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Reinhold, W. C., Varma, S., Sunshine, M., Elloumi, F., Ofori-Atta, K., Lee, S., et al. (2019). RNA sequencing of the NCI-60: Integration into CellMiner and CellMiner CDB. Cancer Res. 79 (13), 3514–3524. doi:10.1158/0008-5472.can-18-2047

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, Y., Miao, J. M., Wang, Y. Y., Fan, Z., Kong, X. B., Yang, L., et al. (2022). Oncolytic viruses combined with immune checkpoint therapy for colorectal cancer is a promising treatment option. Front. Immunol. 13, 961796. doi:10.3389/fimmu.2022.961796

PubMed Abstract | CrossRef Full Text | Google Scholar

Satija, R., Farrell, J. A., Gennert, D., Schier, A. F., and Regev, A. (2015). Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 33 (5), 495–502. doi:10.1038/nbt.3192

PubMed Abstract | CrossRef Full Text | Google Scholar

Sen, T., Della Corte, C. M., Milutinovic, S., Cardnell, R. J., Diao, L., Ramkumar, K., et al. (2019). Combination treatment of the oral CHK1 inhibitor, SRA737, and low-dose gemcitabine enhances the effect of programmed death ligand 1 blockade by modulating the immune microenvironment in SCLC. J. Thorac. Oncol. 14 (12), 2152–2163. doi:10.1016/j.jtho.2019.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, C., Zhang, Z., Tian, Y., Li, F., Zhou, L., Jiang, W., et al. (2021). Sulforaphane enhances the antitumor response of chimeric antigen receptor T cells by regulating PD-1/PD-L1 pathway. BMC Med. 19 (1), 283. doi:10.1186/s12916-021-02161-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Strati, A., Koutsodontis, G., Papaxoinis, G., Angelidis, I., Zavridou, M., Economopoulou, P., et al. (2017). Prognostic significance of PD-L1 expression on circulating tumor cells in patients with head and neck squamous cell carcinoma. Ann. Oncol. 28 (8), 1923–1933. doi:10.1093/annonc/mdx206

PubMed Abstract | CrossRef Full Text | Google Scholar

Sui, Q., Liu, D., Jiang, W., Tang, J., Kong, L., Han, K., et al. (2021). Dickkopf 1 impairs the tumor response to PD-1 blockade by inactivating CD8+ T cells in deficient mismatch repair colorectal cancer. J. Immunother. Cancer 9 (3), e001498. doi:10.1136/jitc-2020-001498

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

Topham, J. T., O'Callaghan, C. J., Feilotter, H., Kennecke, H. F., Lee, Y. S., Li, W., et al. (2022). Circulating tumor DNA identifies diverse landscape of acquired resistance to anti-epidermal growth factor receptor therapy in metastatic colorectal cancer. J. Clin. Oncol. 41, 485–496. doi:10.1200/jco.22.00364

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Zhou, Y., Zhang, Y., Fang, S., Zhang, M., Li, H., et al. (2022). Subtyping of microsatellite stability colorectal cancer reveals guanylate binding protein 2 (GBP2) as a potential immunotherapeutic target. J. Immunother. Cancer 10 (4), e004302. doi:10.1136/jitc-2021-004302

PubMed Abstract | CrossRef Full Text | Google Scholar

Weng, J., Li, S., Zhu, Z., Liu, Q., Zhang, R., Yang, Y., et al. (2022). Exploring immunotherapy in colorectal cancer. J. Hematol. Oncol. 15 (1), 95. doi:10.1186/s13045-022-01294-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, Q., Wang, X., Ren, D., Hu, B., Tang, G., Zhang, Y., et al. (2021). DNA methylation-based signature of CD8+ tumor-infiltrating lymphocytes enables evaluation of immune response and prognosis in colorectal cancer. J. Immunother. Cancer 9 (9), e002671. doi:10.1136/jitc-2021-002671

PubMed Abstract | CrossRef Full Text | Google Scholar