Skip to main content

Single-cell sequencing reveals the immune microenvironment landscape related to anti-PD-1 resistance in metastatic colorectal cancer with high microsatellite instability



The objective response rate of microsatellite instability-high (MSI-H) metastatic colorectal cancer (mCRC) patients with first-line anti-programmed cell death protein-1 (PD-1) monotherapy is only 40–45%. Single-cell RNA sequencing (scRNA-seq) enables unbiased analysis of the full variety of cells comprising the tumor microenvironment. Thus, we used scRNA-seq to assess differences among microenvironment components between therapy-resistant and therapy-sensitive groups in MSI-H/mismatch repair-deficient (dMMR) mCRC. Resistance-related cell types and genes identified by this analysis were subsequently verified in clinical samples and mouse models to further reveal the molecular mechanism of anti-PD-1 resistance in MSI-H or dMMR mCRC.


The response of primary and metastatic lesions to first-line anti-PD-1 monotherapy was evaluated by radiology. Cells from primary lesions of patients with MSI-H/dMMR mCRC were analyzed using scRNA-seq. To identify the marker genes in each cluster, distinct cell clusters were identified and subjected to subcluster analysis. Then, a protein‒protein interaction network was constructed to identify key genes. Immunohistochemistry and immunofluorescence were applied to verify key genes and cell marker molecules in clinical samples. Immunohistochemistry, quantitative real-time PCR, and western blotting were performed to examine the expression of IL-1β and MMP9. Moreover, quantitative analysis and sorting of myeloid-derived suppressor cells (MDSCs) and CD8+ T cells were performed using flow cytometry.


Tumor responses in 23 patients with MSI-H/dMMR mCRC were evaluated by radiology. The objective response rate was 43.48%, and the disease control rate was 69.57%. ScRNA-seq analysis showed that, compared with the treatment-resistant group, the treatment-sensitive group accumulated more CD8+ T cells. Experiments with both clinical samples and mice indicated that infiltration of IL-1β-driven MDSCs and inactivation of CD8+ T cells contribute to anti-PD-1 resistance in MSI-H/dMMR CRC.


CD8+ T cells and IL-1β were identified as the cell type and gene, respectively, with the highest correlation with anti-PD-1 resistance. Infiltration of IL-1β-driven MDSCs was a significant factor in anti-PD-1 resistance in CRC. IL-1β antagonists are expected to be developed as a new treatment for anti-PD-1 inhibitor resistance.

Peer Review reports


Colorectal cancer (CRC) is one of the most common malignant digestive tract tumors and has the second highest mortality rate and the third highest incidence rate among all malignant tumors [1]. Although diagnosis and treatment strategies for CRC have improved rapidly in recent years, the prognosis is unfavorable for many CRC patients. Many patients are not diagnosed until they are already in advanced stages of disease, which is a significant obstacle to effective treatment.

The expression of proteins related to mismatch repair deficiency (dMMR) and microsatellite instability-high (MSI-H) status have been widely recognized as valuable for predicting the efficacy of immunotherapy in CRC patients. Compared to traditional cancer therapies, immunotherapy has improved the objective response rate (ORR) of MSI-H mCRC to some extent [2, 3]. However, the ORR of patients with MSI-H/dMMR mCRC with first-line anti-programmed cell death protein-1 (PD-1) monotherapy was only 43.8% [4]. Thus, more than half of patients do not benefit from immunotherapy or chemotherapy ± targeted therapy because of the MSI-H phenotype [5]. In view of this, this study aims to explore the mechanism of anti-PD-1 resistance in patients with MSI-H mCRC.

ScRNA-seq has contributed to a better understanding of the immune landscape in MSI-H patients, which in turn has offered novel insights into the mechanisms of immunotherapy. Previous studies have mainly focused on the mechanism of immune resistance in microsatellite stable (MSS) mCRC. Several studies have revealed the molecular mechanism of CD40 and CD73 antagonist therapy at the single-cell level, providing a possible theoretical basis for their combined treatment with immune checkpoint inhibitors for MSS mCRC [6,7,8]. However, the mechanism of PD-1 resistance in patients with MSI-H mCRC is still unclear. Therefore, comparing the immune microenvironment between anti-PD-1-resistant and anti-PD-1-sensitive groups by scRNA-seq could help elucidate the mechanisms underlying immunotherapy resistance in MSI-H mCRC patients.

In the present study, tissue samples were collected during colonoscopy from treatment-sensitive and -resistant groups of MSI-H mCRC treated with a PD-1 blocker (tislelizumab). Next, we comprehensively analyzed the cell subtypes and key genes of the two groups using scRNA-seq. Follow-up experiments were performed using clinical samples and mouse models to verify the potential mechanism of anti-PD-1 resistance induced by the key cell types or genes indicated by scRNA-seq.


Sample collection

23 MSI-H/dMMR mCRC patients were treated with anti-PD-1 monotherapy at Yunnan Cancer Hospital (The Third Affiliated Hospital of Kunming Medical University) between August 1, 2020, and May 31, 2022. A PD-1 blocker (200 mg, tislelizumab, Bei Gene Ltd., China) was injected intravenously on the 1st day of each 21-day cycle. Efficacy was evaluated radiologically after every third cycle of treatment. Patients were considered sensitive to anti-PD-1 treatment if they showed a complete response (CR) or partial response (PR), and patients were considered resistant to anti-PD-1 treatment if they had progressive disease (PD) or stable disease (SD). Fresh intestinal tumor tissues (2–4 mm) were taken during colonoscopy. Tissue samples from a total of 23 patients (10 sensitive and 13 resistant) were analyzed by immunohistochemistry (IHC) and immunofluorescence (IF). Six patients (3 PR and 3 PD) were randomly selected for scRNA-seq. All participants provided written informed consent before beginning the study. Additionally, the ethics committee of Yunnan Cancer Hospital approved all study protocols that involved human subjects according to the ethical principles described in the Helsinki Declaration. The inclusion criteria were as follows: primary colorectal adenocarcinoma confirmed by colonoscopy biopsy and distant metastasis confirmed by CT/MR; MSI-H/dMMR status confirmed by multiplex qPCR or IHC; undergoing first-line anti-PD-1 monotherapy. Patients were excluded if we were unable to obtain a sample through colonoscopy due to contraindications or we were unable to perform single-cell sequencing due to insufficient cell activity or quantity.

Radiology and colonoscopy

Siemens (SOMATOM Definition AS +) 128 slice spiral CT was used for plain and enhanced scanning. Patients were fasting and had performed bowel preparation as directed. The scanning range covered the entire abdominal cavity. The scanning layer was 1.0 mm thick with an interval of 0.6 mm. Iohexol (300 mg/ml, 100 ml) was used as the contrast agent. The delay time of arterial phase scanning was 35–40 s, and that of venous phase scanning was 70–80 s. MRI was performed using a Philips Elision 3.0 T, and the scanning sequence included transverse T1WI_Tse, sagittal and coronal T2WI_Tse, high-resolution T2WI_Tse, diffusion-weighted imaging, and multiphase dynamic enhancement. Gadolinium diamine was used as the contrast agent. Colonoscopy was performed using the Olympus CV-290, and tissues were collected by doctors with at least five years of experience.

Single-cell preparation

Fresh tumor tissues obtained from CRC patients were immediately transferred to MACS C-tubes (Miltenyi Biotec) with digestive enzymes. Then, digestion was performed by a gentleMACS Octo Dissociator (Miltenyi Biotec) (30 min at 37 °C). Single cells were processed by Chromium Controller (10X Genomics) according to the manufacturer's protocol. Briefly, cells were washed with 20 mL of RPMI 1640 (Gibco), filtered through a 70-μm nylon strainer (BD Falcon), collected by centrifugation (330 × g, 10 min, 4 °C), and resuspended in a basic solution containing 0.2% fetal bovine serum (FBS; Gibco).

