This article is mentioned in:

Introduction

Pancreatic adenocarcinoma (PAAD) is a highly lethal malignancy, with 80–85% of patients presenting with unresectable disease at advanced stages due to the unique anatomical structure of the pancreas. There were ~60,430 newly diagnosed cases in the United States in 2021. The incidence of PDAC is increasing by 0.5 to 1.0% per year, and it is expected to become the second leading cause of cancer-related deaths by 2030 (1). The incidence of PDAC is increasing by 0.5 to 1.0% per year, and it is expected to become the second leading cause of cancer-related deaths by 2030. The 5-year overall survival rate is only ~4% (2). Treatment approaches for pancreatic cancer have primarily included surgical resection, targeted therapy, immune checkpoint inhibitors (ICIs) and combination therapy. The current first-line standard of care consists of FOLFIRINOX (a combination of 5-fluorouracil, leucovorin, irinotecan and oxaliplatin) or gemcitabine plus albumin-bound (nab) paclitaxel (3). However, as most patients are diagnosed at advanced or metastatic stages, early pancreatic cancer remains nearly asymptomatic. Consequently, >80% of patients are unable to undergo surgical intervention, and available treatments offer limited efficacy. Local or distant recurrence is common, alongside intrinsic drug resistance and poor treatment outcomes (4,5). These challenges highlight the urgent need for more effective therapeutic targets and treatment strategies to improve patient survival and quality of life. Therefore, researchers have increasingly employed bioinformatics approaches to identify potential therapeutic targets. Bioinformatics offers a powerful tool for large-scale data analysis and screening, providing valuable insights into potential treatment targets (6). The present study integrated bioinformatics analysis with experimental validation to explore novel therapeutic targets for pancreatic cancer, aiming to introduce new perspectives and methods for its treatment.

Furthermore, changes in the extracellular matrix (ECM) serve a critical role in pancreatic cancer progression. One of the key functions of the ECM is to maintain tissue homeostasis in response to injury through its immune, vascular and connective tissue components. However, mutations in the oncogene KRAS and the tumor suppressor gene tumor protein 53 (TP53) (7,8) lead to excessive connective tissue proliferation, forming a mechanical barrier around tumor cells. This barrier disrupts proper blood vessel formation, limits chemotherapy exposure and reduces immune cell infiltration (9). Given the significance of the tumor microenvironment (TME) in tumorigenesis, research efforts have increasingly focused on targeting specific TME components. These approaches can be divided into three phases: i) Ablation restricts the evolution of cells that limit the matrix barrier of drug delivery and arrest ECM components; ii) modification strategies for immune responses; and iii) remodeling of cancer cell metabolism (Fig. 1). Early research aimed to break down the stromal barrier that restricts drug delivery and blocks cellular interactions with ECM components (10–15). However, several of these studies overlooked the dual role of connective tissue (16), the heterogeneity of stromal components (17) and the diversity of cancer-associated fibroblasts (CAFs) (16,18). These oversights contributed to discrepancies between preclinical and clinical outcomes, and most therapies failed to demonstrate clear benefits in clinical applications. When combined with conventional chemotherapy, certain treatments even resulted in higher incidences of thrombotic events and other adverse reactions (19). The inability to effectively eliminate the stromal barrier restricting drug delivery prompted a shift in research focus. Although PAAD is often considered an immunologically ‘cold’ tumor, with only a subset of cases exhibiting immune activity (20), numerous studies have established a strong connection between immune processes and PAAD carcinogenesis (21,22). This led to an increased focus on immunotherapeutic strategies (23–36).

Figure 1. Development history of TME treatment strategies. TME, tumor microenvironment; FAK, focal adhesion kinase; FAP, fibroblast activation protein-α; SHH, sonic hedgehog.

