This study utilized scRNA-seq data sourced from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) with accession number GSE166555. The Cancer Genome Atlas (TCGA) portal (https://portal.gdc.cancer.gov/) was used to acquire bulk RNA-seq data. We imported the 10x Genomics data from each sample into R program (v4.3.3) via the Seurat package (v4.3.0). Initially, possible doublets and substandard cells were eliminated utilizing the DoubletFinder algorithm (v2.0.3). Cells were preserved for subsequent studies if they satisfied the following criteria: 500 < nFeature_RNA < 6000 and mitochondrial gene expression being less than 25% of total expression.

The FindAllMarkers program (29, 30) employed the Wilcoxon rank-sum test with default parameters (log fold change > 0.25) to find differentially expressed genes (DEGs) for each cell type and subpopulation. Enrichment analyses were conducted utilizing the clusterProfiler (v4.6.2) and Single Cell Profiler (SCP) (v0.4.8) packages to clarify the functional roles of DEGs within each cell type and subpopulation. All enrichment pathways were obtained from Gene Ontology (GO) (29–34). Furthermore, we utilized AUCell (35, 36) to discern active gene sets and transcription factors at the single-cell resolution.

The InferCNV method (v1.17.0) was employed to ascertain copy number variation (CNV) levels (37). Copy number karyotyping of aneuploid cells was employed during carcinogenesis to distinguish between non-malignant and malignant tumor cells. Endothelial cells served as the reference population for CNV comparison to ascertain whether other tumor cells displayed significant chromosomal copy number variations.

Pseudotime trajectories of CRC tumor cell subpopulations were recreated utilizing the Monocle software (v2.24.0), based on single-cell RNA sequencing data. This methodology simulates the dynamic evolution of individual cells and elucidates the transcriptional alterations linked to the differentiation of tumor cells in CRC progression.

We utilized CytoTRACE to analyze variations in developmental and differentiation states among CRC tumor cell subpopulations, inferring and ranking the differentiation capacity of all tumor cell clusters. Furthermore, the getlineage and getCurves functions were employed to deduce lineage trajectories and evaluate dynamic alterations in gene expression during pseudotime. The Slingshot software (v2.6.0) was utilized to generate developmental trajectories, offering insights into the differentiation status and advancement of each tumor cell subpopulation (38).

We utilized CytoTRACE2 for a comprehensive analysis of the scRNA-seq data, with the objective of predicting the potential subtypes and absolute developmental capacities of tumor cell subpopulations originating from tumors at various clinical stages. In CytoTRACE2, latent classes are delineated according to cellular developmental potential, offering a continuous spectrum from 0 (completely differentiated) to 1 (totipotent) to quantify developmental capacity. This facilitates direct comparisons among datasets in absolute developmental space. CytoTRACE primarily relies on gene expression diversity, while CytoTRACE2 employs an integrated model that combines multiple feature-optimized approaches. We therefore believe that the two methods are complementary, which helps enhance the stability and reliability of the inferred results. Moreover, CytoTRACE2 enables the detection and validation of differentiation status among all tumor cell subpopulations.

Utilizing scRNA-seq data, we employed the CellChat software (v1.6.1) to forecast probable intercellular interactions among several cell types, encompassing tumor cell subpopulations and additional stromal or immune cells. Ligand-receptor interactions were deduced with CellChatDB.human as the reference database. A significance level of P < 0.05 was utilized to discern substantial interactions across different cell types or tumor subpopulations. This research concentrated on visualizing and analyzing the interactions between major tumor subpopulations and relevant immune cells, such as macrophages and mast cells, to identify potential signaling pathways of relevance.

We conducted SCENIC analysis to reconstruct gene regulation networks from scRNA-seq data and to identify stable cellular states. Specifically, we utilized the pySCENIC module (v0.12.1) in Python (v3.9.19) to infer transcription factor activity and generate an AUCell matrix for evaluating regulon enrichment and activity. The outcomes were later visualized using R software (v4.3.3) to enhance the understanding of transcriptional regulation among cell groups.

We evaluated the prognostic significance of essential tumor cell subpopulations as survival indicators in colorectal cancer by utilizing their marker genes as potential predictive attributes. The survival R package was utilized to perform univariate Cox regression and LASSO regression studies to identify additional genes pertinent to prognosis. Using a multivariate Cox regression model, a prognostic signature was created, and the formula was applied to calculate each patient’s risk score: Risk score = (gene1 expression × coefficient1) + (gene2 expression × coefficient2) +… + (geneN expression × coefficientN).

Based on the median risk score, patients were divided into high-risk and low-risk groups, with scores above the median considered high risk and those below considered low risk. We conducted a Kaplan–Meier survival analysis to investigate the disparities in overall survival (OS) across the groups (39). To evaluate the model’s predictive accuracy and calibration, the timeROC package (v0.4) was employed to produce time-dependent ROC curves for 1, 3, and 5 years. A multivariate Cox regression analysis was performed to determine whether the risk score functioned as an independent prognostic factor. A nomogram was created to predict OS at 1, 3, and 5 years, with internal validation conducted through the concordance index (C-index) and calibration plots.

Transcriptomic data from colorectal cancer patients, together with extensive clinical information, were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) for subsequent study. Patients were categorized into two groups according to the expression levels of specific genes: the high MKI67⁺ TCs risk score group (High MTRS group) and the low MKI67⁺ TCs risk score group (Low MTRS group). Kaplan–Meier survival curves were constructed using the survival program (v3.5-5) and displayed with the survminer tool (v0.4.9) to assess survival outcome disparities between these groups.

