N6-methyladenosine methylation regulates the tumor microenvironment of Epstein-Barr virus-associated gastric cancer

Abstract

BACKGROUND

N6-methyladenosine (m6A) methylation modification exists in Epstein-Barr virus (EBV) primary infection, latency, and lytic reactivation. It also modifies EBV latent genes and lytic genes. EBV-associated gastric cancer (EBVaGC) is a distinctive molecular subtype of GC. We hypothesized EBV and m6A methylation regulators interact with each other in EBVaGC to differentiate it from other types of GC.

AIM

To investigate the mechanisms of m6A methylation regulators in EBVaGC to determine the differentiating factors from other types of GC.

METHODS

First, The Cancer Gene Atlas and Gene Expression Omnibus databases were used to analyze the expression pattern of m6A methylation regulators between EBVaGC and EBV-negative GC (EBVnGC). Second, we identified Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment of m6A-related differentially expressed genes. We quantified the relative abundance of immune cells and inflammatory factors in the tumor microenvironment (TME). Finally, cell counting kit-8 cell proliferation test, transwell test, and flow cytometry were used to verify the effect of insulin-like growth factor binding protein 1 (IGFBP1) in EBVaGC cell lines.

RESULTS

m6A methylation regulators were involved in the occurrence and development of EBVaGC. Compared with EBVnGC, the expression levels of m6A methylation regulators Wilms tumor 1-associated protein, RNA binding motif protein 15B, CBL proto-oncogene like 1, leucine rich pentatricopeptide repeat containing, heterogeneous nuclear ribonucleoprotein A2B1, IGFBP1, and insulin-like growth factor 2 binding protein 1 were significantly downregulated in EBVaGC (P < 0.05). The overall survival rate of EBVaGC patients with a lower expression level of IGFBP1 was significantly higher (P = 0.046). GO and KEGG functional enrichment analyses showed that the immunity pathways were significantly activated and rich in immune cell infiltration in EBVaGC. Compared with EBVnGC, the infiltration of activated CD4+ T cells, activated CD8+ T cells, monocytes, activated dendritic cells, and plasmacytoid dendritic cells were significantly upregulated in EBVaGC (P < 0.001). In EBVaGC, the expression level of proinflammatory factors interleukin (IL)-17, IL-21, and interferon-γ and immunosuppressive factor IL-10 were significantly increased (P < 0.05). In vitro experiments demonstrated that the expression level of IGFBP1 was significantly lower in an EBVaGC cell line (SNU719) than in an EBVnGC cell line (AGS) (P < 0.05). IGFBP1 overexpression significantly attenuated proliferation and migration and promoted the apoptosis levels in SNU719. Interfering IGFBP1 significantly promoted proliferation and migration and attenuated the apoptosis levels in AGS.

CONCLUSION

m6A regulators could remodel the TME of EBVaGC, which is classified as an immune-inflamed phenotype and referred to as a “hot” tumor. Among these regulators, we demonstrated that IGFBP1 affected proliferation, migration, and apoptosis.

Key Words: N6-methyladenosine methylation, Tumor microenvironment, Epstein-Barr virus, Gastric cancer, Insulin-like growth factor binding protein 1

Core Tip: In this work, we found that N6-methyladenosine methylation regulators activated the expression of immune cells and inflammatory factors in Epstein-Barr virus-associated gastric cancer (EBVaGC), which could remodel the tumor microenvironment. EBVaGC is classified as an immune-inflamed phenotype and referred to as a “hot” tumor. Insulin-like growth factor binding protein 1 affected the proliferation, migration, and apoptosis of EBVaGC. Thus, our study provided a theoretical basis for future research on the pathogenesis and therapeutic targets of EBVaGC.

INTRODUCTION

Since the 1960s, RNA modification has been fully proven to be a key factor in regulating RNA metabolism in various biological processes[1]. So far, more than 171 known RNA post-transcriptional modifications are known, and N6-methyladenosine (m6A) methylation modification is one of the most abundant and widespread modifications in most eukaryotes[2]. It has been demonstrated that m6A methylation modification and its related enzymes play a pivotal role in different phases of RNA metabolism, including mRNA translation, degradation, splicing, export, and folding[3,4].

