- 1Department of Pathology, Affiliated Hospital, Beihua University, Jilin, China
- 2Department of Pathophysiology, Basic medical college, Beihua University, Jilin, China
- 3Department of Cell Biology, Basic medical college, Beihua University, Jilin, China
- 4Department of Biochemistry and Molecular Biology, Basic medical college, Beihua University, Jilin, China
- 5Department of Physiology, Basic medical college, Beihua University, Jilin, China
- 6Department of Biochemistry Laboratory, Basic medical college, Beihua University, Jilin, China
Introduction: Gastric cancer (GC) is one of the most frequently encountered malignant tumors in the clinic. Because effective early screening techniques are lacking, most patients have advanced disease at first diagnosis. The interleukin (IL)-36 family plays a vital role in regulating the immune system, inflammatory responses, and the occurrence and development of cancer. Hence, this study explored the potential role of IL-36 related genes (IL-36RGs) in GC and built a prognostic risk assessment model for GC based on IL-36RGs, which can help evaluate treatment and prognosis.
Methods: First, relevant datasets were downloaded from public databases. After processing the datasets to remove batch effects, perform differential analysis, and take intersections, IL-36-related differentially expressed genes (IL-36RDEGs) were screened. A prognostic risk model containing nine model genes was constructed based on univariate Cox and least absolute shrinkage and selection operator (LASSO) regression methods. Then, to investigate the potential biological activities of the model genes in GC, we conducted enrichment, PPI interaction network, and immune infiltration analyses. Immunohistochemical staining was conducted to validate the expression of IL-36A in GC.
Results: The prognostic risk model analysis revealed that mortality events in the high-risk group were substantially elevated compared to those in the low-risk group. The model demonstrated excellent predictive capability at 1, 2, and 3 years and showed the best clinical predictive performance at 3 years. Bioinformatics analysis of the model genes indicate that they primarily participate in mechanisms that promote the synthesis and secretion of cytokines in GC. And hub genes may be strongly correlated with host immune response mechanisms. According to the immunohistochemical staining results, IL-36A expression was higher in the STAD group than in the control group.
Conclusions: The results of the above analysis suggest that IL-36RDEGs can serve as independent prognostic biomarkers for GC and provide insights into IL-36RGs from both bioinformatics and experimental validation perspectives.
1 Introduction
Gastric cancer (GC) is one of the most prevalent malignant tumors encountered in clinical practice. Statistical data show that GC ranks fifth as the most frequently occurring cancer and fourth leading cause of death, especially in East Asia and Eastern Europe, where it has the highest incidence and usually a poor prognosis (1). Owing to the lack of obvious symptoms in early-stage patients and the insidious course of the disease, including the current lack of effective early screening methods in clinical practice, the majority of patients are found to be in an advanced stage upon their initial diagnosis (2). Although the pathogenesis of GC is not completely understood, its growth and differentiation are regulated by various factors. Most clinical treatments include surgical resection, targeted therapy, systemic chemotherapy, and radiotherapy. Although these treatment methods are effective (3, 4), the long-term survival outcomes for patients with GC are unsatisfactory, especially for patients with advanced GC (5). Therefore, identification of new prognostic biomarkers and molecular targets is urgently required to better predict the prognosis of patients with GC and guide individualized clinical diagnosis and treatment.
IL-36 belongs to the IL-1 family of cytokines, including three isoforms: IL-36A (IL-36α), IL-36B (IL-36β) and, IL-36G (IL-36γ). It forms a complex by binding to IL-36R (IL-1RL2) and recruiting IL-1 receptor accessory protein (IL-1RAcP), it subsequently triggers the activation of associated signaling pathways, including MyD88-NF-κB (Myeloid differentiation primary response 88) and Mitogen-activated protein kinase (MAPK), which modulate the expression of downstream target genes (6); IL-36Ra and IL-38 are antagonists of IL-36R; upon binding to IL-36R, they prevent the recruitment of IL-1RAcP, thereby inhibiting the assembly of an active IL-36R complex and limiting signal transduction (7, 8). Studies have shown that IL-36 and its related signaling pathways are associated with various inflammatory diseases, autoimmune diseases, and cancers (9–12). The study found that, compared to adjacent tissues, the expression of IL-36 is significantly reduced in human hepatocellular carcinoma (HCC) tissues. HCC cases with IL-36 positive expression have a lower recurrence rate and longer survival time (13). IL-36A inhibits HCC proliferation, survival, and migration, which correlates with a decrease in the expression of cytokines IL-1β and IL-18, suggesting that IL-36A may inhibit pyroptosis (14). In contrast to non-cancerous tissues, the mRNA and protein expression of IL-36 family members is increased in colorectal cancer tissues, and IL-36R can activate the oncogenic phenotype of cancer cells (11, 15). The relationship between IL-36RGs and tumors, whether inhibitory or promotional, may be related to different types of cancer, various samples, and the heterogeneity of the cancer itself. Over the last few years, studies have revealed a significant association between chronic inflammation in the tumor microenvironment and cancer, with IL-36 considered to play a major role in aseptic chronic inflammation (8, 16). A studies has shown that IL-36 is a potent activator of innate immune cells and mediates a strong anti-tumor response through complement and adaptive immunity. IL-36-treated neutrophils can directly kill tumor cells, induce NK cells to generate cytolytic activity, and enhance T cell proliferation. The interactions between these IL-36-treated neutrophils and other immune components in the tumor microenvironment (TME) result in a highly effective anti-tumor response (17). Research has confirmed that IL-36β can promote the activation of CD8+ T cells by activating mTORC1 through PI3K/Akt, IKK, and MyD88 pathways, thereby enhancing the anti-tumor immune response and laying the groundwork for the application of IL-36β in tumor immunotherapy (18). However, the relationship between IL-36RGs and GC has rarely been studied, and their mechanism of action remains unclear.
To reveal the relationship between IL-36RGs and GC, this study integrated and analyzed data from multiple public databases, including TCGA-STAD, GSE19826, and GSE54129. By implementing strict data preprocessing techniques (such as batch effect correction and differential analysis), we obtained IL-36RDEGs. Based on this, we constructed a prognostic risk model using various statistical methods, including univariate Cox regression and LASSO regression, and performed external validation with an independent GEO dataset to ensure the robustness and wide applicability of the study results. Through systematic multi-cohort cross-validation, this research can contribute to a comprehensive assessment of the role of IL-36RDEGs in GC prognosis and their potential mechanisms, thereby providing a theoretical basis for the clinical development of new prognostic biomarkers and molecular therapeutic targets.
2 Materials and methods
2.1 Technical roadmap
The technical roadmap is shown (Figure 1).

