Single-cell atlas of epithelial tissue, intestinal adenoma, and CRC

To characterize colorectal precancerous and malignant lesions, 9 biopsies were performed, comprising 3 pairs of normal intestinal epithelial biopsies (3 cases), adenoma biopsies (3), and CRC biopsies, respectively, from the same patient (Figure S3A, Table S11). All CRC samples were pathologically diagnosed as moderately differentiated (G2), and no other primary tumors were found. For each sampling, individual cells were selected without prior cell filtering, and sequencing data were generated using a 10X chromium platform. After removing low-quality cells, 20,262 high-quality cells were retained for subsequent analysis. To assess different cell types based on gene transcription profiles, downscaling and unsupervised cell clustering were carried out using the Seurat package after removing the batch effect among multiple samples. As shown in t-distributed stochastic neighbor embedding (t-SNE) and UMAP, 26 major cell clusters were finally identified in all samples (Figure S3B), which were then defined as single-cell transcriptome profiles of normal-precancerous lesions and CRC. These clusters were categorized by marker genes into eight known cell lines: T cell (marked by CD3D), plasma B cell (marked by MS4A1), monocyte (marked by S100A9), epithelia cell (EC, marked by EpCAM), follicular B cell (marked by MZB1), Macrophage (marked by CD14, CD163, and CD68), mast cell (marked by KIT), and fibroblast (marked by DCN) (Figure S3C-D). In addition to the typical cell type indices, other genes that marked each cell type were identified (Fig. S3E). The distribution of different samples and different tissue types on UMAP, respectively, is depicted in Fig. S3F and G. Each cell lineage proportion considerably differed among different samples (Fig. S3H), suggesting strong heterogeneity.

Genetic characterization and transcriptional variation of CC

The stability of the tissue microenvironment depends on the mutual communication of different cells, and various mechanisms have evolved for this purpose, of which the most direct and efficient is through channels that directly connect the cytoplasm of neighboring cells [40, 41]. TJs and GJs are the two most classical modes of communication, and their dysfunction is associated with the development of various diseases, especially tumors [42, 43]. Our previous study showed that gap junction protein alpha 4 (Cx43, also named GJA4) is lowly expressed in normal intestinal mesenchymal tissues but highly expressed in CRC mesenchymal tissues and may impair the survival of CRC patients through CAF-related pathways [44]. Building on this finding, we further incorporated TJs into our analysis to explore the significance of CC in CRC from a more comprehensive perspective.

Information on the 47 CC molecules was acquired from the HUGO gene Nomenclature Committee portal, including 24 GJ proteins (connexins and pannexins) and 23 TJ proteins (CLDNs). The localization and regulatory mechanisms of CC in the TME are presented in Fig. 1A. To detect genetic variation in cancer cell communication molecules, 1,482 samples were selected with at least one mutation per 47 CC molecules in the TCGA pan-cancer dataset. The oncoplot revealed the most frequent somatic mutations linked with 10 CC molecules in pan-cancer tissues. Of the 1482 samples, mutations were found in 753 cases, with a mutation frequency of 50.84%. It was revealed that GJA8 exhibited the highest mutation frequency (11%), followed by GJA10 (10%) and GJA1 (7%). Moreover, missense mutation was the most frequent nutation, and ovarian cancer (TCGA-OV) was the most frequent mutated cancer (Fig. 1B). Considering that only 14 types of cancer indicated > 10 paired tumors and normal samples, transcriptome differential expression was compared between these cancers. Further analysis of messenger RNA (mRNA) levels of the 47 CC molecules revealed that the expression of various genes was reduced in several tumor tissues, including CLDN5, CLDN11, GJD3, CLDN19, and CLDN20 (Fig. 1C). These findings indicated that genetic variation is one of the essential factors that influence CC expression. Most cancers revealed a positive relationship between copy number variant (CNV) and mRNA expression levels, especially PANX1 (Fig. 1D). Intriguingly, mRNA expression levels were negatively associated with DNA methylation in a subset of cancer types, with an opposite trend in another subset (Fig. 1E). In addition, the prognostic significance of different CC molecule expressions varied across different types of tumors (Fig. 1F). Meanwhile, transcriptional patterns of CC molecules were markedly heterogeneous in normal and various cancer samples, suggesting a correlation between aberrant expression and transcriptional variants.