Despite the low immunological activity of PAAD and the limitations of overly focusing on inhibition within the TME, the ‘immune response modification strategy’ has produced mixed results in clinical applications. Certain treatments have raised concerns about pulmonary toxicity, yet, overall, the strategy has yielded promising outcomes. With a growing understanding of TME-targeted therapies, the emphasis has shifted from mere inhibition to comprehensive remodeling of the TME (37,38). Consequently, previous research has focused on ‘metabolic remodeling of cancer cells’ (37–41). Research now indicates that the stroma is multifaceted, highlighting the limitations of targeting isolated elements for TME remodeling. Combination therapies addressing multiple tumor phenotypes have shown broad potential. For example, the ‘immune response modification strategy’, which targets CD40 co-stimulatory molecules, has been associated with substantial stromal degradation whilst simultaneously increasing CD86 and major histocompatibility complex class II molecule expression in TAMs (27). Similarly, glutamine antagonism in the ‘remodeling of cancer cell metabolism’ has been associated with reduced hyaluronic acid, increased CD8+ T-cell infiltration and increased sensitivity to anti-programmed cell death protein 1 (PD-1) therapy within the TME (42). In summary, combinatorial and multimodal approaches targeting multiple aspects of the TME offer notable therapeutic potential. Given the limitations of single-target interventions in addressing the complexity of the PAAD microenvironment, the present study developed a multi-gene prognostic model. This integrated framework draws on multidimensional features to create a novel analytical tool for guiding combination therapy.

Moreover, the present study integrated GSE28735 differential expression profiles with The Cancer Genome Atlas (TCGA) multi-omics data to develop a multigene prognostic model for pancreatic ductal adenocarcinoma (PDAC) using the Cox-least absolute shrinkage and selection operator (LASSO) algorithm. The risk score-based stratification system effectively predicted patient survival outcomes, clinicopathological characteristics and immunotherapy responses. Its clinical utility as an independent prognostic indicator was further validated in an independent cohort, demonstrating its potential as a quantitative tool for personalized therapeutic decision-making. By uncovering the heterogeneous characteristics of the tumor microenvironment, this model offers new insights into optimizing ICI treatment strategies and contributes to advancing precision medicine for pancreatic cancer.

Materials and methods

Data sources, differential expression analysis and functional enrichment

