Mito-fission gene prognostic model for colorectal cancer
Chao Liu, Sheng Xu, Yuanyuan Liu, Zhixing Lu, Jianrong Yang
TL;DR
This study identifies key genes related to mitochondrial fission in colorectal cancer and builds a risk model to predict patient survival.
Contribution
A novel risk model based on mitochondrial fission-related genes is developed for predicting colorectal cancer outcomes.
Findings
CCDC68, FAM151A, and MC1R were identified as hub genes in colorectal cancer.
The risk model showed higher risk scores correlated with worse patient survival.
RT-qPCR confirmed gene expression differences between cancer and normal tissues.
Abstract
Dysregulated cellular metabolism is one of the major causes of colorectal cancer (CRC), including mitochondrial fission. Therefore, this study focuses on the specific regulatory mechanisms of mitochondrial dysfunction on CRC, which will provide theoretical guidance for CRC in the future. The Cancer Genome Atlas (TCGA)-CRC dataset, GSE103479 dataset and 40 mitochondrial fission-related genes (MFRGs) were downloaded in this study. The differentially expressed genes (DEGs) were analyzed in TCGA-CRC samples. Using MFRGs scores as traits, key module genes associated with its scores were screened by weighted gene co-expression network analysis (WGCNA). Then, differentially expressed MFRGs (DE-MFRGs) were obtained by intersecting DEGs and key module genes. Next, DE-MFRGs were subjected to univariate Cox, least absolute shrinkage and selection operator (LASSO), multivariate Cox and stepwise…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13- —The Autonomous Region Health Commission Self-funded Research Project
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsGenetic factors in colorectal cancer · Ferroptosis and cancer prognosis · Cancer, Hypoxia, and Metabolism
Introduction
Colorectal cancer (CRC) is one of the most prevalent malignancies of the digestive system and ranks as the second deadliest cancer globally. It is classified into Colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) based on the tumor’s site of origin. Due to their shared etiology, pathogenesis, and histological features, COAD and READ are frequently grouped together under the umbrella term CRC (Islam et al., 2022; Sedlak, Yilmaz & Roper, 2023). The treatment strategy for CRC primarily relies on surgery, supplemented by chemotherapy, radiotherapy, immunotherapy, and targeted therapies (Ghazi et al., 2022). While advances in colectomy, chemotherapy, and immunotherapy have significantly improved the 5-year survival rate for patients with CRC, outcomes for those with advanced CRC remain poor (Yoshino et al., 2022). Consequently, there is an urgent need to unravel the molecular mechanisms underlying CRC and identify novel therapeutic targets to improve patient prognosis.
Mitochondria, integral to multicellular life, generate adenosine triphosphate (ATP) through oxidative phosphorylation, providing the energy required for cellular metabolism. Due to this role, they are often referred to as the “energy factories” of the cell (Harrington et al., 2023). Mitochondria exhibit dynamic structural changes, adjusting their fusion and fission in response to cellular metabolic demands. They also regulate their form, size, quantity, distribution, quality control, and transport within the cell to maintain energy balance (Quiles & Gustafsson, 2022). Disruption of mitochondrial fission leads to alterations in cellular metabolism, proliferation, and apoptosis (Kleele et al., 2021). Recent studies have highlighted that increased mitochondrial fission (MF) can promote tumor growth and metastasis (Colpman, Dasgupta & Archer, 2023). It has been demonstrated that enhanced MF facilitates the metabolic shift from glycolysis to oxidative phosphorylation, thereby supporting tumor cell survival under energy stress, indicating MF’s pivotal role in regulating tumor cell metabolism (Plecitá-Hlavatá et al., 2008; Westermann, 2012). Additionally, mitochondrial fission influences the onset and progression of pancreatic ductal adenocarcinoma, Lung adenocarcinoma, and Hepatocellular carcinoma (Rehman et al., 2012; Kashatus et al., 2015; Zhang et al., 2020). In CRC, mitochondria regulate tumorigenesis through mechanisms such as histone acetylation (Ohshima et al., 2022), and carposide II has been shown to inhibit CRC development by modulating MF and NF-κB pathways (Chen et al., 2019). Despite these findings, the precise mechanism by which mitochondrial fission influences CRC remains unclear. Therefore, this study aims to further explore the underlying mechanisms of mitochondrial fission in CRC.
Utilizing public databases, this study systematically establishes the relationship between mitochondrial fission and CRC by constructing a risk model, investigates its biological functions, and evaluates the prognostic significance of key genes, providing a theoretical foundation for CRC treatment.
Materials & Methods
Data extraction
RNA-seq, clinical, and survival data were retrieved from The Cancer Genome Atlas (TCGA) database, encompassing 606 CRC tissue samples (including survival data) and 51 normal tissue samples, hereafter referred to as TCGA-CRC. Additionally, the GSE103479 dataset was extracted from the Gene Expression Omnibus (GEO) database, which included 155 CRC tissue samples with associated survival data. Two gene sets, GOBP MITOCHONDRIAL FISSION and GOBP POSITIVE REGULATION OF MITOCHONDRIAL FISSION, were selected from MSigDB using “mitochondrial fission” as the keyword, resulting in the identification of 40 mitochondrial fission-related genes (MFRGs) for subsequent analysis (Table 1).
Table 1: List of 40 MFRGs.
Identification of DEGs
In TCGA-CRC samples, differentially expressed genes (DEGs) were identified using the DESeq2 package (version 1.34.0) (Love, Huber & Anders, 2014), with thresholds set at P.adj < 0.05 and log_2_FC > 1.5. The top 20 DEGs, sorted by log_2_FC, were visualized via volcano and heat maps.
Analysis of MFRGs scores
MFRG scores for each sample were calculated using the GSVA package (version 1.42.0) (Hänzelmann, Castelo & Guinney, 2013) to assess differences between sample groups. A rank-sum test was employed to evaluate differences in MFRG scores between the two sample groups. The impact of MFRGs on the survival of patients with CRC was evaluated by stratifying samples based on the optimal cut-off values of MFRG scores. Kaplan–Meier (KM) curves were plotted for high- and low-MFRG score groups using the survminer package (version 0.4.9) (Li et al., 2020). in order to understand the relationship between MFRGs scores and clinical characteristics (age (50-year cut-off), gender, ethnicity and T/N/M stage), CRC patients were divided into different clinical subgroups, and the differences in MFRGs scores between different subgroups were compared using the Wilcoxon test.
Weighted gene co-expression network analysis
The expression matrix of all CRC samples from TCGA-CRC was analyzed using weighted gene co-expression network analysis (WGCNA) (version 1.70.3) (Langfelder & Horvath, 2008). Initially, outlier samples were excluded by clustering. The optimal soft threshold (R^2^ = 0.85) was determined by evaluating the relationship between the soft threshold and the scale-free network evaluation coefficient. A phylogenetic tree of genes was then constructed based on gene similarity. Using MFRG scores as the phenotype, the correlation between the phenotype and gene modules was computed. Correlation coefficients and P-values were calculated, and a correlation heat map was generated. The module with the strongest correlation to the MFRG score (R > 0.5 and P < 0.05) was selected as the key module. Further gene filtering within the module was performed by setting the gene significance (GS) to 0.4 and module membership (MM) to 0.4, thereby identifying key module genes associated with the MFRG score.
Acquisition and enrichment analysis of DE-MFRGs
DE-MFRGs were identified by intersecting DEGs and the key module genes associated with MFRGs using the VennDiagram package (version 1.7.1) (Ito & Murphy, 2013). To explore the biological functions of the DE-MFRGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the clusterProfiler package (version 4.2.2) (Wu et al., 2021), with a significance threshold set at P < 0.05.
Risk model construction
To assess the prognostic impact of DE-MFRGs on patients with CRC in TCGA-CRC, a risk model was constructed. First, univariate Cox regression analysis was performed using the survival package (version 3.5-3) (Deng & Thompson, 2023) to calculate risk scores and identify significant genes (P-value < 0.05). In the second step, Least Absolute Shrinkage and Selection Operator (LASSO) analysis was conducted using the glmnet package (version 4.1-4) (Friedman, Hastie & Tibshirani, 2010) to select genes with strong correlations for constructing the optimal risk model. In the third step, the proportional hazards (PH) assumption was tested for the genes selected from the LASSO analysisd, and those with P > 0.05 were included in further analysis, followed by multivariate Cox regression analysis. Finally, a stepwise regression method was applied to optimize the multivariate Cox regression results by entering variables one by one and checking their significance, removing non-significant variables, thereby constructing the optimal regression equation and identifying the hub genes associated with prognosis.
Validation of risk model
Using the expression levels of hub genes and the risk coefficients obtained from stepwise regression analysis, the risk score formula ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \text{Riskscore}={\mathop{\sum }\nolimits }{\mathrm{i}=1}^{\mathrm{n}}\text{coef}(\mathrm{gen}{\mathrm{e}}{\mathrm{i}})\times \text{expr}(\mathrm{gen}{\mathrm{e}}_{\mathrm{i}})\end{document} ) was used to calculate the risk score for each patient, with patients with CRC categorized into high- and low-risk groups based on the median risk score. The prognostic value of the risk model was assessed by comparing survival statuses between the two risk groups, and KM survival curves were plotted using the survminer package (version 0.4.9) (Li et al., 2020). The model’s validity was further assessed by plotting receiver operating characteristic (ROC) curves. The risk model was then validated in the GSE103479 dataset using the same methodology.
Independent prognostic analysis and clinical characteristic correlation analysis
For univariate Cox analysis, six clinical characteristics (age, gender, race, T/N/M stages, and risk score) were included in the model. The results underwent PH assumption testing and multivariate Cox analysis to identify independent prognostic factors associated with CRC. Additionally, to evaluate the risk model’s applicability for patients with CRC, a nomogram predicting 1-, 2-, and 3-year survival was constructed using the rms package (version 6.5-1) (Sachs, 2017). The nomogram’s validity was verified by calibration curves and ROC analysis. Finally, in TCGA-CRC, differences in risk scores across the six clinical characteristics were analyzed using the rank-sum test (P < 0.05).
GSEA analysis
Differential analysis between the two risk groups in TCGA-CRC was performed using DESeq2 (version 1.34.0) (Love, Huber & Anders, 2014), with —log_2_FC— calculated and ranked. GSEA was conducted using the KEGG background dataset (c2.cp.kegg.v2023.1.Hs.symbols.gmt), with a threshold set at FDR < 0.05. The top six enriched signaling pathways were visualized using the enrichplot package (version 1.18.3) (Zhang et al., 2019).
Immune related analysis
To further investigate the differences in the immune microenvironment between the high- and low-risk groups in CRC samples, the relative abundance of 22 immune cell types across all TCGA-CRC samples was evaluated using the “CIBERSORT” algorithm, and immune cell proportions were compared between high- and low-risk group samples. Differences in immune cell composition between the risk groups were assessed with the “Wilcoxon” test (P < 0.05). The expression levels of 48 immune checkpoints were also compared between these groups, and the TIDE score, immune exclusion score, and immune dysfunction score were calculated using the TIDE database (http://tide.dfci.harvard.edu/) to evaluate immune therapy efficacy.
Drug sensitivity analysis
The half maximal inhibitory concentration (IC50) for 138 common chemotherapeutic agents was calculated for all TCGA-CRC individuals using the “pRRophetic” algorithm (version 0.5) (Geeleher, Cox & Huang, 2014). Differences in IC50 values between the two risk groups were compared using the Wilcoxon test. A Spearman correlation analysis was performed to assess the relationship between the sensitivity to chemotherapeutic agents and hub gene expression (|r| > 0.4).
Pan-cancer analysis
To examine the expression of hub genes across different cancer types, the Wilcoxon test was conducted using the rstatix package (version 0.7.2) with a significance threshold of P < 0.05. Results were visualized with the ggplot2 package (version 3.5.1). Box plots displayed differential expression of hub genes in pan-cancer samples, while violin plots illustrated the differences in hub gene expression between tumor and normal samples in TCGA-CRC.
RNA isolation and quantitative real-time polymerase chain reaction (qRT-PCR)
Twenty tissue samples, including 10 CRC tumor and 10 normal samples, were collected, frozen, and stored at −80 °C. Each 50 mg sample was lysed with TRIzol reagent for total RNA isolation, following the manufacturer’s instructions. RNA concentration was measured by NanoPhotometer N50, and its quality was assessed based on concentration, purity (A260/A280 ratio), and the amplification curve. RNA was reverse transcribed into cDNA using the SureScript First Strand cDNA Synthesis Kit (Servicebio, Wuhan, China). The qRT-PCR reactions consisted of three µL of reverse transcription product, five µL of 2xUniversal Blue SYBR Green qPCR Master Mix, and one µL each of forward and reverse primers. The amplification protocol included an initial denaturation step at 95 °C for 1 min, followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 55 °C for 20 s, and extension at 72 °C for 30 s. Three technical replicates were performed for each sample. Primer sequences are provided in Table 2. GAPDH served as the internal control, and relative gene expression was quantified using the 2^−ΔΔCT^ method. GraphPad Prism 5 was used for graph creation and statistical analysis, with p-values calculated for each comparison.
Table 2: Primer sequence information for RT-qPCR.
Ethics approval and consent to participate
This study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of The People’s Hospital of Guangxi Zhuang Autonomous Region (2024.03.01, KY-ZC-2024-030).
Consent to participate
Written informed consent was obtained from all participants.
Statistical analysis
All statistical analyses were performed using R software. Differential expression analysis was conducted using the DESeq2 package, and survival analysis was carried out using the survminer package to plot KM survival curves. WGCNA was performed using the WGCNA package, and functional enrichment analysis was conducted with the clusterProfiler package. Cox regression analysis was performed using the R survival package, while LASSO regression analysis was carried out with the glmnet package. Additionally, ROC curves were plotted using the survivalROC package, and nomograms were constructed using the rms package. Immune infiltration cell analysis was conducted using the CIBERSORT algorithm, and the IC50 values for common chemotherapy and molecular targeted drugs were calculated using the pRRophetic package. Correlation analysis was performed using Spearman’s correlation. For pairwise comparisons, the Wilcoxon test was applied, and P < 0.05 was considered statistically significant.
Results
Identification of DEGs and characterization of the correlation between MFRGs and CRC
In TCGA-CRC samples, a total of 3,310 DEGs were identified, comprising 1,669 upregulated and 1,641 downregulated genes (Figs. 1A–1B). To investigate the relationship between MFRGs and CRC, differential analysis revealed that the MFRG score was significantly lower in CRC samples compared to normal samples (Fig. 1C). KM survival analysis showed that groups with high MFRG scores exhibited improved survival probabilities, suggesting that MFRGs may influence the survival of patients with CRC (Fig. 1D). Further analysis of clinical characteristics demonstrated significant differences in MFRG scores across various factors, including race, M stage (M0 vs M1), and N stage (N0 vs N1, N0 vs N2) (Fig. 1E).
The expression of DEGs between colorectal cancer group and normal group in TCGACRC samples.(A) Volcanic map of differentially expressed genes distribution between CRC and Normal groups: red dots are up-regulated genes, green dots are down-regulated genes, and gray dots are undifferentiated genes. (B) Heat map of differentially expressed genes between the CRC group and Normal group: the upper part is the heat map of expression quantity density of the Top20 differentially expressed genes in the sample, showing the lines of five quantiles and average values; the next part is the expression heat map of the top 20 differential genes in the sample. (C) MFRGs score difference among different samples in training set. (D) KM survival curve of high and low MFRGs rating groups. (E) Differences in MFRGs scores among different clinical subgroups.
Acquisition of 49 DE-MFRGs
Sample clustering and phenotypic trait heatmaps showed no outlier samples (Fig. 2A). A soft-threshold of 7, with an R^2^ value of 0.85, was selected to construct a scale-free network (Fig. 2B). Following average hierarchical clustering and dynamic tree clipping, 16 modules were identified (Fig. 2C). The correlation heatmap revealed that the MEgreen (Cor = 0.67), MEturquoise (Cor = 0.54), MEred (Cor = −0.52), and MEsalmon (Cor = −0.52) modules were significantly associated with CRC (P < 0.05) (Fig. 2D). A total of 917 genes, 3,143 genes, 907 genes, and 302 genes were derived from these modules, resulting in 1,952 hub module genes (GS = 0.4, MM = 0.4) (Figs. 2E–2H). By intersecting the 3,310 DEGs with the 1,952 hub genes, 49 DE-MFRGs were identified (Fig. 2I).
Acquisition of 49 DE-MFRGs.(A) Sample level clustering after the introduction of sample traits. (B) Soft threshold selection: the determination of the optimal soft threshold mainly refers to the figure on the left, that is, the scale-free fit index (Y-axis) under different soft thresholds (X-axis), where the red line represents the value of the selected scale-free fit index. From the figure on the left, the value when the scale-free fit index is 0.85 is the minimum soft threshold. The figure on the right shows the network connectivity under different soft thresholds. (C) Identification of co-expression modules: the top half is the hierarchical clustering tree of genes, and the bottom half is the gene module, that is, the network module. Corresponding to the top and bottom, it can be seen that the genes that are close to each other (clustering to the same branch) are divided into the same module. (D) Correlation heat map of modules and phenotypes: the leftmost color block represents the module, and the rightmost color bar represents the correlation range. In the middle part of the heat map, the darker the color, the higher the correlation, red indicates positive correlation, blue indicates negative correlation; The numbers in each cell represent relevance and significance. (E–H) Scatterplot of GS and MM in key modules: GS (Gene Significance) represents the correlation between gene expression and phenotypic data. MM (Module membership) refers to the correlation between gene expression and module. (I) Identification of differential MFRGs related genes.
Biological functions and signaling pathways involved in 49 DE-MFRGs
Functional enrichment analysis of the 49 DE-MFRGs revealed 15 significant terms in GO analysis, including caspase binding, helicase activity, kinetochore formation, and cell division site organization (Figs. 3A–3C). KEGG pathway analysis identified three key pathways, such as pentose and glucuronate interconversion and nucleotide sugar biosynthesis (P < 0.05) (Fig. 3D).
Biological functions and signaling pathways involved in 49 DE-MFRGs.(A–D) GO and KEGG enrichment results of differential MFRGS-related genes: color bands on the left represent logFC of genes, and different bands on the right represent different pathways.
Identification of hub gene and risk model construction
The prognostic value of the 49 DE-MFRGs was assessed by constructing a risk model. Initially, univariate Cox analysis was performed on the 49 DE-MFRGs in 606 CRC samples, yielding 10 genes associated with survival (P < 0.05). Of these, HIG1 hypoxia inducible domain family member 1A (HIGD1A), diaphanous related formin 3 (DIAPH3), coiled-coil domain containing 68 (CCDC68) and family with sequence similarity 151 member A (FAM151A) were identified as protective factors, while keratin associated protein 5-1 (KRTAP5.1), heat shock transcription factor 4 (HSF4), melanocortin 1 receptor (MC1R), Zinc finger protein 692 (ZNF692), Lck interacting transmembrane adaptor 1 (LIME1) and tweety family member 3 (TTYH3) were risk factors (Table 3 and Fig. 4A). LASSO analysis, using the minimum lambda value, selected seven genes (HIGD1A, DIAPH3, CCDC68, FAM151A, HSF4, MC1R and TTYH3) to construct the optimal risk model (Figs. 4B–4C). The PH assumption test for these seven genes showed P-values greater than 0.05, and a forest plot from multivariate Cox analysis was generated (Fig. 4D). Finally, stepwise regression analysis identified three hub genes for the construction of a prognostic risk model: CCDC68 (HR = 0.74, 95% CI [0.56–097], P = 0.03) and FAM151A (HR = 0.38 95% CI [0.20–0.71], P = 0.002) were protective factors, and MC1R (HR = 2.22, 95% CI [1.56–3.17], P < 0.001) was a risk factor (Table 4 and Fig. 4E).
Table 3: List of functional information for ten genes from univariate Cox analysis.
Identification of CCDC68, FAM151A and MC1R and risk model construction.(A) Unifactor Cox regression analysis of forest map. (B) Ten cross-validations of the adjusted parameters in LASSO analysis: the horizontal coordinate is the logarithm of the lambdas, and the vertical coordinate is the model error. The optimal lambda value is at the lowest point of the red curve, and the corresponding number of variables is 7. (C) LASSO coefficient spectrum: the horizontal coordinate is the logarithm of the lambdas, and the vertical coordinate is the variable coefficient. As the lambdas increase, the variable coefficient approaches 0. When the optimal lambda is reached, the variable whose culling coefficient is equal to 0. (D) Multivariate Cox regression analysis of forest map. (E) Stepwise regression analysis of forest map.
Risk model had good predictability for CRC patients
In TCGA-CRC, the high- and low-risk groups divided based on the median of the risk score were calculated, with results showing a significant increase in the number of deaths as the risk scores rose (Figs. 5A–5B). KM analysis revealed that patients in the high-risk group had a lower survival rate (Fig. 5C). ROC analysis showed area under the curve (AUC) values of 0.631, 0.626, and 0.663 at 1-, 2-, and 3-year intervals, respectively (Fig. 5D), indicating that the risk model demonstrated reasonable predictive accuracy. The risk model was then validated in the GSE103479 dataset, with results consistent with TCGA-CRC, confirming that high-risk patients with CRC had worse survival and prognosis (Figs. 5E–5H).
RiskScore, age, and N/M stages were independent prognostic factors for CRC
Further investigation of the correlation between clinicopathologic characteristics and the risk model in TCGA-CRC samples revealed that the risk score, along with age and T/N/M stages, showed significant association with the risk model in univariate Cox analysis (P < 0.05) (Fig. 6A), all of which satisfied the PH assumption (P > 0.05). Multivariate Cox analysis demonstrated that, except for T stage, the remaining clinicopathologic characteristics were significantly correlated with the risk model (P < 0.05) (Fig. 6B). A nomogram was developed to visualize the clinical characteristics, showing that higher scores correlated with increased mortality. The validity of the nomogram was confirmed by calibration curves, with an AUC exceeding 0.6, affirming its strong predictive ability for the survival of patients with CRC (Figs. 6C–6E). Further analysis of the risk scores and clinicopathologic features revealed significant differences in risk scores across T/N/M stages (Fig. 6F).
Enrichment pathways of high- and low-risk groups and their effects in the immune micro environment
GSEA identified the top six signaling pathways significantly associated with the risk model, including cardiac muscle contraction, ascorbate and aldarate metabolism, pentose and glucuronate interconversion, retinol metabolism, extracellular matrix (ECM)-receptor interaction, and dilated cardiomyopathy (Fig. 7A).
Differences in the immune microenvironment between the high- and low-risk groups in CRC were also explored. The proportional distribution of 22 immune cell types is shown in Fig. 7B. Notably, except for M0 macrophages, resting dendritic cells, memory CD4+ T cells, activated memory CD4+ T cells, and plasma cells, the low-risk group exhibited more significant differences in immune cell proportions (Fig. 7C), indicating that the immune system of low-risk group patients was more active and more favorable for controlling tumor progression. Furthermore, 16 immune checkpoints showed significant differences between the high- and low-risk groups, including HHLA2, TNFRSF4, CD160, TNFSF4, and TIMP3 (Fig. 7D). These immune checkpoints may be involved in the activation or inhibition of immune cells, suggesting differences in immune treatment responses between the risk groups. In addition, the TIDE score and immune exclusion score of the high-risk group were significantly higher than those of the low-risk group (Fig. 7E). These two scores are more meaningful in evaluating the effectiveness of immunotherapy. The high-risk group patients, due to their higher scores, were associated with poorer prognosis, including shorter survival and higher recurrence rates. The higher the TIDE score, the poorer the response to immune checkpoint inhibitors and the shorter the survival time.
Table 4: Results of stepwise regression analysis.
Risk model had good predictability for CRC patients.(A) Risk curves for the high-low risk group of CRC patients in the training set. (B) Scatterplot of the high-low risk grouping of CRC patients in the training set. (C) Survival curves of CRC high-low risk groups in the training set: red represents high risk group, blue represents low risk group. (D) ROC curves of CRC patients at 1, 2 and 3 years of training set. (E) Risk curves of CRC patients in high-low risk groups were validated. (F) Scatter plots of high-low risk groups of CRC patients were validated. (G) Survival curves of the high-low risk group of CRC patients in the validation set. (H) ROC curves of CRC patients at 1, 2, and 3 years were validated.
RiskScore, age and N/M stages were independent prognostic factors for CRC.(A) Univariate independent prognostic analysis of CRC forest map. (B) Multivariate independent prognostic analysis of CRC forest map. (C) The survival nomogram of CRC patients was constructed based on risk model and clinical features. (D) Calibration curve of clinical feature nomogram: the horizontal axis represents the probability of different clinical outcomes predicted by the model, and the vertical axis represents the probability of actually observed clinical outcomes of patients, which is represented by the form of median plus mean, and an ideal curve with slope of 1 is drawn as a reference. The closer the actual curve is to the ideal curve, the better the calibration degree is, that is, the smaller the deviation between the predicted results of the model and the actual results, the better the model effect. (E) ROC curve of clinical features. (F) Correlation analysis between risk score and clinical features.
Fifty-one common chemotherapeutic drugs were effective in high-risk patients for CRC
The IC50 values of 71 chemotherapeutic drugs showed significant differences between the two risk group samples, with 51 drugs demonstrating better responses in the high-risk group patients (IC50 values lower in the high-risk group than in the low-risk group) (Figs. 8–9), while the remaining 20 drugs performed better in the low-risk group (Fig. 10). It was worth noting that 20 of these drugs had already been clinically applied (Table 5). This suggested a potential significant difference in chemotherapy drug sensitivity between the high-risk and low-risk groups. Spearman correlation analysis revealed that FAM151A was positively correlated with AG.014699, camptothecin, NVP.BEZ235, and TW.37 (r > 0.4), while CCDC68 was negatively correlated with BMS.708163 and X681640 (r < −0.4) (Fig. 11). These findings indicate that the expression levels of certain genes may be associated with the sensitivity to specific drugs, suggesting that gene expression levels could influence tumor cell response to these drugs, or the drugs could impact gene-related functions.
Associated pathways of three hub genes and their effects in the immune micro environment.(A) KEGG enrichment signaling pathway in high-low risk groups: This diagram can be divided into three parts. Part I: The top five lines are the lines of gene Enrichment Score. The vertical axis is the corresponding Running ES, and there is a peak value in the line graph, which is the Enrichemnt score of this gene set, and the genes before the peak value are the core genes under this gene set. The horizontal axis represents each gene under this gene set, corresponding to the bar code-like vertical line in the second part. Part 2: The barcode-like part, called Hits, where each vertical line corresponds to a gene under the gene set. Part 3: Sequencing of genes. (B) Proportion of immune cells in the high-low risk group. (C) Differences in immune cells between high and low risk groups. (D) Differences in 48 immune checkpoints between high and low risk groups. (E) TIDE score difference between high and low risk groups violin chart.
Twenty-seven drugs that perform better in high-risk groups.The horizontal coordinate is the high-low expression group; the ordinate is IC50.
Twenty-four drugs that perform better in high-risk groups.The horizontal coordinate is the high-low expression group; the ordinate is IC50.
Ten drugs that perform better in high-risk groups.The horizontal coordinate is the high and low expression group; the ordinate is IC50.
Table 5: Table of relevant information on 71 discrepant drugs.
Gene correlation between chemotherapy drugs and risk model.
Three hub genes linked in other diseases
Pan-cancer analysis showed that CCDC68, FAM151A and MC1R were significantly different in thyroid cancer (THCA), renal clear cell carcinoma (RCC) and lung squamous cell carcinoma (LUSC) etc. (Fig. 12). In TCGA-CRC, the differential expression between CCDC68 and FAM151A were more significant in the normal group, whereas MC1R had opposite results to the above two genes.
Three hub genes linked in other diseases.(A) Differential analysis of CCDC68 in different samples of pancarcinoma. (B) Differential analysis of FAM151A in different samples of pancarcinoma. (C) Differential analysis of MC1R in different samples of pancarcinoma.
Validation of expression of 3 biomarkers
In TCGA-CRC, the differential expression of CCDC68 and FAM151A was more significant in the normal tissue samples, whereas MC1R exhibited opposite expression patterns (Fig. 13A). RT-qPCR validation of hub gene expression in clinical CRC and normal tissues confirmed these findings. CCDC68 and FAM151A were more highly expressed in normal samples compared to CRC samples, while MC1R expression was higher in CRC samples, consistent with the dataset results (Figs. 13B–13D). Notably, the expression of FAM151A (P < 0.0001) and MC1R (P < 0.0001) exhibited significant differences between clinical CRC and normal samples.
Validation of expression of 3 biomarkers.(A) Difference analysis of risk model genes in different samples of training set. (B) Expression of CCDC68 in normal samples and CRC samples. (C) Expression of FAM151A in normal samples and CRC samples. (D) Expression of MC1R in normal samples and CRC samples.
Discussion
The rising incidence and mortality of CRC have elevated its status as a significant public health issue (Bray et al., 2024). Recent studies have highlighted the pivotal role of mitochondrial fusion–fission dynamics in the initiation and progression of CRC (Wu et al., 2024). Dysregulated mitochondrial fission or fusion can contribute to the metabolic reprogramming of tumor cells, thereby activating oncogenic pathways that drive cell proliferation, invasion, migration, and drug resistance (Wu et al., 2024). This study, leveraging public databases, systematically establishes the association between mitochondrial fission and CRC by constructing a risk model, exploring its biological functions, and evaluating the prognostic efficacy of key genes, providing a theoretical foundation for CRC treatment.
GO and KEGG enrichment analyses were performed on DE-MFRGs. These analyses revealed that biological processes enriched by DE-MFRGs predominantly involve cell division and motility. Cytokinesis, a post-mitosis process, ensures the equal distribution of cell membrane, cytoskeleton, organelles, and soluble proteins to form two daughter cells. The study showed that the cortical protein CTTN upregulates the expression of the cytoplasmic division protein dedicator of cytokinesis 1 (DOCK1), and silencing DOCK1 impairs the migration and invasion capabilities of CTTN. Thus, CTTN promotes CRC metastasis by increasing DOCK1 expression (Jing et al., 2016). Furthermore, mutations in adenomatous polyposis coli (APC), commonly found in CRC, inhibit cell division by preventing mitotic spindle anchorage in the late cortex, thereby obstructing the initiation of the cytokinetic furrow (Caldwell, Green & Kaplan, 2007). Cellular components enriched by DE-MFRGs include chromosomal regions and centromere regions. Chromosomal rearrangements, such as deletions of chromosome arm 8p and amplifications of 8q, are prevalent in CRC (Brueckner et al., 2013). Studies on interchromosomal eight deletions in patients with CRC, using 11 microsatellite markers, revealed that eight markers are located in the centromeric region of chromosome 8p (Chughtai et al., 1999). Tumor suppressor genes may also be involved in a specific site of CRC, notably in a large centromere region between D11S897 and D11S925 (Connolly et al., 1999). These findings underscore the critical role of chromosomal and centromeric regions in CRC initiation and progression, aligning with the results of this study. Enriched molecular functions of DE-MFRGs include ATP-dependent chromatin remodeling enzyme activity, caspase binding, and DNA helicase activity. ATP-dependent chromatin remodeling enzymes, which alter nucleosome structure (Sundaramoorthy & Owen-Hughes, 2020), facilitate the accessibility of DNA sequences to interacting proteins, thereby enabling precise regulation of eukaryotic gene expression (Sundaramoorthy, 2019). Chromatin remodeling is essential for efficient DNA repair, genome stability, and therapeutic responses. Mutations or overexpression of Helicase, lymphoid specific (HELLS), a member of the ATP-dependent chromatin remodeling SNF2 family, have been linked to various cancers, including CRC, hepatocellular carcinoma, and leukemia (Peixoto et al., 2022). Caspases, an evolutionarily conserved cysteine protease family, play pivotal roles in cell death and inflammation (Van Opdenbosch & Lamkanfi, 2019), with caspase activation marking the irreversible point of cell death (Boatright & Salvesen, 2003). The combined expression of Caspase-8 and Caspase-3 exhibits synergistic effects and serves as an effective prognostic indicator for patients with CRC (Jing et al., 2016). Additionally, p20BAP31 has been shown to induce CRC cell apoptosis via the AIF Caspase-independent and ROS/JNK mitochondrial pathways (Jiang et al., 2023). DNA helicases, which catalyze the unwinding of duplex nucleic acids using ATP hydrolysis, are crucial for various DNA-related biological functions (Hidese et al., 2018). These helicases not only maintain genome stability but also play significant roles in cancer. Their involvement in DNA damage, replication stress responses, and repair pathways underscores their critical function in cancer biology (Dhar, Datta & Brosh Jr, 2020). These findings provide valuable insights into the mechanistic understanding of CRC onset and offer potential avenues for its prevention and treatment.
Three key nuclear genes—CCDC68, FAM151A, and MC1R—were identified through the construction of a risk model, with their prognostic significance assessed and validated. CCDC68 may contribute to CRC pathogenesis via multiple signaling pathways, and functional experiments demonstrate its inhibitory effects on CRC cell growth in vitro and tumor formation in vivo (Wang et al., 2021). While CCDC68 acts as a pro-carcinogenic factor in IL-6-stimulated endometrial cancer cells (Li et al., 2021), it also exhibits tumor-suppressor properties in pancreatic ductal adenocarcinoma (Radulovich et al., 2015). Although CCDC68 has not been directly implicated in mitochondrial division, its potential role in cytoskeleton, signaling, and gene expression regulation may indirectly affect mitochondrial morphology and function (Huang et al., 2017). FAM151, a member of the PLC-like phosphodiesterase superfamily, remains enigmatic, with its substrate and function yet to be identified (Findlay et al., 2020). No research has been conducted on the impact of FAM151A in CRC; however, Fam151B homozygous knockout mice develop retinal degeneration, with signs of retinal stress and rapid loss of photoreceptor cells in the eye, while FAM151A homozygous mutant mice have no discernable phenotype, suggesting that Fam151b and FAM151A may be functionally different (Findlay et al., 2020). Although no direct evidence links FAM151A to mitochondrial fission, it may influence mitochondrial dynamics through signal transduction, cellular structural interactions, and gene expression regulation, warranting further investigation. In addition, FAM210B participates in the mitochondrial energy metabolism of erythroblasts and makes a prominent contribution to erythrocyte differentiation by regulating mitochondrial energy metabolism (Suzuki et al., 2022). MC1R, a member of the G protein-coupled receptor (GPCR) subfamily, regulates key physiological and behavioral features via melanocortin binding, making it a potential target for melanoma therapy (Guida, Guida & Goding, 2022). High MC1R expression is significantly associated with microsatellite instability (MSI) (Peng et al., 2021). MC1R signaling accelerates G1/S phase progression and promotes breast cancer progression through the cAMP-CREB/ATF-1 and ERK-NFκB pathways (Chelakkot et al., 2023). Polymorphisms in the MC1R gene, leading to red pigmentation, are linked to a reduced risk of prostate cancer (Weinstein, Virtamo & Albanes, 2013), and MC1R is highly expressed in esophageal squamous cell carcinoma (Zhou et al., 2022). Furthermore, MC1R plays a pivotal role in CRC progression and may serve as a marker of poor prognosis in CRC (Peng et al., 2021). Existing research suggests that MC1R may influence mitochondrial fission through two primary mechanisms. The first involves oxidative stress regulation, as MC1R plays a role in modulating oxidative stress responses in melanocytes (Lu et al., 2024). Mitochondria serve as the primary source of intracellular reactive oxygen species (ROS), and the balance between mitochondrial fission and fusion is essential for maintaining mitochondrial function and cellular redox homeostasis (Westermann, 2012). MC1R activity may influence mitochondrial function, thereby indirectly impacting mitochondrial fission. The second mechanism is mediated through cellular signaling pathways. MC1R activation by α-MSH directly triggers cAMP signaling, leading to AMPK activation (Sun et al., 2023). Under energy stress, AMPK localizes to mitochondria, where it phosphorylates Ser637 of mitochondrial dynamic associated protein 1 (DRP1) and Ser155/Ser172 of its receptor, mitochondrial fission factor (MFF), facilitating DRP1 recruitment to the outer mitochondrial membrane and inducing fission (Hsu et al., 2022). Additionally, MC1R deficiency has been linked to metabolic dysregulation. MC1R-knockout mice exhibit significant hepatomegaly, accompanied by elevated hepatic and plasma cholesterol and triglyceride levels, suggesting that hepatocyte MC1R signaling regulates cholesterol and bile acid metabolism, while its absence promotes hypercholesterolemia (Thapa et al., 2023). Furthermore, some downregulated genes following MC1R knockout have been implicated in cell migration and melanoma metastasis (Seong & Kim, 2014), but its role in CRC cell migration and metastasis remains uncharacterized and warrants further investigation. This study reveals a strong association between CCDC68, FAM151A, and MC1R in CRC development, suggesting a potential therapeutic strategy. The study model demonstrates a correlation with N stage and M stage in CRC, and CCDC68 expression emerges as a promising molecular marker for prognostic evaluation in patients with CRC.
Cancer development and progression are closely linked to alterations in the tumor microenvironment. Cancer cells actively reshape their surroundings by secreting various cytokines, chemokines, and other factors, which reprogram the neighboring cells to support tumor survival and progression (De Visser & Joyce, 2023). Immune cells, key components of the tumor stroma, play a critical role in this dynamic process (Hinshaw & Shevde, 2019). This study assessed the relative proportions of 22 immune cell types in high- and low-risk CRC subgroups, evaluating both immune cell infiltration and immune function scores. To investigate differences in immune infiltration between these subgroups, Wilcoxon tests were performed. The results revealed five immune cell types with significant differences between the high- and low-risk groups: resting dendritic cells, resting memory CD4+ T cells, activated memory CD4+ T cells, plasma cells, and M0 macrophages. Dendritic cells are pivotal in regulating adaptive immune responses and are essential for T cell-mediated cancer immunity (Gardner & Ruffell, 2016). They play a central role in initiating and modulating both innate and adaptive immune responses (Wculek et al., 2020). Memory CD4+ T cells are involved in the immune response and have been associated with better prognosis in various cancers (Liu et al., 2021). Plasma cells, key effectors of adaptive immunity, produce antibodies that protect the body from pathogens (Varlet et al., 2020). Evidence suggests that plasma cells actively contribute to anti-tumor immunity (Wouters & Nelson, 2018). M0 macrophages are predominantly infiltrated in high-risk patients with CRC, where they suppress tumor immune responses and correlate with poor survival (Castrogiovanni et al., 2022; Dong et al., 2024). MC1R has been shown to stimulate HLA-A2-restricted cytotoxic T lymphocytes, enhancing their ability to recognize peptides naturally processed on melanoma cells, further linking immune modulation to tumor immunity (Salazar-Onfray et al., 1997).
This study investigated the relationship between chemotherapy drug sensitivity and gene expression in a risk model. Notably, FAM151A exhibited a positive correlation with four drugs (AG.014699, camptothecin, NVP.BEZ235, TW.37) (r > 0.4), while CCDC68 was positively associated with two drugs (BMS.708163, X681640) (r < −0.4). AG.014699 (rucaparib), a poly(ADP-ribose) polymerase inhibitor, was shown to protect Schwann cells from cell death and reduce glycolysis, though it did not counteract the disruption of the TCA cycle under high-glucose conditions in the absence of pyruvate (Yako et al., 2024). Relevant studies (Bugajova et al., 2024) suggest that NVP.BEZ235 can reduce ATP production. Treatment with NVP.BEZ235 increases the mitochondrial autophagy-related protein BNIP3, promoting mitochondrial fission, enhancing protease activity, and facilitating the degradation of mitochondrial proteins. Camptothecin, a potent topoisomerase inhibitor, was found to increase the proportion of cells with decreased mitochondrial membrane potential, correlating with lower ATP levels. In colorectal carcinoma cell lines (DLD1 and HCT-116), camptothecin elevated apoptosis rates and significantly reduced cell viability after 24 h of treatment by lowering ATP production and pyruvate levels (Liskova et al., 2022). In conclusion, NVP.BEZ235 may influence mitochondrial fission through BNIP3-mediated mitochondrial autophagy, while camptothecin induces mitochondrial dysfunction and fission by reducing mitochondrial membrane potential and ATP levels. The main mechanism of action of AG.014699 (rucaparib) is not directly related to mitochondrial fission, but its regulation of cell metabolism may indirectly affect mitochondrial function.
Although this study has made some progress, there are still some limitations. First, the predictive power of the model (AUC value of about 0.6) has room for further improvement and needs to be optimized and improved by introducing more and more predictive features. Secondly, there is a lack of further in vivo and in vitro experiments for functional verification, especially the exploration of the specific role and mechanism of MC1R gene in colorectal cancer (CRC) is still insufficient. Therefore, future studies should validate the model predictions through clinical data and in vivo/in vitro experiments (such as mouse models, cell models, gene knockout or overexpression, drug treatment, etc.), and further explore the mechanism by which MC1R regulates mitochondrial fission through cAMP/AMPK/DRP1/MFF axis and promotes CRC development. In addition, further exploration of the functional role of these key genes in CRC progression, as well as their potential value in therapy, is of great significance for the discovery of new immunotherapy and targeted therapy strategies.
Conclusions
In this study, a series of bioinformatics methods were used to identify the prognostic value of mitochondrial fission-related genes in CRC. Three genes (CCDC68, FAM151A, MC1R) were identified as potential risk model genes in CRC. The risk model was constructed, and the molecular mechanism of mitochondrial fission-related genes affecting CRC was further analyzed, providing a new direction for CRC treatment.
Supplemental Information
10.7717/peerj.19522/supp-1Supplemental Information 1FPKM (Fragments Per Kilobase of transcript per Million mapped reads) expression data (log-transformed) for colorectal cancer (CRC) samples from the TCGA cohort 1
10.7717/peerj.19522/supp-2Supplemental Information 2FPKM expression data (log-transformed) for CRC samples from the TCGA cohort 2
10.7717/peerj.19522/supp-3Supplemental Information 3FPKM expression data (log-transformed) for CRC samples from the TCGA cohort 3
10.7717/peerj.19522/supp-4Supplemental Information 4FPKM expression data (log-transformed), group information, and survival case data for CRC samples from the TCGA cohort 4
10.7717/peerj.19522/supp-5Supplemental Information 5Expression and survival data for CRC samples from the GSE103479 dataset
10.7717/peerj.19522/supp-6Supplemental Information 6Experimental data for CCDC68, FAM151A, and MC1R, including gene expression and relevant analysis
10.7717/peerj.19522/supp-7Supplemental Information 7MIQE checklist
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Boatright KM Salvesen GS 2003 Caspase activation Biochemical Society Symposia 23324210.1042/bss 070023314587296 · doi ↗ · pubmed ↗
- 2Bray F Laversanne M Sung H Ferlay J Siegel RL Soerjomataram I Jemal A 2024 Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries CA: A Cancer Journal for Clinicians 7422926310.3322/caac.2183438572751 · doi ↗ · pubmed ↗
- 3Brueckner LM Hess EM Schwab M Savelyeva L 2013 Instability at the FRA 8I common fragile site disrupts the genomic integrity of the KIAA 0146, CEBPD and PRKDC genes in colorectal cancer Cancer Letters 336859510.1016/j.canlet.2013.04.00723603433 · doi ↗ · pubmed ↗
- 4Bugajova M Raudenska M Hanelova K Navratil J Gumulec J Petrlak F Vicar T Hrachovinova S Masarik M Kalfert D Grega M Plzak J Betka J Balvan J 2024 Glutamine and serum starvation alters the ATP production, oxidative stress, and abundance of mitochondrial RN As in extracellular vesicles produced by cancer cells Scientific Reports 142581510.1038/s 41598-024-73943-239468126 PMC 11519472 · doi ↗ · pubmed ↗
- 5Caldwell CM Green RA Kaplan KB 2007 APC mutations lead to cytokinetic failures in vitro and tetraploid genotypes in Min mice Journal of Cell Biology 1781109112010.1083/jcb.20070318617893240 PMC 2064647 · doi ↗ · pubmed ↗
- 6Castrogiovanni P Barbagallo I Imbesi R Musumeci G Sanfilippo C Broggi G Caltabiano R Tibullo D Giallongo C Forte S Volti G Li Di Rosa M 2022 Chitinase domain containing 1 increase is associated with low survival rate and M 0 macrophages infiltrates in colorectal cancer patients Pathology, Research and Practice 23715403810.1016/j.prp.2022.15403835932496 · doi ↗ · pubmed ↗
- 7Chelakkot VS Thomas K Romigh T Fong A Li L Ronen S Chen S Funchain P Ni Y Arbesman J 2023 MC 1R signaling through the c AMP-CREB/ATF-1 and ERK-NFκB pathways accelerates G 1/S transition promoting breast cancer progressionnpj Precision Oncology 78510.1038/s 41698-023-00437-137679505 PMC 10485002 · doi ↗ · pubmed ↗
- 8Chen M Ye K Zhang B Xin Q Li P Kong AN Wen X Yang J 2019 Paris Saponin II inhibits colorectal carcinogenesis by regulating mitochondrial fission and NF-κB pathway Pharmacological Research 13927328510.1016/j.phrs.2018.11.02930471409 PMC 8108001 · doi ↗ · pubmed ↗