Fig. 1

Comparison of expression levels of cell communication (CC) molecules. A Patterns mapped on the BioRender website to reveal the regulatory mechanisms of CC in tumors and their functions in the immune microenvironment of tumors. B The waterfall diagram illustrates the most frequent 47 CC molecule’s somatic mutations in The Cancer Genome Atlas (TCGA) pan-cancer data. 50.81% represents the proportion of 753 samples with at least 1 mutation of the top 10 genes among 1482 samples with at least one mutation of 47 CC genes. The percentage value on the right side of each line in the image indicates the number of samples with the specific gene mutation divided by 1,482 samples which had at least one mutation among the 47 CC genes. We label different types of CC molecules in red [gap junctions (GJs)] and blue [tight junctions (TJs)], respectively. C The dot’s color = degree of fold change. Red = high and blue = low expression in cancer tissue. Fold change = mean (tumor)/mean (normal), p-values were used. Field realistic doses (FDR) was utilized for adjusting the t-test and p-value. The size of the bubble indicates FDR; the larger the bubble, the lower the FDR. Genes with > twofold change and significance (FDR > 0.05) were used to plot graphs. If no significant genes are present in a cancer type, that cancer type was not included in the final figure. D Bubble plots display the correlation between the Copy number variant (CNV) (D) as well as DNA methylation (E) and the expression of the mRNA levels. A positive correlation is reflected in red, while a negative correlation is indicated by blue. Darker colors indicate a higher correlation index. The FDR is indicated by the bubble size. F Bubble plots showing the results of a log-rank test of the survival of 47 CC molecules in the TCGA-CRC cohort. Red represents detrimental to survival and blue denotes favorable to survival. The FDR is represented by the bubble size. G The mutation profiles of 44 CC molecules in 544 CRC patients in the TCGA-CRC cohort; co-mutations are shown by the green, mutex-mutations are indicated by the red, and asterisks indicate P values (*P < 0.05,.P < 0.01). H Mutation frequency of 47 CC molecules in 544 CRC patients in the TCGA-CRC cohort. The small graph above is the Tumor Mutational Burden (TMB), and the numbers on the right indicate the mutation frequency of each gene and provide the proportion of each variant. I, J Uniform Manifold Approximation and Projection (UMAP) (I) and violin (J) plot indicates the CC feature level (generated by the “AddModuleScore” function) across different cell types in our single cell RNA data. K Violin plot showing the CC feature level across different tissue types

Due to the previous research base, the present study focused on CRC [44]. The TCGA database revealed substantial co-mutations between CLDN23 and GJB6, CLN23 and GJB7, CLDN14 and CLDN11, CLDN7 and CLDN17, etc. (P < 0.01, Fig. 1G). Of the 544 samples, 144 displayed CC molecular mutations, with a mutation frequency (26.47%). Among the 47 genes, GJA8 had the highest mutation frequency (4%) (Fig. 1H), which was the missense mutation. The correlation among CC molecules in TCGA-CRC is shown in Figure S4. To characterize the localization of CC molecules at the single-cell level, the expression of CC molecules was scored using the “AddModuleScore” function of the “Seurat” package in our scRNA data. The results showed that the CC molecules were up-regulated in ECs (Fig. 1I, J) and significantly down-regulated in cancerous tissues (Fig. 1K), compared with normal intestinal epithelium and intestinal adenomas. These results suggested that an imbalance of CC molecules may lead to CRC.

TJs are involved in the malignant transformation of intestinal epithelial cells (IECs)