Single-cell capture and library preparation

A 10 × Chromium system (10 × Genomics) and library preparation by LC Sciences were utilized to run the single cells according to the recommended protocol for the Chromium Single Cell 30 Reagent Kit (v2 Chemistry). The Illumina HiSeq4000 was used for sequencing, and a 10 × Cell Ranger package (v1.2.0; 10 × Genomics) was used for the postprocessing and quality control of libraries.

Quantitative analysis of the single-cell sequencing data

All single-cell sequencing data were analyzed by Cell Ranger V6.1.1. The results are shown in Additional file 1: Table S1. CRC samples from 3 resistant patients (named R1, R2, R3) and 3 sensitive patients (named S1, S2, S3) contained a total of 56,092 cells; the number of genes detected in each sample ranged from 17,564 to 18,093; the median number of cellular unique molecular identifiers ranged from 3,403 to 6,915; and the sequencing saturation ranged from 47 to 63%. These results indicate that, overall, the sequencing quality was sufficient for use in subsequent correlation analysis.

Processing of CRC single-cell sequencing data

In total, this analysis included 56,092 cells from the sensitive group and resistant group. The Seurat package in R (version 4.0.5) was used for analysis and quality control [9]. Low-quality cells (n = 3,071) were removed if genes were detected only in < 3 cells or if there were < 200 total genes detected by gene-cell matrixes. Next, the global-scaling method LogNormalize was performed to normalize the gene expression values for the remaining 53,021 cells. Then, the FindVariableFeatures function combined with the vst method in R studio was employed to identify the most variable genes, which were used for dimensionality reduction. After principal component analysis, 2000 highly variable genes were identified. Jackstraw and ScoreJackStraw functions were applied to determine the most significant principal components. Finally, graph-based unsupervised clustering was conducted and visualized using a nonlinear t-distributed stochastic neighbor embedding (t-SNE) plot, defined by the FindNeighbors and FindClusters functions.

Identification and characterization of cell subtypes

The identities of cell types were characterized using the SingleR (V1.6.1) package based on the Celldex database. The FindMarkers function in the R package Seurat was used to list the markers of each cell cluster with min.pct = 0.5, logfc.threshold = 1, min.diff.pct = 0.3, and P < 0.05. The markers used in this pipeline are listed in Additional file 2: Table S2.

To investigate the molecular mechanisms involved in each cell subtype, biological process (BP), cell composition (CC), and molecular function (MF) GO enrichment analysis and KEGG pathway analysis were performed using the R package clusterProfiler (version 3.14.3), with the threshold of significance set to adjusted (adj.) P < 0.05.

Screening of differentially expressed genes (DEGs)

The first 2000 highly variable genes were screened using the FindVariableFeatures function combined with the vst method. The FindMarkers function in the R package Seurat was used not only to find the marker genes for different cell subtypes (with screening parameter thresholds of min.pct = 0.5, logfc.threshold = 1, min.diff.pct = 0.3, P < 0.05) but also to identify DEGs of each cell subtype between anti-PD-1-sensitive and anti-PD-1-resistant samples (with the threshold set to |avg_log2FC|≥ 0.5, P ≤ 0.05). The overlapping genes between pseudotime-related genes and PD-1 resistance-related DEGs were considered to be candidate key genes.

Pseudotime analysis

To reveal differences in immune cells in the sensitive and resistant groups, Monocle software (version 2.20.0) [10] was used to analyze sample trajectories and explore the differentiation process. First, a more sophisticated method (dpFeature) was created based on the foundation of cluster and custom developmental marker genes of Monocle. The signature genes with a high degree of dispersion (q < 0.01) were identified among cell subtypes that were selected by dpFeature. Next, DDRTree was applied for dimensionality reduction and pseudotemporal alignment of cells along the trajectory, and finally, the trajectories were visualized as 2D t-SNE maps. Pseudotime-related genes were identified based on q < 0.05 of the abovementioned signature genes.

Protein‒protein interaction analysis of candidate key genes

The Search Tool for the Retrieval of Interacting Genes (STRING; was used to perform protein‒protein interaction (PPI) analysis on candidate key genes. The STRING database can be used to assess the direct (physical) and indirect (functional) associations of proteins [11]. Cytoscape 3.6.1 was used to establish a network model of PPI analysis results. Based on the STRING online tool, the PPI network of the candidate key genes was constructed with medium confidence = 0.4. The top four genes were selected as the key genes in this study based on the connectivity (degree) of each node in the PPI network.

Cell culture and transfection

Colorectal cancer cell line CT26 was selected for this study. The plasmid used for cell transfection was synthesized by Ono Company. Eighteen to twenty-four hours prior to lentivirus infection, 3 × 105/well adherent cells were spread in 6-well plates. The number of cells transfected with lentivirus was approximately 6 × 105/well. When cells adhered to the wells and were 70% confluent, the original culture medium was replaced by 2 ml fresh culture medium containing 8 μg/ml polyzoan and a proper amount of virus suspension. Cells were then incubated at 37 °C for 8 h, after which the virus-containing medium was replaced with fresh medium. If transfection was successful, fluorescent protein was visible after 48–72 h. If no fluorescence was observed, the infection protocol was repeated. Puromycin was added to screen for lentivirus overexpression for one month.

Animal experiments

All animal experiments were approved by the Animal Ethics Committee of Kunming Medical University. In the process of experimental operation, we strictly follow the ARRIVE guidelines. Throughout the experiment, researchers did not know which group the animals taken out of the cage would be assigned to, animal managers and researchers doing the experiment did not know the assignment sequence, and researchers evaluating, testing or quantifying the results of the experiment did not know the means of intervention.

BALB/c mice (male, aged 6–8 weeks old, 20-30 g) were purchased from Beijing Sipeifu Biotechnology Co., Ltd. All mice were maintained in an SPF room in the animal-housing facilities at the Kunming Medical University with food and water provided at will. The experiment was carried out after 1 week of adaptive feeding. Based on the degree of freedom (E) of variance analysis proposed by Mead, we estimated the sample size of the experimental animals we need [12]. The total sample size of mice in this study was 25, and they were randomly divided into 5 groups with 5 mice in each group. All mice were randomly divided into 5 groups using a simple random sampling method, defined as OE-NC, OE-IL-1β, OE-IL-1β + Diacerein, OE-IL-1β + Nivolumab, and OE-IL-1β + Diacerein + Nivolumab groups, respectively. For the control group model (OE-NC group), 1 × 106 empty vector stabilized cells in 100ul PBS were subcutaneously injected into the flank of mice. For the treatment model, 1 × 106 over-expressed IL-1β stably transfected cells in 100ul PBS were subcutaneously injected into the same site of mice. When the tumor volume reached 40mm3 (i.e. the 7th day), the OE-IL-1β + Diacerein group began intraperitoneal (i.p.) injection of the IL-1β antagonist (Diacerein, 0.07 mg/kg), OE-IL-1β + Nivolumab group received i.p. injection of anti-PD-1 antibody (Nivolumab, 2 mg/kg); The OE-IL-1β + Diacerein + Nivolumab group was treated with a combination of diacerein and Nivolumab intraperitoneally. The mice in the OE-NC group and OE-IL-1β group were administered the same volume of PBS. Diacerein and Nivolumab were intraperitoneally injected every 3 days from the 7th day of tumor cell inoculation.

The weight of mice and the volume of tumor were measured in 0d, 7d, 10d, 12d, 14d and 16d in each group, respectively. Tumor volume was calculated as ½ (Length × Width2). On the 16th day after inoculation of tumor cells, all mice were sacrificed with an intraperitoneal injection of sodium pentobarbital (200 mg/kg). Determine mice death based on the disappearance of corneal reflex and the emission of pupils. Tumor tissues were stripped and weighed, and subsequent molecular biological experiment were performed. The selection of 16d as the experimental endpoint is based on the pre- experiment results. The selection of 16d as the experimental endpoint is based on the pre- experiment results. Throughout the entire experiment, the length of tumor should not exceed 20 mm, and the tumor weight should not exceed 10% of the mice weight [13].

Quantitative real-time PCR

Total RNA was extracted from the cultured cell line CT26 or tissues using TRIzol reagent (Ambion) and reverse transcribed into cDNA by the SureScript first strand cDNA synthesis kit (Servicebio) according to the manufacturer’s instructions. 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio) and CFX96 sequence detection system (Bio-Rad, Hercules, CA, USA) were used for qPCR, and the following primers were used: IL-1β (human), forward: 5'- AATCTCCGACCACCACTACA-3' and reverse: 5'-GACAAATCGCTTTTCCATCT-3'; MMP9 (human), forward: 5'-ATGAGCCTCTGGCAGCCCCTGGTCC-3' and reverse: 5'- GGACCAGGGGCTGCCAGAGGCTCAT-3'; GAPDH (human), forward: 5'-CCCATCACCATCTTCCAGG-3' and reverse: 5'-CATCACGCCACAGTTTCCC-3'; IL-1β (mouse), forward: 5'- CCTATGTCTTGCCCGTGG-3' and reverse: 5'- GTGGGTGTGCCGTCTTTC-3'; MMP9 (mouse), forward: 5'-GTGTGTTCCCGTTCATCTTT-3' and reverse: 5'- GCCGTCTATGTCGTCTTTAT-3'; GAPDH (mouse), forward: 5 '- CCTTCCGTGTTCCTACCCC-3' and reverse: 5'-GCCCAAGATGCCCTTCAGT-3'. GAPDH was used as a standardized endogenous control, and 2CT was used to calculate the relative mRNA expression.