We assessed the invasion of 22 immune cell types utilizing the CIBERSORT R program (v0.1.0). Following this, the CIBERSORT, ESTIMATE, and xCell algorithms were employed to evaluate the immune microenvironment in CRC patients, including differences in immune cell infiltration levels and the varied expression of genes related to immune checkpoints. Visual tools were created to show the relationships between immune cells, model genes, OS, and risk scores. The Tumor Immune Dysfunction and Exclusion (TIDE) platform (http://tide.dfci.harvard.edu) was utilized to forecast patients’ responses to immunotherapy. The forecasts of drug immune responses were then validated utilizing data from the TCIA database (https://www.cancerimagingarchive.net/).

Our study included therapeutically pertinent medications to assess drug sensitivity among patient cohorts. pRRophetic program (v0.5) was employed to determine the IC50 for each drug. To assess differences in drug sensitivity, the projected IC50 values for high-risk and low-risk groups were compared.

Procell Life Science & Technology Co., Ltd. (Wuhan, China) supplied the HCT116 and HT-29 human colorectal carcinoma cell lines. The cells were cultured in MEM supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 0.1 mg/mL streptomycin. All cultures were preserved in a humidified incubator at 37°C with 5% CO₂ concentration. Cells in the logarithmic growth phase were harvested for following experimental procedures.

Cells were seeded in 6-well plates at a density of 2 × 10⁵ cells per well for transfection assays. Following a 24-hour incubation period, cells were transfected with HMGA1-targeting small interfering RNAs (siRNAs) obtained from GenePharma (Shanghai, China) at a final concentration of 20 µM. Transfection was performed via Lipofectamine RNAiMAX (Life Technologies, Thermo Fisher Scientific, Brendale, QLD, Australia) in accordance with the manufacturer’s instructions. Cells were collected 24 hours after transfection for further analysis. The sequences of the HMGA1 siRNAs employed were as follows: siRNA1, ACUCCAGGAAGGAAACCAA; siRNA2, AGCGAAGUGCCAACACCUA.

Total RNA was extracted utilizing TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) in accordance with the manufacturer’s guidelines. Subsequently, 500 ng of total RNA was reverse transcribed into cDNA utilizing the PrimeScript™ RT Reagent Kit (TaKaRa, Tokyo, Japan). Quantitative real-time PCR (qRT-PCR) was conducted with SYBR^® Premix Ex Taq™ (TaKaRa) on an ABI ViiA™ 7 Real-Time PCR System (Applied Biosystems, Indianapolis, IN, USA) (40). The primers for HMGA1 were specifically designed as follows: Forward primer: 5′-AGTGAGTCGAGCTCGAAGTC-3′; reverse primer: 5′-GTCTCTTAGGTGTTGGCACT-3′.

The Cell Counting Kit-8 (CCK-8; Dojindo Laboratories, Kumamoto, Japan) was employed to assess cell viability. The cells were inoculated into 96-well plates at a density of 1,000 cells per well and incubated overnight. Subsequently, 100 µL of CCK-8 working solution was added to each well and incubated for 1 hour at 37°C. The optical density at 450 nm was measured daily for four consecutive days utilizing a microplate reader. Cell growth curves were constructed by graphing OD₄₅₀ values over time to assess cell viability and proliferation.

Cells in the logarithmic growth phase were collected, resuspended, and diluted to the specified concentration. One thousand cells per well were inoculated into 6-well plates and cultivated for ten days under conventional conditions, with periodic inspection. Upon the formation of visible colonies, the culture media was discarded, and the cells were meticulously washed twice with ice-cold phosphate-buffered saline (PBS). Cells were subsequently fixed with 4% paraformaldehyde for 20 minutes at ambient temperature, followed by staining with 0.1% crystal violet for 10 minutes. The quantity of colonies was assessed utilizing a gel imaging analysis equipment (G: BOX-F3EE, Syngene, Bangalore, India).

Cells in the logarithmic growth phase were collected, resuspended, and diluted accordingly. 1×10³ cells were inoculated per well in 6-well plates. The EdU incorporation test kit from RiboBio (Guangzhou, China) was utilized to evaluate cell proliferation in accordance with the manufacturer’s guidelines. Following staining, EdU-positive cells were observed utilizing a fluorescence microscope. Quantification involved counting EdU-positive cells and total cells in a minimum of six randomly chosen fields per sample.

