Deciphering the molecular crosstalk between pancreatic ductal adenocarcinoma and type 2 diabetes through multi-dataset integration

Introduction

Now surpassing breast cancer, pancreatic ductal adenocarcinoma (PDAC) ranks as the third most common cause of cancer mortality globally and is projected to overtake colorectal cancer by 2040, ranking second only to lung cancer (1,2). The high lethality of PDAC is primarily due to its typically late diagnosis, often occurring after distant metastasis. Although 10–15% of cases can be linked to germline mutations or established risk factors such as chronic pancreatitis and mucinous cystic neoplasms (3), the majority of patients present without identifiable predisposing factors, further complicating early detection. Early-stage PDAC is typically asymptomatic or presents with non-specific signs, and there are currently no reliable biomarkers for early diagnosis, combined with the pancreas’s deep and inaccessible anatomical location, limiting the effectiveness of routine screening (4). Moreover, the incidence of PDAC continues to rise, driven in part by the global obesity epidemic and increasing life expectancy (5). These challenges emphasize the critical necessity for new diagnostic and therapeutic strategies to enhance outcomes for patients with this highly lethal disease.

Multiple studies have demonstrated that a substantial proportion of PDAC patients develop diabetes months to years prior to their cancer diagnosis (6,7). In such cases, hyperglycemia is associated with the presence of subclinical, undiagnosed pancreatic tumors and is classified as pancreatic cancer-related diabetes (PCRD), or type 3c diabetes (8,9). Although type 2 diabetes mellitus (T2DM) is typically associated with weight gain, pancreatogenic diabetes often presents with weight loss (10). However, the relationship between PDAC risk and the temporal course of diabetes, as well as the degree of weight loss in the general population, remains poorly understood (11). Furthermore, it is unclear whether intentional weight loss affects the accuracy of PDAC risk prediction models (12). Importantly, many patients with pancreatogenic diabetes exhibit localized, resectable tumors at the onset of hyperglycemia, providing a critical window for early cancer detection and intervention.

With the rapid advancement of sequencing technologies and bioinformatics, the application of gene expression profiling has become a powerful tool for identifying key genes and elucidating shared pathogenic mechanisms between diseases (12-15). Herein, we employed weighted gene co-expression network analysis (WGCNA) to determine overlapping genes between the apoptosis- and necrosis-related “lightcyan” module in PDAC and differentially expressed genes (DEGs) in T2DM. Then, we examined the role of these genes in immune modulation and drug sensitivity, aiming to detect occult PDAC in patients with T2DM (16-18). This approach offers valuable insights for the timely detection and clinical management of PDAC. We present this article in accordance with the MDAR and TRIPOD reporting checklists (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2427/rc).

Methods

Data sourcing

We accessed normalized transcriptomic data for pancreatic adenocarcinoma (PAAD) from The Cancer Genome Atlas (TCGA) portal, including entries from both cancer tissues (n=179) and normal tissues (n=4). Meanwhile, microarray expression data were obtained from Gene Expression Omnibus (GEO) database. Specifically, we downloaded the Series Matrix File of dataset GSE25724, which contains expression profiles for 13 samples: 7 from healthy controls and 6 from diabetic patients. The annotation platform for this dataset is GPL96. Furthermore, we obtained single-cell RNA sequencing data from GEO dataset GSE235449 to support downstream single-cell level analyses. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

WGCNA analysis

WGCNA was employed to identify clusters (modules) of highly correlated genes, examine the relationship between gene modules and clinical phenotypes, and determine potential hub genes within key modules. A gene co-expression network was established from the PAAD expression dataset using the WGCNA package in R. To enhance computational efficiency and biological relevance, the top 5,000 genes ranked by variance were used to construct the network. A soft-thresholding power of 6 was chosen to achieve scale-free topology. The resulting adjacency matrix was transformed into a Topological Overlap Matrix (TOM), which provides a robust measure of network connectivity by accounting for shared neighbors between genes. Hierarchical clustering was then applied to the TOM to produce a gene dendrogram, where individual branches correspond to gene modules and are assigned distinct colors. Genes exhibiting comparable expression patterns were clustered into the same module, allowing the complex expression matrix to be simplified into biologically meaningful clusters. This modular organization facilitates the identification of gene groups that may be functionally relevant to PDAC pathogenesis.

Differential expression profiling

Differential gene expression was analyzed with the Limma package in R. Comparisons were made between control and disease groups to identify genes exhibiting significant expression differences. The criteria for DEG selection were defined as P<0.05 and |log₂ fold change (logFC)| >0.585. Volcano plot and heatmap were constructed to visualize the distribution and expression patterns of the identified DEGs.

Functional enrichment analysis

To evaluate the biological role of the DEGs, functional annotation was applied through Metascape (www.metascape.org). Gene Ontology (GO) was applied to assess the biological processes, cellular components, and molecular functions associated with the DEGs. Enrichment terms were considered statistically significant if the minimum gene overlap was ≥3 and the P value was ≤0.01.

Random survival forest analysis

To prioritize prognostically relevant genes, feature selection was carried out using the randomForestSRC package in R. A random survival forest model was implemented with 1,000 Monte Carlo iterations (nrep =1,000) to assess the relative importance of each gene. Genes with a relative importance score >0.3 were chosen as key prognostic markers for subsequent analyses.