Western blotting

Protein samples were isolated from tissues or cells using RIPA lysis buffer (Servicebio, Wuhan, China) containing 1% protease and phosphatase inhibitors (PMSF; Servicebio). A BCA protein assay kit (Biyuntian Biotechnology) was used for protein quantification. Sodium dodecyl sulfate‒polyacrylamide gel electrophoresis was applied to separate proteins of different molecular weights, and proteins were transferred to a polyvinylidene fluoride (Servicebio) membrane. The membrane was blocked with 5% skim milk for 90 min at room temperature and subsequently incubated with primary antibodies (Proteintech; IL-1β 1:1000; MMP9 1:1000; β-actin 1:25,000) at 4 °C overnight, followed by incubation with HRP-conjugated secondary antibodies for 2 h at room temperature.

Immunohistochemistry and immunofluorescence

Paraffin sections of tissues were deparaffinized and rehydrated. Then, sodium citrate buffer was used to extract the antigens to be detected, and antigen retrieval was completed by heat induction. Sectioned tissues were incubated with 3% H2O2 for 15 min to block endogenous peroxidase activity (IF proceeded without this step) and then blocked with PBS containing 5% fetal bovine serum for 30 min. Tissues were incubated with primary antibodies overnight at 4 degrees C, followed by incubation in the dark with the conjugated secondary antibodies at room temperature for another 2 h. DAB was used as nuclear markers for IHC. DAPI (EX:330-380 nm, Em:420 nm) was used to stain the cell nuclei (blue), Alexa Fluor 488 (EX:495 nm, Em:519 nm) was used to stain CD11b (green), Alexa Fluor 555(EX:555 nm, Em:565 nm) was used to stain CD14, CD15 and CD8 (red), Alexa Fluor 594 (EX:590 nm, Em:617 nm) was used to stain CD33 (orange).

Flow cytometry

Polymorphonuclear (PMN)-MDSCs/monocytic (M)-MDSCs were stained with CD11b-FITC (Biolegend, 101,205, USA), Ly-6G-PE (Biolegend, 127,607, USA), and Ly-6C-APC (Biolegend, 128,016, USA), and CD8+ T cells were stained with CD3-FITC (Biolegend, 100,203, USA) and CD8-PE (Biolegend, 100,707, USA) according to the manufacturer’s instructions. Samples were run on a Guava easyCyte 8HT flow cytometer (Millipore). Forward and side scatter gating were performed using FlowJo_V10.

Statistical analysis

Bioinformatics analyses were performed using R software. The varElect online tool was applied to analyze the correlation between genes in the PPI network and CRC. A high score indicated a strong correlation, with a significance threshold of P < 0.05 unless otherwise stated. All data are presented as the mean ± standard error (SE) of independent experiments. Two-tailed one-way analysis of variance (ANOVA) with multiple comparison post hoc analysis was used, and P values < 0.05 (*), P < 0.01 (**), P < 0.001 (***), and P < 0.0001 (****) are indicated as significant. Statistical analysis was performed using GraphPad Prism 9.0.


Efficacy evaluation for MSI-H/dMMR mCRC after anti-PD-1 monotherapy

A total of 23 MSI-H/dMMR mCRC patients were treated with first-line anti-PD-1 monotherapy between August 1, 2020, and May 31, 2022. Treatment response was evaluated by radiological examination after every 3 cycles of PD-1 inhibitor. PR was recorded for seven patients, CR was recorded for three patients, SD was recorded for six patients, and PD was recorded for seven patients. The ORR was 43.48% (10/23), and the disease control rate (DCR) was 69.57% (16/23) (Fig. 1A). The radiological findings for primary and metastatic lesions in six patients pre- and post-immunotherapy can be seen in Fig. 1B and C. Figure 1B shows the changes in primary and metastatic lesions in 3 patients in the resistance group (R1, R2 and R3) after anti-PD-1 treatment. The length of the primary lesions of R1 increased from 2.1 cm to 3.2 cm, while that of the metastatic lesions increased from 1.6 cm to 2.7 cm; the primary and metastatic lesions of R2 increased from 0.5 cm to 1.9 cm and from 0.3 cm to 0.8 cm, respectively; the length of primary lesions of R3 increased from 1.5 cm to 1.9 cm, and the number of metastatic lesions increased from 3 to more than 20. The responses of R1, R2 and R3 were evaluated as PD. Similarly, Fig. 1C shows the changes in primary and metastatic lesions in 3 patients in the sensitive group (S1, S2 and S3) after anti-PD-1 monotherapy. The length of the primary tumor of S1 decreased from 2.0 cm to 1.5 cm and that of the metastatic tumor decreased from 2.1 cm to 1.2 cm; the lengths of the primary and metastatic lesions of S2 decreased from 1.2 cm to 0.7 cm and from 5.3 cm to 4.4 cm, respectively; and the lengths of the primary and metastatic lesions of S3 decreased from 1.3 cm to 0.8 cm and from 2.1 cm to 0.5 cm, respectively. The responses of S1, S2 and S3 were evaluated as PR.

Fig. 1
figure 1

Efficacy evaluation of MSI-H mCRC after first-line PD-1 monotherapy. A Flow chart for screening 23 patients receiving first-line PD-1 monotherapy. Images of primary and metastatic lesions pre- and post-immunotherapy are shown for (B) three resistant patients and (C) three sensitive patients