The study flowchart is presented in Fig. 2. The present study utilized the RNA-sequencing (RNA-seq) dataset of PAAD from the TCGA database (https://portal.gdc.cancer.gov/), along with its corresponding mutation dataset (available on the TCGA official website as of June 30, 2025). The RNA-seq dataset included a total of 181 samples, comprising 4 non-tumor (11A) samples and 177 tumor (01A) samples, which were used for prognostic gene exploration and model construction. Additionally, the mutation dataset contained 180 samples for mutation landscape analysis. From Gene Expression Omnibus, the GSE28735 dataset (11A=45, 01A=45; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE28735) was employed for differential expression analysis and validation of the predictive models constructed using the TCGA dataset.

Figure 2. Flow chart of the present study. LASSO, least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment; RT-qPCR, reverse transcription-quantitative PCR.

Differential expression analysis of mRNAs was performed using the limma package in R software (https://bioconductor.org/packages/limma). Genes with an adjusted P<0.05 and a log2(fold change) >1 or <-1 were considered significantly differentially expressed. The identified differentially expressed genes (DEGs) were subsequently visualized using volcano and heat maps. Functional enrichment analyses were systematically performed on both upregulated and downregulated gene sets using Gene Ontology (GO; http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) pathway analyses.

Construction and evaluation of predictive models

A univariate Cox analysis was first performed to identify genes associated with prognosis (P<0.05). LASSO regression was then applied to refine the selection of key genes influencing patient outcomes. A prognostic model was developed using the training cohort (TCGA dataset), incorporating these genes and their corresponding coefficients. Patients were stratified into high- and low-risk groups based on the median risk score. Survival analyses were performed for both groups, and the predictive accuracy of the model was evaluated.

Construction of nomograms

Nomograms were constructed to visualize the results of Cox regression. A scoring system was developed based on the magnitude of the regression coefficients for all independent variables, allowing for the calculation of a total score for each patient to estimate the probability of clinical outcomes. The nomograms were used to predict the 1-, 3- and 5-year mortality rates in patients (TCGA dataset). Receiver operating characteristic (ROC) curves were then generated to assess the predictive performance of the model.

External validation of predictive models

The predictive model was further validated using the GSE28735 cohort. Patients were categorized into high- and low-risk groups based on the median risk score. Survival analyses were performed, and the accuracy of the model and its value as an independent prognostic indicator were assessed.

Assessment of immune infiltration and prediction of immunotherapy

Immune infiltration and immune function analyses were performed using the ssGSEA algorithm (https://bioconductor.org/packages/GSVA/). The ‘ESTIMATE’ software package (https://bioconductor.org/packages/estimate/) was used to compute immune and stromal scores, with each sample assigned a StromalScore, ImmuneScore and ESTIMATEScore. Box plots were generated to compare immune scores between high- and low-risk groups. To evaluate the potential efficacy of immunotherapy, the Tumor Immune Dysfunction and Exclusion (TIDE) score for each sample was calculated using the TIDE database (http://tide.dfci.harvard.edu/). A scatter plot illustrating the correlation between TIDE scores and risk groups was generated to assess their relationship.

Analysis of mutated landscapes

Mutation data for pancreatic cancer were obtained from the TCGA database using the ‘MAfTools’ package (https://bioconductor.org/packages/maftools), with ‘Mutect2’ selected as the mutation data type. Tumor mutation burden (TMB) for each sample was separately calculated for high- and low-risk groups using the TMB function within the ‘MAfTools’ package.

Cell lines and cell culture

The present study used paired pancreatic ductal models for primary screening: PANC-1 cells (representing classical PDAC) and HPDE-6 cells, an immortalized normal pancreatic ductal epithelial line. PANC-1 was selected due to its extensive use in PDAC research (>4,400 PubMed-indexed studies as of June 4, 2025) and its genomic profile, which harbors key driver mutations (KRAS and TP53) frequently identified in patient-derived PAAD cohorts. This mutational status was specifically documented for PANC-1 by Deer et al (43). The PANC-1 cell line (HLCL-023) and the HPDE-6 cell line (H1-3201) were purchased from iCell Bioscience, Inc., and OriCell Therapeutics. PANC-1 was cultured in DMEM (Provider: Hyclone,) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. HPDE-6 was maintained in the iCell Primary Epithelial Cell Culture System, supplemented with 5 ml Primary Epithelial Cell Culture Additive (iCell Bioscience, Inc.), 2% FBS and 5 ml penicillin/streptomycin. HPDE-6 cells and PANC-1 cells were removed from the liquid nitrogen tank and placed in a 37°C water bath with rapid shaking to thaw. After thawing, the cells were transferred to centrifuge tubes, and 2 ml of corresponding complete medium was added to each of them. The supernatant was removed by centrifugation, 2 ml of complete medium was re-added, gently pipetted and mixed, and the content was transferred to a six-well plate for incubation in a 37°C incubator. Cell passaging was performed when the cell density grew to 80–90%. The passaging ratio was 1:3. When the cell density of the passage reaches ~95%, reverse transcription-quantitative PCR (RT-qPCR) experiments were performed on each cell.

RNA extraction and RT-qPCR

Total RNA was extracted from HPDE-6 cells and PANC-1 cells using Trizol Plus (catalog number 15596026; Takara Bio, Inc.; iCell Bioscience, Inc.), and cDNA (reverse transcribed at 42°C for 15 min, then at 85°C for 5 min, and then placed on ice to cool) was synthesized using a first-strand cDNA synthesis kit (catalog number #K1622; Thero Fisher Scientific, Inc.). The RT-qPCR reaction mixture was prepared according to the standard protocol, consisting of 8.2 µl ddH2O, 10.0 µl SYBR Green Master Mix (catalog number 04913914001; Roche Diagnostics GmbH), 0.4 µl forward primer, 0.4 µl reverse primer and 1.0 µl cDNA template, making a final volume of 20.0 µl. The components were thoroughly mixed in a PCR reaction tube. Primers for RT-qPCR were designed and synthesized using Primer Premier 5.0 software (PREMIER Biosoft) and are listed in Table SI. β-actin was used as the internal reference. Amplification was performed under the following conditions: Initial denaturation at 95°C for 5 min, followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec. The instrument automatically analyzed the results. After PCR amplification, the qPCR system was calibrated using the negative control to set the threshold and baseline, determining the Ct value for each sample. The validity of the Ct value was confirmed through the melting curve analysis. The results were then exported, and differences in target gene expression between the control group and each experimental group were analyzed using the 2−ΔΔCq method (44).

Data analysis

Data were analyzed using R software (version 4.3.1; The R Foundation). The following packages were used for analysis: The limma package was used for variance analysis (Linear Models and Empirical Bayes); the survival package (https://cran.r-project.org/package=survival) was used for Kaplan-Meier curve plotting (survfit function) and cox regression analysis (coxph function); the GLMNet package (https://cran.r-project.org/package=glmnet) was used for regularization analysis (GLMNet function) in LASSO regression; pROC (https://cran.r-project.org/package=pROC) and timeROC (https://cran.r-project.org/package=timeROC) packages were used for ROC analysis and timeROC analysis (ROC function and timeROC function); the RMS package (https://cran.r-project.org/package=rms) was used for nomogram analysis and calibration curve analysis (CPH function and calibrate function); the GSVA package (https://bioconductor.org/packages/GSVA) for single sample Gene Set Enrichment Analysis immune infiltration and immune function analysis (gsva function); the estimate package was used to calculate tumor purity, the immune score and the stromal score (estimateScore function); and the MafTools package was used to construct a mutation landscape analysis (TMB function). Statistical analyses for validation experiments were performed using GraphPad Prism 9 software (Dotmatics). All quantitative data are expressed as mean ± standard error of the mean, based on results from ≥3 independent experiments. P<0.05 was considered to indicate a statistically significant difference.

Results

Analysis of differentially expressed mRNAs in PAAD and their functional enrichment

Based on the criteria of log2(foldchange) >1 or <-1, DEGs were identified by comparing PAAD tissues with normal tissues. A total of 244 upregulated and 170 downregulated DEGs were detected. Volcano and heat maps were then generated to visualize these findings (Fig. 3A and B).

Figure 3. PAAD gene expression analysis. (A) Heat map and (B) volcano plot of differentially expressed RNA in the TCGA dataset (Blue, downregulated expression; red, upregulated expression). (C) Upregulated gene GO enrichment analysis. (D) Upregulated gene KEGG enrichment analysis. (E) Downregulated gene GO enrichment analysis. (F) Downregulated gene KEGG enrichment analysis. PAAD, pancreatic adenocarcinoma; TCGA, The Cancer Genome Atlas; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change; FDR, false discovery rate; BP, Biological Process; CC, Cellular Component; MF, Molecular Function.

Subsequent GO and KEGG enrichment analyses were performed on both upregulated and downregulated DEGs to assess differences in functional pathways. Upregulated DEGs demonstrated significant enrichment in pathways characteristic of pancreatic cancer, especially those involving ECM remodeling. Regarding biological processes, these genes were primarily associated with ECM organization and extracellular structure organization. At the cellular component level, enrichment was observed in the collagen-containing ECM and basement membrane. Molecular function analysis indicated involvement in ECM structural constituent and integrin binding (Fig. 3C). Conversely, downregulated DEGs were associated with normal pancreatic physiological activities. Biological processes included digestion and the disruption of cellular anatomical structure in another organism. Cellular components were associated with platelet α granule and brush border membrane. Molecular functions involved serine-type peptidase activity and serine hydrolase activity (Fig. 3E). KEGG pathway analysis further revealed several key signaling pathways: Upregulated DEGs were enriched in ECM-receptor interaction and cytoskeletal regulation in muscle cells, whilst downregulated DEGs were enriched in pancreatic secretion and protein digestion and absorption pathways (Fig. 3D and F).

Construction and evaluation of a prognostic model for TCGA cohort

To further evaluate genes associated with prognosis, a univariate Cox analysis was performed on 414 DEGs in the TCGA dataset, identifying 18 genes with significant prognostic relevance. The results of the univariate Cox analysis for these 18 genes are presented in a forest plot (Fig. 4A). LASSO regression analysis was then performed (Fig. 4B and C), and multifactorial Cox regression analysis was used to calculate risk scores, resulting in a final prognostic model comprising 12 genes. Patients were categorized into high- and low-risk groups based on the median risk score. The results revealed that patients in the high-risk group had significantly worse prognoses than those in the low-risk group (P<0.05; Fig. 4D). Moreover, the ROC curves demonstrated that the area under the curve (AUC) values for predicting patient prognosis at 1, 3, and 5 years were 0.812, 0.832 and 0.858, respectively. These values exceeded those of clinical characteristics such as sex, age and disease stage (Fig. 4E and F).

Figure 4. Establishment of The Cancer Genome Atlas data set patient prediction model. (A) Univariate Cox analysis of differentially expressed genes. (B and C) Least absolute shrinkage and selection operator regression analysis: (B) Optimal λ selection and (C) Coefficient paths. (D) Kaplan-Meier curve analysis of high- and low-risk groups. (E) Riskscore clinical features timeROC curve. (F) Common clinical features ROC curve. ROC, receiver operating characteristic; AUC, area under the curve; T, tumor; N, node; M, metastasis.

Performance and clinical applicability of constructing predictive models

Given the association between risk score and patient prognosis, a Cox regression analysis was performed to assess the relationship between the risk score and other common clinical parameters (Fig. 5A). A nomogram was then constructed to estimate the 1-, 3- and 5-year overall survival (OS) of patients with PAAD (Fig. 5B). Additionally, the AUC values for several clinical factors used to predict 1-, 3- and 5-year OS were assessed. The results revealed that the nomogram exhibited good prognostic predictive ability, though its performance was not as strong as the riskScore clinical feature (Fig. 5C). Furthermore, calibration curve analysis demonstrated that the predicted values of the nomogram closely aligned with actual observed survival outcomes, indicating high accuracy (Fig. 5D).

Figure 5. Nomogram construction and performance evaluation. (A) Cox regression analysis of clinical features. (B) Construction of the nomogram. (C) Nomogram timeROC curve. (D) Nomogram calibration curve. ***P<0.001. ROC, receiver operating characteristic; T, tumor; N, node; AUC, area under the curve; OS, overall survival.

External validation and evaluation of forecasting models

Similar results were observed in the external validation cohort GSE28735. In this cohort, patients in the high-risk group had significantly worse prognoses than those in the low-risk group (Fig. 6A). Due to the absence of survival data beyond 3 years, ROC curves were analyzed for risk scores at 1, 2 and 3 years. The AUC for these time points were 0.845, 0.921 and 0.960, respectively (Fig. 6B). Univariate and multivariate Cox analyses confirmed that the risk score functioned as an independent prognostic factor (Fig. 6C).

Figure 6. Verification results of external data set GSE28735. (A) Kaplan-Meier curve analysis of high- and low-risk groups in the GSE28735 cohort. (B) TimeROC curve. (C) Risk Score: Univariate and multivariate Cox regression analyses. ***P<0.001. AUC, area under the curve.

Exploration of the immunological properties of PAAD

The immune microenvironment serves a crucial role in determining patient prognosis (45). In the present study, the low-risk group exhibited significantly elevated immune infiltration compared with the high-risk group in myeloid-derived suppressor cells (MDSCs), macrophages and regulatory T cells (Fig. 7A). Additionally, the T-cell co-stimulation immune function pathway demonstrated significantly increased activity in the low-risk group compared with in the high-risk group (Fig. 7B). Furthermore, immune-related scores, namely the ESTIMATE, immune and stromal scores, were significantly higher in the high-risk group than in the low-risk group (Fig. 7C).

Figure 7. Immunoassay of high- and low-risk groups. (A) ssGSEA immune infiltration analysis. (B) ssGSEA immune function analysis. (C) Immune scoring box plot. (D) TIDE box line diagram for the high- and low-risk groups. (E) TIDE and RiskScore scatter plot. *P<0.05; **P<0.01; ***P<0.001. ssGSEA, single sample Gene Set Enrichment Analysis; TIDE, Tumor Immune Dysfunction and Exclusion; MDSC, myeloid-derived suppressor cell; HLA, human leukocyte antigen; MHC, major histocompatibility complex; TME, tumor microenvironment; APC, antigen-presenting cell; CCR, C-C motif chemokine receptor.

The present study also evaluated differences in immunotherapy responses between the two groups (Fig. 7D). The results indicated that TIDE scores were significantly lower in the low-risk group compared with in the high-risk group, suggesting a lower likelihood of tumor immune escape and a greater potential benefit from immunotherapy. Moreover, TIDE analysis revealed a significant positive correlation between the TIDE score and the risk score (Fig. 7E). As the risk score increased, the TIDE score increased, indicating a higher likelihood of tumor immune escape and a reduced response to immunotherapy. These findings suggest that the risk score may serve as a predictive biomarker for ICI efficacy, aiding in the clinical selection of patients who are most likely to benefit from immunotherapy whilst minimizing overtreatment in non-responsive individuals.

Analysis of mutated landscapes

Gene mutations serve a notable role in tumor progression and patient outcomes (46). In the present study, the five most frequently mutated genes in both the high- and low-risk groups were KRAS, TP53, ‘mothers against decapentaplegic homolog 4’, cyclin dependent kinase inhibitor 2A and titin (TTN) (Fig. 8). These findings highlight the need for further investigation into targeted therapies against high-frequency mutations (such as KRAS and TP53) in conjunction with risk-stratification strategies to optimize precision treatment approaches.

Figure 8. Analysis of abrupt landscape changes. Mutation landscape of the (A) high-risk and (B) low-risk groups. TMB, tumor mutation burden.

Expression of prognostic genes in pancreatic cancer cell lines detected by RT-qPCR

To validate the aforementioned findings, RT-qPCR analysis was performed on the pancreatic cancer cell line, PANC-1, and the normal pancreatic cell line, HPDE-6. The results demonstrated significant differences in the expression levels of prognostic genes between the two cell lines (Fig. 9). Compared with the HPDE-6 group, the PANC-1 group demonstrated significantly upregulated expression of chromogranin B, dehydrogenase/reductase 9, endothelial cell specific molecule 1, fibronectin type III domain containing 1, glutaminase 2, 15-hydroxyprostaglandin dehydrogenase, matrix metallopeptidase 9 (MMP9), metallothionein 1M and TTN, whilst insulin growth factor-like family member 2, kynureninase and solute carrier family 16 member 10 were significantly downregulated.

Figure 9. Reverse transcription-quantitative PCR analysis of CHGB, DHRS9, ESM1, FNDC1, GLS2, HPGD, IGFL2, KYNU, MMP9, MT1M, SLC16A10 and TTN expression. *P<0.05; **P<0.01; ***P<0.001. CHGB, chromogranin B; DHRS9, dehydrogenase/reductase 9; ESM1, endothelial cell specific molecule 1; FNDC1, fibronectin type III domain containing 1; GLS2, glutaminase 2; HPGD, 15-hydroxyprostaglandin dehydrogenase; IGFL2, insulin growth factor-like family member 2; KYNU, kynureninase; MMP9, matrix metallopeptidase 9; MT1M, metallothionein 1M; SLC16A10, solute carrier family 16 member 10; TTN, titin.

Discussion

Pancreatic cancer is among the most aggressive and lethal malignancies (47), characterized by a high incidence, frequent recurrence, rapid metastasis and high mortality. It is also associated with a low early detection rate and a poor 5-year survival rate (48,49). According to estimates from the International Agency for Research on Cancer's Global Cancer Observatory, the global incidence of pancreatic cancer is projected to rise by 61.7% by 2040 (50). Given these challenges, early diagnosis is essential for improving treatment outcomes and reducing mortality. In the present study, 12 genes closely associated with the prognosis of patients with PAAD were identified using univariate Cox and LASSO regression analyses following variance analysis. A risk score system was then developed to assess patient prognosis. The findings revealed that a high-risk score was associated with immune competence and could serve as a prognostic indicator within the standard treatment framework for PAAD.

The TME has emerged as a promising immune target for several malignancies, including PAAD, due to its critical role in cancer progression and drug resistance (51). In the present study, a systematic evaluation of immune microenvironment characteristics between high- and low-risk groups revealed significantly elevated infiltration levels of MDSCs, macrophages and regulatory T cells in the low-risk group, along with simultaneous activation of the immune checkpoint and T-cell co-stimulatory pathways (Fig. 7A and B). The abnormal accumulation of MDSCs and Tregs in the high-risk group reflects a key mechanism of immune evasion. MDSCs deplete arginine in the TME through Arginase-1, directly impairing both the metabolic function and proliferation of CD8+ T cells (52). Tregs downregulate CD80/CD86 expression on antigen-presenting cells via cytotoxic T-lymphocyte associated protein 4 (CTLA-4)-dependent trogocytosis, weakening their co-stimulatory function whilst releasing free programmed death-ligand 1 (PD-L1) to further suppress effector T cells (53). At the same time, the elevated macrophage presence likely indicates a predominance of M2-polarized TAMs, which promote immunosuppression through IL-10/TGF-β secretion and ECM remodeling driven by MMP9 overexpression (54). Notably, the simultaneous activation of immune checkpoint pathways (such as PD-1/CTLA-4) and T-cell co-stimulatory pathways (such as inducible T cell costimulator/CD28) suggests an adaptive resistance mechanism: Although co-stimulatory signals may reflect sustained antigen exposure, they are preferentially diverted to expand immunosuppressive Treg populations rather than stimulate effector T cells. This leads to T-cell exhaustion and establishes a self-perpetuating immunosuppressive environment (55). This suggests that the TME in high-risk patients exhibits a complex interplay between immunosuppression and immune activation, likely driven by a dual mechanism in which tumors evade immune responses by recruiting suppressive immune cells whilst upregulating co-stimulatory signals (23). Furthermore, the concurrent increase in ESTIMATE, immune and stromal scores in the high-risk group (Fig. 7C) highlights the complexity of stromal-immune interactions within the TME. High stromal scores may indicate reinforcement of the physical barrier to drug delivery (16), whereas elevated immune scores reflect the presence of immunosuppressive cell infiltration. This paradoxical coexistence may explain the shortcomings of previous single-targeting strategies. For instance, selectively inhibiting CAFs could intensify immunosuppression by disregarding their role as a spatial barrier that restricts immune cell infiltration (18). Conversely, relying solely on ICIs may lead to ineffective drug penetration due to the dense ECM (19). Additionally, the present study identified a stronger immunotherapy response potential in the low-risk group, as indicated by the TIDE analysis (Fig. 7D). The positive correlation between the risk score and TIDE score (Fig. 7E) suggests that the risk score could serve as a potential biomarker for clinical stratification. This finding aligns with emerging research trends: Whilst single-targeted checkpoint blockade (such as PD-1/PD-L1 inhibition) has demonstrated limited efficacy in PAAD (24), multimodal strategies that combine immunotherapy with matrix remodeling [such as targeting hyaluronidase to improve drug penetration (26)] or metabolic interventions [such as regulating amino acid metabolism in MDSCs (32)] have shown synergistic potential. A preclinical study further reported that the simultaneous inhibition of C-X-C motif chemokine receptor 4 (which regulates immune cell chemotaxis) and PD-L1 markedly improves therapeutic response in pancreatic cancer models (34), suggesting that multi-targeted interventions could help overcome the cycle of immunosuppression and physical barriers within the TME.

Although the present study offers a potential theoretical foundation for improving early screening and individualized treatment strategies for pancreatic cancer, several limitations should be considered. First, the relatively small clinical sample size may limit the generalizability of the results, and the stability of the prognostic model needs to be tested in large-scale, multi-center prospective cohorts from diverse geographic regions. The preliminary in vitro validation, which relied mainly on a single cancer cell line, may also restrict the broader applicability of the results and should be expanded to include additional cellular models in future investigations. Second, the conclusions of the present study are primarily based on bioinformatics analysis of retrospective transcriptomic data, lacking multidimensional experimental validation. Future research should incorporate in vitro cell models, in vivo animal studies and prospective clinical cohorts to systematically assess the biological functions and clinical relevance of core genes in pancreatic cancer progression, metastasis and drug resistance. Additionally, the regulatory networks of key differentially expressed genes remain incompletely understood. Their upstream epigenetic modifications (such as DNA methylation or non-coding RNA regulation) and downstream effector pathways (such as molecular mechanisms related to EMT or metabolic reprogramming) require further exploration using high-throughput technologies such as CRISPR screening and chromatin conformation capture.

In conclusion, the present study highlights the covariation patterns of key TME features through multi-omics analysis, addressing the limitations of traditional single-phenotype studies. As demonstrated by the evolution of TME-targeted therapeutic strategies [from ‘matrix ablation’ to ‘immune modification’ to ‘metabolic remodeling’ (Fig. 1)], combined interventions targeting multiple dimensions of tumor biology have become an inevitable trend. Moving forward, based on the risk score model established in the present study, further screening of key molecular targets that regulate immunosuppression (such as MDSCs), matrix remodeling (such as ECM components secreted by CAFs) and metabolic reprogramming (such as lactate metabolism) could provide a theoretical framework for developing integrated ‘immune-matrix-metabolism’ therapeutic strategies.

Furthermore, the present study developed and validated a prognostic model for pancreatic cancer based on 12 key DEGs, demonstrating strong predictive power for patient survival in both the TCGA and GSE28735 cohorts, independent of traditional clinical parameters. Patients in the high-risk group exhibited immunosuppressive microenvironmental characteristics, such as increased MDSC infiltration and elevated TIDE scores, suggesting lower sensitivity to immunotherapy. The differential expression patterns of core genes in pancreatic cancer cell lines were confirmed through RT-qPCR experiments. This model presents a potential tool for prognostic stratification and predicting immunotherapy response in patients with PAAD. However, further refinement incorporating multi-omics data is needed to optimize its clinical applicability.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by the Henan Science and Technology Research Project (grant nos. 242102310490, 242102111077 and 252102310194) and the Excellent Young Teachers Training Program of Sanquan College of Xinxiang Medical University (grant no. SQ2023YQJH02).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

ZLZ and XJL confirm the authenticity of all the raw data. This paper was written by ZLZ, who also performed data acquisition and data analysis. XJL was responsible for the study design, revision of the entire manuscript and interpretation of the data. JLZ, FYL, YLZ and HW contributed their expertise in the field and participated in the bioinformatics analysis. BW, XML and RL performed the experiments, whilst JWSM, JYZ and YSX performed the analysis of experimental data. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References