The raw drug response data of the Cancer Therapeutics Response Portal (CTRP v2) were obtained from the CTD² Data Portal (https://ocg.cancer.gov/programs/ctd2/data-portal/). The RNA-seq data of corresponding pancreatic cancer cell lines are available from the DepMap portal (https://depmap.org/portal/). The “TCGAbiolinks” R package (9) was used to download and process the bulk RNA-seq results of 178 PC patients along with clinical, mutation, and copy number variation (CNV) files in the TCGA database. The TCGA cohort was used as the training cohort to develop a predictive gene signature and nomogram. Meanwhile, a total of 186 patients with complete follow-up information from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database were merged as the validation cohort, including GSE28735, GSE62452, GSE85916 and GSE102238 datasets.

As there are no data regarding normal pancreatic cell lines in the DepMap portal, the transcriptome data of normal pancreatic duct epithelial cell lines (HPNE and HPDE) were obtained from GSE97003 and GSE181625. We also explored the expression levels of key genes regulating gemcitabine resistance between normal and tumor samples at both bulk and single-cell resolutions. Two scRNA-seq datasets were included in our study, GSE212996 and CRA00160 (10). The sequencing data of normal pancreatic tissues from the GTEx portal (https://gtexportal.org/home/) were analyzed in combination with the TCGA cohort. We also compared the expression levels of these genes in matched tumor-normal samples from GSE28735 (45 pairs), GSE102238 (28 pairs), and GSE101448 (18 pairs).

Raw RNA-seq and microarray data of bulk tissues were normalized using transcripts per kilobase million (TPM) transformation and robust multichip average (RMA) normalization, respectively, followed by log2 (x + 1) conversion (11, 12). The “Seraut” package was used to create an S4 object from the unique molecular identifier (UMI) count and barcode matrices of scRNA sequencing in accordance with the standard pipeline (13). After excluding low-quality cells (< 200 genes/cell, > 3000 genes/cell, < 3 cells/gene, hemoglobin genes < 3%, and > 10% mitochondrial genes), linear principal component analysis (PCA) and nonlinear uniform manifold approximation and projection (UMAP) algorithm were adopted to visualize the clustering results. The primary cell type annotations of the “singleR” package (14) were checked and adjusted manually based on marker genes from the CellMarker database (http://xteam.xbio.top/CellMarker/).

A lower gemcitabine resistance score indicates a higher drug sensitivity. According to the scaled gemcitabine resistance score, pancreatic cells were subtyped into gemcitabine-sensitive (Z score ≤ –0.5), intermediate gemcitabine-resistant (–0.5< Z score < 0.5), and gemcitabine-resistant (0.5 ≤ Z score) groups. Subsequently, the “limma” package was used to identify differentially expressed genes (DEGs) between gemcitabine-sensitive and gemcitabine-resistant cells (15). The DEGs with a P value < 0.01 and log₂ fold change (FC) > 1 were defined as the key regulators contributing to gemcitabine resistance.

The least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was applied to develop a risk model from the prognostic genes associated with chemoresistance in the TCGA cohort (16). The optimal penalization parameter was selected by tenfold cross validation. Gene signature with better performance and fewer included genes was selected as the appropriated model from multiple runs. The risk score of each enrolled patient was equal to the sum of each gene expression value multiplied by the corresponding coefficient. In line with other studies, the median risk score was selected as the threshold to equally divide the patients in the training cohort into two groups (17). Utilizing the “survival”, “survminer”, “timeROC”, and “stats” R packages, the prognostic values in the TCGA, GEO, and whole cohort were evaluated by Kaplan-Meier (K-M) curves with log-rank tests, time-dependent receiver operating characteristic (tROC) curves, and principal component analyses (PCA).

Univariate and multivariate analyses were conducted to identify independent prognostic factors with statistical significance for constructing a predictive nomogram (18). The performance of the nomogram was assessed by the concordance index (C index), tROC, decision curve analysis (DCA), and calibration curve. The following R packages were used for analyses and visualization: “survival”, “rms”, “pec”, and “ggDCA”.

Comparison of functional enrichment scores between the high- and low-risk groups was implemented by the “GSVA” and “limma” packages. Gene sets for GSVA analyses are open-access from the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gse), including Hallmark (c2.all.v2022.1.Hs.symbols) and KEGG (c2.KEGG.v2022.1.Hs.symbols) collections. TMB in the TCGA cohort was calculated as the number of mutated bases per million bases utilizing the “TCGAbiolinks” package.

The “calcPhenotype” function of “oncoPredict” R package was used to predict the drug sensitivity for each patient based on bulk RNA-seq result (19). The datasets for fitting the ridge regression model are updated and uploaded by the “oncoPredict” team, including resources from the Genomics of Drug Sensitivity in Cancer (GDSC2) and CTRP V2 databases (https://osf.io/c6tfx/).

The deconvolution module of the “IOBR” R package (20) was used to estimate the landscape of the TME, which integrates six open-source deconvolution methodologies, namely, CIBERSORT (21), MCPcounter (22), EPIC (23), xCell (24), quantiseq (25), and TIMER (26). In addition, the ssGSEA algorithm was also adopted to evaluate the TME score for each patient (27). The reference gene sets for ssGSEA estimation were gathered from previous works (28–32).

First, the ICB resistance score was calculated using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (32). However, the TIDE signature may be biased, as only tumor-infiltrating cell signatures are incorporated for estimation. Generally, the immune response (immune cells, intracellular networks, and intercellular networks) and tumor antigenicity (TMB) are two hallmarks involved in regulating the response to immunotherapy (33). Then, the Estimate Systems Immune Response (EaSIeR) method was introduced to predict ICB response, taking into account the tumor microenvironment as a whole (34). A lower TIDE value and a higher EaSIeR score indicate a better ICB response.

Three human pancreatic cancer cell lines (AsPC-1, CFPAC-1, and PANC-1) were purchased from Procell Life Science & Technology Co., Ltd. (Wuhan, China). Human pancreatic normal cell lines (hTERT-HPNE) were obtained from Meisen Cell Biotechnology Co., Ltd. (Hangzhou, China). AsPC-1 and CFPAC-1 cells were maintained in Roswell Park Memorial Institute-1640 and Iscove’s Modified Dulbecco’s Medium (IMDM; Procell), respectively. PANC-1 and HPNE cells were cultivated in high-glucose Dulbecco’s modified Eagle’s medium (DMEM; Procell). All media contained 10% fetal bovine serum (FBS; Procell) and 100 μg/mL penicillin-streptomycin solution (Procell). The humidity incubator was set to 37°C with 5% CO2.

Cells in six-well plates were transiently transfected with siRNA using Lipofectamine™ 3000 transfection reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. siRNAs were synthesized by GenePharma (Shanghai, China), including ALDH3B1 siRNA (sequence: 5′-GCU GAA GCC AUC GGA GAU UAG tt-3′), NCEH1 siRNA (sequence: 5′-CAA UGA UCG UUA ACA AUC Att-3′), and universal negative control (NC) siRNA.

Total RNA was isolated using RNA Isolater Total RNA Extraction Reagent (Vazyme, Nanjing, China) and dissolved in RNase-free ddH₂O (Vazyme). Reverse transcription and quantitative real-time PCR (qRT-PCR) procedures were performed using Hiscript III Reverse Transcriptase (Vazyme) and ChamQ Universal SYBR qPCR Master Mix (Vazyme). The primers were synthesized by Sangon (Shanghai, China), and the sequences were as follows: ALDH3B1, 5′-GCC CTG GAA CTA TCC GCT G-3′ (forward), 5′-CGT TCT TGC TAA TCT CCG ATG G-3′ (reverse); NCEH1, 5′-GAA TAC AGG CTA GTT CCA AAG-3′ (forward), 5′-TAC TTC TGT AAG ACT TCT GGC-3′ (reverse); GAPDH, 5′-GGA GCG AGA TCC CTC CAA AAT-3′ (forward), 5′-GGC TGT TGT CAT ACT TCT CAT GG-3′ (reverse). The whole experimental process was carried out as described previously by our group (35). Gene expression was compared using the semiquantitative 2-^ΔΔCt method.

Western blotting (WB) was performed according to the previous protocol (36). In brief, total protein was extracted using RIPA buffer (Sevenbio, Beijing, China), denatured for 10 min in a dry bath, and then quantified using BCA kit (Sevenbio).A total of 20ug protein was loaded for 10% SDS-PAGE electrophoresis and transferred to PVDF membrane (Millipore, Hertfordshire, UK), which was then blocked by 5% non-fat milk. After overnight incubation with primary antibodies at 4°C, the PVDF membrane was incubated with secondary antibodies for 2 h at room temperature. Signals were determined using ECL kit (Sevenbio). Primary antibodies against ALDH3B1 (1:500, Zen Bioscience, Chengdu, China), NCEH1 (1:500, Zen Bioscience), and β-actin (1:10000, ProteinTech Group, Rosemont, IL, USA) were used. Gray value analysis of all bands was performed using Image J software (Image J 1.53, NIH).

Cell Counting Kit-8 (CCK-8; GlpBio, CA, USA) was used to measure cell viability and proliferation. The experimental procedure was the same as the official instructions. For the proliferation assay, transfected cells (2000 cells per well) were seeded in 96-well plates, and optical density (OD) 450 nm values were measured at 0, 24, 48, 72, 96, and 120 hours with a microplate reader (BioTek, Vermont, USA). Cell proliferation fold was calculated using the following method: cell proliferation fold = (OD_{0-120 hours} – OD_blank)/(OD_{0 hours} – OD_blank). For the cytotoxicity assay, 4000 cells per well were plated in 96-well plates 24 hours before gemcitabine (GlpBio) treatment. After 48 hours, OD values were determined, and cell viability was obtained using the following algorithm: cell viability = (OD_{gemcitabine-treated group} – OD_blank)/(OD_{control group} – OD_blank).

All statistical analyses and graphical visualizations were implemented using R software (version 4.21) and packages from the Comprehensive R Archive Network (CRAN) or Bioconductor repositories. Two-tailed Student’s t test, nonparametric Wilcoxon, and Spearman tests were used for comparison and correlation analyses, as appropriate. Two-tailed P < 0.05 was considered significant.

The workflow of the current study is presented in Supplementary Figure S1 . To better compare the differences in gemcitabine sensitivity, data from the CTRP database were first Z-scored. As similar responses were observed in seven cell lines (PANC0403, PATU8988S, PANC0813, PANC1, CAPAN1, MIAPACA2, and SUIT2), 0.5 and –0.5 were selected as the grouping thresholds. Of the 30 PC cell lines, 13 were classified as gemcitabine-resistant, 7 as intermediate gemcitabine-resistant, and 10 as gemcitabine-sensitive ( Figure 1A and Supplementary Figure S2 ). Differential expression analysis revealed that the expression levels of 17 genes were significantly higher in the gemcitabine-resistant group than in the gemcitabine-sensitive group ( Figures 1B, C ). Additionally, we summarized the CNV and somatic mutation frequencies in the TCGA cohort ( Figures 1D, E ). The number of CNV gains was much higher than CNV losses in OSBPL1A, NCEH1, DUSP1, and ALDH3B1, while the opposite tendencies were observed in ERAP2, CHST11, GALM, and MSLN. The mutation rates of the DEGs related to gemcitabine chemoresistance were less than 3%, except for OSBPL1A (13%).

K-M analysis revealed that 6 out of 17 DEGs associated with gemcitabine chemoresistance had prognostic value in the TCGA cohort, namely, ALDH3B1, CHST11, EGFR, ERAP2, MSLN, and NCEH1 ( Figure 2A ). The expression levels of the six prognostic DEGs were significantly higher in the TCGA tumor samples than in the nonmatched GTEx normal samples ( Figure 2B ). Consistently, these genes were more highly expressed in gemcitabine-resistant PC cell lines than in normal pancreatic cell lines ( Figure 2C ). However, matched tumor-normal comparisons indicated that EGFR was not relatively upregulated in malignant tissues ( Figures 2D-F ).

After quality control, the number of normal/tumor cells in GSE212966 was 10751/17338 and in CRA00160 was 13771/34974. At the single-cell level, tumor cells expressed low levels of CHST11 compared with ALDH3B1, EGFR, ERAP2, MSLN, and NCEH1 ( Figures 3A-E , 4A-E ). The fibroinflammatory stroma and cancer-associated fibroblast (CAF)-enriched TME are two characteristics of PC (37). Pancreatic stellate cells are the precursors of CAFs, which are essential in tumorigenesis and gemcitabine resistance (37, 38). As shown in the UMAP plots, most of the stellate and fibroblast cells came from the tumor samples ( Figures 3A-C , 4A-C ). The relatively low expression of EGFR in bulk RNA-seq analyses was further confirmed by single-cell analyses ( Figures 3D , 4D ). However, it was still upregulated in malignant and stromal components of tumor tissues, proving that it was indispensable for tumor progression ( Figures 3E , 4E ). In contrast to the results in CRA00160, the mean UMI count of ERAP2 in GSE212966 was lower in tumor samples, which could be accounted for by the low proportion of detected malignant cells ( Figures 3A-D ).

Except for CHST11, the other five chemoresistance-related genes were included in the LASSO Cox model predicting the OS of the TCGA training cohort ( Figures 2G, H ). The risk score formula was as follows: risk score = (0.3731 * EGFR expression) + (0.2364 * ERAP2 expression) + (0.1588 * MSLN expression) + (0.1254 * NCEH1 expression) + (0.2233 * ALDH3B1 expression). On the basis of the median risk score of the patients in the TCGA cohort, all the enrolled patients were separated into high- or low-risk groups. K-M analyses and risk plots indicated that patients with higher risk scores had significantly worse prognoses ( Figures 5A, B ). Log-rank test P values of the TCGA, GEO, and whole cohorts were less than 0.001. The prediction accuracies of the 1-, 3-, and 5-year survival rates were 0.750/0.781/0.741 in the TCGA cohort, 0.596/0.615/0.755 in the GEO cohort, and 0.653/0.668/0.765 in the whole cohort, respectively ( Figure 5C ). Positive correlations were observed between the risk score and expression levels of the five genes ( Figure 5D ). Patients in different risk groups could be clearly distinguished in PCA reduction diagrams ( Figure 5E ).

In the TCGA cohort, the chemoresistance-related risk score was the only independent prognostic factor identified by univariate and multivariate Cox analyses ( Figures 6A, B ). In addition, the risk score was also a significant risk factor in the GEO and whole cohorts ( Figure 6A ). As the clinicopathological information of the GEO patients was incomplete, multivariate Cox regression was not applicable. Therefore, a nomogram predicting 1-, 3-, and 5-year OS rates was built based on the expression levels of the five genes constituting the risk model ( Figure 7A ). C index analysis indicated that the nomogram had the highest prediction performance compared with other factors in the TCGA cohort, which was further supported by tROC and DCA curves ( Figures 7B-D ). Then, the accuracies of the 1-, 3-, and 5-year OS predictions were validated internally and externally, and the deviations between the nomogram-predicted and observed OS were small ( Figure 7E ).

As the five prognostic genes (ALDH3B1, EGFR, ERAP2, MSLN, and NCEH1) were highly expressed in gemcitabine-resistant PC cell lines, this established gene signature was, therefore, a very promising indicator of gemcitabine sensitivity for PC patients. This hypothesis was supported by subsequent analyses, which showed that the risk score increased with the gemcitabine resistance score, and the low-risk group was more sensitive to gemcitabine ( Figures 8A-D ). TMB is a genetic characteristic of cancers. In the TCGA cohort, the risk score was positively correlated with TMB ( Figure 8E ). The top five mutated genes in the high- and low-risk groups were KRAS (84%/46%), TP53 (75%/48%), SMAD4 (29%/17%), CDKN2A (29%/12%), and TTN (14%/10%) ( Figure 8F ). Since this gene signature was closely related to gemcitabine sensitivity, we further explored its potential in predicting the response to other commonly used anti-PC drugs. As shown in Figure 8G , patients with higher risk scores were more likely to be resistant to doxorubicin, olaparib, paclitaxel, and FOLFRINOX components (fluorouracil, irinotecan, oxaliplatin), while the same tendencies were also observed in patients with higher TMB levels.

Moreover, we employed the GSVA method to understand the mechanisms underlying chemoresistance. The pathways enriched in the high-risk group were mainly associated with pancreatic cancer, apoptosis, mismatch repair, cell cycle, spliceosome, hypoxia, glycolysis, P53 signaling, MYC signaling, TGF-beta signaling, and PI3K/Akt signaling ( Figures 9A, B ). The low-risk group showed enrichment for pathways related to metabolism and downregulation of KRAS. Thus, hyperactivation of oncogenic pathways could account for the dismal prognosis and chemoresistance in the high-risk group. Taken together, this gene signature was a versatile tool favoring multidrug resistance prediction.

Seven published methods were integrated to comprehensively evaluate the TME compositions of PC patients in TCGA and GEO cohorts. The risk score exhibits positive correlations with cells favoring immunosuppression, including myeloid-derived suppressor cells (MDSCs), regulatory T (Treg) cells, cancer-associated fibroblasts (CAFs), and T helper type 2 (Th2) cells ( Figure 10A ) (39). Consistently, the anticancer cells were significantly lower in the high-risk group, such as CD8+ T cells, cytotoxic cells, and effector memory T (Tem) cells ( Figure 10B ). The TIDE method was applied to estimate the ICB resistance score, and the results showed that the low-risk group was more sensitive to ICB treatment ( Figure 10C ).

Notably, the high-risk group patients had higher TMB levels and higher expression levels of immune checkpoints (PDCD1 and PD-L1), which was contradictory to the TIDE estimation ( Figures 8E , 10C, D ). Then, the “easier” R package was employed to objectively assess the ICB response by leveraging multiple proxies of the immune response, including cancer type, TMB, pathway activities, cell fractions, transfer factor activities, and intracellular and intercellular communications. Interestingly, TMB exhibited a negative correlation with the EaSIeR score ( Figure 10E ). However, there was no significant difference regarding the EaSIeR scores between the high- and low-risk groups ( Figure 10F ).

Because previous studies reported that EGFR, ERAP2, and MSLN were associated with gemcitabine response, siRNA-mediated knockdowns were carried out to investigate the roles of NCEH1 and ALDH3B1 in cell proliferation and gemcitabine resistance (40–42). First, RT-qPCR analyses demonstrated that the expression levels of ALDH3B1 and NCEH1 were relatively higher in the three PC cell lines (AsPC-1, CFPAC-1, and PANC-1) than in the normal pancreatic epithelial cell line (hTERT-HPNE) ( Figures 11A, B ). At the mRNA protein level, high silencing efficiencies (> 70% for PCR assay, > 40% for WB assay) of both siRNAs were achieved at 48 hours posttransfection ( Figures 11C-H ; Supplementary Figure S3A-C ).

The CCK-8 proliferation assay indicated that the proliferation rates of the three PC cell lines were noticeably diminished after knocking down ALDH3B1 and NCEH1 ( Figures 11I-K , Supplementary Figure S3D-F ). Gemcitabine cytotoxicity assays demonstrated that knockdown of ALDH3B1 and NCEH1 enhanced gemcitabine sensitivity ( Figures 11L-N ). Of note, this effect was concentration-dependent. The concentration threshold for viability differences was 0.01/0.1/10 for CFPAC-1/AsPC-1/PANC1 cells. These findings indicated that NCEH1 and ALDH3B1 were pivotal regulators of cell proliferation and gemcitabine sensitivity in PC.