A total of 30 tumor tissues [10 HP-negative (−) and 20 HP-positive (+)] were collected from patients with gastric cancer who underwent gastrectomy at the Department of Gastrointestinal Surgery, Union Hospital of Huazhong University of Science and Technology (Wuhan, China). All procedures were approved by the Ethics Committee of the Union Hospital of Huazhong University of Science and Technology, and conducted according to Declaration of Helsinki principles. Prior written and informed consent was obtained from each patient or their guardians. The clinicopathological information of the patients is shown in Table SI.

Genes and miRNA expression profiles, as well as the corresponding clinical information, from patients with gastric cancer, were downloaded from The Cancer Genome Atlas (TCGA) database using TCGA biolinks package (19). Samples containing clear HP infection information were included. Samples with unavailable HP infection information, samples of paracancer tissue and samples of normal tissue were removed. The HP status of the tissue samples from the enrolled patients was confirmed using a 13C-urea breath test with 75 mg urea (UREA 13C breath test Heliforce kit; Beijing Richen-Force Science & Technology Co. Ltd.). A total of 141 gastric samples [16 cases were HP(+) and 125 cases were HP(−)] were selected for the following bioinformatics analysis. The differentially expressed genes (DEGs) with significant differential expression between the HP(+) and HP(−) groups were selected. A gene would be selected using the following conditions: i) Expression level 1.2 times higher in 20% of the total number of samples compared with the median expression level in the total number of samples; and ii) the variance in the expression level in every sample was higher than the median level. Subsequently, log2 transformation and Z correction were performed to normalize the expression value of each gene.

Limma package (20) was applied to compute DEGs and DEMs based on the aforementioned normalized mRNA and miRNA expression data. Bayes test was used to identify DEGs with ≥2 fold-change and a P-value cutoff of 0.05 was defined as statistically significant.

Genes and gene products were annotated using GO analysis, which is a common and useful method for identifying characteristic biological attributes of transcriptome data or from high-throughput genomic data. KEGG serves as a knowledge base for systematic analysis of gene functions, and linkage between genomic information and higher-order functional information. Comprehensively mapping a set of genes to the associated biological annotation is an important foundation for the success of the gene functional analysis of any high-throughput data. KEGG pathway analysis and GO enrichment was performed using the ClusterProfiler package (21) to analyze the DEGs at the functional level. P<0.05 was considered to indicate a statistically significant difference.

The PPI network was investigated using the Search Tool for the Retrieval of Interacting Genes (STRING) online database (22). STRING (v11.0) was used, which covers 24,584,628 proteins from 5,090 organisms. The DEGs were mapped using the STRING database to determine the interactive associations among DEGs, and then, the experimentally validated interactions with a significant combined score >0.4 were selected for further analysis. PPI networks were then constructed with the CytoScape software (v3.6.1) (23). Subsequently, the sub-network modules were screened using the plug-in molecular complex detection (MCODE) in CytoScape. The following screening criteria was used: MCODE score >3 and node number >6. In addition, the function and pathway enrichment analysis of DEGs in the sub-module was also determined. P<0.05 was considered to indicate a statistically significant difference.

The TFs and DEMs were selected to construct TF-miRNA regulatory networks using the miRWalk2.0 database (24). The TF-miRNAs were identified as hub interactions, as predicted by the miRWalk2.0 database, TargetScan6.2 database (25), miRanda database (26) and RNA22 database (27) simultaneously. The miRNA-DEG regulatory networks were constructed across the aforementioned method. Finally, TF-miRNA-gene axes were constructed to identify a transcriptional regulation model and reveal significant networks involved in HP-associated gastric cancer.

The gastric cancer MKN45 cell line was purchased from the American Type Culture Collection. The MNK-450 cells were cultured in RPMI 1640 medium (HyClone; Cytiva) with 10% FBS and 1% penicillin-streptomycin under a humidified atmosphere with 5% CO₂ at 37°C.

DE TFs and DEMs were confirmed using reverse transcription-quantitative PCR (RT-qPCR) in samples from enrolled patients. To confirm the results from the integrated bioinformatics analysis, 20 fresh HP(+) and 10 fresh HP(−) gastric cancer tissues were collected and examined by experienced pathologists at the Department of Gastrointestinal Surgery, Union Hospital of Huazhong University of Science and Technology between September 1 and December 30, 2019. The tissue samples were frozen and stored immediately in liquid nitrogen following surgical resection. Total cellular RNA was extracted from the gastric samples using RNAiso Plus (Takara Bio, Inc.). For miRNA quantification, synthesis of cDNA was performed using the PrimeScript™ RT reagent kit (GeneCopoeia, Inc.) according to the manufacturer's protocols. The purity and concentration of RNA was measured using a Nanodrop2000 (Thermo Fisher Scientific, Inc.). qPCR was performed using the Real-time PCR Detection System (SLAN; Shanghai Hongshi Medical Technology, Co., Ltd.) and the SYBR Premix Ex Taq II kit (Takara Bio, Inc.). The conditions of PCR cycling were as following: Activation of TaqMan at 95°C for 10 min, and then 40 cycles of denaturation at 95°C for 10 sec, and annealing/extension at 60°C for 60 sec. Relative miRNA expression levels were normalized to U6, while mRNA was normalized to GAPDH, and quantification was performed using the 2^−ΔΔCq method (28). The primer sequences are shown in Table I. R language (https://www.r-project.org/) was used to analyze the RT-qPCR data. All reactions were performed in triplicate.

