All datasets were obtained from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). The microarray datasets GSE65146 and GSE109227, based on the GPL6246 platform, were used for transcriptome analysis of pancreatic tissues from AP and control mice. Batch effects were corrected, and DEGs between AP and CTRL were identified using the limma package in R (version 4.3.1). Genes were defined as differentially expressed if |(logFC)| > 1 and adjusted P-value < 0.05. Principal component analysis (PCA) was conducted to assess sample clustering and to evaluate the effectiveness of batch effect correction. The datasets GSE169076 (GPL23479) and GSE298193 (GPL25947) served as independent external validation cohorts in this study. The single-cell RNA-seq dataset GSE279876 (GPL19057) was used for cell-level resolution analysis. Detailed dataset characteristics and group assignments are provided in Supplementary Table 1 . To explore the biological functions and pathways associated with the DEGs, Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the clusterProfiler package in R. Terms or pathways with a P-value < 0.05 were considered significantly enriched.The overall workflow of this study is illustrated in Figure 1 .

Glycolysis-related genes were compiled from gene sets in the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/). After harmonizing gene symbols and removing duplicates ( Supplementary Table 2 ), Glycolysis-related genes were intersected with the DEGs identified in our dataset to create a subset for downstream analysis. Three complementary machine learning algorithms including LASSO regression, Boruta feature selection, and SVM-RFE. LASSO logistic regression was performed using the glmnet R package and 10-fold cross-validation to select features based on the optimal lambda values (lambda.min and lambda.1se). Boruta in the Boruta R package was utilized to assess feature importance over 50 iterations, and tentative features were refined using TentativeRoughFix(). SVM-RFE was conducted using the sigFeature package and a custom msvmRFE script with 10-fold cross-validation to select the most predictive genes iteratively. Final key genes were defined by the intersection of the results from all three methods and validated by ROC curve analysis in the validation dataset.

Gene Set Enrichment Analysis (GSEA) was performed using the “HALLMARK_GLYCOLYSIS” gene set obtained from the MSigDB database to explore the functional relevance of glycolysis-related phenotypes. The analysis was based on a ranked list of all expressed genes ordered by their differential expression statistics. GSEA was conducted using the clusterProfiler R package. Enrichment significance was evaluated by normalized enrichment scores (NES) and false discovery rate (FDR) values. To further investigate the functional roles of key glycolysis-related DEGs, we performed GSEA on GO biological process terms for each gene individually. Genes with an absolute Spearman correlation ≥ 0.8 and p-value < 0.05 were selected for downstream analysis. Based on NES, the top five significantly enriched upregulated and downregulated GO biological processes were identified.

