Single-cell mRNA-sequencing data (Che et al., 2021) were collected from 6 CRC patients, patients numbered COL15, COL 17, and COL18 were patients who received chemotherapy and the others who were not. Bulk RNA-seq datasets GSE41258, GSE39582, GSE103479, GSE38832, GSE192667 and GSE15921 as well as TCGA bulk RNA-seq data with corresponding detailed clinical information were included in our analysis. This study adhered to the guidelines set by the TCGA and GEO databases.

scRNA-Seq data were processed by following Seurat pipelines in R (Hao et al., 2021). Briefly, genes expressed in less than 3 cells, as well as cells expressed less than 250 or more than 3000 genes were filtered out. Cells with high mitochondria and rRNA gene proportions were also excluded. Then, log-normalization was conducted to normalize the data from the 6 samples. The highly variable genes were identified using the FindVariableFeatures function, followed by scaling of all genes. PCA dimensionality reduction was performed to identify anchors. The cells were clustered with a resolution of 0.2. After initial sample integration, cell clustering and annotation, we generated a gene expression and phenotype matrix of 1897 CAF cells from all 111,292 cells.

Metabolic activities of CAFs were evaluated by SCmetabolism packages (Wu et al., 2022). Each CAF cell was scored using the VISION algorithm, and finally the activity scores of the cell in different metabolic pathways were obtained. Second round clustering and subtyping of CAFs cells based on SCmetabolism scores were conducted to identify the heterogeneity of CAFs by Seurat pipelines.

Cell-cell communications based on ligand-receptor interactions were inferred by CellphoneDB (Efremova et al., 2020). To gain more critical cell-cell interactions in the colon cancer tumor microenvironment, we selected receptor-ligand pairs associated with hub genes for further analysis, aiming to explore potential interactions between immune cells. Significantly differential expressed ligand-receptor pairs (p < 0.05) were visualized.

Single-cell trajectories and determination of the continuous process of CAFs were analyzed by Monocle 2.0 package (v 2.10.0) (Jin et al., 2021). Pseudo-temporal analysis was applied to classify cells in pseudo-chronological order using the top 1000 differentially expressed genes in fibroblasts. Subsequently, a branch expression analysis model (BEAM analysis) was used to analyze branch fate-related genes.

CAFs-associated gene signatures were generated by identifying the marker genes of all CAFs clusters. The activities of these genes in each sample of all CRC datasets were calculated using GSVA. Log-rank test and Cox proportional hazards regression were performed to explore the relationship between CAFs characteristics and patient prognosis, including overall survival (OS) and recurrence-free survival (RFS) rates. Cutoffs for different cell characteristics in different public datasets were determined by the survminer package and used for plotting Kaplan-Meier curves.

Mutation data of CRC were obtained from the TCGA, and analyzed using the TCGAbiolinks package (Colaprico et al., 2016). Mutation landscape and lollipop plots were generated using maftools (Mayakonda et al., 2018).

We predicted the prognostic characteristics of CRC patients by identifying CAF marker genes from scRNA-seq clusters. Using GSE192667 as training dataset, all CAF marker gene were investigated by univariate Cox regression models for the prognostic evaluation of OS time. Genes with significant prognostic effect (p < 0.05) were determined as candidate prognostic genes. The LASSO regression analysis was then used to identify the feature genes and optimize the model to prevent overfitting. Based on the coefficients generated from the LASSO analysis, a risk score was assigned to each colon cancer patient. Finally, we divided all colon cancer patients into high- and low-risk groups based on their risk score by the median. The association between the risk score and OS was assessed using Kaplan-Meier analysis. Heatmaps were generated to visualize the associations between CAF risk scores and candidate genes. The time-dependent prediction accuracy of our model in the training, internal, and external test datasets was evaluated using AUC.

