SEMA3D Plays a Critical Role in Peptic Ulcer Disease-Related Carcinogenesis Induced by H. pylori Infection

Introduction

Helicobacter pylori (HP) is a Gram-negative gastrointestinal bacterium that infects nearly half of the world population.¹ HP infection is associated with the occurrence and progress of peptic ulcer disease (PUD) and gastrointestinal cancer (GC).^2,3 The causes of GC, and its connection to PUD, remain poorly defined. Recently, HP-induced chronic inflammation has been shown to play an important role in GC occurrence and development.⁴ Prior studies defined several chemokines and cytokines involved in inflammation of the gastric epithelium.⁵

The cell microenvironment is the environment where tumor cells originate and develop. This region also consists of stromal cells, the tumor vascular system, immune cells, the extracellular matrix (ECM), and the acidic and hypoxic environment of the tumor.^6,7 Immune cells are the major cell types in the microenvironment and release many chemokines and cytokines that dictate disease outcomes in response to infection.⁸ During the development of peptic ulcers caused by HP infection, NF-kB signaling is activated by inflammatory factors like IL-8.⁹

In recent years, bioinformatics tools and software have been developed to quickly explore differentially expressed target genes and identify hub genes that contribute to disease progression.¹⁰ The CIBERSORT algorithm is a newly developed tool to assess the association between immune cell landscapes in the cell microenvironment using the existing 22 immune cell signatures.¹¹ This method has been successfully used to identify prognostic immune markers in lung, breast, and gastric cancer.^12,13 The present study used the CIBERSORT algorithm to calculate the proportions of 22 immune cells that infiltrate the ECM during PUD and GC, based on The Cancer Genome Atlas (TCGA) (for early-stage GC) and the Gene Expression Omnibus GEO (for PUD) databases. Inflammation and differentially expressed genes (DEGs) associated with tumors were also assessed and SEMA3D was shown to correlate with the development of peptic ulcers and the immune cell signatures. SEMA3D is a member of Class-3 semaphorins (SEMA3s), which are reported to play pivotal roles in immune response, angiogenesis, apoptosis, cell migration, and local and metastatic cancer spread in pan-cancer.^14,15 SEMA3E, a gene in the same family as SEMA3D, is a lymph node metastasis-related gene expressed in gastric cancer.¹⁶ SEMA3E deficiency dysregulates many immune cell functions both directly and indirectly.^17,18 Prior studies indicated that Semaphorin 3D and 3E have similar cellular functions, however, the exact role of SEMA3D in gastric cancer remains poorly understood.^19,20 Inappropriately regulated gastric immune responses to HP in the cell microenvironment are critical to the development of gastroduodenal disease and responses to treatment.⁵ For example, CD4+T cell-derived IFN-γprovides the key stimulus for the development of gastric premalignant lesions that progress to GC.⁹ SEMA3D also contributes to CD4+T cell infiltration in osteoarthritis joints.²¹ Other immune cells, including macrophages, dendritic cells (DCs), B cells, and gastric epithelial cells (GECs) contribute to the mucosal response to HP infection.^18,22 DCs affect the Treg/Th17 balance induced by HP infection and indirectly activate T cells.^19,23 In the current study, SEMA3D expression was primarily found in DCs from PUD samples. Similarly, prior research indicates that SEMA3E regulates DC function.^17,20 Thus, it was hypothesized that SEMA3D contributes to gastric epithelium carcinogenesis by regulating immune cell infiltration. He findings reported here indicate that SEMA3D may play an essential role in the cell microenvironment and could serve as a promising prognostic biomarker for the malignant transformation of peptic ulcers.

Materials and Methods

Data Preparation