In mammalian cells, m6A methylation modification is a dynamic and reversible process, which is regulated by different types of regulatory proteins called the methyltransferases [writers: Methytransferase-like 3/14 (METTL3/14); Wilms tumor 1-associated protein (WTAP); Vir-like m6A methyltransferase associated (VIRMA/KIAA1429); etc], the demethylases [erasers: Fat mass and obesity-associated protein (FTO); AlkB homolog 5; etc], and binding proteins [readers: YTH domain-containing family protein 1/2/3 (YTHDF1/2/3); YTH domain-containing 1/2 (YTHDC1/2); etc][5]. m6A methylation modification plays a crucial role in multiple human physiological processes and pathological phenomena, such as obesity, diabetes, immune dysregulation, neuronal diseases, tumors, etc[6,7].

It has been found that m6A methylation modification not only exists in eukaryotes but also in viruses, such as human immunodeficiency virus, Zika virus, influenza A virus, hepatitis B virus, hepatitis C virus, Kaposi’s sarcoma-associated herpesvirus, Epstein-Barr virus (EBV), etc[8,9]. Crucially, the m6A methylation modification is reported to regulate viral replication and infectivity. EBV, belonging to the γ-Herpesvirus subfamily, was first discovered in Burkitt lymphoma patients in 1964[10]. EBV is closely related to the initiation and progression of various malignancies, such as Hodgkin’s lymphoma, Burkitt’s lymphoma, natural killer/T cell lymphoma, nasopharyngeal carcinoma, and gastric cancer (GC)[11].

In recent years, some researchers have noticed that the EBV lifecycle is regulated by m6A methylation modification in EBV-related tumors. Previous studies found that m6A methylation modification not only existed in EBV primary infection, latency, or lytic reactivation but also modified EBV latent genes and lytic genes[12,13]. m6A regulators YTHDF1, YTHDF2, and METTL3 further suppress EBV infection and lytic replication by repressing the expression of EBV lytic genes BZLF1 and BRLF1[12,14,15]. YTHDF1 suppresses EBV infection and replication by recognizing m6A modification sites in EBV transcripts and promoting m6A-dependent RNA decay[14]. Zhang et al[16] found that deletion of YTHDF2, METTL3, METTL14, VIRMA/KIAA1429, and WTAP promoted EBV lytic protein expression. EBV-encoded latent protein EBNA3C can activate the transcription of METTL14 and interact with METTL14 to promote its stability and tumorigenesis[12]. Based on the above research results, we partially understand that m6A methylation modification can further regulate EBV-mediated tumorigenesis and development by regulating the lifecycle of EBV.

EBV-associated GC (EBVaGC) is a distinctive molecular subtype of GC that accounts for approximately 10% of cases with specific clinicopathologic and molecular features[17]. At present, some studies have found that m6A methylation plays an important role in the occurrence and development of GC, such as regulating the proliferation, invasion, apoptosis, and clinical prognosis[18,19]. The expression of most m6A methylation regulators in GC tissues is significantly upregulated[19,20]. When GC patients were further divided into EBVaGC and EBV-negative GC (EBVnGC) according to the EBV infection status, the m6A methylation regulators insulin-like growth factor 2 binding protein 1 (IGF2BP1), VIRMA/KIAA1429, leucine rich pentatricopeptide repeat containing (LRPPRC), YTHDF3, and ZC3H13 were significantly downregulated in EBVaGC patients compared to EBVnGC patients. They also observed that FTO was significantly upregulated[20]. Therefore, we hypothesized that EBV and m6A methylation regulators interact with each other in EBVaGC, which differentiates it from other types of GC. At present, there are few studies on the relationship between m6A methylation regulators and EBVaGC.

In order to further understand the relationship between m6A methylation regulators and EBVaGC, we explored the expression pattern of m6A methylation regulators in EBVaGC using The Cancer Gene Atlas (TCGA) and the Gene Expression Omnibus (GEO) database and provided guidance value for further research on the mechanism of EBV in the initiation and progression of GC.

MATERIALS AND METHODS

Data source and processing

Public gene expression data and full clinical annotation were searched in GEO and TCGA databases. The GSE51575 dataset, including 12 EBVaGC and 14 EBVnGC patients, was obtained from the GEO database. TCGA-STAD dataset, including 22 EBVaGC and 224 EBVnGC patients, was obtained from the TCGA database.

Selection of m6A methylation regulators