In 6-well plates, cells were placed at a density of 2 × 10⁵ cells per well and incubated overnight to ensure they adhere. A linear incision was made in the confluent monolayer utilizing a 10 µL pipette tip positioned perpendicularly to the plate surface. The wells were washed thrice with PBS to eliminate unattached cells, after which fresh serum-free media was introduced. The plates were thereafter incubated at 37°C in a humidified incubator containing 5% CO₂. Wound closure was assessed by taking images at 0 and 48 hours utilizing bright-field microscopy.

A Transwell test was performed to evaluate cell migration and invasion. Cells were inoculated into 24-well Transwell insert chambers (BD Biosciences, USA) filled with serum-free media. For the invasion assay, the inserts were pre-treated with 2% Matrigel. The lower chamber housed a medium supplemented with 20% FBS to serve as a chemoattractant. After a 48-hour incubation period, cells adhering to the upper membrane surface were meticulously removed, while those that penetrated the membrane were fixed with paraformaldehyde and stained with crystal violet. The dyed cells were subsequently photographed and quantified using a microscope.

Using R software (v4.3.3) and Python software (v3.9.19), statistical analysis was carried out. Wilcoxon’s test and Pearson’s correlation coefficient were used to evaluate the significance of differences between groups. The subsequent interpretation was utilized for P-values: *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. Insignificant differences were designated as ‘ns’. The statistical techniques and significance levels used evaluated the validity of the experimental findings and offered strong backing for the conclusions.

The overarching concept and substance of our investigation are illustrated in Figure 1 . To begin with, We picked 13 colorectal cancer tumor tissue samples representing various tumor stages from the dataset. Following rigorous quality control and batch effect elimination, a total of 37,236 high-quality cells were preserved. Following dimensionality reduction and clustering investigations, ten different cell groupings were discovered. Clusters were annotated according to the differential expression of established cell-type-specific marker genes as follows: B-plasma cells, T-NK cells, endothelial progenitor cells (EPCs), macrophages, fibroblasts, proliferating cells, mast cells (MCs), endothelial cells (ECs), smooth muscle cells (SMCs), and Schwann cells ( Figure 2A ). Figure 2B displays the five principal marker genes for each cluster, with EPCs identified by PHGR1, PIGR, TFF3, FABP1, and MUC2. EPCs were primarily concentrated in tumor stages III, I, and IV ( Figure 2C ). InferCNV analysis indicated a unique CNV pattern in EPC, suggesting that it may be related to tumor origin ( Supplementary Figure 1A ). Moreover, EPCs and proliferating cells demonstrated markedly increased G2/M and S scores. EPCs exhibited the greatest nFeature_RNA counts, succeeded by macrophages, proliferating cells, ECs, fibroblasts, and SMCs. In terms of cellular stemness, macrophages, MCs, ECs, and SMCs had comparatively high scores, while EPCs ranked marginally lower although remained elevated relative to other cell types ( Figure 2D ). Figure 2E depicts the significant relative density of EPCs in the G2/M score, S score, and nFeature_RNA. EPCs, macrophages, MCs, and T-NK cells demonstrated elevated concentrations of cell stemness AUC. A bubble plot demonstrated that MYC and KLF4 were among the most highly expressed stemness-associated genes in EPCs ( Figure 2F ). UMAP imaging validated the heightened expression of KLF4 and MYC primarily within the EPC population ( Figure 2G ). To enhance the characterization of the TME and EPCs in CRC, we performed a series of enrichment analyses. The heatmap illustrated differential gene expression among all cell types identified genes like KRT8, LGALS4, PIGR, PHGR1, FCGBP, SPINK4, and TFF3 as considerably enriched in EPCs ( Figure 2H ). PINK1 was identified as a significantly enriched gene ( Figure 2I ). The five most upregulated genes in EPCs were PDSS1, PODXL, CXCL3, HSPG2, and CEBPD, while the five most downregulated genes were MEI1, FCRL5, LYZ, IGFBP4, and GLUL ( Figure 2J ). The GO Biological Process (GO-BP) analysis indicates that EPCs are mostly linked to aerobic respiration, oxidative phosphorylation, cellular respiration, ribose phosphate metabolism, ribonucleotide metabolism, and purine ribonucleotide metabolism ( Figure 2K ). GO-GSEA analysis additionally associated EPCs with epithelial growth, epithelial cell differentiation, inorganic anion transport, digestion, digestive system functions, cell–cell junction architecture, and organ morphogenesis in animals ( Figure 2L ). Finally, the enrichment network ( Figure 2M ) identified pathways associated with cell adhesion, cardiac muscle contraction, mechanoreceptor differentiation, tissue homeostasis, digestion, wound healing regulation, intermediate filaments, negative regulation of endopeptidase activity, external encapsulating structures, bacterial defense response, and serine hydrolase activity in relation to EPCs.

