The study design is presented in Figure 1.

We obtained RNA-seq data (TPM) and clinical annotations for TCGA-CRC from the Genomic Data Commons (https://portal.gdc.cancer.gov/) and converted all transcript abundances to TPM (transcripts per kilobase million) for downstream analyses, excluding genes with mean expression ≤ 0.1 TPM. The GEO microarray cohort GSE39582 served as an external validation set (23) (raw or normalized matrices available), with TCGA-CRC as the discovery cohort. We also curated CRC single-cell and spatial transcriptomic datasets from GEO: GSE200997 (16 tumors, 7 adjacent normals) (24)and GSE210038 (two tumors with paired adjacent normals) (25). Single-cell RNA-seq data were processed in Seurat v5.3.0 (26): expression matrices were imported with Read10X and assembled via CreateSeuratObject. Quality control retained cells with >200 and <3,000 detected genes, removed cells with mitochondrial content ≥20%, and excluded those with UMI counts (nCount_RNA) ≥20,000 (Supplementary Figure 1A). Batch effects were corrected with Harmony, followed by normalization and PCA (Supplementary Figure 1B). Using the top 20 PCs, we built a neighbor graph and clustered with FindNeighbors/FindClusters (resolution = 0.1, guided by clustree) (Supplementary Figure 1C), and visualized embeddings with RunUMAP. Cell types were annotated by canonical markers: B cells, endothelial cells, epithelial cells, fibroblasts, mast cells, monocytes, plasma B cells, and T cells. Differential expression was used FindMarkers with a two-sided Wilcoxon test (P < 0.05). NECSO scores were computed with AddModuleScore.For spatial transcriptomics, data were normalized with SCTransform, reduced and clustered with RunPCA/RunUMAP, and scored with AddModuleScore. Spatial distributions were visualized by SpatialFeaturePlot, defining spots with NECSO score > 0 as positive. Immunohistochemistry (IHC) images of colonic tissue were retrieved from the Human Protein Atlas (HPA).

Using TRPM4 as the index gene (9), we identified NECSO-related candidates in colorectal cancer (CRC) transcriptomes with a two-stage “differential expression + correlation” workflow. TPM-scale matrices were standardized (duplicate gene symbols averaged), and low-abundance genes (row-mean ≤ 0.1 TPM) were removed. Samples were stratified into tumor and normal groups by TCGA barcodes. Stage 1 (differential expression): for each gene, two-sided Wilcoxon rank-sum tests (tumor vs normal) were applied; genes with P < 0.05 formed the evaluation set. Stage 2 (tumor-only correlation): Pearson correlations between TRPM4 and each evaluation-set gene were computed, requiring |r| > 0.30, P < 0.001, and tumor expression variability sd > 0.1. Genes satisfying both screens were designated TRPM4-associated (NECSO-related) candidates and carried forward for pathway enrichment, tumor–immune microenvironment analyses, and construction of the prognostic model. All analyses were performed in R.

In the training cohort, we first ran univariable Cox proportional-hazards models to identify genes associated with prognosis (P < 0.05). These candidates were then entered into a LASSO–Cox penalized regression (R package glmnet (27)) to achieve variable shrinkage and derive a parsimonious signature. The risk score was computed as

where Expi is the expression level of gene i and βi its coefficient estimated by the LASSO model. Survival curves were generated with survminer (28), and time-dependent ROC analyses for 1-, 3-, and 5-year overall survival (OS) were performed using survROC (29) to assess discrimination. Model robustness was evaluated in the internal training set and an external independent cohort. To enhance interpretability, we conducted single-feature SHAP dependence analyses to quantify each gene’s marginal contribution to predicted risk, thereby strengthening the model’s potential clinical utility (30).

We used the R package clusterProfiler (31) for functional annotation and pathway analysis of NECSO-related candidates, covering KEGG and the three Gene Ontology (GO) domains—biological process (BP), molecular function (MF), and cellular component (CC). Multiple testing was controlled with the Benjamini–Hochberg procedure; terms with FDR-adjusted P < 0.05 were deemed significant. We then performed gene set enrichment analysis (GSEA) (32, 33) to compare pathway activity between high- and low-risk groups using MSigDB-curated gene sets. Enrichments were ranked by FDR-adjusted P (< 0.05) and normalized enrichment score (NES). Transcriptome-wide, gene-set–level activity was further assessed with GSVA (34) using MSigDB gene sets (35) (http://www.gsea-msigdb.org/gsea/index.jsp); limma was applied in R to compare GSVA scores and delineate functional differences across samples.

To visualize individualized prognostic predictions for CRC, we built a nomogram with the rms package incorporating sex, the model-derived risk score, age, and clinical stage. Each variable was assigned points proportional to its regression coefficient, and the total was mapped to predicted 1-, 3-, and 5-year overall survival (OS) probabilities. Clinical utility was assessed by decision curve analysis (DCA).

To assess immune infiltration in colorectal cancer, we applied six deconvolution tools—xCell (36), EPIC, MCP-counter (37), QUANTISEQ, CIBERSORT (38), and TIMER (39)—to the TCGA-CRC RNA-seq cohort. These methods infer per-sample relative abundance or enrichment of immune-cell subsets from predefined marker signatures. Group differences stratified by the NECSO score were tested with the Wilcoxon rank-sum test, and multi-algorithm outputs were integrated in heatmaps for visualization. We also summarized seven classes of immunomodulators (40, 41), the canonical seven immune subtypes, and the immunophenotype score (IPS) (42) to refine immune profiling within the NECSO framework. Finally, IPS and tumor immune dysfunction and exclusion (TIDE) scores were used to estimate potential benefit from immunotherapy and gauge immune-evasion propensity; IPS and TIDE were retrieved from TCIA (https://tcia.at/home/ and the TIDE server (http://tide.dfci.harvard.edu/), respectively (43).

We analyzed TCGA mutation annotation format (MAF) files with the R package maftools (44) and generated separate oncoplots (waterfall plots) for high- and low-risk groups. Using precomputed mutation frequencies, we displayed the top 15 recurrently mutated genes and overlaid clinical annotations indicating risk stratification.

To delineate crosstalk among epithelial subpopulations, we built cell–cell communication networks using CellChat (45), restricting analyses to “Secreted Signaling.” We estimated ligand–receptor probabilities with computeCommunProb, denoised interactions via filterCommunication, and aggregated signals to pathways using computeCommunProbPathway. Global network structure was summarized with aggregateNet, and pathway-specific patterns were visualized as heatmaps with netVisual_heatmap. To chart developmental trajectories, we applied Monocle2 (46): normalized matrices were imported, cluster DE genes (FDR-adjusted P < 0.05) served as ordering features, and default dimensionality reduction/cell ordering inferred pseudotime. Gene trends along pseudotime were z-score–standardized, smoothed, and visualized to highlight positively or negatively correlated sets. Lineage relationships during subtype bifurcation were inferred with Slingshot (47) using getLineages and getCurves to derive differentiation paths, which were projected onto UMAP for visualization. Cellular differentiation potential was quantified with CytoTRACE (48); the top 100 genes positively and negatively correlated with the CytoTRACE score were designated as immaturity and differentiation markers, respectively. CytoTRACE scores range from 0–1, where 0 denotes a more differentiated state and 1 a less differentiated state.

We used the R package oncoPredict (v0.2) (49)—leveraging drug-response profiles and transcriptomes from the Broad Institute’s CTRP and the Sanger GDSC—to predict each patient’s responsiveness to multiple clinically relevant agents. To prioritize therapeutics, we performed molecular dynamics (MD) simulations to estimate binding free energies and characterize interactions between candidate small molecules and their protein targets, thereby assessing targetability. Two-dimensional ligand structures were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/), and three-dimensional receptor structures from the RCSB Protein Data Bank (PDB) (https://www.rcsb.org/). Initial protein–ligand docking used the CB-Dock2 web platform, and resulting complexes were subjected to MD evaluation of stability and binding interactions.

Molecular dynamics (MD) simulations were performed with GROMACS 2025.2. Small-molecule ligands were prepared in AmberTools 22 with the GAFF force field; geometries were protonated and RESP charges computed in Gaussian 16W. Parameters were written into the topology. Simulations ran at 300 K and 1 bar using the amber99sb-ildn protein force field and TIP3P water. Na⁺ counterions neutralized the system. After steepest-descent minimization, systems underwent NVT and NPT equilibration for 100,000 steps each (τ = 0.1 ps; ~100 ps total). Production trajectories were generated for 50,000,000 steps with a 2 fs timestep (100 ns total). Post-processing with GROMACS tools yielded root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), solvent-accessible surface area (SASA), and radius of gyration (R_g). Binding thermodynamics were estimated by MM/GBSA, and conformational landscapes were profiled via free-energy surface analyses.

Cell lines NCM460, HCT-116, and SW480 were obtained from Haixing Biosciences (Suzhou, China). Cells were cultured in DMEM supplemented with 10% fetal bovine serum (FBS; Life Technologies, USA) and 100 U/mL penicillin–streptomycin at 37 °C in a humidified 5% CO₂ incubator. A lentiviral RNAi vector targeting NEK8 (and a negative control) was provided by GeneChem (Shanghai, China). For infection, HCT-116 and SW480 cells were seeded in 6-well plates at ~2 × 10⁵ cells/well and allowed to adhere overnight. When cultures reached ~30–50% confluence the next day, cells were exposed to lentivirus encoding shNEK8 or the negative control (NC). Stable polyclonal lines were generated by selection in puromycin (2 mg/mL) for 2 weeks. Knockdown efficiency was evaluated 72 h postinfection by Western blotting. shRNA target sequences were: shNEK8#1, 5′-CGGGTGATTGCTACACTTT-3′; shNEK8#2, 5′-TGGTGATCATCAAGCAGAT-3′.

Cells were lysed in RIPA buffer (Solarbio) with ultrasonication. Lysates were clarified by centrifugation at 14,000 × g for 15 min at 4 °C, and supernatants were collected. Protein concentration was measured by BCA. Equal protein amounts were resolved on 12.5% SDS–PAGE and transferred to PVDF membranes. Membranes were blocked in TBST with 5% skim milk and incubated with primary antibodies against NEK8 (Abcam, AB169406) and β-actin (ABclonal, AC026; 1:50,000), the latter as the loading control. After incubation with HRP-conjugated secondary antibody (ABclonal, AS014), signals were detected by enhanced chemiluminescence (ECL).

Total RNA was extracted from cultured cells with TRIzol (Invitrogen) per the manufacturer’s instructions. RNA concentration and purity were assessed by UV spectrophotometry; samples with an A260/A280 of 1.8–2.0 were accepted. Two micrograms of RNA were reverse-transcribed to cDNA using a commercial kit (TOYOBO, Osaka, Japan). Quantitative real-time PCR (qPCR) was performed with SYBR Green chemistry (Sevenbio, China). Primer sequences were: NEK8, forward 5′-CTTCGTGCAGATCCTGCTTG-3′ and reverse 5′-GGAGATGCCGAAATCACCGAT-3′; β-actin, forward 5′-CCATCGTCCACCGCAAAT-3′ and reverse 5′-GCTGTCACCTTCACCGTTCC-3′. Each 20 μL reaction contained 10 μL SYBR Green master mix, 1 μL forward primer (10 μM), 1 μL reverse primer (10 μM), 2 μL cDNA, and 6 μL nuclease-free water. β-actin served as the internal control. All samples were run in technical triplicate, and each experiment was independently repeated at least three times. Relative expression was calculated using the 2−ΔΔCt method.

To assess proliferation, CRC cells transduced with sh-NC or sh-NEK8 were seeded in 96-well plates. Growth was quantified with the Cell Counting Kit-8 (CCK-8; Solarbio, China) per the manufacturer’s protocol, and optical density (OD) was read at 450 nm. For colony formation, transduced cells were plated in 6-well plates at 1,000 cells/well and cultured for 10 days in RPMI-1640 with 10% fetal bovine serum. Colonies were fixed in methanol, stained with 1% crystal violet (Beyotime, China), imaged, and analyzed.

For the wound-healing assay, a linear scratch was made across confluent monolayers with a pipette tip, and cells were cultured in serum-free medium. Wound closure was imaged and quantified at 48 h. Cell migration was assessed in 24-well Transwell chambers with 8-µm polycarbonate membranes. Matrigel (100 μL; 1:8 in serum-free medium) was spread on the upper chamber and allowed to solidify. Cell suspensions were added to the upper chamber, and the lower chamber contained medium with 10% fetal bovine serum. At assay end, membranes were rinsed twice in PBS, fixed in methanol, and stained with 0.5% crystal violet. Cells remaining on the upper surface were gently removed with a cotton swab, and cells that migrated to the underside were counted in predefined fields. All incubations were at 37 °C in a humidified 5% CO₂ atmosphere.

HCT-116 and SW480 cells were fixed in ice-cold 70% ethanol for 12 h, then incubated with RNase at 37 °C for 30 min. Nuclei were stained with a propidium iodide (PI) kit (Sevenbio, China); PI staining was performed at 4 °C per the manufacturer’s instructions. Cell-cycle distribution was analyzed by flow cytometry within 1 h of staining. Data were processed and modeled with ModFit LT (Verity Software, Topsham, ME) and FlowJo v10 (Treestar, Ashland, OR).

For normally distributed variables, between-group differences were tested with the unpaired Student’s t-test; for non-normal data, the Wilcoxon rank-sum test was used. For multi-group comparisons, one-way ANOVA was applied to parametric data and the Kruskal–Wallis test to nonparametric data. Two-sided Fisher’s exact tests were used for contingency tables. All analyses were conducted in R (v4.5.0), GraphPad Prism (v9.5.0), and Python (v3.8).

We leveraged TCGA-CRC expression data to examine the putative role of TRPM4 in sodium-overload cell death (NECSO). TRPM4—a Na⁺-conducting member of the transient receptor potential melastatin subfamily—is broadly expressed in heart, brain, and immune cells, where it modulates membrane potential and ionic homeostasis. Aberrant TRPM4 activation has been linked to hypoxia, energy depletion, and sodium-overload–related cell death (9). Using Pearson correlation, we quantified associations between TRPM4 and all other genes, identifying 359 transcripts significantly correlated with TRPM4 (Supplementary Table 1), which we designated as NECSO-related genes. We then generated a heatmap and gene–gene interaction network for the top 50 correlates (Figure 2A). TRPM4 showed robust positive correlations with PP1R37, NADSY1, and AP1G2, suggesting co-involvement in TRPM4-mediated NECSO and a role in CRC initiation and progression, whereas a subset of genes was negatively correlated, consistent with inhibitory regulation of TRPM4-linked functions. KEGG analysis (Supplementary Table 2) revealed significant enrichment across pathways governing lipid metabolism, the p53 axis, and insulin signaling, indicating that sodium overload may promote cell death by perturbing lipid programs, modulating cell-cycle control, and altering metabolic balance (Figure 2B). Gene Ontology enrichment (Supplementary Table 3) likewise showed overrepresentation across BP, CC, and MF domains: fatty-acid and cellular lipid catabolic processes (BP) linked to lipid metabolism; dendritic spines and neuronal synapses (CC) indicated synaptic interfaces; and protein serine/threonine kinase activity and Hsp90 binding (MF) implicated signaling and stress responses (Figures 2C–E). To prioritize clinically relevant NECSO genes, we performed univariable Cox analysis followed by LASSO regression (Figures 2F–H), yielding a five-gene signature. A NECSO risk index derived from these genes (Figure 2I; Supplementary Table 4) was interpreted with SHAP modeling. Mean SHAP values (Figure 2J) indicated substantial contributions from ACOT11, NEK8, DRD4, CYTH2, and EPHB2; distribution plots (Figure 2K) showed that expression variability in these genes materially influenced the index. The SHAP waterfall plot (Figure 2L) showed positive contributions from ACOT11 (+0.264) and EPHB2 (+0.460) and negative contributions from NEK8 (−0.366), DRD4 (−0.352), and CYTH2 (−0.294). The final prediction was 1.23 relative to an expected baseline E[f(x)]=1.52, suggesting heterogeneous roles of these genes within the NECSO program that may inform prognosis.

Using coefficients for five NECSO-related genes (DRD4, NEK8, EPHB2, CYTH2, ACOT11), we built a prognostic risk-score model and dichotomized patients at the median into high- and low-risk groups. Individual risk was calculated as “ Risk score"=∑_i (β_i×Exp_i ), where Expi is the expression of gene i and βi its LASSO-derived coefficient. TCGA-CRC served as the training set and GSE39582 as the validation set. High-risk patients had significantly shorter survival than low-risk patients (Figure 3A). Concordant patterns of risk-score distributions and survival status in both cohorts supported model stability (Figure 3B). Heatmaps showed marked expression differences across strata for DRD4, NEK8, and EPHB2, among others (Figures 3C, D). We further evaluated performance by Kaplan–Meier analyses in GSE39582, TCGA, and TCGA-PFS: high risk consistently predicted worse outcomes (GSE39582, p = 0.01; TCGA, p < 0.0001) and shorter progression-free survival (TCGA-PFS, p < 0.0001) (Figure 3E). In univariable Cox models, age (HR = 1.030, p = 0.009), overall stage (HR = 2.469, p < 0.001), T stage (HR = 3.222, p < 0.001), M stage (HR = 5.088, p < 0.001), N stage (HR = 2.155, p < 0.001), and the risk score (HR = 520.572, p < 0.001) were significant, whereas sex was not (HR = 1.028, p = 0.911) (Figure 3F). In multivariable models, the risk score (HR = 112.259, p < 0.001) remained the sole independent prognostic factor; age also retained significance (HR = 1.041, p < 0.001), while stage and T/M/N categories were not significant (Figure 3G). Clinical characteristics differed by risk group: age and sex were comparable (p = 0.601 and 0.064), whereas overall stage, T stage, M stage, and N stage were skewed toward advanced disease in the high-risk group (Stage III/IV, T3/T4, M1, N1/N2) and earlier disease in the low-risk group (Stage I/II, T1/T2, M0, N0) (p < 0.001, 0.005, 0.030, < 0.001). A composite heatmap underscored the link between clinical staging and risk (Figure 3H). Stratified Kaplan–Meier curves showed poorer survival for high-risk patients irrespective of age (>65 vs ≤65 years, both p < 0.001), sex (male or female, both p < 0.001), or stage (early vs advanced, both p < 0.001) (Figure 3I). A nomogram integrating sex, risk score, age, and stage (Figure 3J) yielded individualized prognostic scores and survival probabilities; high-risk patients accrued higher points, with predicted probabilities of surviving beyond 5, 3, and 1 years of 0.842, 0.916, and 0.956. Calibration plots showed good agreement between predicted and observed survival at 1, 3, and 5 years (Figure 3K). Time-dependent ROC analyses produced AUCs of 0.651 (1 year), 0.698 (3 years), and 0.655 (5 years). Compared with other variables, the risk score achieved an AUC of 0.651, while age (AUC = 0.823) and overall stage (AUC = 0.746) showed higher discrimination; sex (AUC = 0.441), T stage (AUC = 0.644), M stage (AUC = 0.677), and N stage (AUC = 0.702) contributed less but provided complementary information (Figure 3L).

To assess differences in the tumor immune microenvironment associated with the NECSO-based risk score, we compared immune-cell infiltration and pathway activity between high- and low-risk groups. In the high-risk cohort, we observed reduced infiltration of key effector populations, including CD8⁺ T cells, dendritic cells, and regulatory T cells, which suggests an immunosuppressive microenvironment. This was further supported by the lower Immune/ESTIMATE scores in the high-risk group, which indicate a reduced presence of immune cells. While Immune/ESTIMATE scores were higher in the low-risk group, indicating a greater overall immune presence, this was associated with a more robust immune response, with higher infiltration of effector immune cells (Figure 4A). Pathway analyses further revealed greater activation of antigen-presenting cell (APC) co-stimulation and cytotoxic programs in the low-risk group, with attenuated responses in the high-risk group (Figure 4B). These results are consistent with the idea that low-risk tumors maintain a more active immune response, whereas high-risk tumors exhibit immune evasion. In line with this, immunomodulator profiling demonstrated higher expression of costimulatory molecules (e.g., CD28, CD80) and chemokines (e.g., CCL5, CXCL10) in the low-risk group, supporting the notion of a more favorable immune environment (Figure 4C).Correlative analyses between NECSO signature genes and immune infiltration highlighted DRD4 as strongly positively associated with T cells and dendritic cells, suggesting a role in shaping antitumor immunity (Figure 4D).Immune-subtype stratification further supported these findings, with low-risk tumors predominantly classified as C1 (Wound Healing), a subtype linked to immune activation and tissue repair, while high-risk tumors were enriched for the C2 (IFN-γ Dominant) subtype, which is associated with immune suppression and aggressive disease behavior (Figure 4E).Finally, the risk score correlated positively with the ESTIMATE score (R = 0.62, p = 0.018) and Immune score (R = 0.63, p = 0.025), indicating greater immune suppression in high-risk tumors. The correlation with the Stromal score (R = 0.53, p = 0.19) was weaker, further supporting the idea that immune evasion mechanisms play a key role in the poor prognosis of high-risk tumors (Figure 4F).

Using GSEA and GSVA, we profiled immune and biological pathway differences between high- and low-risk groups. In GSEA, the high-risk group was significantly enriched for tumor-aggressiveness pathways—axon guidance, ECM–receptor interaction, and focal adhesion—suggesting enhanced invasiveness and adaptive responses. By contrast, the low-risk group showed enrichment of metabolism-related programs (drug metabolism—other enzymes; pentose and glucuronate interconversions), indicative of stronger metabolic capacity and a more active immune milieu (Figures 5A, B). Consistently, GSVA revealed preferential enrichment of Hedgehog signaling, myogenesis, and Wnt/β-catenin signaling in the high-risk cohort, whereas the low-risk cohort exhibited higher activity in peroxisome and bile acid metabolism pathways, pointing to heightened metabolic and antioxidant responses (Figure 5C). Together, these results suggest that the low-risk group harbors more robust immune and metabolic activity, whereas the high-risk group aligns with invasive phenotypes and immune evasion. Tumor mutational burden (TMB) analyses supported this divergence. In the high-risk group, APC, TP53, and KRAS showed higher mutation frequencies with predominance of missense and frameshift variants (Figure 5D), whereas the low-risk group had fewer mutations, largely missense (Figure 5E). Boxplots confirmed significantly lower TMB in the low-risk group (p = 0.0016; Figure 5F), concordant with a positive correlation between risk score and TMB (R = 0.62, p = 2 × 10^-4; Figure 5G). Survival analyses indicated shorter overall survival for high-TMB patients (p = 0.031; Figure 5H), with the worst outcomes in the high-TMB/high-risk subgroup (Figure 5I), underscoring the complementary prognostic value of integrating TMB with the NECSO-based risk score.

Analyzing the GSE200997 single-cell RNA-sequencing (scRNA-seq) dataset, we profiled 49,589 cells from 16 colorectal cancer (CRC) specimens and 7 adjacent normal tissues (Figure 6A). A t-SNE embedding (Figure 6B) delineated the cellular landscape, revealing discrete clusters aligned with major lineages. Using canonical markers curated in CellMarker 2.0 (50), we annotated eight principal types: T cells, B cells, epithelial cells, plasma cells, fibroblasts, monocytes, endothelial cells, and mast cells (Figure 6E). Comparison of tumor versus normal composition (Figure 6C) showed pronounced shifts in immune and stromal compartments: proportions of B, T, and plasma cells changed, with notable increases in monocytes and fibroblasts in tumors, consistent with immune remodeling and stromal activation. Spatial distributions of hallmark genes (e.g., CD3D, S100A8, VWF, CD19, CD68, CD79A, COL1A1, EPCAM) matched the corresponding clusters on the t-SNE map, supporting robust cell-type annotation (Figure 6D). To quantify NECSO activity at single-cell resolution, we scored a NECSO gene set (DRD4, NEK8, EPHB2, CYTH2, ACOT11) using AddModuleScore. NECSO scores were elevated predominantly in epithelial cells (Figures 6F, G), implicating malignant epithelia as a principal compartment with heightened NECSO-related transcriptional programs.

Using scRNA-seq, we isolated 7,155 epithelial cells and partitioned them into four transcriptionally distinct subclusters (Figure 7A). Subcluster identities were defined by selective marker expression, and NECSO module scores mapped across subclusters revealed marked inter-subcluster heterogeneity (Figure 7B), implying divergent roles of these epithelial states within the tumor microenvironment. Composition analysis (Figure 7C) further showed altered epithelial ecology between normal and tumor tissues, with a notable expansion of the PCK1⁺ subcluster in tumors. To characterize function, we identified the top five Gene Ontology biological processes enriched among marker genes per cluster, yielding a molecular “fingerprint” for each (Figure 7D). The t-SNE embedding (Figure 7I) delineated a clear boundary between tumor-derived epithelia and their normal counterparts. CytoTRACE (Figure 7J) indicated reduced developmental potential in tumor epithelia, corroborated by boxplots (Figure 7K) showing lower scores—especially in GGH⁺Epi—consistent with a more constrained/differentiated state. We next interrogated epithelial crosstalk with CellChat: GGH⁺Epi displayed frequent interactions with C1_PCK1 and C2_GGH, but fewer with C3_ITLN1 and C4_HTR3E (Figure 7E). Signaling-pattern analysis (Figure 7F) highlighted prominent CypA, MK (midkine), and MIF activity in GGH⁺Epi. Network decomposition (Figure 7G) showed GGH⁺Epi acting chiefly as a sender in MK and as a receiver/mediator in MIF and CypA, underscoring its central role in immune modulation. To improve interpretability, we further prioritized the most significant ligand–receptor interactions that were biologically relevant to immune modulation. Specifically, the top interactions included the MIF axis [e.g., MIF–(CD74+CD44)] and the MDK axis [e.g., MDK–(ITGA6+ITGB1)], both showing high communication probabilities and consistent enrichment in the GGH⁺Epi-centered network (Figure 7H). Together, these findings nominate GGH⁺Epi as a pivotal regulator of the tumor immune microenvironment operating through multiple signaling axes and potentially contributing to immune evasion.

Having established the cellular provenance of NECSO signaling and the hub role of GGH⁺Epi, we next charted epithelial dynamics over pseudotime. Monocle2 revealed a continuous trajectory from early to late states, with distinct terminal clusters corresponding to specific epithelial subpopulations (Figure 8A). Stratifying by tissue of origin showed early stages dominated by normal epithelia and late stages enriched for tumor cells (Figure 8B), indicating that malignant epithelia occupy downstream states with phenotypic remodeling. Slingshot identified two principal differentiation branches leading to divergent terminal phenotypes (Figure 8C), consistent with bifurcating trajectories in CRC. To relate NECSO activity to crosstalk within the tumor immune microenvironment (TIME), epithelial cells were stratified by NECSO expression (NECSO_high vs NECSO_low) and ligand–receptor networks reconstructed with CellChat. Overall, NECSO_high epithelia showed more numerous and stronger interactions with fibroblasts, endothelial cells, and B/T cells (Figures 8D, E), whereas NECSO_low cells interacted preferentially with monocytes and mast cells. Significant ligand–receptor events mirrored these trends: NECSO_high epithelia acted chiefly as senders, while NECSO_low cells functioned as receivers (Figures 8F–H), suggesting that upregulated NECSO programs align with a driver phenotype. Pathway decomposition highlighted the IGFBP and MK (midkine) axes: along both, NECSO_high epithelia assumed the primary sender/influencer role, with fibroblasts, endothelial cells, and B cells as receivers/mediators (Figures 8I, J). Collectively, high NECSO expression is associated with a shift from signal reception to emission in epithelial cells, accompanied by directed remodeling of stromal and humoral compartments—implicating NECSO_high epithelia, via the IGFBP and MK axes, in active modulation of the TIME.

In the TCGA cohort, analysis of the five NECSO-signature genes (Figure 9A) showed NEK8, DRD4, EPHB2, and CYTH2 were broadly upregulated in tumors, whereas ACOT11 was significantly downregulated. Differences in both dispersion and median expression supported a clear tumor–normal separation at the NECSO-related molecular level. To resolve spatial features, we profiled two tumor sections of differing differentiation grades with matched adjacent normals by spatial transcriptomics (Figures 9B, C). NECSO scores were markedly elevated within tumor epithelium and increased in extent and magnitude with poorer differentiation; adjacent normals exhibited uniformly lower scores. Concordantly, NEK8/DRD4/EPHB2/CYTH2 showed patchy high expression in tumor epithelium and at the invasive front, colocalizing with NECSO-high regions, whereas ACOT11 was broadly reduced in tumors and comparatively higher in normal glands/epithelium. This pattern reproduced across both specimens, indicating tumor-restricted activation and regionalization of the NECSO program and further supporting the relatively protective profile of ACOT11.

To functionally validate the NECSO model genes, we prioritized NEK8 and performed genetic perturbation. A pan-cancer analysis with TIMER showed NEK8 upregulation in most tumors relative to adjacent normal tissue (multiple types p < 0.05; several p < 0.001), including colorectal (COAD/READ), gastric (STAD), esophageal (ESCA), lung adenocarcinoma and squamous carcinoma (LUAD/LUSC), breast (BRCA), head and neck squamous carcinoma (HNSC), liver (LIHC), and pancreatic (PAAD) cancers (Figure 10A). In TCGA-CRC, Kaplan–Meier analysis stratified by median NEK8 expression indicated worse overall survival in NEK8-high tumors (p = 0.04; Figure 10B). Human Protein Atlas IHC corroborated elevated NEK8 protein in CRC tissues (Figure 10C), and paired samples confirmed higher NEK8 mRNA in tumors (Figure 10D). In vitro, Western blotting showed higher NEK8 protein in HCT116 and SW480 than in normal NCM460 cells (Figure 10E), with significant quantification (Figure 10F). For functional assays, two independent shRNAs (sh-NEK8#1/#2) were introduced into HCT116 and SW480; both markedly reduced NEK8 protein versus negative control (NC) (Figure 10G; densitometry in Figure 10H) and achieved concordant transcript knockdown by qRT-PCR (Figure 10I).

Having shown that NEK8 is upregulated in CRC and associated with poor prognosis—and after establishing stable knockdown with two independent shRNAs in HCT116 and SW480—we performed phenotypic validation. NEK8 silencing markedly suppressed proliferation (CCK-8, 12–36 h; Figure 11A), reduced migratory and invasive capacities (wound-healing and Matrigel Transwell; Figures 11B, D), and decreased clonogenicity (Figure 11C). Flow cytometry further revealed an increased G1 fraction with concomitant decreases in S/G2, consistent with G1/S arrest (Figure 11E). These effects were concordant across both shRNAs and cell lines, indicating that NEK8 drives CRC cell proliferation, migration, and invasion and acts as a key promoter of malignant phenotypes.

Stratification by NEK8 expression using the immunophenoscore (IPS) revealed divergent immunotherapy-relevant patterns. Under a PD-1 monotherapy context (CTLA4^-/PD-1⁺), the NEK8-high group showed a significantly higher IPS (p = 0.041), suggesting greater putative responsiveness to PD-1 blockade; conversely, under combined checkpoint inhibition (CTLA4⁺/PD-1⁺), IPS was significantly lower in NEK8-high tumors (p = 0.0062) (Figure 12A). However, we acknowledge that these findings are speculative and should be interpreted as hypothesis-generating unless validated in an immunotherapy-treated CRC cohort. Further studies are necessary to verify the predictive value of NEK8 expression for immunotherapy response, including potential confounders like microsatellite instability (MSI) and other CRC-specific factors. In drug prioritization, we focused on vorinostat (SAHA), an FDA-approved oral histone deacetylase inhibitor for cutaneous T-cell lymphoma (51). Within CRC, NEK8-high tumors exhibited reduced in silico sensitivity to vorinostat (higher AUC/IC₅₀; p = 0.00088) (Figure 12B). Molecular docking indicated favorable NEK8–vorinostat binding (AutoDock Vina score = −7.9 kcal/mol) (Figure 12C). Subsequent 100-ns MD simulations supported rapid attainment and maintenance of complex stability, evidenced by plateaued RMSD (~0.6–0.8 nm), decreased then stabilized SASA, constant radius of gyration, persistence of 1–3 hydrogen bonds, and a single deep basin on the free-energy surface (Figures 12D–I). Collectively, these findings nominate NEK8 as a biomarker for immunotherapy stratification. Although NEK8-high samples show lower predicted cytotoxic sensitivity to vorinostat, computational evidence of stable binding supports mechanistic follow-up and structure-guided optimization.