Prognostic value of genes associated with metastasis and propionate metabolism in rectal cancer

Introduction

Colorectal cancer (CRC) is a prevalent malignant tumor of the digestive tract, with peak incidence occurring in individuals over 50 years of age (1,2). CRC development is strongly associated with familial genetic predisposition, dietary patterns, and lifestyle factors. Rectal cancer, as a subtype of CRC, is characterized by high morbidity and mortality, with survival outcomes varying markedly across disease stages. Current treatment strategies for rectal cancer adopt a multimodal approach, including surgical resection, radiotherapy, chemotherapy, and increasingly, targeted therapy and immunotherapy (3,4). However, prognosis remains poor in cases involving distant metastasis, primarily due to challenges in achieving complete response and effective disease control. Despite advances in systemic therapies such as targeted agents and immunotherapy, the management of metastatic disease continues to present substantial challenges that adversely affect long-term survival and quality of life. Established prognostic biomarkers, including KRAS and BRAF mutations, microsatellite instability-high (MSI-H) status, and carcinoembryonic antigen (CEA) levels (5-7), are widely used for prognostic stratification and therapeutic decision-making in rectal cancer. Nevertheless, their predictive performance remains limited; therefore, the identification of novel biomarkers is essential to improve personalized treatment strategies and overall survival (OS) outcomes.

Metabolic reprogramming represents a pivotal strategy by which cancer cells acquire the capabilities required for metastasis during tumor progression. Propionate metabolism encompasses the synthesis and catabolic processes of propionate in biological systems. It is primarily involved in energy metabolism and fatty acid synthesis and plays an important role in maintaining intestinal microecological homeostasis. Evidence suggests that propionate metabolism regulates cellular proliferation and inflammatory responses during tumorigenesis and may influence cancer progression through modulation of the tumor microenvironment. In rectal cancer, propionate metabolism has been implicated in lymphatic metastasis, with associated metabolic alterations potentially enhancing tumor cell invasiveness (8,9). Notably, hyperlipoproteinemia has been reported to promote lymphangiogenesis and lymph node metastasis in CRC by inducing intestinal microbiota dysbiosis and inhibiting the GPR41 signaling pathway (10). However, most existing studies have focused on the macroscopic regulatory effects of propionate metabolism, and substantial gaps remain in the investigation of PMRGs in rectal cancer. On one hand, the specific functional associations and molecular mechanisms linking PMRGs to lymph node metastasis in rectal cancer remain incompletely elucidated. On the other hand, systematic identification of prognostic biomarkers and construction of risk prediction models based on PMRGs are still lacking. Accordingly, the present study aims to address these gaps through bioinformatics analyses, systematically identify prognostic biomarkers associated with both lymph node metastasis and propionate metabolism, and establish a robust risk prediction model for rectal cancer.

In this study, transcriptomic and clinical data related to lymph node metastasis in rectal cancer were retrieved from public databases, together with genes associated with propionate metabolism. Bioinformatics approaches were employed to identify prognostic genes associated with lymph node metastasis and propionate metabolism and to evaluate their prognostic value. In addition, predictive models were constructed and validated to investigate the association between risk stratification and immunotherapy response, as well as to explore the relationship between prognostic genes and drug sensitivity. These findings provide important insights for clinical outcome prediction and therapeutic strategy development in patients with rectal cancer. The analytical workflow of this study is illustrated in Figure 1. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1061/rc).

Figure 1 The research analysis workflow. DEGs, differentially expressed genes; GSEA, Gene Set Enrichment Analysis; LASSO, least absolute shrinkage and selection operator; TCGA-READ, The Cancer Genome Atlas-Rectum Adenocarcinoma.

Methods