Estimation of immune cell infiltration

Single-sample gene set enrichment analysis (ssGSEA) was employed to assess immune infiltration across samples. This approach allows for the quantification of immune cell type enrichment within individual expression profiles. A total of 29 immune cell phenotypes, including subsets of T cells, B cells, and natural killer (NK) cells, were evaluated. The ssGSEA scores were used to infer the relative abundance of each immune cell type within the tumor microenvironment.

Regulatory network analysis of key genes

To investigate the transcriptional regulatory mechanisms of key genes, we employed the RcisTarget package in R to estimate transcription factors according to motif enrichment analysis. All computations within RcisTarget are motif-centric, and the normalized enrichment score (NES) of each motif reflects its relative importance within the gene set, accounting for the total motif content in the database. Beyond motifs that have direct annotations, further annotations were inferred based on motif similarity and gene sequence context. Initially, the area under the curve (AUC) was calculated for every motif-gene set pair, representing the recovery of the gene set across a motif-ranked list. NES values were subsequently derived by normalizing these AUCs relative to the background distribution of all motifs. For this study, the RcisTarget.hg19.motifDBs.cisbpOnly.500bp database was employed for motif-to-gene mapping within a 500-bp upstream window. This approach enabled the identification of putative transcriptional regulators associated with the key genes.

Drug sensitivity analysis

To determine the potential therapeutic relevance of the identified key genes, drug sensitivity prediction was conducted using the pRRophetic package in R, using data from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). A ridge regression model is applied by the package to predict the half-maximal inhibitory concentration (IC₅₀) for each chemotherapeutic drug across individual tumor samples. Model performance was evaluated using 10-fold cross-validation on the GDSC training set to ensure robust prediction accuracy. Default parameters were used throughout the analysis, including the implementation of the “combat” method for batch effect correction and averaging of duplicate gene expression values. This analysis may aid in assessing drug response heterogeneity and potential gene-drug interactions, but its predictive results require further experimental validation.

Gene set enrichment analysis (GSEA)

To investigate the functional implications of gene expression variability, GSEA was implemented by stratifying samples into groups with high or low expression depending on target gene levels. GSEA was conducted using annotated datasets of genes accessed via MSigDB, specifically curated for subtype-specific signaling pathways. The analysis identified statistically significant enrichment of pathways by comparing gene expression distributions between the two groups. Enrichment was considered significant for gene sets with an adjusted P value <0.05, and outcomes were ranked based on NES and consistency across gene sets. GSEA is a commonly adopted approach for interpreting gene expression data in the terms of predefined biological pathways, which provides insights into molecular mechanisms underlying disease phenotypes.

Gene set variation analysis (GSVA)

GSVA is a non-parametric, unsupervised approach for estimating variation in gene set enrichment across samples in a transcriptomic dataset. Unlike traditional enrichment analyses, GSVA converts gene expression data from a focus on individual genes to a pathway-oriented perspective, allowing the assessment of biological process activity at the individual sample level. Herein, curated gene sets were retrieved from the MSigDB v7.0. Enrichment scores for each gene set across all samples were calculated using the GSVA algorithm. This allowed for the identification of functional pathway differences between conditions, providing a refined understanding of the molecular heterogeneity associated with disease progression.

Construction of the nomogram model

A nomogram was developed according to multivariate regression analysis to predict clinical outcomes by integrating gene expression levels and clinical parameters. Each predictor variable in the model was assigned a point value proportional to its regression coefficient, representing its relative contribution to the outcome. These point values were plotted as scaled line segments on a common plane. The cumulative score, calculated by summing the individual scores for each variable, was then used to estimate the predicted probability of the outcome. This visual tool allows for individualized risk estimation and aids in clinical decision-making.

Single-cell RNA sequencing analysis

Single-cell transcriptomic analysis was performed using the Seurat R package. Initially, low-expression genes were filtered out to improve data quality. The expression data were subsequently normalized and scaled. t-distributed stochastic neighbor embedding (t-SNE) was used for dimensionality reduction to visualize the spatial distribution and relationships between cell clusters. Reference data from the Human Primary Cell Atlas (HPCA) were used to annotate cell types in each cluster, enabling the identification of cell populations potentially relevant to disease onset and progression.

Cell lines and transfection

MIA PaCa-2 human pancreatic cancer cells (RRID: CVCL_0428) were obtained from ATCC (CRL-1420; lot/catalog No. BFN60803981). MIA PaCa-2 cells (human pancreatic cancer cells) were seeded in DMEM complete medium and cultured in a 37 ℃, 5% CO₂ incubator, with fresh medium replaced every 2 days. Log-phase MIA PaCa-2 cells (1×10⁵ cells/well) were digested with trypsin in a six-well plate. PCDH-GBP2 and PCDH-vector were transfected separately according to the Lipofectamine 3000 (Invitrogen) protocol. All in vitro validation experiments using commercially available cell lines (MIA PaCa-2) were conducted in accordance with institutional biosafety regulations.

Flow cytometry for cell cycle analysis

