Cancer-Associated Fibroblast-Centric Risk Model Predicts Immunotherapy Resistance in Pancreatic Cancer and Reveals PLOD2 as a Key Stromal Therapeutic Target

Abstract

Background:

Pancreatic neoplasms, particularly pancreatic adenocarcinoma (PAAD), are aggressive malignancies marked by extensive infiltration of cancer-associated fibroblasts (CAFs) and a highly complex tumor immune microenvironment. These pathological features are strongly associated with poor patient survival. However, the precise mechanisms underlying the role of CAFs in PAAD have not been determined.

Methods:

To systematically analyze the functions of CAFs in PAAD and their associations patient outcomes, an integrative approach combining multi-omics data with experimental validation was used.

Results:

Integrated weighted gene co-expression network and protein–protein interaction network analyses revealed CAF-related genes with functional significance. Experimental verification was performed to examine the influence of candidate CAF-related genes identified using multi-database analyses on tumor cell behavior. COL28A1, TGFB2, TGFBI, PLOD2, and COL22A1 were core genes closely associated with CAFs. Patients in the high-risk group demonstrated a higher immune escape ability and lower predictive response rate to immunotherapy than those in the low-risk group. Several potential targeted therapeutic compounds were identified, including dihydrorotenone and sorafenib. Single-cell RNA sequencing and expression analyses confirmed elevated expression of TGFBI and PLOD2 in CAFs. Functional experiments demonstrated that PLOD2 promotes tumor progression by regulating extracellular matrix remodeling.

Conclusions:

This study provides insights into the molecular mechanisms underlying PAAD and establishes a theoretical foundation for the development of CAF-targeting therapeutic strategies.

Graphical abstract

Keywords

pancreatic neoplasms / cancer-associated fibroblast / tumor microenvironment / extracellular matrix

1. Introduction

Pancreatic adenocarcinoma (PAAD), a highly malignant cancer, is responsible for approximately 300,000 deaths globally each year [1]. It is characterized by strong invasiveness and a high propensity for metastasis. Despite available treatment modalities, such as surgery, radiotherapy, chemotherapy, and immunotherapy, the overall therapeutic efficacy remains suboptimal, with a 5-year survival rate below 10% [2]. The primary determinant of the poor prognosis of PAAD is the complex tumor immune microenvironment (TIME), which involves various components, including stromal cells, immune cells, and a dense extracellular matrix (ECM) [3, 4]. The ECM leads to vascular collapse, increased interstitial fluid pressure, and obstruction of lymphatic return, thereby hindering effective drug delivery, deep penetration, and immune cell infiltration [5, 6]. Cancer-associated fibroblasts (CAFs), key cellular components of the ECM, are closely linked to tumor malignancy [7]. Consequently, an in-depth investigation of the functions of CAFs and their associated genes in PAAD holds promise for the development of targeted PAAD therapies.

CAFs arise from the activation and transformation of normal fibroblasts within the TIME by various stimuli, such as transforming growth factor-\beta (TGF-\beta), Hedgehog (Hh), nuclear factor kappa-B (NF-\kappaB) [8, 9]. High \alpha-smooth muscle actin expression is a canonical CAF marker. CAFs regulate the biological behavior of tumors through multiple mechanisms, including promoting tumor cell growth, survival, invasion, angiogenesis, and immune evasion [10]. By modulating the ECM composition and remodeling via paracrine signaling, CAFs induce epithelial-mesenchymal transition (EMT) and increase tumor cell migration and proliferation [11]. CAFs secrete various cytokines and chemokines, including growth differentiation factor 15, CXC chemokine ligand 12, vascular endothelial growth factor, fibroblast growth factor 5, and interleukin 6, which play crucial roles in shaping an immunosuppressive microenvironment, promoting tumor immune escape, and enhancing the invasiveness of pancreatic cancer cells [12, 13, 14].

However, the key pathways by which CAFs contribute to PAAD and their impact on patient survival remain unclear. This research integrates advanced computational frameworks—such as weighted gene co-expression network analysis (WGCNA), protein–protein interaction (PPI) networks, and single-cell RNA sequencing (scRNA-seq)—with targeted functional experiments to deconstruct the role of CAFs in PAAD progression. We hypothesize that a distinct gene expression signature in CAFs drives ECM remodeling and immune suppression in PAAD and that targeting these molecular mediators could overcome immunotherapy resistance. Our findings illuminate the molecular mechanisms underlying PAAD and provide a basis for the development of actionable clinical tools for improving immunotherapeutic interventions and tailoring drug regimens based on individual drug sensitivity profiles.

2. Materials and Methods

2.1 Data Download and Processing

The microarray datasets GSE57495 (GPL15048 plat-form) and GSE62452 (GPL6244), involving 63 and 69 (total n = 69; 4 samples excluded due to missing clinical data) PAAD samples with corresponding clinical information, respectively, were retrieved from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). Each dataset was processed and normalized independently. Data were normalized using the “normalize.quantiles” function of the “preprocessCore” package (version 1.68.0; Bioconductor, Fred Hutchinson Cancer Center, Seattle, WA, USA) in R software (version 4.4.2; R Foundation for Statistical Computing, Vienna, Austria). Gene symbols were assigned to probe sets based on manufacturer annotation files. Ambiguous mappings (where one probe corresponded to multiple genes) were removed. For genes detected by multiple probes, their expression values were averaged to obtain a single representative value per gene. When the two datasets were integrated for a combined analysis (e.g., during model construction and validation), we applied the removeBatchEffect function from the “limma” R package (version 3.62.2; Bioconductor, Fred Hutchinson Cancer Center, Seattle, WA, USA) to adjust for potential technical batch effects arising from different platforms.

2.2 CAF Scoring and WGCNA

