Six publicly available single-cell transcriptomics datasets of colorectal cancer were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), including GSE205506 (Li et al., 2023), GSE178318 (Che et al., 2021), GSE188711 (Guo et al., 2022), GSE225857 (Wang et al., 2023), GSE144735 (Lee et al., 2020), GSE158692 (Giguelay et al., 2022). To minimize batch effects introduced by different sequencing platforms, all six datasets were generated using the GPL24676 platform (Illumina NovaSeq 6,000, Homo sapiens). In total, the scRNA-seq data comprise 112 samples, including tumor tissue, normal tissue, and peripheral blood mononuclear cells (PBMCs) from 48 patients. The clinical information of all samples is summarized in Supplementary Table S1. In addition, one pancreatic ductal adenocarcinoma (PDAC) dataset (GSE197177) and one hepatocellular carcinoma (HCC) dataset (GSE156337) were downloaded from the GEO database.

A comprehensive analysis and visualization of single-cell transcriptomics data were conducted using the R package Seurat (version 4.4.0) (Hao et al., 2021), with all cells integrated according to sample ID by Harmony (version 1.2.0) (Korsunsky et al., 2019). To enhance data quality, a filtering criterion was applied, excluding cells with <400 or >50,000 unique molecular identifiers (UMIs) or >25% mitochondrial genes. To further eliminate doublet data, the scDblFinder R package (1.16.0) (Germain et al., 2022) was used. Log-normalization was then applied to standardize the single-cell data using the NormalizeData function, and the top 2,000 most variable genes were identified with the FindVariableFeatures function. After the gene expression matrices were transformed to the natural log scale by the ScaleData function, principal components analysis (PCA) was performed to reduce the dimensionality of the data. The RunHarmony function was subsequently employed to mitigate batch effects and integrate the datasets from the two samples. The number of included components (PCs) was determined based on the ElbowPlot function, with a total of thirty PCs retained. Cells were clustered using the FindNeighbors and FindClusters functions, and resolutions from 0.1 to 0.8 were explored for optimal cell clustering, with the final resolution set to 0.2. For data visualization, the Two-dimensional Uniform Manifold Approximation and Projection (UMAP) was applied using the RunUMAP function. The cells were then split into nine clusters based on commonly used cell markers: TNK cells expressing CD30 and NKG7 (n = 17,2941), B cells expressing CD79A and MS4A1 (n = 27,560), plasma cells expressing CD79A and JCHAIN (n = 47,326), myeloid cells expressing LYZ (n = 38,387), endothelial cells expressing PECAM1 (n = 14,386), fibroblasts expressing COL1A1 (n = 31,904), epithelial cells expressing EPCAM and KRT8 (n = 51,206), proliferating cells expressing MKI67 (n = 7,739), and mast cells expressing TPSB2 (n = 3,907).

To further resolve fibroblast-related subpopulations, a second round of UMAP clustering was conducted. The fibroblast cells were clustered in the 0.5 resolution. Subclusters were identified based on differentially expressed genes using the FindAllMarkers function in Seurat, with cell type annotations assigned according to the most highly expressed genes in each cluster. Among fibroblast-related subpopulations, some clusters encompassed other stromal cell types in addition to fibroblasts. Fibroblasts were distinguished using the markers DCN and LUM, while pericytes, smooth muscle cells, and glial cells were identified by RGS5, RERGL, MYH11 and GPM6B, respectively. Additionally, proliferating fibroblasts were characterized by the marker STMN1. For the remaining fibroblast subpopulations, cell type annotations were assigned based on differentially expressed genes and highly expressed markers (Du et al., 2024; Gao et al., 2024).

The FindMarkers function was employed to identify differentially expressed genes (DEGs), and the top 10 DEGs ranked by log₂FC were selected for heatmap visualization. The R package clusterProfiler (v4.10.1) (Yu et al., 2012) was used for functional annotation of DEGs through the enrichGO and enrichKEGG functions, with the top pathways or terms ranked by -log10(p value) visualized as bar plots or dot plots. The hypergeometric p-value was adjusted using the Benjamini–Hochberg correction. The GSVA R package (v1.40.1) (Hänzelmann et al., 2013) was also utilized, incorporating the hallmark gene sets obtained from the Molecular Signatures Database (MSigDB). Statistical significance was determined using the criteria p.val <0.05 and |avg. logFC| > 0.5.