Forty-eight hours after transfection, MIA PaCa-2 cells were harvested, trypsinized, and fixed in pre-chilled 70% ethanol at 4 ℃ overnight. Fixed cells underwent washing with phosphate-buffered saline (PBS) and incubated in a staining solution containing 50 µg/mL propidium iodide (PI) and 100 µg/mL RNase A for 30 min in the dark at room temperature. Cell cycle distribution was evaluated via flow cytometry, and the data were analyzed to examine cell proportion in G0/G1, S, and G2/M phases.

Transwell invasion assay

Cell invasion ability was evaluated using Transwell chambers pre-coated with Matrigel to simulate the extracellular matrix. MIA PaCa-2 cells (1×10⁵ cells/mL), transfected either with PCDH-GBP2 or the empty control vector, were seeded into the upper chamber in serum-free medium. The bottom chamber was filled with medium containing 10% fetal bovine serum (FBS) as a chemoattractant. After a 24-hour incubation at 37 ℃, non-invading cells were gently removed from the upper membrane surface, while invading cells on the lower surface were fixed with paraformaldehyde (4%), followed by staining with crystal violet, and counting in five randomly selected fields under a light microscope.

Wound healing assay

Cell migration was assessed using a wound healing assay. MIA PaCa-2 cells were transfected with the PCDH-GBP2 plasmid for GBP2 overexpression or an empty PCDH vector as a control. Once the cell monolayer reached approximately 90% confluence, a linear scratch was created using a sterile 200 µL pipette tip. Non-adherent cells were removed by PBS wash, and adherent cells were cultured in 1% FBS medium to minimize proliferation. Wound closure was monitored by capturing images at 0, 24, and 48 hours using an inverted microscope. The width of the scratch was assessed through ImageJ, and the migration distance was calculated to quantify cell motility.

Statistical analysis

All statistical analyses were conducted using R software (version 4.2). Unless otherwise specified, a P value <0.05 was considered statistically significant. All in vitro assays were conducted with three biological replicates, and each biological replicate included three technical replicates.

Results

Identification of necroptosis-associated gene modules via WGCNA

We retrieved the PAAD gene expression dataset from TCGA database, comprising 183 samples, including 4 normal tissue and 179 tumor specimens. To investigate the role of necroptosis in PDAC, we obtained a list of necroptosis-related genes from GeneCards, and selected the top 200 genes ranked by Relevance Score for further analysis. ssGSEA was employed to calculate a sample-specific quantitative score for necroptosis according to the expression of these genes. Using these scores, WGCNA was used to identify gene modules related to necroptosis activity in PAAD. A soft thresholding power (β=6) was applied to maintain scale-free characteristics (Figure 1A,1B), and a TOM was used to detect gene co-expression modules. Overall, 10 distinct gene modules were identified (Figure 1C), color-coded as follows: black (298 genes), blue [1,297], brown [1,638], cyan [128], gray [77], lightcyan [71], magenta [256], purple [369], red [374], and yellow [492]. Correlation analysis between module eigengenes and necroptosis scores revealed that the lightcyan module demonstrated the strongest positive correlation (correlation coefficient =0.64, P=1e−22), indicating its potential involvement in necroptosis-related regulatory mechanisms in PDAC (Figure 1D). This module contains a total of 71 genes and was selected for subsequent analysis.

Figure 1 Construction and analysis of WGCNA network for PAAD based on necroptosis-related genes. (A) Scale-free topology model fitting plot, showing the relationship between soft threshold power and scale independence, with the selected soft threshold power β set at 6. (B) Mean connectivity plot, illustrating the average connectivity of genes at different soft threshold powers. (C) Gene dendrogram and module colors, generated by hierarchical clustering and dynamic tree cut, with 10 distinct gene modules identified in the data, namely black (298 genes), blue (1,297 genes), brown (1,638 genes), cyan (128 genes), grey (77 genes), lightcyan (71 genes), magenta (256 genes), purple (369 genes), red (374 genes), and yellow (492 genes). (D) Module-trait relationships, showing the correlation between each module and necroptosis, with the lightcyan module exhibiting the highest correlation (cor =0.64, P=1e−22). PAAD, pancreatic adenocarcinoma; WGCNA, weighted gene co-expression network analysis.

Identification of shared genes between T2DM and PDAC

A total of 13 samples from the GSE25724 dataset related to T2DM were obtained, comprising 7 control and 6 disease specimens. Differential gene expression analysis was conducted at P<0.05 and |log₂FC| >0.585. This analysis yielded 2,471 DEGs, involving 997 upregulated and 1,474 downregulated genes (Figure 2A,2B). To explore the shared molecular mechanisms between T2DM and PDAC, these DEGs were overlapped with the genes in the lightcyan module identified in the WGCNA analysis of the PAAD dataset (Figure 2C), resulting in 17 overlapping genes. Functional enrichment analysis of these 17 genes indicated obvious involvement in pathways related to molecular function inhibition and the regulation of the innate immune response (Figure 2D). These 17 genes were subsequently used as a candidate gene set for downstream analyses.