We extracted a total of 24 m6A methylation regulators to identify significant differences in m6A regulators between EBVaGC and EBVnGC, including 7 writers [METTL3, METTL14, WTAP, ZC3H13, RBM15, RNA binding motif protein 15B (RBM15B), and CBL proto-oncogene like 1 (CBLL1)], 15 readers [YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPR, heterogeneous nuclear ribonucleoprotein A2B1 (HNRNPA2B1), insulin-like growth factor binding protein 1 (IGFBP1), IGFBP2, IGFBP3, RBMX, ELAVL1, and IGF2BP1], and 2 erasers (FTO and ALKBH5). We used the “limma” package to screen significant m6A regulators between EBVaGC and EBVnGC. A P < 0.05 and log₂|fold change|> 1 were used as screening criteria. The results were visualized with a heat map and histogram.

Construction of a random forest model and support vector machine

The random forest model (RF) and support vector machine (SVM) model were constructed as a training model to predict the occurrence of EBVaGC. The model was evaluated using the “reverse cumulative distribution of residuals”, “boxplots of residual”, and receiver operating characteristic curve. We established an RF model to select candidate m6A regulators to predict the occurrence of EBVaGC among the seven significant m6A regulators.

Construction of the nomogram model

We constructed a nomogram model to predict the occurrence of EBVaGC based on the selected candidate m6A regulators using the “rms” package. The calibration curve was used to assess the consistency of our predicted values with reality. Decision curve analysis and clinical impact curve were performed to evaluate whether the model was beneficial to the screening of EBVaGC.

Functional enrichment analysis

We used the “limma” package in R to screen for m6A-related differentially expressed genes (DEGs) between distinct GC. P < 0.05 and log₂|fold change| > 1 after adjustment were selected as the screening criterion. Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed in m6A-related DEGs involved in EBVaGC using the “clusterprofiler” package. The results were visualized with an enrichment circle diagram.

Immune cell infiltration

We used the single sample gene enrichment analysis to quantify the relative abundance of immune cells. The gene set for marking each infiltrating immune cell type was obtained from another study (Charoentong, unreferenced). The enrichment scores generated by single sample gene enrichment analysis were utilized to represent the relative abundance of immune cells in each sample, and the results were visualized with a histogram.

Inflammatory factors expression

EBV-infected cells could induce local secretion of proinflammatory, immunosuppressive cytokines, and chemokines. After consulting the literature, we identified 10 inflammatory factors related to EBVaGC, including five proinflammatory factors [interleukin (IL)-1β, IL-17, IL-21, IL-22, and interferon (IFN)-γ] and five immunosuppressive factors (IL-4, IL-6, IL-8, IL-10, and IL-13). We analyzed the expression level of inflammatory factors in EBVaGC. P < 0.05 was selected as the screening criteria. The results were visualized with histogram.

Cell lines and culture conditions

EBVaGC cell lines (SNU719) and EBVnGC cell lines (AGS) were used in our study. SNU719 was purchased from Baiyi Biotechnology Center (Shanghai, China). AGS was purchased from Pricella (Wuhan, China). All cell lines were maintained in Dulbecco’s modified Eagle’s medium (Gibco, Thermo Fisher Scientific, Waltham, MA, United States) containing 10% fetal bovine serum (Gibco) at 37 °C with 5% CO₂.

Plasmid construction and RNA interference

The pcDNA3.1-IGFBP1 plasmids were produced by Jiman Biology (Nanjing, China). The small interfering RNA (siRNA) target sequences for IGFBP1 mRNA were designed and synthesized by Ruibo Bio (Jiangsu, China). Plasmid transfection and RNA interference were performed using Lipofectamine 2000 transfection reagent (Thermo Fisher Scientific) according to the manufacturer’s protocols. A total of 5 μL of plasmid DNA per six-well plate was used for transfection. The cells were harvested after 48 h for downstream analyses.

Western blotting

The cells were extracted using RIPA lysis buffer containing phosphatase inhibitors and phenylmethylsulfonyl fluoride and added to 1 × loading buffer and heated for 10 min at 100 °C. Proteins were separated by sodium-dodecyl sulfate gel electrophoresis electrophoresis at 70 V for 30 min and then increased to 120 V to completion. The proteins were transferred to polyvinylidene fluoride films at 220 mA for 1 h in an ice-water mixture. The membrane was blocked for 2 h in 5% non-fat dry milk and then incubated with the following primary antibodies: Rabbit monoclonal anti-IGFBP1 (1:1000; ab181141, Abcam, Cambridge, United Kingdom) and rabbit monoclonal anti-β-tubulin (1:5000; A12289, Abclonal, Woburn, MA, United States) antibodies. Then, the membranes were incubated with HRP-conjugated secondary antibody (1:10000; AS111, Abclonal). Proteins were visualized using an Immobilon TM Western Chemiluminescent HRP substrate (Millipore, Burlington, MA, United States).

