Xiaoyu He
Wenhao Wang
Li Li3- 1Department of Gastroenterology, People’s Hospital of Deyang City, Deyang, Sichuan, China
- 2Department of Obstetrics and Gynecology, The Second Clinical Medical College, Shanxi Medical University, Taiyuan, Shanxi, China
- 3Department of Geriatrics, The Second Hospital of Shanxi Medical University, Taiyuan, Shanxi, China
- 4Key Laboratory of Cellular Physiology, Ministry of Education, Department of Physiology, Shanxi Medical University, Taiyuan, China
Background: Colorectal cancer (CRC) development and progression are linked to genetic factors, environmental influences, and dysregulated signaling pathways.
Methods: The differentially expressed pan-apoptotic genes (CPAN_DEGs) between CRC and normal colon tissues were screened from bulk RNA-sequencing (RNA-Seq) datasets. The putative biological functions of these CPAN_DEGs were explored through functional enrichment analysis and the protein-protein interaction (PI) network. Unsupervised clustering was used to stratify patients on the basis of CPAN_DEGs, and a prognostic model was constructed using LASSO dimensionality reduction. Based on the CPAN-index score, the patients were divided into the high-risk and low-risk groups, and the survival rates and immunophenotypes were compared. The predictive performance of the CPAN-index model was confirmed in an external validation set. The expression patterns of PANoptosis genes across different cell types in CRC samples, and the distribution of CPAN-index-positive cells within each subpopulation were analyzed using single cell RNA-Seq (scRNA-Seq) datasets. The expression of CDKN2A was confirmed in CRC cell lines, and its functional role was evaluated by gene knockdown.
Results: The expression levels of PANoptosis-related genes showed significant heterogeneity across CRC samples, and the highest percentage (87.4%) was that of apoptosis-related genes. The differentially expressed genes (DEGs) between the CRC and normal tissue samples were significantly enriched in pathways related to metabolism and immune regulation. The CPAN-index constructed using 11 CPAN_DEGs effectively distinguished CRC patients in to the high-risk and low-risk groups, and the high-risk group showed an “invasion-metabolism-immunosuppressive” phenotype, along with immune tolerance and non-classical immune escape. The CPAN-index gene CDKN2A was upregulated in the CRC cell lines, and knocking down the CDKN2A gene inhibited their proliferation and promoted apoptosis in vitro. ScRNA-Seq data revealed a higher proportion of CPAN-index-positive immune cells, and a lower proportion of tumor cells positive for CPAN-index, thus underscoring its critical role in the tumor immune microenvironment.
Conclusions: CDKN2A-mediated PANoptosis signaling network drives CRC progression by reshaping the immune microenvironment and metabolic reprogramming. The CPAN-index provides a new tool for accurate risk stratification of CRC patients, and suggests potential therapeutic strategies targeting the immunometabolism-death interaction network.
1 Introduction
Colorectal cancer (CRC) was the third most prevalent malignancy worldwide in 2020, with 1,931,590 newly diagnosed cases that accounted for 10% of all cancers. In addition, the mortality rate associated with CRC was 9.4%, with approximately 935,173 recorded deaths (1, 2). Given its high incidence and mortality, as well as the rise in incidence among individuals previously considered to be at low risk (3), CRC poses a significant economic burden. Smoking, drinking and weight gain are among the risk factors of CRC (4), whereas a healthy diet can reduce the risk associated with colon cancer (5) (6) (7). Mechanistically, CRC initiation and progression involve aberrant cellular proliferation and differentiation, along with increased resistance to apoptosis (8, 9). At the molecular level, pathogenesis of CRC is associated with the overactivation of phosphoinositol 3-kinase (PI3K)/AKT and transforming growth factor-β (TGF-β) signaling pathways due to multiple oncogenic mutations (10). The clinical manifestation of CRC, including fatigue, anemia, and weight loss have limited predictive power for CRC in older patients. As CRC screening programs become more widely available, many cases are detected prior to the appearance of symptoms (11). For symptomatic patients, the preferred diagnostic method is colonoscopy, although other endoscopic techniques may be selected. The primary radical treatment for patients with non-metastatic CRC is surgery (6, 12). In recent years, immunotherapy has gradually been incorporated for the treatment of advanced cases in combination with radiotherapy (13). Nevertheless, it is essential to elucidate the molecular mechanisms underlying the pathogenesis and development of CRC in order to devise more effective and personalized therapeutic strategies.
Programmed cell death, including pyroptosis, apoptosis, and necrosis, is the central mechanism of most anti-cancer therapies and the immune defense against tumor cells (11, 14, 15). As the three cell death modes overlap each other, Malireddi et al. proposed the concept of PANoptosis, which is stimulated by particular triggers and controlled by the PANoptosome complex (16, 17). Upstream receptors such as ZBP1, RIPK1, and AIM2 detect particular stimuli and initiate the formation of the PANoptosome (18, 19). Recent studies indicate that PANoptosis is significantly linked to the onset and progression of infectious diseases and cancer (20), thereby warranting a thorough exploration of the mechanisms underlying PANoptosis to formulate effective treatment strategies (21). In this study, we analyzed the role of PANoptosis-related genes in the prognosis and immunomodulation of CRC.
2 Methods
2.1 Data retrieval
The Bulk of colorectal adenocarcinoma (COAD) and colorectal adenocarcinoma (READ) were downloaded from NHGRI databases RNA-seq and batch removal via the “sva” packet combat function are finally merged into the CRC data set. We downloaded Bulk RNA-seq (GSE41328, GSE22598, GSE41328, GSE22598) from NCBI and NIH for colorectal cancer. GSE23878 and based on the “IOBR” R package (https://github.com/IOBR/IOBR) remove_batcheffect function merge is used to calculate carcinoma and genetic variations. In addition, Bulk RNA-seq (GSE39582, GSE72970, GSE17536) obtained from this data source contained prognostic information to assess the efficacy of the prognostic model. Single-cell RNA-Seq (scRNA-Seq) data related to CRC was obtained from the TISCH database (http://tisch.comp-genomics.org/home/), and included the GSE146771, EMTAB8107, and GSE166555 datasets (22). The missing samples were eliminated. In case of multiple rows for a gene in the expression matrix, the average was calculated. All data sources used in this study are publicly available and can be accessed without additional ethical approval. Pan-apoptotic genes were obtained from a supplementary dataset from literature (23). The research flow is shown in Figure 1.
Figure 1. The flow chart of this study.
2.2 Biological and pathway analysis
Based on the “limma” R package, we identified cancer and paratancerecancer differential genes (|logFC|> 1&p.Val < 0.05) in the combined CRC dataset and cancer and paratancerecancer differential genes (|logFC|> 1&p. Val < 0.05) in the GEO dataset, and performed heat map visualization. Based on the intersection of differential genes from two datasets, the common differential expression panapoptosis genes (CPAN_DEGs) were obtained, and then the protein interaction network (PPI) of these genes was constructed (24, 25).
2.3 Unsupervised clustering and correlation analysis
The “ConsensusClusterPlus” R package was used to conduct an unsupervised clustering analysis on the CPAN_DEGs. Proportional fuzzy clustering (PAC) and cumulative distribution function (CDF) curves (26) were used to find the most suitable group (k=2 in this analysis), and the prognosis of the two clusters was compared. The DEGs between the two clusters were identified using the “limma” R package (|logFC|>1 & p. Val < 0.05). CIBERSORT and ESTIMATE algorithms were used to analyze the immune microenvironment and immune checkpoint expression in both clusters.
2.4 Construction and verification of prognostic model
The prognostic value of the DEGs was evaluated in the CRC cohort by univariate Cox analysis (27). Through a systematic alteration of the regularization parameter log lambda, we observe the convergence pattern of Binomial Deviance to identify the ideal lambda value, which reflects the best equilibrium between the complexity of the model and its predictive efficiency. Based on this optimization parameter, the coefficient path plot of gene coefficient variation with L1 norm was plotted, and 15 candidate genes with significant prognostic value were incorporated in the multivariate Cox proportional risk model. The coefficients were calculated by the maximum partial likelihood estimation method to quantify their independent predictive power for survival outcomes. Based on the variable compression characteristics of LASSO regression and the results of survival analysis in the Cox model, the “CPAN-index” prognostic model was established. The coefficient distribution of model genes was determined by a lollipop diagram. The predictive performance of the model was tested on the GSE39582, GSE17536, and GSE72970 datasets by receiver operating characteristic (ROC) curve analysis, and the area under the curve was calculated. In addition, the association between the individual risk scores and survival duration was determined for each dataset.
2.5 Analysis of immune characteristics
The immune profile of the CRC dataset was analyzed using the CIBERSORT algorithm (28). To investigate the variations among CRC samples categorized as h-risk and l-risk, we initially conducted an analysis of the merged data set for CRC. Through the data set to submit to the TIDE’s, we get the samples of each CRC score or express quantity on several key indicators, The dysfunction score, CD274 expression, IFNG, TIDE score, Merck18, and CD8 were included and visualized to observe the immunosuppression characteristics difference. Gene-set enrichment analysis (GSEA) was performed using the “GseaVis” package (29). The initial key pathways were delineated to illustrate the enrichment of genes unique to both the high-risk and low-risk categories. Furthermore, the expression of immune checkpoint genes in both risk groups were also analyzed.
2.6 Single cell analysis
Three CRC scRNA-seq datasets were integrated, followed by batch removal through the “harmony” R package, and automatic annotation using “SingleR” data provided by TISCH. The single cell data was reduced to two dimensions using the Uniform Manifold Approximation and Projection (UMAP) algorithm, and the cell distribution of Seurat clusters, cell type and each data set was visualized. The genes with significant variations in expression across various cell types were screened using the FindMarkers function in Seurat. The upregulated and downregulated genes within each primary cell type were visualized using a Volcano plot. The AddModuleScore function in Seurat was used to score the predefined CPAN-index gene modules in the single-cell dataset in order to quantify the expression of CPAN-index genes in each cell. Furthermore, the proportion of CPAN-index positive cells in the different populations were calculated.
2.7 Cell culture and transfection
The normal colon epithelial cell line HIEC-6, and the CRC lines LoVo, HT-29, HCT116 and SW480 were obtained from Cell Bank of the Chinese Academy of Sciences. The HIEC-6 cells were cultured in OGM BulletKit medium (Lonza, Switzerland). HCT116 and SW480 cell lines were transfected with siRNA constructs specific for the CDKN2A gene, and a negative control (NC) construct (Sangon, China). Briefly, once the cells were 80% confluent in the 6-well plate, the lipo3000 transfection reagent and siRNA are diluted using Thermo’s (USA) Opti-MEM reduced serum medium. After 5 minutes, mix the two solutions and leave for another 20 minutes. The prepared transfer solution was introduced into the six-well plate, and the culture medium was refreshed 5 hours after transfection.
2.8 Total RNA extraction and RT-qPCR
The cells were harvested and lysed with 950 μL Trizol reagent (Takara, Japan) in the presence of RNase inhibitors. After incubating for 5 minutes, 150 μL chloroform (Sinopharm Holdings, China) was added to each tube, and the samples were mixed by vortex oscillation and centrifuged for 5 minutes. The supernatant was collected, and isopropyl alcohol (SINOPHARM, China) was added to precipitate the RNA. The samples were centrifuged for 5 minutes, and the RNA pellets were washed with 1 mL 75% or anhydrous ethanol, and air dried. The RNA concentration, and the level of DNA and protein contamination were measured. Following removal of genomic DNA, the RNA samples were reverse transcribed to cDNA using the PrimeScript RT kit (TaKaRa, Japan) as per the manufacturer’s guidelines. RT-qPCR was performed using the SYBR GreenER Supermix (TaKaRa, Japan) on the Roche480 PCR system (Roche, Switzerland).
2.9 Cell counting kit-8 assay
One day post-transfection, the cells were seeded into 96-well plates at the density of 4000 cells per well, and three replicate wells were set up for each group. At pre-determined time points, the CCK8 solution (KeyGEN, China) diluted in complete medium to a final volume of 100 µl was added to each well. The cells were incubated for 1.5 hours, and the absorbance of each well at 450nm was measured using a microplate reader.
2.10 Transwell migration assay
After transfection, we digested cells to harvest them and resuspended in serum-free medium. To ensure result accuracy, we seeded cells at 40,000 per well, with 250 μl suspension in each chamber, and added 550 μl complete culture medium to each 24 well plate well. Then we put chambers into plates and incubated for 24 hours. After that, we removed chamber fluid, wiped off non invaded cells, fixed cells with paraformaldehyde for 30 minutes, washed chambers twice with PBS, stained cells with 0.1% crystal violet for 20 minutes, washed again, and let them air dry. Finally, we took chamber images under a microscope for further analysis.
2.11 Flow cytometry
We employed flow cytometry to measure apoptosis. First, the cells were digested using EDTA free trypsin (Beyotime, China) and then centrifuged at 2,000 rpm to collect the cell pellets. Next, each group of cells was resuspended in 500 μL of Binding Buffer (Biosharp, China) and transferred into a flow tube. Following the manufacturer’s instructions, each group was supplemented with 5 μL of propidium iodide (PI, Biosharp, China) and 5 μL of fluorescein 5-isothiocyanate (FITC, Biosharp, China). After staining the cells in a 37°C water bath for 15 minutes, they were analyzed by flow cytometry. This entire process was repeated three times for each group.
2.12 Wound healing assay
HCT116 and SW480 cells were seeded in a 6-well plate (5×105 cells per well). After 24 h incubation, a vertical scratch lines was generated in each well. After washing with PBS, serum-free medium was added for cell culture. The images were taken at 0 h and 48 h. The experiment was repeated thrice.
2.13 Statistical analysis
All data in this study were based on at least 3 biological replicates, and all statistical studies used R software (version 4.4.1). The t test was used for comparison between groups, and continuous data were expressed as mean ± standard deviation (mean ± SD). P values < 0.05 were deemed statistically significant, and all P values were tested in a two-sided manner.
3 Results
3.1 Identification and functional annotation of differentially expressed PANoptosis genes
The concept of PANoptosis has been outlined in Figure 2a. The apoptosis-related genes comprised 87.4% of the PANoptosis gene list, whereas 9.7% and 2.9% of the genes were respectively associated with pyroptosis and necroptosis (Figure 2b). The top 20 DEGs between the CRC and normal colon tissues in the CRC and combined GSE databases are shown in Figures 2c–d. The overlap of DEGs in the CRC and GSE datasets with the PANoptosis gene list was determined by plotting Venn diagrams, and a total of 908 DEGs were common to both datasets, of which 15 genes were part of the PANoptosis list (Figure 2e). The protein-protein network (PPI) of these genes was constructed, and as shown Figure 2f, TNFRSF10B, FAS, CAV1, and PTPN1 were among the core genes of the network. The CPAN_DEGs were functionally annotated on the basis of Gene Ontology (GO) terms and KEGG pathways. The significantly enriched terms related to biological processes (GO-BP) included “enhancing proteolysis,” “promoting hydrolase activity,” and “facilitating cysteine-type endopeptidase activity”, and in particular “positively regulating proteolysis” and “positively regulating hydrolase activity” (p.adj < 0.0001, Figure 3a). Furthermore, “protein kinase inhibitor activity” and “kinase inhibitor activity” were among the molecular function terms (GO-MF) showing significant enrichment (p.adj < 0.01; Figure 3b). Analysis of the KEGG database revealed that multiple pathways related to apoptosis, immune response and viral infection were significantly enriched in the CPAN_DEGs, and the “apoptosis” pathway showed the highest enrichment. In addition, the “platinum drug resistance”, “p53 signaling pathway”, “viral carcinogenic effect” and other pathways also showed significant enrichment, while those associated with “human papillomavirus infection” and “TNF signaling pathway” were also enriched to some extent (Figure 3c).
Figure 2. Identification of PANoptosis-related DEGs (CPAN_DEGs) in CRC. (A) PANoptosis concept map (created by BioRender). (B) The proportion of apoptosis, necrosis and pyroptosis-related genes in the PANoptosis gene list. (C, D) Heat maps of the top 20 up-regulated and down-regulated genes between CRC and normal tissue in TCGA and combined GEO databases. (E) Venn map showing the overlap of DEGs and PANoptosis gene set. (F) PPI network of differentially expressed PANoptosis genes (CPAN_DEGs).
Figure 3. GO/KEGG/GSEA enrichment analysis of CPAN_DEGs. (A, B) Significantly enriched GO-BP (A) and GO-MF (B) terms in CPAN_DEGs. (C) KEGG pathways significantly associated with CPAN_DEGs.
3.2 Unsupervised clustering and correlation analysis
Unsupervised clustering was performed on the CRC samples based on the CPAN_DEGs. Both intra-cluster consistency (Figure 4a) and smoothness under the CDF curve (Figure 4b) decreased with the increase in k value. As shown in the PAC diagram in Figure 4c, the value is lowest when k=2, indicating that the sample distribution is less affected by random perturbations. Accordingly, the CRC cohort was divided into the C1 cluster (n=418) and the C2 cluster (n=128) (Figure 4d). The Kaplan-Meier curves of the two clusters showed no significant difference in survival (p > 0.05, Figure 4e). However, there were obvious differences in the expression patterns of CPAN_DEGs between C1 and C2 (Figure 4f).
Figure 4. Clustering of CRC samples based on CPAN_DEGs. (A, B) Evaluation of mean consistency in clusters and area under CDF curve when k = 2 ~ 10. (C) The PAC score indicates the optimal number of subtypes in the CRC dataset. (D) The CRC cohort was divided into two CRC subtypes by consensus clustering. (E) Kaplan-Meier survival analysis of CPAN_DEGs clusters. (F) Heat maps showing the expression of CPAN_DEGs in different CRC clusters.
3.3 Analysis of immunophenotypes and mutations
We observed unique immunophenotypes and mutational landscapes across the CRC samples. As shown in Figure 5a, the immune score for the C2 cluster was considerably greater than that for the C1 cluster (p < 0.05); however, no significant difference was observed in the matrix score (Figure 5b, c). This suggests that the C2 cluster may have stronger immune activity. Furthermore, comparison of the ESTIMATE scores revealed higher microenvironment score in the C2 cluster compared to that of C1 cluster (p < 0.05) (Figure 5d), which further supports greater immunoreactivity in the former. However, the expression levels of immune checkpoints such as TIGIT, PDCD1LG2, PDCD1, LAG3, HAVCR2, CTLA4, and CD274 was significantly higher in the C2 cluster than in C1cluster, which may reflect an immunosuppressive state in C2 cluster (Figure 5e). The tumor mutation burden (TMB) of the CPAN_DEGs clusters were similar (Figure 5f).
Figure 5. Immunophenotyping of CPAN_DEGs clusters. (A) Waterfall map of the distribution of 24 immune cells in CRC clusters. (B, C) Immune scores and matrix scores in CRC clusters. (D) Microenvironment scores in the CRC clusters. (E) Ridge map showing differences in immune checkpoint molecules between the CRC clusters. (F) Tumor mutation load in the CRC clusters. (*p < 0.05).
3.4 Establishment and validation of PANoptosis-related prognostic model
GSN, CDKN2A, CAV1 and PDK4 were the risk factors of CRC, whereas the remaining CPAN_DEGs were protective factors (Figure 6a). Through LASSO regression, we selected 11 CPAN_DEGs, including CDKN2A, GZMB, PMAIP1, CASP5, ENDOD1, FAS, CASP7, CAV1, PDK4, CHP2, and PRKACB, to construct a prognostic model (Figure 6b). The coefficients of the CPAN-index genes in overall clinical outcome are shown in Figure 6c. The AUC of the model for predicting 1-, 3-, and 5-year survival of CRC patients in the training cohort were 0.62, 0.6, and 0.62 respectively, indicating good predictive performance (Figure 6d). The distribution of the CPAN index among the CRC patients, and the survival status of individual patients are shown in Figure 6e. Based on the expression patterns of the 11 model genes, the patients were stratified into the high-risk and low-risk groups (Figure 6e). The CPAN index showed exponential distribution across these three GEO validation sets, indicating a significant correlation between patient survival outcomes and the expression levels of the model genes (Figure 7a). Furthermore, patients classified as high-risk had significantly worse overall survival compared to the low-risk patients across all validation sets (p < 0.0001; Figure 7b), thereby confirming the robust predictive power of the CPAN-index model.
Figure 6. Construction of a PANoptosis risk score model (CPAN-index) for CRC based on CPAN_DEGs. (A) Forest maps showing the prognostic value of 15 CPAN_DEGs in the CRC cohort. (B) Selection of model CPAN_DEGs by 10-fold cross-validation Lasso regression analysis. (C) Lollipop chart showing the coefficients for each CPAN_DEG included in CPAN-index. (D) Receiver operator characteristic (ROC) curves of CPAN-index for 1-, 3-, and 5-year survival. (E) Heat maps of CAPN-index distribution, survival status, and expression levels of model genes.
Figure 7. Validation of CPAN-index in multiple CRC cohorts. (A) Thermal maps of CPAN-index distribution, survival status, and expression levels of model genes in the three GEO validation datasets. (B) Kaplan-Meier survival curves of high-risk and low-risk groups in the validation cohorts. (C) Box plots showing immune populations in the CPAN-index groups. (D) Immune-score, Stromal-score and ESTIMATE-score in the CPAN-index groups. (E) Correlation between CPAN-index and immune cell infiltration in CRC. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).
3.5 Immune profiling of CPAN-index risk groups
The risk groups also showed differences in immune infiltration. As shown in Figure 7c, the infiltration of CD4 memory resting T cells was notably diminished in the high-risk cohort compared to the low-risk group (p < 0.0001), M0 macrophages were more abundant in the high-risk samples relative to the low-risk samples (p < 0.001). However, no significant variations were observed in the immune scores, stromal scores, tumor purity, or ESTIMATE scores between the risk groups (Figure 7d). The CPAN index also correlated with multiple immune cell types, including follicular helper T cells, M1 macrophages, CD8 T cells, naïve B cells, regulatory T cells (Tregs), M2 macrophages, CD4 memory resting T cells, M0 macrophages, activated mast cells, resting NK cells, neutrophils, CD4 memory activated T cells, activated dendritic cells (DCs), CD4 naïve T cells, memory B cells, and eosinophils, gamma delta T cells. Furthermore, the monocytes, activated NK cells, plasma cells, resting DCs and resting mast cells were positively or negatively correlated with each other. We also observed a significant positive correlation between the T cell cluster and B cell cluster, specifically the memory B cells (p < 0.05), suggesting that these cells may work synergistically in mounting an immune response against CRC. In addition, M1 macrophages also showed a high positive correlation with M2 macrophages, which likely reflects their coordination in the inflammatory response (p < 0.05). In contrast, the relationship between NK cells and the T cell cluster was mostly negative (p < 0.05), suggesting antagonistic roles in immune response (Figure 7e). Taken together, the CPAN index is closely associated with the immune microenvironment of CRC.
The Dysfunction scores (Figure 8a), CD274 scores (Figure 8b), IFNG scores (Figure 8c), Merck18 scores (Figure 8e), and CD8 scores (Figure 8f) were similar in the two risk groups. Conversely, the TIDE score (Figure 8d) of the high-risk group was significantly higher than that of the low-risk group (p < 0.001), indicating a greater propensity of tumor immune escape in the former. Furthermore, the HALLMARK pathways related to the DEGs in the two risk groups were identified through GSEA (p < 0.0001, Figure 8g). The expression levels of TIGIT, PDCD1LG2, PD-L1, LAG3, HAVCR2, CTLA4 and CD274 were significantly different between the high-risk and low-risk samples (Figure 8h).
Figure 8. CPAN-index correlated with immune landscape and molecular heterogeneity of CRC. (A–F) TIDE, immune dysfunction, CD274, Merck18, CD8 and IFNγ scores in low-risk and high-risk groups. (G) Significantly enriched pathways in high-risk and low-risk groups. (H) Immune checkpoint expression in high-risk and low-risk groups. (*p < 0.05, **p < 0.01, ***p < 0.001).
3.6 Single cell analysis
We performed dimensional-reduction clustering from three different scRNA-seq cohorts (EMTAB107, GSE146771, and GSE166555) using UMAP. The distribution of cell types was largely similar among the cohorts, indicating that there are common characteristics in the immune landscape of the CRC microenvironment. (Figures 9a–c). The up-regulated and down-regulated genes in each major cell type were also screened. MS4A1, BANK1, VPREB3, HLA-DRA, and TNFRSF13C were among the significantly up-regulated genes in B cells, whereas IL7R, CD3D, TRAC, TRBC1, CD3E, etc. showed obvious upregulation in the CD4+ T cells. In addition, CD9, KRT19, ANAX2, KLRD1, and RRBP1 were differentially expressed in various cell types (Figure 9d). The prognostic model in single cell datasets, and the score distribution in the different populations is shown in Figure 9e. We also calculated the proportion of cells positive for CPAN-index in each population. As shown in Figure 9f, NK cells had the highest positive rate of CPAN-index (70.3%), followed by myofibroblasts (66.9%), proliferative T cells (64.6%) and CD8 T cells (64.4%). In contrast, only 25.1% of the malignant cells and 13.5% of the mast cells were positive for CPAN-index (Figure 9f).
Figure 9. CPAN-index in single cell analysis. (A–C) UMAP visualization of cells from three common CRC scRNA-seq cohorts. (D) Volcanic maps showing downregulated and upregulated genes in each population. (E) Characteristic gene expression at the single-cell level measured by AddModuleScore () function in Seurat. (F) Histogram showing the proportion of CPAN-index positive cells in each population.
3.7 CDKN2A regulates the proliferation and migration of CRC cells
Currently, research on the role of CDKN2A in CRC among model proteins is scarce, whereas the functions of other proteins have been well documented. This study focuses on CDKN2A as the research target to investigate its role in CRC. Firstly, we assessed its expression levels in four CRC cell lines and one normal colon epithelial cell line. The findings revealed that CDKN2A mRNA expression was notably higher in CRC cell lines compared to the control cell line (p < 0.01, Figure 10A). Western blot analysis further confirmed that CDKN2A protein expression was also significantly elevated relative to the control group (p < 0.0001, Figure 10B). Next, we constructed two si-CDKN2A nucleic acid systems to knock down CDKN2A expression in HCT116 and SW480 cell lines, laying the groundwork for subsequent experiments. We verified the effective knockdown efficiency of these two systems at both mRNA and protein levels through RTqPCR and Western blot assays (p < 0.0001, Figures 10C, D). CCK8 assay results indicated that after CDKN2A knockdown in the two CRC cell lines, absorbance decreased significantly, suggesting a marked inhibition of cell proliferation (p < 0.0001, Figures 10E, F). Flow cytometry demonstrated that the total apoptosis rate of cells in the two CDKN2A knockdown groups was significantly higher than that in the control group (p < 0.0001, Figure 11A). Transwell experiments confirmed a significant reduction in the invasion and migration capabilities of CRC cells following CDKN2A knockdown (p < 0.0001, Figure 11B). Scratch tests further corroborated the inhibitory effect of CDKN2A knockdown on cell migration (p < 0.01, Figure 11C).
Figure 10. Expression and functional analysis of CDKN2A in colorectal cancer cell lines. (A) Relative expression levels of CDKN2A mRNA in the HIEC-6, LoVo, HT-29, HCT116 and SW480 cell lines. (B) Relative expression levels of CDKN2A protein in the HIEC-6, LoVo, HT-29, HCT116 and SW480 cell lines. (C-D) Relative CDKN2A mRNA and protein expression in the HCT116 and SW480 cells transfected with si-NC, siCDKN2A-1, or siCDKN2A-2 constructs. (E) Proliferation rates of HCT116 cells transfected with si-NC, siCDKN2A-1, or siCDKN2A-2 constructs. The line graph shows OD values at 450 nm over a period of 4 days. (F) Proliferation rates of SW480 cells transfected with si-NC, siCDKN2A-1, or siCDKN2A-2 constructs. The line graph shows OD values at 450 nm over a period of 4 days. ***p <0.01, ****p < 0.0001..
Figure 11. Impact of CDKN2A knockdown on apoptosis, invasion, and migration of colorectal cancer (CRC) cells. (A) Flow cytometry analysis showing the total apoptosis rate of CRC cells in two CDKN2A knockdown groups (transfected with siCDKN2A-1 and siCDKN2A-2) compared to the control group (transfected with si-NC). (B) Transwell assay results depicting the invasion and migration capabilities of CRC cells following CDKN2A knockdown. (C) Wound healing results demonstrating the effect of CDKN2A knockdown on CRC cell migration. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
4 Discussion
A greater understanding of the pathological mechanisms of CRC is essential to design better treatment strategies (30). In recent years, the role of cell death pathways in tumorigenesis and development has received significant attention. In this study, we analyzed the clinical significance of PANoptosis, a complex process involving apoptosis, pyroptosis, and necroptosis, in CRC, and constructed a prognostic model based the PANoptosis gene signature. The apoptosis-related genes comprised the major chunk of the PANoptosis gene set in CRC, indicating a significant role in cancer development and progression. In addition, the differentially expressed PANoptosis genes between the CRC and normal tissues (CPAN_DEGs) were significantly enriched in pathways related to cell death and responses to metabolic stress, as well as autoimmune diseases, thereby presenting novel targets for developing effective therapies.
Based on the expression levels of CPAN_DEGs, the CRC samples were stratified into two distinct clusters, which differed in terms of gene expression profiles and immune cell distribution, but showed similar survival rates and TMB. Thus, unsupervised clustering based on 15 CPAN_DEGs could not distinguish patient outcomes, which can be attributed to the complexity of the PANoptosis pathway, and the co-existence of pro-death and pro-survival signals. Using LASSO dimensionality reduction, we selected 11 core CPAN_DEGs, including CDKN2A, CASP5 and PDK4 among others, to construct a prognostic model. This CPAN-index model accurately stratified CRC patients into the high-risk and low-risk groups, indicating superior predictive efficacy. In addition, the prognostic model showed consistent performance across multiple independent datasets, indicating that its generalization ability and clinical translation potential. Therefore, the expression patterns of these PANoptosis gene markers can help identify high-risk CRC patients and guide individualized treatment strategies.
The high-risk group showed significantly higher TIDE score compared to the low-risk group, which is indicative of stronger immune escape potential due to T cell dysfunction or enhanced matrix barriers in the former (31, 32). Furthermore, the enrichment of resting NK cells and unpolarized macrophages in the high-risk samples may reflect an immunosuppressive state. Although the overall immune score of the microenvironment was not lower, the immune function was inhibited. On the other hand, the elevation of resting CD4 memory T cells and plasma cells in the low-risk samples could be indicative of potential future immune responses. Interestingly, the indicators of immunotherapy responses, such as CD8+ T cell levels and PD-L1 (CD274) expression, remained comparable between the two groups. This suggests that the prognostic significance of the CPAN-index may operate independently of established markers, and is likely related to the functional polarization of immune cells (including M0/M2 equilibrium and suppression of NK cell activity), as well as non-classical immune evasion pathways such as pyroptosis-driven chronic inflammation. In conclusion, CPAN-index can reveal differences in the immune landscape that cannot be captured by traditional microenvironment scores, and provide a new perspective for prognosis assessment and combination therapies targeting the immunometabolism-death interactive network.
The high-risk group showed significant heterogeneity in multiple HALLMARK pathways related to tumor progression and microenvironment remodeling. The activation of epithelial mesenchymal transformation (EMT) and mitotic spindle pathways suggests that the tumor cells in high-risk samples have stronger metastasis ability that may drive aggressive phenotypes through EMT-related transcription factors (such as SNAI1 and ZEB1), along with abnormal mitotic processes that may lead to increased genomic instability (27, 33). Upregulation of interferon alpha response is usually indicative of immune stress induced by chronic viral infection or endogenous retroelement activation. However, when considering the elevated TIDE scores and functional silencing of NK cells in high-risk samples, the interferon signaling pathway may be associated with immune tolerance or T cell depletion due to persistent overactivation rather than an effective anti-tumor immune response. It is worth noting that the enrichment of myogenic pathways in the high-risk samples may be due to the abnormal activation of cancer-associated fibroblasts, which promote stromal sclerosis and an immunosuppressive TME by secreting factors such as TGF-β. On the other hand, downregulation of the KRAS signaling inhibitory pathway suggested continuous activation of KRAS driving signals, which may promote cell proliferation and drug resistance through the MAPK/ERK pathway. The general inhibition of metabolic pathways is consistent with tumor cells relying on glycolysis (Warburg effect) and lipid synthesis reprogramming to meet their high energy demands. This metabolic pressure may exacerbate the acidity in the TME and inhibit immune cell function. In addition, the downregulation of IL6-JAK-STAT3 signaling in the high-risk samples was inconsistent with the enrichment of M2-type macrophages, which indicates that non-classical activation of STAT3, such as through the EGFR or inflammasome pathway, is the dominant driver of cancer. In addition, decreased metabolic capacity may impair the detoxification of chemotherapy drugs and lead to treatment resistance. Taken together, the molecular characteristics of the high-risk samples correlate with the transfer-metabolism-immunosuppression phenotype.
Single-cell analysis showed that a higher proportion of NK cells, CD8+ T cells and myofibroblasts were positive for the CPAN-index compared to tumor cells. This suggests that immune cells could not effectively kill tumors due to overactivation that eventually led to apoptosis or pyroptosis. In addition, the myofibroblasts likely contributed to the formation of a physical barrier and immunosuppressive TME by secreting pro-inflammatory and pro-fibrotic factors, and their metabolic abnormalities related to PANoptosis further aggravated chronic inflammation and tumor metastasis. At the same time, tumor cells escape immune-mediated clearance by downregulating programmed death genes, and the inflammation and metabolic stress in the TME provided an additional survival advantage. This vicious cycle of immune depletion, matrix remodeling, and tumor escape provides a new direction for targeted interventions, such as combining immune checkpoint inhibitors with anti-fibrotic drugs, or regulating PANoptosis.
The CDKN2A gene plays a complex and crucial role in the occurrence, progression, and treatment response of CRC (34). Its functional abnormalities drive malignant tumor phenotypes through multiple mechanisms, including cell cycle dysregulation, epigenetic modification, immune microenvironment remodeling, and metabolic reprogramming. CDKN2A exhibits a double-edged sword effect in CRC. Low expression leads to a loss of its tumor-suppressive function, promoting cell proliferation and invasion; high expression, however, may promote tumor progression through non-classical pathways. Mechanistic studies have found that high doses of CDKN2A can enhance tumor cell migration by activating the Wnt/β-catenin pathway or promoting goblet cell mucus secretion. Furthermore, CDKN2A is involved in immune escape in CRC. For example, the long non-coding RNA SNHG26 inhibits CD8+ T cell infiltration by degrading CDKN2A mRNA, while simultaneously promoting regulatory T cell (Treg) differentiation, creating an immunosuppressive microenvironment. This immunomodulatory effect is closely related to CDKN2A expression levels—high expression allows tumor cells to secrete IL-10, which can induce Treg differentiation, further weakening the anti-tumor immune response. Epigenetic modification is an important mechanism for the dysfunction of CDKN2A in CRC. Abnormal DNA methylation leads to CDKN2A silencing, while histone modification enhances the transcription of pro-inflammatory factors and promotes tumor-associated inflammation.
PANoptosis plays a crucial role in the tumor heterogeneity, immune remodeling, and prognosis in CRC, and the high-risk patients identified by the CPAN-index model present a metastatic-metabolic-immunosuppressive phenotype. The CPAN-index accurately captured the core network driving CRC progression by integrating the heterogeneity in PANoptosis genes, thus providing a theoretical basis for developing combination therapy strategies targeting the immunometabolism-death interaction network. In the future, it is necessary to further analyze the mechanism of action of key genes (such as CDKN2A) and promote clinical translational applications.
Our study offers valuable insights into how the CDKN2A mediated PANoptosis signaling network affects CRC progression, but has limitations. First, results mainly come from in vitro experiments with colon cancer cell lines; lack of in vivo validation limits extrapolation to living organisms’ complex physiological environments, as in vivo studies can consider factors like tissue architecture and tumor microenvironment interactions not replicated in vitro. Second, we haven’t directly measured PANoptosis and immune microenvironment changes after CDKN2A perturbation. Though our design and analysis imply their involvement, direct measurement methods would provide more definitive evidence and enhance understanding of CDKN2A perturbation’s precise mechanisms. Third, our study heavily relies on public datasets, which, while valuable, have inherent limitations like batch effects that may confound results and clinical heterogeneity among samples, introducing variability that could impact the generalizability and robustness of findings.
5 Conclusion
PANoptosis plays a key role in driving the malignant progression of CRC by regulating tumor heterogeneity, immune remodeling, and metabolic reprogramming. The CPAN-index model could stratify the patients based on their prognosis, and the high-risk patients were characterized by the upregulation of pathways related to cancer invasion, metabolic reprogramming, and immunosuppression. While this study provides a new direction for accurate prognosis assessment and targeted intervention strategies in CRC, the mechanisms of CDKN2A and other key genes need to be elucidated further to promote the clinical application of the CPAN-index model.
Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.
Ethics statement
Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.
Author contributions
XH: Conceptualization, Writing – original draft. WW: Formal analysis, Writing – original draft. LL: Writing – original draft, Methodology. YY: Writing – original draft, Methodology. SD: Conceptualization, Writing – review & editing.
Funding
The author(s) declared financial support was received for this work and/or its publication. The authors declare that financial support was received for the research, authorship, and publication of this article. The Shanxi Provincial Health Commission Fund (grant number 2020070) gave funding for the study.
Acknowledgments
We thank the GEO database, KEGG database, and OmicShare tools for supporting the research.
Conflict of interest
The authors declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that Generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1. Ionescu VA, Gheorghe G, Bacalbasa N, Chiotoroiu AL, and Diaconu C. Colorectal cancer: from risk factors to oncogenesis. Med (Kaunas). (2023) 59:1646. doi: 10.3390/medicina59091646
2. Li Z, Mao L, Yu B, Liu H, Zhang Q, Bian Z, et al. GB7 acetate, a galbulimima alkaloid from Galbulimima belgraveana, possesses anticancer effects in colorectal cancer cells. J Pharm Anal. (2022) 12:339–49. doi: 10.1016/j.jpha.2021.06.007
3. Arnold M, Sierra MS, Laversanne M, Soerjomataram I, Jemal A, and Bray F. Global patterns and trends in colorectal cancer incidence and mortality. Gut. (2017) 66:683–91. doi: 10.1136/gutjnl-2015-310912
4. Lawal B, Wang YC, Wu ATH, and Huang HS. Pro-oncogenic c-met/EGFR, biomarker signatures of the tumor microenvironment are clinical and therapy response prognosticators in colorectal cancer, and therapeutic targets of 3-phenyl-2H-benzo[e][1,3]-oxazine-2,4(3H)-dione derivatives. Front Pharmacol. (2021) 12:691234. doi: 10.3389/fphar.2021.691234
5. Wang N, Lu Y, Khankari NK, Long J, Li HL, Gao J, et al. Evaluation of genetic variants in association with colorectal cancer risk and survival in Asians. Int J Cancer. (2017) 141:1130–9. doi: 10.1002/ijc.30812
6. Kuipers EJ, Grady WM, Lieberman D, Seufferlein T, Sung JJ, Boelens PG, et al. Colorectal cancer. Nat Rev Dis Primers. (2015) 1:15065. doi: 10.1038/nrdp.2015.65
7. Vasen HF, Tomlinson I, and Castells A. Clinical management of hereditary colorectal cancer syndromes. Nat Rev Gastroenterol Hepatol. (2015) 12:88–97. doi: 10.1038/nrgastro.2014.229
8. Poursheikhani A, Abbaszadegan MR, Nokhandani N, and Kerachian MA. Integration analysis of long non-coding RNA (lncRNA) role in tumorigenesis of colon adenocarcinoma. BMC Med Genomics. (2020) 13:108. doi: 10.1186/s12920-020-00757-2
9. Barresi V. Colorectal cancer: from pathophysiology to novel therapeutic approaches. Biomedicines. (2021) 9:1858. doi: 10.3390/biomedicines9121858
10. Li Q, Geng S, Luo H, Wang W, Mo YQ, Luo Q, et al. Signaling pathways involved in colorectal cancer: pathogenesis and targeted therapy. Signal Transduct Target Ther. (2024) 9:266. doi: 10.1038/s41392-024-01953-7
11. Veldwijk J, Lambooij MS, Kallenberg FG, van Kranen HJ, Bredenoord AL, Dekker E, et al. Preferences for genetic testing for colorectal cancer within a population-based screening program: a discrete choice experiment. Eur J Hum Genet. (2016) 24:361–6. doi: 10.1038/ejhg.2015.117
12. Jakubowska K, Koda M, Grudzińska M, Kańczuga-Koda L, and Famulski W. Monocyte-to-lymphocyte ratio as a prognostic factor in peripheral whole blood samples of colorectal cancer patients. World J Gastroenterol. (2020) 26:4639–55. doi: 10.3748/wjg.v26.i31.4639
13. Ashrafizadeh M, Farhood B, Eleojo Musa A, Taeb S, Rezaeyan A, and Najafi M. Abscopal effect in radioimmunotherapy. Int Immunopharmacol. (2020) 85:106663. doi: 10.1016/j.intimp.2020.106663
14. Ning J, Wang Y, and Tao Z. The complex role of immune cells in antigen presentation and regulation of T-cell responses in hepatocellular carcinoma: progress, challenges, and future directions. Front Immunol. (2024) 15:1483834. doi: 10.3389/fimmu.2024.1483834
15. Wang Y, Wang F, Pan Q, and Wang J.Integrating multi-omics proteomic approaches in deciphering the tumor microenvironment and therapeutic resistance mechanisms in prostate cancer. Biochim Biophys Acta Rev Cancer (2025) 1880:189426. doi: 10.1016/j.bbcan.2025.189426
16. Huang J, Jiang S, Liang L, He H, Liu Y, Cong L, et al. Analysis of PANoptosis-related lncRNA-miRNA-mRNA network reveals lncRNA SNHG7 involved in chemo-resistance in colon adenocarcinoma. Front Oncol. (2022) 12:888105. doi: 10.3389/fonc.2022.888105
17. Wang Y, Yan K, Guo Y, Lu Y, Su H, and Li H. IP-score correlated to endogenous tumour antigen peptide processing: A candidate clinical response score algorithm of immune checkpoint inhibitors therapy in multiple cohorts. Front Immunol. (2023) 13:1085491. doi: 10.3389/fimmu.2022.1085491
18. Karki R, Sundaram B, Sharma BR, Lee S, Malireddi RKS, Nguyen LN, et al. ADAR1 restricts ZBP1-mediated immune response and PANoptosis to promote tumorigenesis. Cell Rep. (2021) 37:109858. doi: 10.1016/j.celrep.2021.109858
19. Ning J, Wang Y, and Tao Z. The complex role of immune cells in antigen presentation and regulation of T-cell responses in hepatocellular carcinoma: progress, challenges, and future directions. Front Immunol. (2024) 15:1483834. doi: 10.3389/fimmu.2024.1483834
20. Liu F, Li Z, He X, Yu H, and Feng J. Ghrelin attenuates neuroinflammation and demyelination in experimental autoimmune encephalomyelitis involving NLRP3 inflammasome signaling pathway and pyroptosis. Front Pharmacol. (2019) 10:1320. doi: 10.3389/fphar.2019.01320
21. Shi C, Cao P, Wang Y, Zhang Q, Zhang D, Wang Y, et al. PANoptosis: A cell death characterized by pyroptosis, apoptosis, and necroptosis. J Inflammation Res. (2023) 16:1523–32. doi: 10.2147/jir.S403819
22. Zhang S, Yu H, Li J, Zhao L, Tan L, Song Q, et al. Identification of prognostic and tumor microenvironment by shelterin complex-related signatures in oral squamous cell carcinoma. Oxid Med Cell Longev. (2022) 2022:6849304. doi: 10.1155/2022/6849304
23. Song F, Wang CG, Mao JZ, Wang TL, Liang XL, Hu CW, et al. PANoptosis-based molecular subtyping and HPAN-index predicts therapeutic response and survival in hepatocellular carcinoma. Front Immunol. (2023) 14:1197152. doi: 10.3389/fimmu.2023.1197152
24. Zheng JN, Li Y, Yan YM, Yu Y, Shao WQ, and Wang Q. Increased serum calpain activity is associated with HMGB1 levels in systemic sclerosis. Arthritis Res Ther. (2020) 22:110. doi: 10.1186/s13075-020-02195-y
25. Yan H, Gao Y, Ding Q, Liu J, Li Y, Jin M, et al. Exosomal micro RNAs derived from dermal papilla cells mediate hair follicle stem cell proliferation and differentiation. Int J Biol Sci. (2019) 15:1368–82. doi: 10.7150/ijbs.33233
26. Liu Z, Liu L, Weng S, Guo C, Dang Q, Xu H, et al. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat Commun. (2022) 13:816. doi: 10.1038/s41467-022-28421-6
27. Bai Y, Gong Y, Bai J, Liu J, Deng HW, Calhoun V, et al. A joint analysis of multi-paradigm fMRI data with its application to cognitive study. IEEE Trans Med Imaging. (2021) 40:951–62. doi: 10.1109/tmi.2020.3042786
28. Weis-Banke SE, Lisle TL, Perez-Penco M, SChina A, Hübbe ML, Siersbæk M, et al. Arginase-2-specific cytotoxic T cells specifically recognize functional regulatory T cells. J Immunother Cancer. (2022) 10:e005326. doi: 10.1136/jitc-2022-005326
29. Zhang J, Li H, Tao W, and Zhou J. GseaVis: An R Package for Enhanced Visualization of Gene Set Enrichment Analysis in Biomedicine. Med Res 1:131–5. doi: 10.1002/mdr2.70000
30. Ye SP, Qiu H, Liao SJ, Ai JH, and Shi J. Mini-invasive vs open resection of colorectal cancer and liver metastases: A meta-analysis. World J Gastroenterol. (2019) 25:2819–32. doi: 10.3748/wjg.v25.i22.2819
31. Xia Y, Zhang R, Wang M, Li J, Dong J, He K, et al. Development and validation of a necroptosis-related gene prognostic score to predict prognosis and efficiency of immunotherapy in gastric cancer. Front Immunol. (2022) 13:977338. doi: 10.3389/fimmu.2022.977338
32. He X, Wang J, Yu H, Lv W, Wang Y, Zhang Q, et al. Clinical significance for diagnosis and prognosis of POP1 and its potential role in breast cancer: a comprehensive analysis based on multiple databases. Aging (Albany NY). (2022) 14:6936–56. doi: 10.18632/aging.204255
33. Yang L, Jin WQ, Tang XL, Zhang S, Ma R, Zhao DQ, et al. Ginseng-derived nanoparticles inhibit lung cancer cell epithelial mesenchymal transition by repressing pentose phosphate pathway activity. Front Oncol. (2022) 12:942020. doi: 10.3389/fonc.2022.942020
Keywords: CDKN2A, colorectal cancer (CRC), immune microenvironment, PANoptosis, single-cell RNA sequencing
Citation: He X, Wang W, Li L, Yin Y and Ding S (2026) Unveiling the immunometabolic landscape of colorectal cancer through PANoptosis-related gene expression. Front. Immunol. 16:1615022. doi: 10.3389/fimmu.2025.1615022
Received: 20 April 2025; Accepted: 12 December 2025; Revised: 30 November 2025;
Published: 12 January 2026.
Edited by:
Zheng Liu, Beijing Arion Cancer Hospital, ChinaCopyright © 2026 He, Wang, Li, Yin and Ding. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Shunbin Ding, NDIwMzg0MTA5QHFxLmNvbQ==