Figure 2 Identification and analysis of DEGs and their intersection with WGCNA module genes in type 2 diabetes. (A) Volcano plot showing the distribution of DEGs between control and disease groups, with the horizontal axis representing log₂FC and the vertical axis representing −log₁₀ (P value). The red dots indicate upregulated genes, the blue dots indicate downregulated genes, and the gray dots represent non-significant genes. (B) Bar plot showing the number of upregulated and downregulated DEGs, with a total of 2,471 DEGs identified, including 997 upregulated genes and 1,474 downregulated genes. (C) Venn diagram showing the intersection of genes from the lightcyan module in the WGCNA network and DEGs from the type 2 diabetes dataset, with 17 overlapping genes identified. (D) Bubble plot showing the enrichment analysis of the 17 overlapping genes in various pathways, with significant enrichment in pathways such as molecular function inhibitor activity and regulation of innate immune response. DEGs, differentially expressed genes; FC, fold change; GO, Gene Ontology; WGCNA, weighted gene co-expression network analysis.

Identification of key genes via random survival forest analysis

To further pinpoint key genes implicated in both PDAC and T2DM, we conducted a random survival forest analysis on the 17 intersecting genes. Genes with a relative importance score >0.3 were chosen as candidate markers, identifying four genes of high importance (Figure 3A,3B). Subsequent survival analysis of these four genes (Figure 3C,3D) revealed that expression levels of two genes were markedly associated with patient survival (P<0.05). These two genes, GBP2 and LY6E, were prioritized for further investigation.

Figure 3 Identification of key genes influencing PDAC and type 2 diabetes through random survival forest analysis and survival analysis. (A) Variable relative importance plot generated by random survival forest analysis, showing the relative importance of the 17 overlapping genes. (B) Genes with relative importance greater than 0.3 were identified as final markers, and the importance order of the top 4 genes is displayed. (C) Kaplan-Meier survival curve for GBP2, showing the survival probability over time for high and low expression groups. The log-rank test indicates a significant difference in survival (P=0.02). (D) Kaplan-Meier survival curve for LY6E, showing the survival probability over time for high and low expression groups. The log-rank test indicates a significant difference in survival (P=0.007). GBP2 and LY6E were identified as key genes for subsequent research due to their significant survival differences. PDAC, pancreatic ductal adenocarcinoma.

Association of GBP2 and LY6E with immune cell infiltration and the tumor microenvironment in PDAC

The tumor microenvironment profoundly influences disease diagnosis and therapeutic response. To explore the molecular mechanisms by which the identified genes modulate PDAC progression, the association between these genes and immune cell infiltration within the PDAC dataset was analyzed. Our results reveal the proportions of various immune cell types across samples and characterize their interrelationships (Figure 4A,4B). Notably, obvious differences were noted between tumor and control groups in immune features such as APC co-stimulation, immune checkpoint molecules, HLA expression, immature dendritic cells (iDCs), neutrophils, T helper cells, and regulatory T cells (Tregs) (Figure 4C). Further correlation analyses demonstrated strong relationships between the key genes and specific immune cells (Figure 4D). Specifically, GBP2 showed significant positive correlations with parainflammation and major histocompatibility complex (MHC) class I expression, while LY6E was positively correlated with type I interferon response and parainflammation. Further analysis using the TISIDB database identified significant correlations of these genes with diverse immune regulators, chemokines, and receptors (Figure S1). Collectively, these findings suggest that GBP2 and LY6E closely correlate with immune cell infiltration and are crucial for shaping the PDAC immune landscape.

Figure 4 Analysis of the relationship between key genes and immune infiltration in PDAC. (A) Heatmap showing the proportion of immune cells in each sample, with rows representing different immune cell types and columns representing individual samples. The color intensity indicates the relative abundance of each immune cell type in the samples. (B) Correlation heatmap showing the relationships between different immune cell types, with the color scale representing the correlation coefficients. (C) Bar plot showing the differences in immune cell infiltration between normal and tumor groups, with significant differences observed in APC_co_stimulation, Check-point, HLA, iDCs, Neutrophils, T_helper_cells, and Treg (indicated by asterisks). (D) Scatter plot showing the correlation between the key gene GBP2, LY6E and various immune cell types. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. aDC, activated dendritic cell; APC, antigen-presenting cell; CCR, creatinine clearance rate; DCs, dendritic cells; HLA, Human Leukocyte Antigen; iDC, immature dendritic cell; MHC, major histocompatibility complex; NK, natural killer; PDAC, pancreatic ductal adenocarcinoma; pDCs, plasmacytoid dendritic cells; Treg, regulatory T cells.

Chemotherapy sensitivity analysis of key genes GBP2 and LY6E

The combined treatment of surgery and chemotherapy is well-established for early-stage PDAC. In this study, we employed the R package “pRRophetic” to predict the chemotherapy response of individual tumor samples. We further investigated the potential association between the expression of GBP2 and LY6E and sensitivity to commonly used chemotherapeutic agents. Our results indicated that GBP2 expression was may correlated with sensitivity to A.443654, AMG.706, A.770041, and CCT007093 (Figure 5A). Similarly, LY6E expression was showed a potential associated with sensitivity to ABT.263, A.443654, ATRA, AMG.706, A.770041, and CCT007093 (Figure 5B).

Figure 5 Drug sensitivity analysis. (A) Violin plots depicting the distribution of estimated IC₅₀ values for A.443654, AMG.706, A.770041, CCT007093, ABT.263 and ATRA in relation to GBP2 expression levels, categorized into LExp and HExp groups. (B) Violin plots illustrating the distribution of estimated IC₅₀ values for ABT.263, A.443654, ATRA, AMG.706, A.770041, and CCT007093 in relation to LY6E expression levels. HExp, high expression; IC₅₀, half-maximal inhibitory concentration; LExp, low expression.

