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.

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.

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.

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.

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.

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.

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.

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 $\times$g 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.

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).

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).

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.

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).

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.

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).

Fig. 1.

Overall survival based on CAF scores and WGCNA of CAF-related genes. (A) Kaplan–Meier survival curves for patients with PAAD from GSE57495 (n = 63) stratified by CAF scores (MCPcounter) (p = 0.010, log-rank test). (B) Kaplan–Meier survival curves for patients with PAAD from GSE57495 (n = 63) stratified by CAF scores (TIDE) (p = 0.010, log-rank test). (C) Gene dendrogram and modules before merging in the GSE57495 dataset. (D) Correlation analysis between MEcyan module eigengenes and CAF score (left, r = 0.86, p 0.001) and StromalScore (right, r = 0.79, p 0.001) in GSE57495 (n = 63). The color scale represents the correlation coefficient. p-values were calculated for Pearson correlation coefficients. (E) Scatterplot of module membership (MM) vs. gene significance (GS) for the cyan module in GSE57495 (n = 63). The Pearson correlation coefficient (r = 0.89) and p-value (p 0.001) were calculated. CAF, cancer-associated fibroblast; WGCNA, weighted gene co-expression network analysis; PAAD, particularly pancreatic adenocarcinoma; TIDE, Tumor Immune Dysfunction and Exclusion.

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).

Fig. 2.

Key protein–protein interaction (PPI) cluster in the Cyan module identified using the MNC algorithm. A PPI network was constructed using 501 genes from the cyan module, with a minimum confidence score of 0.4. Isolated nodes were removed, and the MNC algorithm was applied to identify the top 50 hub genes. Red indicates core proteins with high connectivity.

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).

Fig. 3.

Establishment and identification of a risk signature. (A) Univariate Cox regression analysis of 50 candidate genes identified using the MNC algorithm. Genes with p 0.05 were selected for further analyses (COL28A1, p = 0.017; TGFB2, p = 0.012; CDH2, p = 0.041; MMP14, p = 0.045; TGFBI, p = 0.025; PLOD2, p = 0.003; COL22A1, p = 0.028; RUNX2, p = 0.039). (B,C) LASSO regression analysis with 10-fold cross-validation. (B) LASSO coefficient profiles for eight genes. (C) Ten-fold cross-validation for tuning parameter ($\lambda$) selection. (D,E) Kaplan–Meier analysis of overall survival in the training cohort (GSE57495, n = 63, p 0.05) (D) and external validation cohort (GSE62452, n = 65, p 0.05) (E) p-values were calculated using the log-rank test. MNC, Maximum Neighborhood Component; LASSO, Least Absolute Shrinkage and Selection Operator.

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).

Fig. 4.

Expression levels of model genes and CAF marker genes. (A) Correlations between the CAF score, stromal score, and risk score in GSE57495 (n = 63). Spearman’s correlation coefficients were used. (B) Correlation analysis between model genes and CAF marker genes in GSE57495 (n = 63). Spearman’s correlation was used, with * indicating p 0.05 (e.g., PLOD2 vs. COL1A1: r = 0.71, p 0.001; TGFBI vs. COL3A1: r = 0.68, p 0.001). Significance levels are denoted as follows: *p 0.05, **p 0.01, ***p 0.001. (C) Heatmap showing the expression levels of the five model genes and CAF marker genes in the high-risk (n = 32) and low-risk (n = 31) groups in GSE57495 (n = 63). The risk groups were defined by the median risk score. Differential expression for all genes between risk groups was significant (p 0.01, Student’s t-test).

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.

Fig. 5.

TIDE analysis and drug sensitivity analysis using risk models. (A) Proportion of samples responding to immunotherapy in the high-risk (n = 31) and low-risk (n = 32) groups in GSE57495 (n = 63) (p = 0.004, Chi-square test). (B) TIDE scores in the high-risk (n = 31) and low-risk (n = 32) groups (p = 0.0035, Student’s t-test). (C–H) Sensitivity analysis of drugs in the risk model. Box plots show the estimated IC50 values for AGI-5198 (C, p 0.001), dihydrorotenone (D, p 0.001), elephantin (E, p 0.001), irinotecan (F, p 0.001), palbociclib (G, p 0.001), and sorafenib (H, p 0.001) in the high-risk and low-risk groups. p-values were calculated using the Student’s t-test. Significance levels are denoted as follows: ***p 0.001.

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).

