IL18 Works Like a Two-Side Coin in Acute Pancreatitis

Introduction

Acute Pancreatitis (AP) is a common acute non-infectious inflammatory disease of the pancreas in clinic. The diagnostic criteria for AP need to consider comprehensively by combining serum amylase, lipase, imaging tests, and clinical symptoms. The causes of AP mainly include biliary calculi, hyperlipidemia, drug-induced, procedural, and non-specific, etc.¹ The incidence of AP worldwide is 34 cases per 100,000 person-years, and it is on the rise globally,² causing a huge burden on medical resources and social economy.³ In addition, about 18% of patients with AP will have a recurrence, and 8% of patients will develop chronic pancreatitis, even periampullary cancers.⁴ Both of these conditions will bring additional economic burdens to the healthcare system and long-term troubles to patients’ lives and work.^5,6

In AP, multiple mechanisms (such as pancreatic duct obstruction, direct alcohol toxicity) disrupt the calcium homeostasis within the acinar cells, leading to premature activation of trypsinogen, destruction of the function of the secretory cells, release of inflammatory and cytokines, triggering a cascading molecular interaction network, and digestion of the pancreas and surrounding tissues.⁷ The main cellular events include calcium ion overload, mitochondrial dysfunction, endoplasmic reticulum stress (ERS), impaired unfolded protein response (UPR), and impaired autophagy, etc.^8,9

Acute pancreatitis (AP) is a clinically common inflammatory disease that is mild in most situations.¹⁰ About 35% of AP patients will progress to moderate severe acute pancreatitis (MSAP) or severe acute pancreatitis (SAP), conditions with poor prognosis and life-threatening.¹¹ AP is mostly self-limiting in the early stage of the disease, thus early identification contributes to timely treatment with good prognosis.¹² With the development of genomics in the last decade, our understanding of the pathophysiological mechanism of AP has gradually shifted from traditional clinical indicators to the gene era.^13–15

High-throughput profiling methods have accelerated basic research and made deep molecular characterization of patient samples routine.¹⁶ Machine learning is superior to traditional statistical models and can build informative and predictive models of the underlying biological process.^17–19 It is becoming an integral part of modern data mining and clinical diagnosis.^20–22 Eolving diagnostic tools (including novel biomarkers and imaging modalities) and innovative therapies (from precision medicines to immunotherapies) are improving survival and quality of life.²³ And it is the computational advancement that brings in a new era of personal health-based therapy.²⁴

Interleukin-18 (IL-18) is a cytokine that shares structural features with the interleukin-1 (IL-1) family.²⁵ As a unique cytokine, it is involved in the activation and differentiation of multiple T cell groups.²⁶ IL-18 expression was increased in AP at an early stage and was proved to be correlated with disease severity.^27–29 At present, the mainstream view is that IL-18 is involved in the deterioration of AP.^27,30 However, some scholars believe that IL-18 seems to have a protective effect on AP.³¹

Mendelian randomization (MR) uses genetic variation “single nucleotide polymorphisms (SNPs)” as an instrumental variable (IV) to determine whether observed associations between risk factors and outcomes are consistent with causal effects.³² Since these genetic variants are not usually associated with confounders, differences in outcomes between those who carry the variant and those who do not can be attributed to differences in risk factors.

In this study, we used high-throughput sequences via machine learning to find the characteristic gene of AP. Then the expression level, clinical correlation, biological function and immune infiltration of IL18 in AP were analyzed comprehensively. In addition, we established murine models to validate the relative mRNA and protein expression of IL18 in experimental AP. Through MR study, we initially explored and speculated that there might be a connection between IL18 and the self-limiting course of AP.

Materials and Methods

Datasets and Mice Models

The series matrix file GSE194331 of AP patients was obtained from GEO database, downloaded from https://www.ncbi.nlm.nih.gov/geo/. The dataset was derived from whole blood samples, including 87 AP patients and 32 healthy controls (total=119). GTEx data is downloaded from UCSC (https://xena.ucsc.edu/). The GEO and GTEx databases are publicly available, which does not require institutional review board approval and informed consent. According to the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects (China, February 18, 2023), Article 32, items 1 and 2, research using publicly available, anonymized data that does not involve personal privacy or re-identification is exempt from IRB approval.