The dynamic changes in molecular signaling during EC malignant transformation (a protagonist in TME) have attracted considerable attention [45, 46]. Our scRNA data showed that CC molecules had the highest expression levels in ECs compared with other cells (Fig. 1I, J). After subclustering all ECs, 11 distinct groups were identified (Fig. 2A). These cell fractions were assigned to four major cell types, including normal, adenoma, cancer, and normal/adenoma cells, based on the type of tissue origin of each cell (Fig. 2B–D). Specific genetic markers for each cell type are displayed in Fig. 2E, further reflecting a high degree of heterogeneity of ECs. Then, cell trajectory analysis was carried out using the Monocle toolkit to further elucidate possible evolutionary routes between cell types. The pseudotime trajectory axis based on transcriptional profiles revealed two distinct trajectories of transdifferentiation (Fig. 2F). Using CytoTRACE, Cluster 2 (high CytoTRACE score) was identified as the normal-lesion transdifferentiation initiation point (Clusters 2, 7, and 8 were predominantly derived from normal tissues, Fig. S5A and Fig. 2G–I). Moreover, normal IECs were further classified into different functional subpopulations based on specific gene markers (Figure S5B). It was found that the two routes had distinct characteristics; specifically, a fraction of normal IECs (Cluster 2, Cluster 7, and Cluster 8) transdifferentiated into adenoma cells, whereas another fraction (mixed with adenoma cells) eventually transdifferentiated into cancer cells (Fig. 2J, K). To characterize the function of CC during transdifferentiation, CLDN3, CLDN4, and CLDN7 were selected to represent TJs, and GJA4, GJB1, and GJC1 were selected to represent GJs, respectively. Pseudotemporal expression dynamics of the representative genes showed that TJ levels increased and then decreased during normal-adenoma transdifferentiation and decreased directly during normal-cancer transdifferentiation. However, no change in both transdifferentiation routes was observed in GJs, which were always maintained at low levels (Figs. 2L, M and S6, S7). Further, multiplex immunofluorescence (mIF) staining of independent normal and cancer resection specimens confirmed that epithelial TJ levels were down-regulated in CRC tissues (as evidenced by an attenuated pattern of co-expression of EC-specific marker proteins epithelial cellular adhesion molecule [EpCAM] and CLDN4), whereas GJ levels did not differ significantly between normal intestinal epithelium and CRC tissues (Fig. 2N–Q). Loss of CLDN expression usually accompanies TJ destruction during tumor progression and leads to the acquisition of a malignant phenotype in cancer cells [47, 48]. Thus, the ultrastructure of CRC cells was observed under TEM. Compared with normal controls, CRC cells exhibited reduced TJ-mediated barrier formation and loose junction structures (Fig. 2R–S). These findings suggest that TJs are involved in the transdifferentiation of IECs to adenoma and carcinoma cells.

Fig. 2

Heterogeneous landscape of CC molecules across different lesions in epithelial cells. A–C UMAP plot of all epithelial cells, color-coded for eleven seurat clusters (A), three tissue types (B), and four cell types (C). D The fraction of three tissue types in four cell types. E Heatmap showing differentially expressed genes among the four cell types (fold change > 1.5, FDR < 0.01). F The trajectories of all epithelial cells constructed by Monocle 3. Each point corresponds to a single cell and is colour coded by pseudotime. G–I Box (G) and t-distributed stochastic neighbor embedding (t-SNE) (H, I) plot demonstrate the degree of differentiation of cluster 2, 7, and 8 (three normal epithelial cell clusters) assessed by CytoTRACE. J–K Monocle 3 demonstrates two trajectories of cellular differentiation present in epithelial cells, including from normal cells to adenoma cells (J) and from normal/adenoma cells to cancer cells (K). L, M Two-dimensional plots showing the dynamic expression of representative CC molecules during the epithelial cell transitions during the pseudotime. N Multiplex immunofluorescence (mIF) staining images of CLDN4 (green), EPCAM (pink), and GJA4 (red) in a resected normal colon specimen (blue, DAPI), and CLDN4 (for TJs) is upregulated in epithelial cells (marked by EPCAM). Scale bars are labeled on the graph. O mIF staining images of CLDN4 (green), EPCAM (pink), and GJA4 (red) in a resected colon cancer specimen (blue, DAPI), And CLDN4 (for TJs) is downregulated in epithelial cells (marked by EPCAM). Scale bars are labeled on the graph. P, Q Co-localization was determined using the Pearson correlation coefficient in normal colon specimens (P, R = 0.9459, P < 0.0001) and colon cancer specimens (Q, R = 0.5616, P < 0.0001), respectively. The co-localization relationship between CLDN4 and EPCAM was weaker in tumor tissue compared to that of normal tissue. The X-axis represents each pixel point on the image, and the Y-axis represents the gray value corresponding to each pixel point. R, S The ultrastructure of junctions was examined using transmission electron microscope (EM) in in normal colon specimen and colon cancer specimen. Orange arrows indicate TJs and green arrows denote GJs