We classified a total of 7,759 tumor cells into separate subpopulations based on their DEGs and designated each cluster according to its most prominently expressed marker gene: C0 NAP1L1⁺ TCs, C1 CA2⁺ TCs, C2 MKI67⁺ TCs, and C3 ITLN1⁺ TCs ( Figure 3A ). The UMAP visualization revealed that the hallmark genes NAP1L1, CA2, MKI67, and ITLN1 were primarily expressed in their specific subpopulations ( Figure 3B ). Moreover, these four marker genes had the greatest expression density within their respective clusters ( Figure 3C ). C0 NAP1L1⁺ TCs comprised the predominant fraction and were significantly represented across all five tumor stages. C1 CA2⁺ TCs exhibited greater abundance in tumor stages 0, I, and III. C2 MKI67⁺ TCs exhibited increased proportions predominantly in stages II and III, ranking second only to C0 NAP1L1⁺ TCs in stage IV. C3 ITLN1⁺ TCs were less prevalent overall but exhibited a comparatively greater presence in stage III ( Figure 3D ). We subsequently analyzed the distribution of G2/M score, S score, nCount_RNA, nFeature_RNA, CNV score, and cell stemness AUC among the tumor cell subpopulations with UMAP plots ( Figure 3E ). We additionally analyzed the disparities in expression density of the G2/M score, S score, and cell stemness AUC among these clusters ( Figure 3F ). C2 MKI67⁺ TCs demonstrated the greatest G2/M and S scores, but C0 NAP1L1⁺ TCs and C2 MKI67⁺ TCs displayed comparatively increased nFeature_RNA and nCount_RNA values. Nonetheless, the variations in cell stemness AUC among the subpopulations were not statistically significant. Violin plots distinctly demonstrated that the G2/M and S scores of C2 MKI67⁺ TCs were significantly elevated compared to those of other clusters, although the difference in cell stemness AUC among subpopulations was minimal ( Figure 3G ). Despite the absence of significant variations in overall cell stemness among clusters, we continued to investigate the expression of genes associated with stemness. Figure 3H ’s bubble plot demonstrated that C0 NAP1L1⁺ TCs exhibited elevated levels of KDM5B; C1 CA2⁺ TCs expressed EPAS1 and KLF4; C2 MKI67⁺ TCs displayed significant expression of EZH2, NOTCH1, MYC, and CD44; whereas C3 ITLN1⁺ TCs mostly expressed KLF4.

We examined the differentiation capability of tumor cell subpopulations with CytoTRACE and CytoTRACE2. The UMAP visualization of CytoTRACE scores indicated that C2 MKI67⁺ TCs and C0 NAP1L1⁺ TCs exhibited comparatively elevated scores ( Figure 4A ). Violin plots further demonstrated that C2 MKI67⁺ TCs exhibited the highest CytoTRACE score, succeeded by C0 NAP1L1⁺ TCs ( Figure 4B ). Moreover, CytoTRACE scores were heightened in cells during the G2/M and S phases, reaching their zenith at tumor stage II, subsequently followed by stages IV, I, and III. Tumor stage 0 displayed the lowest CytoTRACE scores ( Figures 4C, D ). Correspondingly, UMAP visualization of CytoTRACE2 scores ( Figure 4E ) and violin plot analyses ( Figure 4F ) revealed that C2 MKI67⁺ TCs exhibited the highest CytoTRACE2 score, succeeded by C0 NAP1L1⁺ TCs. In accordance with CytoTRACE, CytoTRACE2 scores were elevated during the G2/M and S phases ( Figure 4G ). In contrast to CytoTRACE, CytoTRACE2 scores were elevated in tumor stages II and I, subsequently followed by stage 0 ( Figure 4H ). Moreover, the CytoTRACE2-Potency analysis characterized all tumor cell subpopulations, indicating that the majority of C1 CA2⁺ TCs and C3 ITLN1⁺ TCs were differentiated. Conversely, C2 MKI67⁺ TCs and C0 NAP1L1⁺ TCs—particularly C2 MKI67⁺ TCs—encompassed cells at several differentiation stages, including differentiated, unipotent, oligopotent, and multipotent ( Figure 4I ). CytoTRACE findings together reveal that C2 MKI67⁺ TCs possess the lowest differentiation and the highest differentiation potential, succeeded by C0 NAP1L1⁺ TCs, whereas C1 CA2⁺ TCs demonstrate the most differentiation and the lowest potential ( Figure 4J ). To further examine developmental disparities among tumor subpopulations, we utilized Monocle and Slingshot for trajectory inference. Monocle-based pseudotime ordering on UMAP elucidated the relative pseudotime positions of each subpopulation ( Figure 4K ). The established pseudotime trajectory progressed from right (early) to left (late) ( Figure 4L ). Mapping tumor subpopulations onto this trajectory revealed that C0 NAP1L1⁺ TCs and C2 MKI67⁺ TCs generally occupied the early pseudotime segment, C3 ITLN1⁺ TCs were situated in the intermediate region, and C1 CA2⁺ TCs were primarily found in the late segment ( Figure 4M ). At branching point 1, located at the commencement of the trajectory, we categorized the trajectory into three states: state 1–3 ( Figure 4N ). State1 and state3 denoted the initial pseudotime, whereas state2 encompassed the latter portion of the early segment, the midpoint, and the late pseudotime. State 1 was characterized by a predominance of C0 NAP1L1⁺ TCs and C2 MKI67⁺ TCs; state 2 exhibited an enrichment of C1 CA2⁺ TCs, succeeded by C0 NAP1L1⁺ TCs and C3 ITLN1⁺ TCs; state 3 was mostly constituted of C0 NAP1L1⁺ TCs ( Figure 4O ). Subsequently, we utilized Slingshot to reconstruct a lineage trajectory, designated Lineage1, with the inferred sequence: C1 CA2⁺ TCs → C3 ITLN1⁺ TCs → C0 NAP1L1⁺ TCs → C2 MKI67⁺ TCs ( Figure 4P ). The expression trends of hallmark genes along Lineage1 indicated that NAP1L1 and MKI67 elevated towards the trajectory’s conclusion, CA2 was predominant at the onset, and ITLN1 reached its zenith in the intermediate segment ( Figure 4Q ). Differential expression study across Lineage1 stages indicated that at the terminal stage characterized by C2 MKI67⁺ TCs, genes including CENPF, UBE2C, ASPM, MKI67, and PTTG1 exhibited elevated expression levels. The GO-BP enrichment analysis of these genes underscored their participation in processes such as chromosome segregation, mitotic division, cell cycle regulation, and chromatid organization ( Figure 4R ).

