Cancer-associated fibroblasts shape the formation of budding cancer cells at the invasive front of human colorectal cancer

Abstract

The formation of budding cancer cells at the invasive front of solid tumors is one of the first steps of metastasis. However, this process is still incompletely elucidated. Here, we used spatial molecular imaging to disentangle the complex interactions between cancer cells and the tumor microenvironment at the invasive front of colorectal tumors. Employing a 1000-plex gene panel, we defined all major cell types in tumors and adjacent normal tissue with accurate spatial information. Individual cancer cell clusters were located together, consistent with an expected mutation- and epigenetic-driven clonal evolution. However, cancer cell clusters encompassing budding cells exhibited a markedly different spatial distribution as they also contained cells that were scattered around the periphery of the main cancer cell masses. Moreover, these cells were frequently in contact with cancer-associated fibroblasts (CAFs) and underwent broad gene expression changes, mainly related to epithelial-mesenchymal transition (EMT), remodeling of the extracellular matrix (ECM), and migration. In addition, we defined an 11-gene signature (TYK2, IL2RG, KRT17, HLA-B, NPPC, WIF1, IL32, B2M, CCND1, CRIP1, ITGB1), which characterizes cancer cells en route to metastasis and is associated with inferior outcomes. Collectively, our findings suggest that CAFs induce pro-invasive gene expression changes involved in EMT, ECM remodeling, and migration.

Fig. 1. Spatial gene expression analyses of single cells in colorectal tumors reveal a unique spatial distribution of cell clusters that encompass budding cells.

A Multichannel immunofluorescent (IF) overlay of epithelial (PancK), immune (CD3) and nuclei (DAPI) immunofluorescent markers in green, yellow, and blue, respectively, from a field of view (FOV) containing colorectal adenocarcinoma tissue. Red colorings represent the fiducial beads. B UMAP projection showing the cell type composition of the tissue. C Spatial plot showing the polygonal shape, type, and location of each cell from the colorectal adenocarcinoma tissue. D Spatial plot showing the cellular overlay between the multichannel IF image and each cell detected, represented by a point at each cell centroid and color-labelled by cell type. E Spatial plot showing the overlay of the IF image with the cluster “Cancer cells 4” (light blue), showing a dispersed distribution across the tissue with cells located mainly in the outskirts of the tumor masses and in tumor buds. F Heatmap of the top 10 gene markers specific for each cluster. G Violin plots showing the expression levels of the epithelial markers EPCAM, CDH1, and PIGR for each cluster. Each dot corresponds to a single cell. Significance was determined by unpaired Student’s T tests with “Cancer cells 4” as the reference cluster. H Violin plots showing the expression levels of mesenchymal markers VIM, LUM, and FN1 for each cluster. Each dot corresponds to a single cell. Significance was determined by unpaired Student’s T tests with “Cancer cells 4” as the reference cluster. I Bubble plot showing the scored average expression (blue to red) and the percentage of expression (dot size) per cluster for each EMT state. The red scattered box highlights the “Cancer cells 4” cluster. ns (non-significant), * (P-value < 0.05), ** (P-value < 0.01), *** (P-value < 0.001), **** (P-value < 0.0001).

Fig. 2. A gene signature of early budding is associated with disease free survival.