Source of data

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.The Cancer Genome Atlas-Rectum Adenocarcinoma (TCGA-READ) cohort was obtained from the University of California Santa Cruz (UCSC) Genome Browser repository (https://genome.ucsc.edu/), comprising 163 rectal adenocarcinoma samples, including 82 cases without lymph node involvement (N0) and 77 cases with lymph node involvement (N1–N2), together with RNA sequencing data from 10 normal tissue samples, which served as the training cohort. Among the N1–N2 cases, 73 samples contained available survival information. The GSE87211 cohort was retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds). Based on the GPL13497 platform, this dataset included transcriptomic profiles from 129 rectal cancer samples and was used as the validation cohort. In addition, 604 propionate metabolism-related genes (PMRGs) with correlation scores >7 were obtained from the GeneCards database.

The identification and functional enrichment of potential genes

Differentially expressed genes (DEGs) were identified from the TCGA-READ cohort using two comparisons: DEGs1 (tumor versus normal tissue) and DEGs2 (N0 versus N1–N2 stage), with thresholds of P<0.05 and |log₂ fold change (FC)| >0.5. Differential expression analysis was performed using the edgeR package (version 3.36.0) (11). Volcano plots were generated using the ggplot2 package (version 3.3.5) (12), and heatmaps were constructed using the pheatmap package (version 1.0.12). Overlapping genes between DEGs1 and DEGs2 were defined as crossover genes, and Pearson correlation analysis was applied to evaluate their associations with PMRGs. Candidate genes were selected based on adjusted P<0.05 and |r|>0.4. Functional annotation of candidate genes was performed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, with statistical significance defined as adjusted P<0.05.

Screening biomarkers and constructing risk models

Prognostic biomarkers were identified from candidate genes using univariate Cox regression (13) combined with least absolute shrinkage and selection operator (LASSO) regularization. Patients in the training cohort were stratified into high- and low-risk groups based on the median risk score derived from biomarker expression levels. The risk score was calculated as: Risk score = ∑1ncoef(genei)∗expression(genei). Kaplan-Meier survival curves were generated to assess survival differences between groups. The survivalROC package (version 1.0.3) (14) was used to compute area under the curve (AUC) values for evaluation of model predictive performance. External validation was performed using the GSE87211 dataset to confirm model robustness.

Independent prognostic analysis

The TCGA-READ dataset was further used for independent prognostic analysis of clinical variables, including age, sex, metastasis (M) stage, cancer status, node (N) stage, tumor (T) stage, as well as the risk score. Univariate Cox regression analysis was first conducted, and variables with P<0.05 were retained for further evaluation. These variables were then subjected to proportional hazards (PH) assumption testing using the cox.zph function in the survival package to assess suitability for multivariate Cox regression. Only variables satisfying the PH assumption were included in the multivariate Cox regression model to identify independent prognostic factors. Subsequently, a nomogram was constructed and visualized using the rms package (v 6.2-0) (15). Finally, calibration curves and ROC curves were generated to evaluate model performance.

Screening and functional enrichment of risk-related DEGs (RRDEGs)

Differential gene expression analysis between the two risk stratification groups was performed in the TCGA-READ cohort using the edgeR package (version 3.36.0) (11), with thresholds of P<0.05 and |log₂FC| >0.5 applied to identify RRDEGs. Following identification, RRDEGs were subjected to functional annotation through GO and KEGG pathway enrichment analyses using the clusterProfiler package (16), with statistical significance set at P<0.05.

Immune checkpoint and immunotherapy sensitivity analysis

Expression levels of eight immune checkpoint genes between the two risk subgroups were compared using the Wilcoxon test, and results were visualized as box plots. Spearman’s rank correlation analysis was performed to assess associations between risk scores and differentially expressed immune checkpoints, with significance defined as P<0.05 and |correlation coefficient (cor)| >0.3. In addition, immunotherapy sensitivity analysis was conducted using a SubMap algorithm, and a heatmap illustrating responsiveness to immune checkpoint inhibitors across risk subgroups was generated based on comparative results.

Single-gene Gene Set Enrichment Analysis (GSEA) of biomarkers

Single-gene GSEA was performed using the clusterProfiler package (16) to identify enriched regulatory pathways and biological functions associated with the biomarkers. The top five GO biological process (BP) and KEGG pathways for each biomarker were subsequently ranked and visualized based on adjusted P values.

Drug sensitivity analysis

Drug sensitivity for each sample was evaluated using the pRRophetic algorithm based on the GDSC database. The Wilcoxon test was then applied to compare half-maximal inhibitory concentration (IC₅₀) values between risk groups, and results were visualized using box plots.

Correlation of propionate metabolism pathway score with risk score, clinical outcomes, and immune infiltration in TCGA-READ

To assess overall activity of the propionate metabolism pathway in TCGA-READ, PMRGs were defined as the gene signature of this pathway. Pathway activity scores for each sample were calculated using the single-sample gene set enrichment analysis (ssGSEA) method implemented in the gsva() function of the GSVA package (v 1.50.0)(17). Spearman correlation analysis was conducted using the psych package (v 2.4.6.26) (18) to evaluate associations between PMRG pathway activity scores, risk scores, and immune checkpoints. Correlations with P<0.05 and |cor|>0.3 were considered statistically significant. Furthermore, ssGSEA was applied to all samples in the high- and low-risk groups of the TCGA-READ cohort to estimate the relative infiltration levels of 28 immune cell types. Differences in immune cell infiltration between groups were assessed using the Wilcoxon test (P<0.05). Correlations between differentially infiltrated immune cells and PMRG pathway activity scores were further evaluated (P<0.05 and |cor|>0.3). To assess the association between PMRG pathway activity scores and OS in rectal cancer patients, 73 patients were stratified into high- and low-score groups based on the median PMRG pathway activity score. Survival differences were evaluated using the log-rank test from the survminer package (v 0.4.9)(19), and Kaplan-Meier survival curves were generated (P<0.05).

Statistical analysis

All bioinformatics analyses were performed in R, and the Wilcoxon test was used for differential comparisons.

Results

Screening and functional enrichment of candidate genes

Differential expression analysis identified 5,579 DEGs1 between rectal adenocarcinoma and normal tissue samples, including 2,829 upregulated and 2,750 downregulated genes (Figure 2A,2B, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-1.xls). In parallel, comparison between N1–N2 and N0 staging groups identified 30 DEGs2, comprising 14 upregulated and 16 downregulated genes (Figure 2C,2D, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-2.xls). Intersection analysis between DEGs1 and DEGs2 yielded 22 overlapping genes (Figure 2E, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-3.xls). Subsequently, correlation analysis identified 157 candidate genes for further investigation (table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-4.xls). Functional enrichment analysis demonstrated that these candidates were involved in 1,653 BP terms, 49 cellular components (CC) terms, 112 molecular functions (MF) terms, and 80 KEGG pathways. GO annotation indicated predominant involvement in inflammatory response regulation, lipopolysaccharide response, and related processes (Figure 2F). KEGG pathway enrichment highlighted tumor necrosis factor (TNF) signaling pathways, cytokine-cytokine receptor interaction, among other pathways (Figure 2G, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-5.xls).

Figure 2 Screening and functional enrichment analysis of candidate genes associated with lymph node metastasis and propionate metabolism in rectal cancer. (A,B) Volcano plot and heatmap showing differential gene expression between rectal adenocarcinoma and normal tissue samples (DEGs1: 5,579 genes, with 2,829 upregulated and 2,750 downregulated genes). (C,D) Volcano plot and heatmap displaying differential gene expression between N1–N2 and N0 staging groups (DEGs2: 30 genes, with 14 upregulated and 16 downregulated genes). (E) Venn diagram illustrating the intersection of DEGs1 and DEGs2, yielding 22 overlapping genes. (F) GO enrichment analysis of 157 candidate genes, showing predominant involvement in inflammatory response regulation and lipopolysaccharide response mechanisms across BP, CC, and MF categories. (G) KEGG pathway enrichment analysis highlighting significant enrichment in TNF signaling cascades and cytokine-cytokine receptor interactions. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Biomarker screening and risk model

Through univariate Cox regression combined with LASSO analysis, five biomarkers associated with both lymph node metastasis and propionate metabolism (CCL24, IGFBP3, ODC1, PYGM, and VKORC1) were identified (Figure 3A-3C). Using the median risk score (2.788947269) as the cutoff, participants were stratified into high-risk (n=37) and low-risk (n=36) groups (Figure 3D,3E). Kaplan-Meier analysis showed that patients in the high-risk group had significantly poorer survival outcomes compared with those in the low-risk group within the training cohort (Figure 3F). Receiver operating characteristic (ROC) analysis demonstrated robust prognostic performance, with AUC values exceeding 0.70 for 1-, 3-, and 5-year survival prediction (Figure 3G). To further evaluate predictive robustness, an independent validation cohort was used for external validation. Risk distribution patterns and survival trajectories in the validation dataset were consistent with those observed in the training cohort (Figure 4A-4C). Notably, all AUC values in the validation cohort exceeded 0.63, indicating acceptable predictive performance across different follow-up time points (Figure 4D).

Figure 3 Identification of prognostic biomarkers and construction of risk prediction model in the training cohort. (A,B) Univariate Cox regression analysis and LASSO coefficient profiles identifying five biomarkers associated with lymph node metastasis and propionate metabolism. (C) Cross-validation curves for optimal lambda selection in LASSO regression. (D,E) Risk score distribution and survival status of patients in high-risk (n=37) and low-risk (n=36) groups stratified by median risk score (2.788947269). (F) Kaplan-Meier survival curves demonstrating significantly poorer OS in the high-risk group compared to the low-risk group. (G) Time-dependent ROC curves showing robust predictive performance with AUC values exceeding 0.70 for 1-, 3-, and 5-year survival prediction. AUC, area under the curve; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic.

Figure 4 External validation of the risk prediction model in the GSE87211 cohort. (A,B) Risk score distribution and survival status of patients in the validation cohort, demonstrating patterns consistent with the training data. (C) Kaplan-Meier survival curves confirming significantly different survival outcomes between high-risk and low-risk groups in the validation cohort. (D) Time-dependent ROC curves showing adequate predictive precision with all AUC values surpassing 0.63 across various follow-up intervals (1-, 3-, and 5-year survival). AUC, area under the curve; ROC, receiver operating characteristic.

Construction of the nomogram and assessment of prognostic models

Following analysis and PH assumption testing, two independent prognostic factors (cancer status and risk score) were identified (Figure 5A-5C). A nomogram based on these prognostic factors was constructed to predict OS (Figure 5D). The calibration curves indicated good predictive agreement (Figure 5E). ROC analysis further supported the predictive performance of the nomogram, with all AUC values ≥0.85 (Figure 5F).

Figure 5 Independent prognostic analysis and nomogram construction for OS prediction. (A,B) Forest plots showing univariate and multivariate Cox regression analyses of clinical characteristics and risk scores. (C) PH assumption testing results for significant prognostic factors. (D) Nomogram integrating cancer status and risk score for predicting 1-, 3-, and 5-year OS probabilities. *, P<0.05; **, P<0.01. (E) Calibration curves demonstrating good agreement between nomogram-predicted and observed survival rates. (F) Time-dependent ROC curves validating the effectiveness of the nomogram model with all AUC values ≥0.85. AUC, area under the curve; M, metastasis; N, node; OS, overall survival; PH, proportional hazards; ROC, receiver operating characteristic; T, tumor.

Screening and functional enrichment of RRDEGs

A total of 563 RRDEGs were identified between the two risk groups, including 504 upregulated and 59 downregulated genes (Figure 6A,6B, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-6.xls). GO and KEGG enrichment analyses indicated that these RRDEGs were primarily associated with cell adhesion and angiogenesis-related processes (Figure 6C, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-7.xls). KEGG pathway analysis further revealed significant enrichment in cell adhesion molecules, complement pathways, and coagulation cascades (Figure 6D, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-7.xls).

Figure 6 Screening and functional enrichment analysis of RRDEGs. (A,B) Volcano plot and heatmap showing 563 RRDEGs between high-risk and low-risk groups, comprising 504 upregulated and 59 downregulated genes. (C) GO enrichment analysis revealing primary association of RRDEGs with cell adhesion and angiogenesis processes across BP, CC, and MF categories. (D) KEGG pathway enrichment analysis indicating significant involvement of cell adhesion molecules, complement pathways, and coagulation cascades. BP, biological process; CC, cellular component; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; RRDEGs, risk-related differentially expressed genes.

Immune-related analysis between two risk subgroups

The differential analysis identified six immune checkpoints (LAG3, HAVCR2, TIGIT, PDCD1LG2, CTLA4, and PDCD1) that were differentially expressed between the two risk subgroups (Figure 7A-7F). In addition, a significant positive correlation was observed between risk scores and five of these immune checkpoints (HAVCR2, LAG3, PDCD1, PDCD1LG2, and TIGIT) (Figure 7G). SubMap analysis further indicated that the high-risk group exhibited increased sensitivity to CTLA4 inhibitors (Figure 7H). However, CTLA4 blockade is not a standard therapeutic option for rectal cancer, and this computational prediction lacks clinical validation; therefore, its clinical relevance remains uncertain and should be interpreted as exploratory.

Figure 7 Immune checkpoint expression and immunotherapy sensitivity analysis between risk subgroups. (A-F) Box plots showing differential expression of six immune checkpoints (LAG3, HAVCR2, TIGIT, PDCD1LG2, CTLA4, and PDCD1) between high- and low-risk groups. (G) Correlation heatmap demonstrating significant positive correlations between risk scores and five differential immune checkpoints (HAVCR2, LAG3, PDCD1, PDCD1LG2, and TIGIT). *, P<0.05; **, P<0.01; ***, P<0.001. (H) Submap analysis revealing increased sensitivity to CTLA4 inhibitors in the high-risk group.

Single-gene GSEA of biomarkers

Single-gene GSEA was performed to investigate the regulatory pathways and molecular functions associated with each biomarker. CCL24 and ODC1 showed significant enrichment in GO-BP terms related to ribonucleoprotein complex biogenesis and ribosome biogenesis (Figure 8A,8B). IGFBP3, PYGM, and VKORC1 were enriched in multiple BP, including mitochondrial respiratory chain complex assembly (Figure 8C-8E, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-8.xls). Notably, ODC1 and VKORC1 were predominantly enriched in KEGG pathways such as ribosome-related processes (Figure 9A,9B), whereas CCL24, IGFBP3, and PYGM were mainly enriched in cell adhesion molecule pathways (Figure 9C-9E, table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-9.xls).

Figure 8 GSEA of GO biological processes for individual biomarkers. (A,B) GSEA results showing CCL24 and ODC1 enrichment in ribonucleoprotein complex biogenesis and ribosome biogenesis pathways. (C-E) GSEA results demonstrating IGFBP3, PYGM, and VKORC1 enrichment in various processes, including mitochondrial respiratory chain complex assembly and related biological processes. GO, Gene Ontology; GSEA, gene set enrichment analysis; ncRNA, non-coding RNA.

Figure 9 GSEA of KEGG pathways for individual biomarkers. (A,B) GSEA results showing ODC1 and VKORC1 prominent enrichment in ribosome-related pathways. (C-E) GSEA results demonstrating CCL24, IGFBP3, and PYGM primary enrichment in cell adhesion molecules and related KEGG pathways. COVID-19, coronavirus disease 2019; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Chemotherapy drug sensitivity analysis

Drug sensitivity analysis revealed significant differences in susceptibility to 51 drugs (AP.24534, AZD7762, BIBW2992, CGP.60474, VX.680, etc.) between the two risk groups (table available at https://cdn.amegroups.cn/static/public/jgo-2025-1-1061-10.xls). Among these, temsirolimus, WO2009093972, CCT007093, MK.2206, and DMOG showed the lowest P values (Figure 10). Notably, the low-risk group exhibited lower IC₅₀ values for these five drugs compared with the high-risk group, indicating greater sensitivity to these therapeutic agents. Figure S1 provides additional results.

Figure 10 Chemotherapy drug sensitivity analysis between high-risk and low-risk groups. (A-E) Box plots showing IC₅₀ values for five drugs with the lowest P values: Temsirolimus, WO2009093972, CCT007093, MK.2206, and DMOG. The low-risk group demonstrated significantly lower IC50 values compared to the high-risk group, indicating higher sensitivity to these chemotherapy agents. IC₅₀, half-maximal inhibitory concentration.

Prognostic and immunological significance of the PMRGs pathway activity score in rectal cancer

Correlation analysis demonstrated no significant association between the PMRG pathway activity score and the risk score (Figure 11A). In contrast, the PMRG pathway activity score showed significant positive correlations with most immune checkpoints, with the strongest correlation observed for PDCD1LG2 (cor =0.41, P<0.001) (Figure 11B). Furthermore, 17 immune cell types were differentially infiltrated between the high- and low-risk groups, including activated B cells, central memory CD4 T cells, and activated dendritic cells (Figure 11C). Among these, 16 immune cell types showed significant positive correlations with the PMRG pathway activity score, with type 1 helper cells exhibiting the strongest correlation (Figure 11D). Based on the median PMRG pathway activity score, TCGA-READ patients were stratified into a high-score group (37 samples) and a low-score group (36 samples). Kaplan-Meier survival analysis indicated that OS was lower in the high-score group compared with the low-score group (Figure 11E). These findings suggest that increased activity of the propionate metabolism pathway may contribute to an immunosuppressive tumor microenvironment and facilitate immune evasion in rectal cancer, thereby promoting disease progression. Targeting this pathway may represent a potential immunotherapeutic strategy for patients with high PMRG activity.

Figure 11 Prognostic and immunological significance of the PMRGs pathway activity score in rectal cancer. (A) Analysis of the correlation between the risk score and the PMRGs pathway activity score. (B) Correlation analysis between PMRGs pathway activity score and immune checkpoints. (C) Differences in immune cell infiltration between different risk groups. (D) Correlation analysis between PMRGs pathway activity score and differentially infiltrated immune cells. (E) Association between PMRGs pathway activity score and clinical outcomes. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. PMRGs, propionate metabolism-related genes.

Discussion

Rectal cancer represents a major malignancy of the gastrointestinal tract with substantial morbidity and mortality worldwide. The disease is commonly characterized by aggressive local invasion and early lymph node metastasis, both of which significantly impact patient prognosis (20,21). An increasing number of studies have demonstrated the critical role of propionate metabolism in regulating the tumor microenvironment and cellular processes in rectal cancer, potentially influencing tumor progression and metastasis (22,23). Elucidation of the relationship between propionate metabolism and rectal cancer may provide novel insights into tumor biology and facilitate the identification of innovative prognostic biomarkers.

Through DEG screening and functional enrichment analysis, candidate genes were found to be predominantly enriched in BP, including inflammatory response and lipopolysaccharide response, as well as pathways such as TNF signaling and cytokine-cytokine receptor interaction. Inflammatory responses and the TNF signaling pathway play a critical role in regulating the tumor microenvironment. They facilitate the initiation and progression of rectal cancer by promoting inflammatory responses, accelerating cell proliferation, and mediating immune evasion in rectal cancer (24-26). Cytokine-cytokine receptor interactions drive malignant tumor progression by modulating immune cell functions within the tumor microenvironment and amplifying pro-inflammatory signals (27). Notably, TNF/nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) axis-mediated inflammatory signaling promotes epithelial-mesenchymal transition (EMT) and lymphangiogenesis in tumor cells, thereby inducing lymph node metastasis in rectal cancer (28-30). Meanwhile, as a pivotal regulator of short-chain fatty acids (SCFAs) in the intestinal microenvironment, propionate metabolism modulates local inflammatory activity and interacts with inflammatory signaling pathways to jointly remodel the tumor metastatic microenvironment (31-33). Accordingly, the functional characteristics of these candidate genes not only confirm their involvement in metastasis-related pathways but also uncover mechanistic crosstalk between propionate metabolism and tumor progression, providing valuable insights into the molecular mechanisms linking propionate metabolism to rectal cancer metastasis.

Five biomarkers associated with both propionate metabolism and lymph node metastasis were identified: CCL24, IGFBP3, ODC1, PYGM, and VKORC1. These biomarkers participate in distinct BP and signaling pathways and may influence rectal cancer initiation and progression through heterogeneous molecular mechanisms. CCL24, a chemotactic factor, promotes tumor progression and lymphatic metastasis by recruiting immune cells, including eosinophils, monocytes, and immunosuppressive M2-type macrophages, into the tumor microenvironment. This enhanced immune cell infiltration facilitates tumor cell migration and invasion (34). IGFBP3, a member of the insulin-like growth factor (IGF) binding protein family, is implicated in rectal cancer development and progression through involvement in the IGF signaling pathway, wherein interactions between IGFBP3 and IGF regulate cell proliferation and apoptosis (35). ODC1 is a key enzyme in the polyamine biosynthesis pathway and plays a pivotal role in cellular growth and tumor progression. Polyamines, including spermidine and spermine, exert important regulatory effects on these processes. In high-risk populations, elevated ODC1 expression is associated with enhanced polyamine synthesis, potentially promoting tumor cell proliferation and progression (36). PYGM is an isoenzyme of glycogen phosphorylase. Although its specific role in rectal cancer remains incompletely defined, available evidence suggests an association with patient prognosis, potentially contributing to disease progression through regulation of tumor cell metabolism. In high-risk populations, PYGM may facilitate tumor cell survival and proliferation by promoting glycogenolysis (37). VKORC1 encodes vitamin K epoxide reductase complex 1 and primarily participates in the vitamin K metabolic pathway. Although its direct association with rectal cancer has not been extensively investigated, it has been implicated in anticoagulant-related processes and may indirectly influence tumor angiogenesis and growth. In high-risk groups, VKORC1 may promote tumor growth and metastasis by modulating vitamin K metabolism and related coagulation pathways, while also affecting vascular supply (38).

In this study, six differential immune checkpoints were identified, namely LAG3, HAVCR2, TIGIT, PDCD1LG2, CTLA4, and PDCD1. As classical inhibitory immune checkpoints, LAG3, HAVCR2, TIGIT, and PDCD1 (PD-1) have been widely reported to enhance tumor cell survival by suppressing T cell function and facilitating immune evasion across multiple cancer types. Results demonstrated that patients in the high-risk group exhibited elevated CTLA4 expression, which showed a positive correlation with risk scores. The high-risk group also presented increased expression of additional co-inhibitory molecules, including LAG3 and TIGIT, indicating a more pronounced immunosuppressive state and T cell exhaustion phenotype. Accordingly, a mechanistic hypothesis is proposed. Upregulated CTLA4 has been confirmed to inhibit T cell activation by competing with CD28 for binding to CD80/CD86 on antigen-presenting cells, thereby promoting tumor immune evasion (39,40). In addition, prognostic biomarkers identified in this study, such as CCL24, may further exacerbate T cell exhaustion and upregulate CTLA4 expression through recruitment of immunosuppressive M2 macrophages. This process may amplify the anti-tumor effects of restored T cell function following CTLA4 blockade, thereby rendering high-risk patients more susceptible to CTLA4 blockade therapy. This hypothesis is consistent with observations in other solid tumors, where higher CTLA4 expression has been associated with improved therapeutic responses to ipilimumab (40). Although direct clinical evidence in rectal cancer remains limited and CTLA4 blockade is not a standard therapeutic option, these findings provide bioinformatics-based theoretical support for further investigation of CTLA4 blockade in molecularly stratified high-risk patients. The proposed hypothesis requires further experimental validation in rectal cancer models and subsequent clinical confirmation.

Drug susceptibility analysis revealed significant differences in sensitivity to multiple therapeutic agents across risk groups, providing potential insights for individualized treatment strategies. Temsirolimus, a well-established mammalian target of rapamycin (mTOR) inhibitor with documented anti-tumor activity, showed increased efficacy in the low-risk subgroup, suggesting that mTOR pathway dysregulation may represent a dominant driver in these tumors (41). WO2009093972 and CCT007093, functioning as pan-CDK inhibitors and Chk1 inhibitors, respectively, imply that cell cycle dysregulation constitutes a key mechanism underlying tumor proliferation in the low-risk group (42). Increased sensitivity to these cell cycle inhibitors indicates distinct proliferative dependencies in low-risk tumors, potentially associated with altered cell cycle checkpoint regulation. Overall, the drug susceptibility analysis not only reveals differential therapeutic responses across risk groups but also further supports the feasibility of risk score-based individualized treatment strategies. These agents may provide more effective therapeutic options, particularly for patients in the low-risk group.

In this study, five prognostic biomarkers associated with both lymph node metastasis and propionate metabolism were systematically identified through comprehensive bioinformatics analyses. Based on these biomarkers, a robust prognostic model was constructed. The model effectively stratified patients according to risk levels and demonstrated a significant association between the risk score and the tumor immune microenvironment. Further analyses suggested that these biomarkers may play important roles in tumor initiation, progression, and immune evasion, particularly in high-risk patients characterized by elevated immunosuppression.

However, certain limitations should be acknowledged. Constrained by time and experimental conditions, this study primarily relied on bioinformatics analyses and prognostic model construction based on publicly available datasets, without experimental validation of the functional roles of the five prognostic genes in rectal cancer cell proliferation, migration, and chemotherapy sensitivity through in vitro experiments. Future studies will employ rectal cancer cell lines (e.g., HCT116, SW480) to systematically validate the functions of these genes and their association with propionate metabolism through Cell Counting Kit-8 (CCK-8), Transwell, and IC₅₀ assays. Continued investigation will focus on elucidating the roles of these genes in rectal cancer to support the development of more personalized and precise therapeutic strategies.

Conclusions

Five biomarkers (CCL24, IGFBP3, ODC1, PYGM, and VKORC1) were identified as being associated with both lymph node metastasis and propionate metabolic pathways.

Acknowledgments

We would like to express our sincere gratitude to all individuals and organizations that supported and assisted this research. Special thanks are given to Hao-Tang Wei and Ting-Yu Mou. Appreciation is also extended to all contributors whose support made this work possible.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1061/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1061/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-1-1061/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Huo Z, Liu G, Li J. Recent research progress and clinical status of immunotherapy for colorectal cancer. J Adv Res 2026;82:729-50. [Crossref] [PubMed]
2. Chen P, Gao YH, Xie SX, et al. Lactylation in colorectal cancer: From mechanism to clinical application. Biochem Pharmacol 2026;244:117621. [Crossref] [PubMed]
3. Maria-Alexia P, Cristina R, Andreea-Ramona T, et al. Familial colorectal cancer: risk factors, screening strategies and personalized medicine. Cancer Genet 2026;302-303:166-75. [Crossref] [PubMed]
4. Yang L, Chen N, Zhao F, et al. Local control treatment for pulmonary metastases of colorectal cancer. Crit Rev Oncol Hematol 2026;217:105062. [Crossref] [PubMed]
5. Zhang Y, Chen J, She Y, et al. Paneth-like transition drives resistance to dual targeting of KRAS and EGFR in colorectal cancer. Cancer Cell 2026;44:77-93.e8. [Crossref] [PubMed]
6. Weiss L, Stintzing S, Modest DP, et al. First-line treatment efficacy of anti-EGFR versus anti-VEGF antibodies in BRAF(V600E)-mutated metastatic colorectal cancer according to primary tumor sidedness: A pooled analysis of seven clinical trials performed in the first-line treatment of mCRC (German AIO Study Group). Eur J Cancer 2026;232:116105. [Crossref] [PubMed]
7. Cataldi C, Karaoğlan BB, Liotta E, et al. Decoding Immunotherapy Response in Colorectal Cancer: Translational Insights Beyond MSI. Cancers (Basel) 2026;18:852. [Crossref] [PubMed]
8. Nshanian M, Gruber JJ, Geller BS, et al. Short-chain fatty acid metabolites propionate and butyrate are unique epigenetic regulatory elements linking diet, metabolism and gene expression. Nat Metab 2025;7:196-211. [Crossref] [PubMed]
9. Zhao S, Zheng X, Yang C, et al. Gut microbiota causally affects ulcerative colitis by potential mediation of plasma metabolites: A Mendelian randomization study. Medicine (Baltimore) 2025;104:e42791. [Crossref] [PubMed]
10. Wang J, Tan Q, Ni M, et al. High-Fat Diet-Induced Gut Microbiota Disruption Promotes Colorectal Cancer Lymphatic Metastasis via Propionate/GPR41 Signaling. Digestion 2025;106:573-90. [Crossref] [PubMed]
11. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010;26:139-40. [Crossref] [PubMed]
12. Ito K, Murphy D. Application of ggplot2 to Pharmacometric Graphics. CPT Pharmacometrics Syst Pharmacol 2013;2:e79. [Crossref] [PubMed]
13. Ye Y, Dai Q, Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov 2021;7:71. [Crossref] [PubMed]
14. Yang CC, Tsai CG, Yeh JA. Miniaturization of dielectric liquid microlens in package. Biomicrofluidics 2010;4:43006. [Crossref] [PubMed]
15. In J, Lee DK. Survival analysis: part II - applied clinical data analysis. Korean J Anesthesiol 2019;72:441-57. [Crossref] [PubMed]
16. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
17. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
18. Schober P, Boer C, Schwarte LA. Correlation Coefficients: Appropriate Use and Interpretation. Anesth Analg 2018;126:1763-8. [Crossref] [PubMed]
19. Liu TT, Li R, Huo C, et al. Identification of CDK2-Related Immune Forecast Model and ceRNA in Lung Adenocarcinoma, a Pan-Cancer Analysis. Front Cell Dev Biol 2021;9:682002. [Crossref] [PubMed]
20. Yeh C, Wheless MC, Ciombor KK, et al. New Frontiers in the Treatment of Rectal Cancer. Annu Rev Med 2026;77:177-92. [Crossref] [PubMed]
21. Choi YH, Huh JW, Lee WY, et al. Risk factors for lymph node metastasis in p/ypT1-2 rectal cancer. Am J Surg 2026;253:116782. [Crossref] [PubMed]
22. Zhang C, Shen W, Leng Z, et al. A novel radio-immunotherapy strategy: Gut microbiota metabolite combined with radioactive hydrogel for the treatment of low rectal cancer. Biomaterials 2025;322:123386. [Crossref] [PubMed]
23. Wang X, Hu Y. High‑dose sodium propionate contributes to tumor immune escape through the IGF2BP3/PD‑L1 axis in colorectal cancer. Oncol Lett 2025;29:303. [Crossref] [PubMed]
24. Arroyo Almenas F, Törő G, Szaniszlo P, et al. Cystathionine Gamma-Lyase Regulates TNF-α-Mediated Injury Response in Human Colonic Epithelial Cells and Colonoids. Antioxidants (Basel) 2024;13:1067. [Crossref] [PubMed]
25. Wang Q, Zhong W, Xiao Y, et al. Pan-immune-inflammation value predicts immunotherapy response and reflects local antitumor immune response in rectal cancer. Cancer Sci 2025;116:350-66. [Crossref] [PubMed]
26. Ma Q, Zheng L, Cheng H, et al. PDCD4-induced oxidative stress through FGR/NF-κB axis in rectal cancer radiotherapy-induced AKI. Int Immunopharmacol 2024;132:111779. [Crossref] [PubMed]
27. De Simone V, Franzè E, Ronchetti G, et al. Th17-type cytokines, IL-6 and TNF-α synergistically activate STAT3 and NF-kB to promote colorectal cancer cell growth. Oncogene 2015;34:3493-503. [Crossref] [PubMed]
28. de Carvalho TG, Lara P, Jorquera-Cordero C, et al. Inhibition of murine colorectal cancer metastasis by targeting M2-TAM through STAT3/NF-kB/AKT signaling using macrophage 1-derived extracellular vesicles loaded with oxaliplatin, retinoic acid, and Libidibia ferrea. Biomed Pharmacother 2023;168:115663. [Crossref] [PubMed]
29. Long X, Hu Y, Duan S, et al. MRGBP promotes colorectal cancer metastasis via DKK1/Wnt/β-catenin and NF-kB/p65 pathways mediated EMT. Exp Cell Res 2022;421:113375. [Crossref] [PubMed]
30. Li Y, He J, Zhong D, et al. High-mobility group box 1 protein activating nuclear factor-κB to upregulate vascular endothelial growth factor C is involved in lymphangiogenesis and lymphatic node metastasis in colon cancer. J Int Med Res 2015;43:494-505. [Crossref] [PubMed]
31. Hou H, Chen D, Zhang K, et al. Gut microbiota-derived short-chain fatty acids and colorectal cancer: Ready for clinical translation? Cancer Lett 2022;526:225-35. [Crossref] [PubMed]
32. Dong Y, Zhang K, Wei J, et al. Gut microbiota-derived short-chain fatty acids regulate gastrointestinal tumor immunity: a novel therapeutic strategy? Front Immunol 2023;14:1158200. [Crossref] [PubMed]
33. Du Y, He C, An Y, et al. The Role of Short Chain Fatty Acids in Inflammation and Body Health. Int J Mol Sci 2024;25:7379. [Crossref] [PubMed]
34. Cheadle EJ, Riyad K, Subar D, et al. Eotaxin-2 and colorectal cancer: a potential target for immune therapy. Clin Cancer Res 2007;13:5719-28. [Crossref] [PubMed]
35. Murphy N, Carreras-Torres R, Song M, et al. Circulating Levels of Insulin-like Growth Factor 1 and Insulin-like Growth Factor Binding Protein 3 Associate With Risk of Colorectal Cancer Based on Serologic and Mendelian Randomization Analyses. Gastroenterology 2020;158:1300-1312.e20. [Crossref] [PubMed]
36. Liu N, Zhang T, Steer CJ, et al. MicroRNA-378a-3p prevents initiation and growth of colorectal cancer by fine tuning polyamine synthesis. Cell Biosci 2022;12:192. [Crossref] [PubMed]
37. Newgard CB, Littman DR, van Genderen C, et al. Human brain glycogen phosphorylase. Cloning, sequence analysis, chromosomal mapping, tissue expression, and comparison with the human liver and muscle isozymes. J Biol Chem 1988;263:3850-7.
38. Yang X, Wang Z, Zandkarimi F, et al. Regulation of VKORC1L1 is critical for p53-mediated tumor suppression through vitamin K metabolism. Cell Metab 2023;35:1474-1490.e8. [Crossref] [PubMed]
39. Lin S, Cheng L, Ye W, et al. Chimeric CTLA4-CD28-CD3z T Cells Potentiate Antitumor Activity Against CD80/CD86-Positive B Cell Malignancies. Front Immunol 2021;12:642528. [Crossref] [PubMed]
40. Xu X, Dennett P, Zhang J, et al. CTLA4 depletes T cell endogenous and trogocytosed B7 ligands via cis-endocytosis. J Exp Med 2023;220:e20221391. [Crossref] [PubMed]
41. Faivre S, Kroemer G, Raymond E. Current development of mTOR inhibitors as anticancer agents. Nat Rev Drug Discov 2006;5:671-88. [Crossref] [PubMed]
42. Thoma OM, Neurath MF, Waldner MJ. Cyclin-Dependent Kinase Inhibitors and Their Therapeutic Potential in Colorectal Cancer Treatment. Front Pharmacol 2021;12:757120. [Crossref] [PubMed]