Pathway enrichment analysis of key genes GBP2 and LY6E

To further investigate the signal transduction cascades associated with the key genes and to elucidate their molecular mechanisms underlying disease progression, GSEA and GSVA analyses were performed. GSEA revealed that GBP2 was significantly enriched in the B cell receptor and tumor necrosis factor (TNF) pathways (Figure 6A,6B). For LY6E, enrichment was observed in the interleukin (IL)-17, nucleotide-binding oligomerization domain (NOD)-like receptor, and other pathways (Figure 6C,6D). GSVA analysis further demonstrated that GBP2 was enriched in the P53 pathway and IL2/STAT5 axis (Figure 6E), while LY6E was enriched in Wnt/β-catenin and mTORC1 pathways (Figure 6F). These findings imply that GBP2 and LY6E may influence PDAC progression via distinct but critical immune and oncogenic signaling networks.

Figure 6 Pathway enrichment analysis of the two key genes using GSEA and GSVA. (A) GSEA enrichment plot showing that GBP2 is significantly associated with the B cell receptor signaling pathway, with a NES indicating the strength of the association. (B) GSEA enrichment plot showing that GBP2 is significantly associated with the TNF signaling pathway, highlighting its potential role in immune regulation. (C) GSEA enrichment plot showing that LY6E is significantly associated with the IL-17 signaling pathway, suggesting its involvement in inflammatory responses. (D) GSEA enrichment plot showing that LY6E is significantly associated with the NOD-like receptor signaling pathway, indicating its role in innate immunity. (E) GSVA analysis showing that GBP2 is associated with pathways such as P53_PATHWAY and IL2_STAT5_SIGNALING, with t-values comparing HExp and LExp groups. (F) GSVA analysis showing that LY6E is associated with pathways such as WNT_BETA_CATENIN_SIGNALING and MTORC1_SIGNALING, with t-values comparing HExp and LExp groups. These analyses provide insights into the potential molecular mechanisms through which the key genes influence disease progression. GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; IL, interleukin; NES, normalized enrichment score; NOD, nucleotide-binding oligomerization domain; TNF, tumor necrosis factor.

Construction and evaluation of a nomogram based on core gene expression

To integrate the prognostic significance of the key genes with clinical indicators, we constructed a nomogram based on multivariate regression data. This model includes the expression levels of GBP2 and LY6E, along with clinical parameters, to predict overall survival (OS) in PDAC patients. The nomogram visually represents the relative contribution of each variable to the total prognostic score (Figure S2A). Furthermore, we evaluated the model’s predictive accuracy through a comparison of the predicted and observed OS at 1 and 3 years. Calibration curves demonstrated strong concordance between predicted and actual outcomes, indicating that the nomogram possesses robust predictive accuracy for patient survival (Figure S2B).

Single-cell expression analysis of key genes in the tumor microenvironment

We retrieved the GSE235449 dataset from the NCBI GEO database to assess the cell-type-specific expression patterns of the key genes GBP2 and LY6E within the pancreatic tumor microenvironment. The expression levels of these genes were analyzed across seven major immune and stromal cell types, including T cells, monocytes, macrophages, B cells, dendritic cells (DCs), NK cells, and tissue stem cells (Figure 7). To further explore potential interactions between key genes and tumor progression pathways, we identified three well-established PDAC progression-related genes, such as ATM, CCND1, and MET, from GeneCards. Co-expression analyses were conducted between these progression genes and the two key genes across the seven cell types, and the results were visualized to highlight potential cellular and molecular associations (Figure S3).

Figure 7 Single-cell RNA sequencing analysis of key gene expression across various immune cell types in the GSE235449 dataset. (A) t-SNE plot showing the expression levels of GBP2 and LY6E across different cell types, including T cells, monocytes, macrophages, B cells, DCs, NK cells, and tissue stem cells. (B) t-SNE plot showing the expression levels of GBP2 across the same cell types. (C) t-SNE plot showing the expression levels of LY6E across the same cell types. Similar to panel A, the color gradient represents the average expression levels of GBP2 and LY6E, and the dot size indicates the percentage of cells expressing GBP2 and LY6E. DCs, dendritic cells; NK, natural killer; t-SNE, t-distributed stochastic neighbor embedding.

Differential expression analysis of key genes in pancreatic cancer patients and diabetic conditions

To investigate whether key genes exhibit expression differences associated with diabetic status in pancreatic cancer patients, we analyzed clinical information from the TCGA-PAAD database. Among these patients, 38 had a history of diabetes, while 111 did not. Box plots were generated by grouping samples based on the presence of T2DM to visualize expression differences between groups. As shown in Figure S4, expression levels of LY6E and GBP2 in pancreatic cancer patients fluctuated between the diabetic (yes) and non-diabetic (no) groups but did not exhibit statistically significant differences.

Functional validation of GBP2 in PDAC cell migration, invasion, and cell cycle regulation