Tissue lysates were extracted using the RIPA lysis buffer (Beyotime Institute of Biotechnology) supplemented with 1% phenylmethylsulfonyl fluoride (Sigma-Aldrich; Merck KGaA) according to the manufacturer's instructions. The protein concentration was measured using the Pierce bicinchoninic acid protein assay kit (Thermo Fisher Scientific, Inc.). The same amount of protein (30 µg) from all patients with gastric cancer was separated on 10% gels using SDS-PAGE (Bio-Rad Laboratories, Inc.) and transferred onto a PVDF membrane (EMD Millipore) After that, the membrane was blocked for 1 h with 5% non-fat milk and 1X TBST (0.05%) at room temperature. Next, the primary antibodies were used to incubate the membrane at 4°C overnight followed by incubation with secondary antibodies for 2 h at room temperature. β-actin was used as a loading control. The following primary antibodies were used: Rabbit anti-ZEB1 (1:1,000; cat. no. 021544-1-AP; ProteinTech Group, Inc.), rabbit anti-LHX3 (1:1,000; cat. no. 20745-1-AP, ProteinTech Group, Inc.), rabbit anti-PAX2 (1:2,000; cat. no. ab79389; Abcam) and mouse anti-β-actin (1:5,000; cat. no. 60008-1-Ig; ProteinTech Group, Inc.). Anti-mouse IgG (1:5,000; cat. no. ab208001; Abcam) and anti-rabbit IgG (1:5,000; cat. no. ab207999; Abcam) were used as the secondary antibodies. Proteins bands were visualized using electrochemiluminescence (ECL Western Blotting Substrate; cat. no. 32106; Pierce; Thermo Fisher Scientific, Inc.) and analyzed using the ChemiDoc™ XRS Molecular Imager 3.0 system (Bio-Rad Laboratories, Inc.).

ZEB1-small interfering (si)RNA (siZEB1), PAX2-siRNA (siPAX2) and corresponding negative control (NC-siRNA) (Table II) were purchased from Guangzhou RiboBio Co., Ltd. ZEB1 and PAX2 pcDNA3.1 plasmid (ZEB1 and PAX2 overexpression plasmid) and the control vector were purchased from Shanghai GeneChem Co., Ltd. All siRNAs were transfected at a final concentration of 50 nmol/l, and the plasmid was transfected at final concentration using 1.6 µg for a 12-well plate. Lipofectamine^® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) was used for cell transfection at 37°C for 48 h according to the manufacturer's instructions and then subsequent experiments were immediately performed. The promoter reporter construct (ΔE-box) was constructed with the QuikChange site-directed mutagenesis kit (Stratagene; Agilent Technologies, Inc.).