GJs are involved in the malignant transformation of fibroblasts

CAFs are important tumor stroma components [49, 50]. The signaling crosstalk that develops between CAFs and other cells is necessary to maintain the TME [51, 52] (Fig. 3A, estimated using CellChat in our scRNA data). Here we further established the expression atlas of CC features in fibroblasts. Cluster analysis of all fibroblasts was performed, and two significant cell subpopulations (Cluster 0 and Cluster 1) were identified (Fig. 3B). Based on the proportion of single-cell tissue sources in the two clusters (Fig. 3C, D), CAFs were mainly represented by C0, indicating high expression of BGN, POSTN, and ACTA2, and NFs were mainly represented by C1, indicating high expression of OGN, POSTN, and MFAP5 (Fig. 3E, F). Next, gene set enrichment analysis (GSEA) was performed to decipher the differences in molecular characteristics between CAFs and NFs. Compared with NFs, CAFs were associated with signaling pathways such as cell-substrate junction, focal adhesion, actin filament binding, and regulation of the actin cytoskeleton (Fig. 3G, H). Notably, CAFs were mostly derived from NFs, which were recruited to the tumor region and reprogrammed into the former by cancer cells secreting cytokines. To further investigate this ongoing process, a trajectory analysis of CAFs and NFs was conducted (Fig. 3I). Contrary to the phenomenon observed in ECs, GJs but not TJs were involved in the transdifferentiation of NFs to CAFs. Specifically, GJ genes were all up-regulated during this differentiation, while the levels of TJ genes remained unchanged (Figs. 3J and S8). mIF staining of independent samples also confirmed that GJ protein expression was markedly increased in tumor mesenchymal tissues (high expression of actin alpha 2 [ACTA2] and low expression of cytokeratin) than in normal mesenchymal tissues, whereas no difference was found in TJs (Fig. 3K–N). Meanwhile, a higher density of GJ channels was observed in the ultrastructures of tumor stroma (Fig. 3O, P, Green arrow). To spatially confirm these co-localization and co-expression relationships, normal and cancerous were further studied. Spatial transcriptomics sections from GEO database were distinguished between epithelium and stroma by pathological sections and specific gene markers (EpCAM for epithelium and ACTA2 for stroma). Similarly, the co-localization of TJs and parenchyma was weakened and the co-localization of GJs and stroma was strengthened in CRC compared with normal intestinal tissue (Fig. 3Q, R). Collectively, these data imply that GJs may be important in the activation of CAFs.

Fig. 3