Gene expression data from 118 early-stage GC samples and 32 PUD gastric mucosa (uninfected or HP-infected) were downloaded from TCGA and GEO databases, respectively (Table 1. GC patients who were diagnosed with stages I or II according to the 6th and 7th editions of the AJCC Cancer Staging Manual were included in this study. Clinical information for each patient was obtained from the TCGA database following TCGA publication guidelines and data access policies. Patients were excluded if they had recurrent GC, therapies performed before admission, other observed clinical disorders, or other GC clinical stages. PUD gene expression information was obtained from the GSE60427 dataset. The microarray platform for GSE60427 was GPL1707. Eight mucosa tissue samples were included in the normal group (GSM1479654, GSM1479655, GSM1479656, GSM1479657, GSM1479670, GSM1479671, GSM1479672 and GSM1479673) and 24 samples were included in the HP+ group (GSM1479658, GSM1479659, GSM1479660, GSM1479661, GSM1479662, GSM1479663, GSM1479664, GSM1479665, GSM1479666, GSM1479667, GSM1479668, GSM1479669, GSM1479674, GSM1479675, GSM1479676, GSM1479677, GSM1479678, GSM1479679, GSM1479680, GSM1479681, GSM1479682, GSM1479683, GSM1479684 and GSM1479685). The 32 samples used for microarray analysis were selected from 293 patient subjects. All the patients provided written informed consent and the protocols were approved by the ethics committees of Oita University (Japan). Patients with PUD and GC were identified by endoscopy. Gastritis was defined as HP gastritis in the absence of peptic ulcers or gastric malignancy. Patients with a history of partial gastric resection or who had received HP eradication therapy or treatment with antibiotics, bismuth-containing compounds, H2-receptor blockers, or proton pump inhibitors within four weeks prior to the study were excluded. All the clinical information was obtained from the previous papers.¹⁴ The number of included and excluded subjects in the study was summarized in Flow Chart (Figure S1). The protocols described above were approved by the ethics committees of Affiliated Hospital of Hebei University (AHHU20211029).

Table 1 The KEGG Enrichment Pathway List

DEG Identification

Gene expression profiles were screened using the R package, and DEGs were identified in both groups. Based on PUD and GC integrated analysis, a common gene set was identified for the two groups. DEGs were determined based on an absolute value of log2 fold change (|log2FC|) >1 and a false discovery rate (FDR) <0.05. Heatmaps of DEGs were drawn using the “pheatmap” package in R and the common differentially expressed genes between the two datasets were determined using the “Venn diagrams” in R. Immune-related genes were extracted from the DEGs after KEGG enrichment analysis. Immune-related genes were downloaded from the IMMPORT database (https://www.immport.org/home).

Pathway Enrichment and Annotation

Enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway for DEGs was performed using the KOBAS online tool (http://kobas.cbi.pku.edu.cn/index.php). KEGG analysis showed DEG enrichment in the signaling pathways.

PPI Network Analysis

The immune-associated DEGs were then used for PPI analysis. The DEG PPI network was constructed using the STRING database. Nodes with the confidence of interactive relationship > 0.7 were defined as the threshold. Subsequently, CytoHubba was utilized to identify the top ten hub genes and Venn diagrams were used to visualize common genes between the top 20 hub genes from GC and PUD patients, respectively.

Gene Set Enrichment Analysis (GSEA) and CIBERSORT Algorithm

To explore the role of the ten hub genes in regulating the cell immunological micro-environment, the CIBERSORT algorithm was applied to assess the proportion of immune cells that infiltrated the ECM using the LM22 signature. The LM22 signature consisting of 547 genes was used to identify 22 types of infiltrating immune cells. The cell fraction of the PUD and GC datasets was identified. Difference and correlation analyses were performed to explore the correlation between SMA3D expression and the types of immune cells that infiltrated the GC and PUD microenvironment (P < 0.05).

Clinicopathological Characteristics Analysis and Survival Analysis

GC patients were classified into a high- and low-group based on SMA3D and NRP2 expression levels. The association between clinicopathological characteristics and SMA3D and NRP2 expression was evaluated. K-M plotter was used to plot survival curves, which were compared using the Log rank test. Univariate and multivariate analyses were performed using the Cox proportional hazards model to investigate the influence of genotypes on cancer risk. P<0.05 was regarded as statistically significant.

Statistical Analysis

Statistical analyses were performed using R software (version 4.0.2) and GraphPad Prism. All statistical methods and appropriate R packages were described throughout the study. Hypergeometric test/Fisher’s exact test was employed to perform KEGG enrichment analysis. Method proposed by Benjamini & Hochberg was used to control for the false discovery rate (FDR).²⁴ ANOVA (one-way) and the Tukey’s test were used to compare multiple groups. Chi-square was used to analyze the correlation between SEMA3D and NRP2 expression in GC tissues and the patient clinical data. P and q<0.05 were considered statistically significant.

Results

DEGs Identification and Enrichment Analysis

Differential analysis was performed to determine immune-related DEGs. Heatmaps showed the differential gene expression profiles of GC and PUD patients (Figure 1A and B). A total of 6032 and 2032 DEGs were identified between the stromal low-score and high-score groups in each data set, respectively. The threshold of difference was |log2FC| >1 and FDR <0.05. In both the GC and PUD groups, 437 DEGs were identified as common DEGs using Venn diagrams (Figure 1C). The common DEGs were overlapped with immune-related gene sets from the IMMPORT database and 49 immune-related DEGs were grouped to identify their primary functions (Figure 1D). The top three KEGG enrichment scores indicated that these DEGs were enriched in the cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, and Rap1 signaling pathway pathways (Figure 1E). The complete pathways list of the KEGG enrichment is shown in Table 1.

Figure 1 DEGs commonly expressed in three datasets. (A and B) Hierarchical clustering heatmap of upregulated and downregulated DEGs in the PUD or GC groups from each dataset (green points) and genes without significance (black points). The differences threshold was set as |log2FC| >1.0 and adjusted P-value <0.05. The Venn diagram shows 437 DEGs (C) and 49 immune-related DEG (D) commonly expressed across the two datasets. (E) KEGG enrichment analysis of immune-related DEGs.

PPI Network Construction

The STRINGs database was used to assess the interaction between immune-associated DEGs (Figure 2A). A PPI network was then constructed with a confidence of interactive relationship >0.7 as the threshold. CytoHubba, a plug-in of the Cytoscape software, was used to screen the hub gene through three terms of degrees, closeness, and betweenness. Seven hub genes were identified as described in the methods (Figure 2B and C).

Figure 2 PPI networks and hub gene analysis of commonly expressed DEGs in the immune-related dataset. (A) PPI networks constructed by the STRINGs. (B) Major PPI network analysis of the top 10 hub genes using Cytohubba software by three methods. The node color reflects the degree of connectivity. (C) The Venn diagram of the three methods.

Correlation of the Survival and Clinicopathological Characteristics with Hub Gene Expression

Of the seven hub genes, VEGFA, EPO, SPP1, IFNG, and PLXNB3 were closely related to GC progression. NRP2 and SEMA3D were selected for survival comparative analysis in the GC group. The Kaplan-Meier survival curve showed that GC patients with low expression of SEMA3Dlow had a better overall survival rate than those with high expression (Figure 3; P<0.05; adjust HR=2.446, 95% CI 1.225–4.882). In contrast, differences in NRP2 expression did not have much effect on overall survival (Figure 3B). The correlation between SEMA3D expression and clinical characteristics was assessed using COX analysis. SEMA3D expression was closely correlated with advanced disease stages but not with TNM classification, indicating that SEMASD merits a higher clinical prognostic value (Table 2, *p<0.05).

Table 2 Cox Regression Analysis of Many Clinical-Pathological Characteristic in GC Dataset with SEMAD3

Figure 3 Correlation of NRP2 (A) and SEMA3D (B) expression with GC patient survival. COX analysis was performed to get an adjusted HR: (SEMA3D: P=0.01292; adjust HR=2.446, 95% CI 1.225–4.882), (NRP2: P=0.19287; adjust HR=1.313, 95% CI 10.6540–2.635).

Immune Cell Infiltration Analysis and Correlation Analysis

The CIBERSORT algorithm with 22 immune cells signature was employed to perform immune cell infiltration analysis using GC and PUD tissue samples. A higher fraction of active mast cells were found in the GC group than in the PUD group (Figure 4A). As shown in Figure 4B, SEMA3D expression and active mast cells showed a significant positive correlation in the GC group, but SEMA3D expression and DCs were positively correlated in the PUD group.

Figure 4 Immune cell infiltration analysis and correlation analysis. (A) Violin plot showing significant changes in immune cell infiltration in GC compared with PUD groups (P-value <0.05). (B) Correlation between gene expression and the relative percentages of immune cells in PUD and GC tissue.(C) The expression value of the two hub genes in different immune cells.

Discussion

The primary goal of this study was to explore significant hub genes associated with the malignant transformation of PUD into early GC. Seven significant hub genes were identified using bioinformatics. The SEMA3D gene was found to correlate with advanced clinicopathological stages of GC and patient survival.

KEGG pathway enrichment analysis showed that DEGs regulate many types of immune response in clinical tissue. SEMA3s always require additional neuropilin (NRP) receptors to bind VEGF, and the VEGF/SEMA3s balance is a prognostic marker of disease.²⁵ Results from this study showed that high expression of SEMA3D and NRP2 correlated with activation of the non-canonical VEGF pathway, while VEGFA signaling was inhibited. Given the TLR4 was shown to mediate CD8+T cell activation during particular innate immune responses to disease, VEGF pathway was regarded as a critical regulator in pro-inflammatory responses.^21,26,27

SEMA3D, which encodes a semaphorin III family secreted protein, is a critical regulator of neuron development and diverse tumorigenic processes like proliferation, invasion, and angiogenesis.^19,28,29 Abnormal SEMA3D expression is associated with a poor prognosis in many nervous system diseases and cancers.^{20,25,30–32} Decreased SEMA3D expression in gastrointestinal tumors correlates significantly with colorectal cancer progression while overexpression is a favorable prognostic factor for survival.³³ SEMA3D is reported to participate in the recruitment of immune cells to the disease site.³⁴ Results from this study showed that SEMA3D was significantly correlated with PUD patient outcomes. These findings imply that SMA3D could serve as a potential biomarker for early diagnosis of GC.

PPI and KEGG analyses showed that SEMA3D was involved in regulating immune-related pathways and the extracellular microenvironment of GC. The proportions of immune cells that infiltrated the cellular microenvironment were estimated using the CIBERSORT algorithm and NRP2 and SEMA3D were expressed in similar cell types. SEMA3D expression was primarily correlated with three infiltrating immune cell types, CD4+T cells, DCs, and mast cells. Interestingly, SEMA3D and NER2 expression were enriched in CD4+T cells and DCs from the PUD samples but primarily expressed in CD4+T cells and mast cells from the GC samples. Different DC subsets can differentially regulate T cell function. In PUD samples, DCs primarily functioned to induce T and B cell activation.²³ Mast cell function during cancer remains unclear, however. Recent studies show that mast cells promote gastric tumor cancer by releasing angiogenic cytokines.³⁵ Results from this study suggest that while SEMA3D expression in DCs from PUD tissue samples may help them to subvert the host immune response by activating T cells, SEMASD expression in mast cells may promote tumorigenesis. Similar heterogeneous functions of other SEMA3s are reported in other cancers.³⁶ A diagram that summarizes the findings of this study and hypothesizes how SEMA3D expression impacts DC and mast cell function is shown in Hypothetic Diagram (Figure S2). The detailed molecular mechanism of how this occurs requires additional study.

Results from this study defined seven hub genes associated with PUD-related carcinogenesis, provided strong evidence that SEMA3D correlates with tumor-related immune activation or dysfunction, and provided a new direction to study how hub gene functions during PUD and GC. However, this study does not describe the detailed mechanism by which hub genes participate in DC and mast cell function during PUD inflammation or the potential relationship between these genes and HP infection. HP-induced PUD is associated with gastric cancer, but there are few biomarkers that aid disease prognosis in clinical practice. In this study, SEMASD was defined as a potential prognostic molecule for PUD and GC, though its mechanism of action and clinical value require further research.

Conclusions

Using comprehensive bioinformatics, this study found that the hub gene, SEMA3D, was associated with the infiltration of immune cells, in particular DCs and mast cells, into PUD and GC tissue samples. Additional research on how SEMASD impacts immune cell function in the PUD and GC dataset will help to elucidate the mechanism of malignant transformation during PUD.