The JASPAR database (http://jaspar.genereg.net) was used to predict the binding region of the miRNA-144 according to the manufacturer's protocols (29). To determine whether ZEB1 targeted the promoter region of the miR-144 gene, using Lipofectamine 3000, the gastric adenocarcinoma MKN-45 cell line was transfected with the recombinant pGL4.1 plasmid (Promega Corporation), including the wild-type or mutant miR-144 gene promoter, empty vector, ZEB1 overexpression plasmid and siZEB1 and NC (siNC). The empty vector and siNC were used as the negative controls for the luciferase reporter gene assay. After 48 h, firefly and Renilla luciferase activities were detected using the dual-luciferase reporter assay system (Promega Corporation). Renilla luciferase activity was used for normalization.

Statistical analyses were performed applying R software (4.0.0; http://www.r-project.org/) unless otherwise noted. All results are presented as the mean ± SD and were assessed via Graphpad Prism 7 (Graphpad Software, Inc.). The significance of differences between two groups was evaluated via Student's unpaired t-test, and the comparison among multiple groups was conducted via one-way ANOVA followed by the LSD post hoc test. All experiments were performed 3 times. P<0.05 was used to indicate a statistically significant difference.

There were 24,991 gene and 1,881 miRNA expression profiles from 141 gastric cancer samples obtained from TCGA mRNA and miRNA datasets, respectively. A total of 1,050 genes and 13 miRNAs were found to be DE in the HP-associated gastric cancer samples, as shown in Fig. 1. The top 20 DEGs are shown in Fig. 2A, while the 13 DEMs are listed in Fig. 2B.

All of the 1,050 DEGs identified, which were based on the normalized mRNA expression data, were screened for gene set enrichment analysis. A total of 180 significant enrichments were identified using the GO analysis, including in the three main categories of biological process, cell component and molecular function. The top three GO terms were ‘T cell activation’, ‘leukocyte proliferation’ and ‘leukocyte migration’. The top 12 GO terms according to the P-values and the top 5 GO enrichments with the gene linkages are displayed in Fig. 3A and B, respectively. Pathway analysis using the KEGG database found that these genes were significantly enriched in ‘Staphylococcus aureus infection’, ‘systemic lupus erythematosus’ and ‘hematopoietic cell lineage’. The top 12 KEGG pathways and the top 5 pathways with the gene linkages are displayed in Fig. 3C and D, respectively.

The top 400 DEGs were imported into the STRING database to identify the interactive association between the proteins. Only the genes with a combined score >0.4 were selected to create the network. Finally, 1,434 pairs of protein associations were identified following removal of the unmatched genes. In addition, genes with >60 edges were considered as hub genes. A total of 6 significant genes were identified based the annotated information from the STRING database: TYROBP, CCR5, C3AR1, LCP2, CCL5 and SPI1 (Fig. 3E).

A total of 252 nodes and 1,434 edges were analyzed using the plug-in, MCODE. The top 2 significant modules were selected for the gene functional annotation involved in the modules. Enrichment analysis showed that the genes in those 2 modules were primarily associated with phytochelatin metabolic process and positive regulation of apoptosis (Fig. 3F).

Using the mRNA and miRNA expression data from TCGA database, 1,050 DEGs and 13 DEMs were identified. From these, 31 DE TFs were identified from the Trrust database (Fig. 4A). Next, the regulatory associations were constructed using the 31 TFs and 13 miRNAs and the miRWalk database, and finally, 7 TF-miRNA pairs (including 7 genes and 4 miRNAs) were constructed (Fig. 4B). Subsequently, 107 pairs of miRNA-DEG interactions were also constructed using the 4 miRNAs and 1,050 DEGs, and the miRWalk database (Fig. 4C). Combining the aforementioned regulatory associations, 181 TF-miRNA-DEG regulatory associations were finally identified (Fig. 4D) (Table III).

RT-qPCR was used to investigate the expression level of the 7 TF-miRNAs (Table III). It was revealed that the expression level of 3 genes and 2 miRNAs were dysregulated [Figs. 5 and 6; HP(+) is represented by the case group and HP(−) is represented by the control group]. Among these genes and miRNAs, the expression levels of PAX2, LHX3, miRNA-144-3p, miRNA-592 were downregulated, while that of ZEB1 was upregulated, which was consistent with the screening results. Furthermore, western blot analysis was used to determine the protein expression levels of the associated proteins in the normal (HP(−) tissues, as well as the tumor groups, and the clinicopathological information of the patients with gastric cancer is shown in Table SI. LHX3 was found to be downregulated and ZEB1 was upregulated, while there was no difference in PAX2 expression levels (Fig. 6H). To verify whether ZEB1 binds to the promoter region of miR-144-3p, to suppress transcriptional activity, the JASPAR database was used to predict the binding region of the miRNA-144, which was subsequently amplified and cloned into the luciferase reporter vector (Fig. 6I). As shown in Fig. 6J, ZEB1 knockdown enhanced the luciferase activity, while the overexpression of ZEB1 markedly inhibited transcriptional activity in the presence of the wild-type E-box construct only, which confirmed that ZEB1 directly targeted the E-box of the miR-144-3p gene promoter to trigger its expression. In order to further investigate whether there was significant regulatory association between ZEB1 and miR-144-3p, an siRNA for ZEB1 (siZEB1) and an overexpression plasmid (ZEB1-OE) for ZEB1 were designed and applied to attenuate or increase the expression of ZEB1 (Fig. 6J). Fig. 6K shows that miR-144-3p expression was significantly decreased in ZEB1 overexpression cells, while being significantly increased in ZEB1 knockdown cells, which further confirmed the regulatory relationship between ZEB1 and miR-144-3p. Additionally, Fig. 6L shows that luciferase activity in the WT miR-144-3p promoter was upregulated under ZEB1 deficiency, while the opposite result was observed when ZEB1 was overexpressed. No significant difference was observed when the binding site was mutated. After silencing and overexpressing PAX2 (Fig. 6M), it was found that there was no significant regulatory association between PAX2 and miR-592 (Fig. 6N). In addition, the results of the luciferase reporter assay showed that there was no significant transcriptional regulation between PAX2 and miR-592 (Fig. 6O).