Identification of 23 cell clusters based on single-cell sequencing data from CRC samples

A total of six patients (three PR and three PD) were randomly selected for scRNA-seq. 10 × Genomics scRNA sequencing datasets were obtained from fresh CRC tissues obtained from three resistant (named R1, R2, R3) and three sensitive (named S1, S2, S3) patients. After removing 3071 low-quality cells, a total of 53,021 cells were used in the final analysis (Fig. 2A); specifically, there were 7679 cells from R1, 10,797 cells from R2, 9020 from R3, 7880 cells from S1, 9369 cells from S2, and 8276 cells from S3. Figure 2B shows the top 2000 highly variant genes, with the ten most variable genes labeled. PCA of the 2000 highly variable genes demonstrated no significant separation in CRC cells between resistant and sensitive groups (Fig. 2C). Twenty principal components were selected based on linear dimensionality reduction analysis (Fig. 2D). Nonlinear dimensionality reduction of the data according to the dimension value 20 combined with the RunMap function revealed a more uniform distribution of cells in each sample (Fig. 2E). Subsequently, these cells were classified into 23 cell clusters based on gene expression levels by t-SNE, and the distribution of these cellular taxa was found to be nearly identical across the sensitive and resistant groups (Fig. 2F).

Fig. 2
figure 2

Cell clusters in six CRC samples were identified by scRNA sequencing. A Number of cells in different samples. B Display of the top 2000 variant genes. C PCA of the top 2000 variant genes. D Linear dimensionality reduction analysis of principal components. E Analysis of cell distributions by nonlinear dimensionality reduction combined with the RunMap function. F Analysis of cell cluster distributions in the sensitive and resistant groups through t-SNE

Characterization of the nine cell subtypes

Using the R packages SingleR and celldex, 23 cell clusters were annotated into nine cell subtypes: CD8 + T cells, epithelial cells, B cells, dendritic cells (DCs), hematopoietic stem cells, monocytes, fibroblasts, myocytes, and endothelial cells (Fig. 3A). Combined with Fig. 2F, the aggregation of CD8+ T cells and monocytes in the sensitive group was significantly higher than that in the resistant group. The number of cells of each cellular subtype in each sample is shown in Additional file 3: Table S3. A total of 679 marker genes were obtained for the nine cell subtypes (Additional file 2: Table S2): there were 176 markers for fibroblasts, 143 markers for endothelial cells, 131 markers for myocytes, 67 markers for DCs, 65 markers for monocytes, 40 markers for hematopoietic stem cells, 27 markers for epithelial cells, 17 markers for CD8+ T cells, and 13 markers for B cells (Fig. 3B). The top gene for each cell subtype is shown in Fig. 3C. The heatmap of these genes illustrated that each cluster displayed distinct gene expression features (Fig. 3D).

Fig. 3
figure 3

Annotation of cell subtypes and functional enrichment analysis of marker genes. A Annotation of cell subtypes according to the R package SingleR and celldex. B Number of marker genes for each cell subtype. C Scatter plot showing the top gene marker for each subtype. D Heatmap showing the top gene maker for each subtype

GO analysis showed that these marker genes were mainly enriched in immune-related biological processes such as differentiation of lymphocytes and monocytes and positive regulation of leukocytes. Moreover, they were also significantly associated with “cytokine-mediated signaling pathway” (Additional file 4: Figure S1A-C). KEGG enrichment analysis (Additional file 4: Figure S1D) indicated that the marker genes were significantly correlated with immune/inflammatory responses (e.g.“cytokine − cytokine receptor interaction”, “Fc epsilon RI signaling pathway”, “AGE − RAGE signaling pathway”, “cell adhesion molecules”). Cancer-related pathways, including “PI3K-Akt signaling pathway” and “proteoglycans in cancer”, were also significantly enriched.

Simulation of cell cluster developmental trajectories and screening of key genes

Developmental trajectories were simulated for the nine cell subtypes, and 1623 feature genes were found to be differentially expressed among the cell subpopulations (Additional file 5: Table S4). Figure 4A illustrates the relatively high dispersion of these feature genes. Individual cells were then classified by these feature genes using the R package Monocle, and a tree structure of the entire spectrum of differentiation trajectories was constructed (Fig. 4B). From a cell typing perspective, overall, the proposed temporal evolutionary trend of the cell subtypes was a gradual transition from epithelial cells to intermediate cells (e.g., fibroblasts, hematopoietic stem cells, monocytes, B cells) and eventually to CD8+ T cells (Fig. 4C). A total of 1454 pseudotime-related genes were identified from the 1623 feature genes (P < 0.05; Additional file 6: Table S5). KEGG analysis (Additional file 7: Figure S2A) revealed that these genes were involved in “cytokine‒cytokine receptor interaction”, “viral protein interaction with cytokine and cytokine receptor”, and “cell cycle”; they were also inextricably linked to the tumor-associated pathways “NF-kappa B signaling pathway”, “p53 signaling pathway”, and “PPAR signaling pathway”. GO-BP analysis revealed that these pseudotime-related genes were related to the immune response, inflammatory response, immune cell differentiation, proliferation, migration, and chemotaxis (Additional file 7: Figure S2B). Moreover, these genes also performed molecular functions such as “antigen binding”, “immunoglobulin receptor binding”, and “extracellular matrix structural constituent” (Additional file 7: Figure S2C) in cellular fractions such as “immunoglobulin complex”, “external side of plasma membrane”, and “immunoglobulin complex, circulating” (Additional file 7: Figure S2D).

Fig. 4
figure 4

Construction of cell clusters developmental trajectories and identification of key genes. A Relative dispersion of feature genes. B Construction of the tree structure of the differentiation trajectory spectrum. C Temporal evolutionary trend of cell subtypes. D Identification of common genes through overlap analysis. E A PPI network based on 130 common genes

To observe differences in gene expression for each cell subtype between the sensitive and resistant groups, we performed differential analysis using the FindMarkers function in the Seurat package. A total of 155 DEGs were obtained; 97 were upregulated, and 98 were downregulated (Additional file 8: Table S6). Functional enrichment analysis indicated that these genes were related to immune and inflammatory responses. DEGs in the CD8+ T-cell subtype were mainly involved in the IFN-related response and oxygen transport, while DEGs in the DC and fibroblast subtypes were associated with the chemokine-related response and neutrophil migration (Additional file 9: Figure S3A and B). Enrichment results for GO-CC and GO-MF are presented in Additional file 9: Figure S3B and C. KEGG enrichment analysis showed that DEGs in the B-cell and monocyte subtypes were involved in similar pathways and were enriched in antigen processing and presentation; DEGs in the DC and fibroblast subtypes were involved in cancer-related pathways, including chemokine and IL-17 signaling (Additional file 9: Figure S3D).

Through overlap analysis, 130 common genes were identified among 1454 pseudotime-related genes and 155 (deduplicated) PD-1 resistance-related genes (Fig. 4D; Additional file 10: Table S7). A PPI network containing 109 nodes and 435 edges was drawn based on the 130 common genes using the STRING database (Fig. 4E). VarElect analysis of the 130 common genes identified the two genes with the highest connectivity in the PPI network as IL-1β (score = 17.72) and MMP9 (score = 13.45), indicating that these two genes are most closely associated with CRC (Additional file 11: Table S8; Additional file 12: Figure S4). KEGG pathway enrichment analysis of these two key genes showed that they are involved in cytokine‒cytokine receptor interaction (IL-1β) and the MAPK and PI3K-Akt signaling pathways (IL-1β and MMP9). The specific signaling pathway maps are shown in Additional file 12: Figure S5.