Figure 1. Work flow diagram of this study. TCGA, The Cancer Genome Atlas; STAD, Stomach adenocarcinoma; DEGs, Differentially Expressed Genes; GSEA, Gene Set Enrichment Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ExpDiff&ROC, Expression Differential & Receiver Operating Characteristic; LASSO, Least Absolute Shrinkage and Selection Operator; DCA, Decision Curve Analysis; GSVA, Gene Set Variation Analysis; ssGSEA, single-sample Gene-Set Enrichment Analysis; PPI, Protein-Protein Interaction; IHC, Immunohistochemistry.
2.2 Data download
The TCGA-STAD dataset was sourced from the TCGA database (https://portal.gdc.cancer.gov/) with the aid of the R package TCGAbiolinks (19) and was used as the test set, excluding samples that did not have clinical details. After normalization, the associated clinical data were extracted from the UCSC Xena database (20) (https://xena.ucsc.edu/) and are summarized in Table 1. The GSE19826 and GSE54129 datasets were obtained from the GEO database (21) (https://www.ncbi.nlm.nih.gov/geo/) using the GEOquery (22) package in R for further validation purposes (23). The sva package in R was used to eliminate batch effects from the data (24). The GSE19826 and GSE54129 datasets were processed for probe annotation, standardization, normalization, and other treatments using the R package limma (25). All selected STAD and control samples were used in this study. Further information is provided in Table 2. A total of 44 IL-36RGs were sourced from the GeneCards (26) (https://www.genecards.org/) and MsigDB (27) (Molecular Signatures Database). After using the term “Interleukin-36” as a search keyword and keeping only “protein-coding” IL-36RGs, a total of 44 IL-36RGs were obtained in the GeneCards. In addition, we also got seven IL-36RGs in MsigDB database. IL-36RGs obtained in the above way were merged and deduplicated to obtain a total of 44 IL-36RGs.
2.3 IL-36-related differentially expressed genes
We employed the limma package in R to perform differential analysis of genes in the STAD and Control groups across the TCGA-STAD, GSE19826, and GSE54129 datasets. We set the thresholds for DEGs as |log2FC| > 0 and p value < 0.05, and the outcomes of the differential analysis were plotted as a volcano plot using the ggplot2 package in R. We intersected the DEGs obtained from the differential analysis in the TCGA-STAD dataset with IL-36RGs and obtain IL-36RDEGs to draw a Venn diagram, which was then displayed as a heatmap using the pheatmap package in R.
2.4 Creation of a risk model for predicting prognosis related to GC
Using the R package survival (28), the prognostic risk model in the TCGA-STAD dataset was created to analyze the effect of IL-36RDEGs on prognosis through single and multivariate Cox regression analyses based on clinical information and to ascertain whether IL-36RDEGs are an independent prognostic factor. The process began with univariate Cox regression analysis of genes showing p value < 1, after which LASSO regression analysis with family = “cox” applying the package glmnet in R (29) was conducted to ascertain the model genes for the prognostic risk model. We adopted a 10-fold cross-validation approach to determine the variables in the LASSO regression model and identify the optimal penalty parameter λ. After completing the training for all 10 folds, we calculated the average of C-index to ascertain the optimal λ that produced the highest mean C-index. This optimal λ was employed to construct the LASSO regression model, facilitating the selection of the most significant genes.
Finally, using the LASSO risk score and clinical information, we conducted a multivariate Cox regression analysis. The following approach was used to compute risk scores:
Furthermore, a risk factor graph was created using the LASSO risk score and the package ggplot2 in R.
This study sought to explore the differences in overall survival (OS) between high-risk (High) and low-risk (Low) patients in the STAD group of the E dataset. Using the R package survival, we performed a Kaplan–Meier (KM) (30) curve analysis and plotted the KM curves according to the LASSO risk score. Then, the survivalROC package in R (31) was employed to generate time-dependent ROC curves utilizing the LASSO risk score and OS data, and the area under the curve (AUC) was computed to estimate the survival outcomes for 1, 2, and 3-year periods within the STAD group from TCGA-STAD dataset.
2.5 Analysis of enrichment in GO and KEGG pathways
To further investigate the biological and signaling pathways associated with the genes linked to the risk score, we performed enrichment analysis using Gene Ontology (GO) (32) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (33). The model genes were subjected to GO and KEGG enrichment analysis using the clusterProfiler (34) package in R, utilizing entry screening criteria that considered a p value < 0.05 and an FDR value (q value) < 0.25 as statistically significant. Finally, the GO and KEGG enrichment analysis results were visualized using the R package PathView (35).
2.6 Gene set enrichment analysis
In the TCGA-STAD dataset, the STAD group was sorted into high- and low-risk groups using the median value of the LASSO risk score, and GSEA was performed for all genes in the STAD group using the clusterProfiler package in R (34). The GSEA screening criteria were set as adj.p < 0.05 and an FDR value (q value) < 0.25, correcting the p value using the Benjamini–Hochberg method.
2.7 Gene set variant analysis
The h.all.v7.4.symbols.gmt gene set was sourced from the MSigDB database, and GSVA (36) was conducted on all genes in the TCGA-STAD dataset to evaluate the differences in enrichment of functions among the high- and low-risk groups, using a GSVA selection criterion of p value < 0.05. This was performed to assess whether the various samples were enriched in different pathways.
2.8 Assessment of the prognostic risk model for GC
To illustrate the findings of the multivariate Cox regression analysis, a forest plot was created to represent the expression of the LASSO risk score and clinical variables considered in the analysis. A nomogram (37) was created based on the outcomes of the multivariate Cox regression analysis using the rms package in R, illustrating the correlation between clinical information and the LASSO risk score inside the multivariate Cox regression model. Calibration analysis was conducted to draw calibration curves and evaluate the accuracy and discrimination capacity of the prognostic risk model based on the LASSO risk score. Applying the R package ggDCA, a decision curve analysis (DCA) (38) diagram derived from the LASSO risk score was generated to assess the accuracy and discrimination of the prognostic risk model for GC.
2.9 PPI interaction network and functional similarity analysis
Using the STRING database (39), hub genes were screened, a hub gene-related protein-protein interaction network (PPI Network) was created, and the PPI network model was visualized using Cytoscape (40). Using the GeneMANIA database (41) (http://genemania.org), we predicted genes with functions analogous to the hub genes online and downloaded the interaction network. The inner circle in the figure represents the hub genes in our study, and the outer circle represents functionally similar genes.
2.10 Validation of differential expression of hub genes and ROC curve analysis
To further explore the differences in hub gene expression in the TCGA-STAD dataset and in STAD and Controls in the GSE19826 and GSE54129 datasets, group comparison plots were constructed based on the expression of hub genes. Finally, the package pROC in R was employed to generate the ROC curve for the hub genes and to compute the AUC to assess the diagnostic impact of hub gene expression levels on the occurrence of GC. Subsequently, to investigate the connections among hub genes, the Spearman algorithm was employed to assess the association between hub gene expression and TCGA-STAD dataset. A heat map generated using the pheatmap package in R represents the outcomes of the relationship analysis.
2.11 Immune infiltration examination of the high- and low-risk groups
First, various infiltrating immune cell subtypes were identified. Subsequently, the enrichment scores derived from ssGSEA (42) represented the comparative immune cell infiltration abundance in each sample, which produced an immune cell infiltration matrix for the STAD group in the TCGA-STAD dataset. The package ggplot2 in R was used to demonstrate the differences in immune cell expression between the high- and low-risk groups within the STAD group of TCGA-STAD dataset. The immune cells showing considerable differences between the two groups were filtered for additional analysis; the correlation between immune cells and hub genes was assessed using the Spearman algorithm; and the findings were depicted through a correlation bubble map produced with the R package ggplot2.
2.12 Patient and tissue samples
This study was approved by the Ethics Committee of the Affiliated Hospital of Beihua University in Jilin City, Jilin Province, China (Approval No. 20240084). All patients signed an informed consent form. This study used GC and adjacent tissue samples collected by our research group with clear pathological diagnoses. All patients with GC underwent radical surgery and did not receive endocrine or radiation therapy prior to surgery.
2.13 Immunohistochemical staining and evaluation methods
Immunohistochemistry was conducted using an IL-36A antibody (24173-1-AP; Proteintech, China) diluted at 1:100, following the manufacturer’s instructions and based on preliminary experiments. Two pathologists independently assessed the pathological grouping of each slice. The positive expression rate of the tissue slices was calculated using the ImageJ software. The IL-36A positive expression rate was calculated as follows: Positive area/Total area × 100%.
2.14 Statistical analysis
The R software was used to analyze all data in this study (Version 4.4.1). Unless otherwise specified, statistical significance for comparisons of two categories of continuous variables was examined by applying the independent Student’s t-test for normally distributed variables and the Mann–Whitney U test, commonly known as the Wilcoxon rank-sum test, for variables that were not normally distributed. The Kruskal–Wallis test was used to assess differences across three or more groups. Correlation coefficients between different molecules were calculated using Spearman’s correlation analysis. All statistical p-values were calculated using a two-tailed method, with a p value < 0.05, considered statistically meaningful. Quantitative data are represented as mean ± standard deviation (SD).
3 Results
3.1 Processing of GC datasets
First, batch effects in the GSE19826 and GSE54129 datasets were eliminated using the R package sva. Subsequently, boxplots were generated to compare the differences in illustration values of the datasets before and after batch-effect removal (Figures 2A–D). The outcomes from the box plots indicate that the batch effects among the samples in the datasets were largely eradicated after batch removal.

Figure 2. Batch Effects Removal of GSE19826 and GSE54129. (A) Distribution boxplot of GSE19826 dataset before going batch. (B) Distribution boxplot of the post-batch GSE19826 dataset. (C) The distribution of GSE54129 dataset boxplot before batch processing. (D) The distribution boxplot of the post-batch GSE54129 dataset.
3.2 Differentially expressed genes related to IL-36
Differential analysis was conducted on STAD and control samples from the TCGA-STAD dataset, together with the GSE19826 and GSE54129 datasets, applying the R package limma to identify DEGs, and volcano plots were created from the results of each dataset’s differential analysis (Figures 3A–C). The intersections of all DEGs from the TCGA-STAD dataset, which fulfilled the thresholds of |log2FC| > 0 and p < 0.05, and IL-36RGs were selected. A Venn diagram was plotted showing 16 IL-36RDEGs (Figure 3D). Following intersection analysis, the differences of IL-36RDEGs in various sample groups within the TCGA-STAD dataset were evaluated, and a heatmap was generated using the package pheatmap in R to illustrate the results (Figure 3E).

Figure 3. Differential Gene Expression Analysis. (A-C) Volcano plot of DEGs analysis between STAD and Control groups in the TCGA-STAD, GSE19826 and GSE54129 datasets. (D) The DEGs in the TCGA-STAD dataset and IL-36RGs Wayne figure. (E) Heatmap of the expression of IL-36RDEGs in the TCGA-STAD dataset.
3.3 Construction of a prognostic risk model for GC
To develop a prognostic risk model for GC, we conducted univariate Cox regression analysis using clinical information from the STAD group in the TCGA-STAD dataset combined with IL-36RDEGs. The LASSO regression analysis included all variables with a p value < 1 from the univariate analysis, and a forest plot was constructed to present the outcomes (Figure 4A). To further determine the prognostic value of the genes from the univariate Cox regression model in GC, LASSO regression analysis was conducted and a LASSO regression model was constructed. The C-index was used as the criterion for selecting LASSO variables, and the results were visually presented in the LASSO regression model (Figure 4B) and LASSO variable trajectory (Figure 4C) plots. The LASSO regression model included 9 LASSO regression model genes: IL-36A, AP1S3, IL1RAP, CARD14, IL1A, TNXB, CD276, SLC44A4 and IRAK1. The risk factors derived from the LASSO risk score were plotted utilizing the package ggplot2 in R (Figure 4D). The RiskScore was calculated as follows:

Figure 4. LASSO and Cox Regression Analysis. (A) Forest plot of the univariate Cox regression analysis of the prognosis of IL-36RDEGs. (B, C) Prognostic risk model plot and variable trajectory plot of the LASSO regression model. (D) Risk factor plot of Model Genes prognostic LASSO model. (E) Prognostic KM curves between high-risk and low-risk groups of LASSO risk score and OS in the STAD group. (F) Time-dependent ROC curve of STAD group in the TCGA-STAD dataset. p value < 0.001, extremely statistically significant.
The data demonstrated that the death rate was higher in the high-risk group than in the low-risk group, and CD276, TNXB, IL1A, CARD14 and IL1RAP were expressed at elevated levels in the high-risk group.
Subsequently, we conducted a prognostic KM curve analysis based on the LASSO risk score combined with the OS of the STAD group in the TCGA-STAD dataset using median value grouping (Figure 4E). The findings revealed a statistically significant difference in OS between the high- and low-risk groups within the STAD group of TCGA-STAD dataset (p < 0.001). Furthermore, we plotted a time-dependent ROC curve for the STAD group in TCGA-STAD dataset (Figure 4F). The study revealed that the GC prognostic risk model had an effective diagnostic capability (0.5 <AUC <0.7).
3.4 GO and KEGG enrichment analysis
Through GO and KEGG enrichment analyses, we further examined the biological processes (BP), molecular functions (MF), cellular components (CC), and the relationship between biological KEGG and GC of the nine model genes. These nine model genes were used for the GO and KEGG enrichment analyses. The findings indicated that the nine model genes were primarily abundant in BP-like cytokine-mediated signaling pathways within STAD samples and in MF, such as interleukin (IL)-1 receptor binding and growth factor receptor binding. Additionally, they were enriched in biological pathways (Kyoto Encyclopedia of Genes and Genomes), including the MAPK and NF-kappa B signaling pathways. The results of the GO and KEGG enrichment analyses are illustrated in a bubble chart (Figure 5A). Concurrently, network diagrams of BP, MF, and biological pathways were outlined according to the results obtained from the GO and KEGG enrichment analyses (Figures 5B–D).

Figure 5. GO and KEGG Enrichment Analysis for Model genes. (A) The bubble chart displays the results of the GO and KEGG enrichment analyses for the model genes: BP, MF, and KEGG. (B-D) The network diagrams show the results of the GO and pathway (KEGG) enrichment analyses for the model genes: BP (B), MF (C), and KEGG (D).
3.5 Gene set enrichment analysis
To assess how the expression levels of all genes in STAD samples influenced the high and low risk of GC, GSEA was employed to examine the connections among gene expression, the biological processes in which they participate, the affected cellular components, and the molecular functions they exert, as represented through a mountain plot (Figure 6A). The analysis demonstrated that every gene in the TCGA-STAD dataset showed significant enrichment in biological functions and signaling pathways related to cellular metabolism, signal transduction, and regulation of gene expression (Figures 6B–E).

Figure 6. GSEA for TCGA-STAD. (A) Mountain map of GSEA 4 term in the TCGA-STAD dataset. (B-E) GSEA shows that all genes are significantly enriched in HNF3B pathway (B), metabolism of fat soluble vitamins (C), folate metabolism (D) and metabolism of porphyrins (E).
3.6 Gene set variant analysis
To investigate the differences between h.all.v7.4. symbols.gmt gene set among the high- and low-risk groups in the E dataset, we conducted GSVA for all genes in the TCGA-STAD dataset. Subsequently, we selected the top 10 positively and negatively enriched pathways with p values < 0.05 and log2FC rankings and analyzed the differential expression of the 20 pathways among the high- and low-risk groups, visualizing the results via a heatmap (Figure 7A). We then performed differential validation based on the Mann–Whitney U test and created a grouped comparison chart to display the results (Figure 7B).

Figure 7. GSVA Analysis. (A, B) Heatmap (A) and group comparison of GSVA results between the high-risk and low-risk groups in the TCGA-STAD dataset (B). *p value < 0.05, statistically significant; ** p value < 0.01, highly statistically significant; *** p value < 0.001, extremely statistically significant.
3.7 Prognostic analysis of the prognostic risk model for GC
Based on the findings of the LASSO regression analysis, LASSO regression was performed to explore the correlation between the LASSO risk score and clinical prognosis. The results of the multivariate Cox regression analysis were visualized using a forest plot (Figure 8A). To further support the significance of the GC prognostic risk model, a nomogram was built using the LASSO risk score and clinical details to display the interconnections among the genes (Figure 8B). The results revealed that the LASSO risk score was significantly more effective in the GC prognostic risk model than other factors. Moreover, we conducted a prognostic calibration assessment for the GC risk model across 1, 2, and 3 years and plotted the calibration curves (Figures 8C–E). The outcomes showed that the prognostic risk model for GC had the best clinical predictive performance at the 3-year follow-up. Finally, we evaluated the clinical utility of GC prognostic risk model utilizing DCA for 1, 2, and 3 years (Figures 8F–H). The findings indicated that the ranking of the clinical predictive performance of our established multivariate Cox regression model was as follows: 3 years > 2 years > 1 year.

Figure 8. Prognostic Analysis. (A, B) Forest plot (A) and Nomogram (B) of the LASSO risk score and clinical information in the multivariate Cox regression model. (C-E) Calibration Curve for the GC prognostic risk model at 1 year (C), 2 years (D), and 3 years (E). (F-H) DCA plots for the GC prognostic risk model at 1 year (F), 2 years (G), and 3 years (H).
3.8 PPI interaction network and functional similarity analysis
Using the STRING database to analyze the PPI interaction network of the nine model genes, we retained only the genes that were connected to other nodes and designated them as hub genes for subsequent analysis. This resulted in the construction of a PPI network consisting of six hub genes (IL-36A, AP1S3, IL1RAP, CARD14, IL1A, and IRAK1), which were visualized using the Cytoscape software (Figure 9A).

Figure 9. PPI interactions network and functional similarity analysis. (A) PPI network diagram. (B) The GeneMANIA website predicts the interaction network of functionally similar genes for the Hub genes.
The correlation of the six hub genes with other genes was analyzed using the GeneMANIA database (Figure 9B). The findings indicated that the six hub genes primarily exhibited co-expression, predicted interactions, physical interactions, pathways, shared protein domains, and co-localization with other genes.
3.9 Differential expression validation of hub genes and ROC curve analysis
The differences in the expression of the six hub genes among the STAD and control groups in the TCGA-STAD dataset were analyzed and presented in a grouped comparison chart (Figure 10A). Analysis of differences demonstrated that the expression of the six hub genes were statistically significant (p < 0.001). Correlation analysis was conducted on the expression of the six hub genes in the TCGA-STAD dataset, and a correlation heatmap was generated (Figure 10B). Among these, IL1RAP and IL1A showed the greatest positive association (r = 0.43, p < 0.05). Finally, using the R package pROC, ROC curves were created from the expression of hub genes in the TCGA-STAD dataset (Figures 10C–E). The ROC curves indicated that the four hub genes (AP1S3, IL1RAP, CARD14, and IRAK1) showed a notable level of accuracy in distinguishing the STAD group from the control group based on their expression (0.7 < AUC < 0.9).

Figure 10. Differential Expression Validation and ROC Curve Analysis of Hub genes. (A) Grouped comparison chart of hub genes between STAD and Control groups in the TCGA-STAD dataset. (B) Correlation heatmap of hub genes in the TCGA-STAD dataset. (C-E) ROC curves for hub genes IL-36A and AP1S3 (C), IL1RAP and CARD14 (D), IL1A and IRAK1 (E) in the TCGA-STAD dataset. (F, G) Grouped comparison chart of hub genes between STAD and Control groups in the GSE19826 (F) and GSE54129 (G) datasets. * p value < 0.05; ** p value < 0.01; *** p value < 0.001.
The validation results are presented through group comparison graphs, demonstrating the expression difference analysis results of six hub genes in the STAD and control groups of the GSE19826 and GSE54129 (Figures 10F, G). The differential results indicated that the representation differences of the six hub genes in the dataset GSE54129 were significant and statistically meaningful.
3.10 Immunohistochemical results
Among the six hub genes, previous studies have been conducted on the roles of the other five hub genes in GC (43–51), except for IL-36A. Hence, IL-36A was selected for immunohistochemistry (IHC) experiments to detect its expression in GC. In this study, a sum of 46 STAD samples and 22 adjacent normal tissue samples (controls) were collected. The specimens were embedded with paraffin and sectioned into 0.5-μm thick pieces. The expression of IL-36A protein was detected (Figure 11A), with positive expression indicated by brownish-yellow or brown cytoplasmic granules. The positive expression rate was assessed using the ImageJ software.

Figure 11. The expression levels of IL-36A protein. (A) The expression level of IL-36A in the Control and STAD group tissues was detected using the IHC staining method (×200). (B) Quantitative analysis of the average optical density (IOD/Area) of IL-36A IHC images in the Control and STAD groups was calculated using Image-J software, and the data are presented as mean ± standard deviation. * p < 0.05.
The IHC results showed that IL-36A was primarily expressed in gastric foveolar glands and stroma in the control group tissues (0.152 ± 0.04), while in the STAD tissues, IL-36A was mainly expressed in the cytoplasm (0.181 ± 0.05). For IL-36A, a statistically significant difference was found in comparison with the control group (p = 0.0186, p < 0.05) (Figure 11B).
3.11 Immune infiltration analysis
In the expression matrix of the STAD group in the TCGA-STAD dataset, the ssGSEA algorithm was employed to assess the immune infiltration abundance of 28 types of immune cells within the high- and low-risk groups. First, a grouped comparison chart (Figure 12A) was presented to display the discrepancies in immune cell infiltration abundance across various groups. The grouped comparison chart indicated that the 20 categories of immune cells showed statistically significant differences (p < 0.05).

Figure 12. Risk Group Immune Infiltration Analysis by ssGSEA Algorithm. (A) Comparison of immune cells in the TCGA-STAD dataset between the high-risk and low-risk group. (B, C) Correlation bubble plots of immune cell infiltration abundance and hub genes in the high-risk (B) and low-risk (C). ns p value ≥ 0.05, non statistically significant; * p value < 0.05; ** p value < 0.01; *** p value < 0.001.
Finally, the correlation bubble plot showed the relationship between hub genes and levels of immune cells (Figures 12B, C). The results of the correlation bubble plot indicated that certain hub genes showed strong positive associations with specific immune cell types in the high- and low-risk groups. In the low-risk group, IL1RAP shows the strongest positive relationship with immature dendritic cells (r = 0.337, p < 0.05), whereas in the high-risk group, IL1A shows the highest positive association with neutrophils (r = 0.551, p < 0.05).
4 Discussion
GC has a complex etiology, is difficult to detect in its early stages, has a poor prognosis after treatment, and remains a malignant tumor with a high mortality rate, imposing a heavy economic and psychological burden on patients and their families. Patients with GC generally have a poor prognosis and low 5-year survival rate (1, 52). Important advancements have been made in recent years regarding prognostic biomarkers and molecular targets for GC. For example, prognostic risk models have been established using apoptosis-related molecules like p53, BCL-2, and Caspases-3, along with seven identified Anoikis-Related Long Non-Coding RNAs (ar-lncRNAs) (53–56). Additionally, such as the immune checkpoint PD-L1 (57), and some newer biomarkers like Human Epidermal Growth Factor Receptor 2 (HER2) (58, 59) and Fibroblast Growth Factor Receptor 2 (FGFR2) (60) have also emerged. And others existing prognostic models for GC that involve iron death-related genes (61), copper death-related genes (62), and mitochondrial-related gene models (63), these biomarkers not only serve as prognostic predictors but also provide new avenues for molecular targeted therapy. However, the efficacy and cytotoxic effects of these biomarkers in relation to targeted drugs remain unclear, there is a necessary to explore other prognostic model for GC. Our study presents a novel prognostic model associated with IL36-RGs, which has not been reported before. The IL-1 family is known to have connections with inflammation, immunity, and cancer (64); however, there has been no previous report documenting a GC prognostic model specifically focused on IL-36. As previously mentioned, some studies have confirmed a connection between the IL-36 family and the occurrence and development of tumors. However, research has typically focused on the regulatory mechanisms or predictive models of a specific factor corresponding to a single pathway. Disease occurrence and development result from an interplay between multiple factors. Therefore, this study constructed a prognostic model using IL-36RDEGs through various bioinformatics analysis methods. Firstly, the study identified 16 IL-36RDEGs. We utilized univariate Cox regression and LASSO regression analyses to identify nine feature genes, including IL-36A, AP1S3, IL1RAP, CARD14, IL1A, and IRAK1 to build a prognostic risk model for GC, in which patients with high-risk scores had a significantly reduced survival duration. Calibration analysis and DCA showed that the accuracy of this prognostic risk model was sufficiently good, especially in clinical predictions over 3 years, indicating that over time, the model’s predictive ability in clinical applications gradually improved.
Regarding these nine model genes, we employed various enrichment analyses to reveal their participation in important biological processes and signaling pathways. GO and KEGG functional enrichment analyses indicated that the nine model genes primarily participate in BP that promote cytokine synthesis and cytokine-mediated signaling pathways in GC. Cytokines activate a series of intracellular signaling pathways by binding to their receptors, thereby triggering changes in the cell phenotype, metabolism, and function. This study showed that the main MF of the model genes was binding to interleukin-1 receptors and growth factor receptors, and their enrichment was predominantly observed in the cytokine-cytokine receptor interaction and the MAPK signaling pathway (KEGG). The MAPK signaling pathway can regulate interactions between cells and the microenvironment, promoting the generation of new blood vessels in tumors and interactions between tumor and immune cells (65). These observations suggest that the model genes potentially influence the tumor microenvironment of GC. Studies have shown that IL-36A is an important predictor of unfavorable prognosis in patients with non-small cell lung cancer (66). IL-36A is expressed in all types of immune and non-immune cells including T cells, neutrophils, and epithelial cells. Evidence reveals that IL-36A plays a significant pro-inflammatory biological role in the communication between different cells, such as dendritic cells, neutrophils, and epithelial cells, in the course of initiating, sustaining, and amplifying inflammation. IL-36A can also activate MAPKs and the NF-κB pathway, as reflected in the GO and KEGG analyses of the model genes in this analysis. The function of IL-36A in the tumor microenvironment is also receiving increasing attention, with one of the main mechanisms being the enhancement of immune cell infiltration through upregulation of the expression of various chemokines (67). The connection between hub genes and immune cell infiltration abundance in this study showed that IL36A positively correlated with neutrophils in the low-risk group and negatively correlated with regulatory T cells in the high-risk group. Earlier research has indicated that IL-36A could be closely linked to the tumor immune microenvironment (68). Thus, IL-36A is a potentially effective target for clinical immunotherapy of GC. Our study revealed that IL-36A, as a factor in the GC prognostic risk model, is being reported for the first time. The validation of the differential analysis results demonstrated statistical significance for the Hub genes in the GSE54129 dataset. However, the differential analysis of IL-36A in the TCGA-STAD dataset showed no statistical significance, and there is limited research on its expression in GC tissues. Given the existing research background, it is crucial to further clarify its expression in GC tissues. This study detected the expression of IL-36A in GC using immunohistochemistry. The analysis demonstrated a statistically significant difference in the expression between the STAD and Control groups. In summary, this study provides a theoretical foundation for studying the mechanism of immune action of IL-36A in GC.
Research in GC indicated that AP1S3 plays a role in enhancing the development of GC (43, 44). IL1RAP shows higher expression in GC tissues than in non-cancerous tissues and is involved in the occurrence of GC and regulation of the inflammatory process (45). The findings on AP1S3 and IL1RAP were consistent with the analysis of differentially expressed hub genes in this study, and the ROC curve analysis showed good accuracy. In the diffuse GC group, elevated CARD14 expression was significantly associated with worse patient prognosis (46). Malignant cells, tumor-infiltrating immune cells, and stromal cells can express IL1A (47, 69), which is overexpressed in GC and is closely related to the clinical features of patients with GC (48, 49, 70). IRAK1 is associated with negative prognosis and invasiveness of GC (50, 51). Research findings on CARD14, IL1A, and IRAK1 suggest that they can serve as biomarkers for GC prognosis, which concurs with the outcomes of this investigation. Western blotting and immunohistochemical experiments have indicated that TNXB is overexpressed in patients with gastric adenocarcinoma with lymph node metastasis (71), which is linked to a lower survival rate in these patients and can serve as a prognostic marker. CD276 and SLC44A4 have also been confirmed as prospective targets for immune checkpoints, diagnosis, prognosis, and treatment of tumors (72–74). In summary, the study reveals that all nine model genes were associated with tumors or inflammation, highlighting their potential functions and prognostic abilities in GC. This provides strong support for the reliability of this study and is consistent with the outcomes of the prognostic risk model.
To better elucidate the interactions among model genes, we created a PPI network and identified six hub genes. Hub genes exhibit physical interactions, pathways, and co-expression. By analyzing the relationship between hub genes and differential immune cell infiltration abundance, it was observed that in the low-risk group, IL1RAP exhibited a positive relationship with the levels of immature dendritic cell infiltration. The innate immune system, especially immature dendritic cells, can be activated by recognizing DAMPs released by tumor cells through pattern recognition receptors, triggering metabolic changes in dendritic cells and playing a role in cancer immune surveillance (75). IL1A levels were positively correlated with the abundance of neutrophil infiltration in high-risk samples. Neutrophils are a major component of the tumor microenvironment. Studies have shown that neutrophils interact with GC cells, promoting their invasion and migration of GC cells through the induction of the epithelial-mesenchymal transition (EMT) and activation of the ERK pathway (76). Neutrophils may also support the development and progression of GC by facilitating angiogenesis and reducing antitumor T cell activity (77, 78), suggesting that neutrophils contribute to the occurrence and development of GC. Hub genes may influence the tumor microenvironment or patient prognosis by affecting the abundance or activity of specific immune cells. These findings provide important clues for personalization of immunotherapy.
GSEA showed that all genes within the STAD group were primarily focused on biological processes and signaling pathways related to cell development and metabolism. Changes in cell development and metabolism, which can trigger alterations in the tumor microenvironment, have gained increasing attention (79–83). These pathways may be directly connected to and closely related to the occurrence, development, and changes in the GC microenvironment, thus providing important clues for mechanistic research and potential therapeutic strategies for GC. The outcomes of the enrichment analysis of these pathways through GSVA can reveal biological processes or signaling pathways that exhibit notable variations between the high- and low-risk groups, such as angiogenesis, EMT, and inflammatory response pathways, which could be linked to the onset and progression of GC. These findings provide guidance for future research on these mechanisms.
5 Conclusion
The prognostic risk model constructed using IL-36RDEGs is an independent prognostic risk factor for GC and can be used to assess the prognosis of patients with GC. This model provides a powerful tool for clinical practice, as it can identify high-risk patients and optimize treatment strategies, thereby leading to a better overall prognosis for patients with GC. The expression of critical genes was initially validated through immunohistochemical staining. However, some limitations remain unresolved. For instance, the sample size was relatively small, and there is a deficiency of extensive, multicenter clinical samples for further external validation. Although we validated the expression of core genes and their ability to distinguish GC from normal tissues using two independent datasets, GSE19826 and GSE54129 from the GEO databases, the public GEO datasets generally lack comprehensive clinical prognostic follow-up data. Consequently, a systematic assessment of the overall prognostic predictive ability of the model using independent external cohorts is not yet possible. Therefore, the generalizability and clinical applicability of the model require further validation in independent cohorts with complete prognostic information. Additionally, some potential influencing factors may have been inadequately considered due to limitations in data sources and the analytical methods employed. The model achieved an AUC value of just 0.584 for the 3-year overall survival rate, suggesting a constrained predictive capability. It is necessary to optimize and enhance the prognostic predictive efficacy of the model in the future by increasing sample sizes and refining multi-center clinical cohort data. We intend to gather and analyze independent cohorts with clinical prognostic information to improve the model’s utility and generalizability. Overall, this study provides new evidence for the prognostic direction of GC, which can act as a potential therapeutic target for further mechanistic validation and clinical evaluation.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The studies involving humans were approved by the Ethics Committee of the Affiliated Hospital of Beihua University in Jilin City, Jilin Province, China. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
YZ: Formal analysis, Resources, Visualization, Writing – original draft. YaL: Conceptualization, Methodology, Project administration, Supervision, Writing – review & editing. XG: Conceptualization, Methodology, Supervision, Writing – review & editing. MQ: Conceptualization, Methodology, Supervision, Writing – review & editing. DW: Data curation, Formal analysis, Investigation, Writing – review & editing. NL: Data curation, Investigation, Resources, Writing – review & editing. ZL: Data curation, Investigation, Writing – review & editing. YuL: Formal analysis, Investigation, Software, Visualization, Writing – review & editing. HW: Formal analysis, Investigation, Software, Visualization, Writing – review & editing. LY: Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the Natural Science Foundation of Jilin (YDZJ202401033ZYTS, YDZJ202201ZYTS194). Department of Science and Technology of Jilin Province (20210402015GH). Education Department of Jilin Province (JJKH20220068KJ, JJKH20230080KJ). The Health Commission of Jilin Province (2022JC026). The Beihua University College Student Innovation Program (202410201212).
Acknowledgments
The authors thank the TCGA network and the GEO network for their contributions.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1566993/full#supplementary-material
References
1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36cancersin 185 countries. CA Cancer J Clin. (2021) 71:209–49. doi: 10.3322/caac.21660
2. Machlowska J, Baj J, Sitarz M, Maciejewski R, and Sitarz R. Gastric cancer: epidemiology, risk factors, classification, genomic characteristics and treatment strategies. Int J Mol Sci. (2020) 21:4012. doi: 10.3390/ijms21114012
3. Li P, Huang CM, Zheng CH, Russo A, Kasbekar P, Brennan MF, et al. Comparisonof gastric cancer survival after R0 resection in the US and China. J Surg Oncol. (2018) 118:975–82. doi: 10.1002/jso.25220
4. Chen K. Totally laparoscopic gastrectomy for gastric cancer: a systematic reviewandmeta-analysis of outcomes compared with open surgery. World J Gastroenterol. (2014) 20:15867. doi: 10.3748/wjg.v20.i42.15867
5. Li K, Zhang A, Li X, Zhang H, and Zhao L. Advances in clinical immunotherapy for gastriccancer. Biochim Biophys Acta Rev Cancer. (2021) 1876:188615. doi: 10.1016/j.bbcan.2021.188615
6. Dinarello CA. The IL-1 family of cytokines and receptors in rheumatic diseases. Nat Rev Rheumatol. (2019) 15:612–32. doi: 10.1038/s41584-019-02778
7. van de Veerdonk FL, Stoeckman AK, Wu G, Boeckermann AN, Azam T, Netea MG, et al. IL-38 binds to the IL-36 receptor and has biological effects on immune cells similar to IL-36receptor antagonist. Proc Natl Acad Sci. (2012) 109:3001–5. doi: 10.1073/pnas.1121534109
8. Neurath MF. IL-36 in chronic inflammation and cancer. Cytokine Growth Factor Rev. (2020) 55:70–9. doi: 10.1016/j.cytogfr.2020.06.006
9. O’Reilly S. Interleukin-36α is elevated in diffuse systemic sclerosis and may potentiate fibrosis. Cytokine (Philadelphia Pa). (2022) 156:155921. doi: 10.1016/j.cyto.2022.155921
10. Byrne J, Baker K, Houston A, and Brint E. IL-36 cytokines in inflammatory and Malignant diseases: not the new kid on the block anymore. Cell Mol Life Sci. (2021) 78:6215–27. doi: 10.1007/s00018-021-03909-4
11. Baker KJ, Brint E, and Houston A. Transcriptomic and functional analyses reveal a tumour-promoting role for the IL-36 receptor in colon cancer and crosstalk between IL-36 signallingand665 the IL-17/IL-23 axis. Br J Cancer. (2023) 128:735–47. doi: 10.1038/s41416-022-02083-z
12. Le N, Luk I, Chisanga D, Shi W, Pang L, Scholz G, et al. IL-36G promotes cancer-cell intrinsic hallmarks in human gastric cancer cells. Cytokine (Philadelphia Pa). (2022) 155:155887. doi: 10.1016/j.cyto.2022.155887
13. Hu M, Tong Y, Fang H, Tang J, Liu L, Hu Y, et al. IL36 indicating good prognosis inhuman Hepatocellular Carcinoma. J Cancer. (2020) 11:6248–55. doi: 10.7150/jca.47106
14. Song Y, Chu H, Liu F, Guo W, Gao N, Chen C, et al. The pro-tumor biological functionof IL- 36α Plays an important role in the tumor microenvironment of HCC. Cancer Manag Res. (2023) 15:895–904. doi: 10.2147/CMAR.S407123
15. Baker K, O’Donnell C, Bendix M, Keogh S, Byrne J, O’Riordain M, et al. IL-36signalling enhances a pro-tumorigenic phenotype in colon cancer cells with cancer cell growth restrictedby administration of the IL-36R antagonist. Oncogene. (2022) 41:2672–84. doi: 10.1038/s41388-022-67702281-2
16. Sullivan GP, Henry CM, Clancy DM, Mametnabiev T, Belotcerkovskaya E, Davidovich P, et al. Suppressing Suppressing IL-36-driven inflammation using peptide pseudosubstrates for neutrophil proteases. Cell Death Dis. (2018) 9:378. doi: 10.1038/s41419-018-0385-4
17. Roy S, Fitzgerald K, Lalani A, Lai C, Kim A, Kim J, et al. Autonomous IL-36Rsignalingin neutrophils activates potent antitumor effector functions. J Clin Invest. (2023) 133:1–17. doi: 10.1172/JCI162088
18. Zhao X, Chen X, Shen X, Tang P, Chen C, Zhu Q, et al. IL-36β Promotes CD8+TCell activation and antitumor immune responses by activating mTORC1. Front Immunol. (2019) 10:1803. doi: 10.3389/fimmu.2019.01803
19. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. (2016) 44:e71. doi: 10.1093/nar/gkv1507
20. Goldman MJ, Craft B, Hastie M, Repečka K, McDade F, Kamath A, et al. Visualizingand interpreting cancer genomics data via the Xena platform. Nat Biotechnol. (2020) 38:675–8. doi: 10.1038/s41587-020-0546-8
21. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. (2012) 41:D991–5. doi: 10.1093/nar/gks1193
22. Davis S and Meltzer PS. Geoquery: a bridge between the gene expression omnibus (GEO) and Bioconductor. Bioinformatics. (2007) 23:1846–7. doi: 10.1093/bioinformatics/btm254
23. Wang Q, Wen Y, Li D, Xia J, Zhou C, Yan D, et al. Upregulated INHBA expression is associated with poor survival in gastric cancer. Med Oncol (Northwood London England). (2012) 29:77–83. doi: 10.1007/s12032-010-9766-y
24. Leek JT, Johnson WE, Parker HS, Jaffe AE, and Storey JD. The sva package for removingbatch effects and other unwanted variation in high-throughput experiments. Bioinformatics. (2012) 28:882–3. doi: 10.1093/bioinformatics/bts034
25. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007
26. Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, et al. The geneCards suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinf. (2016) 54:1.30.1–1.30.33. doi: 10.1002/cpbi.5
27. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, and Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. (2011) 27:1739–40. doi: 10.1093/bioinformatics/btr260
28. Therneau T, Lumley T, and Elizabeth A. A package for survival analysis in R. In: R package version (2024). p. 34–0.
29. Engebretsen S and Bohlin J. Statistical predictions with glmnet. Clin Epigenet. (2019) 11:123. doi: 10.1186/s13148-019-0730-1
30. Rich JT, Neely JG, Paniello RC, Voelker CCJ, Nussenbaum B, and Wang EW. Apractical guideto understanding Kaplan-Meier curves. Otolaryngology–Head Neck Surg. (2010) 143:331–6. doi: 10.1016/j.otohns.2010.05.007
31. Park SH, Goo JM, and Jo CH. Receiver operating characteristic (ROC) curve: practical reviewfor radiologists. Korean J Radiol. (2004) 5:11–8. doi: 10.3348/kjr.2004.5.1.11
32. Mi H, Muruganujan A, Ebert D, Huang X, and Thomas PD. Panther version 14: more genomes, PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res. (2019) 47:D419–26. doi: 10.1093/nar/gky1038
33. Kanehisa M and Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic AcidsRes. (2000) 28:27–30. doi: 10.1093/nar/28.1.27
34. Yu G, Wang L, Han Y, and He Q. Clusterprofiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
35. Luo W and Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization. Bioinformatics. (2013) 29:1830–1. doi: 10.1093/bioinformatics/btt285
36. Hanzelmann S, Castelo R, and Guinney J. GSVA: gene set variation analysis for microarrayand RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7
37. Wu J, Zhang H, Li L, Hu M, Chen L, Xu B, et al. A nomogram for predicting overall survival in patients with low-grade endometrial stromal sarcoma: a population-based analysis. Cancer Commun (Lond). (2020) 40:301–12. doi: 10.1002/cac2.12067
38. Van Calster B, Wynants L, Verbeek JFM, Verbakel JY, Christodoulou E, Vickers AJ, et al. Reporting and interpreting decision curve analysis: a guide for investigators. Eur Urol. (2018) 74:796–804. doi: 10.1016/j.eururo.2018.08.038
39. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. String v11: protein-protein association networks with increased coverage, supporting functional discovery ingenome-wide experimental datasets. Nucleic Acids Res. (2019) 47:D607–13. doi: 10.1093/nar/gky1131
40. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
41. Franz M, Rodriguez H, Lopes C, Zuberi K, Montojo J, Bader GD, et al. GeneMANIAupdate 2018. Nucleic Acids Res. (2018) 46:W60–4. doi: 10.1093/nar/gky311
42. Xiao B, Liu L, Li A, Xiang C, Wang P, Li H, et al. Identification and verification of immune-related gene prognostic signature based on ssGSEA for osteosarcoma. Front Oncol. (2020) 10:607622. doi: 10.3389/fonc.2020.607622
43. Xu L, Liu D, and Wang X. Linc00562 drives gastric cancer development by regulating miR-4636-AP1S3 axis. Korean J Physiol Pharmacol. (2023) 27:197–208. doi: 10.4196/kjpp.2023.27.3.197
44. Ye T, Cheng Y, and Li C. Adaptor protein complex 1 sigma 3 is highly expressed in gliomaand could enhance its progression. Comput Math Methods Med. (2021) 2021:1–11. doi: 10.1155/2021/5086236
45. Lv Q, Xia Q, Li A, and Wang Z. The potential role of IL1RAP on tumor microenvironment-related inflammatory factors in stomach adenocarcinoma. Technol Cancer Res Treat. (2021) 20:1180463896. doi: 10.1177/1533033821995282
46. Carino A, Graziosi L, Marchiano S, Biagioli M, Marino E, Sepe V, et al. Analysis of gastric cancer transcriptome allows the identification of histotype specific molecular signatureswith prognostic potential. Front Oncol. (2021) 11:663771. doi: 10.3389/fonc.2021.663771
47. Malik A and Kanneganti TD. Function and regulation of IL-1α in inflammatory diseases andcancer. Immunol Rev. (2018) 281:124–37. doi: 10.1111/imr.12615
48. Kemik O, Kemik AS, Begenik H, Erdur FM, Emre H, Sumer A, et al. The relationshipamong acute-phase responce proteins, cytokines, and hormones in various gastrointestinal cancer types patients with cachectic. Hum Exp Toxicol. (2012) 31:117–25. doi: 10.1177/0960327111417271
49. Uefuji K, Ichikura T, and Mochizuki H. Increased expression of interleukin-1αand cyclooxygenase-2 in human gastric cancer: a possible role in tumor progression. Anticancer Res. (2005) 25:3225–30. doi: 10.1007/s10120-022-01297-7
50. Yamanaka N, Morisaki T, Tsuneyoshi M, Tanaka M, Katano M, Nakashima H, et al. Interleukin 1beta enhances invasive ability of gastric carcinomathrough nuclear factor-kappaB activation. Clin Cancer Res. (2004) 10:1853–9. doi: 10.1158/1078-7730432.CCR-03-0300
51. Kogo R, Mimori K, Tanaka F, Komune S, and Mori M. Clinical significance ofmir-146a ingastric cancer cases. Clin Cancer Res. (2011) 17:4277–84. doi: 10.1158/1078-0432.CCR-10-2866
52. Smyth EC, Nilsson M, Grabsch HI, van Grieken NC, and Lordick F. Gastric cancer. Lancet. (2020) 396:635–48. doi: 10.1016/S0140-6736(20)31288-5
53. Huang KH, Fang WL, Li AF, Liang PH, Wu CW, Shyr YM, et al. Caspase-3, a keyapoptotic protein, as a prognostic marker in gastric cancer after curative surgery. Int J Surg. (2018) 52:258–63. doi: 10.1016/j.ijsu.2018.02.055
54. Min K, Kim D, Son BK, Kim DH, Kim E, Seo J, et al. A high ki67/BCL2 index couldPredict lower disease-free and overall survival in intestinal-type gastric cancer. Eur Surg Res. (2017) 58:158–68. doi: 10.1159/000448945
55. Qiu Z, Qin R, Zhang Z, Zhang T, Zhang Z, Qiao C, et al. Expression of p53 as a biomarkerin determining response to apatinib for advanced gastric cancer. Front Oncol. (2023) 13:1203980. doi: 10.3389/fonc.2023.1203980
56. Meng W, Guo J, Huang L, Zhang Y, Zhu Y, Tang L, et al. Anoikis-related long non-codingrna signatures to predict prognosis and immune infiltration of gastric cancer. Bioengineering. (2024) 11:893. doi: 10.3390/bioengineering11090893
57. Böger C, Behrens H, Mathiak M, Krüger S, Kalthoff H, and Röcken C. PD-L1 is an independent prognostic predictor in gastric cancer of Western patients. Oncotarget. (2016) 7:24269–83. doi: 10.18632/oncotarget.8169
58. Kurokawa Y, Matsuura N, Kimura Y, Adachi S, Fujita J, Imamura H, et al. Multicenter large-scale study of prognostic impact of HER2 expression in patients with resectable gastriccancer. Gastric Cancer. (2015) 18:691–7. doi: 10.1007/s10120-014-0430-7
59. Kim KC, Koh YW, Chang H, Kim TH, Yook JH, Kim BS, et al. Evaluation of HER2protein expression in gastric carcinomas: comparative analysis of 1,414 cases of whole-tissue sections and 595 cases of tissue microarrays. Ann Surg Oncol. (2011) 18:2833–40. doi: 10.1245/s10434-011-7991695-2
60. Ahn S, Lee J, Hong M, Kim ST, Park SH, Choi MG, et al. FGFR2 in gastric cancer: protein overexpression predicts gene amplification and high H-index predicts poor survival. ModPathol. (2016) 29:1095–103. doi: 10.1038/modpathol.2016.96
61. Liu D, Peng J, Xie J, and Xie Y. Comprehensive analysis of the function of helicobacter-associated ferroptosis gene YWHAE in gastric cancer through multi-omics integration, molecular docking, and machine learning. Apoptosis. (2024) 29:439–56. doi: 10.1007/s10495-023-01916-3
62. Wang J, Qin D, Tao Z, Wang B, Xie Y, Wang Y, et al. Identification of cuproptosis-related subtypes, construction of a prognosis model, and tumor microenvironment landscape ingastric cancer. Front Immunol. (2022) 13:105693247. doi: 10.3389/fimmu.2022.105693247
63. Chang J, Wu H, Wu J, Liu M, Zhang W, Hu Y, et al. Constructing a novel mitochondrial-related gene signature for evaluating the tumor immune microenvironment and predicting survival in stomach adenocarcinoma. J Transl Med. (2023) 21:191. doi: 10.1186/s12967-023-04033-6
64. Sun R, Gao DS, Shoush J, and Lu B. The IL-1 family in tumorigenesis and antitumor immunity. Semin Cancer Biol. (2022) 86:280–95. doi: 10.1016/j.semcancer.2022.05.002
65. Brägelmann J, Lorenz C, Borchmann S, Nishii K, Wegner J, Meder L, et al. MAPK-pathway inhibition mediates inflammatory reprogramming and sensitizes tumors to targeted activationof innate immunity sensor RIG-I. Nat Commun. (2021) 12:5505. doi: 10.1038/s41467-02125728-8
66. Xie X, Hu H, He J, Liu Y, Guo F, Luo F, et al. Interleukin-36alpha suppresses growthof non-small cell lung cancer in vitro by reducing angiogenesis. FEBS Open Bio. (2021) 11:1353–63. doi: 10.1002/2211-5463.13141
67. Wei X, Yao Y, Wang X, Sun J, Zhao W, Qiu L, et al. Interleukin-36α inhibits colorectal cancer metastasis by enhancing the infiltration and activity of CD8+ T lymphocytes. Int Immunopharmacol. (2021) 100:108152. doi: 10.1016/j.intimp.20
68. Vela-Patiño S, Salazar MI, Taniguchi-Ponciano K, Vadillo E, Gomez-Apo E, Escobar-España A, et al. The immune microenvironment landscape of pituitary neuroEndocrine tumors, a transcriptomic approach. Genes (Basel). (2024) 15:531. doi: 10.3390/genes15050531
69. Sakurai T, He G, Matsuzawa A, Yu G, Maeda S, Hardiman G, et al. Hepatocyte necrosis induced by oxidative stress and IL-1α release mediate carcinogen-induced compensatory proliferationand liver tumorigenesis. Cancer Cell. (2008) 14:156–65. doi: 10.1016/j.ccr.2008.06.016
70. Zhou L, Niu Z, Wang Y, Zheng Y, Zhu Y, Wang C, et al. Senescence as a dictator of patient outcomes and therapeutic efficacies in human gastric cancer. Cell Death Discov. (2022) 8:13. doi: 10.1038/s41420-021-00769-6
71. Lu X, Fu Y, Gu L, Zhang Y, Liao AY, Wang T, et al. Integrated proteome and phosphoproteome analysis of gastric adenocarcinoma reveals molecular signatures capable of stratifyingpatient outcome. Mol Oncol. (2023) 17:261–83. doi: 10.1002/1878-0261.1336133
72. Liu S, Liang J, Liu Z, Zhang C, Wang Y, Watson AH, et al. The role of CD276 in cancers. Front Oncol. (2021) 11:654684. doi: 10.3389/fonc.2021.654684
73. Cheng M, Chen S, Li K, Wang G, Xiong G, Ling R, et al. CD276-dependent efferocytosisby tumor-associated macrophages promotes immune evasion in bladder cancer. Nat Commun. (2024) 839:15. doi: 10.1038/s41467-024-46735-5
74. Kang W, Zhang M, Wang Q, Gu D, Huang Z, Wang H, et al. The SLC family are candidate diagnostic and prognostic biomarkers in clear cell renal cell carcinoma. BioMed Res Int. (2020) 2020:1–17. doi: 10.1155/2020/1932948
75. Hernandez C, Huebener P, and Schwabe RF. Damage-associated molecular patterns incancer: a double-edged sword. Oncogene. (2016) 35:5931–41. doi: 10.1038/onc.2016.104
76. Zhang W, Gu J, Chen J, Zhang P, Ji R, Qian H, et al. Interaction with neutrophils promotes gastric cancer cell migration and invasion by inducing epithelial-mesenchymal transition. Oncol Rep. (2017) 38:2959–66. doi: 10.3892/or.2017.5942
77. Wang T, Zhao Y, Peng L, Chen N, Chen W, Lv Y, et al. Tumour-activated neutrophils ingastric cancer foster immune suppression and disease progression through GM-CSF-PD-L1 pathway. Gut. (2017) 66:1900–11. doi: 10.1136/gutjnl-2016-313075
78. Li TJ, Jiang YM, Hu YF, Huang L, Yu J, Zhao LY, et al. Interleukin-17-producing neutrophils link inflammatory stimuli to disease progression by promoting angiogenesis in gastric cancer. Clin Cancer Res. (2017) 23:1575–85. doi: 10.1158/1078-0432.CCR-16-0617
79. Shibu MA, Huang CY, and Ding DC. Comparison of two hepatocyte differentiation protocolsin human umbilical cord mesenchymal stem cells: in vitro study. Tissue Cell. (2023) 83:102153. doi: 10.1016/j.tice.2023.102153
80. Jeyakumar SM. Micronutrient deficiency in pulmonary tuberculosis-perspective onHepatic drug metabolism and pharmacokinetic variability of first-line anti- tuberculosis drugs: special reference to fat-soluble vitamins A, D, & E and nutri-epigenetics. Drug Metab Lett. (2021) 14:166–76. doi: 10.2174/1872312814999211130093625
81. Avanzino BC, Prabhakar K, Dalvi P, Hartstein S, Kehm H, Balasubramani A, et al. AT-cell engaging bispecific antibody with a tumor-selective bivalent folate receptor alpha bindingarmfor the treatment of ovarian cancer. Oncoimmunology. (2022) 11:2113697. doi: 10.1080/2162402X.2022.2113697
82. Davis CD and Hord NG. Nutritional “omics” technologies for elucidating the role(s) of bioactive food components in colon cancer prevention. J Nutr. (2005) 135:2694–7. doi: 10.1093/jn/135.11.2694
Keywords: GC, IL-36RDEGs, prognosis model, risk score, prediction
Citation: Zhang Y, Liu Y, Guan X, Qu M, Wu D, Liu N, Lin Z, Liu Y, Wang H and Yang L (2025) IL-36-related genes predict prognosis of gastric cancer. Front. Oncol. 15:1566993. doi: 10.3389/fonc.2025.1566993
Received: 17 February 2025; Accepted: 30 May 2025;
Published: 18 June 2025.
Edited by:
Alessandro Mangogna, University of Udine, ItalyCopyright © 2025 Zhang, Liu, Guan, Qu, Wu, Liu, Lin, Liu, Wang and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Lijuan Yang, eWFuZ2xpanVhbkBiZWlodWEuZWR1LmNu