Identification of mitochondrial-related subtypes and development of a prognostic model for pancreatic ductal adenocarcinoma

Introduction

Pancreatic cancer is the third leading cause of cancer-related deaths worldwide, accounting for nearly 5% of all cancer-related deaths (1,2). Pancreatic ductal adenocarcinoma (PDAC) is the major histologic subtype of pancreatic cancer, accounting for more than 95% of all pancreatic cancer cases (3). Chemotherapy and surgery are the primary treatment options for PDAC; however, only 15–20% of PDAC patients are eligible for surgery at diagnosis due to the frequent presence of metastasis (4). Thus, improved early detection and therapeutic strategies urgently need to be established. Due to the deep anatomical location of the pancreas and the asymptomatic progression of PDAC, early detection remains challenging (5), and biomarkers for effective diagnosis and therapeutic intervention are lacking (6). Thus, the molecular mechanisms driving PDAC progression and metastasis urgently need to be identified to improve patient outcomes.

Mitochondria, often described as the energy factories of cells, are indispensable for key cellular functions, including energy production, redox balance, reactive oxygen species (ROS) generation, and calcium homeostasis (7). Dysregulated mitochondrial function is a hallmark of cancer, characterized by abnormal morphology, mitochondrial copy number defects, altered energy metabolism, ROS accumulation, and biogenesis and mitophagy imbalances (8). Notably, while severe mitochondrial dysfunction may suppress tumor growth, mild dysfunction can promote tumor proliferation and invasion (9).

Mitochondria play a vital role in the progression of pancreatic cancer. Despite their complex metabolism, the mitochondria of pancreatic cancer cells function well, efficiently producing ATP energy and providing the building blocks required for rapid cell proliferation. Pancreatic cancer cells use mitochondria-related functions to transform surrounding stromal cells and the extracellular matrix to optimize the tumor microenvironment (TME) (5). In addition to serving as metabolic hubs, dynamic changes in mitochondrial shape (particularly the tilt towards division) play a key role in the progression of PDAC. PDAC may benefit more from mitochondrial division than fusion. The division process is often associated with cell proliferation, metabolic reprogramming, and the clearing of damaged mitochondria, which may provide advantages to fast-growing tumor cells (10).

In addition, mitochondrial DNA-related genetic information and regulatory networks are key factors driving the malignant behavior of pancreatic cancer. Mutations in mitochondrial DNA are directly related to chemotherapy resistance and the metastatic process of pancreatic cancer. Alterations in the genetic material of mitochondria are a key factor in the aggressiveness and treatment resistance of pancreatic cancer (11). In PDAC, where low oxygen and nutrient levels are common, mitochondrial dysfunction may drive metabolic programming, remodel the tumor environment, and facilitate immune evasion. Thus, mitochondria could serve as prognostic biomarkers for PDAC. The role of mitochondria has been widely researched in other cancer types, but little is known about their role in PDAC, particularly in terms of the effects of mitochondria on the TME.

Public databases, such as The Cancer Genome Atlas (TCGA), along with advancements in immunotherapy and the development of risk models, have emerged as powerful tools for analyzing cancer outcomes (12-14). Typically, researchers leverage large-scale genomic data from TCGA to identify novel molecular biomarkers, such as specific gene mutations, expression patterns, and immune cell infiltration scores. These biomarkers are then used to construct sophisticated risk models that can accurately predict patient prognosis, therapeutic response, and immunotherapy efficacy. This data-driven approach is revolutionizing our understanding of tumor biology and paving the way for personalized cancer treatment.

This study integrated analyses of mitochondrial-related genes (MitoRGs) and publicly available data to identify the molecular mechanisms underlying PDAC. By constructing risk models, analyzing the TME, and leveraging single-cell data to explore immune-related features, the research aimed to identify key prognostic and therapeutic biomarkers. Our findings provide novel insights into the role of mitochondrial dysfunction in PDAC oncogenesis and progression, potentially improving early detection, risk stratification, and the development of precision treatment strategies. We present this article in accordance with the TRIPOD and MDAR reporting checklists (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-613/rc).

Methods

Data and clinical information

