Introduction

Immunotherapy with immune checkpoint inhibitors (ICI) targeting PD-1/PD-L1 has revolutionized the treatment for cancer (1). A pioneering clinical trial in 2015 investigating PD-1 inhibitor efficacy highlighted high microsatellite instability (MSI-H) as a pan-solid tumor biomarker (1). However, only approximately 5% of metastatic colorectal cancer are MSI-H, while the remaining 95% low microsatellite instability (MSI-L)/patients with microsatellite-stable (MSS) metastatic colorectal cancer barely respond to ICIs monotherapy (2, 3).

Recent advances identified the reasons for MSS colorectal cancer immunotherapy resistance: inadequate antigen release and presentation, poor immunogenicity, and an immune-cold or immune-excluded tumor microenvironment (TME; ref. 2). The abundant immunosuppressive cells, including tumor-associated macrophages (TAM), myeloid-derived suppressor cells (MDSC), regulatory T cells, and microenvironment factors like VEGFA and TGFβ (2), could exclude CD8+ T-cell infiltration and impair antitumor response. ICIs combined with VEGFR2 blockade has shown enhanced efficacy in MSS colorectal cancer preclinical models (4). Thus, elucidating crucial immune evasive factors in the TME may facilitate solving the obstacles of ICI treatment in MSS colorectal cancer.

Immunometabolism represents the interaction between metabolism and the immune TME (5). Metabolic flux in the TME, including deprivation of metabolic substrates, accumulation of metabolic waste and metabolism activity in various cell types, has crucial effects on the antitumor response (5). Immune regulatory metabolic genes, enzymes, and metabolites are recognized as “immunometabolic checkpoints” (6). Recently, targeting an immunometabolic checkpoint through glutamine enzyme inhibition was found to suppress glutamine circuits of tumor cells while upregulating oxidative metabolism in effector T cells (7). Metabolic antagonists may sensitize the TME to immunotherapy by triggering divergent metabolic programs that amplify the function of immune cells and reshape the immune TME with low toxicity and high efficiency (6). At present, immunometabolism strategies, such as glucose metabolism inhibitor Metformin, glutamine pathway inhibitor CB-839, adenosine inhibitor oleclumab, and IDO1/IDO2 inhibitors, combined with anti–PD-1/PD-L1 are under investigation in preclinical or clinical research (5, 8). Metabolism interventions have great potential in tackling the difficulties of poor MSS colorectal cancer immunotherapy response.

We identified a metabolic-immune subtype with the worst prognosis in MSS colorectal cancer, characterized by abundant stroma and high-level immune evasion. Furthermore, we discovered chondroitin sulfate (CS), specifically chondroitin-6-sulfate (C-6-S), as the key factor of poor prognosis. CS, a type of glycosaminoglycans (GAG), is a linear acidic polysaccharide comprised of repeating disaccharides (9). CS has shown anti-inflammatory activity in osteoarthritis and proinvasive effect in tumors (9, 10). However, the role of CS in the TME and immunotherapy is unclear. Here, we targeted C-6-S as an immunometabolic checkpoint through a metabolism-immune microenvironment exploration workflow, which may contribute to improving MSS colorectal cancer ICIs strategies.

Materials and MethodsHuman and mouse cell lines

We obtained murine colon carcinoma cell line CT26 from the National Collection of Authenticated Cell Cultures (NCACC). Cells were cultured in RPMI1640 medium (Solarbio, 31800) supplemented with 10% FBS (Solarbio, S9030). Murine MC38 colon cancer cell line was from BNCC and cultured in DMEM (Solarbio, 12100) with 10% FBS (Solarbio, S9030). Human monocytic cell line THP-1 was from NCACC and cultured with RPMI1640 medium (Solarbio, 31800) supplemented with 10% FBS (Solarbio, S9030). When indicated, THP-1 monocytes were incubated with phorbol 12-myristate 13-acetate (PMA; 100 ng/mL; Selleck, S7791) for 24 hours to differentiate into macrophages. Murine macrophage cell line RAW264.7 was from NCACC and cultured in DMEM (Solarbio, 12100) with 10% FBS (Solarbio, S9030). All cell lines were cultured at 37°C in a humidified atmosphere with 5% CO2. Authentication for all cell lines was acquired from the manufacturers and reauthenticated with short tandem repeat test within the past year. Cells were cultured for approximately 3 months for experiments and were regularly tested for Mycoplasma with PCR method. No more than 10 passages for THP-1 and RAW264.7 cells and no more than 20 passages for CT26 and MC38 cells were used.

Study cohorts