We employed volcano plots to ascertain the five most upregulated and five most downregulated genes for each tumor cell subpopulation ( Figure 5A ). In C0 NAP1L1⁺ TCs, the five most upregulated genes were RPL12, RPS6, RPS7, RPL7, and RPL6, whereas the five most downregulated genes were FCGBP, PIGR, GUCA2A, ITM2C, and PHGR1. In C1 CA2⁺ TCs, the five most upregulated genes were S100A6, RHOC, TMEM54, GPA33, and GUCA2A, while the five most downregulated genes were NUCKS1, YBX1, RPS8, RPL11, and PTMA. In C2 MKI67⁺ TCs, the five most upregulated genes were NCAPH, TTK, NCAPG, ASF1B, and ARHGAP11A, whereas the five most downregulated genes were TSPAN1, FXYD3, NEAT1, LGALS4, and S100A6. In C3 ITLN1⁺ TCs, the five most upregulated genes were NEURL1, RAP1GAP, REG4, CLCA1, and ITLN1, while the five most downregulated genes were CD24, SFN, LDHA, PRSS3, and GPRC5A. Word cloud analysis indicated that C0 NAP1L1⁺ TCs were significantly associated with terms including localization, mRNA, telomere, and telomerase; C1 CA2⁺ TCs correlated with terms such as ion, transport, lipid, stress-activated, and metabolic; C2 MKI67⁺ TCs were connected to cycle, mitotic, and transition-related terms; and C3 ITLN1⁺ TCs were associated with unfolded protein response, glycosylation, and stress ( Figure 5B ). GO enrichment analysis of DEGs within each subpopulation indicated that C2 MKI67⁺ TCs were predominantly linked to ribonucleoprotein complex formation, rRNA metabolic processes, RNA localization, rRNA processing, and ribonucleoprotein complex assembly ( Figure 5C ). The enrichment network for C2 MKI67⁺ TCs emphasized pathways including the regulation of the metaphase/anaphase transition in the cell cycle, the positive regulation of protein localization to chromosomes and telomeric regions, and the regulation of the adaptive immune response ( Figure 5D ). GO-GSEA analysis revealed that C2 MKI67⁺ TCs were significantly linked to chromosome organization, mitotic cell cycle processes, chromosome segregation, nuclear division, organelle fission, and ribonucleoprotein complex biogenesis ( Figure 5E ). Two sample enriched phrases, nuclear division and nuclear chromosome segregation, were illustrated in detail ( Figure 5F ).

A heatmap illustrated the DEGs across subpopulations, revealing that C2 MKI67⁺ TCs exhibited elevated expression levels of CENPF, TPX2, TOP2A, ATAD5, BRCA2, MKI67, PTTG1, TUBA1B, HMGB2, and HMGB1. GO-BP analysis of these genes demonstrated significant enrichment in nuclear division, DNA recombination, somatic diversification of immune receptors through germline recombination at a single locus, somatic cell DNA recombination, and somatic diversification of immune receptors ( Figure 5G ).