The heterogeneous landscape of CC molecules across different lesions in fibroblasts. A Circle plot showing the possible ligand-receptor pairs between fibroblasts and other type cells (predicted by CellChat). B The UMAP plot of all fibroblasts, color-coded for two seurat clusters. C, D The fraction of three tissue types in two cell types showed by histogram (C) and UMAP plot (D). E All fibroblasts were defined as Normal fibroblasts (NFs) and Cancer-associated fibroblasts (CAFs), respectively, according to tissue origin. F Heatmap showing differentially expressed genes between the two cell types (fold change > 1.5, FDR < 0.01). G, H The bubble plots indicate the up-regulated gene set in NFs (G) and CAFs (H), differently. I Trajectory of fibroblasts constructed by Monocle 3. Each point corresponds to a single cell and is colour-coded by pseudotime. J Two-dimensional plots showing the dynamic expression of representative CC molecules during the fibroblast transitions along the pseudotime. K mIF staining images of ACTA2 (pink), CLDN4 (green), and GJA4 (red) in a resected normal colon specimen (blue, DAPI), and GJA4 (for GJs) is lowly expressed in normal mesenchymal tissues (marked by ACTA2). Scale bars are labeled on the graph. L mIF staining images of CK (orange), ACTA2 (pink), CLDN4 (green), and GJA4 (red) in a resected colon cancer specimen (blue, DAPI), and GJA4 is highly expressed in cancer mesenchymal tissues (marked by ACTA2 and CK). Scale bars are labeled on the graph. M, N Co-localization was evaluates based on the Pearson correlation coefficient in normal colon specimens (M, R = 0.2692, P < 0.0001) and colon cancer specimen (N, R = 0.8806, P < 0.0001), respectively. The co-localization relationship between GJA4 and ACTA2 was stronger in tumor tissue compared to normal tissues. The X-axis represents each pixel point on the image, and the Y-axis represents the gray value corresponding to each pixel point. O, P The ultrastructure of junctions in normal colon specimens and colon cancer specimen was examined using a transmission electron microscope (EM). Orange arrows represent TJs, and green arrows represent GJs. Q, R Spatial transcription sections indicate the spatial expression of EPCAM, ACTA2, TJs markers (CLDN3/4/7), and GJs markers (GJA4, GJB1, and GJC1) in normal colonic tissue (Q) and colon cancer tissue (R). The dot color indicates the expression level of the markers. Green boxes for the parenchyma and pink boxes for the mesenchyme

Unsupervised learning identified two different CC-level patterns

The above results suggested that a spatiotemporal heterogeneity of the CC molecules existed during malignant transformation. To comprehensively understand the integrated mechanism of the CC feature in CRC, unsupervised clustering was performed on 620 samples from TCGA-CRC. Two unique modification patterns were identified, named Cluster 1 (C1, 525 cases) and Cluster 2 (C2, 95 cases) (Fig. 4A, B). PCA confirmed that the two clusters were distinguishable by the 47 CC molecule expression levels (Fig. 4C). The thermogram showed that TJs had lower levels but GJs had higher levels in C2 than in C1 (Fig. 4D). Considering previous findings (Figs. 2 and 3), we speculated that the stroma of C2 might be in a more active state, while its epithelial component might have a weaker malignant feature. We analyzed survival prognosis differences between the two CC subtypes, and results showed that C1 had a distinct and significant survival advantage, whereas C2 had a poorer prognosis (OS, P = 0.047; DFS, P = 0.000407; DSS, P = 0.00242; PFS, P = 0.535; Fig. 4E). Next, the CC phenotype was compared with several commonly used clinical indicators. It was found that most C1 patients were in the advanced stages, whereas most C2 patients were in the early stages (Fig. 4F). Notably, the abundance of stromal components was lower in the tumor of C1 patients and higher in the tumor of C2 patients (TCGA pathology slides, Fig. 4G). The “ESTIMATE” method showed that C2 had high stromal and immune signals (Fig. 4H), both of which have been shown to correlate with poor outcomes [53, 54]. In addition, immunosuppressive stromal cells or abundant stroma hinder the success of immunotherapy. A computational method was applied to model the two main mechanisms of tumor immune evasion and tumor immune dysfunction and exclusion (TIDE) to predict the ICB response based on transcriptional profiles. The results showed that C2 had a significantly higher TIDE score than C1, suggesting that C2 patients may not benefit from immunotherapy (Fig. 4I, P < 0.0001). Immunological estimates indicated that M2 macrophages were highly enriched in C2 and that M2 macrophages are a major modulator of immune tolerance in cancer cells andconfer resistance to immunotherapy (Fig. 4J).

Fig. 4