We performed transcriptomic analysis in The Cancer Genome Atlas (TCGA) colorectal cancer cohort. The cohort was filtered with tissue type (primary tumor). MSI-L/MSS samples are similar to proficient mismatch repair (pMMR) status, while MSI-H are similar to deficient mismatch repair status (11). After filtering samples with MSI status (MSI-L/MSS), 383 patients were selected and the cohort is referred to TCGA MSS colorectal cancer cohort (Supplementary Table S1). We assembled a clinical cohort with 112 colorectal cancer patient's paraffin-embedded surgical tumor specimens from Nanfang Hospital (Guangzhou, P.R. China; n = 35), Guangzhou First People's Hospital (Guangzhou, P.R. China; n = 74), and Zhujiang Hospital (Guangzhou, P.R. China; n = 3). Patients were diagnosed with colorectal cancer and received surgery in 2014–2018. Signed informed consents were obtained from all the patients for the use of clinical information and tissue samples. Paraffin-embedded tissue sections from primary tumors were collected and stored at 4°C until analysis. Tissue sections were observed under microscopy and tumor area was calculated as tumor cells–infiltrated area/total area × 100%. Tissue sections with more than 20% tumor area were qualified to perform further analysis. MMR status was defined on the basis of MSH2, MSH6, PMS2, and MLH1 protein expression via immunohistochemistry (IHC). pMMR status was considered as MSS, according to previous report (11). A total of 96 patients with pMMR were included for further investigation, among which 30 patients were from Nanfang Hospital (Guangzhou, P.R. China), 63 from Guangzhou First People's Hospital (Guangzhou, P.R. China), and 3 from Zhujiang Hospital (Guangzhou, P.R. China). The study was conducted in accordance with the Declaration of Helsinki. The use of human tissue samples and clinical data was approved by the ethics committee of Nanfang Hospital (Guangzhou, P.R. China).

DNA MMR status classification

DNA MMR system is constituted by four MMR genes and their proteins (MLH1, MSH2, MSH6, and PMS2). For collected colorectal cancer samples, IHC of the four MMR proteins was performed as described in the IHC assays. One pathologist blind to data analysis interpreted the results. Normal colonic crypt epithelium adjacent to tumor, lymphoid cells, and stroma cells were used as internal positive control. Any positive staining of tumor cells is considered positive. If the expression of all four MMR indicators is normal, the sample is classified as proficient MMR/MSS; if the expression of one or more indicators is missing, then the sample is classified as deficient MMR/MSI-H.

Calculation of metabolism and immune characteristics score