We depicted intercellular communication across all cell types with CellChat. A chord diagram illustrated the total intensity and quantity of interactions among all cells ( Figure 6A ). In the chord diagram illustrating the signaling from the important subpopulation C2 MKI67⁺ TCs to other cells, the interaction intensity from C2 MKI67⁺ TCs to macrophages and ECs was significantly greater. Concurrently, macrophages acquired enhanced signaling from fibroblasts, ECs, SMCs, and Schwann cells ( Figure 6B ). The heatmap illustrating receptor-ligand expression patterns across all cells is presented for both outgoing and incoming signaling patterns ( Figure 6C ). MIF was significantly expressed in C2 MKI67⁺ TCs within the outgoing signaling pathways. Consequently, we concentrated on the MIF-associated signaling pathway for enhanced visualization. Within the MIF signaling network, C2 MKI67⁺ TCs demonstrated a significant sender importance score, while macrophages displayed a notably high receiver important score ( Figure 6D ). Subsequent comprehensive investigation indicated that MIF was significantly expressed in C2 MKI67⁺ TCs, CD74 was mostly found in macrophages, and CD44 was detected in both MCs and macrophages ( Figure 6E ). Through chord diagrams and circular plots, we illustrated that the communication strength along the MIF-(CD74+CD44) axis was much greater between C2 MKI67⁺ TCs and macrophages ( Figures 6F, G ). A hierarchy diagram ultimately validated that across both the MIF signaling network and the MIF-(CD74+CD44) axis, C2 MKI67⁺ TCs exert a significant paracrine influence on macrophages ( Figure 6H ).

We depicted the five principal transcription factors (TFs) across all tumor subpopulations utilizing a heatmap ( Figure 7A ). The five most prominent transcription factors in the C2 MKI67⁺ T cell subpopulation were E2F8, E2F1, E2F7, HMGA1, and BPTF. Furthermore, we evaluated all TFs in C2 MKI67⁺ TCs according to their specificity scores, producing results that align with those presented in the heatmap ( Figure 7B ). Subsequently, we depicted the expression of these five predominant TFs on a UMAP plot, revealing elevated expression levels solely within the C2 MKI67⁺ T cell subpopulation ( Figure 7C ). The relative expression levels of these five TFs were elevated in C2 MKI67⁺ tumor cells ( Figures 7D-H ). Additionally, violin plots were employed to effectively compare the expression levels of the five TFs across various tumor cell subpopulations and cell cycle phases. All five TFs exhibited peak expression in C2 MKI67⁺ tumor cells. HMGA1 expression was most elevated in C2 MKI67⁺ TCs, followed by a comparatively high expression in C0 NAP1L1⁺ TCs. E2F8, HMGA1, and BPTF exhibited increased expression in cells during the G2/M phase; E2F1 was primarily expressed in the S phase; and E2F7 demonstrated no substantial difference between the G2/M and S phases, although it was marginally elevated in the S phase. Subsequently, cells in the G1, G2/M, and S phases were distinctly seen on UMAP plots ( Figure 7I ), corroborating prior results that cells in the G2/M phase predominantly clustered within the C2 MKI67⁺ TCs group. Concentrating on G2/M phase cells, we ranked all TFs by their specificity scores, identifying E2F8, E2F1, HMGA1, E2F7, and SOX9 as the five most prominent. HMGA1 was distinctly represented on a UMAP plot, exhibiting significant expression in C2 MKI67⁺ TCs.

Utilizing the CSI matrix, we categorized the TFs of all tumor cell subpopulations into three regulatory modules, labeled M1, M2, and M3 ( Figure 8A ). Subsequently, we illustrated the distribution of M1, M2, and M3 throughout all subpopulations utilizing UMAP plots ( Figure 8B ). M1 was mostly concentrated in C1 CA2⁺ TCs; M2 was chiefly dispersed throughout C0 NAP1L1⁺ TCs, C1 CA2⁺ TCs, and C3 ITLN1⁺ TCs; whereas M3 was predominantly linked to C2 MKI67⁺ TCs. Violin plots effectively illustrated the expression levels of these regulatory modules among tumor subpopulations, with C2 MKI67⁺ TCs exhibiting the highest expression in M3, succeeded by C0 NAP1L1⁺ TCs ( Figure 8C ). Furthermore, we executed pySCENIC analysis and implemented a novel dimensionality reduction and clustering based on the regulatory activity of TFs in CRC tumor cells, culminating in a new UMAP visualization depicted in Figure 8D . Facet plots illustrated the distributions of tumor subpopulations and cells in the G1, G2/M, and S phases. M1, M2, and M3 were subsequently displayed on the new UMAP plot ( Figure 8E ). The findings revealed that C1 CA2⁺ TCs constituted a greater percentage in M1; C0 NAP1L1⁺ TCs and C1 CA2⁺ TCs were more prevalent in M2; and C2 MKI67⁺ TCs and C0 NAP1L1⁺ TCs were predominant in M3. These results aligned with the earlier visual representations derived from the original UMAP graphic. Subpopulations of tumors within M1, M2, and M3 were subsequently graded based on their regulon activity scores ( Figure 8F ). In M1, C1 CA2⁺ TCs were ranked highest, whereas C2 MKI67⁺ TCs were ranked lowest. In M2, C3 ITLN1⁺ TCs achieved the highest ranking, while C2 MKI67⁺ TCs rated the lowest once more. In M3, C2 MKI67⁺ TCs were the most prominent. Ultimately, we prioritized TFs inside M3 according to the proportion of variance between subtypes, highlighting E2F8, E2F1, and HMGA1 as the foremost candidates ( Figure 8G ).