A Venn diagram showing the unique and common genes being upregulated when comparing the “Cancer cells 4” cluster (CC4) and all the other cancer cell clusters (CCs), and when comparing CC4 and the CAFs cluster. B Dot plot showing the average expression and percentage of expression of each of the 11 genes that were upregulated in both comparisons in A. The red scattered box highlights the “Cancer cells 4” cluster. C Box plots showing the interquartile range of the average expression data for the 11 genes identified in A, extracted from TCGA (COAD). The horizontal line in the box indicates median value and the whiskers indicate 10^th and 90^th percentiles. Colorectal adenocarcinoma (T), separated by Microsatellite instability status: MSI high (MSI-H), MSI low (MSI-L), and Microsatellite stable (MSS) were compared to adjacent normal tissue (N). Significance was determined by one-way ANOVA. D Kaplan-Meier plot showing the disease-free survival based on the average expression levels of the 11 genes identified in A. The two groups were defined based on the median. Significance was determined using a Mantel-Cox test. ns (non-significant), * (P-value < 0.05), ** (P-value < 0.01), *** (P-value < 0.001), **** (P-value < 0.0001). E Schematic of the patient cohort derived from Pelka et al. single-cell RNA-seq dataset (GSE178341) (Upper panel). UMAP plots showing the integrated single-cell data, stratified by tissue type: adjacent normal colon (N) and tumor (T) tissues. Cells are colored by tumor stage (pT1-4) and annotated with major cell type labels (Lower panel). F Bubble plot showing the scored average expression (blue to red) and the percentage of expression (dot size) per cluster for the Budding Signature (11 genes found in A) and each EMT state for the epithelial subcluster Stem/transition amplifying (TA)-like cells. The red scattered box highlights pT3 and pT4 stages. G Bubble plot showing the scored average expression (blue to red) and the percentage of expression (dot size) per cluster for the each of the 11 genes found in A from the epithelial subcluster Stem/transition amplifying (TA)-like cells. The red scattered box highlights pT3 and pT4 stages. H IHC on human colorectal cancer tissue sections extracted from the Protein Atlas database showing the protein staining corresponding to 10 of the 11 genes identified in A. Green arrows highlight invasive front areas and budding cells.

Fig. 3. Individual FOVs from high, mid, and low-density budding areas share upregulated and downregulated genes.

Volcano plots depicting P-adjusted values relative to log2 fold changes (FCs) for comparisons between the cancer clusters that encompass most budding cells and all other cancer clusters for each field of view (FOV). The top panel in red shows the two FOVs annotated as high-density budding areas. The mid panel in light yellow shows the six FOVs annotated as mid-density budding areas. The bottom panel in light blue shows the six FOVs annotated as low-density budding areas. Red dots depict genes with a P-Adjusted (Benjamini-Hochberg correction) value < 0.05 and a log2 FC > 1. Blue dots depict genes with a P-Adjusted (Benjamini-Hochberg correction) value < 0.05 and a log2 FC < 1. ns: non-significant.

Fig. 4. Enrichment analysis on high, mid, and low-density budding areas revealed a progressive adaptation of the budding cells to the surrounding stroma.

A Venn diagrams showing the overlap of upregulated genes in cancer cell clusters containing most budding cells relative to all other cancer cell clusters for the two high-density budding (HDB) field of views (FOVs) (left panel), the six mid-density budding (MDB) FOVs (middle panel), and the six low-density budding (LDB) FOVs (right panel). For the two HDB FOVs, 129 genes were commonly upregulated. For the six MDB FOVs, 13 genes were commonly upregulated, 16 genes were upregulated in five FOVs, and seven genes were upregulated in four FOVs. For the six LDB FOVs, five genes were commonly upregulated, 17 genes were upregulated in five FOVs, and 57 genes were upregulated in four FOVs. B Bar plots showing pathway enrichment analyses performed on the 129, 36, and 79 genes identified in the HDB, MDB, and LDB areas, respectively. The Top 10 terms ranked by P-Adjusted values for “Molecular Function” from Gene Ontology, “WikiPathways”, and “REACTOME” databases are shown. C Venn diagram highlighting the shared and unique pathways from the enrichment analyses shown in B.

Fig. 5. Cancer cells from the budding clusters are frequently in contact with CAFs.

Heatmaps showing z-scores of the cell-cell interaction values for each of the cell types per field of view (FOV). The FOVs are divided into three sub-groups based on their budding density, high, mid, and low-density budding areas. Green boxes highlight the budding clusters.

Fig. 6. Cancer cells in contact with CAFs change the ECM composition to facilitate migration and avoid apoptosis.

Dot plots showing the “ECM Organization Score” (A), “Collagen ECM Score” (B), “Focal Adhesion Score” (C), “Apoptotic process Score” (D), and “miRNA Targets Score” (E). The dots represent mean values, and the error bars indicate standard deviation. FOVs are separated into high-density budding (HDB), mid-density budding (MDB), and low-density budding (LDB) areas. Blue dots correspond to the mean score of the budding cells that are in contact with CAFs, whereas light green dots correspond to the mean score of budding cells not in contact with CAFs. The box plots to the right show the mean values for the HDB, MDB, and LDB areas. Error bars indicate standard deviations.