To evaluate the functional role of GBP2 in PDAC, we investigated the effects of its overexpression on the migratory, invasive, and proliferative properties of MIA PaCa-2 cells. GBP2 was overexpressed via transfection with the PCDH-GBP2 construct, and corresponding vector controls were used for comparison. Wound healing assays demonstrated a significantly reduced closure of the scratch area in the PCDH-GBP2 group compared to the PCDH-vector group at both 24 and 48 hours (Figure 8A). Quantitative analysis confirmed that GBP2 overexpression significantly reduced cell migration distance (P<0.05), suggesting an inhibitory effect on migratory capacity. Transwell invasion assays further supported these findings (Figure 8B), exhibiting a notable reduction in invading cell counts in the GBP2-overexpressing group compared to controls (P<0.01). This indicates that GBP2 also suppresses the invasive potential of MIA PaCa-2 cells. To investigate the potential mechanism by which GBP2 affects proliferation, we conducted flow cytometry-based cell cycle analysis (Figure 8C). Overexpression of GBP2 led to a significant cell deposition into the G0/G1 phase and a reduction in the G2/M and S phases compared to vector controls (P<0.05), indicating that GBP2 induces G0/G1 phase cell cycle arrest. Altogether, these results demonstrate that GBP2 overexpression inhibits migration and invasion, and impairs cell cycle progression in pancreatic cancer cells, suggesting that GBP2 holds potential as a therapeutic target in pancreatic cancer treatment.

Figure 8 Analysis of the effects of GBP2 overexpression on cell migration, invasion, and cell cycle progression in MIA PaCa-2 cells. (A) Representative images of the wound healing assay at 24 and 48 hours post-wounding. The scratch area in the PCDH-GBP2 group was significantly less filled compared to the PCDH-vector group, indicating reduced migration distance (P<0.05). Scale bar =100 μm. (B) Representative images of the transwell invasion assay. The number of invading cells in the PCDH-GBP2 group was significantly lower than in the PCDH-vector group (P<0.01), demonstrating that GBP2 overexpression markedly reduces the invasive capacity of MIA PaCa-2 cells. Cells were fixed with 4% paraformaldehyde for 15 min and stained with 0.1% crystal violet for 20 min; invading cells were counted in ten randomly selected fields under a light microscope at 200× magnification; scale bar =100 μm. (C) Cell cycle analysis by flow cytometry. Compared to the PCDH-vector group, the PCDH-GBP2 group exhibited an increased proportion of cells in the G0/G1 phase (P<0.05) and decreased proportions in the S and G2/M phases (P<0.05). This suggests that GBP2 overexpression induces cell cycle arrest in the G0/G1 phase, thereby inhibiting the proliferation of MIA PaCa-2 cells.

Discussion

Diabetes is a heterogeneous metabolic condition defined by chronic hyperglycemia caused by impaired insulin production, insulin action, or both mechanisms, leading to peripheral insulin resistance (19). Diabetes is associated with a significantly increased risk of complications, including stroke, blindness, limb amputation, cardiovascular disease, kidney failure, and various forms of cancer (20). Diabetes is commonly classified into three primary types: type 1 diabetes mellitus, T2DM, and gestational diabetes mellitus (21). Another less common form, known as type 3c diabetes mellitus, is associated with exocrine PDAC (21). Diabetes may serve as both a risk factor for PDAC and a consequence of it (22). Notably, a subset of PDAC patients are diagnosed with diabetes within 3 to 5 years prior to their cancer diagnosis. This form of newly-onset diabetes is referred to as PCRD, which is increasingly recognized as a potential early marker or harbinger of PDAC (23). The pathogenesis of PCRD is thought to involve several mechanisms, including tumor-secreted mediators that induce β-cell dysfunction and peripheral insulin resistance, altered interactions between pancreatic tumors and adipose tissue, dysregulation of diabetes-associated oncogenic pathways via microRNAs, and enhanced transforming growth factor-beta (TGF-β) signaling that promotes β-cell exhaustion through increased TGF-β secretion (23). Clinically, PCRD presents with features that are distinct from those of T2DM, suggesting a different underlying pathophysiology (24). Due to the bidirectional relationship between PDAC and DM, the association between the two conditions is complex. Distinguishing PCRD from T2DM is critical, especially in new-onset diabetes, as it may serve as an early indicator of occult PDAC. Therefore, the identification of novel, clinically applicable biomarkers with high sensitivity and specificity is essential for the early detection of PDAC, which could significantly improve therapeutic outcomes. Emerging evidence suggests that early-stage pancreatic neoplasms, such as pancreatic intraepithelial neoplasia, may secrete bioactive molecules within exosomes, small extracellular vesicles, which disrupt islet function and peripheral glucose metabolism (25). These tumor-derived exosomes may contribute to β-cell dysfunction and insulin resistance, positioning them as a potential mechanistic link between early tumorigenesis and the development of PCRD. Accordingly, the role of exosomal cargo in mediating metabolic reprogramming has become a growing area of interest in PCRD research.

Our WGCNA analysis identified the lightcyan module as being strongly associated with necroptosis in PDAC (correlation =0.64, P=1e−22). Cross-referencing the genes within this module with T2DM-related DEGs revealed 17 overlapping genes. Among these, GBP2 and LY6E emerged as key candidates due to their significant prognostic value in survival analyses. These genes were enriched in immune-related pathways, such as TNF and IL-17 signaling, as well as metabolic and oncogenic pathways, including Wnt/β-catenin and mTORC1 signaling. These results imply a dual role for GBP2 and LY6E in both tumor progression and dysregulation of glucose metabolism.