The distribution of fibroblast subsets across various tissues was analyzed using the ratio of observed to expected cell numbers (Ro/e) (Zhang et al., 2018). This metric was calculated to determine the tissue-specific preference of each subset in three tissues groups. Expected cell numbers for each fibroblast subset-tissue combination were derived using the chi-square test. A fibroblast subset was considered to be enriched in a particular tissue if the Ro/e ratio exceeded 1.

Monocle2 (v2.28.0) (Qiu et al., 2017) was employed for pseudotime analysis to infer the developmental trajectories of cells within specific cell types based on their gene expression profiles. The UMI matrix was used as input, and variable genes that were detected by dispersionTable function were used for a building trace. A pivotal feature of monocle2, the plot_genes_in_pseudotime function, facilitated the visualization of gene expression changes along pseudotime, representing the inferred developmental trajectory of individual cells. For branch site differential genes analysis of fibroblasts, the BEAM function was used to detect genes that contributed most significantly when cells branched (q < 1 × 10⁻⁴) and the plot_genes_branched_heatmap function was employed for heatmap visualization.

Cell cycle phase assignment was performed using the CellCycleScoring function in Seurat (v4.3.0) with the updated human cell cycle gene sets (cc.genes.updated.2019). For each cell, S phase scores were calculated as the average normalized expression of the corresponding gene sets, subtracted by the average expression of the remaining genes.

Cell-cell communication analysis was conducted using the CellChat (v1.6.0) software to model the intercellular pathway network in colorectal cancer. Standardized single-cell data were used as input, incorporating prior knowledge of interacting ligands, receptors, and auxiliary factors for modeling. The interactions between fibroblast subgroups and other cell populations, such as immune cells and epithelial cells, were focused on, with specific emphasis on the expression patterns of ligand-receptor combinations across these cell types. The significance of interactions was calculated using permutation tests. The CellChatDB.human database was employed by loading the “Secreted Signaling Pathways” to infer the cell-cell communication network. The analysis followed the procedures outlined on the CellChat official website (Jin et al., 2021). Centrality measures for each cell type within the network were calculated using the netAnalysis_computeCentrality function to assess their relative importance in cellular communication, while the netAnalysis_signalingRole_network function was used to identify the roles of cells in signal transduction, such as signal senders, receivers, or mediators, within the network.