Cell counting kit-8 assay

Cell proliferation was detected using the cell counting kit-8 (Qidongzi, China) according to the manufacturer’s protocol. Three hours after the siRNAs or plasmids were transfected into cells, absorbance was tested using a Soft-Max apparatus (ELx808; BioTek, Winooski, VT, United States) at a wavelength of 450 nm.

Transwell assay

The migration ability of the cells was measured using transwell plates. Cells were plated in the upper chamber, and 0.6 mL Dulbecco’s modified Eagle’s medium containing 10% fetal bovine serum was added to the lower chamber. After 24 h, cells that migrated through the membrane to the lower surface were fixed with paraformaldehyde and stained with crystal violet. Then, photos were taken through the microscope, and the numbers were counted.

Apoptosis analysis

Apoptotic cells were measured by a FITC Annexin V Apoptosis Detection Kit I (BD Pharmingen, San Diego, CA, United States) according to the manufacturer’s protocol.

Statistical analysis

R software (version 4.1.2) was used for statistical analysis. All cell experiments were repeated at least three times. GraphPad Prism 9.3 software was used for statistical analysis. All parametric analyses were based on two-tailed tests, and statistical significance was set at P < 0.05.

RESULTS

Landscape of the 24 m6A regulators in EBVaGC

m6A methylation was one of the most common RNA modifications involved in the progression of numerous cancers. To explore the expression pattern of m6A methylation regulators in EBVaGC, we first extracted and analyzed the expression levels of m6A regulators from the GSE51575 database using the “limma” package. As shown in Figure 1A and B, the expression levels of WATP, RBM15B, CBLL1, LRPPRC, HNRNPA2B1, IGFBP1, and IGF2BP1 were significantly downregulated in EBVaGC tissues compared with EBVnGC (P < 0.05). In addition, the 24 m6A regulatory factors could interact with each other during the occurrence and development of GC, showing a positive or negative correlation (Figure 1C).

Figure 1 Landscape of the 24 N6-methyladenosine regulators in Epstein-Barr virus-associated gastric cancer. A: Expression heat map of the 24 N6-methyladenosine (m6A) regulators in Epstein-Barr virus (EBV)-associated gastric cancer (EBVaGC) and EBV-negative gastric cancer (EBVnGC); B: Differential expression histogram of the 24 m6A regulators identified between EBVaGC and EBVnGC; C: The interaction between m6A regulators. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001. EBV: Epstein-Barr virus.

Then, we explored the prognostic role of m6A regulation in EBVaGC. We evaluated the correlation of m6A regulators with corresponding clinical follow-up information through Kaplan-Meier analysis from the TCGA-STAD dataset. As shown in Figure 2, the overall survival rate of EBVaGC patients with low IGFBP1 expression was significantly higher (P = 0.046), while the expression level of other regulators were not correlated with the overall survival rate of EBVaGC patients (P > 0.05). The above results suggested that EBV infection status in GC might lead to differences in the expression levels of m6A regulators, thus affecting the survival prognosis of GC patients.

Figure 2 Kaplan-meter survival curves of 7 N6-methyladenosine regulators with differences in Epstein-Barr virus-associated gastric cancer. A: Heterogeneous nuclear ribonucleoprotein A2B1; B: Insulin like growth factor 2 binding protein 1; C: Insulin like growth factor binding protein 1; D: Leucine rich pentatricopeptide repeat containing; E: RNA Binding Motif Protein 15B; F: CBL proto-oncogene like 1; G: Wilms tumor 1 associated protein. HNRNPA2B1: Heterogeneous nuclear ribonucleoprotein A2B1; IGF2BP1: Insulin like growth factor 2 binding protein 1; IGFBP1: Insulin like growth factor binding protein 1; LRPPRC: Leucine rich pentatricopeptide repeat containing; RBM15B: RNA Binding Motif Protein 15B; CBLL1: CBL proto-oncogene like 1; WTAP: Wilms tumor 1 associated protein.

Construction of the RF model and SVM model

To forecast the incidence of EBVaGC, we developed the RF and SVM model to screen prospective m6A regulators from the seven significant m6A regulators. Both “reverse cumulation distribution of residual” (Figure 3A) and “boxplots of residual” (Figure 3B) demonstrated that the RF model had minimum residuals, making it the best model for predicting the incidence of EBVaGC. The area under the curve value of the receiver operating characteristic curve indicated that the RF model was more accurate than the SVM model. We sorted and visualized the DEGs of EBVaGC according to the Gini index of genes (Figure 3C). The Gini index of IGF2BP1 and IGFBP1 was significantly greater than that of other m6A regulators, indicating that IGF2BP1 and IGFBP1 play a key role in the occurrence and progression of EBVaGC (Figure 3D).