To investigate the possible connections between CRC tumor cell subpopulations and cardiac dysfunction, which may result in heart failure, we assembled a gene collection associated with heart failure and myocardial fibrosis. This gene set was compiled from current research and established consensus, encompassing genes expressed in colorectal tumor tissues that are recognized to be linked to heart failure and myocardial fibrosis. Utilizing this gene set, we evaluated each tumor cell subpopulation for its association with heart failure and myocardial fibrosis. Initially, we categorized tumor subpopulations based on cell cycle stages ( Figure 9A ). Heart failure scores were elevated in C0 NAP1L1⁺ TCs (G1), C0 NAP1L1⁺ TCs (G2/M), C0 NAP1L1⁺ TCs (S), C2 MKI67⁺ TCs (G2/M), and C2 MKI67⁺ TCs (S). Concurrently, cardiac fibrosis scores were heightened in C0 NAP1L1⁺ TCs (G2/M), C0 NAP1L1⁺ TCs (S), C1 CA2⁺ TCs (G2/M), and C1 CA2⁺ TCs (S). Subsequently, we evaluated tumor subpopulations across various tumor stages ( Figure 9B ). The cardiac failure scores were significantly elevated in C0 NAP1L1⁺ TCs (stage II), C1 CA2⁺ TCs (stage II), and C2 MKI67⁺ TCs (stage II). Myocardial fibrosis scores were heightened in C0 NAP1L1⁺ TCs (stage I), C1 CA2⁺ TCs (stage 0), C1 CA2⁺ TCs (stage I), and C0 NAP1L1⁺ TCs (stage II). We additionally illustrated the expression pattern and density of heart failure scores among subpopulations by UMAP plots ( Figures 9C–E ), indicating that peak values are primarily concentrated in C0 NAP1L1⁺ TCs and C2 MKI67⁺ TCs. Violin plots comparing heart failure ratings throughout tumor stages indicated that stage II tumors had the highest scores, succeeded by stage IV ( Figure 9F ). Violin plots illustrating inter-subpopulation disparities indicated that C0 NAP1L1⁺ TCs and C2 MKI67⁺ TCs possessed considerably elevated heart failure scores, whereas C1 CA2⁺ TCs demonstrated moderate values, and C3 ITLN1⁺ TCs displayed the lowest scores ( Figure 9G ).

We conducted pertinent in vitro research centered on the HMGA1 gene to investigate its function as a principal transcription factor in C2 MKI67⁺ tumor cells. HMGA1 knockdown was performed on two CRC cell lines, HCT116 and HT-29. After knockdown, the mRNA and protein expression levels of HMGA1 were markedly diminished in both cell lines ( Figure 10A ). Moreover, in comparison to control groups, cell viability was significantly reduced in both cell lines following HMGA1 knockdown ( Figure 10B ). Colony formation experiments revealed a marked decrease in colony numbers following HMGA1 silencing ( Figure 10C ). The EDU labeling experiments demonstrated that decreased HMGA1 levels impeded cell growth ( Figure 10D ). Furthermore, wound healing studies demonstrated that HMGA1 knockdown not only inhibited cell migration ( Figure 10E ) but also markedly reduced the wound closure rate and cellular proliferation capability ( Figure 10F ). Transwell experiments demonstrated that the migratory and invasive capacities of HCT116 and HT-29 cells were significantly diminished after HMGA1 knockdown ( Figures 10G, H ).

To study the potential clinical relevance of the MKI67⁺/HMGA1 regulatory network, we created a risk prediction model based on the major tumor cell subpopulation C2 MKI67⁺ TCs. First, univariate Cox regression analysis was done to identify significantly differentially expressed genes and evaluate their prognostic significance ( Figure 11A ). Subsequently, LASSO regression analysis was employed to identify genes significantly associated with prognosis ( Figure 11B ). Multivariate Cox regression analysis found MAFK, IRF7, and HEYL as independent negative prognostic variables (hazard ratio [HR] > 1) ( Figure 11C ). Subsequently, we computed the coefficient values of the chosen genes to assess their correlation with survival outcomes ( Figure 11D ). Utilizing the relevant formula, we calculated an MKI67⁺ TCs risk score (MTRS) for each patient and conducted a differential gene expression study. Patients in the TCGA cohort were categorized into two groups, High MTRS and Low MTRS, utilizing the best cutoff point for MTRS. Survival analysis through scatter plots indicated that the High MTRS group demonstrated markedly reduced survival rates and a worse prognosis ( Figure 11E ). A heatmap depicted the differential expression of five model-related genes across the two groups, revealing increased expression of HEYL, MAFK, and IRF7 in the High MTRS group, and heightened expression of ETS2 and XBP1 in the Low MTRS group, corroborating prior findings ( Figure 11F ). The Kaplan-Meier survival curves showed that the High MTRS group had significantly worse survival outcomes compared to the Low MTRS group ( Figure 11G ). ROC curves demonstrated the model’s predictive accuracy and consistency, exhibiting AUC values at 1-, 3-, and 5-year intervals within the TCGA cohort ( Figure 11H ). Furthermore, risk ratings exhibited a negative correlation with OS, suggesting that elevated risk scores forecasted diminished OS ( Figure 11I ). The model’s accuracy was further validated by calculating the C-index, which resulted in values exceeding 0.5 at 1, 3, and 5 years, confirming its strong predictive ability ( Figure 11J ). Correlation heatmaps and scatter plots were created to thoroughly analyze the correlations among the five model genes, risk scores, and OS ( Figure 11K ). Additionally, subpopulation analyses based on race, age, tumor stages (I, II, III, IV), and risk scores were performed to assess their varying effects on OS ( Figures 11L, M ). In the end, Figure 11 F’s heatmap provided a detailed visualization of the expression differences of the five model genes between the two risk groups ( Figure 11N ). Kaplan-Meier survival analyses revealed significant disparities in survival between groups exhibiting high and low expression of these genes ( Figure 11O ). Elevated expression of ETS2 and XBP1 connected with improved survival, aligning with their protective functions, while increased expression of MAFK, HEYL, and IRF7 was linked to lower survival outcomes, reinforcing their status as negative prognostic indicators.