After obtaining differentially expressed genes between CAF subtypes, Metascape (https://Metascape.org) was used for gene set enrichment analysis. To estimate the infiltrating immune cells in the tumor microenvironment, CIBERSORT package was used to infer the relative abundance of immune cells in each sample (Yoshihara et al., 2013). Gene sets of tumor-associated canonical pathways were obtained from previous study (Sanchez-Vega et al., 2018). Activities of these gene sets were generated by single cell gene set enrichment analysis (ssGSEA) for cell state assessment of each tumor sample.

Tumor and adjacent tissue samples were fixed in formalin and embedded in paraffin. For IHC staining, sliced samples were deparaffinized and rehydrated. After that, antigen retrieval was performed, and normal goat serum was used for 10 min at room temperature to block non-specific binding site. Each slide was treated with mouse monoclonal anti-human NNMT antibody (1E7, diluted 1:1400) and incubated in 37°C for 40 min. Then, slides were incubated with biotinylated goat anti-mouse antibody for 30 min, and chromogenic reaction was carried out using a diaminobenzidine (DAB) Substrate Kit. Finally, Digital slide scanning system (KF-PRO-005) was used to capture images of IHC. The staining scores of NNMT protein expression were evaluated by two independent pathologists based on their clinical information. The protein expression levels were classified as 0 (no staining), 1+ (weak staining), 2+ (moderate staining), or 3+ (intense staining), and the staining score was calculated by integrating the percentage of positive cells and the respective intensity scores. The staining score ranged from a minimum value of 0 to a maximum value of 300.

Chemosensitivities of high and low CAF score groups were evaluated by oncoPredict (Maeser et al., 2021). Briefly, a ridge regression model with 10-fold cross-validation was built to infer half-maximal inhibitory concentrations (IC₅₀) Value. Pharmacogenomics database Genomics of Cancer Drug Sensitivity (GDSC; https://www.cancerrxgene.org) (Yang et al., 2013) was used to assess the response of CRC patients to chemotherapy. In addition, the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (Jiang et al., 2018) was implemented to predict immune checkpoint blockade treatment response between the two groups.

Statistical analyses were conducted using R software (v 4.1.0) and data visualizations were generated using appropriate R packages. Non-parametric tests such as the Wilcoxon rank-sum test were used for comparing two groups of continuous variables that did not follow normal distributions, while the Kruskal–Wallis test was used for testing three or more groups. Cox regression was employed to calculate hazard ratios (HR), and Kaplan-Meier analysis was used for prognostic evaluation. Statistical significance was set at a two-sided p-value of <0.05. Spearman’s correlations were determined for correlation analysis (p < 0.05, **p < 0.01).

Six samples obtained from the SMC cohort (GSE178318) were included in our study, which underwent quality control based on cell characteristics and mitochondrial and ribosomal gene expression. Subsequently, dimensionality reduction was performed to classify all cells. T-Distributed Random Neighborhood Embedding (t-SNE) was used to divide the cells into 9 major clusters and 32 more detailed minor clusters (Figure 1A). To further explore the metabolic signature of tumor fibroblasts, we calculated their metabolic signature pathways scores using SCmetabolism and clustered these CAFs based on these scores. Two fibroblast subclusters were identified (Figure 1C), and the important differentially expressed pathways for each cluster were visualized by heatmap for 11 metabolic pathways, including carbohydrate metabolism, energy metabolism, lipid metabolism and other nutrients metabolism (Figure 1D). According to the expression activities of metabolic pathways, two fibroblast subclusters were determined as hypermetabolism CAFs (hyperCAFs) and hypometabolism CAFs (hypoCAFs).

In primary CRC, total CAFs were significantly more abundant than in liver metastases (Figure 1B). Proportions of both CAFs subtypes also were greater in primary tumors than in metastases (Figure 1E). CAF in CRC patients showed significant heterogeneity in different sample, COL12 having a more even proportion in CAFs cells than others (Supplementary Figure S1). However, the proportion of hypoCAFs in other patients was significantly lower than hyperCAFs.

We identified differential expressed genes between clusters (Figure 1F). Hypometabolic CAFs shows higher expression of genes involved in myogenesis and pericyte-associated markers, such as MALAT1. While hypermetabolic CAFs highly expressed THY1, COL1A2, and some metabolism-related genes, such as PLA2G2A and NNMT. We focused on nicotinamide N-methyltransferase (NNMT), which is a cytosolic enzyme that has been identified as a significant metabolic regulator of cancer-associated fibroblasts (Eckert et al., 2019). Our findings were supported by immunohistochemical staining, which revealed that NNMT is predominantly expressed in fibroblasts (Figure 1G). The expression of other metabolism-related genes, including CRABP2, PLA2G2A, OGN, were all highly expressed in hyperCAFs subtypes (Supplementary Figure S2).

To describe and explain the functional heterogeneity of the two CAF subpopulations, several sets of genes characterizing the related functions of CAFs were used. Heatmap showed that different metabolic CAF subpopulations were characterized by significant differences in the expression of collagen genes, angiogenesis genes, smooth muscle-related contractile genes, and members of the RAS superfamily (Figure 2A).

Trajectory analysis of CAFs was performed based on the Monocle 2 algorithm to infer the maturation process of CAFs (Figure 2B, Supplementary Figure S3). In particular, we dissected gene patterns involved in CRC cell state transitions. Cell-to-cell communication analysis revealed large-scale interactions between the two CAF subpopulations and other cell types. Hypometabolism CAFs had cellular interactions among hypometabolism CAFs, tumor cells, and endothelial cells (Figure 2C), While hypermetabolism CAFs had the strongest interaction on tumor cells and mast cells (Figure 2C).

We compared different CAF subtypes between CRC and metastasis with genes related to tumor proliferation, metastasis, and progression pathways, to explore whether significant interactions were observed among different cell subsets. The results showed that CAFs, tumor cells, and B cells can participate in a series of functional interactions involving CXCL12 receptor-mediated APP, COPA, and MIF signaling (Figure 2D).

Activity of each transcription factor (TF) and its regulated genes were also inferred in both CAF subtypes. By comparing regulator specificity scores (RSS), we examined key regulators for each cell type and visualized the top 5 regulators (Figures 2E,F). E2F1, RUNX3, and ZNF224 were identified as the top regulators for hypometabolism CAFs, while CREB3L1, NFIB, ARNT, and GTF3A were identified as the top regulators for hypermetabolism CAFs.

Based on differentially expressed genes (DEG) (Supplementary Table S1) between two CAF subtypes, we performed the functional enrichment analysis using the online gene ontology (GO) enrichment analysis tool Matascape. Highly expressed genes of hypometabolism CAFs were enriched in the regulation of RNA splicing and muscle structure development GO terms, while genes of hypermetabolism CAFs were enriched in the transport of small molecules, cGMP-PKG signaling pathway, and some metabolism-related pathways such as fatty acid degradation (Figures 3A–D).

To investigate the association between CAF subtypes signatures and overall survival (OS) and recurrence-free survival (RFS) of CRC patients, we computed metabolic subtype scores for CAF subtypes by GSVA. Prognostic analysis was performed on all differentially expressed genes (DEGs) across nine publicly available cohorts that were classified into three types according to the metastasis status (i.e., primary tumor cohort, tumor cohort with metastases, and tumor cohort with liver metastases). We conducted a meta-analysis to obtain stable prognostic results for CAF subtypes, and compared our metabolic subtypes with Pan-CAF signatures derived from previous studies (Galbo et al., 2021). Our analysis revealed significant differences between subgroups in hypometabolism CAF scores compared to hypermetabolism CAF scores, in relation to RFS and OS. As a result, we defined hypometabolism CAF-type cells as a class of cells that are specific to colorectal cancer patients with LM.

To explore the prognostic genes associated with hypoCAF, we selected 20 genes based on the univariate Cox regression analysis in the GSE159216 dataset (Figure 4A). These 20 genes underwent Lasso-Cox regression analysis with 10-fold cross-validation to generate the optimal model, which highlighted 13 genes with the smallest partial likelihood deviation and optimal regression efficiency, including TINAGL1, ADIRF, ELP6, CSRP2, PPP1R15A, PABPN1, PHLDA1, ID3, KNOP1, DSTN, PPP1R10, CCDC107, and UBALD2. The risk score was calculated using the formula, and applied to GSE159216, TCGA, and GSE72970 datasets. Results indicated that colorectal cancer patients with high risk score had a higher mortality rate (Figures 4C–E). Heatmap results showed significant differences in the 13 genes expression between two risk score groups.

According to KM curves, patients in the high risk score group had lower survival rates than those in the low risk score group (Figures 5A–C). We also performed time-dependent ROC analysis, with the AUC values of our model in the training set for 1-year, 3-year, and 5-year overall survival being 0.85, 0.78, and 0.80, respectively (Figure Figure5A). In TCGA, the AUC for our 1-, 3-, and 5-year survival models were 0.66, 0.67, and 0.64, respectively, while in the GSE72970 dataset, the AUC values for 1-, 3-, and 5-year survival models were 0.74, 0.75, and 0.74, respectively. Cox regression analysis in three cohorts showed CAF score model would be an independent prognostic marker for CRC with LM (Figure 5G). We developed a nomogram based on CAF score to predict overall survival in CRC patients at 1, 3, and 5 years (Figure 5H). The calibration curves for each time point showed excellent predictive performance.

We analyzed the correlation between the hypoCAFs signature score and several pathways, the score was negatively correlated with most of the metabolic-related pathways. It was also negatively correlated with some cancer-related pathways such as cell cycle, gene duplication, homologous recombination, and p53 signaling (Figure 6A).

After deconvolution analysis, samples with low CAF scores exhibited high levels of activated memory CD4 T cells and memory B cells, while samples with high CAF scores showed high expression levels of regulatory T cells and resting memory CD4 T cells. To evaluate immune competence, we examined the expression of immune checkpoints (CD274, CTLA4, HAVCR2, IDO1, LAG3) and immune competence factors (CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF) (Figure 6B). The CAF score demonstrated a negative correlation with 14 out of 75 immunomodulators and 24 immune cells (Figure 6B). Finally, we compared the immune score (Figure 6D), ESTIMATE (Figure 6E), and stroma score (Figure 6F) between samples with high and low CAF scores, and observed that high CAF score samples exhibited elevated matrix and ESTIMATE scores.

To explore the role of CAF score in immunotherapy, we investigated the correlation between risk score and TMB. Our findings revealed that TMB expression was significantly higher in the low-risk subgroup than in the high-risk subgroup (Supplementary Figure S4). Moreover, to gain further insights into the nature of immunity in different risk subgroups, we analyzed genetic mutations. The top 20 genes with the highest mutation rates were identified in both the high-risk and low-risk subgroups (Supplementary Figure S4).

For ICB response prediction, we determined the correlations of CAF scores with TIDE, dysfunction, exclusion, and MSI Expr signature (Figures 7A–D). The results showed that CAF scores were positively associated with TIDE, dysfunction and exclusion, and negatively associated with MSI Expr sig. It was found that CRC patients with lower hypoCAFs scores had a higher possibility of responding to immunotherapy and may had better prognosis after immunotherapy, indicating that patients with lower CAFs scores were more likely to benefit from immune checkpoint therapy (Wilcoxon test, p = 0.0001, Figures 7E–H).

The expressions of five immune checkpoint molecules (PD1, PD-L1, CTLA4, LAG3 and HAVCR2) were compared between groups with high and low CAF score, the results showed that four immune checkpoint molecules (PD1, PD-L1, LAG3, and HAVCR2) were significantly upregulated in the high-risk group (Figure 7G). Furthermore, CAF score was strongly associated with the desert and inflamed immunophenotype (Kruskal–Wallis, p = 0.0029, Figure 7I). Our study also found that CAF score was positively correlated with PD-L1 expression in tumor cells and PD-1 expression in immune cells (Figure 7J).

Chemotherapy is the main treatment options for colorectal cancer liver metastases, therefor, whether CAF scores can accurately predict chemotherapy outcomes in colorectal cancer patients was investigated. GDSC is used to predict response to conventional chemotherapy in patients with colorectal cancer liver metastasis. A ridge regression model was built to predict IC₅₀ of different drugs. We found that the IC₅₀ of cisplatin, gemcitabine and other chemotherapeutics in the high CAF score group were significantly lower than those in the low CAF score group, suggesting that CAF score was positively correlated with chemotherapeutic drug sensitivity of colon cancer liver metastasis. In addition, we used the database to predict small molecule drugs (Figure 8B). The drugs vemurafenib, PLX-4720, dasatinib, and PI-103 were found to be negatively associated with the CAF score, with lower estimated AUC values in samples with high CAF score. These findings suggest that the predicted small-molecule drugs may be more effective in patients with high CAF score.