Relationship between IL-1β and MDSCs or CD8+ T cells in CRC patients

IL-1β and MMP9 were identified as the top two genes with the highest correlation with anti-PD-1 resistance through scRNA-seq. We used IHC to detect the expression of IL-1β and MMP9 in 10 sensitive and 13 resistant tumor tissues from MSI-H/dMMR patients. The IHC grade of IL-1β in IL-1β-positive tissues was significantly higher than that in IL-1β-negative tumor tissues (P < 0.001; Fig. 5A). Next, all 23 patients with were categorized into IL-1β-negative or IL-1β-positive groups based on the IHC grade of IL-1β. Of the 13 tissues from resistant patients, eleven tissues had high expression of IL-1β; of the 10 tissues from sensitive patients, only one tissue had high expression (Fig. 5B). A similar expression trend was observed for another key gene, MMP9 (P < 0.0001; Fig. 5C). MMP9 expression was positively correlated with IL-1β expression (R = 0.6945, P < 0.0001; Fig. 5D).

Fig. 5
figure 5

Relationship between IL-1β and MDSCs or CD8+ T cells in MSI-H/dMMR CRC patients. A IHC staining for IL-1β in CRC tissues. B Proportion of IL-1β-positive and IL-1β-negative tissues in the sensitive and resistant groups. C IHC staining for MMP9 in IL-1β-positive and -negative groups. D Correlation between MMP9 and IL-1β according to IHC scores. IF examination and quantification of (E) PMN-MDSCs, (F) M-MDSCs, and (G) CD8+ T cells in the IL-1β-negative and -positive groups. Pearson correlation analysis of IF scores between IL-1β and (H) PMN-MDSCs, (I) M-MDSCs, and (J) CD8.+ T cells. ****P < 0.0001

MDSCs are heterogeneous cells derived from bone marrow that lead to the inactivation of CD8+ T cells and immune resistance. Multiple studies have shown that IL-1β might play crucial roles in the aggregation and differentiation of MDSCs. We further examined the relationship between IL-1β expression and MDSCs using IF. The fluorescence intensity of makers in PMN-MDSCs and M-MDSCs in the tissues of IL-1β-positive patients was significantly higher than that in IL-1β-negative patients (P < 0.0001; Fig. 5E-F). CD8+ T cells were identified as one of the most significantly different cell clusters between the sensitive group and resistant group by scRNA-seq. CD8+ T cells in the tumor microenvironment (TME) are essential for the antitumor effects of immunotherapy. We used IF to detect the expression of CD8+ T cell marker in tissues and assessed the correlation between IL-1β and CD8+ T cells. The fluorescence intensity of CD8 in IL-1β positive tissue was significantly lower than that in IL-1β negative tissue (P < 0.0001; Fig. 5G). There was a significant correlation between IL-1β and PMN-MDSCs (R = 0.8168, P < 0.0001; Fig. 5H), M-MDSCs (R = 0.7604, P < 0.0001; Fig. 5I) or CD8+ T cells (R = 0.7684, P < 0.0001; Fig. 5J).

Combining an IL-1β antagonist with an anti-PD-1 inhibitor overcame anti-PD-1 resistance in a xenograft model

To further demonstrate the influence of IL-1β on PD-1 resistance in colorectal cancer, an in vivo xenograft model was used. The human IL-1β gene was stably transfected into CT26 cell lines to increase the expression of IL-1β, as verified by qPCR and western blotting (Additional file 12: Figure S6). IL-1β-overexpressing or untransfected control CT26 cells were then subcutaneously injected into male BALB/c mice (n = 5 each group; Fig. 6A and B). After 7 days of tumor cell inoculation, progressive tumor growth was observed. As shown in Fig. 6C-D, IL-1β upregulation resulted in a substantial increase in tumor volume, and there was no significant difference in mouse weight among the five groups. We treated the three groups with diacerein, nivolumab or both and observed the growth of tumors to determine whether the IL-1β antagonist cooperates with the immune agent to inhibit the growth of tumors. IL-1β-induced elevation of CRC tumor volume (745.74 ± 188.34) and weight (0.73 ± 0.16) was greatly attenuated by diacerein (619.59 ± 127.03, 0.49 ± 0.09) or nivolumab (540.47 ± 90.92, 0.49 ± 0.10, P < 0.05), and treatment with both drugs led to the greatest attenuation of tumor growth (355.49 ± 39.89, 0.32 ± 0.08, P < 0.001). Diacerein and nivolumab showed marked synergistic effects. Tumor sections from the inoculated mice are shown in Fig. 6E.

Fig. 6
figure 6

IL-1β antagonist and PD-1 inhibitor have a synergistic antitumor effect. A Diagram of the experimental timeline. B Xenograft tumors of BALB/c mice 16 days after inoculation with CT26 cells containing an empty vector or IL-1β overexpression vector (n = 5). C Tumor weight (n = 5). D Average body weight and average tumor volume of BALB/c mice in each treatment group (n = 5). E Hematoxylin–eosin staining of tumors from inoculated BALB/c mice (n = 5). *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001

IL-1β antagonism attenuates resistance to PD-1 by regulating MDSCs and CD8+ T cells

To examine whether the effectiveness of diacerein treatment against IL-1β-induced nivolumab resistance was mediated by MDSCs, flow cytometry was performed using biomarkers for PMN-MDSCs, M-MDSCs, and CD8+ T cells. As shown in Fig. 7, IL-1β resulted in a significant increase in the number of PMN-MDSCs (17.62 ± 5.56 vs. 4.64 ± 0.84, P < 0.0001) and M-MDSCs (5.61 ± 1.35 vs. 1.26 ± 0.33, P < 0.0001) in the TME, decreased the proportion of CD8+ T cells (2.94 ± 1.17 vs. 4.35 ± 1.37, P < 0.01) in the tumor tissues of mice, and enhanced immunosuppression. Next, we further assessed whether diacerein and nivolumab showed similar synergistic effects in reversing the IL-1β-mediated effect in PMN-MDSCs, M-MDSCs and CD8+ T cells. Monotherapy with either diacerein or nivolumab suppressed the differentiation of PMN-MDSCs (11.15 ± 2.90 vs. 17.621 ± 5.561, 12.0 ± 1.88 vs. 17.621 ± 5.561, P < 0.01) and M-MDSCs (3.21 ± 0.77 vs. 1.26 ± 0.33, 4.09 ± 0.90 vs. 1.26 ± 0.33, P < 0.01) and increased the number of CD8+ T cells (4.16 ± 0.76 vs. 2.94 ± 1.17, 4.81 ± 1.07 vs. 2.94 ± 1.17, P < 0.01). Figure 7B and C suggest that compared with nivolumab alone, the combined treatment significantly reduced the number of PMN-MDSCs (8.79 ± 2.98 vs. 12.0 ± 1.89, P < 0.01) and M-MDSCs (2.26 ± 0.72 vs. 4.09 ± 0.90, P < 0.001). Figure 7D shows that the degree of CD8+ T cell aggregation in the monotherapy groups with nivolumab (4.81 ± 1.07 vs. 5.695 ± 0.90, P < 0.05) was substantially reduced compared to that in the combination therapy group.

Fig. 7
figure 7

IL-1β stimulates the aggregation of MDSCs and inhibits the proliferation of CD8+ T cells in tumors. A Flow cytometry was performed to detect PMN-MDSCs, M-MDSCs, and CD8 + T cells (n = 5). Quantitative statistical diagram for expression of (B) PMN-MDSCs, (C) M-MDSCs, (D) CD8.+ T cells in each group (n = 5). *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001

