The datasets GSE194331 and GSE279876 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, accessed on 15 October 2024). GSE194331 is an RNA-seq dataset of peripheral blood from 87 clinical AP patients (Mild = 57, Moderately-Severe = 20, Severe = 10). We extracted peripheral blood gene expression data from 32 healthy individuals and 10 SAP patients for analysis. GSE279876 is a scRNA-seq dataset of pancreatic tissue from mice with AP, comparing normal diet and high-fat diet conditions. AP was induced by intraperitoneal injection of caerulein (50 μg/kg), every hour for 12 consecutive times. Then, pancreatic tissues from each group were collected and subjected to 10x single-cell sequencing. We extracted the data of the normal control group and the acute pancreatitis group induced by a normal diet for analysis.

The GeneCards database (https://www.genecards.org/, accessed on 20 October 2024) was searched for “mitophagy” to obtain MRGs, and genes with relevance score >2 were selected. We search “mitophagy” in PubMed to obtain MRGs from the literature (30). After removing duplicates, 218 MRGs were finally obtained ( Supplementary Material 1 ).

The raw count data were normalized and then analyzed for differential expression using the “DESeq2” package in the R software. Genes with |logFC| > 1 and adjusted p < 0.05 were considered DEGs, and volcano and heat maps were plotted using the ggplot2 and “pheatmap” packages.

Gene symbol conversion was performed using the “org.Hs.eg.db” and “org.Mm.eg.db” packages, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “clusterProfiler” package, and the results were visualized using the “ggplot2” package.

The “WGCNA” package was used to construct the gene co-expression matrix. Samples were first clustered and the optimal soft threshold (β) was determined to be 11. Using parameters minModuleSize = 30 and MergeCutHeight = 0.25, we constructed a scale-free co-expression network, thereby converting the adjacency matrix into a topological overlap matrix (TOM). We performed cluster analysis to identify gene modules, and constructed a dendrogram via hierarchical clustering to calculate the correlation between module eigengenes and disease phenotypes.

To further screen for hub genes, three machine learning algorithms were used: LASSO logistic regression, RF, and SVM-RFE. Previous studies have demonstrated the effectiveness of these algorithms for gene screening (31). The LASSO regression model was built using the “glmnet” package, with the minimum λ value was selected as the optimal parameter for prediction. We constructed multiple decision trees using the “randomForest” package and aggregated their results to perform classification, regression, and feature selection. For the SVM-RFE algorithm, we used the “e1071” and “caret” packages, applying five-fold cross-validation to obtain the results. Finally, venn diagrams were plotted to show the intersection of the three algorithms.

The “ggpubr” package was used to visualize the expression level of MAPK14. Receiver operating characteristic (ROC) analysis was performed using “pROC” package to evaluate its diagnostic ability in SAP.

The GO, KEGG, REACTOME, and HALLMARK gene sets were downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb, accessed on 30 October 2024). Single-gene GSEA was performed using the “clusterProfiler” package to explore the potential functions of MAPK14, with results visualized via the “enrichplot” package.

Single sample gene set enrichment analysis (ssGSEA) was performed using the “GSVA” package. Immune cell infiltration and immune function activity for each sample were calculated based on sample gene expression. Correlations between immune cells abundance and immune pathways activity were calculated using the “corrplot” package. The “ggplot2” package was used to draw boxplots to compare the differences in immune cells infiltration and immune pathways activity between the SAP and control groups. A subset of MAPK14 expression was extracted from the original expression matrix, and the immune cell abundance matrix was combined with the target gene expression matrix. A lollipop plot was created using the “ggplot2” package to demonstrate correlation between MAPK14 expression and immune cells abundance as well as immune pathways activity.