The public gene expression data and detailed clinical information were obtained from USUC Xena (https://xena.ucsc.edu/). Survival analysis was performed by R package survival (v3.5–7) and overall survival from the TCGA COADREAD datasets was used. Hazard ratio was calculated by the Cox proportional hazards model and 95% confidence interval was reported, and the Kaplan–Meier survival curve was modeled by survfit function. The Kaplan–Meier survival curves were plotted by using ggsurvplot function and compared using the two-sided long-rank test. Specifically, to test the correlation of each stromal subset and patients’ survival, we used top10 DEGs as gene set for each subset. GSVA enrichment scores for each gene set were calculated using the R package GSVA (v1.50.5) from Bioconductor. The mean of the enrichment scores was then used as a threshold to group samples into high and low groups for downstream survival analysis.

The single-cell regulatory network inference and clustering (SCENIC) package (v1.1.2.1) (Aibar et al., 2017) was used to analyze the TF activity of fibroblasts from mCC and mLC samples. The super cells were generated by combining the data from every five single cells in each cluster to reduce computing resource consumption. The mean values of normalized counts from five single cells were calculated as the raw input data for SCENIC. The transcription factor prediction website KnockTF2.0 (https://bio.liclab.net/KnockTFv2/index.php) was also used to predict key transcription factors (TFs) of PRELP⁺ CAF clusters in mLC samples.

The formalin-fixed paraffin-embedded (FFPE) tissue sections of colorectal cancer and its liver metastasis sites were collected from the Pathology Department of the First Affiliated Hospital of Anhui Medical University. The PRELP⁺ CAFs, Myo Fibroblasts, T cells, B cells and Plasma were characterized by a multi-marker panel in FFPE tissue sections using dual immunofluorescence. Briefly, the sections were deparaffinized in xylene, followed by rehydration in a graded ethanol series. Antigen retrieval was carried out by heating the slides in citrate buffer (pH 6.0) at 95 °C for 30 min. After cooling, the slides were incubated with blocking buffer (5% BSA in PBS) for 1 h at room temperature to reduce nonspecific binding. Rabbit anti-human PRELP (Abcam, ab229719, 1:100) was incubated overnight at 4 °C with mouse anti-human COL1A1 (CST, #6648, 1:200), mouse anti-human aSMA (Fluidgm, 3141017D, 1:200), mouse anti-human CD3 (Santa Cruz, sc-59010, 1:80), mouse anti-human CD138 (Santa Cruz, sc-12765, 1:100), and mouse anti-human CD20 (Servicebio, GB14030-50, 1:150) individually in separate reactions. The following day, secondary antibodies conjugated to different fluorophores (Abcam, ab150077, ab150115,1:100) were applied and incubated for 1 h at room temperature. After washing, the DAPl Staining Solution (biosharp, BL105A) was used to stain the nuclei in the sections. Images were observed with Olympus microscopy and were analyzed with OlyVIA Version 4.1.

To investigate fibroblast heterogeneity in primary CRC and liver metastases, we retrieved all relevant single-cell RNA sequencing (scRNA-seq) data from public databases (Lee et al., 2020; Che et al., 2021; Giguelay et al., 2022; Guo et al., 2022; Li et al., 2023; Wang et al., 2023) and conducted a comprehensive scRNA-seq analysis. The dataset comprised 75 tumor samples, 30 adjacent normal tissue samples, and seven peripheral blood mononuclear cell (PBMC) samples from 48 patients (Figure 1a; Supplementary Table S1). We collected clinical information from all samples as much as possible, including age, treatment information, MSI-H/MSS, and AJCC stage information (Supplementary Table S1). Based on histological classification at diagnosis, the samples were categorized as follows: primary CRC without liver metastasis (nCC) and its normal control (nCN), primary CRC with liver metastasis (mCC) and its normal control (mCN), and metastatic liver cancer (mLC) with its normal control (mLN) (Figure 1a; Supplementary Figure S1a).

After quality control and batch effect correction, a total of 440,592 cells were retained for clustering, identifying nine major clusters. The cell types included 75.1% immune cells (47.9% TNK cells, 6.7% B cells, 10.8% plasma cells, and 9.7% myeloid cells), 11.6% epithelial cells, and 7.2% fibroblasts (Figure 1b). These clusters were annotated based on the expression of canonical cell-type marker genes (Supplementary Figure S1c; Supplementary Table S2). Analysis of the cellular composition across samples revealed that PBMC samples were predominantly composed of TNK cells, B cells, and myeloid cells, confirming high-quality cluster classification. Notably, we observed significant heterogeneity in cellular composition between tumor and normal samples, reflecting the inherent variability contributed by individual patient characteristics (Supplementary Figure S1b).

To further investigate the heterogeneity of fibroblasts, we retrieved 31,904 cells from the fibroblast clusters expressing the COL1A1 marker gene and subjected them to additional clustering, resulting in 12 subgroups. Based on marker gene expression, three subgroups were identified as non-fibroblast cells, including glial cells, pericytes (PCs), and smooth muscle cells (SMCs), while the remaining nine subgroups were classified as distinct fibroblast populations (Figure 1c; Supplementary Figure S1d; Supplementary Table S3). We then analyzed the sources of these fibroblast subtypes. MMP1⁺, PRELP⁺, and COL11A1⁺ fibroblasts were derived from tumor tissues and were not found in normal tissues, demonstrating a strong association with tumors and thus being identified as cancer-associated fibroblasts (CAFs). In contrast, ADAMDEC1⁺ Fibro, C3⁺ Fibro, CFD⁺ Fibro, DPT⁺ Fibro, Proliferating Fibro, and Myo Fibro were present in both tumor and normal tissues and were therefore classified under the Fibro group (Figure 1c). Myofibroblasts are typically defined by high expression of MYH11 and low expression of RERGL. In multiple tumor models, fibroblasts exhibiting elevated MYH11 expression have been reported to play roles in immune regulation, promotion of tumor cell proliferation and migration, as well as inhibition of apoptosis(Grout et al., 2022; Lin et al., 2025). Consistent with previous reports, ADAMDEC1⁺ fibroblasts, characterized by high expression of ADAMDEC1 and ADAM28, are involved in tissue remodeling and healing induced by inflammation (Jasso et al., 2022) and were predominantly found in colorectal samples (Figure 1d). In contrast, a specific fibroblast subgroup, the PRELP⁺ CAFs, was significantly enriched in liver metastasis samples compared to both nCC and mCC samples (Figure 1d; Supplementary Figure S1e,f; Supplementary Table S4). A total of 87 samples in five sample groups contained fibroblasts, and most samples included more than eight cell groups (Supplementary Figure S1f; Supplementary Table S8).

The top 10 upregulated genes from each fibroblast subtype were considered marker genes, and these marker genes exhibited subtype-specific expression patterns (Figure 2a). We identified four representative marker genes for each fibroblast subtype, consistent with previously characterized fibroblast populations reported in the literature, including ADAMDEC1⁺ fibroblasts (Jasso et al., 2022), CFD⁺ fibroblasts (Du et al., 2024), DPT⁺ fibroblasts (Peng et al., 2022; Peng et al., 2023), COL11A1⁺ CAFs and proliferating fibroblasts (Du et al., 2024). Notably, the expression pattern of the top 10 marker genes in PRELP⁺ CAFs showed some similarity to those in COL11A1⁺ CAFs; however, their overall transcriptional profiles were distinct.

Prolargin (proline/arginine-rich terminal-rich leucine repeat protein), encoded by PRELP gene, is a protein that bind type I collagen to basement membranes and type II collagen to cartilage (Lewis, 2003). In most tumors, prolargin is considered a tumor suppressor, but studies in colorectal cancer have shown that it promotes epithelial-mesenchymal transition (EMT) in colorectal cancer cells, which in turn promotes the growth and metastasis of colorectal cancer (Chen et al., 2015; Hong et al., 2020; Dozen et al., 2022; Gui et al., 2024). Other genes highly expressed in PRELP⁺ CAFs also encoded extracellular matrix proteins, including lectins (COLEC11), latent transforming growth factor (TGF)-beta binding proteins (LTBP2), and integrins (ITGBL1) (Figure 2a). GSVA enrichment analysis revealed that many genes in pathways involved in proliferation and differentiation, such as TGF-β signaling and Wnt-β signaling were highly expressed (Figure 2b; Supplementary Figure S2a). Similar to PRELP⁺ CAFs, COL11A1⁺ CAFs were characterized by high expression of extracellular matrix proteins, including COL11A1, COMP, CTHRC1, and MMP11 (Figure 2a; Supplementary Figure S2b). Additionally, COL11A1⁺ CAFs showed increased activity in pathways related to EMT and angiogenesis (Figure 2b). This finding further supports the idea that both PRELP⁺ and COL11A1⁺ CAFs play important roles in the tumor microenvironment, particularly in regulating tumor metastasis and angiogenesis. The pathway enrichment results for PRELP + CAFs in mLC and mCC tissues show that this subpopulation is energetically active in mCC tissues, suggesting that the cells are producing and storing large amounts of energy in preparation for upcoming energy-consuming life activities. In mLC tissues, this subpopulation is enriched in pathways related to extracellular matrix remodeling and adhesion and attachment functions, indicating the functions it performs (Supplementary Figure S3).

Distinct fibroblast distributions were observed across tumor sites. ADAMDEC1⁺ fibroblasts and MMP1⁺ CAFs were more abundant in primary colorectal cancer (mCC and nCC) than in metastatic liver cancer (mLC) (Figure 2c; Supplementary Figure S2c). In contrast, PRELP⁺ CAFs were predominantly found in mLC, with lower proportions in mCC and minimal presence in nCC (Figure 2c). To validate the scRNA-seq analysis results, we performed immunofluorescence assays and examined histological slides of both mCC and mLC samples. We conducted dual immunofluorescence staining using antibodies against PRELP, a protein highly expressed in PRELP⁺ CAFs, and COL1A1, a common fibroblast marker protein. We observed significantly fewer PRELP⁺ CAFs in mCC samples, consistent with scRNA-seq analysis (Figure 2d). This suggests that PRELP⁺ CAFs are associated with advanced colorectal cancer and may contribute to metastasis. Using the top 10 upregulated genes as a signature for PRELP⁺ CAFs, we calculated GSVA score for patients in the TCGA colorectal cancer cohort. Stratification based on these scores revealed that high PRELP⁺ CAF levels were significantly associated with poorer overall survival (Figure 2e), suggesting that PRELP⁺ CAFs play a significant role in tumor progression and advancement.

To further investigate whether CAFs resembling PRELP⁺ CAFs in CRCLM are present in other liver lesions, we analyzed single-cell datasets from pancreatic ductal adenocarcinoma liver metastases (PDADLM) and hepatocellular carcinoma (HCC) (Sharma et al., 2020; Zhang et al., 2023). Since previous studies did not report clusters with features similar to PRELP⁺ CAFs, we re-analyzed fibroblast subpopulations at higher resolution in datasets GSE197177 (PDADLM) and GSE156337 (HCC). In the PDADLM dataset, fibroblasts were clustered into 10 subpopulations. Among these, the C5 cluster showed high expression of PRELP⁺ CAF marker genes—including COLEC11, PRELP, and COL10A1 (Supplementary Figure S4a,b,e). However, unlike PRELP⁺ CAFs in CRCLM, which are enriched in liver metastases, the C5 cluster was similarly distributed between primary and metastatic sites (Supplementary Figure S4c). KEGG enrichment analysis further revealed that C5 was primarily associated with immune-related and cellular lifespan pathways—such as apoptosis, TNF signaling, and FoxO signaling—which only partially overlapped with the functional profile of PRELP⁺ CAFs (Supplementary Figure S4g). Thus, we propose that C5 represents a fibroblast subset in PDADLM that shares certain phenotypic features with PRELP⁺ CAFs but differs in tissue distribution and functional emphasis. In contrast, among the five fibroblast subpopulations identified in HCC, none exhibited marker gene expression similar to that of PRELP⁺ CAFs (Supplementary Figure S4d,f). These findings suggest that PRELP⁺ CAFs are not a universal feature of liver-associated fibroblasts, but rather a distinct population specifically associated with colorectal cancer liver metastasis.

To investigate the phylogenetic relationships among fibroblast subtypes, we performed pseudotime trajectory analysis using Monocle2, revealing their developmental trajectories. The analysis results for all sample groups indicated a complex differentiation pathway, where fibroblasts transition from initial progenitor-like stages into three distinct terminal states: nCC state, mCC state, and mLC state, as determined by the proportional distribution of fibroblast subtypes across sample groups (Figure 3a; Supplementary Figure S5a-c). The distribution of fibroblast subtypes across these states exhibited distinct preferences. CFD⁺ fibroblasts and DPT⁺ fibroblasts were highly concentrated in the early developmental stages, with their density gradually declining along the trajectory. In contrast, PRELP⁺ CAFs showed a progressive increase in density, reaching their peak at the mLC terminal stage (Figure 3a). Additionally, MMP1⁺ CAFs predominantly occupied the terminal branch of the nCC state, while ADAMDEC1⁺ fibroblasts clustered at the terminal branches of both the nCC state and the mCC state. Meanwhile, COL11A1⁺ CAFs and C3⁺ Fibroblasts were positioned in intermediate transitional states (Figure 3a; Supplementary Figure S5a,b).

Additionally, we examined the expression of proliferative genes across all subpopulations and conducted cell cycle analysis, which further corroborated the pseudotime trajectory results. With the exception of Proliferating Fibro, none of the other subpopulations exhibited high expression of proliferative genes (Supplementary Figure S6a). S-phase scoring revealed that cell clusters at the initial differentiation stage—such as CFD⁺ Fibro and DPT⁺ Fibro—displayed overall higher S-phase scores, suggesting enhanced proliferative and differentiation potential. These populations showed significant differences compared to terminal-stage PRELP⁺ CAFs (Supplementary Figure S6b). The differentiation trends identified in separate analyses of the mCC and mLC sample groups were consistent with those from the integrated analysis. PRELP⁺ CAFs were primarily localized at the terminal end of the differentiation trajectory, although a subset was also observed at intermediate stages (Supplementary Figure S6c). Together, these results indicate a dynamic differentiation process among fibroblast subpopulations, highlighting their distinct functional contributions to colorectal cancer progression and metastasis.

To elucidate the molecular mechanisms underlying these transitions, we identified DEGs along the pseudotime trajectory. These DEGs were clustered into three distinct expression patterns corresponding to the initial developmental stage (pre-branch), mCC state, and mLC state (Figure 3b). The heatmap clearly highlights three distinct gene sets highly expressed in each state, including marker genes for fibroblast subtypes such as PRELP, C3, and CFD, along with numerous functional genes (Figure 3b). To further investigate the biological significance of these genes, KEGG pathway enrichment analysis was performed. The results revealed that genes in the mLC state cluster were primarily enriched in pathways related to ECM organization, ECM-receptor interactions, the PI3K-Akt signaling pathway and the p53 signaling pathway (Figure 3b). The PI3K-Akt signaling pathway is known to promote tumor cell proliferation, survival, and metastasis (Leiphrakpam and Are, 2024). And the p53 signaling pathway mediates the TGFβ–SMAD signaling pathway activation of fibroblasts, which can contribute to an immunosuppressive TME and EMT (Efe et al., 2024). These findings suggest that PRELP⁺ CAFs represent the terminal differentiation state of fibroblasts within liver metastatic tumors, highlighting their role in shaping and maintaining the TME in liver metastases.

CAFs contribute to tumor progression by interacting with immune cells through cytokines, chemokines, metabolites, and exosomes (Li et al., 2021). To determine whether PRELP⁺ CAFs are more closely adjacent to immune cells than other subtypes, we performed immunofluorescence staining and compared with Myofibroblasts. The results showed that PRELP⁺ CAFs co-localized spatially with multiple immune cell types (including T cells, B cells, and plasma cells) and were closer to immune cells than Myofibroblasts (Figure 4a). This finding suggests potential interactions between PRELP⁺ CAFs and immune cells within the TME. To further elucidate how PRELP⁺ CAFs modulate the TME, we conducted a cell-cell communication analysis using CellChat on all identified subsets in mLC. To enhance resolution, TNK cells were subdivided into 12 subtypes, while both myeloid and epithelial cells were further clustered into 14 subtypes each (Supplementary Figure S7a; Supplementary Tables S5–S7). This refined classification allowed for a comprehensive evaluation of intercellular signaling dynamics within the TME.

Cellular communication analysis showed that PRELP⁺ CAFs were a highly signaling-active cell type in mLCs, suggesting that they play an important role through strong interactions with other cell types in the tumor microenvironment (Figure 4b). The amyloid precursor protein (APP) is upregulated in colorectal cancer and is linked to increased tumor cell proliferation, migration, and invasion (Pandey et al., 2016). The immunosuppressive effects of the APP-CD74 axis primarily impact B cells and myeloid cells, as observed in testicular tumors (Chen et al., 2024). Additionally, the collagen-CD44 axis plays a key role in mediating interactions between fibroblasts and various immune cells, facilitating the directional migration of immune cells away from the tumor parenchyma and toward the tumor stroma (Yang et al., 2025). Our findings indicate that the APP and collagen signaling pathways mediated by PRELP⁺ CAFs may inhibit immune cell function, promote immune evasion by tumor cells, remodel the ECM, and influence immune cell infiltration and activity through the APP-CD74 and collagen-CD44 axes (Chen et al., 2024; Cheung et al., 2025) (Figure 4c).

Among the 12 TNK cell subsets, PRELP⁺ CAFs exhibited strong signaling, particularly toward CD8⁺ Tem, CD8⁺ Trm, and CD4⁺ Tregs (Supplementary Figure S7b). For example, PRELP⁺ CAFs acted as the primary signal senders in Galectin signaling pathway, with PRELP⁺ fibroblasts generating LGALS9 (Galectin-9). This ligand binds to the glycosylated structure of the CD45 receptor, expressed on all TNK cell subsets (Peng et al., 2024) (Figure 4d; Supplementary Figure S7c). This interaction impairs anti-tumor immune responses and suppresses CD8⁺ Tem proliferation. The TIGIT-NECTIN2 immune checkpoint axis is another key interaction in liver cancer, promoting the creation of an immunosuppressive microenvironment that supports cancer growth (Ho et al., 2021). TIGIT, an inhibitory receptor highly expressed by CD4⁺ Tregs, is activated by NECTIN-2, expressed on PRELP⁺ CAFs (Figure 4d; Supplementary Figure S7c). This activation enhances the immunosuppressive function of Tregs, promoting the secretion of IL-10 and TGF-β(Riquelme et al., 2018). In colorectal cancer research, it has also been reported that cancer-associated fibroblasts (CAFs) can directly suppress effector T cell function, impair proliferation, and promote exhaustion through the TIGIT–NECTIN2 axis, thereby driving immunosuppression (Agorku et al., 2024; Stary et al., 2024). In addition to acting as signal senders in interactions with TNK cells, PRELP⁺ CAFs also function as signal receivers, engaging in crosstalk with myeloid cells. For instance, LAMP3⁺ DCs express SEMA7A, which may regulate gene expression in PRELP⁺ CAFs through SEMA7A-ITGB1/A1 (integrin β1 and α1) interactions. Notably, ITGA1 has been identified as a pro-malignant biomarker in pancreatic cancer, promoting drug resistance and metastatic potential (Gharibi et al., 2017) (Figure 4e). Furthermore, PRELP⁺ CAFs interact specifically with BATF3⁺ DCs, characterized by high expression of CLEC9A, through multiple signaling pathways, including CD226, SEMA4, THY1 and NECTIN (Figure 4f; Supplementary Figure S7d). These findings suggest that the phenotype and function of PRELP⁺ CAFs are shaped by intercellular communication and that they actively contribute to the establishment of an immunosuppressive tumor microenvironment.

As PRELP⁺ CAFs represent the terminal differentiation state of fibroblasts in colorectal cancer liver metastases, we sought to identify the transcription factors (TFs) driving their distinct phenotype. SCENIC analysis revealed NR2F2, NFIC, JUND, JUN, and NFIA as the top 5 TFs with significantly elevated regulatory activity in PRELP⁺ CAFs compared to other fibroblast subtypes in metastatic liver tumors (Figure 5a; Supplementary Figure S7e). Among these, NR2F2 emerged as the most prominent regulator, with its target genes showing the highest level of activation in PRELP⁺ CAFs (Figure 5b). NR2F2, a ligand-inducible nuclear receptor of the steroid/thyroid hormone receptor superfamily, is mechanistically linked to colorectal cancer progression through its role in TGF-β-dependent EMT, a process critical for tumor cell invasion, metastasis, and poor clinical outcomes (Wang et al., 2015; Wang et al., 2017). To validate the SCENIC findings, we further analyzed the top 50 upregulated DEGs in PRELP⁺ fibroblasts using the KnockTF database to predict potential TF regulators. Remarkably, 90% of these DEGs in PRELP⁺ CAFs were predicted to be regulated by NR2F2, including key markers such as PRELP, COLEC11, ITGBL1, LTBP2, and collagen gene (Figure 5C). These ECM-related genes are essential for remodeling the tumor microenvironment and reinforcing the matrix-forming characteristics of PRELP⁺ CAFs. Furthermore, NR2F2 expression strongly correlated with PRELP levels, the defining marker of this CAF subset (Supplementary Figure S7f), solidifying NR2F2 as a key driver of both the molecular identity and functional specialization of PRELP⁺ CAFs.

In addition to NR2F2, the transcription factors JUN (encoding c-JUN) and JUND, members of the AP-1 transcription factor complex, exhibited significant transcriptional activity in PRELP⁺ CAFs (Figure 5b). The AP-1 complex, composed of JUN, FOS, and ATF family proteins, is a well-established mediator of pro-tumorigenic immune modulation (Eferl and Wagner, 2003). Mechanistically, c-JUN shapes the immunosuppressive phenotype of PRELP⁺ CAFs by modulating the expression of key immune-related proteins. For example, c-Jun not only directly binds to the AP-1 regulatory element in the APP promoter, activating APP expression (Davis, 2015), but has also been reported to enhance NECTIN2 transcription, a key molecule involved in immune evasion (Lui et al., 2006). Moreover, AP-1 signaling also amplifies TGF-β1 production, which induces galectin-9 expression in colorectal cancer cells (Selnø et al., 2020). Galectin-9 interacts with TIM-3 on immune cells, exacerbating T-cell exhaustion and myeloid cell dysfunction—key mechanisms of immune evasion in metastatic niches.

Together, these findings highlight the intrinsic regulatory mechanisms underlying the diverse functions of PRELP⁺ CAFs, emphasizing their dual roles in stromal remodeling and immunosuppressive regulation within the metastatic tumor microenvironment.