Correlation between IL-1β and MMP9 expression levels in a mouse xenograft model

As one of the key genes related to resistance identified by single-cell sequencing, MMP9 is considered to have strong antivascular and immunosuppressive functions. However, the correlation between MMP9 and IL-1β remains unclear. MMP9 was detected in mouse CRC tissue using qPCR, western blotting, and IHC. As shown in Fig. 8, IL-1β treatment significantly increased MMP9 expression in mouse tumor tissues on the basis of the results of qPCR (9.57 ± 0.26 vs. 1.08 ± 0.44, P < 0.0001) and western blotting (1.13 ± 0.07 vs. 0.21 ± 0.04, P < 0.0001). Next, we assessed whether diacerein could cooperate with nivolumab to reduce the expression of MMP9. Monotherapy with either diacerein (4.21 ± 0.60 vs. 9.57 ± 0.26, P < 0.01) or nivolumab (6.36 ± 1.10 vs. 9.57 ± 0.26, P < 0.01) attenuated the IL-1β-induced upregulation of MMP9 through qPCR analysis. Combination therapy reduced MMP9 expression more potently than monotherapy with diacerein (1.62 ± 0.59 vs. 4.21 ± 0.60, P < 0.0001) or nivolumab (1.62 ± 0.59 vs. 6.36 ± 1.10, P < 0.0001). The expression of MMP9 in tumor tissue was almost positively correlated with that of IL-1β (the value of R2 is shown in Fig. 8I).

Fig. 8
figure 8

Correlation between IL-1β and MMP9 in the mouse xenograft model. qPCR of (A) IL-1β and (B) MMP9 mRNA expression in each treatment group (n = 3). C Western blotting analysis (n = 3). Grayscale analysis of (D) IL-1β and (E) MMP9 protein (n = 3). IHC staining and scores of (F, G) IL-1β and (H, I) MMP9 expression (n = 5). (J) Pearson correlation analysis of MMP9 and IL-1β expression (n = 5). *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001


Compared with traditional chemotherapy or targeted therapy, immunotherapy has a unique molecular mechanism by which it shrinks tumors. It is generally considered that immune checkpoint inhibitors enhance the body’s general immunity by inducing T-cell activation and “releasing the immune brake” in the TME [14, 15]. In a previous study from our center, 29 patients with MSI-H/dMMR locally advanced colorectal cancer received neoadjuvant immunotherapy with a single-agent PD-1 inhibitor, and the ORR was 100% (29/29), consistent with the results of the NICHE study [16, 17]. In the present study, 23 patients with MSI-H/dMMR mCRC received first-line PD-1 inhibitor treatment, but the ORR was only 43.5% (10/23). These divergent responses of MSI-H mCRC to anti-PD-1 monotherapy are likely caused by differences in TME composition [18,19,20].

In view of this, through scRNA-seq, we defined various cellular subtypes, screened key genes, and revealed the TME signaling pathways that differ between the anti-PD-1 treatment-sensitive and -resistant groups. Together, the results suggested that the proportions of CD8+ T-cell subsets in the sensitive group were significantly higher than those in the resistant group [21,22,23,24]. Multiple CD8+ T-cell subsets in the TME may affect the response rate of new therapies targeting the immune system [21, 22]. MDSCs are heterogeneous cells derived from bone marrow that promote angiogenesis and immunosuppression [25]. Recent studies have shown that decreasing the aggregation of MDSCs in the TME can significantly increase the infiltration of CD8+ T cells and enhance the antitumor effects of PD-1 inhibition [26, 27]. The immunosuppressive effect of MDSCs might be related to the secretion of various cytokines, such as inducible nitric oxide synthase, arginase 1, interleukin-6, and transforming growth factor-β [28].

In this study, IL-1β ranked first in terms of connectivity in the overlap of pseudotime-related genes and PD-1 resistance-related genes. Monocyte-derived macrophages are a major source of IL-1β production during the immune response to pathogen infection [29, 30]. Previous studies have suggested that IL-1β could promote the invasion and growth of human colon cancer cells [31], and IL-1β polymorphisms are associated with CRC recurrence [32]. IL-1β blockade has been shown to reverse immunosuppression and has been shown to exhibit synergy with PD-1 inhibitors to promote the elimination of several tumor types, including breast [32], pancreatic [18], and renal [33] tumors. Nevertheless, the precise role of IL-1β in CRC immunotherapy has yet to be elucidated. Several studies have suggested that IL-1β can induce the infiltration of MDSCs by producing granulocyte colony stimulating factor, various CXC chemokines, and vascular adhesion molecules [34]. However, the role of IL-1β-mediated MDSC aggregation in CRC immunotherapy remains unclear.

MMP9 is a matrix metalloproteinase that has been identified as a component of the angiogenic switch during carcinogenesis [35, 36]. MMP9 mediates tumor invasion, metastasis, and immune escape through a pro-oncogenic signaling pathway [37,38,39] and is associated with relapse and poor prognosis in patients with CRC [40]. The combination of MMP9 inhibition with immune checkpoint inhibition enhanced the efficacy of immunotherapy in mouse models of melanoma [41] and breast cancer [42]. MMP9 has also shown promise in the stratification of prognosis and immune checkpoint treatment responsiveness in patients with hepatocellular carcinoma [43]. The present KEGG analysis results for 130 key genes revealed that IL-1β and MMP9 were highly correlated with the MAPK and PI3K-Akt signaling pathways. Recent studies have also indicated that IL-1β upregulates MMP9 expression through different signaling pathways in different diseases [44,45,46]. The infiltration of MDSCs could further stimulate the secretion of MMP9 [47], and IL-1β-driven accumulation of MDSCs might be one of the main sources of MMP9. However, the function of MMP9 upregulation mediated by IL-1β in CRC immunotherapy needs further study.

In the present study, we described the cellular landscape of immunotherapy-resistant and immunotherapy-sensitive groups at the single-cell level. The degree of aggregation of CD8+ T cells and monocytes in the sensitive group was significantly higher than that in the resistant group. IL-1β and MMP9 were identified as the two genes with the highest correlation with anti-PD-1 resistance. Moreover, IL-1β-driven infiltration of MDSCs enhanced anti-PD-1 resistance in MSI-H/dMMR CRC.


In the present study, the ORR and DCR of MSI-H/dMMR mCRC treated with first-line PD-1 monotherapy were 43.75% (7/16) and 68.75% (11/16), respectively. IL-1β and CD8+ T cells were found to have the highest correlation with anti-PD-1 resistance among genes and cell types, respectively. Mouse experiments demonstrated that IL-1β-driven MDSC infiltration suppressed the accumulation of CD8+ T cells, which enhanced the anti-PD-1 resistance of MSI-H/dMMR CRC. Together, these findings suggest that IL-1β antagonists may prove promising as new drugs to reverse resistance to PD-1 inhibitors.

Availability of data and materials

The scRNA-seq datasets generated and analysed during the current study are available in the NCBI repository [].



Biological process


Cell composition


Complete response


Colorectal cancer


Disease control rate


Dendritic cells


Differentially expressed genes


Deficient mismatch repair






Myeloid-derived suppressor cells


Metastatic colorectal cancer


Monocytic myeloid-derived suppressor cells


Microsatellite instability-high


Microsatellite stable


Objective response rate


Progressive disease


Programmed cell death protein-1


Polymorphonuclear myeloid-derived suppressor cells


Partial response


Perform protein–protein interaction


Molecular function


Single-cell RNA sequencing


Stable disease


Search Tool for the Retrieval of Interacting Genes


Tumor microenvironment