Protein-protein interaction (PPI) networks were constructed to explore potential functional relationships among the glycolysis-related DEGs. Interaction data were obtained from the STRING database (https://cn.string-db.org/). Network visualization and analysis were carried out using Cytoscape (version 3.9.1). Hub genes were identified using the Maximal Clique Centrality (MCC) algorithm implemented in the CytoHubba plugin.

ScRNA-seq data from control and AP samples (GSE298193) were processed using Scanpy (Python 3.8). Cells with < 200 genes, < 1,000 UMIs, or > 5% mitochondrial reads were excluded to remove low-quality or dying cells, and genes detected in fewer than 3 cells were removed to reduce noise. Putative doublets were identified using Scrublet (31), which simulates artificial doublets from the observed data and classifies cells based on transcriptome similarity. Scrublet was run separately for each sample with automatic thresholding, and cells predicted as doublets were removed before downstream analyses. Counts were normalized with a delta-method shifted log transformation, and the 2,000 most highly variable genes (HVG) were selected via the Pearson-residual method. PCA was performed on the scaled HVG matrix (top 50 PCs retained), and Leiden clustering was applied to the PCA space (resolution = 0.25). Cluster identities were assigned manually based on marker genes curated from CellMarker 2.0 (http://117.50.127.228/CellMarker/) and PanglaoDB (https://panglaodb.se/index.html). For visualization, a modified distance embedding was computed from the PCA space.

The levels of immune cell infiltration were estimated using single-sample Gene Set Enrichment Analysis (ssGSEA) based on immune cell marker gene sets (32)( Supplementary Table 4 ). The enrichment scores were normalized to calculate the relative proportions of immune cells across samples. Differences in immune infiltration between groups were evaluated using Student’s t-test with multiple testing correction. Spearman correlation analysis was conducted to assess relationships among immune cell types as well as between immune infiltration and the key glycolysis-related DEGs. The results were visualized using heatmaps and boxplots to depict patterns of immune cell infiltration and their associations with core gene expression.

NetworkAnalyst (https://www.networkanalyst.ca/) was used to generate the regulatory networks. TF-gene interactions were derived from ENCODE ChIP-seq data applying the BETA Minus algorithm with a peak intensity signal threshold of <500 and a predicted regulatory potential score of <1. MiRNA–gene interactions were obtained from miRTarBase v9.0. The integrated TF–miRNA–mRNA regulatory relationships were visualized using Cytoscape.

C57BL/6 mice (6–8 weeks old, 18–24 g weight) were purchased from Jiangsu GemPharmatech Co., Ltd. and housed under specific pathogen-free conditions (22–24 °C, 50–60% humidity, 12 h light/dark cycle). The AP model was induced by 10 hourly intraperitoneal (i.p.) injections of caerulein (100 μg/kg; Tocris, UK, Cat No.6264). Lipopolysaccharide (LPS; 10mg/kg; Sigma, USA, Cat No.2880) was administered (i.p.) immediately after the 10th caerulein injection. Mice were sacrificed 12 h after LPS injection. Amylase and lipase in serum and myeloperoxidase (MPO) activity in pancreatic and lung tissues were measured. Specimens of pancreatic and lung tissue were fixed in 10% formalin at room temperature for 24 hours and embedded in paraffin. Sections were stained with hematoxylin and eosin (H&E). Pathological severity was semi-quantitatively scored based on edema, inflammatory infiltration, and acinar necrosis, each rated on a 0–3 scale, as previously described (33, 34). Lung injury assessment was conducted by quantifying edema and the accumulation of inflammatory cells in 10 randomly selected fields per section at 100× magnification, as described in previous studies (33). Whole-slide imaging was performed using the Olympus SLIDEVIEW VS200 system (Olympus, Japan), and Fiji software was used for digital image analysis.

Total RNA was extracted from pancreatic tissues using TRIzol reagent and was reverse transcribed into cDNA. qPCR was conducted using the Bio-Rad CFX96 Real-Time PCR System (Bio-Rad, USA) and SYBR Green Master Mix (Vazyme, China). Gene expression was quantified by the 2^^−ΔΔCt method, normalized to 18S rRNA. Reagents and primer sequences are provided in the Supplementary Materials .

Western blot and immunohistochemistry (IHC) were performed to detect CLDN3 and CD44 expression in pancreatic tissues following standard protocols. The detailed methods were described in the Supplementary Materials .

Statistical analyses were conducted using R software (version 4.3.1) and GraphPad Prism 9.5 (GraphPad Software Inc., USA). Normally distributed data were presented as mean ± standard deviation and were analyzed using Student’s t-test or one-way analysis of variance. Continuous data were analyzed using the Mann–Whitney U test. Correlations between variables were assessed using Pearson’s or Spearman’s correlation coefficients. All tests were two-tailed, and p < 0.05 was considered statistically significant.

After batch correction, normalized expression values were evenly distributed across the 10 Ctrl (Control) and 9 AP samples ( Figure 2A ; Supplementary Figure 1 ). Three-dimensional PCA confirmed a clear separation between AP and Ctrl samples ( Figure 2B ). Differential expression analysis using the limma package (|log₂FC| ≥ 1, adjusted P < 0.05) identified 1,943 DEGs, comprising 1,156 upregulated and 787 downregulated genes ( Figure 2C ). Hierarchical clustering of the top 50 DEGs clearly distinguished AP from Ctrl samples, further supporting the differential expression patterns ( Figure 2D ). GO enrichment revealed that these DEGs are predominantly involved in biological processes (BP) such as wound healing, actin-filament organization, and positive regulation of cytokine production; associated with cellular components (CC) including the cell–substrate junction, focal adhesion, and actin cytoskeleton; and exhibited molecular functions (MF) such as cadherin binding, actin binding, and protein serine/threonine/tyrosine kinase activity ( Figure 2E ). KEGG pathway analysis showed significant enrichment in the PI3K–AKT, and MAPK signaling pathways ( Figure 2F ). These pathways reflect key pathophysiological features of AP, including immune activation, inflammatory signaling, cellular stress responses, and tissue repair mechanisms.

GSEA revealed significant enrichment of the glycolysis pathway in AP samples, with the NES of 2.48 and FDR below 0.25 ( Figure 3A ). By intersecting 1,483 DEGs with 313 glycolysis-related DEGs using a Venn diagram, 43 glycolysis-related DEGs were identified ( Figure 3B ). Boxplots comparing expression levels of the 43 glycolysis-related DEGs between AP and Ctrl groups showed distinct expression differences ( Figure 3C ). The chromosomal distribution revealed that glycolysis-related DEGs are clustered on chromosomes 5 and 6( Figure 3D , Supplementary Figure 2 ). A PPI network of the glycolysis-related DEGs was constructed based on the MCC algorithm to investigate potential interactions ( Figure 3E , Supplementary Figure 3 ). Node colors represent MCC scores, highlighting hub genes within the network.

Glycolysis-related key DEGs of AP were identified by utilizing machine learning approaches. LASSO regression selected key features with optimal lambda values determined by cross-validation ( Figures 4A, B ). SVM-RFE models were evaluated through accuracy and error rate analyses across multiple parameters ( Figures 4C, D ). Boruta feature selection further refined candidates by assessing feature importance and Z-scores over multiple iterations ( Figures 4E, F ). Integrating results from these methods identified 2 overlapping key glycolysis-related DEGs: Claudin-3 (CLDN3) and CD44 molecule ( Figure 4G ). Validation in an independent dataset confirmed significant upregulation of these genes in the AP samples ( Figure 4H ). In the validation datasets, ROC analysis yielded an Area Under the Curve (AUC) of 1.0 for both Cldn3 and Cd44 ( Figure 4I , Supplementary Figure 4 ). Rather than indicating diagnostic applicability, these results primarily serve as evidence that the differential expression of these genes is robust and consistently observed across independent datasets.

GSEA was performed on GO terms to explore the biological roles of CLDN3 and CD44 ( Figure 5 ). The upregulated GO terms of these two genes were related to cell adhesion, extracellular vesicle formation, and RNA processing, which were essential for immune cell communication, migration, and activation. Conversely, the downregulated GO terms were enriched in sensory perception, ion transport, and transcription factor activity, which reflected a suppression of neuronal-like signaling pathways. These results suggest that CLDN3 and CD44 may contribute to immune-mediated structural remodeling and transcriptional regulation within the inflammatory microenvironment of AP.

Single-cell quality control confirmed high data integrity ( Figure 6A ). Unsupervised clustering and marker-gene annotation identified the main cell populations of the pancreas—acinar, ductal, endothelial, fibroblast, B-cell, T-cell, NK-cell, macrophage, and neutrophil clusters ( Figures 6B, C ). There was an increasing number of immune cell populations and a remarkable loss of acinar cells in AP samples ( Figures 6D, E ). Cell-type–resolved mapping revealed the expression changes of glycolysis-related key genes in AP (Figures 6F-I). CD44 was enriched in both ductal and immune cells, which is consistent with the enhanced activation and migration of the immune cells. Notably, CLDN3 was significantly downregulated in acinar cells, indicating a disruption in epithelial integrity. This finding was inconsistent with bulk RNA-seq results, thereby prompting subsequent protein-level validation.

The ssGSEA analysis revealed distinct immune infiltration patterns of the pancreas in AP ( Figure 7A ). An increase in the expression of immune cell subsets was observed, indicating an alteration from a steady-state immune surveillance in normal status to an activated state. Correlation analysis among immune cells demonstrated the coordinated infiltration patterns, which indicated synergistic immune responses to the microenvironment in AP ( Figure 7B ). Notably, glycolysis-related key DEGs showed strong correlations with the majority of immune cells ( Figure 7C , Supplementary Table 5 ), suggesting their broad regulatory roles in shaping the immune landscape during AP.

To elucidate the upstream regulation of the glycolysis-related key DEGs, TFs, and miRNA networks were assembled from curated datasets ( Figure 8 ). A total of 33 TFs were predicted to target these genes. At the post-transcriptional level, 138 miRNAs were identified to regulate the glycolysis-related key DEGs. Focusing on miRNAs with a degree of ≥1, 21 miRNAs were shared between CD44 and CLDN3( Supplementary Figure 6 ). Notably, SAP30, KDM5B, and PHF8 were identified as TFs regulating CD44, while TRIM24, SMARCA5, E2F5, and BCOR were predicted to regulate CLDN3. The identification of these shared miRNAs suggests post-transcriptional regulation of CD44 and CLDN3.

To validate the bioinformatic findings, we established a mouse AP model by the combination of caerulein and LPS, as confirmed by serum amylase, lipase, and trypsin activity, and MPO activity in pancreatic tissue ( Figures 9A-D ), as well as H&E histopathological assessment ( Figures 9F, G ). The elevation of lung MPO activity and H&E scores indicated the involvement of the lungs ( Figure 9E, H, I ). We next assessed Cd44 and Cldn3 expression at both mRNA and protein levels by qPCR and western blotting ( Figures 9J, K ). Cd44 was upregulated at both transcriptional and protein levels in AP mice. Correspondingly, IHC staining showed that CD44 was predominantly expressed in inflammatory cells within pancreatic tissues of AP mice ( Figure 9L ). Interestingly, bulk qPCR indicated increased Cldn3 mRNA, whereas Western blotting revealed a marked protein reduction. This discrepancy likely reflects cell-type compositional shifts in AP, whereby increased Cldn3 expression in infiltrating immune and ductal cells masks its downregulation in acinar cells(p<0.0001), as revealed by single-cell transcriptomics ( Figure 6G ). IHC further showed a redistribution of CLDN3 from acinar membranes in controls to diffuse cytoplasmic staining in AP mice( Figure 9M ).