Figure 3 Random forest model and support vector machine model. A: Reverse cumulative distribution of residual was plotted to show the residual distribution of random forest (RF) and support vector machine (SVM) model; B: Boxplots of residual was plotted to show the residual distribution of RF and SVM model; C: Receiver operator characteristic curves indicated the accuracy of the RF and SVM model; D: Gini index of the 7 N6-methyladenosine regulators with differences in Epstein-Barr virus-associated gastric cancer based on the RF model. RF: Random forest; SVM: Support vector machine; HNRNPA2B1: Heterogeneous nuclear ribonucleoprotein A2B1; IGF2BP1: Insulin like growth factor 2 binding protein 1; IGFBP1: Insulin like growth factor binding protein 1; LRPPRC: Leucine rich pentatricopeptide repeat containing; RBM15B: RNA Binding Motif Protein 15B; CBLL1: CBL proto-oncogene like 1; WTAP: Wilms tumor 1 associated protein.

Construction of the nomogram model

We constructed a nomogram model based on the two candidate m6A regulators IGF2BP1 and IGFBP1 by using the “rms” package to predict the occurrence of EBVaGC (Figure 4A). The calibration curve showed that the nomogram model predictivity was accurate (Figure 4B). The red line in the decision curve analysis remained above the gray and black lines from 0 to 0.8, indicating that the decision of the nomogram model may be beneficial in the screening of EBVaGC (Figure 4C). The clinical impact curve revealed that the predictive power of the nomogram model was remarkable (Figure 4D).

Figure 4 Construction of the nomogram model. A: Construction of the nomogram model based on the insulin like growth factor binding protein 1 and insulin like growth factor 2 binding protein 1; B: Predictive calues of the model by the calibration curve; C: Decision curve analysis to evaluate the model; D: Clinical impact of the model. IGF2BP1: Insulin like growth factor 2 binding protein 1; IGFBP1: Insulin like growth factor binding protein 1; m6A: N6-methyladenosine.

Enrichment analysis of genes participating in EBVaGC

A total of 356 m6A-related DEGs were selected between EBVaGC and EBVnGC (Figure 5A). To further explore the potential mechanism of these DEGs in EBVaGC, GO and KEGG functional enrichment analyses were applied. The GO functional enrichment analyses showed that “biological processes” mainly enriched in T cell activation, leukocyte cell-cell adhesion, and regulation of T cell activation. “Cellular components” were mainly enriched in external side of plasma membrane, endocytic vesicle membrane, and MHC protein complex. “Molecular function” was mainly enriched in receptor ligand activity, signaling receptor activator activity, and immune receptor activity (Figure 5B). In addition, the KEGG enrichment analysis was mainly involved in cell adhesion molecules, hematopoietic cell lineage, and cytokine-cytokine receptor interaction (Figure 5C). The results showed that enrichment of biological processes remarkably related to m6A modification and immunity, which confirmed that m6a modification played a positive role in immune regulation.

Figure 5 Enrichment analysis of genes participating in Epstein-Barr virus-associated gastric cancer. A: N6-methyladenosine (m6A)-related differentially expressed genes (DEGs); B: Gene ontology analysis that explores the potential mechanism underlying the effect of the 365 m6A-related DEGs on the occurrence and development of Epstein-Barr virus-associated gastric cancer; C: Kyoto Encyclopedia of Genes and Genomes analysis that explores the potential signal pathway. FC: Fold change.

Comparison of the immune activity between subgroups

EBV could induce a highly variable composition of tumor microenvironment (TME), which consists of both innate and adaptive immune cells. The TME on which tumor cells depend for growth and survival also plays an essential component in the occurrence and development of EBV-related malignancies[21]. Therefore, we calculated the abundance of immune cell infiltration and inflammatory factors in two subgroups. We found that the EBVaGC group was remarkably rich in innate and adaptive immune cell infiltration, including B cells, activated CD4+ T cells, activated CD8+ T cells, regulatory T cells, natural killer cells, macrophages, eosinophils, mast cells, myeloid-derived suppressor cells, and plasmacytoid dendritic cells (P < 0.05). The expression of activated CD4+ T cells, activated CD8+ T cells, monocytes, and dendritic cells were significantly higher in EBVaGC (Figure 6A). Compared with EBVnGC, the expression level of proinflammatory factors IL-17, IL-21, and IFN-γ and immunosuppressive factor IL-10 were significantly increased in EBVaGC (P < 0.05) (Figure 6B). The results showed that EBVaGC was classified as an immune-inflamed phenotype, characterized by adaptive immune cell infiltration and immune activation, and m6A modification played a non-negligible regulation role in shaping different TME landscapes.