CAF scores were calculated via the MCPcounter.probesets and MCPcounter.genes functions of the “MCPcounter” R package (version 1.2.0; developed by the Tumor Microenvironment and Cancer Immunology research team, INSERM U1138, Cordeliers Research Centre, Paris, France), based on a fibroblast-specific gene signature (e.g., COL1A1, COL3A1, COL6A1, COL6A2, DCN, PAMR1, TAGLN and GREM1) [15]. To identify gene expression modules linked to CAF abundance, a Weighted Gene Co-expression Network Analysis was executed using the WGCNA package (version 3.5.2; Peter Langfelder and Steve Horvath, University of California, Los Angeles, CA, USA) [16]. In this method, a co-expression network is constructed based on pairwise correlations between gene profiles. Genes are then grouped into discrete modules according to a topological overlap measure [17]. During this process, an automatically selected soft-thresholding power was used to ensure a scale-free network architecture. Finally, relationships between each module and the CAF score were assessed by calculating Pearson correlation coefficients, identifying key regulatory modules.

2.3 Construction and Validation of a CAF Signature

To pinpoint genes central to CAFs, we initially built a PPI network using the STRING database (version 12.0, https://string-db.org). After filtering interactions with a minimum confidence score of 0.4 and eliminating isolated nodes, we employed the CytoHubba plugin within Cytoscape (version 3.8.2; The Cytoscape Consortium, University of California, San Diego, CA, USA). Specifically, by applying the Maximum Neighborhood Component (MNC) algorithm, we identified the top 50 most influential genes [18]. These candidate genes were then subjected to a survival analysis via Cox proportional hazards regression and variable selection using Least Absolute Shrinkage and Selection Operator (LASSO) regression in R. The optimal penalty parameter (\lambda) was determined through 10-fold cross-validation, selecting the value that minimized the mean cross-validated error. Based on coefficients from the LASSO model, we developed a risk score formula: Risk Score = \mathrm{\Sigma } (Coefficient \times Gene Expression Level). For external validation, the formula was applied to the GSE62452 dataset using the coefficients derived from the training cohort, without any retraining or parameter adjustment. Subsequently, patients were dichotomized into high- and low-risk groups using the median risk score. Finally, survival differences between these groups were assessed with the Kaplan–Meier method, and the prognostic accuracy of the final model was evaluated through calibration plots.

2.4 Tumor Immune Dysfunction and Exclusion (TIDE) Prediction

The TIDE score assesses how well a tumor can evade the immune system, helping predict the effectiveness of immunotherapy in different patient risk groups. High-scoring tumors exhibit robust immune evasion, whereas low scores indicate a heightened likelihood of responding positively to immune checkpoint inhibitor (ICI) therapy. As PAAD is not a distinct category within the TIDE database, expression data were submitted with the tumor type designated as “Other” and the “prior immunotherapy” option set to “No” to exclude the confounding effects of prior treatments.

2.5 Drug Sensitivity Analysis

To investigate whether distinct risk strata exhibit differential drug responses, we utilized computational models trained on the Genomics of Drug Sensitivity in Cancer (GDSC) repository to forecast therapeutic sensitivity. GDSC cell line expression data served as the training set, whereas samples from the constructed risk model served as the test set. Batch effect correction was applied using the “eb” (Empirical Bayes) method to control for systematic bias. The prediction parameters were set as variation cut-off = 0.2 and minimum sample size = 10 to ensure model robustness. Sensitivity prediction scores for each test sample were obtained for various drugs.

2.6 Cancer Cell Line Encyclopedia (CCLE) Analysis

To confirm distinct transcriptional activity of preselected candidate genes between tumor cells and fibroblasts, RNA-seq datasets were acquired from the CCLE repository (https://depmap.org/portal/ccle). After curating focused expression data for our genes of interest from available PAAD cell lines and fibroblast samples, we performed a comparative analysis. This was executed in the R computing environment, employing the moderated t-test in the limma framework (version 3.62.2; Bioconductor), alongside the ggpubr add-on (version 0.6.0; available via the Comprehensive R Archive Network) for data presentation. Given potential deviations from a normal distribution, we used rank-based, non-parametric tests. p-values were subsequently adjusted using established false discovery rate (FDR) control procedures to ensure statistical reliability.

2.7 Single-Cell RNA Sequencing Analysis

The scRNA-seq data were analyzed using the SCAR_Atlas_0799 dataset and the ScarAtlas online platform (http://www.scaratlas.com/). Cells were initially filtered based on platform-preset quality control criteria (number of detected genes: 200–6000; mitochondrial gene percentage: 20%). An automated workflow was then used for dimensionality reduction and clustering, including a principal component analysis, UMAP visualization, and cell clustering using the Louvain algorithm. Cell-type annotation was performed by actively comparing expression profiles against established marker genes within the platform. For PLOD2 and TGFBI expression features, the visualization tools of the platform were utilized directly; their spatial distributions were mapped onto a clustered UMAP plot and bar charts were generated to display the average expression level and frequency of positive cells (pctExpress) within each cell population.

2.8 Tumor Cell Line and CAF Culture

PAAD cell lines, CAPAN-1, SW1990 and Miapaca-2, were obtained from Zhang DX Laboratory (Hefei, China). Cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Gibco, Thermo Fisher Scientific, Waltham, MA, USA) containing 10% fetal bovine serum at 37 °C in a humidified atmosphere comprising 5% CO₂. CAFs (Nanjing, China, Cat# JY-J883) were purchased from Jinyuan Biotechnology Co. (https://www.ssrcc.com.cn/). A short tandem repeat (STR) analysis was leveraged to validate the purity and ensure cell line stability. To control for potential off-target effects of CRISPR/Cas9, two independent sgRNAs (sg1-PLOD2 and sg2-PLOD2) targeting distinct sequences of the PLOD2 gene were used for knockout. Detailed information on the plasmid vectors and the precise sgRNA sequences (sg1-PLOD2 and sg2-PLOD2) are provided in Supplementary Table 1. To prepare conditioned media from CAFs (CAF-CM) and CAFs with PLOD2 knockout (sgCAF-CM), CAFs or sgCAFs were cultured in serum-free medium. Cells in the logarithmic growth phase were harvested and resuspended to a concentration of 10⁶ cells/mL. Twenty milliliters of the cell suspension were seeded into culture flasks. Upon reaching 80% confluence, the supernatant was collected, centrifuged at 300 \timesg for 20 min to remove cellular debris, and stored at –20 °C until use. All sequences were provided by Tsingke Biotech (Beijing, China) and are detailed in Supplementary Table 1. All cell lines were validated by STR profiling and tested negative for mycoplasma.

2.9 Western Blotting

Cells were lysed using RIPA buffer (Beyotime Biotechnology, Shanghai, China, Cat# P0013B). Total proteins (10 µg per lane) were separated by SDS-PAGE and then transferred onto PVDF membranes. The membranes were blocked with 5% non-fat milk in TBST for 1 hour at room temperature and then incubated overnight at 4 °C with primary antibodies against TGFBI (rabbit polyclonal, ZEN BIO, Chengdu, China, Cat# 380025; 1:2000 dilution) and PLOD2 (mouse monoclonal, Proteintech, Wuhan, China, Cat# 66342-1-Ig, 1:2000 dilution). After washing, the membranes were incubated with HRP-conjugated secondary antibodies, including goat anti-rabbit IgG (Abbkine, Redlands, CA, USA, Cat# A21020; 1:5000) and goat anti-mouse IgG (Abbkine, Redlands, CA, USA, Cat# A21010; 1:5000), for 2 hours at room temperature. Protein bands were visualized using an ECL detection kit (Vazyme, Nanjing, China, Cat# E412-01). The chemiluminescent signals were captured and quantified using ImageJ software (version 1.52a; National Institutes of Health, Bethesda, MD, USA).

2.10 Colony Formation Assay

The capacity for clonogenic survival was assessed by seeding approximately 1000 treated cells into each well of a 6-well plate. Cultures were maintained at 37 °C in a humidified atmosphere containing 5% CO₂ for approximately 14 days or until distinct colonies became evident. Colony visualization was achieved by fixing with 4% paraformaldehyde (Beyotime Biotechnology, Shanghai, China, Cat# P0099) for 15 min and staining with 0.1% crystal violet (Solarbio, Beijing, China, Cat# C8470) for 20 minutes. Upon air-drying, digital images of the stained colonies were captured, and the total number of colonies per well was determined using ImageJ software (version 1.52a; National Institutes of Health, Bethesda, MD, USA).

2.11 Wound Healing Assay

Starting with a density of 3 \times 10⁵ cells per well, cultures were established in 6-well plates and grown to complete monolayer confluency. At this stage, a uniform scratch was created across the center of each well using a sterile 10 µL pipette tip. Detached cells were then removed by aspirating the medium and performing two successive washes with phosphate-buffered saline. Fresh, serum-free DMEM was added to each well. Digital images of the wound were captured at the initial time point (0 hours) and after a 48-hour incubation period. Cell migration was quantified by calculating the percentage decrease in the scratch width over time, using the following equation: % Migration = [(Scratch Width_0⁢h – Scratch Width_48⁢h)/Scratch Width_0⁢h] \times 100.

2.12 Collagen Contraction Assay

The contractile capacity of CAFs was evaluated using a collagen-based assay kit (Cell Biolabs, Inc., San Diego, CA, USA, Cat#CBA-201) following the manufacturer’s two-step model. Control CAFs (CON) and PLOD2-knockout CAFs (sg1-PLOD2, sg2-PLOD2) were prepared. Cells were mixed with cold Collagen Gel Working Solution to a final density of 0.5–1 \times 10⁶ cells/mL. Then, 0.5 mL of the mixture was added to each well of a 24-well plate for polymerization (1 hour, 37 °C). After adding 1 mL of medium, gels were incubated for 48 hours. Gels were then released from the well walls to initiate contraction. Images were taken at various time points. The contraction index was measured as the reduction in gel area using ImageJ software (version 1.52a; National Institutes of Health, Bethesda, MD, USA).

2.13 Statistical Analysis

Statistical analyses were performed using R (version 4.4.2; R Foundation for Statistical Computing, Vienna, Austria) and GraphPad Prism (version 9.0.0; GraphPad Software, San Diego, CA, USA). Data are presented as means \pm SD, and p 0.05 was considered significant. The Shapiro–Wilk test was used to assess normality. For two-group comparisons, the unpaired t-test (normal data) or Mann–Whitney U test (non-normal data) was used. Multiple groups were compared using one-way ANOVA with Tukey’s test (normal) or Kruskal–Wallis with Dunn’s test (non-normal). For survival analyses, the Kaplan–Meier method with the log-rank test was used. Correlations were analyzed with Pearson (normal) or Spearman (non-normal) methods. The prognostic performance of the risk model was assessed using a time-dependent ROC analysis. Univariate and multivariate Cox regression analyses were used to evaluate its independence from and performance relative to common clinical variables. For multiple comparisons, the FDR was controlled using the Benjamini–Hochberg method. Exact statistical tests, p-values, and sample sizes are provided in the figures and legends.

3. Results

3.1 CAF Scoring, Survival Analysis, and Identification of Key CAF-Related Genes

A survival analysis based on CAF MCPcounter scores revealed that patients with higher CAF scores exhibited significantly worse overall survival (OS) than that of patients with lower CAF scores (p 0.05) in the GSE57495 dataset (Fig. 1A). Similar results were obtained using the CAF TIDE score (Fig. 1B). WGCNA partitioned all genes into 16 distinct modules (Fig. 1C). Given the typical enrichment of CAFs in the tumor stroma, we analyzed the correlations between the eigengene of each module and CAF and stromal scores. The Cyan module exhibited a strong positive correlation with CAF (r = 0.86, p 0.001) and stromal scores (r = 0.79, p 0.001) (Fig. 1D). Within the Cyan module, a similarly strong positive correlation with the CAF score was observed (r = 0.89, p 0.05), as illustrated in the MM-GS scatter plot (Fig. 1E).

This description is inaccurate and should be replaced with: A PPI network was constructed using 501 genes from the Cyan module. After removing 185 isolated nodes that lacked interactions with other proteins, a connected network was obtained for subsequent analysis.(Supplementary Fig. 1) Subsequently, the MNC algorithm was applied, yielding a network cluster containing 50 key genes, including LOXL2, MMP3, PLOD2, and P4HA3, involved in stromal remodeling (Fig. 2).

3.2 Univariate Cox Regression, LASSO Regression Analysis of Candidate Genes, and Prognostic Risk Model

Univariate Cox regression analyses of the 50 candidate genes identified using the MNC algorithm were performed. COL28A1, TGFB2, CDH2, MMP14, TGFBI, pLOD2, COL22A1, and RUNX2 demonstrated significant associations with overall survival (p 0.05). All other genes with p > 0.05 were excluded from further analyses (Fig. 3A). To eliminate collinearity among candidate genes and prevent overfitting of the prognostic model, a LASSO regression analysis of the eight prognostically significant genes was performed (Fig. 3B,C). The risk coefficient was calculated for each gene, and five genes (COL28A1, TGFB2, TGFBI, pLOD2, and COL22A1) were ultimately selected for model construction and risk scoring. A survival analysis demonstrated that patients in the high-risk group had significantly worse OS than that of patients in the low-risk group (p 0.05) (Fig. 3D). This indicates that the constructed CAF-associated risk model is closely related to CAF levels within the TIME and effectively predicts the prognosis of patients with PAAD. Our risk model demonstrated strong generalizability, as evidenced by its performance using the external validation set (GSE62452). In alignment with the results for the training cohort, the model accurately identified a high-risk group exhibiting significantly worse survival rates than those of the low-risk group (Fig. 3E), thereby substantiating its reliability and clinical utility across different patient populations. A time-dependent ROC analysis demonstrated the predictive accuracy of the model for 1-, 3-, and 5-year survival in both the training (AUC = 0.691, 0.725, 0.612) and validation cohorts (AUC = 0.649, 0.569, 0.706) (Supplementary Fig. 2A,B). To evaluate the independent prognostic value of the risk score, we performed Cox regression analyses using both cohorts. In the training cohort, the risk score was a strong and independent prognostic factor in both univariate (HR = 4.047, 95% CI: 2.040–8.027, p 0.001) and multivariate analyses (HR = 3.994, 95% CI: 1.855–8.599, p 0.001), outperforming clinical variables, such as stage, CA19-9, and CA125 (Supplementary Fig. 2C,D). In the independent validation cohort, the risk score showed a trend in the univariate analysis (HR = 1.051, 95% CI: 0.888–1.245, p = 0.563) but did not retain statistical significance in multivariate analysis when adjusted for sex, CA19-9, and CA125 (HR = 1.059, 95% CI: 0.849–1.319, p = 0.613; Supplementary Fig. 2E,F).

3.3 Expression Differences and Correlations Between Model Genes and CAF Marker Genes

The putative CAF marker genes served as a foundation, informed by prior investigations in the field. For visualization, a comprehensive heatmap was employed to juxtapose the expression dynamics of the selected genes with those of recognized CAF markers within the final risk estimator. Spearman’s method demonstrated significant, positive correlations (p 0.05) between the CAF-associated risk score and metrics representing CAF infiltration and stromal abundance (Fig. 4A). Among the five genes incorporated into the model, four (COL22A1, TGFB2, TGFBI, and PLOD2) were positively correlated with CAF marker expression, whereas COL28A1 showed a negative correlation (Fig. 4B). Additionally, CAF marker genes sourced from the literature showed higher expression levels in the high-risk group than in the low-risk group. Among the model genes, COL22A1, TGFB2, TGFBI, and pLOD2 were expressed at high levels in the high-risk group, whereas COL28A1 was expressed at high levels in the low-risk group (Fig. 4C).

3.4 Association of Risk Stratification Models With Differential Immunotherapy Response and Chemotherapy Sensitivity Among Risk Subgroups

To evaluate the value of the CAF-associated risk model in predicting the immunotherapy response, samples from the GSE57495 datasets were analyzed using the TIDE platform. The TIDE scores quantified the immune escape potential and predicted the sensitivity to ICIs for each sample. The immunotherapy response differed significantly among risk groups (p 0.05); the high-risk group contained more non-responders, while the low-risk group had more responders (Fig. 5A). TIDE scores were significantly higher in the high-risk group than in the low-risk group (p 0.05) (Fig. 5B), indicating a stronger immune escape capacity and lower immunotherapy response in high-risk samples. To identify potential therapeutic agents for PAAD, 21 drugs exhibiting differential sensitivities (p 0.05) between risk groups were identified. AGI-5198, dihydrorotenone, elephantin, irinotecan, palbociclib, and sorafenib showed higher sensitivity in the high-risk group than in the low-risk group (Fig. 5C–H), suggesting that these drugs are effective against CAF-associated PAAD.

3.5 Validation Using Cell Lines and scRNA-seq

To determine whether the CAF-related genes included in the signature originated primarily from CAFs, multidimensional validation was performed, including cell line and scRNA-seq analyses. RNA-Seq data were collected from 47 fibroblast and 37 PAAD cell lines. As determined using the “limma” package (Fig. 6A) and Wilcoxon tests (Fig. 6B), COL28A1 expression levels were significantly lower and TGFB2, TGFBI, pLOD2, and COL22A1 levels were significantly higher in fibroblasts than in PAAD cell lines. The two genes (TGFBI and PLOD2) showing the most significant differences were selected for validation. The scRNA-seq analysis successfully resolved four major cell populations (Fig. 6C): CD8+ T cells, fibroblasts, malignant epithelial cells, and monocytes/macrophages. The fibroblast population exhibited a prominent fibrotic gene expression signature. Specifically, the collagen cross-linking enzyme gene PLOD2 showed a higher pctExpress in fibroblasts than in other cell types (Fig. 6D,F). Similarly, the TGF-\beta pathway effector gene TGFBI exhibited a higher pctExpress in fibroblasts than in malignant epithelial cells (Fig. 6E,G). Western blotting confirmed that protein levels of PLOD2 and TGFBI are significantly higher in CAFs than in pancreatic cancer cell lines (Fig. 6H–J).

3.6 Knockout of PLOD2 in CAFs Attenuates Their Pro-Proliferative Effect on Pancreatic Cancer Cells

To validate the effect of CAF-derived PLOD2 on pancreatic cancer cells, sgRNAs targeting PLOD2 (sg1-PLOD2 and sg2-PLOD2) were transduced into CAFs for PLOD2 gene knockout. Fig. 7A shows the transfection efficiency. CAFs were cultured in serum-free DMEM for 48 hours, and the collected supernatant (CM) was co-cultured with pancreatic cancer cells (Fig. 7B). Colony formation assays demonstrated that the CM from control CAFs (CON-CM) resulted in significantly greater clonogenic survival and proliferation of CAPAN-1 and SW1990 pancreatic cancer cells than those for the negative control (NC) medium (p 0.05). However, CM from CAFs with PLOD2 knockout (sgPLOD2-CM) attenuated this pro-tumorigenic effect significantly; the colony formation efficiency of both Capan-1 and SW1990 cells treated with sgPLOD2-CM was significantly lower than that of cells treated with CON-CM (p 0.05) (Fig. 7C–E). Wound healing assays showed that CON-CM promoted pancreatic cancer cell migration. In Capan-1 cells, the wound healing percentage was significantly higher in the CON-CM group than in the NC group (p 0.05); a similar increase was observed in SW1990 cells (p 0.05). CM from CAFs with PLOD2 knockout (sg1-PLOD2-CM and sg2-PLOD2-CM) maintained a pro-migratory effect compared with results for the NC in Capan-1 (p 0.05 for sg2-PLOD2-CM) and SW1990 (p 0.05 for sg1-PLOD2-CM) cells (Fig. 7F–I). A collagen contraction assay revealed that PLOD2 is critical for CAF-mediated matrix remodeling. Control CAFs (CON) significantly contracted collagen gels, while PLOD2 knockout (sg1- and sg2-PLOD2-CM) severely impaired this contractile function.Assessment of collagen gel contraction demonstrated that gels in the NC group had the largest diameter. The diameter was significantly reduced in the CON-CM group, and it was significantly larger in both the sg1-PLOD2-CM and sg2-PLOD2-CM groups compared to the CON-CM group (p 0.01) (Fig. 7J,K). These results confirm the critical role of PLOD2 in CAF-mediated ECM remodeling.

4. Discussion

Pancreatic cancer, a malignancy with a dismal prognosis, is characterized by a TIME composed of dense ECM and extensive inflammatory cell infiltration. CAF enrichment within the TIME, which promotes tumor progression, is linked to poor patient outcomes [19]. The prognostic value of CAFs has been established in some cancers, such as hepatocellular carcinoma and ovarian cancer [20], highlighting their potential as therapeutic targets in PAAD.

4.1 Prognostic Significance of the CAF-Centric Risk Model and Its Association With the Immunotherapy Response

Our results demonstrated that a high CAF-associated risk score and stromal score are correlated with worse OS in patients with PAAD and predict lower response rates to immunotherapy. This is likely attributable to the dense desmoplastic stroma and abundant immune cell infiltration typically accompanying high CAF scores, which collectively foster immune escape [21]. Although immunotherapy has achieved remarkable success in the treatment of various malignancies, patients with pancreatic cancer generally exhibit poor immune responses, partly because the dense ECM acts as a physical barrier [22]. As a crucial cellular component within the TIME, CAFs contribute to angiogenesis as well as tumor cell infiltration and metastasis [23]. Recent research underscores the key role of CAFs in shaping the immune landscape of the TIME by modulating the ECM-mediated anchoring and trafficking of immune cells and suppressing immune activation [24]. Yu et al. [25] revealed that CAFs promote tumor immune escape by affecting immune checkpoint activity through high CD73 expression. We stratified patients into high- and low-risk groups based on CAF-associated risk scores, revealing that the high-risk group possessed stronger immune escape potential and lower immunotherapy efficiency, confirming that CAFs affect immune evasion through multiple pathways.

4.2 Key CAF-Related Genes and Their Functional Implications

We identified a panel of potential therapeutic agents targeting CAFs. For instance, dihydrorotenone, an inhibitor of mitochondrial electron transport chain complex I, effectively inactivates CAFs, thereby slowing the progression of transplanted gastric cancer in mice [26]. In triple-negative breast cancer, Zhang et al. [27] showed that dasatinib increases the permeability and therapeutic sensitivity of anti-PD-1 agents by modulating CAFs. These studies offer new perspectives for targeted PAAD therapy; however, the clinical efficacy of these agents requires further validation.

Through multidimensional analyses, we systematically elucidated the expression characteristics and functional mechanisms of key genes in PAAD CAFs. Using WGCNA, PPI analysis, and Cox/LASSO regression algorithms, we identified a signature based on five key genes (COL28A1, TGFB2, TGFBI, PLOD2, and COL22A1). Expression levels of TGFB2, TGFBI, PLOD2, and COL22A1 were significantly positively correlated with levels of established CAF markers, whereas COL28A1 showed a negative correlation. PLOD2 and TGFBI were expressed at particularly high levels in CAFs, providing novel insight into the molecular basis of CAF regulation in the PAAD microenvironment.

Among collagen family members, COL28A1 and COL22A1 exhibited distinct expression patterns. Although COL28A1 expression was relatively low, its positive correlation with CAF markers suggests an indirect regulatory role in TME remodeling [28]. TGFBI was overexpressed in CAFs. As an ECM protein downstream of the TGF-\beta signaling pathway, TGFBI activates key signaling pathways (e.g., FAK/PI3K/Akt) by binding integrin receptors (e.g., \alphav\beta3 and \alphav\beta5), thereby promoting tumor cell adhesion, migration, and invasion [29, 30, 31]. Multiple studies have linked high TGFBI expression to the invasive and metastatic properties of various malignancies [32, 33], potentially involving EMT regulation and remodeling of the TIME [34]. The expression profile of PLOD2 was equally striking. As a key enzyme for collagen cross-linking, PLOD2 promotes the formation of stable collagen fibrils, creating a fibrotic microenvironment conducive to tumor progression [35]. PLOD2 expression was also higher in CAFs than in tumor cells. Therefore, PLOD2 has emerged as a potential novel therapeutic target for PAAD.

Based on our CAF-associated risk stratification, we identified several potential therapeutic agents targeting CAFs. Dihydrorotenone (DHR), a mitochondrial complex I inhibitor, effectively inactivates CAFs and delays transplanted gastric cancer progression in mice [26]. This effect may be mediated by metabolic reprogramming, wherein DHR-induced mitochondrial dysfunction promotes CAF reversion to a quiescent state [36]. Similarly, dasatinib has CAF-targeting potential via ECM remodeling rather than direct cytotoxicity [27]. By reducing ECM production and collagen cross-linking, dasatinib alleviates matrix stiffness, enhances drug penetration, and promotes cytotoxic T lymphocyte infiltration into tumor cores [37]. These findings provide a mechanistic rationale for targeting CAFs to reprogram the fibrotic tumor microenvironment and overcome immunotherapy resistance in PAAD.

4.3 Experimental Validation of PLOD2 Function and Therapeutic Potential

In functional experiments, the CM increased the proliferation and invasion of pancreatic cancer cells significantly. Crucially, this pro-tumorigenic effect was significantly attenuated upon PLOD2 knockout in CAFs, indicating a pivotal role for PLOD2 in CAF-mediated tumor progression. PLOD2, a member of the lysyl hydroxylase family, catalyzes the hydroxylation of lysine residues in collagen, promoting collagen fibril cross-linking and stabilization [38]. Our results provide experimental support that inhibiting PLOD2 in CAFs can disrupt their role in cancer cell aggressiveness. PLOD2 knockout reduces collagen cross-linking, altering the physical properties of the ECM and, consequently, diminishing its tumor-promoting effects, consistent with previous findings [39]. High PLOD2 expression enhances ECM rigidity and density and activates mechanosignaling pathways (e.g., YAP/TAZ) via increased matrix tension, promoting tumor cell infiltration and metastasis [40, 41]. Thus, PLOD2 promotes tumor progression through multiple mechanisms (ECM remodeling, mechanosignaling activation, and secretory regulation) and its knockout effectively inhibits these effects. Future studies should further investigate the downstream molecular networks.

4.4 Limitations and Future Directions

This study clarifies the role of CAFs in pancreatic cancer; however, several limitations warrant acknowledgment. The use of public datasets with limited sample sizes may affect the reliability of the findings. Furthermore, the reliance on cell-based experiments, without in vivo confirmation, represented a constraint. The mechanisms by which PLOD2 remodels collagen remain unclear, and the efficacy and safety of promising drugs should be verified in living organisms. Future work should include larger patient cohorts, develop animal models for PLOD2 targeting, and advance promising treatments into clinical trials.

5. Conclusion

In conclusion, our study establishes a CAF-centric risk model that predicts patient survival and immunotherapy resistance in pancreatic cancer and identifies PLOD2 as a key stromal therapeutic target. The functional link between CAF-derived PLOD2 and tumor progression underscores the value of targeting the stromal compartment. This work provides a framework for developing stroma-directed combination therapies and emphasizes the need for further research focused on PLOD2 inhibition to overcome treatment resistance in PAAD.

Availability of Data and Materials

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

References

[1] Hu ZI, O’Reilly EM. Therapeutic developments in pancreatic cancer. Nature Reviews. Gastroenterology & Hepatology. 2024; 21: 7–24. https://doi.org/10.1038/s41575-023-00840-w.

[2] Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA: A Cancer Journal for Clinicians. 2024; 74: 12–49. https://doi.org/10.3322/caac.21820.

[3] N Kalimuthu S, Wilson GW, Grant RC, Seto M, O’Kane G, Vajpeyi R, et al. Morphological classification of pancreatic ductal adenocarcinoma that predicts molecular subtypes and correlates with clinical outcome. Gut. 2020; 69: 317–328. https://doi.org/10.1136/gutjnl-2019-318217.

[4] Hwang WL, Jagadeesh KA, Guo JA, Hoffman HI, Yadollahpour P, Reeves JW, et al. Single-nucleus and spatial transcriptome profiling of pancreatic cancer identifies multicellular dynamics associated with neoadjuvant treatment. Nature Genetics. 2022; 54: 1178–1191. https://doi.org/10.1038/s41588-022-01134-8.

[5] Liu X, Shao Y, Li Y, Chen Z, Shi T, Tong Q, et al. Extensive Review of Nanomedicine Strategies Targeting the Tumor Microenvironment in PDAC. International Journal of Nanomedicine. 2025; 20: 3379–3406. https://doi.org/10.2147/IJN.S504503.

[6] Ho WJ, Jaffee EM, Zheng L. The tumour microenvironment in pancreatic cancer - clinical challenges and opportunities. Nature Reviews. Clinical Oncology. 2020; 17: 527–540. https://doi.org/10.1038/s41571-020-0363-5.

[7] Katarkar A, Bottoni G, Clocchiatti A, Goruppi S, Bordignon P, Lazzaroni F, et al. NOTCH1 gene amplification promotes expansion of Cancer Associated Fibroblast populations in human skin. Nature Communications. 2020; 11: 5126. https://doi.org/10.1038/s41467-020-18919-2.

[8] Mucciolo G, Araos Henriquez J, Jihad M, Pinto Teles S, Manansala JS, Li W, et al. Egfr-activated myofibroblasts promote metastasis of pancreatic cancer. Cancer Cell. 2024; 42: 101–118. https://doi.org/10.1016/j.ccell.2023.12.002.

[9] Fang Z, Meng Q, Xu J, Wang W, Zhang B, Liu J, et al. Signaling pathways in cancer-associated fibroblasts: recent advances and future perspectives. Cancer Communications (London, England). 2023; 43: 3–41. https://doi.org/10.1002/cac2.12392.

[10] Lunardi S, Muschel RJ, Brunner TB. The stromal compartments in pancreatic cancer: are there any therapeutic targets? Cancer Letters. 2014; 343: 147–155. https://doi.org/10.1016/j.canlet.2013.09.039.

[11] Tew BY, Legendre C, Gooden GC, Johnson KN, Martinez RA, Kiefer J, et al. Isolation and characterization of patient-derived CNS metastasis-associated stromal cell lines. Oncogene. 2019; 38: 4002–4014. https://doi.org/10.1038/s41388-019-0680-2.

[12] Elkabets M, Gifford AM, Scheel C, Nilsson B, Reinhardt F, Bray MA, et al. Human tumors instigate granulin-expressing hematopoietic cells that promote malignancy by activating stromal fibroblasts in mice. The Journal of Clinical Investigation. 2011; 121: 784–799. https://doi.org/10.1172/JCI43757.

[13] Cheteh EH, Sarne V, Ceder S, Bianchi J, Augsten M, Rundqvist H, et al. Interleukin-6 derived from cancer-associated fibroblasts attenuates the p53 response to doxorubicin in prostate cancer cells. Cell Death Discovery. 2020; 6: 42. https://doi.org/10.1038/s41420-020-0272-5.

[14] Bulle A, Lim KH. Beyond just a tight fortress: contribution of stroma to epithelial-mesenchymal transition in pancreatic cancer. Signal Transduction and Targeted Therapy. 2020; 5: 249. https://doi.org/10.1038/s41392-020-00341-1.

[15] Han C, Liu T, Yin R. Biomarkers for cancer-associated fibroblasts. Biomarker Research. 2020; 8: 64. https://doi.org/10.1186/s40364-020-00245-w.

[16] Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9: 559. https://doi.org/10.1186/1471-2105-9-559.

[17] Tran TM, Guha R, Portugal S, Skinner J, Ongoiba A, Bhardwaj J, et al. A Molecular Signature in Blood Reveals a Role for p53 in Regulating Malaria-Induced Inflammation. Immunity. 2019; 51: 750–765.e10. https://doi.org/10.1016/j.immuni.2019.08.009.

[18] Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Systems Biology. 2014; 8: S11. https://doi.org/10.1186/1752-0509-8-S4-S11.

[19] Takamatsu H, Yamamoto KI, Tomonobu N, Murata H, Inoue Y, Yamauchi A, et al. Extracellular S100A11 Plays a Critical Role in Spread of the Fibroblast Population in Pancreatic Cancers. Oncology Research. 2019; 27: 713–727. https://doi.org/10.3727/096504018X15433161908259.

[20] Song L, Li Q, Lu Y, Feng X, Yang R, Wang S. Cancer Progression Mediated by CAFs Relating to HCC and Identification of Genetic Characteristics Influencing Prognosis. Journal of Oncology. 2022; 2022: 2495361. https://doi.org/10.1155/2022/2495361.

[21] Timperi E, Croizer H, Khantakova D, Rana M, Molgora M, Guerriero JL, et al. At the Interface of Tumor-Associated Macrophages and Fibroblasts: Immune-Suppressive Networks and Emerging Exploitable Targets. Clinical Cancer Research: an Official Journal of the American Association for Cancer Research. 2024; 30: 5242–5251. https://doi.org/10.1158/1078-0432.CCR-24-1690.

[22] Xu X, Wu Y, Qian X, Wang Y, Wang J, Li J, et al. Nanomedicine Strategies to Circumvent Intratumor Extracellular Matrix Barriers for Cancer Therapy. Advanced Healthcare Materials. 2022; 11: e2101428. https://doi.org/10.1002/adhm.202101428.

[23] Peng C, Xu Y, Wu J, Wu D, Zhou L, Xia X. TME-Related Biomimetic Strategies Against Cancer. International Journal of Nanomedicine. 2024; 19: 109–135. https://doi.org/10.2147/IJN.S441135.

[24] Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nature Reviews. Immunology. 2015; 15: 669–682. https://doi.org/10.1038/nri3902.

[25] Yu M, Guo G, Huang L, Deng L, Chang CS, Achyut BR, et al. CD73 on cancer-associated fibroblasts enhanced by the A_2⁢B-mediated feedforward circuit enforces an immune checkpoint. Nature Communications. 2020; 11: 515. https://doi.org/10.1038/s41467-019-14060-x.

[26] Lee E, Yeo SY, Lee KW, Lee JA, Kim KK, Kim SH. New screening system using Twist1 promoter activity identifies dihydrorotenone as a potent drug targeting cancer-associated fibroblasts. Scientific Reports. 2020; 10: 7058. https://doi.org/10.1038/s41598-020-63996-4.

[27] Zhang Y, Zhou J, Wang Y, Wu Y, Li Y, Wang B, et al. Stimuli-responsive polymer-dasatinib prodrug to reprogram cancer-associated fibroblasts for boosted immunotherapy. Journal of Controlled Release: Official Journal of the Controlled Release Society. 2025; 381: 113606. https://doi.org/10.1016/j.jconrel.2025.113606.

[28] Guo Y, Lu G, Mao H, Zhou S, Tong X, Wu J, et al. miR-133b Suppresses Invasion and Migration of Gastric Cancer Cells via the COL1A1/TGF-β Axis. OncoTargets and Therapy. 2020; 13: 7985–7995. https://doi.org/10.2147/OTT.S249667.

[29] Shelton L, Rada JAS. Inhibition of human scleral fibroblast cell attachment to collagen type I by TGFBIp. Investigative Ophthalmology & Visual Science. 2009; 50: 3542–3552. https://doi.org/10.1167/iovs.09-3460.

[30] Zhan S, Bai X, Zhao Y, Tuoheti K, Yisha Z, Zuo Y, et al. TGFBI promotes proliferation and epithelial-mesenchymal transition in renal cell carcinoma through PI3K/AKT/mTOR/HIF-1α pathway. Cancer Cell International. 2024; 24: 265. https://doi.org/10.1186/s12935-024-03454-7.

[31] Wang H, Xu YH, Guo Y. Novel prognostic marker TGFBI affects the migration and invasion function of ovarian cancer cells and activates the integrin αvβ3-PI3K-Akt signaling pathway. Journal of Ovarian Research. 2024; 17: 50. https://doi.org/10.1186/s13048-024-01377-5.

[32] Sato T, Muramatsu T, Tanabe M, Inazawa J. Identification and characterization of transforming growth factor beta-induced in circulating tumor cell subline from pancreatic cancer cell line. Cancer Science. 2018; 109: 3623–3633. https://doi.org/10.1111/cas.13783.

[33] Huang H, Tang Q, Li S, Qin Y, Zhu G. TGFBI: A novel therapeutic target for cancer. International Immunopharmacology. 2024; 134: 112180. https://doi.org/10.1016/j.intimp.2024.112180.

[34] Lecker LSM, Berlato C, Maniati E, Delaine-Smith R, Pearce OMT, Heath O, et al. TGFBI Production by Macrophages Contributes to an Immunosuppressive Microenvironment in Ovarian Cancer. Cancer Research. 2021; 81: 5706–5719. https://doi.org/10.1158/0008-5472.CAN-21-0536.

[35] Li SS, Lian YF, Huang YL, Huang YH, Xiao J. Overexpressing PLOD family genes predict poor prognosis in gastric cancer. Journal of Cancer. 2020; 11: 121–131. https://doi.org/10.7150/jca.35763.

[36] Folmes CDL, Terzic A. Energy metabolism in the acquisition and maintenance of stemness. Seminars in Cell & Developmental Biology. 2016; 52: 68–75. https://doi.org/10.1016/j.semcdb.2016.02.010.

[37] Liu L, Zhang B, Wu X, Cheng G, Han X, Xin X, et al. Bioresponsive nanocomplex integrating cancer-associated fibroblast deactivation and immunogenic chemotherapy for rebuilding immune-excluded tumors. Nanomedicine: Nanotechnology, Biology, and Medicine. 2024; 58: 102743. https://doi.org/10.1016/j.nano.2024.102743.

[38] Tong Y, Qi Y, Xiong G, Li J, Scott TL, Chen J, et al. The PLOD2/succinate axis regulates the epithelial-mesenchymal plasticity and cancer cell stemness. Proceedings of the National Academy of Sciences of the United States of America. 2023; 120: e2214942120. https://doi.org/10.1073/pnas.2214942120.

[39] Qi Y, Xu R. Roles of PLODs in Collagen Synthesis and Cancer Progression. Frontiers in Cell and Developmental Biology. 2018; 6: 66. https://doi.org/10.3389/fcell.2018.00066.

[40] Lewis DM, Pruitt H, Jain N, Ciccaglione M, McCaffery JM, Xia Z, et al. A Feedback Loop between Hypoxia and Matrix Stress Relaxation Increases Oxygen-Axis Migration and Metastasis in Sarcoma. Cancer Research. 2019; 79: 1981–1995. https://doi.org/10.1158/0008-5472.CAN-18-1984.

[41] Yin M, Weng Y, Qi T. PLOD2 exacerbates cervical squamous cell carcinoma by suppressing p53 by binding to YAP1. Molecular Medicine Reports. 2025; 31: 23. https://doi.org/10.3892/mmr.2024.13388.

Funding

Anhui Province Science and Technology Innovation Tackle Plan Project(202423k09020009)

Graduate Research and Practice Innovation Project of Anhui Medical University(YJS20240022)

National Natural Science Foundation of China(82073372)