Overexpression of GBP2 in MIA PaCa-2 cells significantly inhibited cell migration, invasion, and cell cycle progression (P<0.05), aligning with its previously reported tumor-suppressive functions in colorectal and glioblastoma cancers through modulation of the WNT and STAT3 signaling pathways (26,27). In contrast, LY6E (a pro-oncogenic member of the LY6 family) was found to promote tumor cell adhesion and epithelial-mesenchymal transition via activating the TGF-β1/Smad pathways (28). These findings highlight the opposing roles of GBP2 and LY6E in PDAC, suggesting their context-dependent functions within the TME.

Immune infiltration analysis indicated remarkable correlations between GBP2/LY6E expression and immune cells, including Tregs, neutrophils, and MHC class I molecules. GBP2 was strongly associated with parainflammatory responses and antigen presentation, while LY6E was related to type I interferon signaling (P<0.01). This immune dysregulation may contribute both to insulin resistance in T2DM and to the establishment of an immunosuppressive TME in PDAC, potentially forming a bidirectional feedback loop (24). Drug sensitivity profiling further identified GBP2 and LY6E as potential biomarkers of chemotherapeutic response. Elevated GBP2 expression correlated with increased sensitivity to A-443654 and AMG-706, whereas LY6E expression was associated with responsiveness to ABT-263 and ATRA (P<0.05). These results indicate the possible clinical applicability of GBP2 and LY6E as predictive markers for personalized therapy, particularly in PDAC patients with comorbid diabetes.

The nomogram model incorporating GBP2, LY6E, and clinical variables demonstrated strong predictive performance for 1- and 3-year OS, with a C-index of 0.78, indicating its utility as a practical tool for patient risk stratification. Notably, single-cell RNA-seq analysis showed elevated levels of GBP2 and LY6E in macrophages and T cells, implicating these genes in immune evasion and stromal remodeling. These findings are consistent with emerging evidence that tumor-derived exosomes contribute to PCRD, in which exosomal cargo may disrupt β-cell function and exacerbate metabolic dysfunction (29).

Despite offering novel insights, there are several limitations in this study. First, the limited sample size of the T2DM dataset (GSE25724, n=13) may restrict the generalizability of the findings. Second, experimental validation was restricted to in vitro models; further studies employing genetically engineered mouse models or patient-derived xenografts are warranted to elucidate the underlying mechanisms. Third, the functional involvement of GBP2 and LY6E in exosome-mediated communication between pancreatic cancer cells and islet β-cells remains to be investigated. Furthermore, the drug sensitivity analysis in this study is based on computational prediction models. While the results suggest a potential association between gene expression and chemotherapy response, they have not yet been directly validated through in vitro or in vivo drug sensitivity experiments. Therefore, further experimental confirmation is required before clinical translation.

Conclusions

In conclusion, GBP2 and LY6E appear to serve as key molecular mediators linking PDAC and T2DM, influencing immune cell infiltration, metabolic reprogramming, and treatment response. Their paradoxical roles in both tumor suppression and promotion suggest the complexity of PDAC-T2DM comorbidity. Future work should prioritize the validation of these biomarkers in large patient populations and explore their roles in exosomal signaling, with the ultimate goal of advancing early detection and precision therapies for high-risk patient populations.

Acknowledgments

The authors thank the Department of Nuclear Medicine, The Third Affiliated Hospital of Sun Yat-sen University, for providing technical assistance and scientific guidance.

Footnote