Figure 6 Tumor microenvironment cell infiltration characteristics in Epstein-Barr virus-associated gastric cancer. A: Single sample gene set enrichment analysis; B: the expression of inflammatory factors between Epstein-Barr virus-associated gastric cancer and Epstein-Barr virus-negative gastric cancer. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001. EBV: Epstein-Barr virus; IL: Interleukin; IFN: Interferon; MDSC: Myeloid-derived suppressor cell.

IGFBP1 expression verification

To verify the expression of IGFBP1 in EBVaGC cell lines (SNU719), we used western blotting to detect the protein expression level of IGFBP1 in EBVaGC cell lines (SNU719) and EBVnGC cell lines (AGS). According to western blotting, the expression of the IGFBP1 protein in the EBVaGC cell line (SNU719) was significantly lower than in the EBVnGC cell line (AGS) (P < 0.05) (Figure 7A).

Figure 7 Role of insulin like growth factor binding protein 1 in Epstein-Barr virus-associated gastric cancer and Epstein-Barr virus-negative gastric cancer cell lines. A: the expression of insulin like growth factor binding protein 1 (IGFBP1) in Epstein-Barr virus-associated gastric cancer (EBVaG) and Epstein-Barr virus-negative gastric cancer (EBVnGC) cell lines; B: Role of IGFBP1 in proliferation in EBVaG and EBVnGC cell lines; C: Role of IGFBP1 in migration in EBVaG and EBVnGC cell lines; D: Role of IGFBP1 in apoptosis in EBVaG and EBVnGC cell lines. ^aP < 0.05, ^bP < 0.01, ^cP < 0.001. IGFBP1: Insulin like growth factor binding protein 1; NC: Normal control; PI: Propidium iodide.

Next, we investigated the potential effects of IGFBP1 on cell proliferation, migration, and apoptosis in the EBVaGC cell line (SNU719). We transfected PCDNA3.1-IGFBP1 and PCDNA3.1-NC plasmids into the EBVaGC cell line (SNU719). The results showed that the overexpression of IGFBP1 could inhibit proliferation of the EBVaGC cell line (SNU719) (Figure 7B). The transwell results showed that the overexpression of IGFBP1 could inhibit the migration of the EBVaGC cell line (SNU719) (Figure 7C). The Annexin V staining assay showed that the overexpression of IGFBP1 induced a significant increase in the percentage of late apoptosis cells compared with PCDNA3.1-NC (Figure 7D). Meanwhile, we also used siRNA in the EBVnGC cell line (AGS). The results showed that the proliferation and migration of AGS cells were significantly increased (Figure 7B and C), and the apoptosis of AGS cells was significantly decreased compared with the siRNA-NC (Figure 7D).

DISCUSSION

m6A methylation modification is the most common modification in human mRNA[2], which is considered a new layer of epigenetic regulation on mRNA processing, translation, and stability[3]. Mounting evidence has shown that dysregulation of m6A modification is closely related to various physiological and pathological phenomena in humans, including obesity, immune dysregulation, carcinogenesis, etc[7]. m6A modification could promote or restrict herpesvirus infection and replication by regulating the life cycle of many viruses[9,16,22].

Recent studies have shown that both EBV latent and lytic transcripts are modified by m6A during initial infection, the incubation period, or lytic reactivation[12,14]. EBV uses m6A methylation to promote cell survival and progeny virus production during the lytic cycle. The depletion of YTHDF1 and YTHDF2 has been shown to promote EBV lytic protein expression[12,13]. YTHDF2 modifies viral RNA by binding m6A methylation, which could restrict EBV replication by promoting RNA decay[12]. The recent findings demonstrated that YTHDF 1-3 proteins control the viral RNA decay and replication through PIASI-mediated SUMOylation (E3 small ubiquitin-like modifier), which is known as an EBV restriction factor[23].