Quality control and preprocessing of data were performed using the “Seurat” package. The “PercentageEigenSet” function was used to remove mitochondrial genes (“^mt-”) and red blood cell genes (Hba1, Hba2, Hbb, Hbd, Hbe1, Hbg1, Hbg2, Hbm, Hbq1, and Hbz). Final cell quality control was carried out based on the following parameters: nFeature_RNA > 300 & nFeature_RNA < 7000 & nCount_RNA > 1000 & percent.mt < 20 & percent.HB < 1. The “NormalizeData” function was used for data normalization, the “FindVariableFeatures” function was used to find 3000 highly variable genes, and the “ScaleData” function was used to standardize the data. After PCA dimensionality reduction, the “Harmony” package was used for data integration. Then, the “FindNeighbors” and “FindClusters” functions were employed for clustering, and the “DoubletFinder” package was used to remove double cells. The results were visualized through t-distributed stochastic neighbor embedding (tSNE). By referring to the CellMarker2.0 (http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.html, accessed on 10 November 2024), PanglaoDB (https://panglaodb.se/index.html, accessed on 10 November 2024), and Cell Taxonomy (https://ngdc.cncb.ac.cn/celltaxonomy/, accessed on 10 November 2024) databases for cell annotation, 20643 cells were annotated as 11 cell types, including acinar cells, fibroblasts, neutrophils, macrophages, T cells, B cells, endothelial cells, mesothelial cells, beta cells, alpha cells, and pericytes.

The cluster-specific genes of macrophages were identified using the “FindAllMarkers” function. Subsequently, adjusted p-values were calculated using the Wilcoxon rank sum test. These differential genes were then used for GO and KEGG enrichment analysis.

C57BL/6 mice were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. Newly arrived mice were allowed free access to food and water for one week. Then they were fasted 12 h before the experiment. To induce SAP, mice were intraperitoneally injected with cerulein (MCE, China) at a dose of 100 μg/kg every 1 h for a total of 8 times. The last injection of cerulein was accompanied by 15 mg/kg of Lipopolysaccharide (MCE, China). For the SB203580 + SAP group, mice were intraperitoneally injected with SB203580 (10 mg/kg) 30 min before the first injection of cerulein. After the treatment, the mice were sacrificed by cervical dislocation, and pancreatic tissues were collected for subsequent experiments. This study was approved by the Ethics Committee of Xuanwu Hospital of Capital Medical University (XW20211223-1).

The tissue sections were deparaffinized in xylene and then hydrated through a series of alcohols with decreasing concentrations to facilitate dye penetration. They were stained with hematoxylin, rinsed under running water to remove excess dye, and then differentiated with 0.1% hydrochloric acid in ethanol. After another rinse to remove the acid, the sections were stained with eosin, followed by a further rinse under running water to remove excess eosin. Subsequently, the sections were dehydrated using a series of alcohols with increasing concentrations, mounted with a neutral resin, and photographed for observation. Two pathologists independently scored the severity of pancreatitis in mice according to Kusske et al.’s pancreatic pathology scoring criteria (32).

The tissue sections were placed in a thermostat at 60°C and baked for 2h. they were dewaxed in xylene and hydrated in gradient ethanol. Subsequently, antigen retrieval was performed using sodium citrate. The sections were then blocked with 5% goat serum. After blocking, they were incubated overnight at 4°C with primary antibodies against p38 MAPK (ABclonal, A14401) and phospho-p38 MAPK (ABclonal, AP0057). After the incubation with primary antibodies, the sections were incubated with secondary antibodies. Subsequently, the expression of target proteins was detected using 3,3′-diaminobenzidine (DAB) solution. Two professional pathologists independently scored the results using the immunoreactivity score (IRS) as follows: The intensity of cellular staining was divided into 4 levels, no positive staining (negative) scored 0, yellow (weakly positive) scored 1, tawny (positive) scored 2, and brown (strongly positive) scored 3. The percentage of positive cells was divided into 4 levels, ≤25% scored 1, 26%-50% scored 2, 51%-75% scored 3, >75% scored 4. The final score was obtained by multiplying these two scores.

The pancreatic tissue was lysed on ice using RIPA lysis buffer and the supernatant was collected by centrifugation. Proteins were separated via sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to polyvinylidene fluoride (PVDF) membranes. After blocking the PVDF membranes with 5% skim milk, the membranes were incubated overnight at 4°C with primary antibodies targeting β-actin (Proteintech, 20536-1-AP), Pink1 (Proteintech, 23274-1-AP), Parkin (Proteintech, 14060-1-AP), and Bnip3l/Nix (Selleck, F0469). Then, the membranes were incubated with secondary antibodies at room temperature for 1 h. Finally, protein bands were visualized using a western blot imaging system.

Bioinformatics-related statistical analyses were performed using R software (R-4.4.1, 64-bit). For molecular biology experiments, data from at least three independent experimental replicates were analyzed using two-sample t-tests or Wilcoxon rank-sum tests using Prism 8 (GraphPad), A p-value < 0.05 was considered statistically significant.

In this study, we used the GSE194331 dataset, which included sample data from 10 SAP patients and 32 healthy individuals. DEGs were screened based on the criteria of |logFC|> 1 and an adjusted p-value < 0.05. A total of 4110 DEGs were identified, comprising 2290 up-regulated genes and 1820 down-regulated genes ( Figure 1A ). The top 50 up-regulated and top 50 down-regulated genes were respectively selected for heatmap visualization ( Figure 1B ).

To further explore the potential biological functions regulated by these DEGs, we conducted GO and KEGG enrichment analyses on the screened DEGs. The results indicated that these DEGs were primarily involved in various biological processes, including the immune response-activating signaling pathway, positive regulation of cytokine production, mononuclear cell differentiation, leukocyte mediated immunity, and immune response-regulating cell surface receptor signaling pathway. Moreover, the DEGs were also involved in multiple signaling pathways, such as cytokine-cytokine receptor interaction, Neutrophil extracellular trap formation, and NF-kappa B signaling pathway ( Figures 1C, D ).

To identify core genes associated with SAP, we constructed a gene co-expression network using WGCNA. Sample clustering analysis revealed distinct clustering of samples ( Figure 2A ), indicating reliable data quality. A soft threshold of β = 11 was selected, yielding a scale-free topology fit index (R²) of 0.85 ( Figure 2B ), which confirmed the construction of a scale-free co-expression network. Using gene correlation matrices, we constructed a hierarchical clustering dendrogram of genes, identifying 15 distinct gene modules ( Figure 2C ). Among these, the “blue” module, containing 924 genes, exhibited the strongest correlation with SAP phenotypes (Cor = 0.79, p = 5 × 10^-10), making it the most clinically relevant module ( Figure 2D ). A scatter plot revealed a significant positive correlation between gene salience (GS) and module members (MM) within the “blue” module (Cor = 0.84, p < 1 × 10^-200), indicating high consistency between module genes and SAP relevance ( Figure 2E ).

To investigate the role of mitophagy in SAP pathogenesis, we performed an intersection analysis among DEGs, MRGs, and genes in the “blue” module, identifying 8 overlapping signature genes ( Figure 3A ). These genes were subjected to GO and KEGG enrichment analyses to characterize their potential biological functions and pathways associated with mitophagy in SAP pathogenesis. GO functional enrichment revealed that the signature genes mainly regulate biological functions such as autophagy of mitochondrion, organelle disassembly, macroautophagy, autophagosome organization, and vacuole organization ( Figure 3B ). KEGG pathway enrichment revealed their involvement in autophagy, NOD-like receptor signaling pathway, mitophagy, FoxO signaling pathway, efferocytosis, sphingolipid metabolism, and cholesterol metabolism ( Figure 3C ).

We utilized three machine learning algorithms to further screen for core genes. Using LASSO regression, we identified five core genes that best characterize the MRGs in SAP: ATG3, MAPK14, CAMKK2, TSPO, and GABARAPL2 ( Figures 4A, B ). RF identified two significant genes: MAPK14 and MFF ( Figures 4C, D ), SVM-RFE identified three candidate genes with minimal error and maximal accuracy: MAPK14, MFF, and ATG3 ( Figures 4E, F ). By intersecting the results from the three algorithms, we identified MAPK14 as the mitophagy-related gene most strongly associated with SAP ( Figure 5A ). Further analysis revealed that MAPK14 expression was significantly higher in SAP patients than in healthy individuals ( Figure 5B ). Moreover, ROC curve analysis showed that MAPK14 had a good predictive value for SAP, with an AUC value of 0.900 ( Figure 5C ).

Our previous analysis identified MAPK14 as the mitophagy-related gene most strongly associated with SAP pathogenesis, potentially regulatory SAP progression. To further explore its molecular functions and regulated biological pathways, we performed single-gene GSEA. MAPK14 was enriched in GO terms including ATP synthesis coupled electron transport, macroautophagy, mitochondrial RNA metabolic process, oxidative phosphorylation, and positive regulation of inflammatory response ( Figure 5D ). GSEA of the KEGG gene set revealed that MAPK14 is involved in regulating pathways such as electron transfer in complex IV, IL1-IL1R-p38 signaling pathway, mismatch repair, PRNP-PI3K-NOX2 signaling pathway, and TGFR-EGFR-RAS-ERK signaling pathway ( Figure 5E ). GSEA of the REACTOME gene set revealed that MAPK14 is involved in regulating pathways such as antigen processing cross presentation, apoptosis, cellular response to hypoxia, DNA repair, and eukaryotic translation initiation ( Figure 5F ). And GSEA of the HALLMARK gene set revealed that MAPK14 is involved in regulating pathways such as apoptosis, hypoxia, complement, IL6-JAK-STAT3 signaling, and inflammatory response ( Figure 5G ). This further highlights the potential of targeting MAPK14 to modulate SAP acinar cell death and inflammation.

Immune mechanisms play an important role in the development of SAP. Therefore, we explored the differences in immune infiltration between SAP patients and normal individuals using ssGSEA. The heatmap shows the distribution of 28 immune cells in 32 healthy individuals and 10 SAP patients in the GSE194331 dataset ( Figure 6A ). Figure 6B shows the correlation between these 28 immune cells. We observed that the infiltration of neutrophil, monocyte, activated dendritic cell, immature dendritic cell, regulatory T cell, macrophage, gamma delta T cell, and mast cell was significantly higher in SAP, while the infiltration of immature B cell, central memory CD4 T cell, activated CD8 T cell, effector memory CD8 T cell, and activated B cell decreased in SAP. This suggesting that these immune cells may play a crucial role in the development of SAP ( Figure 6C ). In addition, we further explored the distribution of 11 immune functions between SAP patients and healthy individuals, as well as the correlation between these 11 immune functions ( Figures 7A, B ). Differential analysis of immune function revealed that T cell co-stimulation, inflammation-promoting, APC co inhibition, parainflammation, and check-point were closely associated with the development of SAP ( Figure 7C ). Correlation analysis revealed that MAPK14 was statistically correlated with 23 immune cells and 8 immune functions. Among them, regulatory T cell, gamma delta T cell, macrophage, neutrophil, mast cell, activated dendritic cell, type 17 T helper cell, immature dendritic cell, and APC co inhibition immune function were positively correlated with MAPK14 ( Figures 7D, E ).

To investigate the potential role of MAPK14 in AP at the single-cell level, we analyzed the scRNA-seq data from the NFD group (control) and the AP group in the GSE279876 dataset. A total of 20643 cells from C57BL/6 mouse pancreatic samples were processed, and the dataset underwent dimensionality reduction, clustering, and tSNE visualization ( Figure 8A ). Fifteen major cell clusters were identified, and the “FindAllMarker” package was used to identify marker genes ( Figure 8B ). Based on these markers, 20643 cells were annotated into 11 cell types: acinar cells, fibroblasts, neutrophils, macrophages, T cells, B cells, endothelial cells, mesothelial cells, beta cells, alpha cells, and pericytes ( Figure 8C ). Figure 8D displays the marker genes for the 15 cell clusters. We found that Mapk14 was expressed across multiple cell types, with the highest expression observed in macrophages ( Figures 8E, F ). To explore its role in macrophages, we extracted macrophage subpopulations for differential gene expression analysis and performed GO and KEGG enrichment analysis. Differentially expressed genes were primarily involved in biological processes such as RNA splicing, immune response activating signaling pathway, and regulation of apoptotic signaling pathway ( Figure 8G ), as well as pathways such as endocytosis, MAPK signaling pathway, autophagy, oxidative phosphorylation, lysosome, and apoptosis ( Figure 8H ). These findings indicate that Mapk14 may play a regulatory role in these biological processes, highlighting its potential importance in macrophage-mediated functions during AP.

To investigate the role of p38α in SAP pathogenesis in vivo, we established SAP mouse models. IHC staining of mouse pancreatic tissues revealed that p38α expression was significantly higher in SAP mice than in control mice, with phosphorylated p38α exhibiting a more pronounced increase ( Figure 9A ). This indicates activation of the p38α in the pancreas of SAP mice. To validate the role of p38α in SAP-associated inflammation, we treated mice with the p38-specific inhibitor SB203580 and observed significant reductions in pancreatic inflammation and lung injury ( Figure 9B ). In addition, we found that the expression of mitophagy markers Pink1, Parkin, and Bnip3l/Nix was reduced in the pancreatic tissue of SAP mice, indicating impaired mitophagy in SAP ( Figures 9C, D ). The expression of these markers was restored after p38α inhibition ( Figures 9C, D ). This suggests that inhibiting p38α may alleviate pancreatitis severity by promoting mitophagy in SAP.