Unsupervised Machine Learning algorithms used to identify 2 molecular subtypes in TCGA-CRC. A Heat map showing the sample clustering at K = 2 (the optimal cluster number) in TCGA-CRC. B Left: The cumulative distribution function (CDF) curve in consensus cluster analysis. The consensus score’s CDF curves with various subtype numbers (k = 2, 3, 4, 5, and 6) are shown. Right: Relative change in area under the CDF curve for k = 2–6. C The TCGA-CRC samples were classified via Principal Component Analysis (PCA) based on the CC molecules expression profile. Different colors = C1 and C2 subtypes, respectively. Each point is a single sample. D The distribution of 47 CC molecules between two subtypes in TCGA-CRC. GJs molecules are upregulated in C2 and TJs molecules are upregulated in C1. E Survival analysis in terms of Overall Survival (OS), Disease-Specific Survival (DSS), Progression-Free Survival (PFS), and Disease-Free Survival (DFS) based on 2 subtypes (TCGA-CRC, Logrank test, n = 620). F The Sankey diagram completely showing the association between the subtypes and clinicopathological attributes. G Representative images of pathological Hematoxylin–eosin (HE) staining of 2 CC phenotypes (above, scale bars = 500 μm; below, scale bars = 50 μm). C2 contained a more abundant matrix component than C1. H Violin plots showing the immune score and stromal score of different CC patterns (Wilcoxon test). I The box plot indicating the Tumor Immune Dysfunction and Exclusion (TIDE) score of different CC patterns (Wilcoxon test). J Comparison of TME infiltrating cells between the two CC phenotypes (Wilcoxon test). ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05

Identification of LMOD1 as CC phenotype-associated factors