We first assessed the amounts of immune cell infiltration in the High MTRS and Low MTRS groups ( Figures 12A, B ). The immune cell types with notable abundance comprised M0 macrophages, CD4⁺ naïve T cells, M2 macrophages, and activated CD4⁺ memory T cells. The High MTRS group demonstrated a markedly elevated Tumor Immune Dysfunction and Exclusion (TIDE) score ( Figure 12C ), indicating enhanced immune suppression and an increased probability of unfavorable clinical outcomes in this cohort. Moreover, the Low MTRS group had a greater proportion of cell subpopulations lacking the immunological checkpoint markers CTLA-4 and PD-1, in contrast to their diminished presence in the High MTRS group ( Figure 12D ). The reduced expression of CTLA-4 and PD-1, pivotal immune checkpoint proteins linked to T cell exhaustion and immunological suppression, in the Low MTRS group may indicate a more active or less repressed immune condition. Consequently, in comparison to the High MTRS group, the Low MTRS group likely possesses a more vigorous antitumor immune response and a TME characterized by diminished immune suppression.

Among the five model genes, XBP1, ETS2, and HEYL exhibited increased CNV gains ( Figure 12E ). The High MTRS group consistently demonstrated increased ESTIMATE, Immune, and Stromal scores ( Figure 12F ). Correlation analyses demonstrated positive relationships between the risk prediction model and M0, M1, and M2 macrophage phenotypes, with notably robust positive correlations between HEYL and both M2 and M0 macrophages ( Figures 12G-I ). Additionally, the relationships among the five model genes, risk score, OS, and immune checkpoint-related genes were examined ( Figure 12J ). HEYL and IRF7 exhibited a positive correlation with the majority of immunological checkpoint genes, whereas ETS2 shown a negative correlation with most of these genes. Significantly, drug sensitivity analyses demonstrated that the High MTRS group displayed reduced half-maximal IC50 values for A.770041, AG.014699, WH.4.023, TW.37, Sunitinib, and GSK269962A, in comparison to the Low MTRS group, signifying heightened sensitivity to these compounds ( Figure 12K ). In contrast, Sorafenib and Mitomycin C exhibited reduced IC50 values in the Low MTRS group, indicating these medications may be more efficacious for patients within this subpopulation.

We initially showed the DEGs between the high- and low-risk groups using a volcano plot, highlighting both upregulated and downregulated genes ( Figure 13A ). Subsequently, many enrichment analysis methodologies were employed, including Gene Ontology Cellular Component (GO-CC), GO-BP, Molecular Function (GO-MF), Kyoto Encyclopedia of Genes and Genomes (KEGG), and GSEA. In the GO-CC analysis, the upregulated genes were predominantly associated with the collagen-rich extracellular matrix, whereas the downregulated genes were primarily connected to the glutamatergic synapse and postsynaptic membrane ( Figure 13B ). In GO-BP, upregulated genes exhibited considerable enrichment in response to estradiol, whereas downregulated genes were concentrated in lipid catabolic and carbohydrate biosynthetic activities ( Figure 13C ). The GO-MF results demonstrated that elevated genes were primarily enriched in extracellular matrix structural components and heparin binding ( Figure 13D ). KEGG pathway analysis indicated that genes exhibiting elevated expression were predominantly associated with the PPAR signaling pathway, whereas genes demonstrating reduced expression were connected to pathways related to human papillomavirus infection ( Figure 13E ). GSEA corroborated these results, indicating that upregulated genes were considerably enriched in collagen fibril organization, positive control of excitatory postsynaptic potential, postsynaptic specialization organization, and modulation of excitatory postsynaptic potential. In contrast, downregulated genes exhibited enrichment in uronic acid metabolic processes, cytochrome complex assembly, NADH dehydrogenase complex assembly, and mitochondrial RNA metabolic processes ( Figure 13F ). Finally, GSVA was applied to gene sets from both the High and Low MTRS groups, with the outcomes illustrated in Figure 13G .