The metabolism pathways were obtained from Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/; ref. 12) and curated on the basis of a previous report (13). Metabolism-related pathways with gene set size of 3–200 were selected, resulting in 113 metabolism pathways with 1,784 metabolic genes. GAGs metabolism could be divided into CS, heparan sulfate, and keratan sulfate metabolism according to the KEGG database. Metabolism pathway gene sets, including subsets from GAGs metabolism, are provided in Supplementary Table S2. Considering cell type–specific expression signature and dominant cell types in colorectal cancer, we applied curated immune and microenvironment cell gene sets of 25 cell types reported by previous study (14), which were filtered from CIBERSORT and MCP-counter. We integrated gene sets from Msigdb (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) and previous research to generate TME gene sets [including hypoxia, reactive oxygen species, angiogenesis, lymphangiogenesis (15) and stroma (16) pathways] and immune phenotype regulatory factor gene sets including MHC machinery BioCarta, IFNG signaling (BioCarta), CD8 T effector, immune checkpoints (17), cytokines related to immune response (Reactome), TGFβ (17), angiogenesis, stroma, and DNA replication (Gene Ontology) pathways. In addition, immune response and evasion-related mechanic pathways (18), including nonredundant 24 HALLMARK pathways, 21 IPA pathways (http://www.ingenuity.com), adenosine, immungentic cell death, and NOS1 pathways were obtained. Tumor molecular function portrait (19), including oncogenic pathways (20), pro- and antitumor microenvironmental factors and stromal network, based on single-sample gene set enrichment analysis (ssGSEA) scores was performed. Transcripts per kilobase million (TPM) values were used, unless specified, to perform ssGSEA via GSVA R package (21). Enrichment scores of each pathway for each patient were calculated for further investigation.

Clustering of metabolism subtypes

We used K-means clustering algorithms and the ConsensusClusterPlus R package (22) to identify metabolism subtypes. To achieve stable clustering, we performed 1,000 iterations with 80% resampling (k = 2–6). Cluster stability was evaluated with cumulative distribution function and delta area plot. k = 3 was chosen on the basis of clustering stability and clinical significance. Samples were reordered according to K-means clustering results. Metabolism pathway scores were scaled before heatmap graphing using ComplexHeatmap R package (23).

Metabolism subtype–specific pathway and category activation ratio

To select subtype-specific pathways, we first applied Boruta algorithms based on a random forest model. With metabolism subtype as classification indicator, featured pathway selection was carried out using Boruta R package (24) (maxRuns = 200). Comparisons of metabolism pathway scores [stroma metabolism (SM) subtype vs. Others, nucleotide metabolism (NM) subtype vs. Others, energy metabolism (EM) subtype vs. Others] were conducted using the Limma R package (25). log2 (fold change) > 0 and FDR < 0.05 was considered significant upregulation. The top 10 significantly upregulated metabolism pathways that passed the Boruta importance test in each subtype were defined as subtype-specific metabolism pathways (3 for SM subtype, 9 for NM subtype, and 10 for EM subtype). On the basis of information in the KEGG database, pathways were classified into 11 metabolism categories (Supplementary Table S2). Categories with three or more pathways (n = 9) were used in metabolism category activation ratio analysis. Metabolism category activation ratios were calculated as the proportion of subtype-specific upregulated pathways in each metabolism category for each subtype.

Immune response and evasion score

T cell–inflamed gene expression profile (GEP; ref. 26) reflects IFNγ signaling related antitumor immune response. High GEP was reported to predict immunotherapy efficacy in melanoma and non–small cell lung cancer (26). In this study, we calculated the GEP score as the average expression of 18 genes in the signature using TPM values with log2 transformation. T-cell dysfunction and exclusion (TIDE) score (27) describes immune evasion; high TIDE scores predict immunotherapy resistance. Standardized data were generated using the gene expression median as the normalization control and uploaded to obtain TIDE score for tumor samples (http://tide.dfci.harvard.edu/log/). Cytotoxicity score was defined by the average expression values of GZMA and PRF1 (28).

Digital pathology analysis

Image analysis processing included the following steps: (i) Two regions of interest (ROI) in each histologic slide were manually outlined, namely the tumor core (TC) and the invasive margin (IM). The TC is the main area of the tumor tissue, and the IM is a 500 μm wide band-shaped area at the interface between tumor and non-tumor tissue, according to previous study (29). (ii) Open-source software QuPath (v0.2.2, https://github.com/qupath/qupath; ref. 30) was used to automatically detect all cell structures in the ROI, and the staining intensity threshold was set to recognize positive-staining cell. In general, dark brown is considered as strong positive, brown-yellow as moderate positive, light yellow as weak positive, and blue nuclei as negative. The detection results were manually checked. (iii) Output the density of positive cells (positive cell number per mm2) or H-Score [1× percentage of weak positive cells (%) + 2× percentage of moderate positive cells (%) + 3× percentage of strong positive cells (%)] for further analysis. (iv) Use QuPath's built-in random forest cell recognition module when necessary, and train the random forest model to identify cells in the ROI, using annotated tumor cells and stromal cells regions as input. Each slide was manually checked and a suitable classification model was applied. The steps above were supervised by two pathologists. As for immune phenotype, high CD8+ TC and high CD8+ IM is described as “hot”; low CD8+ TC and low CD8+ IM is “cold,” and high CD8+ IM but low CD8+ TC is “excluded” phenotype. A total of 70th percentile of CD8+ cell density was used as the cutoff to determine high/low grouping. The digitized tissue slides at 40× magnification were acquired using Aperio GT450 digital pathology scanner and analysis system (Leica). Because of automated focusing errors, insufficient staining and limited clinical specimen tissues, the tissue sections that failed to produce qualified digital slides were marked as NA values.

IHC assays

Tissue samples were collected and formalin fixed and paraffin embedded. A total of 4-μm tissue sections were used for immunostaining as reported previously (31). Briefly, after antigen retrieval, samples were blocked with BSA (Solarbio, 9048-46-8) for 1 hour, and then incubated with primary antibodies overnight at 4°C. Tissue sections were then washed in PBS and incubated with horseradish peroxidase–labeled goat anti-rabbit IgG antibody (Bioss, bs-0295G) for 1 hour at room temperature. 3,3′-Diaminobenzidine staining and hematoxylin counterstaining were applied and mounted. Primary antibodies were used as follows: C-6-S (Millipore, MAB2035, RRID: AB_11214309, 1:300), CD8 (Proteintech, 66868-1-Ig, RRID: AB_2882205, 1:3,000), CD163 (Proteintech, 16646-1-AP, RRID: AB_2756528, 1:800), CD206 (Proteintech, 60143-1-Ig, RRID: AB_2144924, 1:8,000), αSMA (Proteintech, 14395-1-AP, RRID: AB_2756528, 1:8,000), GLI1 (Proteintech, 66905-1-Ig, RRID: AB_2882232, 1:1,000), pSTAT3 (Cell Signaling Technology, 9145, RRID: AB_2491009, 1:800), MLH1 (Proteintech, 11697-1-AP, RRID: AB_2145604, 1:200), MSH2 (Proteintech, 60161-1-Ig, RRID: AB_10666855, 1:200), MSH6 (Proteintech, 66172-1-Ig, RRID: AB_2881567, 1:400), PMS2 (Proteintech, 66075-1-Ig, RRID: AB_11182595, 1:400).

Primary cultured fibroblasts

Fresh human colorectal cancer specimens and the adjacent non-tumor tissues were collected to isolate cancer-associated factors (CAF) and NFs, respectively. Murine CAFs were isolated from CT26 subcutaneous tumors. Tumor or non-tumor tissues were diced into approximately 1 mm3 with a razor blade and were digested with a solution of 5 mg/mL DNase I (Meilunbio, MB3069) + 20 mg/mL collagenase IV (Biofroxx, 2091MG100) + 20 mg/mL hyaluronidase (Solarbio, H8030) for 3–4 hours at 37°C on a rotating platform. Then cells were resuspended, filtered through a cell strainer (75 μm), and rinsed with PBS. Cells were collected, seeded into 24-well plates, and were cultured in DMEM (Solarbio, 12100) supplemented with 10% FBS (Solarbio, S9030) and 1% antibiotic-antimycotic (Solarbio, P8420). Primary fibroblasts were purified from other cell populations by differential adhesion and serial passage. Cell identity was confirmed by αSMA immunofluorescence (IF) staining. All primary fibroblasts in experiments were below 10 passages.

Primary mouse bone marrow–derived macrophages, mouse peritoneal macrophages, and human colorectal cancer TAMs

Six- to 10-week-old C57BL/6 mice were euthanatized to obtain femurs. The femurs were dissected with scissors and the muscles attached to the bone were removed. Bone marrow was flushed out using DMEM (Solarbio, 12100) and placed into sterile tube. The bone marrow was then homogenized with plastic pipette and the primary progenitor cell suspension was generated. Then the cell suspensions were incubated in DMEM (Solarbio, 12100) with addition of 10% FBS (Solarbio, S9030), 100 U/mL penicillin-streptomycin (Solarbio, P1400) and 10% L929 cell-conditioned medium, for 7 days. L929 cell-conditioned medium was generated through cultivation of L929 cells in RPMI1640 medium (Solarbio, 31800) with 10% FBS (Solarbio, S9030) for 10 days, which contains macrophage colony-stimulating factor. Differentiation of macrophages was determined using anti-F4/80 (MultiSciences, 70-AM048010-20, Clone:BM8.1) with flow cytometry as described previously (32). Primary peritoneal macrophages (PM) were obtained by flushing the mouse peritoneal cavity with PBS. Cells were washed with PBS twice and cultured in DMEM (Solarbio, 12100) for 1 hour at 37°C and 5% CO2. PBS was then used to remove nonadherent cells. The purity of macrophages was tested with anti-F4/80 (MultiSciences, 70-AM048010-20, Clone:BM8.1) by flow cytometry (purity >80%). Human TAMs were isolated from fresh colorectal cancer tissues as described previously (33). After digestion, cell suspensions were obtained and placed in a 15 mL tube with 5 mL 45% Percoll (Solarbio, P8370) in the middle and 5 mL 60% at the bottom, which were then centrifuged at 800 × g for 30 minutes. Cells from the interphase were then isolated with CD14 isolation kit (Miltenyi Biotec, 130-097-052) based on the manufacturer's instructions.

Three-dimensional cocultivation assays

Collagen type I (3 mg/mL; Solarbio, C8062), 1 mol/L NaOH solution (Acmec, S41251), 10× PBS (Servicebio, G0002) solution, and dH2O were mixed in a volume ratio of 40:1:6:13 to prepare a pH neutral collagen working solution. Aliquots of 300 μL were used to prepare three-dimensional (3D) collagen. After the primary CAF cells and macrophages were counted, 3.0 × 105 cells of each were mixed in equal proportions, dissolved in 3D collagen, and placed in 37°C and 5% CO2 atmosphere. After gelation, the 3D collagen system was covered with culture medium to construct a 3D coculture model. IF was used to detect the expression and localization of related indicators.

Coculture system for CAFs and macrophages

In the co-culture model, fibroblasts and macrophages were cultured in a chamber (JET Biofil, TCS016012). A total of 2 × 105 adherent fibroblasts were added to the upper layer, and 2 × 105 macrophages were cultured in the lower layer. Co-cultivation lasted for continuous 48 hours. The following reagents were added as indicated: Surfen (Millipore, S6951, 20 μmol/L), chondroitinase ABC (Ch-ABC; Yuanye Bio-Technology, S31309, 0.2 U), Stattic (Selleck, S7024, 2.5 μmol/L), and Vismodegib (Targetmol, T2590, 0.05 μmol/L). The lower layer of macrophages was collected for RNA extraction, and M1 and M2 phenotype polarization indicators were detected via qRT-PCR and flow cytometry.

IF assays

The IF was performed as described previously (31). The samples were incubated with the primary antibody at 4°C overnight. Then the secondary antibody Alexa fluor 647-labeled goat anti-mouse IgG(H+L) (Beyotime, A0473, RRID: AB_2891322, 1:500), Alexa Fluor 488-labeled goat anti-rabbit IgG(H+L) (Beyotime, A0423, RRID: AB_2891323, 1:500), Cy3-labeled goat anti-rat IgG(H+L) (Beyotime, A0507, 1:500) or Alexa Fluor 594 AffiniPure donkey anti-goat IgG (H+L) (Yeasen, 34312ES60, 1:200) was used for incubation for 1 hour. Next, the samples were incubated with the methanol dilution of DAPI (Beyotime, C1002, 1:1,000) at room temperature for 5 minutes. The images were taken with a fluorescence and laser confocal microscope (Nikon ECLIPSE Ti2) and were analyzed using ImageJ software (34). Radial fluorescence intensity analysis was performed using Plot Profile, a plug-in of ImageJ software, to observe fluorescence colocalization. The following antibody concentrations were used: C-6-S (Millipore, MAB2035, RRID: AB_11214309, 1:200), EPCAM (Proteintech, 21050-1-AP, RRID: AB_10693684, 1:100), αSMA (Novus Biologicals, NB300-978SS, 1:200), CD163 (Proteintech, 16646-1-AP, RRID: AB_2144924, 1:200), CD8 (Novus Biologicals, NB200-578, RRID: AB_10003082, 1:200), CD206 (Cell Signaling Technology, 91992, RRID: AB_2800175, 1:200), GLI1 (Proteintech, 66905-1-Ig, RRID: AB_2882232, 1:200), pSTAT3 (Cell Signaling Technology, 9145, RRID: AB_2491009, 1:200).

qRT-PCR

According to the manufacturer's instructions, total RNA from tissue samples was extracted with TRIzol reagent (Invitrogen). The reverse transcription kit HiScript II Q RT SuperMix for qPCR (R222-01, Vazyme, Nanjing) was used to synthesize cDNA from total RNA. Real-time qPCR was performed using the LightCycler 480 system Version 1.5 (Roche). The indicator gene expression was scaled using GAPDH expression as control. The 2–ΔΔCt method was used to calculate the expression fold change. Each qRT-PCR experiment was independently repeated in triplicate. The primer sequences used for qRT-PCR in this study are shown in the Supplementary Table S3.

Western blotting and coimmunoprecipitation

Western blotting (WB) was performed as described previously (31). Immunoblots were detected with fluorophore-conjugated goat anti-rabbit or anti-mouse secondary antibodies by an Odyssey imaging system (LI-COR). Antibodies used were as follows: GLI1 (Proteintech, 66905-1-Ig, RRID: AB_2882232, 1:1,000), pSTAT3 (Cell Signaling Technology, 9145, RRID: AB_2491009, 1:1,000). For the immunoprecipitation (IP) experiment, FLAG-GLI plasmid (GeneCopoeia, EX-F0407-Lv242) was transiently transfected into macrophages using Lipofectamine 2000 (Invitrogen) as described previously (35). Following 24 hours, the medium was replaced and treated with C-6-S (1 mg/mL for RAW264.7 and 2.5 mg/mL for THP-1) for another 24 hours. FLAG tag antibody (Cell Signaling Technology, 14793, AB_2572291) was used to pull down FLAG-GLI1 protein and its binding protein complex. pSTAT3 antibody (Cell Signaling Technology, 9145, RRID: AB_2491009) was used to pull down the complex of its binding protein to determine the interaction between the two transcription factors. Briefly, after rinsing with cold PBS, cells were lysed in IP lysis buffer (Meilunbio, MB9900). Then, Protein A+G Sepharose Beads (7Sea biotech) and the primary IP antibody were added and placed on a vibration platform overnight at 4°C to precipitate the immune complexes. Samples were rinsed with IP lysis buffer five times, boiled and eluted in SDS-PAGE buffer (LEAGENE, PE0025) to carry out further WB analysis. Equal amounts of protein were electrophoresed on SDS-PAGE and then transferred to nitrocellulose membrane. The immunoblots were blocked with 5% skim milk powder and detected by using enhanced chemiluminescence reagent. Antibodies were used as follows: GLI1 (Proteintech, 66905-1-Ig, RRID: AB_2882232) pulldown (1:50), WB (1:1,000); pSTAT3 (Cell Signaling Technology, 9145, RRID: AB_2491009) pulldown (1:50), WB (1:1,000); FLAG (Cell Signaling Technology, 14793, RRID: AB_2572291) pulldown (1:50), WB (1: 1,000).

In vitro and in vivo flow cytometry analysis

In the in vitro experiments, cells were incubated with the antibody conjugated with fluorescence in 100 μL of FACS staining buffer [1×PBS (Servicebio, G0002) containing 1% BSA (Solarbio, 9048-46-8)] and were protected from light and incubated for 30 minutes. In the in vivo flow cytometry analysis, the subcutaneous tumor was removed and then mechanically separated with scissors in sterile PBS. The tumor tissues were passed through a 75 μm cell strainer to obtain a single-cell suspension. After resuspending in PBS containing 0.5% BSA, mouse tumor tissue monocyte separation medium kit (Solarbio, P3970) was used to extract monocyte macrophages according to the manufacturer's instructions. Cells were incubated with appropriate antibodies used for cell labeling for 30 minutes. Use the following antibodies: PE rat monoclonal CD206 antibody (eBioscience, 12-2061-82, MR6F3, RRID: AB_2637421, 0.125 μg/test), Super Bright 436 rat monoclonal PD-L1 antibody (eBioscience, 62-5982-80, MIH5, RRID: AB_2637417, 0.25 μg/test), PE-Cy7 mouse monoclonal F4/80 antibody (MultiSciences, 70-AM048010-20, Clone:BM8.1, 0.25 μg/test), APC mouse monoclonal CD11c antibody (MultiSciences, 70-AM011C05-20, 0.125 μg/test), FITC anti-human CD14 antibody (BioLegend, 301803, M5E2, RRID: AB_314185, 5 μL per million cells), PerCP/Cyanine5.5 anti-human CD86 antibody (BioLegend, 305419, IT2.2, RRID: AB_1575070, 5 μL per million cells), PE anti-human CD163 antibody (BioLegend, 333605, GHI/61, RRID: AB_1134005, 5 μL per million cells). A FACS Aria II (BD Biosciences) was used to detect fluorescence, and FlowJo X (v10.6.2) software was used to analyze the data.

In vivo mouse studies

All animal experiments were performed in accordance with the protocol approved by the Ethics Committee of Nanfang Hospital of Southern Medical University (Guangzhou, P.R. China). Female BALB/c mice (6–7 weeks of age) and C57BL/6 mice (7–8 weeks of age) were obtained from the Experimental Animal Center, Nanfang Hospital (Guangzhou, P.R. China), Southern Medical University (Guangzhou, P.R. China). The mice were maintained at 22°C–24°C temperature, 60 ± 10% humidity, with the 12-hour light/dark cycle, under pathogen-free conditions. Standard rodent laboratory diet and water ad libitum were provided. Murine CT26 colorectal cancer cell line was obtained from the NCACC (Shanghai, P.R. China), which was reported to be MSS (36). A total of 1 × 106 CT26 cells were injected into the right thigh of the mice (marked as day 0), and the tumor size was measured every other day with calipers. Using the following formula to calculate the tumor volume (mm3): tumor volume = (π)/6 × L × W, where L is the long axis size and W is the vertical size. Seven days after tumor implantation, the mice were randomly divided into experimental groups. The following 2-week treatments were given: Surfen (Selleck, S6951), 1 mg/kg, intraperitoneal injection, days 7–13 for 7 consecutive days; anti–PD-1 (BioXcell, BE0146), 12.5 mg/kg, i.p., twice a week, that is, on days 10, 13, 17, 20; mFOLFOX6, oxaliplatin (Sanofi) 6 mg/kg followed 2 hours later by 5-fluorouracil (Xudong) 5 mg/kg and leucovorin (Yaoyou) 90 mg/kg, all injected intraperitoneally, once a week, that is on day 8 and day 15; Regorafenib (Bayer, BAY73-4506), dissolved in 0.5% methylcellulose (aladin, 9004-65-3), administered orally, daily for 2 weeks. In addition, azoxymethane (Sigma-Aldrich, A5486)/dextran sulfate sodium (MP Biomedicals, 9011-18-1) [azoxymethane/dextran sulfate sodium (AOM/DSS)] induced colorectal cancer mice model were established as reported previously (37). Briefly, C57BL/6 mice were treated with a single intraperitoneal injection of 10 mg/kg AOM, followed by 6 days 2% DSS drinking water and then 16 days normal water for three cycles. After that, mice were randomized into experimental groups and given a 2-week treatment. Surfen or anti–PD-1 was applied when indicated (Surfen: 1 mg/kg, i.p., everyday; anti–PD-1: 12.5 mg/kg, i.p., twice a week). When the treatment completed, mice were then sacrificed for analysis. The investigators were blind to the treatment groups during the experiment and outcome assessment. The mice were monitored daily and euthanized by cervical dislocation when showing any sign of discomfort.

Statistical analysis

The χ2 test was used to test the relationship between clinical information and metabolic subtypes. The Kolmogorov–Smirnov test was used to confirm whether the data follow a normal distribution. Independent-sample t tests, paired-sample t tests, one-way ANOVA, Wilcoxon tests, and Kruskal–Wallis tests were used to compare continuous variables where appropriate. Correlation coefficients were calculated using Pearson and Spearman rank test. The GSVA R package (21) was used to generate ssGSEA scores for indicated gene set in each patient. Wald statistical test and generalized linear model was used for differential gene analysis via the DESeq2 R package (38). Differential ssGSEA scores were analyzed using the Limma package (25). The Kaplan–Meier method was used for survival analysis to generate survival curves, and the log-rank test was applied to determine the statistical significance. A univariate and multivariate Cox proportional hazard regression model was used to calculate the HR and 95% confidence intervals. Multivariate logistic regression was performed with treatment response as a binary outcome. When indicated, the survminer package was used to determine the optimal threshold based on the maximum rank statistic, which was then used to divide patients into high and low expression groups. R packages ggplot2, ComplexHeatmap (23), clusterProfiler (39), and survival were used to analyze data and generate plots. Gene set enrichment analysis (GSEA) java software was used for gene enrichment analysis. The Benjamini-Hochberg test was used to adjust the P value to reduce the false positive rate. All tests are bilateral, and P < 0.05 is considered significant. Statistical tests were performed using R software (version 3.6.2, http://www.R-project.org) or GraphPad Prism 8.0.

Data and code availability

The datasets used in the current study, including TCGA cancer cohorts, GSE39582 (40), GSE17536 (41), GSE33113 (42), GSE81861 (43), and GSE35602 (44), are available in TCGA database (tcga-data.nci.nih.gov/tcga) or the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo). Datasets referring to cohorts receiving ICI treatment were acquired through public database or appropriate request to authors of previous research. The urothelial cancer cohort (17), receiving anti–PD-L1 drugs (atezolizumab), can be downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies. The data and code that support the findings of this study are included in the article and Supplementary Data.

ResultsComprehensive transcriptomic analyses identify an SM subtype with poor prognosis in MSS colorectal cancer

To systematically interpret the interaction of metabolism and the immune TME, we profiled the transcriptome of 14 types of solid tumors in TCGA database. There was a general correlation between metabolic pathways activation and immune infiltration (Supplementary Fig. S1A; Supplementary Tables S1, S2, and S4). Within colorectal cancer, the MSI-L/MSS subgroup had a higher proportion of metabolic-immune correlation (MSI-L/MSS vs. MSI-H: 54% vs. 18.7%; Supplementary Fig. S1A), indicating key roles of metabolic factors in regulating the immune TME.

To unravel metabolic heterogeneity, we performed unsupervised clustering of TCGA MSS colorectal cancer cohort and identified three clusters with distinct metabolism characteristics (Fig. 1A; Supplementary Fig. S1B–S1D). On the basis of importance score from the Boruta algorithm (24) and subtype-specific expression, 22 feature pathways were profiled. Cluster 1 displayed high GAGs biosynthesis, cyclooxygenase arachidonic acid metabolism and prostaglandin biosynthesis; Cluster 2 highly expressed pyrimidine synthesis and metabolism; Cluster 3 was characterized by nitrogen metabolism and urea cycle pathways (Fig. 1A; Supplementary Table S5). Metabolic activation ratio analysis revealed a prominent activation (100%) of NM in Cluster 2 and EM and xenobiotic biodegradation in Cluster 3 (Fig. 1B). Thus, we named Cluster 2 and Cluster 3 as NM subtype and EM subtype, respectively. In contrast, Cluster 1 exhibited minimum activation (≤10%) in any of metabolic categories (Fig. 1B). We found Cluster 1 overlapped with the subtypes with abundant stroma infiltration and poor prognosis in published CCMS, CRCA, CCS, Stroma Contribution Subtype, and CMS (45), suggesting that Cluster 1 was associated with stroma remolding (Supplementary Fig. S1E; Supplementary Table S6). Moreover, the feature metabolic pathways of Cluster 1 were more highly enriched in colorectal cancer stroma than epithelial tissue [GSE35602 (44)], based on which we named it SM subtype (Supplementary Fig. S1F). Metabolic subtypes showed significant prognosis difference, among which SM subtype was the worst and EM subtype was better [disease-free survival (DFS): EM vs. SM, P = 0.0056; overall survival (OS): EM vs. NM, P = 0.016; OS: EM vs. SM, P = 0.0022; Fig. 1C]. The reproducibility of our clustering and prognosis value were validated by external cohorts [GSE39582 (40), GSE17536 (41); Supplementary Fig. S1G–S1I].

Figure 1.

Comprehensive transcriptomic metabolism pathway-based clustering in MSS colorectal cancer. A, Metabolic pathway scores and unsupervised k-means clustering were performed in TCGA colorectal cancer cohort (n = 472). Patients with MSS colorectal cancer (n = 383) were filtered and three metabolic subtypes are shown with the specifically upregulated metabolic pathways. B, The proportions (%) of significantly upregulated pathways (log2 FC > 0, FDR < 0.05) in each metabolism categories among metabolic subtypes (a, Carbohydrate metabolism; b, Lipid metabolism; c, Amino acid metabolism; d, NM; e, EM; f, Metabolism of other amino acids; g, Glycan biosynthesis and metabolism; h, Xenobiotics biodegradation and metabolism; i, Metabolism of cofactors and vitamins). C, Kaplan–Meier plot showing DFS (left) and OS (right) among metabolic subtypes. P value was evaluated by log-rank test. D, Molecular function portrait (including oncogenic pathways, pro- and anti-TME factors and stromal network) based on ssGSEA score difference among subtypes (EM vs. Others, NM vs. Others, SM vs. Others). Antitumor effect is marked with blue and protumor factors with red, and the absolute value of log2 FC was shown with color and size changes. **, P < 0.01; *, P < 0.05.

In terms of clinicopathologic characteristics, there were no significant differences in sex, tumor site, stage, or KRAS, BRAF and EGFR mutations among the three subtypes (Supplementary Fig. S2A). Furthermore, we applied molecular function portrait (ref. 19; including oncogenic pathways, protumor and antitumor microenvironmental factors and stromal network) to characterize the three subtypes. Most oncogenic pathways were downregulated in the EM subtype, while cell cycle, TP53, and proliferation scores were upregulated in the NM subtype (Fig. 1D). Multiple oncogenic pathways, including Hippo, NOTCH, and Wnt pathways, were amplified in the SM subtype, indicating a high degree of malignancy (Fig. 1D). The SM subtype also featured simultaneous activation of both protumor and antitumor microenvironment factors, along with stromal network dysregulation (Fig. 1D). GSEA showed that TNFα inflammatory and angiogenesis pathways were enriched in the SM subtype (Supplementary Fig. S2B and S2C). Taken together, we identified distinct metabolism subtypes in MSS colorectal cancer, among which the SM subtype, with distinct TME and stromal remolding, had the worst prognosis and requires further investigation.

CS metabolism and immune microenvironment cross-talk in the MSS colorectal cancer SM subtype

Potential immunotherapy response can be marked by a high cytolytic activity (CYT) score (represented by PRF1 and GZMA) and PD-L1. We observed 65.1% of SM subtype patients highly expressed CYT and PD-L1, implying the SM subtype could potentially benefit from immunotherapy despite the poor prognosis (Fig. 2A). We next applied GEP and TIDE scores (26, 27) and found that the SM subtype exhibited high GEP and TIDE score simultaneously (Fig. 2B). Both antitumor immune response and immune evasion factors may be activated in the SM subtype.

Figure 2.

Chondroitin sulfate metabolism and immune microenvironment crosstalk in MSS colorectal cancer SM subtype. A, Classification of patients with TCGA MSS colorectal cancer based on median CYT score and CD274 (PD-L1) expression (log2 (TPM+1)). The proportions (%) of the four groups in each metabolic subtype were shown on the right. I, CYT-high, CD274-high; II, CYT-low, CD274-high; III, CYT-low, CD274-low; IV, CYT-high, CD274-low. B, Violin plots presenting GEP (left) and TIDE (right) scores among metabolic subtypes. The Wilcoxon test was performed to assess significance. C, Radar chart presenting the ssGSEA score of immune regulatory factors. Scores were rescaled into 0%–100% with simple linear conversion. D, Differential analysis of infiltration scores of immune and TME cells (n = 25) between metabolic subtypes (SM vs. Others, NM vs. Others, EM vs. Others). E, Univariate Cox regression analysis for OS of metabolic pathways, including 25 TME cell types, 10 oncogenic pathways and 5 TME factors. F, Volcano plot, the X axis represents the Pearson correlation coefficient between factors and TIDE score. G, Top 10 metabolic pathways with highest HR in univariate Cox model for OS and their correlation with immune regulatory factors. ****, P < 0.0001; ***, P < 0.001; **, P < 0.01; *, P < 0.05; ns, P > 0.05.

In general, successful antitumor immune response requires seven steps in the cancer-immunity cycle (27, 46), in which immune evasion could occur at any step. We assessed seven regulatory factors: (i) antigen presentation molecules and costimulators; (ii) IFNγ signaling; (iii) cytolytic activity; (iv) immune checkpoint expression; (v) cytokine and chemokines; (vi) tumor proliferation marker Ki-67, cell cycle and DNA replication; (vii) TME TGFβ signaling, stroma remodeling and angiogenesis. Most MHC molecules and costimulatory molecules were upregulated in EM and SM subtypes (Supplementary Fig. S3A and S3B). IFNγ signaling, CYT and most immune checkpoints were elevated in the SM subtype (Supplementary Fig. S3C–S3E). As for cytokine and chemokines, the EM subtype highly expressed CCL28 (Supplementary Fig. S3F), which is associated with transportation of lymphocytes (47). We noticed CSF1/CSF1R axis and CCL2 were enriched in the SM subtype, indicating the relevance to macrophage recruitment and M2 polarization (Supplementary Fig. S3G and S3H; ref. 48). Immune exclusion-related TME factors also enhanced in the SM subtype included TGFβ signaling, stroma remodeling, and angiogenesis score (Fig. 2C; Supplementary Fig. S3I). In addition, we evaluated the variety and abundance of tumor-infiltrated immune cells in the TME. The EM subtype displayed increased infiltration of plasma cells, most immune cell infiltration was decreased in the NM subtype, and the SM subtype displayed observed an increase of stromal cells (CAFs, endothelial cells) and immunosuppressive myeloid cells (M2 macrophages, MDSCs), which is closely related to immune exclusion (Fig. 2D; Supplementary Fig. S3J; Supplementary Table S7). Therefore, we extended its definition to be the SM–immune excluded (SM-IE) subtype.

To determine key factors causing immune exclusion in the SM-IE subtype, we comprehensively evaluated oncogenic pathways, metabolic pathways, and TME cells and factors through univariate Cox model (top five factors with highest or lowest HR value selected; Fig. 2E). Correlation coefficient values with TIDE score were regarded as another parameter (Fig. 2F). We found that none of the TME factors or cells were significant in both tests, while the oncogenic NOTCH pathway was a risk factor for prognosis and positively correlated with TIDE score. In terms of metabolism, only CS metabolism gained significance, with a higher effect size compared with the NOTCH pathway (Fig. 2E and F; Supplementary Table S8). Multivariate Cox regression analysis confirmed that CS metabolism was a risk factor for OS, independent of age, stage, and oncogenic pathways, which was also supported by colorectal cancer external cohorts [GSE33113 (42); Supplementary Table S9].

We then evaluated the correlations between the top 10 metabolic pathways with highest HR and immune TME factors. Results demonstrated a close relationship between CS metabolism and M2 macrophages, CAFs, and immune exclusion-related TME factors (angiogenesis and stroma; Fig. 2G; Supplementary Table S10). Collectively, these results suggested immune exclusion in MSS colorectal cancer SM-IE subtype, and CS metabolism may be intimately intertwined with an immune suppressive TME.

An immune exclusion barrier is defined by C-6-S and M2 macrophages in the MSS colorectal cancer invasive margin

It is reported that CS subtype C-6-S, not chondroitin-4-sulfate (C-4-S), markedly increases in colorectal cancer (49). To determine spatial patterns of C-6-S and immune cells in the TME, we carried out digital pathology analysis on 96 patients with MSS colorectal cancer in a multicenter clinical cohort (see clinicopathologic characteristics in Supplementary Table S11). QuPath (30) was used for cell detection and cell intensity classifications (Supplementary Tables S12 and S13). The distribution of C-6-S, CD8+ T cells (CD8), and M2 macrophages (CD163) were analyzed in TC and IM (Fig. 3A). Interestingly, C-6-S+ IM (P = 0.008), but not C-6-S+ TC (P = 0.311), predicted poor prognosis of MSS colorectal cancer (Fig. 3B). IM is the boundary compartment between tumor and normal tissue, reshaped by chemokines, cytokines, and TME cells like CAFs and macrophages (29). These results suggested the interplay between C-6-S and immune cells located in the IM.

Figure 3.

An “exclusion barrier” was constructed by C-6-S and M2 macrophages in MSS C