Differences in the gene expression between CC subtypes were assessed to identify potential regulators. A total of 157 DEGs were identified between C1 and C2 (Fig. S9A). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses revealed that the extracellular matrix (ECM) organization and collagen-comprising ECM, ECM structural constituent, and other collagen-matrigel matrix-related signals were significantly activated in C2 (Fig. S9B, Table S12). For hallmark enrichment, C2 was enriched for a large set of cancer-related genes, such as PI3K_AKT_MTOR_SIGNALING, MTORC1 SIGNALING, WNT_BETA_CATRNIN_SIGNALING, and NOTCH SIGNLING (Figure S9C). Afterward, a MEGENA network was established based on the DEGs (Figure S9D), and 8 modules and 157 module genes were obtained, with the largest module C1_5 encompassing 46 genes, followed by C1_2 with 33 genes and module C1_11 with 27 genes (Fig. S9E, Table S13). A Cox multivariable proportional hazards model constructed based on hub genes indicated that COL1A2, LMOD1, and MYH11 were significant adverse factors for CRC (Fig. S9F–H, TCGA-CRC, PFI, P < 0.05). Considering the very limited number of studies on LMOD1 in CRC, the current study focused on LOMD1, and Kaplan–Meier curves confirmed that LMOD1 was associated with a worse DSS and PFI in CRC (Fig. S9I–K). Interestingly, there was a low expression of LMOD1 in the CRC samples compared with the controls (Fig. S9L, P < 0.001). We further analyzed the difference in LMOD1 gene expression among different clinical stages of CRC cases in the TCGA-CRC cohort and found a positive correlation (Figure S9M–P; T stage, P < 0.05; N stage, P < 0.01; M stage, P < 0.05; Pathologic stage, P < 0.05). Next, a network of LMOD1-related genes was generated in the GeneMANIA (http://www.genemania.org) and found a close link between LMOD1 and mesenchymal-related genes, such as ACTA2, ITGA1, and ACTG2 (Figure S9Q). Overall, these results indicated that LMOD1 may be an important oncogene.

LMOD1 resembles TJs in the epithelium and GJs in the stroma

Furthermore, we explored the possible connection between LMOD1 and CC features. Trajectory analysis (our scRNA data) showed that LMOD1 was not involved in the transformation of normal epithelium into adenomas (Figs. 5A and S10A). Meanwhile, for the malignant transformation of ECs, the sequence of LMOD1 expression changes was similar to that of TJs; for malignant transformation of fibroblasts, the sequence of LMOD1 expression changes was in line with that of GJs and indicators of CAFs activation (Figs. 5B, C, S10B, and S11). This interesting pattern is summarized in Fig. 5D. mIF staining of independent samples confirmed that for the epithelium, LMOD1 expression was lower in cancer tissues than in healthy tissues (Fig. 5E–H); the opposite result was observed for the stroma (Fig. 5I–L). In addition, we described the spatial distribution of LMOD1 in the parenchyma and stroma of CRC tissues in the ST data (GEO database). As expected, LMOD1 expression disappeared in the parenchymal component of the tumor (compared with normal parenchyma), whereas it was relatively increased in the tumor stromal component (Fig. 5M, N), as summarized in Fig. 5O. Subsequently, the Tumor Immune Estimation Resource (TIMER) tool and “ESTIMATE” package were employed for further quantification in TCGA-CRC. The results revealed a marked positive association between LMOD1 expression and CAF infiltration level (colon adenocarcinoma [COAD], R = 0.719, P = 5.96e−45; rectum adenocarcinoma [READ], R = 0.516, P = 2.00e−07), and a significant negative correlation with tumor purity (R = − 0.354, P = 1.98e−13) (Fig. 5P). As depicted in Fig. 5Q, LMOD1 expression was significantly positively correlated with all the scores, especially the stromal score (R = 0.749, P < 0.001). Subsequently, we selected ACTA2 (Fig. 5R, TCGA-CRC, R = 0.884, P < 0.001) and fibroblast activation protein (FAP, Fig. 5S, TCGA-CRC, R = 0.677, P < 0.001) as markers of CAF activation. Human primary CAFs were successfully isolated, and it was confirmed in vitro that LMOD1 overexpression significantly promoted the expression of ACTA2 and FAP, whereas down-regulation of LMOD1 demonstrated the opposite result (Fig. 5T, U, P < 0.001).

Fig. 5

Transcriptional landscape heterogeneity of LMOD1. A Two-dimensional plots illustrating the invariabilities in LMOD1 expression during the transitions (from normal cells to adenoma cells) along the pseudotime. B Two-dimensional plots showing the variations (decrease) in LMOD1 expression during the transitions (from normal/adenoma cells to cancer cells) along the pseudotime. C Two-dimensional plots indicating the variations (increase) in LMOD1 expression and fibroblasts activation markers during the transitions (from NFs to CAFs) along the pseudotime. D LMOD1 is involved in transdifferentiation from normal epithelium to colorectal cancer but not from normal epithelium to intestine adenomas (grey⊥stands for no change; purple↓stands for decrease). E, F Double immunofluorescence (dIF) staining images of EPCAM (green) and LMOD1 (red) in a resected normal colon specimen (E) and a resected colon cancer specimen (F) (blue, DAPI). LMOD1 is upregulated in normal epithelial cells but downregulate in cancer cells. Scale bars are provided on the graph. G, H Co-localization was determined using the Pearson correlation coefficient in normal colon specimen (G, R = 0.5367, P < 0.0001) and colon cancer specimen (H, R = − 0.5860, P < 0.0001), respectively. The X-axis represents each pixel point on the image, and the Y-axis represents the gray value corresponding to each pixel point. The co-localization relationship between LMOD1 and EPCAM was weaker in tumor tissue compared to that of normal tissue, just like TJs. I dIF staining images of ACTA2 (green) and LMOD1 (red) in a resected normal colon specimen (blue = DAPI), and LMOD1 downregulated in NFs. Scale bars are labelled on the graph. J mIF staining images of ACTA2 (green), LMOD1 (red), and CK (pink) in a resected colon cancer specimen (blue = DAPI), and LMOD1 is upregulated in cancer mesenchymal tissues. Scale bars are provided on the graph. K, L Co-localization was determined using the Pearson correlation coefficient in normal colon specimen (K, R = 0.7930, P < 0.0001) and colon cancer specimen (L, R = 0.9290, P < 0.0001), respectively. The co-localization relationship between LMOD1 and ACTA2 was stronger in tumor tissue compared with that in normal tissue, simile to GJs. The X-axis represents each pixel point on the image, and the Y-axis represents the gray value corresponding to each pixel point. M, N Spatial transcription sections showing the spatial expression of LMOD1 in normal colonic tissue (M) and colon cancer tissue (N). The dot color represents the expression level of the markers. Green boxes for the parenchyma and pink boxes for the mesenchyme. O LMOD1 exhibited similar behaviors as TJs during the malignant transformation of epithelial cells and to GJs in malignant transformation of fibroblasts. P Spearman correlation between LMOD1 expression and the tumor purity (left) as well as infiltration level of fibroblasts in Colon adenocarcinoma (COAD) (middle) and Rectum adenocarcinoma (READ) (right) was analyzed on TIMER 2.0 (TCGA-CRC). Q Spearman association of LMOD1 expression with stromal score (left), immune score (middle), and estimate score (right) was analyzed by “ESTIMATE” package (TCGA-CRC). R, S Spearman association of LMOD1 with fibroblast activation markers ACTA2 (R) and FAP (S) expression (TCGA-CRC). T, U Double staining technique by ACTA2 (green), and FAP (red) staining in the primary CAFs. Representative images of staining are shown. All assays were conducted thrice, independently (scale bars = 20 µm). ANOVA was applied. ****P < 0.0001

Fibroblast-expressing LMOD1 promotes cancer invasion through FGF1

Since LMOD1 may play distinct roles in the epithelium and stroma, we investigated LMOD1 function in fibroblasts and ECs, respectively. A correlation assay indicated that LMOD1 was positively correlated with most fibroblast growth factors (FGFs) (Fig. 6A, TCGA-CRC), especially FGF1 (Fig. 6B, TCGA-CRC, R = 0.730, P < 0.001). FGF1 has been shown to play an important role in the expansion of CAFs by suppressing the transcription of tumor protein p53 (TP53) and triggering tumor fibrosis [55]. Importantly, CAFs form signaling crosstalk with cancer cells to support tumor growth. NFs co-evolve with cancer cells and transform into CAFs, and CAFs secrete cytokines to induce cancer cell survival [56]. Since we confirmed that LMOD1 promotes the activation of CAFs, we further investigated whether LMOD1 could stimulate cancer cell invasion through the above mechanism. A CAF-cancer cell co-culture model was established as shown in Fig. 6C. Under co-culture with oe-LMOD1-treated CAFs, RKO and SW480 cells exhibited more potent migration ability, and this effect was partially reversed by si-FGF1 (Fig. 6D, E, wound healing, P < 0.0001). Next, transwell experiments were performed in the model (Fig. 6F), which showed that LMOD1-treated CAFs affected the invasion ability of cancer cells (Fig. 6G, H, P < 0.01). Correlation analysis revealed that LMOD1 may have a regulatory function on EMT-related markers (Fig. 6I: TCGA-CRC; CDH2, R = 0.750, P < 0.001; CDH1, R = 0.081, P = 0.040; MMP9, R = 0.550, P < 0.001; MMP2, R = 0.713, P < 0.001; SNAI1, R = 0.386, P < 0.001; SNAI2, R = 0.601, P < 0.001). The WB assay confirmed that LMOD1 up-regulation increased the levels of N-cadherin, matrix metalloproteinase 2 (MMP2), MMP9, Slug, and Snail in RKO and SW480 cells (co-cultured with CAFs in the pattern of Fig. 6C) but reduced the expression of E-cadherin; however, this effect was abolished by ectopic FGF1 expression (Fig. 6J, K). Together, these results suggest that LMOD1, localized on fibroblasts, is a proto-oncogene.

Fig. 6

LMOD1/FGF1 in CAFs promotes CRC cell invasion and metastasis by regulating the EMT process. A Correlation analysis between LMOD1 and fibroblast growth factors (FGFs) based on TCGA-CRC. B Pearsons’s correlation coefficient between LMOD1 and FGF1 expression based on TCGA data. C Non-contact co-culture unit of CAFs and CRC cells for cell migration (Wound healing assay). A 1:1 ratio of the cells was employed. D, E Cell migration (Wound healing assay) in RKO and SW480 cells with the intervention of CAFs over-/under-expressing LMOD1 in the presence or absence of si-FGF1. (Scale bars = 100 μm, magnification, × 200). F The cell invasion (Transwell assay) assay for non-contact co-culture unit of CAFs and CRC cells. A 1:1 ratio of the cells was employed. G, H Cell invasion (Transwell assay) for RKO and SW480 cells with the intervention of CAFs over-/under-expressing LMOD1 in the presence or absence of si-FGF1. (Scale bars = 100 μm, magnification, × 200). I Association of LMOD1