Male C57BL/6J mice weighing 20.76–22.12 g were purchased from SLAC Animal Corporation (Shanghai, China). The mice were raised in an animal laboratory with suitable humidity and temperature, adequate water resources and feed. Ten male C57BL/6J mice were randomly divided into two groups, the wild type group and the caerulein-induced AP group (5 mice per group). The AP induction protocol was 100ug/kg of caerulein injected intraperitoneally every 1 hour, repeated for 10 times. The control group was given the same amount of normal saline intraperitoneal injection for the same time period. Mice were sacrificed under pentobarbital sodium after 3 hours of the last injection of caerulein, and pancreas were collected for subsequent analysis. All animal experiments were performed in accordance with the ARRIVE guidelines and the National Institutes of Health Guide for the Care and Use of Laboratory Animals (8th edition). All procedures performed have been approved by the Experimental Animal Ethics Committee of Shanghai First Hospital Affiliated to Shanghai Jiao Tong University (IACUC: 2023AWS208).

Differentially Expressed Genes and Protein-Protein-Interaction (PPI)

R-package of “edgeR” was used to analyze differentially expressed genes in whole blood between normal controls and AP patients. We took mRNA expression LogFC absolute value >1 and false discovery rate (FDR) <0.05 as the threshold points of differential genes. PPI analysis was performed in the STRING database (https://www.string-db.org/). The connectivity diagram is drawn in R language, and the core gene is defined according to the connectivity ≥5. The threshold of 5 was chosen based on the distribution of interaction counts in the PPI network; the top 10% of genes by degree centrality had ≥5 interactions. This is a commonly used cut-off in protein-protein interaction network analyses to identify hub genes. Cytoscape presented the core genes and their interactions.

Biological Role and Disease Analysis

The R packages of “clusterProfiler”, “enrichplot”, “org.Hs.egdb”, “ggplot2”, “GSEABase” and “DOSE” was applied to analyze the function, pathway and disease of core genes. Gene Ontology (GO) functional enrichment was used to analyze the biological significance of core genes, including Biological Process (BP), Cellular Components (CC) and Molecular Function (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was used to analyze the pathways of core genes. Disease Ontology (DO) enrichment analyzed main diseases of core gene.^33,34 The above results were considered statistically significant with p<0.05 and p<0.05 after adjustment (Benjamini–Hochberg false discovery rate). The visualization of DO is achieved through the R package of “GOplot”.³⁵

Machine Learning

In order to reduce bias, we use three different machine learning algorithms to screen potential characteristic genes. The R packages of “glmnet”, “e1071” and “randomForest” completed the Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression, Support Vector Machine Recursive Feature Elimination (SVM_RFE) and Random Forest (RF) respectively. Finally, the venn diagram shows the overlapping characteristic gene in AP after performing three algorithms.

Expression of the Characteristic Gene in Human Tissues

The expression levels of the characteristic gene in human organs were extracted from GTEx database. Then R packages of “dplyr” and “ggpubr” were used to graph the expression of the characteristic gene in various tissues of human body. Moreover, the mRNA expression of the characteristic gene in AP tissues and normal controls was also mapped.

Clinical Value and ROC Curve

We explored the clinical relevance and diagnostic sensitivity of the characteristic gene. The R-package of “ggpubr” compared the expression level of the characteristic gene under different clinical condition. The R package of “pROC” generated the receiver operating characteristic curve (ROC) and calculated the area under the curve (AUC) value to assess the specificity and sensitivity of the characteristic gene to predict AP.

Function and Pathway Analysis of Co-Expressed Genes

We look for correlated genes that are associated with the characteristic gene in R software. The threshold value of correlated genes was set as |Pearson correlation coefficient| p>0.60 and p<0.001. The R packages of “igraph” and “reshape2” were used to visualize all co-expressed results. R packages of “clusterProfiler”, “org.Hs.egdb”, “DOSE”, “ggplot2”, “GOplot”, “R.Utls” and “pathview” were used for GO and KEGG analyses.

Immune Infiltration and Principal Component Analysis (PCA)

The CIBERSORT algorithm can show the relationship between the expression level of the characteristic gene and 22 types of immune cells. Boxplot exhibited the situation of immune cells in AP and normal controls. PCA maps were used for dimensionality reduction to evaluate the feasibility of distinguishing AP among immune cells. Bar charts visualized the correlation between characteristic gene and immune cells. The ssGSEA algorithm is based on 29 immune gene sets (infiltration scores of 16 immune cells and activity of 13 immune-related pathways) to comprehensively quantify the relative abundance of immune cell types, pathways, functions, and checkpoints in each patient. GSVA R package was used to analyze the differences in immune function between AP patients with low- and high-expression of the characteristic gene. The R packages of “limma”, “preprocessCore”, “GSVA”, “GSEABase”, “reshape2”, “corrplot”, “ggpubr”, and “ggplot2” were used to complete the above tasks.

Two-Sample Mendelian Randomization (MR) Analysis

Two-sample MR was used to investigate the causal relationship between the characteristic gene and the risk of AP, and SNP was defined as IVs. The characteristic gene information was acquired from the Genome-Wide Association Study (GWAS). We performed MR analysis of IL18 and AP. The GWAS ID was ebi-a-GCST90010141 for IL-18 (individuals of European ancestry from the UK Biobank and deCODE cohorts), and finn-b-K11_ACUTPAN for AP. MR Analysis was performed based on the R package of “TwoSampleMR” and the relationship between IL18 and AP was evaluated using inverse variance weighting (IVW). Additional sensitivity analysis was performed by MR-Egger.³⁶

Biological Experiments

Total RNA and protein from murine pancreas were extracted using TRIzol and RIPA according to the instructions. Reverse transcription and SYBR kits were used for RT-qPCR (EnzyArtisan). CT values (2−ΔΔCT) were calculated to determine the relative expression level of IL18 referred to mouse Rplp0. The primer sequences were as follows: IL18, forward: 5’-AACTTTGGCCGACTTCACTGTA-3’, reverse: 5’-TATCAGTCATATCCTCGAACACAGG-3’; Rplp0, forward: 5’-TTATAACCCTGAAGTGCTCGAC-3’, reverse: 5’-CGCTTGTACCCATTGATGATG-3’. Western blotting was performed using IL18 (Proteintech, Cat No. 10663-1-AP) and β-actin (Proteintech, Cat No. 66009-1-Ig,) primary antibodies incubated at 4°C overnight. The secondary antibody of the corresponding species (Boster, Cat. No. BA1054; Cat No. BA1050) was incubated for 60 min the next day. The ECL photograph shows the WB results.

Pancreatic tissue was fixed, embedded and sliced to 4 mm. Sections were dewaxed, hydrated and incubated with citrate antigen retrieval solution (Beyotime Biotechnology, Shanghai, China) for 1 hour. The slides were then incubated overnight with primary antibodies against IL18 (Proteintech, Cat No. 10663-1-AP). Fluorescent secondary antibody (Servicebio, GB25303) combined with DAPI (Servicebio, G1127) was added the next day for tissue immunofluorescence analysis. Half of the slides were incubated for 1 hour with a biotinylated secondary antibody (Servicebio, G1213) for immunohistochemistry (IHC) staining using a microscope (Leica, IL, USA, DFC550).

Statistics

The t-test is used to analyze differences between groups for variables with a normal distribution. Otherwise, the Mann–Whitney U-test is applied. Chi-square tests are used to compare quantities. The Pearson correlation method was used to analyze the correlation between two different genes. All statistical analyses were performed by R software version 4.3.1 (https://www.r-project.org/). p<0.05 was considered statistically significant. We have marked * in the results, where * means p<0.05, ** means p<0.01, and *** means p<0.001.

In order to ensure the stability of the results, the following methods were performed: (i) performance metrics (AUC, accuracy, sensitivity, specificity) are reported in Results, (ii) parameter tuning for the LASSO and random forest models (alpha = 1 for LASSO, ntree = 500 for RF), (iii) we used feature selection within cross-validation loops to prevent data leakage, and (iv) the 10-fold cross-validation strategy with 5 repeats, (v) all three machine learning algorithms were performed 30 repeats.

Results

Differentially Expressed Genes and Function of Core Gene

The transcriptome of whole blood samples from patients with AP and normal controls was analyzed. A total of 1356 differential genes were obtained according to the threshold value (Table S1). The PPI interaction network of 1356 differential genes was constructed via STRING online database (http-string-db.org) to better understand the interactions between these differential genes. And number of interactions ≥5 is considered as the threshold for core genes (Figure S1). The cytoHubba module of Cytoscape software computed and visualized these 100 core genes (Figure 1A).

Figure 1 Core genes and related functions, pathways and diseases in AP. (A) Visualize the core genes using Cytohubba plugin in Cytoscape software. (B) The circular scatter plot of GO enrichment to exhibit the biological process. (C) The circular scatter plot of KEGG enrichment. (D) Circos plot illustrates the relationship between core genes and the top 6 kinds of diseases with DO enrichment. Genes are involved in connection with the DO term through colored connecting lines. The outer ring color represents logFC.

To explore the role of these core genes in the AP process, we focused on function, pathway and disease analysis of them. The results of GO analysis are mainly positive regulation of cytokine production, defense response to bacterium and regulation of inflammatory response. This suggests that a variety of cytokines and inflammatory responses take part in the process of AP (Figure 1B). The pathways of core genes are mainly enriched in Tuberculosis, Inflammatory bowel disease and Staphylococcus aureus infection (Figure 1C). DO analysis showed that the onset of AP may be associated with diseases such as periodontal disease, arteriosclerosis, bacterial infectious disease, etc. The circle graph presented the six diseases corresponding to core genes (Figure 1D).

Machine Learning

Three machine learning methods, namely LASSO regression, SVM-RFE and RF, were used to learn core genes in AP (Table S2). LASSO regression obtained 17 core genes in AP (Figure 2A). The SVM-REF algorithm got 13 core genes in AP (Figure 2B). The RF algorithm gained 8 core genes in AP (Figure 2C). A Venn diagram was drawn to indicate the intersection of gene subsets from three machine learning algorithms. Finally, IL18 was identified as the characteristic gene for pancreatitis (Figure 2D).

Figure 2 Identification of potential biomarkers in AP based on machine learning. (A) The best Lambda value in LASSO algorithm corresponds to 17 genes. (B) Root mean square error (RSME) curve of 100 core genes in SVM-RFE algorithm. The blue dots represent the lowest error rate and correspond to 13 genes. (C) Screening diagnostic genes for AP by random forest algorithm. (D) Venn diagram intersects the characteristic gene in AP after LASSO, SVM-REF and RF algorithms.

Expression Levels and Clinical Value of IL18

Expression level of IL18 in all organs was depicted based on the GTEx database. The red box suggested the expression of IL18 in the pancreas (Figure 3A). According to GSE194331, we compared the expression level of IL18 in whole blood samples in AP and control group (Figure 3B). The ROC curve indicated that IL18 could predict AP with good ability, with an AUC value of 0.917 (95% CI=0.859–0.961) (Figure 3C). The whole blood expression of IL18 in all pathological types of AP was higher than that in normal controls (p<0.001, Figure 3D).

Figure 3 Expression and clinical value of IL18. (A) Expression of IL18 in various organs of normal human body. The red box illustrated the pancreas. (B) Violin plot of IL18 mRNA expression in the GEO dataset (GSE194331). (C) ROC curve to evaluate the diagnostic value of IL18 in GEO cohort. (D) Expression of IL18 in different pathological states of AP. **p<0.01, ***p<0.001. Abbreviation: ns, no significance.

Meanwhile, we conducted AP model in mice, HE staining showed that AP animal models were constructed successfully (Figure 4A). The relative protein and mRNA expression level of IL18 in AP was significantly increased than that in wild type controls based on WB, IHC, PCR and IF (Figure 4B–E).

Figure 4 Experimental results indicating the expression of IL18 in normal wild type and AP murine model. (A) HE staining of pancreas in mice models. (B) WB result of relative protein expression of IL18 in normal controls (n=5) and AP mice (n=5). (C) IHC staining result of expression of IL18 in normal controls and AP mice. (D) PCR result of relative expression of IL18 in normal controls (n=5) and AP mice (n=5). (E) IF result of expression of IL18 in normal controls and AP mice (IL18 in red). **p<0.01.

Correlation and Function of IL18

Genes with a correlation coefficient over 0.6 are considered as IL18 related genes. A total of 964 associated genes are listed in Table S3. The corNetwork showed the correlation among the 12 genes with the largest correlation coefficient with IL18 (Figure 5A). GO analysis suggested that BP mainly involves activation of immune response, immune response-regulating signaling pathway and immune response-activating signaling pathway. CC is mainly secretory granule membrane, secretory granule lumen and cytoplasmic vesicle lumen. MF is mainly immune receptor activity, phosphortyrosine residue binding and oxidoreductase activity, acting on a sulfur group of donors (Figure 5C, details in Table S4). KEGG showed that IL18 was mainly associated with Salmonella infection, Osteoclast differentiation, NF-kappa B signaling pathway, Th1 and Th2 cell differentiation and T cell receptor signaling pathway (Figure 5E, details in Table S4). All these indicated that IL18 was closely related to immune response and status in AP. The pathway map exhibited the correlated genes with altered expression in the pathways of NF-kappa B, Th1 and Th2 cell differentiation, and T cell receptor (Figure 5B, D and F).

Figure 5 Function and pathway of co-expressed genes of IL18. (A) The corNetwork diagram of IL18. (B) Expression of IL18 and co-expressed genes in NF−kappa B signaling pathway. (C) Bubble map of GO enrichment analysis for IL18 co-expression genes. (D) Expression of IL18 and co-expressed genes in Th1 and Th2 cell differentiation signaling pathway. (E) Bubble map of KEGG enrichment analysis for IL18 co-expression genes. (F) Expression of IL18 and co-expressed genes in T cell receptor signaling pathway. Red represents positive regulation, while green represents negative regulation.

Immune Infiltration

We used the CIBERSORT method to measure the status of 22 types of immune cells in normal controls and AP (Figure 6A). Differences existed with regard to 12 kinds of immune cells such as B cells naive, Macrophages, T cells CD8 and T cells CD4 memory resting. PCA results showed that AP and controls could be well distinguished based on immune cells (Figure 6B). IL18 was positively associated with 6 types of immune cells and negatively correlated to 4 types of immune cells (Figure 6C). Except Macrophages, Neutrophils and Treg cells, the other immune pathways were downregulated when the expression of IL18 was high, which were assessed via ssGSEA. This might because the data of the AP group covered all stages of mild to severe conditions, and the characteristics of macrophages and neutrophils undergo dynamic changes at different disease progression stages. Overall, most of immune cells and signaling pathways were inhibited under the condition of high expression of IL18 (Figure 6D). Hence we propose that IL18 works like a coin which has two sides, affecting immune cells and immune pathways in different ways.

Figure 6 Analyses of the association between immune cells and IL18 in normal controls and AP. (A) 22 kinds of immune cells status in normal controls and AP (Blue represents controls, while red represents AP samples). (B) PCA of 22 types of immune cells, revealing differences of immune-phenotype in normal controls and AP. (C) Lollipop diagram of correlation between IL18 and 22 kinds of immune cells. (D) Comparison of ssGSEA scores between low- and high-expression levels of IL18. *p<0.05, **p<0.01, and ***p<0.001.

MR Result

Last but not least, we estimated the causal relationship between IL-18 level and AP. Scatter plots showed the causal effect of SNP on AP: the higher IL18 level was, the lower risk for AP became (Figure 7A). The forest map reflects the results of every SNP via wald ratio method, while the bottom red line reflects the IVW result (Figure 7B). We found that IL18 level was associated with the risk of AP by IVW approach, with an OR of 0.908 (95% CI = 0.843–0.978, p = 0.011) (Five MR methods employed along with sensitivity analyses to comprehensively evaluate the causal relationship between IL18 and AP were displayed in Table S6). The IVW heterogeneity test indicated that heterogeneity did not exist (p=0.777). The causal effects of the funnel plot are roughly symmetrical (Figure 7C). Horizontal pleiotropy showed no confounding factors (p=0.514), suggesting that IL18level was reliable to predict the risk of AP. Leave-one-out sensitivity analysis showed that there was no dominant SNP of IL18 level in AP (Figure 7D), demonstrating significant causal association between all the calculated results of SNPs.

Figure 7 Mendelian randomization study on IL18 and AP. (A) Scatter plot showing the causal effect of IL18 on the risk of AP. (B) Forest plot showing the causal effect of each SNP on the risk of AP. (C) Funnel plots to visualize overall heterogeneity of MR estimates for the effect of IL18 on AP. (D) Leave-one-out plot to visualize causal effect of IL18 on AP when leaving one SNP out.

Discussion

AP is an inflammatory disease caused by the activation of intracellular trypsinogen due to various reasons.^37,38 In 2022, Maryam Nesvaderani provided the first high-throughput sequencing data of blood samples in AP patients.³⁹ This raw dataset is remarkable for researchers to probe for machine learning of novel markers in AP.^40,41 With the help of this dataset, we integrated bioinformatics analysis and machine learning to explore meaningful biomarkers in AP. And IL18 was identified to be the characteristic gene of AP. The expression of IL18 increased significantly in AP compared to normal controls, and murine model confirmed the result. ROC curve illustrated that IL18 was a good indicator for clinical diagnosis of AP. Intriguingly, MR analysis suggested that high serum levels of IL18 could be causally associated with a reduced risk of AP.

Although the triggers for AP may vary, the immune response to cell damage is similar.⁴² IL18 was a member of the cytokine IL-1 family and was recognized as an important regulator of inflammation, immune response, and tissue damage.^32,43 On one hand, IL-18 has been shown to induce the production of cytokines and chemokines by neutrophils.⁴⁴ Neutrophil infiltration in the SAP-damaged pancreas was significantly reduced when IL-18 was knocked out, along with reduced T cell activation.⁴⁵ On the other hand, IL18 is involved in Th1 and Th2 immune responses and activation of M2 macrophages.^26,46 Ueno et al reported that compared with AP wild-type mice, serum amylase, lipase and the numbers of acinar cells with parenchyma vacuolization were significantly increased in AP mice with IL18 knockout, and the above parameters could be improved after pretreatment with recombinant mouse IL-18. Therefore, they believe that IL18 has a protective effect in AP.³¹ Based on the results in our study, IL18 serves as an indicator for predicting and reflecting the course of AP; and from the perspective of genetic predisposition, the high expression level of IL18 is a protective factor for the occurrence of AP in human. Therefore, the bidirectionality of the research results obtained by previous studies is reasonable. Because the intervention of genes in the modeling process is equivalent to a genetic change. At the same time, we also speculate that after IL18 increases, there may be a similar protective mechanism similar to genetic predisposition.

In summary, the level of IL-18 (p<0.05) was correlated positively with the psoriasis area and severity index at an early immunosuppressive state.⁴⁷ Since functions of cell immunity and humoral immunity in primary nephrotic syndrome (PNS) patients remain untreated were suppressed and disordered, researchers found that IL-18 increased, while IL-10 and Foxp3+Treg cells decreased.⁴⁸ Similarly, in metastatic colorectal cancer, PBMCs exhibited significantly lower production of IL-18 cytokines and lower auto-MLR responses, whereas Treg frequency, IL-10 cytokines were increased compared to healthy donors.⁴⁹

Although our study find IL18 is vital in AP, and performed animal experiments to verify our results, there are still some limitations. In the future, we plan to recruit more real-world human samples in our clinical work and use detection methods such as ELISA to further validate the diagnostic efficacy of IL18 in AP. Additionally, we will collect more detailed information about AP patient, aiming to analyze whether IL18 is associated with the scores in the assessment of AP severity. Furthermore, larger-scale animal study including male and female models with more inflammatory markers (IL-1β, IL-6, TNF-α, and pancreatic histology scoring) needs to be done to further elucidate the mechanism of IL18 in AP.

Conclusions

Our study suggests that IL18 promoted several immune cells at the same time suppressed some immune cells and pathways. But there are too many confounding factors. Therefore, we also conducted a two-sample MR Analysis of GWAS data to explore the association between IL18 and the risk of AP.⁵⁰ The results showed that there was a causal relationship between increased IL18 level and reduced risk of AP. In conclusion, IL-18 is like a coin with two sides: genetic speaking, IL18 can reduce the risk of AP; while during the course of AP, IL18 can serve as a biomarker. We hypothesize that the elevated level of IL18 induced by AP inflammation might also have a protective effect similar to genetics at a certain stage, that’s perhaps why we could see AP is a self-limiting disease in its early stage in clinic.