T-distributed stochastic neighbor embedding


  1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–49.

    Article  PubMed  Google Scholar 

  2. Overman MJ, McDermott R, Leach JL, Lonardi S, Lenz HJ, Morse MA, et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (CheckMate 142): an open-label, multicentre, phase 2 study. Lancet Oncol. 2017;18(9):1182–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Andre T, Amonkar M, Norquist JM, Shiu KK, Kim TW, Jensen BV, et al. Health-related quality of life in patients with microsatellite instability-high or mismatch repair deficient metastatic colorectal cancer treated with first-line pembrolizumab versus chemotherapy (KEYNOTE-177): an open-label, randomised, phase 3 trial. Lancet Oncol. 2021;22(5):665–77.

    Article  CAS  PubMed  Google Scholar 

  4. André T, Shiu KK, Kim TW, Jensen BV, Jensen LH, Punt C, et al. Pembrolizumab in microsatellite-instability-high advanced colorectal cancer. N Engl J Med. 2020;383(23):2207–18.

    Article  PubMed  Google Scholar 

  5. Cercek A, Dos Santos Fernandes G, Roxburgh CS, Ganesh K, Ng S, Sanchez-Vega F, et al. Mismatch repair-deficient rectal cancer and resistance to Neoadjuvant chemotherapy. Clin Cancer Res. 2020;26(13):3271–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Zhang L, Li Z, Skrzypczynska KM, Fang Q, Zhang W, O’Brien SA, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell. 2020;181(2):442-459.e429.

    Article  CAS  PubMed  Google Scholar 

  7. Vonderheide RH. CD40 agonist antibodies in cancer immunotherapy. Annu Rev Med. 2020;71:47–58.

    Article  CAS  PubMed  Google Scholar 

  8. Kim M, Min YK, Jang J, Park H, Lee S, Lee CH. Single-cell RNA sequencing reveals distinct cellular factors for response to immunotherapy targeting CD73 and PD-1 in colorectal cancer. J Immunother Cancer. 2021;9(7).

  9. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36(5):411–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32(4):381–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(Database issue):D447-452.

    Article  CAS  PubMed  Google Scholar 

  12. Dell RB, Holleran S, Ramakrishnan R. Sample size determination. Ilar J. 2002;43(4):207–13.

    Article  CAS  PubMed  Google Scholar 

  13. Han Y, Sun J, Yang Y, Liu Y, Lou J, Pan H, et al. TMP195 exerts antitumor effects on colorectal cancer by promoting M1 macrophages polarization. Int J Biol Sci. 2022;18(15):5653–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Wu T, Wu X, Wang HY, Chen L. Immune contexture defined by single cell technology for prognosis prediction and immunotherapy guidance in cancer. Cancer Commun (Lond). 2019;39(1):21.

    Article  PubMed  Google Scholar 

  15. Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019;16(6):361–75.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Zhang X, Wu T, Cai X, Dong J, Xia C, Zhou Y, et al. Neoadjuvant Immunotherapy for MSI-H/dMMR locally advanced colorectal cancer: new strategies and unveiled opportunities. Front Immunol. 2022;13:795972.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Chalabi M, Fanchi LF, Dijkstra KK, Van den Berg JG, Aalbers AG, Sikorska K, et al. Neoadjuvant immunotherapy leads to pathological responses in MMR-proficient and MMR-deficient early-stage colon cancers. Nat Med. 2020;26(4):566–76.

    Article  CAS  PubMed  Google Scholar 

  18. Zhang L, Yu X, Zheng L, Zhang Y, Li Y, Fang Q, et al. Lineage tracking reveals dynamic relationships of T cells in colorectal cancer. Nature. 2018;564(7735):268–72.

    Article  CAS  PubMed  Google Scholar 

  19. Lin A, Zhang J, Luo P. Crosstalk between the MSI status and tumor microenvironment in colorectal cancer. Front Immunol. 2020;11:2039.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Zhang Y, Xu J, Zhang N, Chen M, Wang H, Zhu D. Targeting the tumour immune microenvironment for cancer therapy in human gastrointestinal malignancies. Cancer Lett. 2019;458:123–35.

    Article  CAS  PubMed  Google Scholar 

  21. St Paul M, Ohashi PS. The roles of CD8(+) T cell subsets in antitumor immunity. Trends Cell Biol. 2020;30(9):695–704.

    Article  CAS  PubMed  Google Scholar 

  22. Zebley CC, Gottschalk S, Youngblood B. Rewriting history: epigenetic reprogramming of CD8(+) T cell differentiation to enhance immunotherapy. Trends Immunol. 2020;41(8):665–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Ugel S, Canè S, De Sanctis F, Bronte V. Monocytes in the tumor microenvironment. Annu Rev Pathol. 2021;16:93–122.

    Article  CAS  PubMed  Google Scholar 

  24. Schetters STT, Rodriguez E, Kruijssen LJW, Crommentuijn MHW, Boon L, Van den Bossche J, et al. Monocyte-derived APCs are central to the response of PD1 checkpoint blockade and provide a therapeutic target for combination therapy. J Immunother Cancer. 2020;8(2).

  25. Peranzoni E, Zilio S, Marigo I, Dolcetti L, Zanovello P, Mandruzzato S, et al. Myeloid-derived suppressor cell heterogeneity and subset definition. Curr Opin Immunol. 2010;22(2):238–44.

    Article  CAS  PubMed  Google Scholar 

  26. de Coaña YP, Wolodarski M, Poschke I, Yoshimoto Y, Yang Y, Nyström M, et al. Ipilimumab treatment decreases monocytic MDSCs and increases CD8 effector memory T cells in long-term survivors with advanced melanoma. Oncotarget. 2017;8(13):21539–53.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Kim W, Chu TH, Nienhüser H, Jiang Z, Del Portillo A, Remotti HE, et al. PD-1 signaling promotes tumor-infiltrating myeloid-derived suppressor cells and gastric tumorigenesis in mice. Gastroenterology. 2021;160(3):781–96.

    Article  CAS  PubMed  Google Scholar 

  28. Zhang MdJ, Zhang MdL, Yang MdY, Liu MdQ, Ma MdH, Huang MdA, et al. Polymorphonuclear-MDSCs facilitate tumor regrowth after radiation by suppressing CD8(+) T cells. Int J Radiat Oncol Biol Phys. 2021;109(5):1533–46.

    Article  Google Scholar 

  29. Guo B, Fu S, Zhang J, Liu B, Li Z. Targeting inflammasome/IL-1 pathways for cancer immunotherapy. Sci Rep. 2016;6:36107.

    Article  PubMed  PubMed Central  Google Scholar 

  30. van de Veerdonk FL, Netea MG, Dinarello CA, Joosten LA. Inflammasome activation and IL-1β and IL-18 processing during infection. Trends Immunol. 2011;32(3):110–6.

    Article  CAS  PubMed  Google Scholar 

  31. Li Y, Wang L, Pappan L, Galliher-Beckley A, Shi J. IL-1β promotes stemness and invasiveness of colon cancer cells through Zeb1 activation. Mol Cancer. 2012;11:87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Kaler P, Augenlicht L, Klampfer L. Macrophage-derived IL-1beta stimulates Wnt signaling and growth of colon cancer cells: a crosstalk interrupted by vitamin D3. Oncogene. 2009;28(44):3892–902.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Das S, Shapiro B, Vucic EA, Vogt S, Bar-Sagi D. Tumor cell-derived IL1β promotes Desmoplasia and immune suppression in pancreatic cancer. Cancer Res. 2020;80(5):1088–101.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Tannenbaum CS, Rayman PA, Pavicic PG, Kim JS, Wei W, Polefko A, et al. Mediators of inflammation-driven expansion, trafficking, and function of tumor-infiltrating MDSCs. Cancer Immunol Res. 2019;7(10):1687–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Kessenbrock K, Plaks V, Werb Z. Matrix metalloproteinases: regulators of the tumor microenvironment. Cell. 2010;141(1):52–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Bergers G, Brekken R, McMahon G, Vu TH, Itoh T, Tamaki K, et al. Matrix metalloproteinase-9 triggers the angiogenic switch during carcinogenesis. Nat Cell Biol. 2000;2(10):737–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Bruno A, Bassani B, D’Urso DG, Pitaku I, Cassinotti E, Pelosi G, et al. Angiogenin and the MMP9-TIMP2 axis are up-regulated in proangiogenic, decidual NK-like cells from patients with colorectal cancer. Faseb J. 2018;32(10):5365–77.

    Article  CAS  PubMed  Google Scholar 

  38. Chen H, Ye Y, Yang Y, Zhong M, Gu L, Han Z, et al. TIPE-mediated up-regulation of MMP-9 promotes colorectal cancer invasion and metastasis through MKK-3/p38/NF-κB pro-oncogenic signaling pathway. Signal Transduct Target Ther. 2020;5(1):163.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Shao L, Zhang B, Wang L, Wu L, Kan Q, Fan K. MMP-9-cleaved osteopontin isoform mediates tumor immune escape by inducing expansion of myeloid-derived suppressor cells. Biochem Biophys Res Commun. 2017;493(4):1478–84.

    Article  CAS  PubMed  Google Scholar 

  40. Chu D, Zhao Z, Zhou Y, Li Y, Li J, Zheng J, et al. Matrix metalloproteinase-9 is associated with relapse and prognosis of patients with colorectal cancer. Ann Surg Oncol. 2012;19(1):318–25.

    Article  PubMed  Google Scholar 

  41. Ye Y, Kuang X, Xie Z, Liang L, Zhang Z, Zhang Y, et al. Small-molecule MMP2/MMP9 inhibitor SB-3CT modulates tumor immune surveillance by regulating PD-L1. Genome Med. 2020;12(1):83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Li M, Xing S, Zhang H, Shang S, Li X, Ren B, et al. A matrix metalloproteinase inhibitor enhances anti-cytotoxic T lymphocyte antigen-4 antibody immunotherapy in breast cancer by reprogramming the tumor microenvironment. Oncol Rep. 2016;35(3):1329–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Liu F, Qin L, Liao Z, Song J, Yuan C, Liu Y, et al. Microenvironment characterization and multi-omics signatures related to prognosis and immunotherapy response of hepatocellular carcinoma. Exp Hematol Oncol. 2020;9:10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Du M, Wang Y, Liu Z, Wang L, Cao Z, Zhang C, et al. Effects of IL-1β on MMP-9 expression in Cementoblast-derived cell line and MMP-mediated degradation of Type I collagen. Inflammation. 2019;42(2):413–25.

    Article  CAS  PubMed  Google Scholar 

  45. Huang Q, Lan F, Wang X, Yu Y, Ouyang X, Zheng F, et al. IL-1β-induced activation of p38 promotes metastasis in gastric adenocarcinoma via upregulation of AP-1/c-fos, MMP2 and MMP9. Mol Cancer. 2014;13:18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Cheng CY, Kuo CT, Lin CC, Hsieh HL, Yang CM. IL-1beta induces expression of matrix metalloproteinase-9 and cell migration via a c-Src-dependent, growth factor receptor transactivation in A549 cells. Br J Pharmacol. 2010;160(7):1595–610.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Abrams SI, Waight JD. Identification of a G-CSF-Granulocytic MDSC axis that promotes tumor progression. Oncoimmunology. 2012;1(4):550–1.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