EBV immediate-early protein BZLF1 interacts with the promoter of m6A methyltransferase METTL3, inhibiting its expression[15]. Knockout of METTL3 decreased the expression of viral lytic proteins and reduced progeny virion production[24]. Dai et al[15] provided solid evidence that m6A modification is involved in the crosstalk between EBV and host cells via a positive feedback loop to promote EBV infection. Lang et al[12] revealed that EBV-encoded latent protein EBNA3C upregulated the expression of METTL14 and directly interacts with and stabilizes MRTTL14 to promote the growth and proliferation of tumor cells. Lang et al[12] found that the expression of METTL14 was significantly upregulated in EBV-positive post-transplant lymphoproliferative disorders, and it promoted the proliferation and colony formation of EBV-positive cells. In addition, in our study, we also found that m6A methylation regulators drove the occurrence and development of EBVaGC.

More and more studies have shown that abnormal expression of many molecules are related to m6A methylation modification, which affects the occurrence and progression of GC[25,26]. Most m6A methylation regulators in GC tissues were significantly upregulated compared with normal gastric tissues, and high expression predicted a poor prognosis for GC patients[26,27].

In mammals, m6A methylation regulators interacted with EBV to drive and promote the occurrence and development of tumors. Therefore, it was necessary to further study the impact of m6A methylation modifications on EBVaGC and provide a theoretical basis for further exploring the pathogenesis of EBVaGC. Previous studies found that the m6A total RNA methylation level of the EBVaGC cell line was significantly lower than that of the EBVnGC cell line[28]. Xiao et al[28] found that EBV may activate the nuclear factor-kappa B signaling pathway to inhibit the expression of WTAP to achieve the regulation of m6A methylation. Zhang et al[20] showed that the expression of m6A methylation regulators IGFBP1, VIRMA/KIAA1429, LRPPRC, YTHDF3, and ZC3H13 were remarkably downregulated in EBVaGC patients, while FTO was significantly upregulated. In our study, we also found that the expression of WTAP, RBM15B, CBLL1, LRPPRC, HNRNPA2B1, IGFBP1, and IGF2BP1 was significantly downregulated in EBV-related GC compared with EBVnGC (P < 0.05). There were some differences between our research results and those of the study by Zhang et al[20]. The main reason was that the databases used by them are inconsistent. However, we could find that there were differences in the expression of WTAP, IGF2BP1, and LRPPRC in EBVaGC. In future studies, it will be beneficial to investigate the mechanism of m6A modifications in EBVaGC.

At present, we knew that m6A methylation modifications affect the occurrence and development of EBV-related tumors. In our study, we found that IGFBP1 could effectively predict the occurrence of EBVaGC. IGFBPs belong to the insulin-like growth factor family identified as m6A eraser with a common RNA-binding domain that inhibit the degradation of m6A modified mRNA by recognizing the UGGAC sequence of m6A mRNA and recruiting RNA stabilizers. IGFBP1 played an important role in regulating cell proliferation, differentiation, apoptosis, and carcinogenesis[29-31]. IGFBP1 could inhibit cancer cell invasion and proliferation in hepatocellular carcinoma.

In Helicobacter pylori-associated GC, overexpression IGFBP1 could reduce the promoting effect of matrix metalloproteinase 9 on BGC823 and AGS cell migration. However, previous research found that higher expression of IGFBP1 was associated with hematogenous metastasis, poor prognosis, and recurrence-free survival of GC patients[32]. Serum IGFBP1 had potential clinical value for the early diagnosis of upper gastrointestinal tumors[33]. We also found that higher IGFBP1 of expression had a worse prognosis in EBVaGC patients (P < 0.05). In our study, we found that overexpression of IGFBP1 could inhibit the proliferation and migration of the EBVaGC cell line (SNU719) and promote apoptosis.

Kim et al[34] explored the role of IGFBP1 in liver metastasis from colorectal cancer, and the overexpression of IGFBP1 significantly decreased cell proliferation and invasiveness of SW480 cell lines in the in vitro assays. However, in vivo experiments showed liver metastasis in mice transplanted with IGFBP1-overexpressing SW480 cells, which seemed to contradict with the in vitro results. The study by Kim et al[34] found that β-catenin and c-Myc were continuously expressed in mice with IGFBP1 overexpression with liver metastases. The abnormally activated Wnt/β-catenin signaling axis could promote the renewal and aberrant proliferation of cancer stem cells through downstream molecules such as c-Myc and survivin[35]. FAK is overexpressed and activated in many cancer types. The activation of FAK/p-FAK could promote the invasion and migration of cancer cells[36]. IGFBP1 potentiates cell proliferation, adhesion, and migration through activating the β1/Src/FAK pathway[37]. Based on the above studies, we believe that IGFBP1 may have potential unknown behaviors or alternative pathways that affect the biological function of tumor cells.