Reporting Checklist: The authors have completed the MDAR and TRIPOD reporting checklists. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2427/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2427/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2427/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All in vitro validation experiments using commercially available cell lines (MIA PaCa-2) were conducted in accordance with institutional biosafety regulations.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Khasawneh B, Esmail A, Al-Najjar E, et al. An 81-Year-Old Geriatric Patient with Metastatic Pancreatic Cancer Demonstrating Excellent Response and Well Tolerance to NALIRIFOX: A Case Report and Literature Review. Reports 2025;8:69. (MDPI). [Crossref] [PubMed]
2. Rahib L, Wehner MR, Matrisian LM, et al. Estimated Projection of US Cancer Incidence and Death to 2040. JAMA Netw Open 2021;4:e214708. [Crossref] [PubMed]
3. Demcsák A, Sahin-Tóth M. Misfolding CPA1 mutation accelerates precancerous pancreas lesions in KC mice. Am J Physiol Gastrointest Liver Physiol 2025;329:G546-56. [Crossref] [PubMed]
4. Byeon S, McKay MJ, Molloy MP, et al. Novel serum protein biomarker panel for early diagnosis of pancreatic cancer. Int J Cancer 2024;155:365-71. [Crossref] [PubMed]
5. Zhang T, Moore SC, Fu S, et al. Association between prediagnostic serum metabolites and pancreatic ductal adenocarcinoma risk in two prospective cohorts. Int J Cancer 2025;157:1304-15. [Crossref] [PubMed]
6. Sapoor S, Nageh M, Shalma NM, et al. Bidirectional relationship between pancreatic cancer and diabetes mellitus: a comprehensive literature review. Ann Med Surg (Lond) 2024;86:3522-9. [Crossref] [PubMed]
7. Hill JLE, Hill TG, Parslow D, et al. Pancreatic Cancer and Diabetes: Insights, Hypotheses, and Next Steps. Int J Mol Sci 2025;26:10245. [Crossref] [PubMed]
8. Wayne CD, Benbetka C, Besner GE, et al. Challenges of Managing Type 3c Diabetes in the Context of Pancreatic Resection, Cancer and Trauma. J Clin Med 2024;13:2993. [Crossref] [PubMed]
9. ElSayed NA, Aleppo G, Aroda VR, et al. 2. Classification and Diagnosis of Diabetes: Standards of Care in Diabetes-2023. Diabetes Care 2023;46:S19-40.
10. Sharma A, Kandlakunta H, Nagpal SJS, et al. Model to Determine Risk of Pancreatic Cancer in Patients With New-Onset Diabetes. Gastroenterology 2018;155:730-739.e3. [Crossref] [PubMed]
11. Diagnosis and classification of diabetes mellitus. Diabetes Care 2010;33:S62-9.
12. Pelaez-Luna M, Takahashi N, Fletcher JG, et al. Resectability of presymptomatic pancreatic cancer and its relationship to onset of diabetes: a retrospective review of CT scans and fasting glucose values prior to diagnosis. Am J Gastroenterol 2007;102:2157-63. [Crossref] [PubMed]
13. Paternoster S, Falasca M. The intricate relationship between diabetes, obesity and pancreatic cancer. Biochim Biophys Acta Rev Cancer 2020;1873:188326. [Crossref] [PubMed]
14. Yan S, Sun M, Gao L, et al. Identification of Key LncRNAs and Pathways in Prediabetes and Type 2 Diabetes Mellitus for Hypertriglyceridemia Patients Based on Weighted Gene Co-Expression Network Analysis. Front Endocrinol (Lausanne) 2021;12:800123. [Crossref] [PubMed]
15. Wang S, Tian J, Wang J, et al. Identification of the Biomarkers and Pathological Process of Heterotopic Ossification: Weighted Gene Co-Expression Network Analysis. Front Endocrinol (Lausanne) 2020;11:581768. [Crossref] [PubMed]
16. Yuan C, Babic A, Khalaf N, et al. Diabetes, Weight Change, and Pancreatic Cancer Risk. JAMA Oncol 2020;6:e202948. [Crossref] [PubMed]
17. Yang W, Jiang W, Guo S. Regulation of Macronutrients in Insulin Resistance and Glucose Homeostasis during Type 2 Diabetes Mellitus. Nutrients 2023;15:4671. [Crossref] [PubMed]
18. Pathak AK, Tiwari V, Kulshrestha MR, et al. Impact of essential metals on insulin sensitivity and fetuin in obesity-related type 2 diabetes pathogenesis. J Diabetes Metab Disord 2023;22:703-12. [Crossref] [PubMed]
19. Cloete L. Diabetes mellitus: an overview of the types, symptoms, complications and management. Nurs Stand 2022;37:61-6. [Crossref] [PubMed]
20. Kim H, Lee H. Risk of Stroke and Cardiovascular Disease According to Diabetes Mellitus Status. West J Nurs Res 2023;45:520-7. [Crossref] [PubMed]
21. Krause M, De Vito G. Type 1 and Type 2 Diabetes Mellitus: Commonalities, Differences and the Importance of Exercise and Nutrition. Nutrients 2023;15:4279. [Crossref] [PubMed]
22. Roy A, Sahoo J, Kamalanathan S, et al. Diabetes and pancreatic cancer: Exploring the two-way traffic. World J Gastroenterol 2021;27:4939-62. [Crossref] [PubMed]
23. Park MN. Therapeutic Strategies for Pancreatic-Cancer-Related Type 2 Diabetes Centered around Natural Products. Int J Mol Sci 2023;24:15906. [Crossref] [PubMed]
24. Hu Y, Zeng N, Ge Y, et al. Identification of the Shared Gene Signatures and Biological Mechanism in Type 2 Diabetes and Pancreatic Cancer. Front Endocrinol (Lausanne) 2022;13:847760. [Crossref] [PubMed]
25. Zhao C, Gao F, Weng S, et al. Pancreatic cancer and associated exosomes. Cancer Biomark 2017;20:357-67. [Crossref] [PubMed]
26. Wang J, Min H, Hu B, et al. Guanylate-binding protein-2 inhibits colorectal cancer cell growth and increases the sensitivity to paclitaxel of paclitaxel-resistant colorectal cancer cells by interfering Wnt signaling. J Cell Biochem 2020;121:1250-9. [Crossref] [PubMed]
27. Yu S, Yu X, Sun L, et al. GBP2 enhances glioblastoma invasion through Stat3/fibronectin pathway. Oncogene 2020;39:5042-55. [Crossref] [PubMed]
28. Zhang Q, Jia Y, Pan P, et al. α5-nAChR associated with Ly6E modulates cell migration via TGF-β1/Smad signaling in non-small cell lung cancer. Carcinogenesis 2022;43:393-404. [Crossref] [PubMed]
29. Binang HB, Perera CJ, Apte MV. Role of Pancreatic Tumour-Derived Exosomes and Their Cargo in Pancreatic Cancer-Related Diabetes. Int J Mol Sci 2023;24:10203. [Crossref] [PubMed]