Not applicable.


This subject was funded by grants from the Province Innovation Team Technology Project of Yunnan (No. 2019HC009), Provincial Colorectal Cancer Center Construction Project of Yunnan (No. 2020YZ001), Scientific Research Fund of Yunnan Provincial Education Department (No. 2022J0227) and Province Science and Technology Planning Project of Yunnan (No. 202201AY070001-149).

Author information

Authors and Affiliations



YL and LX2 designed the article. TW, XZ and XL drafted the manuscript and were responsible for analyzing the scRNA-seq data. XC was responsible for animal experiment. DP and RD were responsible for enrolling patients and collecting samples. RL and TS conducted immunohistochemistry and immunofluorescence. RH and JD modified the manuscript. FL and JL performed the colonoscopy. CW conducted the radiologic examination. LX1 and SG performed the statistical analysis. ZY completed the qPCR and western blotting. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Lu Xing or YunFeng Li.

Ethics declarations

Ethics approval and consent to participate

All participants had written informed consent before the study. Additionally, the Ethics Committee of Yunnan Cancer Hospital approved study protocols which were designed based on the medical research involving human subjects ethical principles of the Helsinki Declaration (approval number: KYLX202127). The mice applied in this study are all from the Experimental Animal Department of Kunming Medical University, and the experimental scheme has been approved by the Experimental Animal Ethics Committee of Kunming Medical University (approval number: kmmu20221065).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no no potential conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

 Table S1. Single cell sequencing data of clinical samples from 6 patients with MSI-H/dMMR.

Additional file 2:

Table S2. Markers used in cell subtype analysis pipeline.

Additional file 3: Table S3.

Number of cells of 9 cell subtypes in each sample.

Additional file 4:

Figure S1. GO and KEGG analysis of marker genes.

Additional file 5:

 Table S4. 1623 differentially expressed genes in 9 cell subsets.

Additional file 6:

 Table S5. 1454 pseudotime-related genes were identified from 1623 feature genes.

Additional file 7:

Figure S2. KEGG and GO analysis of pseudotime-related genes.

Additional file 8:

 Table S6. 155 DEGs in gene expression of each cell subtype between sensitive group and resistant group.

Additional file 9:

Figure S3. GO and KEGG analysis of PD-1 resistance-related DEGs.

Additional file 10:

 Table S7. 130 common genes among 454 pseudo-time related genes and 155 (de- duplication) PD-1 resistance related genes.

Additional file 11:

 Table S8. VarElect analyse of 130 common genes.

Additional file 12:

 Figure S4. VarElect analysis results of 130 common genes in 454 pseudotime-related genes and 155 (de-duplication) PD-1 resistance related genes DEGs. Figure S5. (A) The cytokine‒cytokine receptor interaction and (B) the MAPK and PI3K-Akt signaling pathway maps of IL-1β and MMP9. Figure S6. Construction of stable CT26 cell line overexpressing IL-1β. (A) Colorectal cancer cell line CT26 used in this experiment. (B) Plasmid map of IL-1β overexpression vector. (C) Screening positive colonies by ampicillin after plasmid transformation. (D) Gel map of plasmid electrophoresis. (E) Western blotting was applied to detect the expression of IL-1β in stably transformed cell lines. (F) qPCR was used to detect the expression of IL-1β in stably transformed cell lines.

Additional file 13.

 Images of the original blots.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wu, T., Zhang, X., Liu, X. et al. Single-cell sequencing reveals the immune microenvironment landscape related to anti-PD-1 resistance in metastatic colorectal cancer with high microsatellite instability. BMC Med 21, 161 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Metastatic colorectal cancer (mCRC)
  • Programmed cell death protein-1(PD-1)
  • Interleukin-1 beta (IL-1β)
  • Myeloid-derived suppressor cells (MDSCs)
  • microsatellite instability-high (MSI-H)
  • Deficient mismatch repair (dMMR)