The CRC Single‐Cell RNA Sequencing (scRNA-seq) data were obtained from the GEO website (https://www.ncbi.nlm.nih.gov/geo/), GSE number was GSE161277. The dataset includes patient-matched blood, normal, para-cancer, polyp (adenoma), and cancer tissue. As we used publicly available data from databases, we did not need to assign an ethical approval number for this study.

Gene expression data was imported into the R software (version 4.2.0) and analyzed using the Seurat R package (version 4.3.0) (40). Low-quality cells were rejected: (1) 300 < nFeature < 7,500; (2) 500 < nCount < 100,000; (3) mitochondrial gene expression not exceeding 25% of the total number of genes in the cell; (4) erythrocyte gene expression not exceeding 5% of the total number of genes in the cell.

The harmony R package (v.0.1.1) (41, 42) was used to eliminate the batch effect between samples (43–45). The dimensionality reduction parameter (dim value) was set to 30 and the resolution parameter was set to 1.2. All gene expression data were normalized and scaled using the NormalizeData and ScaleData functions in the Seurat R package (v.3.1.4) (46). The “FindVariableFeautres” function (47) was used to select the top 2,000 most variable genes for PCA (48–51). Cells were divided into clusters based on the top 30 principal components (PCs) using FindClusters at a resolution of 1.0. The top 30 essential PCs were selected for uniform manifold approximation and projection (UMAP) downscaling and gene expression visualization (52, 53).

We conducted differential expression analysis of all genes in cell clusters using Seurat’s FindAllMarkers function to identify the marker genes in each cluster (54–56). An adjustment of P-value< 0.05, expression percentage > 0.25, and | log2 [fold change (FC)] | > 0.25 were taken as threshold standards for identifying marker genes. Afterwards, different cell clusters were identified and annotated by the singleR package based on the composition patterns of the marker genes and were then manually confirmed and corrected with the CellMarker database.

To assess the predilection of TCs subtypes for cancer, odds ratios (OR) were calculated using the following computational methods (57).

AUCell (58) is a method for identifying cells with active genes in scRNA-seq data. It takes a gene set as input and outputs the “activity” of that gene set in each cell. It was used in this paper to rate the stemness of subtypes of TCs. Cell stemness was assessed through the utilization of the CytoTRACE R package (version 0.3.3) (59), allowing the temporal order of cell differentiation to be inferred (60, 61).

Single-cell pseudotime trajectories were constructed using the Monocle 2 algorithm (53, 62), an R package designed for single-cell trajectories by Qiu et al (63). This algorithm employs a machine learning technique to narrow down the high-dimensional expression spectrum to a low-dimensional space, visualized as a UMAP plot. Single cells were cast into this space and sorted into a trajectory with branching points. A dynamic expression heatmap was constructed using the plot_pseudotime_heatmap function. Furthermore, differential expression analysis between branches was performed using the plot_genes_branched_heatmap function.

The slingshot R package (version 2.6.0) was used to infer cell lineages and pseudotimes (64). It defined lineage structure using clustering-based minimum spanning trees (MSTs) and applications synchronized master curves to fit branching curves for these lineages. The getCurves function was used to obtain smooth trajectory curves.

Using the Gene Ontology (GO) (52, 65–67), Genome Enrichment Analysis (GSEA) (http://software.broadinstitute.org/gsea/msigdb) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (68) tools in the ClusterProfiler R package (version 4.6.0) (63), differentially expressed genes (DEGs) were enriched and analyzed. Significance of GO terms was determined based on an adjusted P < 0.05 (69, 70).

To investigate the top 5 transcription factors (TFs) with the most prominent expression changes in each TCs subtypes, single-cell regulatory network inference and clustering analyses were performed using the pySCENIC R package (version 0.10.0) in Python (version 3.7). Initially, GRNBoost was leveraged to identify potential target genes for each TFs. Then, DNA-motif analysis helped to determine potential direct binding targets. Lastly, the activity of regulon in the cells was scored by AUCell, and the top 5 TFs with the highest score were identified. The human gene-motif rankings were derived from https://resources.aertslab.org/cistarget/. Regulon modules were then identified according to the Connection Specificity Index (CSI) to confirm specific associating partners (71). Hierarchical clustering with Euclidean distance was then conducted to identify different regulon modules. We then constructed a regulator association network using 0.65 as a cutoff to investigate the relationship between different regulators.

The CellChat R package (version 1.6.1) was utilised to infer sophisticated intercellular interactions and construct regulatory networks based on ligand-receptor levels (72). It utilizes the “netVisualDiffInteraction” function to indicate differences in the strength of intercellular communication and the “IdentifyCommunicationPatterns” function to indicate an estimate of the number of communication patterns. We used a significance threshold setting with a P-value cutoff of 0.05 to predict cell-cell interactions between different cell types.

To elucidate the impact of CRC-associated TCs on predicting patient survival, we employed marker genes of critical tumor cell subtypes as predictive factors for the construction of a prognostic model. Employing univariate Cox analysis and LASSO regression analysis (73–75), we identified the most significant prognostic genes. Multivariate Cox regression analysis was then conducted to calculate the risk coefficients for each prognostic gene (76–78), enabling the establishment of a risk score model:

Risk score = ∑ i n Xi × Yi , where X represents the coefficient and Y represents the gene expression level.

Patients were categorized into low-risk and high-risk groups based on the optimal cutoff value calculated using the “surv_cutpoint” function. To observe the prognostic outcomes in different patient cohorts, we further utilized the R package “Survival” (version 3.3.1) for survival analysis of the constructed risk score model and visualized survival curves using the “ggsurvplot” function. The accuracy and calibration of the predictive model were evaluated and calibrated using the “timeROC” package (version 0.4.0) to plot ROC curves. The AUC values for the time-dependent ROC curve were computed at different time intervals to evaluate the model’s predictive accuracy over time (79).

CIBERSORT R package is a machine learning approach, which estimates the abundance of cell clusters in bulk RNA-seq data (80). We used this tool to digitally purify the transcriptomes of individual cell clusters from bulk data without isolating single cells. We extracted the transcripts per million (TPM) normalization datasets of selected cell types including CD8+ T cells, fibroblasts, myeloid cells and epithelial cells to create the signature matrix in 1000 permutations and without batch correction. We then divided the CRC patients available in the TCGA database into training and testing cohorts at a 1:1 ratio using survival-based randomisation and estimated the proportion of each cell cluster in the training and testing cohorts, respectively, using CIBERSORT. Notably, the bulk RNA-seq data from TCGA was initially normalized to TPM values.

NIC-H716 and SW837 cell lines were obtained from American Type Culture Collection (ATCC). NIC-H716 cell line was cultured in RPMI 1640 medium containing 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin under standard conditions (37°C, 5% CO2, 95% humidity). SW837 cell line was cultured in RPMI 1640 medium containing 10% fetal bovine serum (FBS), 1% penicillin-streptomycin under standard conditions (37°C, 5% CO2, 95% humidity) (81).

DLX2 knockdown was achieved using small interfering RNA (siRNA) constructs from GenePharma (Suzhou, China). Transfection was performed following the manufacturer’s instructions of Lipofectamine 3000 RNAiMAX (Invitrogen, USA). Cells were seeded at 50% confluency in 6-well plates and transfected with negative control (si-NC) and knockdown constructs (Si-DLX2-1 and Si-DLX2-2). Lipofectamine 3000 RNAiMAX (Invitrogen, USA) was used for each transfection.

Cell viability of transfected NIC-H716 and SW837 cells was assessed using CCK-8 assay. Cells were seeded in 96-well plates at a density of 5×10³ cells per well and cultured for 24 hours. After adding 10μL of CCK-8 reagent (A311-01, Vazyme) per well, plates were incubated at 37°C in the dark for 2 hours. Absorbance at 450nm was measured using a microplate reader (A33978, Thermo) on days 1, 2, 3, and 4 post-transfection. Mean OD values were calculated and plotted (82).

Stably transfected cells were seeded in 6-well plates and grown to confluence. A sterile 200μL pipette tip was used to scratch each well, which was then washed with PBS to remove cell debris and incubated in serum-free medium. Images of the scratches were captured at 0 hours and 48 hours, and the width of the scratch was measured using Image-J software (21).

Cells were starved in serum-free medium for 24 hours prior to the experiment. Following treatment with Matrigel (BD Biosciences, USA), cell suspension was added to the upper chamber of Costar plates, with serum-containing medium in the lower chamber. Cells were incubated for 48 hours in a cell culture incubator. Following cultivation, cells were fixed with 4% paraformaldehyde and stained with crystal violet to assess their invasive capacity.

Transfected NIC-H716 and SW837 cells were seeded at 5×10³ cells per well in 6-well plates and cultured overnight. Subsequently, a 2x EdU working solution was prepared by mixing 10 mM EdU with serum-free medium. After a 2-hour incubation at 37°C, cells were washed with PBS, fixed with 4% paraformaldehyde for 30 minutes, permeabilized with 2 mg/mL glycine and 0.5% Triton X-100 for 15 minutes, and then stained with a mixture of 1X Apollo and 1X Hoechst 33342 for 30 minutes at room temperature. Cell proliferation was quantified using fluorescence microscopy.

To identify the major cell types during CRC progression, we collected 13 tissue samples containing adenoma, blood, carcinom, normal and para-cancer from 4 CRC patients on GEO for scRNA-seq. In addition, we checked the quality and completeness of the raw data. We then excluded genes in the samples that did not reach the minimum expression threshold. After initial quality control and removal of batch effects, we maintained a total of 42,133 cells ( Figure 1A ). We classified these 42,133 cells into 32 cell clusters by dimensionality reduction clustering. Based on cellular genetic profiles and typical markers, we finally identified the 32 cell clusters into 8 cell types, including T and NK cells (CD3D), fibroblasts (DCN), mast cells (MCs, TPSB2), epithelial cells (EPCs, PHGR1), B cells (CD74), plasma cells (JCHAIN), proliferating cells (TUBA1B) and myeloid cells (LYZ). From the UMAP plots and the bar plots we could learn that from normal to para-cancer to adenoma and carcinoma, the proportion of EPCs and proliferating cells roughly increases gradually, in contrast to a gradual decrease in the proportion of T and NK cells, and the proportion of B cells from normal to para-cancer to adenoma was gradually increasing but decreased in carcinoma; for cell cycle, most EPCs were in the G1 phase, which was the pre-DNA synthesis period with active material metabolism, aiming to prepare the material and energy for DNA replication in the next stage, the S phase (83), whereas T and NK cells were evenly distributed in all three periods, and proliferating cells were only present in S and G2M phases, synthesizing DNA as well as making sufficient preparations for mitosis, in line with their proliferative properties ( Figures 1B–D ). The bubble plot showed us the first five markers used to distinguish cell types, and the results were consistent with immunology ( Figure 1E ). The bar graphs showed the expression levels of nCount-RNA, nFeature-RNA, G2M.Score and S.Score for different cell types, tissue sources and cell cycle, respectively. From the results, it could be seen that proliferating cells had the highest expression levels of all four items, while EPCs, fibroblasts and plasma cells ranked behind them. In the tissue source, from normal to para-cancer to adenoma to carcinoma, the expression of both nCount-RNA and nFeature-RNA gradually increased, and there were statistically significant differences between all groups ( Figures 1F–I ).

Among the above cell types, EPCs attracted our attention. EPCs not only have relatively abundant RNA counts, but also tend to aggregate gradually during the developmental stages of cancer, whereas the tumor cells of CRC originate from the intestinal epithelial cells of the colon and rectum (84). The tumor cells are the most fundamental cause of tumors, so we next turned our attention to the study of CRC’s tumor cells.

Next, we identified malignant tumor cells from epithelial cells using inferCNV ( Supplementary Figure S1 ). Due to the absence of blood vessels in the epithelial tissues, we removed the blood tissues and compiled scRNA-seq data from the remaining four tissues in the CRC for further subclustering analysis against tumor cells. The analysis yielded seven clusters, resulting in seven different tumor cell subtypes: C0 FXYD5+ TCs, C1 APCDD1+ TCs, C2 MUC2+ TCs, C3 HEPACAM2+ TCs, C4 OTOP2+ TCs, C5 SLC26A3+ TCs, and C6 AVIL+ TCs, and showed the distribution of the various tissue origins and cellular period distribution in the subtypes ( Figure 2A ). Ro/e preference plots combined with box line plots showed that C6 AVIL+ TCs preferred the G2M phase, while C1 APCDD1+ TCs preferred the S phase, and were favored by adenoma tissues, whereas in contrast, C2 MUC2+ TCs, C4 OTOP2+ TCs, and C5 SLC26A3+ TCs were found overwhelmingly in the normal tissue and para-cancer tissues, whereas C0 FXYD5+ TCs preferred carcinoma tissues ( Figures 2B, C ). The above tissue distribution features suggested that C1 APCDD1+ TCs may be a watershed for colorectal progression from adenoma to malignancy, whereas C0 FXYD5+ TCs was a subtype that had become cancerous. The top five marker genes obtained by enrichment in different TCs subtypes were next visualized by bubble plots. We could see that S100A11 was significantly enriched in C0 FXYD5+ TCs, whereas overexpression of S100A11 in cytoplasmic and nuclear subcellular compartments was associated with tumor metastasis and poor prognosis in patients with CRC and meta-analyses demonstrated that the expression level of S100A11 in CRC tissues significantly elevated in CRC tissues, a marker gene for clinical use in identifying CRC (85, 86). Meanwhile, SELENBP1 was also noteworthy, as it was highly expressed in C1 APCDD1+ TCs, and studies had shown that it was down-regulated in CRC but not in normal tissue or adenoma tissue (87), a typical oncogene that inhibits colorectal cancer progression by suppressing the EMT (88), whereas it was low-expressed in C0 FXYD5+ TCs, reflecting its intermediary role in the progression from C1 APCDD1+ TCs to C0 FXYD5+ TCs ( Figure 2D ). From the UMAP plots and bar graphs, we could tell that the C0 FXYD5+ TCs and C1 APCDD1+ TCs subtypes had an abundant number of RNA, as well as a higher CNVscore for C0 FXYD5+ TCs and C2 MUC2+ TCs ( Figures 2E–G ).

Finally, the violin plots showed the marker genes of each TCs subtype in turn ( Figures 2H–N ), which had shown corresponding expression advantages in their respective subtypes. FXYD5 in C0 FXYD5+ TCs deserved our attention, as it had been shown to activate TGF-β/SMAD signaling and drive EMT to promote ovarian cancer progression, and to form a positive feedback loop to drive EMT during the progression of ovarian cancer and to promote tumor growth and metastasis (89); meanwhile, recent study also showed that FXYD5 was upregulated in various tumor types and positively correlates with tumor progression (90).

To explore the expression of stemness genes in TCs subtypes and to comprehend their proliferative differentiation potential, we used bubble plot to illustrate the different expression of stemness genes in TCs subtypes ( Figure 3A ). The differential expression analysis was intended to provide a holistic view of stemness-related gene expression, including both high and low-expression markers, to fully capture the diversity of stemness-associated features across tumor cell subtypes. Results showed the corresponding expression of stemness genes KDM5B, EPAS1, CTNNB1, MYC, KLF4 and CD44 in TCs subtypes and different tissue sources, with MYC highly expressed in C0 FXYD5+ TCs and PROM1 in C1 APCDD1+ TCs. Following this, UMAP plot showed the AUC scores of cells stemness in different tissue sources, which were significantly higher in adenoma and carcinoma than in normal and para-cancer ( Figure 3B ). Next, bar graphs showed the AUC scores of stemness genes in different TCs subtypes and in the cell cycle, respectively ( Figures 3C, D ). We next enumerated the stemness genes enriched in different TCs subtypes and different tissue sources using bubble plots, respectively, and the results remained that MYC was significantly enriched in C0 FXYD5+ TCs, and for the tissues, PROM1, NOTCH1, and CD44 were enriched in adenoma, whereas KDM5B and MYC were enriched in carcinoma. Stemness genes were also enriched in normal and para-cancer, but at lower levels, which suggested that adenoma and carcinoma tissues had a higher proliferative potential ( Figure 3E ). Finally, we visualized and compared the stemness genes with higher expression in C0 FXYD5+ TCs and C1 APCDD1+ TCs using box line plots, and the results further validated the previous results, showing a stronger proliferative capacity in C0 FXYD5+ TCs and C1 APCDD1+ TCs as well as adenoma and carcinoma ( Figures 3F–I ).

Based on the above results, in order to further clarify the differentiation status of TCs subtypes, we first used CytoTRACE to clustered all tumor cells, and obtained UMAP plots and box line plots of cell stemness prediction ordering ( Figures 4A, B ). From these three plots, we could learn that the C0 FXYD5+ TCs and C1 APCDD1+ TCs subtypes were in a high stemness and low differentiation state, which was consistent with the above results ( Figure 3 ). Next, the bar graphs showed the genes that correlated with the CytoTRACE results, where greater than 0 was a positive correlation and vice versa was a negative correlation ( Figure 4C ). We then performed Monocle analysis of tumor cells to obtain a pseudotime atlas of tumor cell subtypes, colored by the order of pseudotime, different subtypes, and different tissue sources ( Figures 4D–G ). First, the pseudotime order speculated that the developmental trajectory started at the bottom left corner and branches out to two more trajectories when it reached state 1. Also in conjunction with the subtypes of TCs, C2 MUC2+ TCs, C4 OTOP2+ TCs, and C5 SLC26A3+ TCs subtypes were concentrated in early developmental trajectories, whereas C1 APCDD1+ TCs, C3 HEPACAM2+ TCs, and C6 AVIL+ TCs were concentrated in late developmental trajectories, and unlike C0 FXYD5+ TCs subtype, it was found in mid- to late-developmental trajectories; for the tissue origin, normal and para-cancer tissues were present at the initial developmental stage, which later i.e., developed into the adenoma and carcinoma.

Furthermore, we then performed Slingshot analysis on tumor cells, which resulted in a pseudotime speculative map colored by pseudotime order and based on the UMAP plot obtained earlier ( Figure 4H ). Combined with the ridge plot, the analysis showed that, again, C2 MUC2+ TCs, C4 OTOP2+ TCs, and C5 SLC26A3+ TCs subtypes were at the anterior end of the developmental trajectory, while the rest of the subtypes were located in the middle to posterior ( Figure 4I ). By observing how Slingshot’s spectral differentiation trajectories were plotted, we could learn that tumor cells go from state 1 to state 3 to state 2; at the same time, this result yielded four spectral trajectories based on cell subtypes, as follows: Lineage1: C5 to C2 to C1 to C0; Lineage2: C5 to C2 to C1 to C3; Lineage3: C5 to C2 to C1 to C6; Lineage4: C5 to C2 to C4 ( Figures 4J, K ). The differences between the four trajectories were mainly in the late stage, with lineage1 ending at C0, lineage2 at C3, lineage3 at C6, and lineage4 at C4, where C0 not only had a preference for, but also a high percentage of tumor tissue. Hence, we predicted that lineage1 represented the lineage of differentiation of TCs associated with CRC tumorigenesis. Furthermore, we noted that lineage1 passed through C1, which was predominantly organized as adenoma, a precancerous, benign tumor, and previous studies had reported that participants with positive diagnostic colonoscopy results and advanced adenomas had a significantly increased risk of colorectal cancer compared to participants without adenomas (91). Therefore, studying how C1 progresses to C0 and the key factors that drive the epithelium to change from benign to malignant has important implications for how to stop CRC from occurring and worsening. Finally, the dynamic trends plots displayed the expression changes and distribution of marker genes for TCs subtypes across the four differentiation trajectories in pseudotime ( Figure 4L ).

Then, we performed GO-BP DEGs in TCs subtypes to reveal their enrichment in biological processes. Heatmap displayed the results of the top five enriched items in the seven TCs subtypes ( Figure 5A ). The C0 subtype was mainly associated with pathways such as cytoplasmic translation, aerobic respiration, oxidative phosphorylation, cellular respiration and energy derivation by oxidation of organic compounds; The C1 subtype was enriched in pathways such as cytoplasmic translation, ribonucleoprotein complex biogenesis, ribosome biogenesis, ribosome assembly and ribonucleoprotein complex; The C2 subtype revealed their close association with antigen processing and presentation of endogenous antigen, endogenous peptide antigen or peptide antigen; The C3 subtype showed enrichment in pathways such as protein targeting to endoplasmic reticulum (ER), establishment of protein localization to endoplasmic reticulum, cotranslational protein targeting to membrane and response to endoplasmic reticulum stress; The C4 subtype was enriched in pathways related to stress response to metal and copper ion, detoxification of copper ion and inorganic compound and regulation of cell morphogenesis; The C5 subtype mainly exhibited enrichment in pathways such as alcohol metabolic process, viral life cycle, negative regulation of response to external stimulus, cellular hormone metabolic process and primary alcohol metabolic process; The C6 subtype revealed pathways related to actin filament organization, regulation of actin filament-based process, supramolecular fiber organization and actin cytoskeleton organization and mononuclear cell differentiation. The enrichment items derived from the C0 subtype caught our attention, and the last four items were all related to cellular energy metabolism, which corresponds to the increased energy metabolism that occurred after cancerous transformation, e.g., in order to maintain proliferative and metastatic capacity, tumor cells increased their energy-producing pathways (92), and which we conjecture mediate cytoplasmic and mitochondrial translational homeostasis in order to modulate the respiratory capacity of the mitochondria to provide the tumor cells to provide energy to adapt to different environmental stresses and growth demands (93). This implied that C0 subtype was strongly correlated with tumor tissue development and progression. Next, we utilized volcano plots to demonstrate the spectrum of DEGs among TCs subtypes ( Figure 5B ). The results showed that the main DEGs upregulated in C0 subtype were TMSB10, S100A10, CLIC1, GSTP1 and RPS21. Among them, TMSB10, S100A10, CLIC1 and GSTP1 have been shown in previous studies to have an oncogenic mechanism in the colorectum as well as to contribute to the poor prognosis of colorectal cancer, and TMSB10 could be used as a minimally invasive serum tumor marker for the detection of CRC, while CLIC1 and GSTP1 may be potential prognostic biomarkers for patients with CRC (94–97).

In addition, the bubble plot showed the results of GSEA ( Figure 5C ). The results showed that the C0 subtype was significantly expressed in the pathways of translation, peptide biosynthetic process, amide biosynthetic process and ribonucleoprotein complex biogenesis. All these pathways indicated that energy metabolism as well as ribose and protein synthesis were actively taking place in the C0 subtype, providing a material basis for the rapid proliferation of tumor cells. The word cloud map showed the pathways enriched in DEGs in the TCs subtypes ( Figures 5D–J ). Consistent with the above results, it indicated that the C0 subtype was undergoing enriched metabolic and oxidative reactions and substance synthesis. Finally, we performed a GSEA of DEGs in the C0 subtype according to GO-BP ( Figure 5K ). We observed that pathways associated with cytoplasmic translation, amide biosynthetic process, peptide biosynthetic process and ribonucleoprotein complex subunit organization, which were consistent with the previous enrichment items, were upregulated in the C0 subtype. While pathways associated with defense response to gram-positive bacterium, positive regulation of cyclase activity and lyase activity and regulation of cyclase activity were down-regulated in the C0 subtype. Where for positive regulation of lyase activity decreases, it may lead to increased levels of degradative enzymes, which in turn are able to degrade proteoglycan basement membrane components, promoting increased tumor cell detachment from the primary tumor and local invasion, resulting in a poor prognosis (98). Thus, the items enriched by various methods in the C0 subtype were suggestive that the C0 subtype might be transforming into tumors with a more aggressive and poorer prognosis, and that we could incorporate this feature into risk group of the tumor population.

TFs act directly on the genome to regulate gene transcription by binding to specific nucleotide sequences upstream of genes, thereby affecting the biological functions of cells. Here, we analyzed the gene regulatory network of the C0 subtype using scenic analysis. First, we performed a dimensionality reduction clustering of TCs based on the activities of regulatory factors and listed the results from different tissue sources ( Figures 6A, B ). Among them, the C0 subtype still mainly consisted of carcinoma tissues. Next, we observed the correlation between the TCs subtypes. Considering that different TFs can co-regulate the expression of certain genes, based on the connectivity specificity index (CSI), we categorized the TCs subtypes into three modules of regulatory factors: M1, M2 and M3 ( Figure 6C ). These three regulatory submodules showed TFs that may be cooperative with each other in gene regulation. As shown in the bar graphs, the M3 regulatory submodule may take a major regulatory role in the biological functions of the C0, C1, and C3 subtypes ( Figure 6D ). Meanwhile, the M3 regulatory submodule was also the dominant regulatory module for the C0 subtype ( Figure 6E ). By further analyzing the key regulators of different TCs subtypes, it was found that three of the five major regulators of the C0 subtype, ETV4, MYC and XBP1, were present in the M3 regulatory submodule ( Figures 6F, G ). The expression of TFs in the C0 subtype was labeled in the UMAP plot ( Figures 6H, I ). It is worth noting that previous studies have indicated that CEBPB regulates the bile acid receptor FXR to accelerate colon cancer progression by modulating aerobic glycolysis (99), corresponding to the analysis of its enrichment pathway; whereas, the expression of ETV4 was significantly correlated with the depth of infiltration, lymphovascular and venous infiltration, lymph node and distant metastasis, and the stage of lymph node metastasis of the pathologic tumors in terms of progression and recurrence, which suggests an unfavorable prognosis (100); while the expression of XBP1 could also promote the proliferation of colon cancer cells (101). These were all factors of poor prognosis of CRC.

In order to specifically elucidate complex cellular responses, we intended to explore intercellular relationships and ligand-receptor communication networks for a better understanding of intercellular interactions. Through analyses via CellChat, we preliminarily established intercellular communication networks involving a variety of cells, such as proliferating cells, fibroblasts, T and NK cells, and different subtypes of TCs ( Figures 7A, B ). We also established a communication network with C0 FXYD5+ TCs as the source cells and target cells ( Figures 7C, D ). After building the intercellular communication network using CellChat analysis, we accounted for the number of interactions (expressed in terms of the thickness of the connecting lines between the two types of cells) and the strength of the interactions (expressed in terms of the weight of the connecting lines, with thicker lines indicating stronger interactions). This method helped to quantify the complexity and strength of the communication pathways between different cell types in the network. From the intricate communication network diagram, we could see that C0 FXYD5+ TCs sent positive communication signals to almost all cells, while fibroblasts and proliferating cells are the cell types that give it the most responses. Next, we used the gene expression pattern analysis method provided by CellChat to examine how cells and signaling pathways interact with each other ( Figures 7E, F ). First, we assessed the relationship between the inferred potential communication patterns and the population of cells secreting signaling molecules to decipher the outgoing communication patterns. Through our analysis, we identified three distinct signaling patterns: pattern 1 (tumor cells and proliferating cells), pattern 2 (fibroblasts), and pattern 3 (myeloid cells, B cells, plasma cells, MCs, and T and NK cells). Whereas for the incoming communication mode, the distinct signaling patterns: pattern 1 (C1 APCDD1+ TCs, proliferating cells, fibroblasts, C0 FXYD5+ TCs, C6 AVIL+ TCs and C3 HEPACAM2+ TCs), pattern 2 (myeloid cells, B cells, T cells and NK cells, MCs and plasma cells), and pattern 3 (C2 MUC2+ TCs, C4 OTOP2+ TCs and C5 SLC26A3+ TCs). Next, the heatmaps showed the signaling molecules of various types of cells under outgoing signaling patterns and under incoming signaling patterns, respectively ( Figure 7G ). Finally, to identify the major afferent and efferent signals associated with the seven TCs subtypes, we quantified the ligand-receptor network using CellChat ( Figures 7H, I ). With this approach, we could predict the major outgoing signals from secretory cells (signal senders) that release various cytokines or ligands. In addition, we evaluated on which cell types act as target cells (signal receivers) and how ligand-receptor-mediated communication between different cell types contributes to the progression of CRC. This analysis helps to illustrate how receptors on these cells are targeted by ligands released by the same or other cell types. In summary, we visualized signaling molecules between different cell types using multiple communication analyses that provide support for describing cellular communication.

Since we focused on the C0 FXYD5+ TCs subtype, we identified 9 genes that could be used as prognostic features based on the top 100 marker genes of the FXYD5+ TCs subtype by univariate Cox regression analysis ( Figure 8A ), of which ETS2 and ATOH1 were protective factors (HR < 1), and the other genes were risk factors. To address the problem of multicollinearity among these genes, we further screened them using LASSO regression analysis and multivariate regression analysis and finally identified 8 genes associated with prognosis ( Figures 8B–D ). The one that caught our attention was DLX2, which was the most dangerous prognostic factor among them with the largest HR. To further investigate the impact of FXYD5 highly expressed TCs on CRC patients, we categorized patients in the TCGA cohort into high and low FTRS (FXYD5+ TCs risk score) groups based on the identified 8 prognostic genes. Subsequently, we created a Nomogram survival prediction model for OS in CRC patients using autonomous prognostic factors to predict the prognosis of CRC patients ( Figure 8E ). The model was validated and performed well in predicting the OS C index ( Figure 8F ). Subsequent prognostic gene versus risk correlation plots, as well as ridge and box line plots, illustrated the clear correlation of prognostic genes with risk and variability across risk groups ( Figures 8G–H ). It could be seen that the difference between the two groups for our gene of interest, DLX2, was statistically significant and was expressed higher in the high FTRS group than in the low FTRS group. Next, to assess the accuracy of risk scores in predicting 1-, 3- and 5-year survival in patients with CRC, we performed ROC curve analyses on the three cohorts. The results showed high predictive accuracy with AUC (1-year) = 0.641, AUC (3-years) = 0.658, and AUC (5-years) = 0.691 ( Figure 8I ). Moreover, calibration curves were plotted to show the consistency between predicted and observed values for 1-, 3- and 5-year OS in both the training and validation cohorts ( Figure 8J ), which showed good agreement. The PCA plot also provided further evidence of the discrete and differential character of the high and low FTRS groups ( Figure 8K ). With definitive evidence, we next analyzed the calculation of coefficient values for these prognostic genes ( Figure 8L ). Using Kaplan-Meier survival curves, we first confirmed the poor prognosis corresponding to high expression of DLX2, and then we further concluded that the survival outcome was worse in the high FTRS group ( Figures 8M, N ). DLX2 was demonstrated to be associated with EMT in osteosarcoma in previous studies and its increased expression was associated with advanced gastric adenocarcinoma, but no article defining the relationship between DLX2 and poor prognosis in colorectal cancer could be found yet (102, 103).

The CIBERSORT analysis revealed the immune infiltration in the tumor samples. With the heatmap, we could observe the estimated proportion of immune cells in the high and low FTRS groups ( Figure 8O ). From the box line plot, we could easily see that there was a statistically significant difference in immune cell expression between the different groups ( Figure 8P ). Also, we were able to find that the relatively high expression of immune cells in the high FTRS group was positively correlated with the risk score ( Figures 8Q, R ). For example, T cells regulatory (Treg) relatively high expression in the high FTRS group and positive correlation with risk score represented a poor prognosis. Tregs have been documented in malignant tumors, where they promote tumor progression by suppressing effective anti-tumor immunity and are associated with a poor prognosis in various types of human cancers (104). In addition, we used ESTIMATE to calculate the stromal score, immune score, and ESTIMATE score for the high FTRS group and low FTRS group, and found that all the scores of the high FTRS group were significantly higher than those of the low FTRS group, which indicated that the level of immune cell infiltration in the tumor samples of the high FTRS group was higher ( Figure 8S ). In addition, tumor purity was higher in the low FTRS group than in the high FTRS group ( Figure 8T ). Since low tumor purity was associated with extensive infiltration of stromal and immune cells, it was independently associated with shorter survival time and faster recurrence, and significantly correlated with mesenchymal, invasive, and metastatic phenotypes suggesting a poor prognosis that corresponds to the prognostic curve. Finally, TIDE is a tool used to comprehensively assess tumor immune escape mechanisms and predict response to immune checkpoint inhibitor (ICI) therapy, with higher TIDE scores correlating with poorer immune checkpoint inhibition therapy (105), and violin plot of high FTRS group scores being significantly greater than lower groups ( Figure 8U ).

Immediately following the previous section, we proceeded to show the differential expression patterns of the genes used to construct the model using heatmap ( Figure 9A ). In addition, the curve and scatter plots showed the differences in risk scores and survival results between the two risk groups, suggesting that the high FTRS group was associated with a poorer prognosis, corresponding to the previous section ( Figure 9B ).

Next, volcano plot showed up- and down-regulated DEGs ( Figure 9C ). The filtering criteria for fold change were based on a log2 fold change threshold of 2, which meant that only genes with a fold change greater than or equal to 4 (or less than or equal to -4) were considered as differentially expressed. In addition to this fold change criterion, a p-value cutoff of 0.05 was applied to further refine the selection of significant DEGs. We then used these DEGs to gain insight into the biological processes involved by employing various enrichment methods. All GO enrichment results and their terms were shown, with different colors representing different subclasses ( Figure 9D ). GO-CC enrichment analysis showed that DEGs were predominantly enriched in lamellar body and rough endoplasmic reticulum, the latter of which suggested that there may be abundant protein synthesis in high-risk groups ( Figure 9E ). In GO-MF enrichment analysis, we observed enrichment in hormone and receptor ligand activity, suggesting that there may be an abundance of communication ( Figure 9F ). As for GO-BP enrichment analysis, the items mainly including regulation of lipid localization, intestinal lipid absorption and lipid transport ( Figure 9G ). All of the above enriched items are associated with lipid metabolism, and the former has been richly studied for its association with cancer. Dysregulated lipid metabolism is one of the most prominent metabolic alterations in cancer, and cancer cells utilize lipid metabolism for energy, biofilm components, and signaling molecules required for proliferation, survival, invasion, metastasis, and response to tumor microenvironmental influences and cancer therapy (106, 107). For another enrichment method analysis KEGG, the pathways were mainly AMPK, adipocytokine and PPAR signaling pathway ( Figure 9H ). Among them, AMPK plays a crucial role in maintaining energy homeostasis and metabolism, and its aberrant activation can lead to pro-carcinogenic effects, which is a double-edged sword (108); recent studies also show that AMPK signaling can also reversibly regulate highly active MAPK signaling in cancer cells by phosphorylating its key component, the RAF/KSR family of kinases, which not only affects carcinogenesis, but also the outcome of targeted cancer therapy against MAPK signaling. And there was an experiment showing that adipocytokines were involved in the carcinogenesis process (109), and an adipocytokine, LCN-2, has been reported to have the ability to destroy the extracellular matrix, which may lead to cancer progression and metastatic spread (110). At the same time, previous article had shown that PPAR-d upregulation increases susceptibility to colon tumorigenesis, which may impact the development of strategies to molecularly target PPAR-d in cancer and non-cancer diseases (111). Through the GSEA results we found that up-regulated DEGs were mainly enriched in the processes of antigen processing and presentation of exogenous peptide antigen, collagen fibril organization, intermediate filament-based process and regulation of intestinal absorption, while down-regulated DEGs were enriched in xenobiotic glucuronidation, midgut development, regulation of chromosome condensation and uronic acid metabolic process ( Figure 9I ).

Finally, in order to study the sensitivity of the high and low FTRS groups to different drugs, we predicted the drug sensitivity of each patient based on the drug sensitivity data in the GDSC database using the “pRRophetic” R software package. The results showed that the drugs AZD6244, CGP.082996, and FTI.277 had lower IC50 values and higher drug sensitivities in tumor cells of the high FTRS group, whereas the other drugs had lower IC50 values in tumor cells of the low FTRS group ( Figures 9J–O ).

Through the previous construction of CRC prognostic model, we obtained the prognostically strongly related gene DLX2. To further investigate its role in CRC, we performed in vitro experiments using NCI-H716 and SW837 cell lines. First, we knocked down DLX2 and measured mRNA and protein expression levels before and after knockdown. Compared with the control group, we observed a significant decrease in both mRNA and protein expression levels in both cell lines ( Figure 10A ). Subsequently, CCK-8 assay showed that tumor cell viability was significantly decreased after DLX2 knockdown ( Figure 10B ). Colony formation assay confirmed that DLX2 gene knockdown inhibited tumor cell aggregation ( Figure 10C ). In addition, after DLX2 knockdown, the migratory and invasive abilities of tumor cells were assessed using the scratch and the transwell method, and the results showed that the levels of migration and invasion were significantly reduced ( Figures 10D, E ). Finally, EDU staining method also confirmed these results ( Figure 10F ). Together, these results indicated that knockdown of DLX2 inhibited tumor cell activity, proliferation, migration, and invasion, thereby suppressing tumor growth.