Fig. 6.

Multidimensional expression validation. (A) Heat map of candidate gene expression in CAFs (n = 47) and PAAD cell lines (n = 37) from CCLE. (B) Box plot of candidate gene expression in CAFs and PAAD cell lines. p-values were calculated by the Wilcoxon test (COL28A1, p 0.001; TGFB2, p 0.05; TGFBI, p 0.001; PLOD2, p 0.001; COL22A1, p = 0.71). (C) UMAP plot showing the major cell types in single-cell RNA sequencing data. (D,E) Differential distribution of PLOD2 (D) and TGFBI (E) expression at the single-cell level. (F,G) Bar charts show the average expression level and percentage of positive cells (pctExpress) in each cell population. (F) Western blot analysis of PLOD2 and TGFBI protein levels in CAFs and pancreatic cancer cell lines (CAPAN-1 and SW1990). Group names are indicated below the lanes (H). Bar graphs show the densitometric analysis results for PLOD2 (I) and TGFBI (J) protein levels from three independent experiments. Data are presented as means $\pm$ SD. **p 0.01, ***p 0.001, ****p 0.0001 by one-way ANOVA with Tukey’s post hoc test (CAFs vs. CAPAN-1 for PLOD2: p = 0.0003; CAFs vs. SW1990 for PLOD2: p 0.0001; CAFs vs. CAPAN-1 for TGFBI: p = 0.0012; CAFs vs. SW1990 for TGFBI: p = 0.0008). CCLE, Cancer Cell Line Encyclopedia.

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.

Fig. 7.

In vitro analysis of PLOD2. (A) Western blot showing the PLOD2 knockout efficiency in CAFs. Two sgRNAs (sg1-PLOD2 and sg2-PLOD2) were used. (B) Schematic diagram of the co-culture experiment. CAFs were cultured in serum-free DMEM for 48 hours, and the conditioned medium (CM) was collected to treat pancreatic cancer cells. (C,D) Colony formation assay of CAPAN-1 (n = 3) (C) and SW1990 (D) cells treated with CM from control CAFs (CON-CM) or PLOD2-knockout CAFs (sgPLOD2-CM). NC, negative control (fresh serum-free DMEM). *p 0.05, **p 0.01 The bar graph (E) shows the quantitative analysis of colony numbers from three independent experiments. Data are presented as means $\pm$ SD. *p 0.05, ***p 0.001 by one-way ANOVA with Tukey’s post hoc tests (for CAPAN-1: NC vs. CON-CM, p = 0.009; CON-CM vs. sg1-PLOD2-CM, p = 0.025; CON-CM vs. sg2-PLOD2-CM, p = 0.017). (F–I) Wound healing assay of CAPAN-1 (n = 3) (F) and SW1990 (H) cells treated with CM. The migration rate was calculated after 48 hours. Data are from three independent experiments and presented as means $\pm$ SD. *p 0.05, **p 0.01 by one-way ANOVA with Tukey’s post hoc tests (for CAPAN-1: NC vs. CON-CM, p = 0.012; CON-CM vs. sg2-PLOD2-CM, p = 0.031). Scale bars = 200 μm. (J,K) PLOD2 knockout impairs CAF-mediated collagen contraction (n = 3). (J) Representative images of collagen gels after 48 hours of contraction. Control CAFs (CON) show significant gel contraction, while PLOD2-knockout CAFs (sg1-PLOD2-CAF, sg2-PLOD2-CAF) exhibit reduced contractility. The red circles indicate the boundary of the gels, highlighting the extent of contraction. (K) Quantitative analysis of gel contraction. Data are presented as mean $\pm$ SD from three independent experiments. *p 0.05, ***p 0.001 vs CON group by one-way ANOVA with Tukey’s post hoc test.

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.

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., $\alpha$v$\beta$3 and $\alpha$v$\beta$5), 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.

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.

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.