EBVaGC represents a distinct subtype, which harbors deregulation of immune response genes, strong evidence of immune infiltration, and a high-level activation of immune checkpoint pathways[38]. Indeed, several studies reported that EBVaGC patients are more sensitive and are beneficiaries of immune checkpoint inhibitors. For example, EBVaGC patients had a good response rate to anti-programmed death protein 1 (PD-1) therapy as well as a longer progression-free survival time after anti-PD-1 therapy[39]. TME plays a critical role in GC progression, immune escape, and immunotherapy[40,41]. Recently, several studies revealed the special correlation between m6A methylation modification and TME infiltrating immune cells[42,43]. In our study, we found that m6A methylation modification is related to some human immune functions. Bose et al[44] also found that m6A plays an important role in overcoming the innate immune response during EBV reactivation. EBV targets the host’s RNA methylation system by inhibiting IFN signaling, weakening the innate immune response and promoting successful lytic replication of the virus[44]. EBV infection could induce innate and adaptive immune responses that inhibit virus replication and bring the infection under control. EBV lytic and latent proteins were rich in antigens that activate natural killer cells as well as specific CD4+ T cells and CD8+ T cells, as anti-EBV immunity. Previous studies reported that high immune cells infiltrating in EBV-positive tumor cells, and the infiltration of CD4+ T cells, CD8+ T cells, dendritic cells, macrophages, and natural killer cells were significantly higher in EBV-positive samples[45,46].

In our study, we found that EBVaGC patients had increased infiltrating immune cells, mainly composed of CD4+ T cells, CD8+ T cells, and dendritic cells. Tumor-infiltrating lymphocytes in the immune microenvironment are associated with the prognosis of GC. CD8+ T cells and natural killer cells predict a better prognosis in EBVaGC[46,47]. EBV could induce local secretion of proinflammatory and immunosuppressive cytokines, which in turn could affect the function of immune cells. In our study, we found that the expression levels of proinflammatory IL-17, IL-21, and IFN-γ and immunosuppressive cytokine IL-10 were significantly increased in EBVaGC. Therefore, we speculated that IGFBP1 might affect the immunotherapeutic effect and prognosis by regulating the TME. These findings suggested that the potential mechanisms for the change of m6A modification patterns are mediated by EBV.

In EBV-related malignant tumors, T cells are an important component of the TME, which are actively recruited by chemokines that are produced by the tumor cells. T cells are an indispensable component of adaptive immunity against EBV. These tumor infiltrating lymphocytes produce an immune-activated TME by eliminating EBV-infected malignant cells. Dendritic cells are professional for antigen presentation that present antigenic peptides to T cells[48]. We found a large amount of immune cell infiltration in EBVaGC. GO and KEGG enrichment analysis showed that the immune activation pathway of EBVaGC was significantly enhanced, including T cell activation, cell adhesion molecules, the interaction between cytokines, cytokine receptors, MHC molecules, receptor ligand activity, signal receptor activator activity, and immune receptor activity. The above results revealed that EBVaGC is a subtype of GC characterized by immune activation.

Zhang et al[20] built a scoring system based on m6A regulatory factors to quantify individuals of different GC patients. The results showed that the lower m6A score is concentrated in the EBVaGC patients, which is classified as an immune-activation phenomenon. These results explained why EBVaGC patients are sensitive to immunotherapy and benefit from immune-checkpoint inhibitor therapy.

In this study, we identified that m6A methylation regulators are involved in the occurrence and development of EBVaGC by using TCGA and GEO databases. However, our study also had some limitations. First, although we had reviewed the literature and determined 24 m6A methylation regulators, a series of newly identified regulators need to be incorporated into the study. Second, the cases of EBVaGC in the GSE51575 and TCGA databases are not rich enough, and further expansion of the sample size is needed for analysis and demonstration in the future. Finally, the potential functions and effects of m6A-associated regulators on the immune microenvironment of EBVaGC should be further explored both in vivo and in vitro.

CONCLUSION

The expression level of m6A methylation regulators were downregulated in EBVaGC, which activated the expression of immune cells and built a “hot” immune microenvironment. m6A shows potential as a target of lytic induction therapy for EBVaGC. IGFBP1 may mediate m6A methylation modification to affect the occurrence and development of EBVaGC.

ACKNOWLEDGEMENTS

We thank the authors who provided the Gene Expression Omnibus and The Cancer Gene Atlas databases.