Deciphering the heterogeneity of epithelial cells in pancreatic ductal adenocarcinoma: implications for metastasis and immune evasion

Abstract

Objective

This study examines the cellular heterogeneity of epithelial cells within pancreatic ductal adenocarcinoma (PDAC) and their contributions to tumor progression, metastasis, and immunosuppressive interactions using single-cell RNA sequencing.

Methods

Single-cell RNA-sequencing data from two datasets (GSE154778 and GSE158356) were integrated using the Harmony algorithm, followed by quality control, clustering, and differential gene expression analysis. Distinct subpopulations of epithelial cells were identified, and their gene expression profiles were analyzed. To assess the malignancy of these subpopulations, single-cell copy number variation (CNV) analysis and trajectory analysis were conducted. Additionally, intercellular communication was examined using the CellChat platform.

Results

The analysis revealed pronounced heterogeneity among PDAC epithelial cells, with specific subpopulations exhibiting distinct roles in tumor proliferation, extracellular matrix remodeling, and metastatic dissemination. Subpopulations 4 and 6 were characterized by increased CNV levels and a more malignant phenotype, suggesting an enhanced capacity for metastasis. Single-cell trajectory analysis, along with CellChat, mapped the temporal evolution of epithelial cells, identifying key regulatory genes such as DCBLD2 and JUN. A prognostic model incorporating five key genes, including KLF6, was developed and demonstrated strong predictive accuracy for patient outcomes. Notably, KLF6 emerged as a critical prognostic marker associated with immune modulation, particularly through interactions with M2 macrophages.

Conclusion

The study highlights the pronounced heterogeneity of epithelial cells in PDAC and their distinct contributions to tumor progression, metastasis, and immune modulation. Through single-cell transcriptomic and CNV analyses, we identified epithelial subpopulations with varying malignant potentials and distinct interactions with the tumor microenvironment. Among these, KLF6 emerged as a key regulator associated with immune modulation and metastasis. Our findings emphasize the significance of epithelial cell heterogeneity in shaping pancreatic cancer progression. These insights provide a foundation for future investigations into novel prognostic markers and therapeutic strategies.

Introduction

Pancreatic ductal adenocarcinoma (PDAC) is among the most lethal forms of cancer, with a five-year survival rate below 12%, placing it in the lowest survival category across all cancers [1]. The high mortality rate of PDAC largely arises from its frequent diagnosis at advanced stages, challenges in surgical resection, rapid metastatic spread, and significant resistance to existing treatments such as chemotherapy [2, 3]. Currently, only a small proportion of patients are eligible for potentially curative surgery, with most presenting with unresectable locally advanced disease or distant metastases [4]. Despite advances in therapeutic approaches, including adjuvant and neoadjuvant chemotherapy, the introduction of immunomodulatory drugs, such as checkpoint inhibitors, has not significantly improved PDAC treatment outcomes [5, 6]. This underscores the urgent need for a deeper understanding of PDAC biology to facilitate the development of more effective therapies.

Epithelial cells play integral roles in pancreatic cancer progression, contributing actively to tumor development through a complex array of genetic and epigenetic changes that drive their transformation from normal to malignant phenotypes [7, 8]. These epithelial cells interact with immune cells, fibroblasts, vascular endothelial cells, and the extracellular matrix within the tumor microenvironment (TME) to collectively promote tumor growth, invasion, and resistance to therapy [9]. The heterogeneity of distinct epithelial cell subpopulations further complicates treatment, as these subpopulations vary in their therapeutic responses and may adopt invasive properties via mechanisms such as epithelial–mesenchymal transition (EMT) and metabolic pathway reprogramming, which support rapid growth and enhance immune evasion [10,11,12]. Expanding our understanding of these characteristics is essential for developing new and more effective therapeutic strategies, including targeted and immunotherapeutic approaches.

In this context, single-cell RNA sequencing (scRNA-seq) has emerged as an invaluable tool in precision oncology, enabling comprehensive characterization of the pancreatic cancer TME and identification of new therapeutic targets [13]. Recent studies have shed light on the heterogeneity of epithelial cells in PDAC using scRNA-seq, but the specific roles of these subpopulations in tumor progression, metastasis, and immune evasion remain unclear. While Werba et al. and Zhang et al. have explored the TME and identified key immune-related pathways, their focus has primarily been on immune cell interactions rather than the functional roles of epithelial cells [14, 15]. In contrast, our study uniquely integrates scRNA-seq data from primary and metastatic tumors, providing a comprehensive analysis of epithelial cell heterogeneity and its implications for metastasis and immune modulation. Given the critical role of epithelial cells in tumor progression, metastasis, and resistance to treatment—a relatively underexplored area—this study integrates single-cell data from two datasets, complemented with bulk RNA-seq data, to investigate the role of epithelial cells within the TME. The findings underscore the central involvement of epithelial cells in the development of pancreatic cancer and suggest potential molecular markers and therapeutic strategies for future clinical applications.

Methods

Data collection