The clinical and genetic data used for the training set were obtained from the TCGA database. The RNA-sequencing transcriptomic data and corresponding clinical information of 178 pancreatic cancer patients were obtained from the TCGA database (https://cancergenome.nih.gov/). The TCGA database includes detailed basic information, as well as clinical, pathological, and follow-up data, all of which were included in this study. Additionally, RNA-sequencing transcriptomic data from 167 normal human pancreatic tissues were downloaded from the Genotype-Tissue Expression (GTEx) project (https://xenabrowser.net/), and the GSE28735 dataset, comprising 45 tumor samples and 45 matched normal tissue samples, was downloaded from the Gene Expression Omnibus (GEO). Additionally, the MitoCarta 3.0 database (https://www.broadinstitute.org/mitocarta/mitocarta30-inventory-mammalian-mitochondrial-proteins-and-pathways) was referenced to access a curated list of 1,136 human protein-coding genes localized to mitochondria. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Functional enrichment analysis of MitoRGs in normal and tumor samples

A differential expression analysis of the MitoRGs between normal and tumor samples was performed using the limma package in R (version 4.2.1). Differentially expressed genes (DEGs) were defined as those with a |log₂ fold change| >0.585 and a false discovery rate <0.05. The results were visualized using a heatmap. Functional enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were then conducted on the DEGs.

Development and validation of risk-score models

The prognostic genes were identified from the DEGs using univariate Cox regression (P<0.05). Mutation and copy number variation analyses were conducted using “maftools” and University of California, Santa Cruz (UCSC) Xena, and the results were visualized in histograms and circle plots. The risk-score model was established by least absolute shrinkage and selection operator (LASSO) Cox regression (“glmnet”) using the following formula:

Riskscore=∑i=1n(Coefi×ExpGenei)[1]

where “Coef” and “ExpGene” were generated by LASSO Cox regression analysis. Survival differences and predictive accuracy were assessed using Kaplan-Meier plots, log-rank tests, time-dependent receiver operating characteristic (ROC) curves (“survivalROC”), and the corresponding areas under the curves (AUCs).

Immune cell infiltration analysis

The CIBERSORT algorithm was used to analyze immune cell infiltration in the TME. The ESTIMATE algorithm was then applied to calculate stromal, immune, and ESTIMATE scores, which indicate tumor purity. Finally, the differences in the microenvironment scores between the risk groups were compared.

ScRNA-seq analysis

The single-cell RNA-sequencing (scRNA-seq) dataset GSE212966, comprising six tumor samples and six normal samples, was used in this study. Data preprocessing and analysis were conducted using the “Seurat” package in R. After initial quality control to remove low-quality cells based on mitochondrial and ribosomal gene content, the gene expression data were log-normalized. The top 2,000 highly variable genes were then identified for downstream dimensionality reduction. Principal component analysis (PCA) was performed on these genes, and the resulting principal components were further integrated using the Harmony package to correct for batch effects. Finally, the cell clusters were visualized in two dimensions using uniform manifold approximation and projection (UMAP) (15). The cell clusters were manually annotated by identifying marker genes.

Cell culture

Human pancreatic ductal epithelial cells of normal origin (HPDE6-C7) were obtained from Immocell Biotechnology Ltd. (Xiamen, China), while pancreatic carcinoma cells (PANC-1) were obtained from Sunncell Biotechnology (Wuhan, China). Both cell lines were authenticated by Short Tandem Repeats (STR) analysis. The HPDE6-C7 and PANC-1 cells were cultured in Dulbecco’s Modified Eagle Medium supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin, under controlled conditions of 37 ℃, 5% carbon dioxide, and 100% relative humidity. Routine passaging and mycoplasma testing were performed on all cell lines, with all tests yielding negative results for contamination.

Real-time quantitative polymerase chain reaction (RT-qPCR)

Total RNA was isolated employing Trizol reagent (catalog number DP419, TIANGEN, China), followed by reverse transcription into complementary DNA (cDNA) using a cDNA synthesis kit (catalog number TSK314S, TSINGKE, China). The cDNA was used for the polymerase chain reaction (PCR) analysis (TSE201, TSINGKE), targeting genes such as ACACB, IFI27, BCL2L1, SUGCT, SLC25A27, SLC25A24, AK4, HTATIP2, ACADVL, MCU, GPD2, MLYCD, and SMDT1, with glyceraldehyde-3-phosphate dehydrogenase as the reference gene. The PCR reaction used a three-step protocol: 95 ℃ pre-denaturation for 1 minute, followed by 40 cycles of 95 ℃ denaturation, 60 ℃ annealing, and 72 ℃ extension, each for 10 seconds. A melting curve analysis was conducted to assess the amplification products, and the relative expression levels were calculated using the 2^-△△Ct method.

Statistical analysis

The statistical analyses were conducted using R software (version 4.2.1). The Wilcoxon rank-sum test was used for comparisons between two groups, while the Kruskal-Wallis test was used for comparisons among multiple groups. Categorical variables were analyzed using the chi-square test or Fisher’s exact test, depending on data distribution. Spearman’s correlation coefficient was used for the correlation analysis, and statistical significance was defined as P<0.05 unless stated otherwise.

Results

Identification of MitoRGs in PDAC through bulk and single-cell transcriptomic analyses

The sequencing data of 178 pancreatic cancer tissue samples were collected from the TCGA database, along with 167 normal pancreatic tissue samples from the GTEx database. The two datasets were merged, and a differential expression analysis comparing the PDAC and normal tissues identified 17,104 DEGs (Figure 1A). Additionally, a differential expression analysis of the GEO sequencing data for pancreatic cancer and adjacent non-tumor tissues revealed 18,776 DEGs (Figure 1B). After combining TCGA expression data with survival information, a univariate Cox regression analysis was conducted to identify the molecular markers related to pancreatic cancer prognosis. Additionally, 1,136 MitoRGs were obtained from the MitoCarta 3.0 database. By intersecting these genes with the DEGs, 13 MitoRGs-DEGs significantly associated with pancreatic cancer prognosis were identified (Figure 1C).

Figure 1 Identification of MitoRGs in PDAC through bulk and scRNA transcriptomic analysis. (A) DEGs in PDAC tissues compared to normal pancreatic tissues using TCGA data. (B) Differential expression analysis of PDAC and adjacent non-tumor tissues using GEO data. (C) Venn diagram showing the intersection of mitochondria-related DEGs and prognostic molecular markers identified in the survival analysis. (D) ScRNA-sequencing analysis clustering cells into 19 groups, highlighting the cell-type heterogeneity in the PDAC samples. (E-G) Annotation of cell clusters into 10 cell types with expression validation of the mitochondrial-related DEGs across clusters using violin plots. (H) Comparison of mitochondrial-related gene expression between HPDE6-C7 (normal cells) and PANC-1 (cancer cells). *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. DEGs, differentially expressed genes; FC, fold change; GEO, Gene Expression Omnibus; MitoRG, mitochondria-related gene; NS, not significant; PDAC, pancreatic ductal adenocarcinoma; scRNA, single-cell RNA; TCGA, The Cancer Genome Atlas.

To examine the expression patterns of these 13 genes, the scRNA-seq data from the GSE212966 dataset were analyzed. A clustering analysis at a resolution of 0.5 grouped the cells into 19 clusters (Figure 1D), which were further classified into 10 cell types, including T cells, fibroblasts, ductal cells, myeloid cells, endothelial cells, stellate cells, proliferating cells, acinar cells, B cells, and cancer stem cells (Figure 1E). Marker genes were used to label the clusters, and violin plots were generated to examine the distribution of the expression levels of the 13 MitoRGs across the various cell types (Figure 1F,1G). Additionally, the gene expression validation in the HPDE and PANC-1 cell lines revealed that most of the genes were more highly expressed in the PANC-1 cells, with the exception of SLC25A27 and ACADVL, which were more highly expressed in the HPDE cells (Figure 1H).

Identification and functional analysis of mitochondria-associated subtypes in PDAC

Using the expression data of the 13 MitoRGs, we conducted unsupervised consensus clustering and identified two distinct pancreatic cancer clusters (k=2) (Figure 2A, Figure S1A-S1L). The clustering analysis revealed that the highest within-cluster correlation was achieved when k=2. The Kaplan-Meier survival analysis showed that patients in cluster A had a significantly better prognosis than those in cluster B (Figure 2B). The PCA revealed a clear distinction between the two clusters (Figure 2C). To explore the functional differences between these subtypes, we conducted a gene set variation analysis (GSVA) using the GEO database, which identified key pathway differences: cluster A was upregulated in pathways such as the regulation of protein kinase A, while cluster B was downregulated in pathways such as the fucosylation pathway and neuroactive ligand-receptor interaction (Figure 2D-2F).

Figure 2 Identification and functional analysis of mitochondria-associated subtypes in PDAC. (A) Unsupervised consensus clustering identified two distinct mitochondrial subtypes (clusters A and B) in PDAC based on 13 MitoRGs. (B) Kaplan-Meier survival analysis revealed better survival outcomes for cluster A compared to cluster B. (C) PCA confirmed distinct separation between the two clusters. (D-F) GSVA compared pathway activity between the clusters. (G,H) GO and KEGG pathway enrichment analyses of mitochondrial-related DEGs. (I) Hallmark pathway enrichment analysis identified key differences between the clusters. DEGs, differentially expressed genes; GO, Gene Ontology; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MitoRG, mitochondria-related gene; PCA, principal component analysis; PDAC, pancreatic ductal adenocarcinoma; TCGA, The Cancer Genome Atlas.

The GO enrichment analysis revealed that these mitochondria-related DEGs were involved in biological processes such as epidermis development, keratinization, keratinocyte differentiation, skin development, and cornified envelope formation (Figure 2G). The KEGG pathway analysis further revealed that these DEGs were significantly enriched in pathways such as calcium signaling, the cell cycle, neuroactive ligand-receptor interaction, p53 signaling, and proteasome (Figure 2H). Finally, the Hallmark gene set enrichment analysis showed that these mitochondrial subtypes were significantly associated with pathways such as E2F targets, G2M checkpoint, glycolysis, interferon (IFN) alpha response, and mitotic spindle (Figure 2I).

Immune microenvironment and key prognostic genes identified through differential expression and LASSO regression analyses in PDAC clusters

We then characterized the immune microenvironment, which revealed stark contrasts between clusters A and B. As shown in the immune cell infiltration landscape (Figure 3A), cluster A was significantly more immunogenic, exhibiting significantly elevated levels of key adaptive and innate immune cells, including activated B cells and CD8⁺ T cells, as well as monocytes and neutrophils. In contrast, the immune profile of cluster B was dominated by CD56 natural killer cells and T helper cells. This clear dichotomy revealed the existence of two distinct immune phenotypes. To explore the molecular basis of these differences, a differential expression analysis was performed, which identified a substantial number of 18,312 DEGs between the clusters (Figure 3B; thresholds: |log₂ fold change| >0.585 and adjusted P<0.05).

Figure 3 Immune microenvironment and key prognostic genes identified by differential expression and LASSO regression analyses. (A) Immune cell infiltration analysis showed higher levels of activated immune cells in cluster A compared to cluster B. (B) Volcano plot of DEGs between the clusters. (C,D) GO and KEGG pathway enrichment analyses for DEGs. (E,F) LASSO regression analysis identified seven key prognostic genes and their optimal lambda values for model construction. ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001. BP, biological process; CC, cellular component; DEG, differentially expressed gene; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; MF, molecular function; NS, not significant.

Subsequent functional enrichment analyses of these DEGs provided insights into their potential biological roles. The GO analysis (Figure 3C) indicated that the most significantly enriched terms were associated with wound healing (biological process), the collagen-containing extracellular matrix (cellular component), and receptor ligand activity (molecular function). Complementing this, the KEGG pathway analysis (Figure 3D) identified the PI3K-AKT signaling pathway as the most significantly enriched pathway, followed by pathways related to the cytoskeleton in muscle cells and human papillomavirus infection. Together, these findings suggest that the distinct immune microenvironments are underpinned by widespread transcriptomic differences influencing processes such as tissue remodeling and key oncogenic signaling.

Notably, the PI3K-AKT signaling pathway, which plays a critical role in regulating cell survival, proliferation, and metabolism, may underlie the distinct biological behaviors observed between the two clusters. To identify the key prognostic genes, we conducted a LASSO regression analysis, with the variable trajectories shown in Figure 3E. As the lambda value increased, the number of non-zero coefficients gradually decreased to zero. Using 10-fold cross-validation, the model achieved optimal performance when the number of genes was reduced to seven, as indicated by the confidence intervals for each lambda value (Figure 3F). A formula was then derived from the coefficients and expression levels of genes to calculate the risk score. The formula was expressed as follows: risk score = (−0.6693 * MTURN expression) + (−0.6345 * TLE2 expression) + (0.3274 * CAV2 expression) + (0.2580 * IL1R2 expression) + (−0.7065 * ATP1A3 expression) + (0.4463 * CLDN1 expression) + (0.3567 * FOXM1 expression).

Validation and evaluation of the PDAC prognostic risk-score model

To evaluate the prognostic significance of the risk-score model, we combined data from TCGA and GEO databases (Figure 4A), and split the merged dataset into a training set (Figure 4B) and validation set (Figure 4C). The patients were divided into high- and low-risk groups based on the median risk score. In the combined model, the patients in the high-risk group showed significantly worse overall survival (OS) than those in the low-risk group (P<0.001). A heatmap illustrated the expression patterns of seven core genes, highlighting clear differences between the two groups (Figure 4D-4F). A distribution plot revealed a significantly higher number of deaths in the high-risk group (Figure 4G-4I).

Figure 4 Validation and evaluation of the prognostic risk-score model. (A-C) PDAC patients’ characteristics by high- and low-risk categories and survival time. (D-F) Heatmaps showing expression patterns of core genes and survival outcomes in risk groups. (G-I) PDAC patients’ characteristics by risk-score categories. (J-L) Survival probability of combined, training, and validation sets. (M-O) ROC curves indicating the predictive accuracy of the risk-score model across multiple datasets, with AUC values for 1-, 3-, and 5-year OS predictions. AUC, area under the curve; OS, overall survival; PDAC, pancreatic ductal adenocarcinoma; ROC, receiver operating characteristic.

The survival analysis further confirmed a significant prognostic difference between the two groups, as shown by the survival curve (P<0.001) (Figure 4J-4L). The AUC values for predicting 1-, 3-, and 5-year OS were 0.781, 0.862, and 0.860, respectively (Figure 4M). Our model demonstrated highly competitive predictive performance compared to those reported in previous studies (16,17). In the training set, the AUC values demonstrated good predictive accuracy (Figure 4N). To further confirm the robustness of our model, we tested it on a validation set. Consistent with the results from the training set, the patients in the high-risk group exhibited significantly worse OS. The ROC curve analysis indicated acceptable predictive performance, with AUC values for predicting 1-, 3-, and 5-year OS of 0.690, 0.759, and 0.766, respectively (Figure 4O), indicating that the model had acceptable reliability and robustness.

Immune cell infiltration and TME characteristics associated with risk groups

Building on the prognostic risk model, we investigated the relationship between risk scores and the TME, focusing on immune cell infiltration and activity. Using the CIBERSORT algorithm, we assessed the expression profiles of 23 immune cell types in the TCGA-PAAD cohort. Statistically significant differences in immune cell infiltration were observed between the high- and low-risk groups, highlighting distinct immune states in their TMEs (Figure 5A).

Figure 5 Immune cell infiltration and TME characteristics associated with risk groups. (A) Immune cell subtype analysis across high- and low-risk groups using CIBERSORT. (B) Immune cell fraction differences highlighting the immunosuppressive microenvironment in the high-risk group. (C) Comparison of TME, stromal, and immune scores between risk groups. (D-M) Correlations between risk scores and immune cell abundance, showing reduced anti-tumor immune activity in the high-risk groups. ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001. TME, tumor microenvironment.

In the high-risk group, anti-tumor immune cells (e.g., activated B cells, activated CD8⁺ T cells, and monocytes) exhibited lower infiltration levels, suggesting a diminished anti-tumor immune response. Conversely, immune-suppressive or tumor-promoting cells (e.g., activated CD4⁺ T cells and Th2 cells) were more abundant, indicating an immunosuppressive microenvironment that may promote tumor progression. These observations suggest that high-risk patients may benefit from immune-targeting therapies aimed at reversing immune suppression.

An analysis of the immune cell fractions in the two groups (Figure 5B) further substantiated the immunosuppressive characteristics of the high-risk group. Key anti-tumor components, including naive CD8⁺ T cells and B cells, were lacking, while immunosuppressive cell types, such as M0 macrophages, were more prevalent.

We further assessed the TME composition using TME, stromal, and immune scores, all of which were consistently higher in the low-risk group (Figure 5C). These results indicated that the low-risk group possessed a more immune-active microenvironment, whereas the high-risk group exhibited a TME more conducive to tumor growth. In the high-risk group, abundant stromal content and suppressed immune infiltration collectively formed an immunosuppressive environment that limit anti-tumor immunity.

The high-risk patients also displayed increased scores in antigen-presenting cell (APC) co-inhibition, alongside lower cytolytic activity and T-cell co-stimulation scores, collectively indicating a suppressive immune microenvironment. Additionally, higher scores in parainflammation and type 1 interferon response in the high-risk group imply a chronic inflammation state that supports tumor growth rather than an effective anti-tumor response. Although the major histocompatibility complex and APC co-stimulation scores were higher in the high-risk group, the elevated co-inhibition signals indicate that antigen presentation was significantly suppressed, thereby restricting T-cell activation. For the high-risk patients, therapeutic approaches that alleviate immune suppression and enhance T-cell activation may be necessary to stimulate a more effective anti-tumor immune response (Figure 5D).

The trend analysis of immune cell abundance in relation to the risk score revealed a gradual decrease in the abundance of monocytes, CD4 memory T cells, and CD8 T cells as the risk score increased (Figure 5E-5M). This indicated that in the high-risk patients, the activity or quantity of these key anti-tumor cells was relatively low, potentially weakening the anti-tumor immune response. Conversely, follicular helper T-cell abundance showed a rising trend with increasing risk scores, implying that in the TME of patients with elevated risk, these cells may facilitate tumor progression and immune escape through heightened immunoregulation or immunosuppressive mechanisms. These shifts in immune cell composition indicate that the TME of high-risk patients is more inclined to suppress anti-tumor responses, potentially facilitating tumor progression and survival.

Gene mutation correlation analysis and nomogram development

The heatmap revealed potential associations between specific genes and immune cell infiltration in the TME. For example, ATP1A3 and TLE2 showed a strong positive correlation with CD8⁺ T cells, indicating that these genes might contribute to promoting anti-tumor immune responses (Figure S2A). Each panel corresponds to a specific gene, illustrating its expression distribution across all cell types (Figure S2B). Notably, IL1R2 was highly expressed in myeloid cells, suggesting that IL1R2 may exert immunosuppressive effects in PDAC.

An enrichment analysis was performed using the GO, KEGG, and HALLMARK databases to identify the key pathways associated with risk scores. The enrichment results for the GO (Figure S3A), KEGG (Figure S3B), and HALLMARK (Figure S3C) datasets are presented for the high-risk group. Similarly, the enrichment results for the GO (Figure S3D), KEGG (Figure S3E), and HALLMARK (Figure S3F; datasets are presented for the low-risk group). Finally, the differentially enriched pathways between these two groups were examined using the GO (Figure S3G), KEGG (Figure S3H), and HALLMARK (Figure S3I) datasets.

To expand on the analysis of the TME characteristics, we further explored the relationship between somatic mutation profiles, MitoRG typing, and risk scores in the pancreatic cancer patients. We retrieved the integrated somatic mutational profiles from the TCGA database for pancreatic cancer patients categorized into high- and low-risk groups. Subsequently, we performed an analysis and visualization of the mutation spectra. We also examined the genetic mutations present in each sample, using distinct colors to represent different mutation types. The analysis of the somatic mutation data from the TCGA database revealed that KRAS, TP53, and SMAD4 were the most commonly mutated genes, especially in the high-risk group, which also exhibited a greater tumor mutation burden (TMB) (Figure 6A-6C). The patients were further stratified into two clusters, A and B, based on MitoRG typing, represented in blue and orange, respectively. A positive correlation was found between risk scores and TMB, with the majority of the cluster B samples clustering in the high-risk, high-TMB region (Figure 6D).

Figure 6 Correlation between risk scores, TMB, and prognosis. (A,B) Mutation analysis of high- and low-risk groups, identifying key mutated genes such as KRAS, TP53, and SMAD4, and (C) their association with TMB. (D-F) Stratification of patients into MitoRG-based clusters, showing higher risk scores and poorer prognosis in cluster B. (G,H) Nomogram model integrating clinical features and risk scores to predict patient outcomes. *, P<0.05; ***, P<0.001. AIC, Akaike information criterion; CI, confidence interval; MitoRG, mitochondria-related gene; TMB, tumor mutation burden.

Risk scoring showed that the cluster B patients had significantly higher risk scores, which were associated with worse prognoses (Figure 6E,6F). A forest plot identified age, N1 stage, and high-risk scores as unfavorable prognostic factors (Figure 6G). To enhance prognosis prediction, a nomogram model incorporating gender, age, stage (T/N), and risk scores (Figure 6H) was established. These findings emphasize the importance of somatic mutations, TMB, and MitoRG typing in assessing pancreatic cancer risk and prognosis.

Discussion

This study examined the molecular characteristics of PDAC and highlighted the critical role of mitochondria in cancer progression. The seminal work of Collisson et al. in 2011, which categorized pancreatic cancer into classical, quasi-mesenchymal, and exocrine types using transcription profiling, has prompted further investigations (18). These subsequent studies, relying on specific genomic and transcriptomic expression signatures, have led to the development of a more detailed molecular classification system for pancreatic cancer (19,20). Through gene expression and survival analyses, we classified PDAC patients into two molecular subtypes (clusters A and B). The patients in cluster A had significantly better survival outcomes. Additionally, we developed a risk-score model based on seven MitoRGs, which effectively stratified the patients based on risk scores (AUC for the training cohort: 0.781–0.860) and demonstrated robust performance in the external validation dataset. These findings provide key insights into the role of mitochondrial function in PDAC and offer a foundation for novel therapeutic strategies.

Mitochondria are not only central to cellular energy metabolism but also play essential roles in apoptosis, oxidative stress, and signal transduction. Our study revealed that MitoRG expression significantly influences both PDAC patient survival and molecular subtype classification. The patients in cluster A exhibited higher immune activity and better prognoses, potentially due to increased immune cell infiltration in their TME. Conversely, the immunosuppressive microenvironment observed in cluster B, characterized by the reduced infiltration of activated CD8 T cells, likely promotes tumor progression. These findings underscore the central role of mitochondrial dysfunction in PDAC tumorigenesis and immune regulation.

Mitochondria-mediated abnormalities are not unique to PDAC and have also been reported in various cancers, including breast cancer (21,22), hepatocellular carcinoma (23), and glioblastoma (24). Such abnormalities are particularly evident in lung adenocarcinoma (LUAD), where the high expression of oxidative phosphorylation (OXPHOS)-related genes has been linked to a poor prognosis. As a core mitochondrial metabolic pathway, OXPHOS is critical for ATP production via the electron transport chain. The variation in OXPHOS levels among tumors and in different regions of the same tumor suggests that it plays a complex role in cancer progression. Based on these findings, targeting OXPHOS metabolism has been proposed as a novel therapeutic strategy for LUAD, underscoring the broader relevance of mitochondrial metabolism as a promising area for therapeutic intervention across multiple cancer types (25). These shared characteristics across cancer types underscore the centrality of mitochondrial dysfunction in tumor biology and reinforce its relevance in PDAC progression and prognosis.

The clinical relevance of MitoRGs has been demonstrated not only in PDAC but also in other diseases, such as idiopathic pulmonary fibrosis (IPF). A recent study successfully established and validated a prognostic model for IPF based on MitoRGs, emphasizing their pivotal role in predicting disease outcomes (26). Similarly, our study developed a mitochondria-based risk-score model for PDAC, which revealed that high-risk patients exhibit an immunosuppressive TME characterized by reduced anti-tumor immune activity and elevated T-cell co-inhibition scores. These features likely limit PDAC response to conventional immunotherapies.

Research has also highlighted the therapeutic potential of targeting mitochondrial pathways in various cancers. For instance, in laryngeal squamous cell carcinoma (27), MitoRGs have been used to develop prognostic models that reveal significant differences in immune function and drug sensitivity between risk groups.

Furthermore, the immunosuppressive microenvironment in PDAC is influenced by cancer-associated fibroblasts (CAFs), particularly the myofibroblastic CAF phenotype defined by MAPK signaling. This phenotype has been shown to reduce CD8⁺ T-cell infiltration and promote immunoregulatory programs, further contributing to tumor progression and therapy resistance (28).

These findings underscore the importance of integrating mitochondrial dysfunction and CAF heterogeneity in understanding PDAC progression. Our mitochondria-based risk-score model offers a promising tool to guide personalized treatment strategies, such as combining immune checkpoint inhibitors with agents targeting mitochondrial pathways or CAF phenotypes, to improve therapeutic outcomes.

A key strength of this study lies in the use of multi-database validation. We constructed and validated the risk-score model using TCGA and GEO datasets, and employed various bioinformatics tools (e.g., LASSO regression and CIBERSORT immune analysis) to ensure the robustness of our results. Additionally, by incorporating single-cell data (GSE212966), we further explored the expression patterns of MitoRGs across different cells, uncovering their complex relationships with the TME. Our findings suggest that the risk-score model could serve as an important tool for predicting PDAC prognosis and informing treatment strategies. For high-risk patients, therapeutic approaches targeting the immunosuppressive TME, such as immune checkpoint inhibitors combined with metabolic targeting agents, may hold promise. For low-risk patients, conventional therapies combined with strategies to enhance anti-tumor immunity may be more effective. Moreover, the identification of MitoRGs as potential therapeutic targets opens up new avenues for drug development. Research has highlighted the critical role of mitochondrial dysfunction in cancer treatment. For instance, gambogic acid was shown to induce pyroptosis in ovarian cancer cells via the ROS/p53/mitochondria signaling pathway (29). These findings further underscore the potential of targeting mitochondrial pathways to regulate cell death mechanisms and enhance therapeutic outcomes in cancers such as ovarian cancer.

The findings of this study align with previous research on PDAC subtyping and immune microenvironments. Molecular markers such as HMGA2 and GATA6 have been found to be associated with distinct PDAC subtypes, with HMGA2⁺/GATA6⁻ tumors showing lower CD8⁺ T-cell infiltration, increased FAP⁺ fibroblasts, and poorer responses to gemcitabine-based chemotherapy (30). These markers highlight the role of the immune microenvironment in shaping patient outcomes. Although Huang et al. (31) also developed risk-score models, their AUC values for 1-, 3-, and 5-year predictions were 0.587, 0.635, and 0.770, respectively. Conversely, those for our model were 0.781, 0.862, and 0.860 for the same time points. The risk-score models developed by Shiye et al. for PDAC, based on stemness and ferroptosis, demonstrated high predictive accuracy for OS, achieving AUC values of 0.784, 0.787, and 0.793 for 1-, 2-, and 3-year predictions, respectively (17). Peng et al. constructed a novel lactylation-related signature for prognostic prediction in pancreatic adenocarcinoma. The model demonstrated robust predictive accuracy for OS in the training cohort, with AUC values of 0.722, 0.793, and 0.782 for 1-, 3-, and 5-year OS, respectively (16). Overall, our model demonstrated superior predictive performance and stability compared to previous studies. Collectively, these results underscore the importance of mitochondrial dysfunction in PDAC progression and highlight the distinct contributions of this study to the field.

Furthermore, while the primary goal of this study was to establish a prognostic tool for stratified management, we must emphasize the critical importance of shifting the paradigm toward early detection and prevention for PDAC. Given that the dismal prognosis of PDAC is largely attributable to late-stage diagnosis, the ultimate clinical utility of molecular markers may lie in their ability to identify the disease in its pre-symptomatic or early stages. The molecular signature underpinning our risk score holds promise for such a translational application. Future studies should seek to validate the key components of our model in high-risk cohorts, such as individuals with a family history of PDAC, hereditary syndromes, or new-onset diabetes. By developing a minimally invasive, blood-based assay (e.g., to detect circulating tumor DNA or proteins related to our signature), a predictive molecular diagnosis strategy could be implemented to screen these populations. The successful identification of individuals with early, potentially curable lesions (e.g., high-grade pancreatic intraepithelial neoplasia) would represent a monumental advance in PDAC management, moving the intervention point from incurable late-stage disease to a preventable or curable early stage.

Despite the significant findings, our study had several limitations. First, our analysis was primarily based on retrospective public datasets. While this approach allows for robust bioinformatic exploration and hypothesis generation, it is inherently constrained by the variables and cohorts available in these datasets. Second, and more importantly, our findings were derived from computational analyses and lack experimental validation in functional models, such as PDAC organoids or in vivo animal studies. Thus, the causal relationship between the identified mitochondrial dysfunction and PDAC progression remains to be experimentally confirmed. Finally, the cohorts used in this study, while substantial, might not fully represent the diversity of PDAC across different ethnicities and healthcare settings, potentially limiting the immediate generalizability of our biomarker signatures to all patient populations. Future research should seek to validate these findings prospectively in multi-center cohorts and employ advanced functional models to delineate the precise mechanistic roles of the key MitoRGs identified in this study.

Conclusions

Our findings have important transformative significance for the development of PDAC precision medicine. Currently, several types of mitochondria-targeting drugs are in preclinical and clinical research stages. These include metabolic regulators that disrupt mitochondrial metabolism (e.g., metformin), compounds that leverage increased mitochondrial membrane potential in cancer cells (e.g., mitochondria-targeted tamoxifen), drugs that inhibit key mitochondrial enzymes, and drugs that target mitochondrial anti-apoptotic proteins (e.g., Bcl-2 inhibitors). Crucially, our research provides a strategic framework for improving the effectiveness of such therapies. By identifying specific subgroups of PDAC patients with significant mitochondrial dysfunction and an immunosuppressive TME, our risk-score model provides a powerful tool for patient selection. We hypothesize that “high-risk” patients, as defined by our model, may be particularly sensitive to mitochondria-targeted drugs. Future research should seek to test this hypothesis in clinical trials, where our model can be used to screen patient populations most likely to respond.

Acknowledgments

We appreciate the TCGA databases for providing their platforms, and thank contributors for their valuable public dataset.

Footnote

Reporting Checklist: The authors have completed the TRIPOD and MDAR reporting checklists. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-613/rc

Data Sharing Statement: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-613/dss

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2025-613/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-613/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. Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA Cancer J Clin 2024;74:12-49. [Crossref] [PubMed]
2. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
3. Mizrahi JD, Surana R, Valle JW, et al. Pancreatic cancer. Lancet 2020;395:2008-20. [Crossref] [PubMed]
4. Halbrook CJ, Lyssiotis CA, Pasca di Magliano M, et al. Pancreatic cancer: Advances and challenges. Cell 2023;186:1729-54. [Crossref] [PubMed]
5. Padinharayil H, Rai V, George A. Mitochondrial Metabolism in Pancreatic Ductal Adenocarcinoma: From Mechanism-Based Perspectives to Therapy. Cancers (Basel) 2023;15:1070. [Crossref] [PubMed]
6. Garg SK, Chari ST. Early detection of pancreatic cancer. Curr Opin Gastroenterol 2020;36:456-61. [Crossref] [PubMed]
7. Li Y, Zhang H, Yu C, et al. New Insights into Mitochondria in Health and Diseases. Int J Mol Sci 2024;25:9975. [Crossref] [PubMed]
8. Bao X, Zhang J, Huang G, et al. The crosstalk between HIFs and mitochondrial dysfunctions in cancer development. Cell Death Dis 2021;12:215. [Crossref] [PubMed]
9. Li Y, Sundquist K, Zhang N, et al. Mitochondrial related genome-wide Mendelian randomization identifies putatively causal genes for multiple cancer types. EBioMedicine 2023;88:104432. [Crossref] [PubMed]
10. Carmona-Carmona CA, Dalla Pozza E, Ambrosini G, et al. Divergent Roles of Mitochondria Dynamics in Pancreatic Ductal Adenocarcinoma. Cancers (Basel) 2022;14:2155. [Crossref] [PubMed]
11. Moro L. Mitochondrial DNA and MitomiR Variations in Pancreatic Cancer: Potential Diagnostic and Prognostic Biomarkers. Int J Mol Sci 2021;22:9692. [Crossref] [PubMed]
12. Gössling GCL, Zhen DB, Pillarisetty VG, et al. Combination immunotherapy for pancreatic cancer: challenges and future considerations. Expert Rev Clin Immunol 2022;18:1173-86. [Crossref] [PubMed]
13. Liang G, He J, Chen T, et al. Identification of ALDH7A1 as a DNA-methylation-driven gene in lung squamous cell carcinoma. Ann Med 2025;57:2442529. [Crossref] [PubMed]
14. Ren H, Du MZ, Liao Y, et al. Deciphering the Significance of Platelet-Derived Chloride Ion Channel Gene (BEST3) Through Platelet-Related Subtypes Mining for Non-Small Cell Lung Cancer. J Cell Mol Med 2024;28:e70233. [Crossref] [PubMed]
15. Becht E, McInnes L, Healy J, et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol 2018; Epub ahead of print. [Crossref] [PubMed]
16. Peng T, Sun F, Yang JC, et al. Novel lactylation-related signature to predict prognosis for pancreatic adenocarcinoma. World J Gastroenterol 2024;30:2575-602. [Crossref] [PubMed]
17. Ruan S, Wang H, Zhang Z, et al. Identification and validation of stemness-based and ferroptosis-related molecular clusters in pancreatic ductal adenocarcinoma. Transl Oncol 2024;41:101877. [Crossref] [PubMed]
18. Collisson EA, Sadanandam A, Olson P, et al. Subtypes of pancreatic ductal adenocarcinoma and their differing responses to therapy. Nat Med 2011;17:500-3. [Crossref] [PubMed]
19. Moffitt RA, Marayati R, Flate EL, et al. Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma. Nat Genet 2015;47:1168-78. [Crossref] [PubMed]
20. Bailey P, Chang DK, Nones K, et al. Genomic analyses identify molecular subtypes of pancreatic cancer. Nature 2016;531:47-52. [Crossref] [PubMed]
21. de Oliveira RC, Dos Reis SP, Cavalcante GC. Mutations in Structural Genes of the Mitochondrial Complex IV May Influence Breast Cancer. Genes (Basel) 2023;14:1465. [Crossref] [PubMed]
22. Fang Y, Zhang Q, Guo C, et al. Mitochondrial-related genes as prognostic and metastatic markers in breast cancer: insights from comprehensive analysis and clinical models. Front Immunol 2024;15:1461489. [Crossref] [PubMed]
23. Zhang X, Guo J, Jabbarzadeh Kaboli P, et al. Analysis of Key Genes Regulating the Warburg Effect in Patients with Gastrointestinal Cancers and Selective Inhibition of This Metabolic Pathway in Liver Cancer Cells. Onco Targets Ther 2020;13:7295-304. [Crossref] [PubMed]
24. Cunha de Oliveira R, Gouvea de Souza F, Bispo AG, et al. Differential gene expression analysis supports dysregulation of mitochondrial activity as a new perspective for glioblastoma's aggressiveness. Heliyon 2024;10:e40414. [Crossref] [PubMed]
25. Fujiwara M, Mimae T, Kushitani K, et al. Mitochondrial Metabolism as a Potential Novel Therapeutic Target for Lung Adenocarcinoma. Anticancer Res 2024;44:5241-52. [Crossref] [PubMed]
26. Wang X, Yang L, Wang Y, et al. Establishment and validation of a prognostic model for idiopathic pulmonary fibrosis based on mitochondrial-related genes. J Thorac Dis 2024;16:7427-45. [Crossref] [PubMed]
27. Wu Z, Chen Y, Jiang D, et al. Mitochondrial-related drug resistance lncRNAs as prognostic biomarkers in laryngeal squamous cell carcinoma. Discov Oncol 2024;15:785. [Crossref] [PubMed]
28. Veghini L, Pasini D, Fang R, et al. Differential activity of MAPK signalling defines fibroblast subtypes in pancreatic cancer. Nat Commun 2024;15:10534. [Crossref] [PubMed]
29. Zhang D, Chen Y, Sun Y, et al. Gambogic acid induces GSDME dependent pyroptotic signaling pathway via ROS/P53/Mitochondria/Caspase-3 in ovarian cancer cells. Biochem Pharmacol 2025;232:116695. [Crossref] [PubMed]
30. Yamamoto N, Dobersch S, Loveless I, et al. HMGA2 Expression Predicts Subtype, Survival, and Treatment Outcome in Pancreatic Ductal Adenocarcinoma. Clin Cancer Res 2025;31:733-45. [Crossref] [PubMed]
31. Huang Z, Li M, Gu B, et al. Ferroptosis-related LINC02535/has-miR-30c-5p/EIF2S1 axis as a novel prognostic biomarker involved in immune infiltration and progression of PDAC. Cell Signal 2024;123:111338. [Crossref] [PubMed]

(English Language Editor: L. Huleatt)