The TCGA-PAAD dataset, including gene expression data (FPKM) and corresponding survival information, was obtained from the UCSC Xena platform (https://xenabrowser.net/datapages/). No specific inclusion or exclusion criteria, such as tumor purity thresholds, were applied to the samples. To ensure robustness in model validation, the dataset was randomly divided into training and validation cohorts at a 1:1 ratio. The scRNA-seq data were retrieved from the GEO database under accession numbers GSE154778 and GSE158356 [16, 17]. In total, 12,707 cells from GSE154778 and 4,890 cells from GSE158356 were included in the analysis, with both datasets sequenced on the Illumina HiSeq 4000 platform.

Data processing

Preprocessing, quality control, normalization, and clustering of single-cell data were conducted using Seurat v4.3.0. DoubletFinder v2.0.3 was applied to detect and process potential doublets within the single-cell data. Quality control criteria required each cell to express at least 400 genes and limited mitochondrial gene expression to a maximum of 20%. Subsequent procedures, including data normalization, selection of highly variable genes, and dimensionality reduction clustering, followed Seurat's default parameters and standard workflow. Harmony v1.2.3 was used to integrate data from multiple samples [15]. To assess integration quality, we visualized the uniform manifold approximation and projection (UMAP) distributions before and after batch correction, as shown in Fig. 1A. The correction successfully mitigated batch effects while preserving meaningful biological variations, indicating that Harmony effectively integrated the data. Batch effect correction was validated using the Local Inverse Simpson’s Index (LISI). LISI quantifies both batch mixing (higher scores indicate better integration) and cell type specificity (lower scores indicate preserved biological clusters). Cell clusters were annotated by referencing known marker genes from the literature and through manual annotation. The FindAllMarkers function identified differentially expressed genes among cell subgroups, using the Wilcoxon test with an adjusted p-value threshold of < 0.05 and an absolute log2 fold change (|log2FC|) of > 0.5.

Fig. 1

Overview of single-cell RNA-seq data analysis from GSE154778 and GSE158356. A UMAP visualization illustrates the distribution of single-cell RNA-seq data before (left) and after (right) batch effect correction. B t-SNE plot showing the clustering of cells into 30 distinct clusters. C Bubble plots depict the expression levels of marker genes across each cell type. D t-SNE plot revealing the distribution of nine identified cell types, along with a stacked bar graph depicting the sample distribution proportions for each cell type. E Heatmap illustrating cell type annotation specificity, with labels shown along the bottom axis. F t-SNE visualization of cell type distributions in metastatic versus in situ cancers. G Stacked bar graph displaying the percentage distribution of each cell type in metastatic and in situ cancer samples

Cell type identification

Marker genes for each cell cluster were identified through differential expression analysis using Seurat’s FindAllMarkers function. Criteria for marker gene identification included an adjusted p-value of < 0.05, an expression percentage > 0.25, and an absolute log2 fold change (|log2FC|) > 0.25. Cell clusters were subsequently identified and annotated using the singleR package, based on the compositional patterns of marker genes. These annotations were manually reviewed and refined with reference to the CellMarker database.

Cellular communication

The "CellChat" R package v1.1.3 was used to analyze cellular communication within the tumor microenvironment, with a focus on receptor-ligand interactions [18]. CellChat quantifies interaction links and compiles communication probabilities, constructing interaction networks that represent the duration and cumulative strength of interactions between any two cell populations. Scatter plots were generated to visualize primary sender (signal source) and receiver cells (targets) within a two-dimensional space, facilitating the identification of significant contributors to outgoing or incoming signals among immune cell populations. A global communication model based on pattern recognition was used to examine the collaborative functions of multiple immune cell types and signaling pathways.

Pseudotemporal trajectory analysis

Cellular pseudotemporal trajectories were constructed using the Monocle 2 algorithm, an R package developed by Qiu et al. for single-cell trajectory analysis [19]. This algorithm applies machine learning techniques to condense high-dimensional expression data into a lower-dimensional space, arranging cells into trajectories that include branch points. Dynamic expression heatmaps were generated using the plot_pseudotime_heatmap function.

Immune infiltration assessment

The CIBERSORT algorithm was used to quantitatively assess immune cell infiltration levels in patients with pancreatic adenocarcinoma (PAAD), focusing on differences in cell abundance between high-risk and low-risk groups [20]. Additionally, Pearson’s correlation was calculated to examine the relationship between immune cell abundance and patient risk scores. To identify potential differences in immune function, enrichment scores were derived from single-sample gene set enrichment analysis (ssGSEA) [21]. The Wilcoxon test was subsequently applied to compare immune function between high-risk and low-risk patient groups.

Statistical analyses

All statistical analyses and data visualizations were conducted using R software (version 4.1.3). Pearson's correlation coefficient was applied to assess correlations between continuous variables. For quantitative data, subgroup comparisons were conducted using two-tailed, unpaired Student's t-tests or one-way Analysis of Variance (ANOVA) with Tukey’s multiple comparison test, as appropriate. To control for false positives arising from multiple comparisons, p-values were adjusted using the Benjamini-Hochberg (BH) method to maintain a controlled false discovery rate (FDR). A p-value of < 0.05 (or an adjusted p-value < 0.05 for multiple comparisons) was considered indicative of statistical significance.

For single-cell and bulk RNA-seq analyses, standard cutoff values were applied to ensure data quality and biological interpretability: a mitochondrial gene threshold of 20% to filter out low-quality cells, a |log2FC|> 0.5 to identify differentially expressed genes, and an expression percentage > 0.25 to focus on genes with meaningful expression levels. These thresholds are widely supported in literatures and were chosen to balance sensitivity and specificity in our analyses [15, 16].

Results

Single-cell sequencing data reveals cellular heterogeneity in primary and metastatic pancreatic cancer

In this study two single-cell sequencing datasets from the GEO database, GSE154778 and GSE158356, were integrated and analyzed, with quality control measures implemented to ensure data integrity, including the screening of low-quality cells and the identification and removal of doublet contamination. The Harmony algorithm was applied to correct for batch effects between samples (Fig. 1A), followed by data normalization, centering, and principal component analysis (PCA) dimensionality reduction. The first 30 principal components were selected for UMAP and t-distributed stochastic neighbor embedding (T-SNE) calculations. Cells were categorized into 30 clusters (Fig. 1B), and the sources of single-cell samples are revealed in Figure S1A, B. LISI scores for batch mixing ranged from 1.000 to 2.000 (mean = 1.333), indicating effective but non-overcorrection of technical variation.

Based on cell type-specific genes reported in the literature and singleR analysis, the 30 clusters were further categorized into 9 distinct cell types (Fig. 1C, Figure S1) [14]. The heterogeneity of cell types within each sample is depicted in Fig. 1D, highlighting epithelial cells, fibroblasts, and T cells as the three major cell subgroups. The top 50 most highly expressed genes were calculated for each cell type, allowing individual cell types to be clearly distinguished and demonstrating the reliability and specificity of the cell annotations (Fig. 1E). Figure 1F, G and Figure S1C–D illustrate the T-SNE distributions of each cell type in both primary and metastatic cancers, as well as the corresponding ratio maps. It was observed that the proportions of NK/T cells and myeloid cells were significantly lower in primary cancers compared to metastatic cancers, consistent with the immunosuppressive environment and reduced immune cell infiltration commonly reported in primary pancreatic cancer. [22,23,24]

Association of epithelial cell subtypes with biological processes related to pancreatic cancer progression

Epithelial cells serve as the primary source of malignant cells in pancreatic cancer and comprise a large proportion of the single-cell data analyzed. Hierarchical clustering was performed on epithelial cells to analyze intratumoral epithelial cell heterogeneity, resulting in the identification of 10 distinct epithelial cell subtypes, as illustrated in Fig. 2A and Figure S2A. Notable, subgroups 4, 5, 6, and 9 exhibited significantly higher representation in metastatic cancer samples compared to in situ cancer samples. Subgroup 1 was observed only in metastatic samples, while subgroups 8 and 10 were found exclusively in in situ cancer samples (Fig. 2B). The distribution of epithelial cell subtypes across individual samples is depicted in Figure S2B–C.

Fig. 2

Characterization of epithelial cell populations in pancreatic cancer. A t-SNE plot illustrating the segmentation of epithelial cells, with distinct colors representing various cell types; the left plot corresponds to metastatic cancer, while the right plot represents carcinoma in situ. B Stacked bar graph displaying the proportion of cells within each group. C Left panel: Graphs depicting the dynamic expression patterns of representative differentially expressed genes (DEGs) across epithelial cell populations. Middle panel: Heatmap highlighting representative DEGs between cell clusters. Right panel: Representative enriched gene ontology (GO) terms for each cell cluster

To analyze gene expression differences among these subtypes and their unique gene signatures, relative expression levels of each gene were calculated in single cells, with scores assigned accordingly. An unsupervised clustering approach then categorized these genes, yielding distinct gene clusters and identifying nine expression trends. These trends highlighted potential roles for these genes in various biological processes (Fig. 2C). For instance, the IGF2 gene, known to be involved in cell proliferation and malignant progression, exhibited increased expression in subgroup 1 (Figure S2D) [25]. Genes associated with developmental processes and inflammatory responses, such as SOX4 and IL8, were highly expressed in subgroup 2 (Figure S2). This subgroup was linked to stress and immune responses, including the regulation of RNA polymerase II promoter transcription in response to stress, modulation of DNA template transcription during stress, and reactions to tumor necrosis factor [26, 27]. Additionally, subgroups 4, 5, 6, and 9 in metastatic cancer samples were involved in processes such as cell proliferation, division, cellular structure maintenance, and immune response-related pathways.

These findings indicate that differences among epithelial cell subtypes may contribute to the growth and metastasis capabilities of the tumor, including processes related to cell proliferation and extracellular matrix remodeling. Specific gene expression patterns may therefore provide a biological foundation for the development and progression of pancreatic cancer.

Single-cell CNV analysis to assess the malignancy of epithelial cell subpopulations in pancreatic cancer

To further assess the malignant status of epithelial cell subpopulations, single-cell CNV profiles were inferred using inferCNV, with T cells, B cells, and NK cells serving as references for comparison [28]. As anticipated, CNVs exhibited considerable heterogeneity across epithelial cell subgroups, with varying degrees of CNV accumulation observed. For instance, significant copy number amplifications were identified on chromosome 8, while notable deletions were detected on chromosomes 3 and 6. Amplification on chromosome 1 was particularly pronounced in subgroup 1, whereas amplification on chromosome 12 was observed across all subgroups except subgroup 1 (Fig. 3A).

Fig. 3

Analysis of (CNV) and gene expression in epithelial cell subtypes. A inferCNV analysis depicting inferred copy number variation for each epithelial cell subtype, with T cells, B cells, and NK cells serving as reference populations. B CNV scores mapped onto the epithelial cell t-SNE plot, where the color intensity represents the CNV score level. C Box plots depicting CNV scores across 10 distinct epithelial cell subtypes. D, E Expression profiles of TOP2A and ITGB1 within epithelial cells. F Differential gene analysis across cell types, conducted using the Wilcoxon rank-sum test; “Highly” indicates adjusted p-values below 0.05, while “Lowly” denotes adjusted p-values below 0.1. G GSEA heatmap displaying the 50 marker gene sets from the MSigDB database, arranged according to subgroups of epithelial cells from in situ and metastatic cancers

CNV scores were calculated for each cell and visualized using a T-SNE plot (Fig. 3B). Quantification of the CNV scores within each subgroup revealed that epithelial cells in subgroups 4 and 6 displayed significantly higher CNV levels compared to other subgroups, indicating a more malignant phenotype. Subgroups 1, 2, and 7 displayed lower CNV levels relative to these subgroups (Fig. 3C).

Genes with elevated expression in each epithelial cell subgroup were identified, with DNA topoisomerase II alpha (TOP2A) significantly overexpressed in subgroup 4 (Fig. 3D). Aberrant expression of TOP2A has been associated with errors in DNA replication and issues in chromosome segregation, which contribute to genomic instability, including CNVs [29]. Subgroup 6 revealed elevated expression of ITGB1, a factor linked to cell adhesion and migration [30].

Differentially expressed genes (DEGs) were characterized in each epithelial cell subgroup when comparing in situ and metastatic cancers, with the top five genes exhibiting significantly high or low expression in each subgroup depicted in a scatter plot using the FindMarkers function (Fig. 3). In subgroup 4, the genes PLAT, LY6D, and ASPH were significantly upregulated in metastatic cancers, indicating heightened invasiveness and metastatic potential of subgroup 4 epithelial cells in metastatic cancers. This finding aligns with the malignancy assessment of subgroup 4 based on CNV analysis.

To further characterize the biological pathway activities across the epithelial subtypes, hallmark gene sets from the Molecular Signatures Database (MSigDB) were utilized in conjunction with the Gene Set Variation Analysis (GSVA) algorithm (Fig. 3G). mTORC1_Signaling, E2F_Targets, and G2M_Checkpoint pathways were notably upregulated, indicating abnormal cell growth and metabolic activity as well as dysregulated cell cycle control, with these effects being more pronounced in metastatic cancers. Additionally, metastatic cancer epithelial cells in subgroup 2 exhibited increased activity in the P53 pathway and heme metabolism. The CNV analysis revealed that subgroups 4 and 6 exhibited the highest genomic instability, suggesting a strong link to malignancy. To further explore the functional consequences of these genetic alterations, we conducted single-cell trajectory analysis to map their dynamic transitions during cancer progression.

Single-cell trajectory analysis of epithelial cell subsets in pancreatic cancer

To further analyze dynamic shifts within epithelial cell subsets, single-cell trajectory analysis was performed using the Monocle package, identifying three distinct cell states as shown in Fig. 4A and B. States 1 and 2 were associated with metastatic carcinoma in pancreatic cancer, while state 3 represented the carcinoma in situ state (Fig. 4C). Epithelial cells from subgroups 2 and 3 predominantly mapped to the in situ cancer endpoints within the trajectory, as depicted in Fig. 4D. We calculated the proportion of each epithelial cell subtype in each state to validate our observations (Fig. 4E).

Fig. 4

Monocle trajectory analysis of epithelial cell subgroups and gene expression dynamics. A–D Monocle trajectory analysis of epithelial cell subgroups, with cells colored by (A) inferred pseudotime, (B) developmental status, (C) group designation, and (D) sample origin. (E) The proportion of each epithelial cell subtype in each state. F Pseudotemporal heatmap displaying expression changes in the top 30 genes identified by BEAM analysis before and after node 1. G Pseudotemporal heatmap revealing the top 30 genes with significant expression shifts along the trajectory. H Gene expression patterns of DCBLD2, GLUL, JUN, and TAGLN2, with color coding based on expression levels throughout the pseudotime progression. I Expression levels of DCBLD2, GLUL, JUN, and TAGLN2 before and after node 1, displayed along pseudotime and color coded by group.(J) Violin plots illustrating the expression of GLUL, JUN, and TAGLN2 in in situ versus metastatic cancers

Dynamic gene expression changes along these trajectories were examined (Fig. 4F) and differential gene trajectories before and after node 1 were analyzed, highlighting the top 30 genes by significance (Fig. 4G) and extending to the top 200 genes (Figure S3A, B). TAGLN2, DCBLD2, JUN, and GLUL emerged as key genes in this analysis. Trajectories for these genes across the projected timeline are shown in Fig. 4H–J.

DCBLD2, a transmembrane receptor implicated in the regulation of tumor neovascularization and metastasis via EMT, displayed low expression in in situ cancers but increased significantly toward the terminal phase of metastatic carcinoma, indicating a time-dependent role for DCBLD2 in metastatic progression. GLUL is involved in glutamate and amino acid metabolism, JUN contributes to cell proliferation, differentiation, and apoptosis as part of the AP-1 transcription factor complex; and TAGLN2, associated with the cytoskeleton, is thought to regulate tumor invasiveness and metastatic potential [31,32,33].

The expression of these genes increased during the early stages of the temporal trajectory, declined gradually, and partially recovered toward the endpoint associated with carcinoma in situ, as illustrated in Fig. 4H and depicted in the violin plot (Fig. 4I). These findings suggest a progressive transformation in epithelial cell states during pancreatic cancer metastasis, with invasion- and metastasis-related genes like TAGLN2, DCBLD2, JUN, and GLUL following specific expression patterns. This continuity is essential for the metastatic regulation of tumor cells and indicates a regulatory role that transcends individual gene expression levels. Further experimental validation and functional studies are necessary to elucidate the role of these genes in pancreatic cancer progression.

Analysis of intercellular communication networks in in situ and metastatic cancers using CellChat

Given that epithelial subpopulations exhibit distinct malignant properties and trajectory patterns, understanding their interactions with immune and stromal cells is crucial for elucidating their role in PDAC progression. To clarify intercellular communication among epithelial cell subpopulations and other cell types in both in situ and metastatic cancers, CellChat analysis was employed (Figure S3), leveraging known ligand–receptor pairs and their cofactors [18]. This analysis focused on four key subgroups identified through Monocle’s temporal sequence analysis. Subgroups 2 and 3 were linked to the in-situ cancer trajectory endpoint, while subgroups 4 and 6 were linked to the metastatic cancer trajectory. Subgroups with low cell numbers were excluded from this analysis. Notably, subgroups 4 and 6 exhibited slightly higher afferent and efferent signaling capacities within the immune communication network (Fig. 5A). The distinct expression patterns of epithelial cells, analyzed via Cophenetic- and Silhouette-based clustering, are shown in Figure S4.

Fig. 5

Interaction analysis between epithelial cell subtypes and immune cells within the tumor microenvironment. A Circle plot illustrating the number (left) and strength (right) of interactions between specified epithelial cell subtypes and immune cells. B Bubble plot revealing interactions between overexpressed ligands and receptors, where bubble size represents the p-value from the alignment test, and color indicates interaction likelihood, with epithelial cell subtypes acting as signal senders. C Visualization of receptor-ligand interactions generated by epithelial subgroups 4 and 6 as signal emitters. D Receptor-ligand interactions received by epithelial subgroups 4 and 6 as signal recipients. E Contribution of each cell type to reciprocal interactions within the PLAU pathway. F Algorithm-based inference of the role of each cell type in the PLAU pathway. G Expression levels of PLAU-related molecules across cell types

The role of epithelial subsets as both signal recipients and initiators was examined, focusing on their ligand–receptor interactions with immune cell subsets (Fig. 5B, Figure S3D). Notably, LGALS9 signaling, involving interactions with P4HB, CD45, and CD44, was specific to subgroups 2 and 3. As LGALS9 is known to regulate immune responses by inducing T-cell apoptosis and modulating T-cell activity, this indicates that subgroups 2 and 3 may facilitate immune evasion by modulating T-cell activity within the pancreatic cancer immune environment [34, 35]. In contrast, subgroups 4 and 6, characterized by high malignancy, predominantly signaled through immune escape-related MIF, metastasis-associated MK, and invasion-linked LAMININ pathways. The primary target cells for these pathways were NK and myeloid cells, further aligning with the aggressive profile of subgroups 4 and 6. Additionally, fibroblast-derived COLLAGEN signaling was notably prominent (Fig. 5C, D).

In subgroups 4 and 6, the interaction between PLAU and its receptor PLAUR was identified, a pathway notably absent in subgroups 2 and 3 (Fig. 5B). PLAU, a serine protease from the S1 family, facilitates fibrinolytic enzyme activation from fibrinogen, directly degrading extracellular matrix (ECM) components (e.g., collagen, laminin), contributing to the progression of various cancers, including squamous cell carcinoma, colorectal cancer, and pancreatic cancer [36, 37]. PLAU and its receptor, PLAUR, play key roles in tumor invasion, cell migration, and immune regulation within the tumor microenvironment, with established associations with lymph node metastasis and the formation of tumor-associated fibroblasts [38, 39]. This indicates that the metastatic potential of subgroups 4 and 6 may be linked to PLAU-mediated molecular pathways and interactions with myeloid cells (Fig. 5E).

Our CellChat analysis further demonstrated that subgroups 4 and 6 emerged as primary Senders, Mediators, and Influencers of the PLAU signaling pathway, with these cells exhibiting significantly elevated PLAU expression levels (Fig. 5F, G). These findings suggest that PLAU-PLAUR signaling is predominantly active in metastatic samples, underscoring its critical role in promoting tumor invasiveness and metastatic dissemination. When assessing signaling patterns, all four epithelial cell subgroups exhibited similar afferent and efferent signaling trends (Figure S5), further supporting the broad relevance of this pathway across different cellular contexts.

Prognostic marker screening using monoclonal trajectory analysis and TCGA-PAAD to differentiate patient risk

To further examine the biological functions of epithelial cell subgroups, Monocle analysis was combined with prognostic data from the TCGA-PAAD dataset. This integrated approach, aligning gene expression profiles with survival information, allowed univariate regression analysis to identify 3,333 genes significantly linked to patient prognosis (p < 0.05). These were cross-referenced with 8,426 significant genes (p < 0.05) through BEAM analysis, resulting in a final set of 1,329 key genes (Fig. 6A). These genes were subsequently refined using a random forest algorithm, identifying five genes—MEIG1, EDN2, DGCR5, CXCL11, and KLF6—as core prognostic markers. These genes were incorporated into a comprehensive prognostic scoring model tailored for pancreatic cancer (Figure S6A–C).

Fig. 6

Prognostic model evaluation and immune infiltration analysis in pancreatic cancer. A Venn diagram depicting the overlap between significant genes identified by one-way Cox regression and trajectory genes with significant changes identified by Monocle. B Kaplan–Meier (K-M) survival curves for the 5-gene prognostic model in the training set. C Kaplan–Meier survival curves for the 5-gene prognostic model in the validation set. D Kaplan–Meier survival curves for the 5-gene prognostic model across the entire dataset. E, F ROC curves assessing model performance with AUCs at 1, 5, 7, and 10 years for both the training set and the entire dataset. G Immunological infiltration in high- and low-risk groups, as determined by the model, assessed using the CIBERSORT approach. H Violin plots revealing immune cell infiltration across different categories. I Correlation analysis between KLF6 expression and immune infiltration levels of M2 macrophages, neutrophils, and Tregs

To assess the predictive performance of the model, the TCGA-PAAD group was divided into two equal groups, with survival curves generated for each. The results revealed significant prognostic differences between these groups (p < 0.0001 for the training set, p = 0.0077 for the validation set, and p = 0.00047 for the overall sample; Figs. 6B–D). For the training set, the model achieved one-, three-, and five-year AUC values of 0.81, 0.71, and 0.77, respectively (Fig. 6E), while in the overall sample, the AUCs were 0.74, 0.76, and 0.76 (Fig. 6F), indicating robust predictive validity.

Immune cell infiltration levels in the TCGA-PAAD high- and low-risk groups were further analyzed using the CIBERSORT algorithm (Fig. 6G, Figure S6D). Notably, patients in the high-risk group exhibited significantly increased infiltration of immune cells, particularly M1-type macrophages (p < 0.0001) and resting dendritic cells (p = 0.004) (Fig. 6H). The relationship between the expression of each of the five prognostic genes and immune cell infiltration in single-cell data was also assessed (Figure S6E). Among these, KLF6, a gene with limited prior exploration in pancreatic cancer research emerged as noteworthy, due to its strong association with increased infiltration of M2-type macrophages, neutrophils, and regulatory T cells in the TCGA-PAAD dataset (Fig. 6I), indicating a potential role in modulating the tumor microenvironment and facilitating immune evasion. To identify specific molecular drivers of these aggressive phenotypes, we next examined the role of KLF6, which emerged as a key regulator in our trajectory and prognostic analyses.

KLF6 expression in pancreatic cancer and its association with poor patient prognosis

In a detailed single-cell analysis of pancreatic cancer, KLF6 expression was found to be prevalent in epithelial cells, myeloid cells, and fibroblasts (Fig. 7A). Within the epithelial cell population, KLF6 expression patterns varied based on cancer type: in situ cancer cells displayed nearly uniform KLF6 expression, while metastatic cancers showed a more balanced distribution of KLF6-positive and KLF6-negative epithelial cells (Fig. 7B, Figure S7A, B). Further gene expression comparisons between KLF6-positive (KLF6( +)) and KLF6-negative (KLF6( −)) epithelial cells identified distinct biological pathways associated with each subset. Notably, gene set enrichment analysis (GSEA) indicated significant activation of immune regulation, cell adhesion, migration, and proliferation pathways in KLF6( +) cells (Fig. 7C).

Fig. 7

KLF6 expression and its implications in pancreatic cancer epithelial cells. A KLF6 expression levels across different cell types. B UMAP plot of epithelial cells, colored by KLF6 expression, where KLF6( +) indicates KLF6-expressing epithelial cells and KLF6( −) indicates non-expressing cells. C GSEA results revealing relevant entries based on differentially expressed genes between KLF6( +) and KLF6( −) cells. D Receptor-ligand interactions generated by KLF6( +) and KLF6( −) epithelial cells as signal senders. E, F Bubble plots indicating significant signaling roles of individual cell types in both outgoing (efferent) and incoming (afferent) signaling modes. G–I Expression levels of EDN1, OCLN, and EPHA2 in KLF6( +) versus KLF6( −) epithelial cells. J Correlation between elevated KLF6 levels and poorer patient prognosis in the TCGA-PAAD dataset

Analysis using the CellChat tool revealed that both KLF6( +) and KLF6( −) epithelial cells received COLLAGEN signals from fibroblasts (Fig. 7D). However, only KLF6( +) cells displayed outgoing signaling activity through the OCLN and EDN pathways. OCLN, a protein essential for cellular tight junctions, is involved in barrier function and cell polarity, while EDN is involved vascular tension, cellular proliferation, and neointima formation, potentially enhancing tumor aggressiveness [40,41,42]. KLF6( +) epithelial cells emerged key players in the OCLN, EDN, and EPHA pathways (Figure S7D–F), suggesting that these cells may enhance invasive and metastatic potential in pancreatic cancer by activating OCLN and EDN signaling, fostering a malignant cellular phenotype.

In terms of signal reception, KLF6( +) epithelial cells uniquely received EPHA signaling, a pathway that influences cell adhesion, morphology, and migration [43]. Additionally, KLF6( +) epithelial cells specifically delivered HLA-F signaling, distinct from KLF6( −) cells (Figure S7), while both KLF6( +) and KLF6( −) subsets displayed similar afferent and transmitting modes (Figure S8).

Comparative analysis of EDN1, OCLN, and EPHA2 expression patterns in epithelial cells revealed that these genes were expressed at significantly higher levels in KLF6-positive (KLF6( +)) cells than in KLF6-negative (KLF6( −)) cells (Fig. 7G–I). These observations support the hypothesis that KLF6( +) epithelial cells may contribute to increased tumor invasiveness and metastatic potential by activating pathways such as EDN, OCLN, and EPHA. These pathways are likely linked to abnormal angiogenesis and intercellular communication within the tumor microenvironment. This hypothesis is further corroborated by analysis of data from the TCGA-PAAD dataset (Fig. 7).

KLF6 influences tumor cell responses to a range of drugs

Previous analyses highlighted the role of KLF6 in promoting malignant transformation in pancreatic cancer epithelial cells, as well as its prognostic significance. To further examine its potential as a therapeutic target, a comprehensive drug sensitivity analysis was performed. Drug sensitivity data were collected from the NCI-60 cell line panel along with RNA sequencing data from the National Cancer Institute (NCI) database. Our analysis included 75 clinical trials and 188 FDA-approved drugs to assess the relationship between KLF6 expression levels and drug sensitivity, specifically examining IC50 values. Pearson's correlation coefficients were calculated to assess the relationships between KLF6 expression and each drug (Fig. 8A). The results indicated that high KLF6 expression was associated with enhanced tumor cell sensitivity to several chemotherapeutic agents, including the periwinkle alkaloid antitumor drugs vinblastine and vinorelbine, the bifunctional alkylating agent antibiotic mitomycin, the antimetabolite fluorouracil, the estrogen receptor modulators tamoxifen and raloxifene, the protein kinase inhibitor staurosporine, and the topoisomerase II inhibitor dexrazoxane.

Fig. 8

Investigation of KLF6 expression and its association with drug sensitivity and gene coexpression networks in pancreatic cancer. A Analysis of the relationship between KLF6 expression levels in tumor cells and drug sensitivity for eight compounds. B High-dimensional Weighted Gene Coexpression Network Analysis (hdWGCNA) based on KLF6( +) cells. C UMAP plot of gene coexpression, displaying 11 modules identified by hdWGCNA. D hdWGCNA module network graph illustrating the top 25 genes from modules 5, 7, and 9, based on gene linkage. Each edge represents a co-expression relationship between two genes within the network. The 10 most central genes are positioned at the center, with the remaining 15 genes arranged in an outer circle

To further examine the co-expression networks specific to KLF6( +) epithelial cells, high-dimensional weighted gene co-expression network analysis (hdWGCNA) was applied, a robust tool for examining single-cell RNA sequencing data through co-expression network modeling [44]. This approach identified 11 distinct gene co-expression modules within KLF6( +) cells (Fig. 8), each representing potential functional groupings of genes.

Module 9 included genes heavily associated with cell cycle regulation and DNA replication, such as BIRC5, RRM2, CDK1, and TOP2A. The dysregulated expression of these genes could influence the cellular response to chemotherapeutic agents like vinblastine, vinorelbine, mitomycin, and fluorouracil, which target proliferative pathways. In Module 5, genes central to cell adhesion and signaling, including ITGB1, ITGB4, and MET, were prevalent. These genes are implicated in pathways that may affect the efficacy of hormone-related drugs like tamoxifen and raloxifene, indicating a potential mechanism by which KLF6( +) cells might alter drug responsiveness.

Furthermore, Module 7 comprised of genes involved in stress responses and apoptosis regulation, such as JUN, ATF3, HSPA1A/B, and DUSP1. This indicated that stress and apoptosis pathways may account for the enhanced sensitivity observed in drugs like staurosporine and dexrazoxane, as depicted in Fig. 8A.

Discussion

The present study involved the integration and analysis of 21 samples sourced from two distinct datasets to investigate potential common mechanisms associated with pancreatic cancer. The intricate dynamics of the TME in PDAC were meticulously examined by using a single-cell RNA sequencing (scRNA-seq) data. This methodology facilitated an in-depth exploration of cellular interactions, transcriptional heterogeneity, and the immune contexture within PDAC, thereby providing a comprehensive understanding of how epithelial cell subpopulations contribute to disease progression.

A notable finding was the significant reduction of NK/T cells and myeloid cells in primary PDAC tumors compared to metastatic sites. This observation aligns with the immunosuppressive characteristics of the primary PDAC microenvironment, as documented in prior research [45, 46]. Liudahl et al. (2021) identified a correlation between reduced immune infiltration in primary PDAC and a microenvironment that promotes immune evasion and tumor persistence [47]. The reduced presence of NK and T cells in primary tumors implies that these cancers may establish early immunosuppressive niches, diminishing the efficacy of immune surveillance and allowing tumor proliferation. In contrast, metastatic sites often exhibit a more intricate immune landscape, where immune cells, such as macrophages, are recruited to facilitate the development of pre-metastatic niches (PMNs) [48]. While previous studies have explored similar immune landscapes, our study uniquely highlights the dynamic shifts in immune composition between primary and metastatic sites, emphasizing the adaptive nature of immune modulation in PDAC and identifying potential therapeutic windows for immunomodulatory treatments.

The study also identified specific epithelial cell subpopulations exhibiting high levels of CNV, notably in subgroups 4 and 6, which were more prevalent in metastatic samples. Previous studies have also employed CNV-based malignancy assessments in PDAC and other cancers. Oketch et al. [49] demonstrated that CNV amplifications in pancreatic cancer organoids correlate with tumor aggressiveness and treatment resistance. Similarly, Usman et al. [50] highlighted the stability of CNV patterns in metastatic progression. However, most of these studies relied on bulk sequencing or organoid models, which may obscure intratumoral heterogeneity. Unlike previous studies that focused primarily on CNV alterations themselves, we integrate CNV profiles with single-cell transcriptomics and immune interaction analyses, revealing that high-CNV epithelial subgroups (subgroups 4 and 6) not only exhibit increased metastatic potential but also display unique immune evasion signatures.This evidence supports the notion that targeting CNV-associated pathways may represent an effective strategy to limit metastatic spread and enhance patient outcomes.

The role of EMT in promoting PDAC metastasis has been well established [51, 52]. EMT equips epithelial cells with increased migratory and invasive properties, facilitating their detachment from the primary tumor site and adaptation within distant organs. In this study, high CNV subpopulations also displayed gene expression patterns linked to EMT, including upregulation of markers such as ITGB1, which plays a role in cell adhesion and migration. Our findings extend beyond prior studies by integrating single-cell trajectory analysis to map the temporal evolution of EMT-related genes, offering new insights into the molecular mechanisms underlying metastasis. These findings contribute to the growing body of evidence indicating that therapeutic interventions targeting EMT-associated pathways may provide a valuable approach for disrupting the metastatic process in PDAC.

The trajectory analysis identified some key genes that exhibited dynamic expression patterns during pancreatic cancer progression. For example, DCBLD2 displayed low expression in early-stage carcinoma in situ but significantly increased toward the terminal phase of metastatic carcinoma. This time-dependent upregulation suggests that DCBLD2 may play a pivotal role in promoting EMT and enhancing the invasive potential of epithelial cells during metastasis. Similarly, JUN exhibited a dynamic expression pattern along the pseudotemporal trajectory, with elevated levels in metastatic subpopulations. This observation aligns with prior research indicating that JUN promotes EMT and enhances the migratory capacity of cancer cells by activating downstream effectors such as MMPs (matrix metalloproteinases) and cytokines involved in extracellular matrix remodeling. Future functional studies are needed to validate their roles and explore therapeutic potential.

A key focus of this study was the transcription factor KLF6, a member of the Sp1/KLF family involved in various biological processes, including cell cycle regulation, vascular remodeling, and programmed cell death [53,54,55,56]. Recent evidence indicates that the role of KLF6 may vary across cancer types; it has been investigated as a tumor suppressor in cancers such as oral squamous cell carcinoma [57]. In breast and ovarian cancers, for example, KLF6 and its splice variants have been associated with increased metastatic potential and poorer patient outcomes [58, 59]. Within pancreatic cancer research, Xiong et al. reported that KLF6 inhibits PDAC progression; however, its influence on the malignant transformation of pancreatic epithelial cells has not been thoroughly examined [60].

This study revealed a critical role for KLF6 in promoting malignant transformation in pancreatic cancer epithelial cells. KLF6 high expression in specific epithelial subpopulations correlates with a more aggressive phenotype and poorer prognosis, indicating that KLF6 may exert dual roles in PDAC by advancing tumor progression in specific cellular contexts, particularly through interactions with the immune microenvironment. Our data indicate that KLF6 expression is associated with the recruitment of immunosuppressive cell types, including M2 macrophages, which support tumor growth and metastasis by suppressing anti-tumor immune responses [48]. Notably, the co-existence of KLF6 and PLAU pathways may reflect a cooperative mechanism where KLF6 modulates immune evasion while PLAU promotes extracellular matrix degradation and stromal remodeling, collectively driving metastasis. The main limitation is the lack of direct experimental validation and clinical sample analysis. While our bioinformatics analyses provide strong correlative evidence linking KLF6 to PDAC progression and immune modulation, functional studies are necessary to confirm its mechanistic role. For example, in vitro experiments involving KLF6 knockdown or overexpression in PDAC cell lines could provide insights into its impact on tumor proliferation, invasion, and immune cell interactions. Additionally, in vivo models—such as patient-derived xenografts (PDX) or genetically engineered mouse models (GEMMs)—could help elucidate whether KLF6 directly influences tumor growth, metastatic potential, and immune evasion [61, 62]. In future studies, we will validate these results through in vitro experiments and clinical sample analysis, further evaluate whether targeting KLF6 can be used as a therapeutic strategy to modulate the tumor microenvironment of pancreatic cancer.

Another limitation of this study is that publicly available scRNA-seq datasets limit sample size and may not fully represent the heterogeneity of PDAC in different patient populations. While scRNA-seq provides high-resolution profiling of cellular states and interactions, it is inherently limited by the depth and coverage of the datasets used. As a result, rare cell populations may be underrepresented, potentially restricting our capacity to thoroughly characterize immune and epithelial cell interactions within the PDAC microenvironment.

Future research could benefit from integrating scRNA-seq with spatial transcriptomics to obtain a more comprehensive view of how spatial organization within the tumor influences cellular interactions and gene expression. Such integration would enhance our understanding of how malignant subpopulations, including those expressing KLF6, interact with specific immune cell subsets within the TME.

Given the emerging role of KLF6 in modulating both epithelial cell behavior and immune interactions, examining its efficacy in combination with existing therapies could offer significant clinical value. Prior studies, including those by Bear et al. (2020), have underscored the challenges posed by the dense stromal architecture and immunosuppressive environment of PDAC [5]. Targeting a transcription factor like KLF6 may provide a novel strategy for circumventing these barriers, enhancing therapeutic responses.

In conclusion, this study provides a comprehensive characterization of epithelial cell heterogeneity in PDAC and reveals its implications for metastasis and immune evasion. Through single-cell analysis, we identify distinct epithelial subpopulations with varying malignancy potential and immune interactions. Future studies should focus on experimental validation and mechanistic investigations to further elucidate the functional contributions of these epithelial subpopulations in PDAC progression and immune modulation.