DNA methylation-regulated YTHDF2 correlates with cell migration and immune cell infiltration in glioma

Analysis of the expression and prognosis of YTHDF2

We employed the TCGA (https://www.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi), CGGA (http://www.cgga.org.cn/) [6], Rembrandt (http://gliovis.bioinfo.cnio.es/) [7] and Gravendeel (http://gliovis.bioinfo.cnio.es/) [8] to examine the expression pattern, clinical significance and prognosis of YTHDF2 in LGG.

GSEA analysis and immune cell infiltration analysis

In this study, we used GSEA software to explore the signaling pathway involved by YTHDF2 in glioma. TIMER (https://cistrome.shinyapps.io/timer/) [9]. TIMER was used to determine the relationship between YTHDF2 expression and immune cell infiltration in glioma.

DNA methylation analysis

MEXPRESS (https://mexpress.be/) is an online tool for DNA methylation database, In our study, we use the SMART (http://www.bioinfo-zs.com/smartapp/), EXPRESS (https://mexpress.be/), and MethSurv (https://biit.cs.ut.ee/methsurv/) databases to examine the relationship between the YTHDF2 expression and CNV, DNA methylation in LGG [10–12].

Drug sensitivity analysis

In this study, we used the GDSC and CTRP databases to explore the correlations between YTHDF2 and various drug sensitivity [13, 14].

Cell culture and real-time PCR

The glioma cell line was purchased from the cell bank of Kunming Institute of Zoology and cultured in DMEM medium (Corning, USA) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. The RT-PCR primers shown in the following, YTHDF2-F: AGCCCCACTTCCTACCAGATG, YTHDF2-R: TGAGAACTGTTATTTCCCCATGC; β-actin-F: CTTCGCGGGCGACGAT, β-actin-R: CCATAGGAATCCTTCTGACC. The expression quantification was obtained with the 2^−ΔΔCt method. The shRNA primer used in this study following: shYTHDF2-1: CCTACTTACCCAGTTACTACA, shYTHDF2-2: GCTCTGGATATAGTAGCAATT.

Cell proliferation, cell migration and IHC assay

Cell proliferation, cell migration and IHC assay was performed as previously described [15]. Cell proliferation and cell migration assay used to determine the function of YTHDF2 on glioma cell migration and proliferation. IHC assay used to explore the expression level of YTHDF2 in glioma tissues.

Expression and prognostic values of YTHDF2 in human cancers

To examine the expression pattern of YTHDF2 in pan-cancer, we used the TIMER database analysis it’s expression pattern, the result demonstrated that YTHDF2 was significantly over-expression in BLCA, CHOL, COAD, ESCA, HNSC, KICH, KIRP, KIRC, LIHC, LUAD, PRAD, SKCM, STAD, THCA and UCEC (Figure 1A). Similar to YTHDF2 mRNA, the protein level of YTHDF2 was highly in BRCA, COAD, OV, KIRC, and UCEC (Figure 1B). Moreover, we employ the CCLE databases to examine YTHDF2 expression in various cancer cells lines. Results suggested that YTHDF2 was elevated in different cancer cell lines (Figure 1C).

Figure 1. The expression of YTHDF2 in human cancer. (A) The expression of YTHDF2 in pan-cancer examine by TIMER database. (B) The protein level of YTHDF2 in diverse cancer examine by the UALCAN database. (C) The expression of YTHDF2 in diverse cancer cell lines examined by CCLE database.

We employed KM plot to explore the prognostic value of YTHDF2 expression in human cancer. Our results confirmed that higher level of YTHDF2 overexpression closely correlates to the poor clinical prognosis of LIHC, SARC, and LGG. On the contrary, up-regulation of YTHDF2 correlated with a better prognosis in BLCA, CESC, KIRC, LUAD, and OV (Figure 2).

Figure 2. The prognosis value of YTHDF2 in human cancer. (A) The prognosis value of YTHDF2 in BLCA, CESC, and KIRC was examined by the km plot database. (B) The prognosis value of YTHDF2 in LIHC, LUAD, and OV was examined by the km plot database. (C) The prognosis value of YTHDF2 in SARC and LGG was examined by the km plot database.

YTHDF2 was up-regulated in LGG

To determine the expression of YTHDF2 in LGG, we adopt the public database to examine the expression of YTHDF2 in glioma tissues and normal group. Results found out that YTHDF2 was overexpressed in the LGG tissues (Figure 3A–3C). Similar results were also obtained from the diverse GEO datasets (Figure 3D–3F). Furthermore, the results demonstrated that the YTHDF2 protein level was increased in glioma tissues (Figure 3G). Collectively, these data indicate that YTHDF2 was up-regulated in glioma tissues and high expression of YTHDF2 was significantly correlated with the tumor grade of LGG.

Figure 3. YTHDF2 was the high expression in LGG. (A–F) The YTHDF2 expression is significantly regulated in glioma examined by the TCGA (A), Rembrandt (B), Gravendeel (C), and GEO datasets (D–F). (G) The immunohistochemistry detection of YTHDF2 in Normal brain tissue, LGG, and HGG.

Relationship between YTHDF2 expression and clinical features

Next, we examine the relationship between YTHDF2 expression and clinical features in LGG. Results confirmed that YTHDF2 differential expression in different histology subtypes (Figure 4A). Additionally, the YTHDF2 expression level was increased with the tumor grade elevated (Figure 4B). We also found that YTHDF2 was markedly reduced in the IDH mutation group and the 1p/19q chromosome co-deletion group (Figure 4C, 4D). Interestingly, YTHDF2 was down-regulation in patients over 40 years old (Supplementary Figure 1A, 1B).

Figure 4. The correlation between the YTHDF2 expression and clinical information in LGG. (A) The expression of YTHDF2 in diverse histology of glioma based on TCGA, CGGA, Rembrandt, and Gravendeel databases. (B) The expression of YTHDF2 in diverse tumor grades of glioma based on TCGA, CGGA, Rembrandt, and Gravendeel databases. (C) The expression of YTHDF2 in diverse IDH mutations of glioma based on TCGA, CGGA, Rembrandt, and Gravendeel databases. (D) The expression of YTHDF2 in diverse 1p/19q codeletion of glioma based on TCGA, CGGA, Rembrandt, and Gravendeel databases.

Prognostic value of YTHDF2

We found that up-regulation of YTHDF2 was correlated with the unfavorable prognosis of glioma patients in different datasets (Figure 5A–5F). Additionally, we perform the ROC curves analysis showed that the area under the curve (AUC) of YTHDF2 were 0.773, 0.815, and 0.789 for the TCGA-LGG datasets (Figure 5G), 0.797, 0.832, and 0.702 for the CGGA datasets (Figure 5H) in 1, 3, and 5 years, respectively. In GBM, the overall of patients no significant difference between high and low YTHDF2 group (Figure 5I). ROC curve of YTHDF2 show an AUC value of 0.982 and 0.931 based on the TCGA-GBM and TCGA-LGG dataset (Figure 5J).

Figure 5. The prognosis of YTHDF2 in LGG. (A–F) The prognosis of YTHDF2 in LGG was examined by TCGA, CGGA, Rembrandt, and Gravendeel databases. (G, H) ROC analyses revealed the predictive value of YTHDF2 in glioma based on TCGA-LGG, CGGA, and TCGA-GBM. (I) The prognostic value of YTHDF2 in GBM. (J) The diagnostic value of YTHDF2 in GBM and LGG.

DNA methylation analysis

DNA methylation plays a crucial role in the regulation of gene expression. To elucidate the abnormal up-regulated mechanisms of YTHDF2 in LGG tissues, we further examine the relationship between DNA methylation, CNV, and the expression of YTHDF2 in gliomas. Firstly, we found that CNV of YTHDF2 was positively correlated with the expression (r = 0.798, p < 0.0001) (Figure 6A). Furthermore, we confirmed that various methylation sites in the promoter region of YTHDF2, among these methylation sites, the cg28734000 was negative correlated with the expression of YTHDF2 in LGG (r = 0.–0.556, p < 0.0001). CGGA dataset also obtained same results (Figure 6B, 6C). 5-Azacytidine, as an inhibitor of DNA methyltransferase [16], was used to treat the U251 cells, the results show that inhibiting the DNA methylation could significantly enhance YTHDF2 RNA levels in U251 cells (Figure 6D). Finally, we uncovered that lower methylation levels on the cg28734000 site correlate with adverse clinical outcomes in the TCGA-LGG dataset examined by the measure database (Figure 6E).

Figure 6. Analysis of the YTHDF2 methylation level. (A) The correlation between CNV, DNA methylation, and YTHDF2 expression is examined by EXPRESS. (B) The DNA methylation level of YTHDF2 in LGG. (C) The DNA methylation level of YTHDF2 in diverse tumor grades of LGG based on CGGA dataset. (D) The expression of YTHDF2 in U251 cells after using 5-AZA treatment was examined by qRT-PCR assay. (E) The prognosis for the methylation level of YTHDF2 in the TCGA-LGG dataset.

YY1 activates YTHDF2 expression in glioma cells

In order to determine the aberrant up-regulation of YTHDF2 in LGG, we employed JASPAR, Knock TF, and PROMO algorithms” to predict the potential transcription factors that may bind to YTHDF2 promoters. The results found that there are 4 common transcription factors in both prediction results, including the STAT3, JUN, MAF, and YY1. Then, we used the GEPIA database to analyze the correlation between these transcription factors and YTHDF2 expression in glioma samples based on the TCGA and CGGA dataset and found that only YY1 was positively associated with YTHDF2 expression. We also found that among these transcription factors, only YY1 was up-regulated in LGG specimens (Figure 7A). Elevated YY1 expression related to poor OS in TCGA-LGG (Figure 7B). Additionally, we found that YY1 expression was significantly increased in diverse WHO grade (p = 0.0006), (Figure 7C), and histology (p = 0.02, Figure 7D). ROC curve analysis of YY1 showed an AUC value of 0.843 in LGG patients (Figure 7E). By analyzing the TCGA LGG and CGGA datasets, we found that YY1 was strongly correlated with YTHDF2 expression (Figure 7F, 7G). Finally, to validate the above results, we performed the qRT-PCR assay, we found out that knock-down of YY1 was markedly reduced the YTHDF2 levels in U251 cells (Figure 7H). These results confirmed that YY1 can increase YTHDF2 expression via modulating the transcription process of YTHDF2 in glioma cells.

Figure 7. YY1 activates YTHDF2 expression in glioma cells. (A) The expression of YY1in LGG tissues and normal tissues. (B) The prognosis of YY1 in TCGA-LGG was examined by GEPIA. (C) The expression of YY1 in diverse tumor grades of LGG examined by CGGA databases. (D) The expression of YY1 in diverse histology of in CGGA-LGG databases. (E) The ROC curve value of YY1 in TCGA-LGG. (F) YY1 expression was strongly positively associated with YTHDF1 expression in TCGA-LGG examined by GEPIA. (G) YY1 expression was strongly positively associated with YTHDF1 expression in LGG examined by CGGA. (H) The expression of YTHDF2 in U251 cells upon depletion of YY1 was examined by qRT-PCR assay.

GSEA analysis

To determine the signaling pathway affected by YTHDF2 in LGG, we performed the GSEA enriched analysis using the GSEA software, results showed that YTHDF2 was mainly involved in the cell cycle, Focal adhesion, Wnt signal pathway, JAK- STAT3-signal pathway (Figure 8).

Figure 8. KEGG signaling pathway explored by GSEA software. The Cell Cycle, Focal adhesion, Wnt signal pathway JAK- STAT3- signal pathway, T cell receptor signal pathway, and B cell receptor signal pathway in LGG were examined by GSEA software.

Association between YTHDF2 expression and immune cell infiltration

Firstly, we analysis the immune infiltration levels based on the different tumor grade, the result demonstrated that different immune infiltration levels observed in the different grade tumor samples (Figure 9A). Furthermore, we found that CNV of YTHDF2 was markedly affect the immune infiltration levels of different immune (Figure 9B). Finally, we found that YTHDF2 was markedly correlated with B cells (cor = 0.577, P = 7.52e-44), CD8+ T cell (cor = 0.382, P = 5.17e-18), CD4+ T cells (cor = 0.373, P = 4.00e-17), macrophages (cor = 0.369 P = 1.06e-16), neutrophils (cor = 0.513, P = 2.68e-33), and dendritic cells (cor = 0.51, P = 7.25e-33) (Figure 9C). We also uncovered that YTHDF2 was positive correlated with the abundances of Th2 cells, T helper cells Tgd, Eosinophils Macrophages, and Neutrophils, negatively associated with the abundances of DCNK, CD56bright cells, and pDC in LGG (Supplementary Figure 1C). In GBM dataset, YTHDF2 was positive correlated with the abundances of Th2 cells, T helper cells, Tgd, and NK cells, negatively associated with the abundances of T cells, DC Macrophages, DC, Neutrophils, and Cytotoxic cells (Supplementary Figure 1D). The Cox model confirmed that B cells and YTHDF2 were significantly related to the poor OS in LGG patients (Figure 9D).

Figure 9. Analysis of the correlation between YTHDF2 expression and diverse immune cell infiltration. (A) The level of diverse immune infiltration in different tumor grades of LGG. (B) The correlation between YTHDF2 expression and somatic copy number alterations. (C) The correlation between YTHDF2 expression and diverse immune cell infiltration. (D) The B cells, CD4+ T cells, CD8+ T cells, dendritic cells, Macrophages, and Neutrophils are correlated with the cumulative survival rate in LGG.

We analyzed the immune checkpoints related gene expression distribution in WHO grade II, III and IV gliomas tissues, the result demonstrated that immune checkpoints related gene differential expression in diverse samples (Figure 10A). In addition, we also analyzed the relationship between expression of YTHDF2 and immune checkpoint-related genes by Pearson correlation analysis in LGGs. Result shown that YTHDF2 expression was significantly positively correlated with CD274 (cor = 0.456, P = 6.65e-28), CTLA4 (cor = 0.464, P = 5.74e-29), HAVCR2 (cor = 0.445, P = 1.9e-26), LAG3 (cor = 0.262, P = 1.59e-09), PDCD1 (cor = 0.594, P = 1.49e-50), PDCD1LG12 (cor = 0.524, P = 8.48e-38), TIGIT (cor = 0.365, P = 9.87e-18), SIGLEC15 (cor = 0.203, P = 3.27e-06) (Figure 10B).

Figure 10. The association between YTHDF2 expression and immune checkpoints. (A) The expression of immune checkpoints-related genes is different from the WHO grade of LGG. (B) The correlations between YTHDF2 expression and diverse immune checkpoints related genes were examined by the TIMER database.

YTHDF2 expressions were correlated with diverse immune cell type markers

We explored the relationship between YTHDF2 expression and various immune cells as well as immune markers in LGG by employing the TISIDB and TIMER database, the TISIDB analysis result indicated that YTHDF2 expression was significantly correlated with 28 tumor-infiltrating lymphocytes (Table 1). Furthermore, the results indicated that YTHDF2 expression was positively correlated with the immune cell marker, in LGG (Table 2).

Drug sensitivity analysis

Given that YTHDF2 has a potential driving effect on the progression of LGG, we used the GSCA tools to analyze the relationship between the YTHDF2 expression and drug sensitivity. Results indicated that YTHDF2 expression positive correlated with the drug sensitivity of XAV939, 17-AAG, Cetuximab, TGX221, Erlotinib, and Trametinib (r > 0.16, p < 0.0001), negative related to the drug sensitivity of NPK76-II-72-1, THZ-2-102-1, PHA-793887, WZ3105, YM201636, GSK1070916, Tubastatin A, I-BET-762, KIN001-102, AT-7519, PIK-93, TPCA-1, BX-912 and UNC0638 in GDSC and CTRP database (r < −0.25, p < 0.0001) (Figure 11A, 11B and Tables 3, 4).

Figure 11. The correlation between YTHDF2 expression and drug sensitivity in LGG. (A) The correlation between the YTHDF2 expression and drug sensitivity was examined by the GDSC database. (B) The correlation between the YTHDF2 expression and drug sensitivity was examined by the CTRP database.

Knock down of YTHDF2 inhibited cell proliferation and migration

To examine the function of YTHDF2 in LGG, we used the qRT-PCR assay to detect the YTHDF2 expression in glioma cell lines. The results show that YTHDF2 was the high expression in glioma cell lines, especially in U251cells (Figure 12A). Owing to the higher expression of YTHDF2 in glioma cell lines, we construction the YTHDF2 knockdown cells and using the qRT-PCR assay verify the knockdown efficiency (Figure 12B). By performing the loss of function for YTHDF2 in U251cell, we found that YTHDF2 knock down was significantly inhibited the cell proliferation of U251 cells by growth curve assay (Figure 12C). Additionally, the decreased migration ability of U251 cells was observed in the knockdown of YTHDF2 (Figure 12D, 12E). We also uncover that YTHDF2 was strongly correlated with Ki-67 expression (r = 0.617, P = 4.21e-25) and vimentin (r = 0.705, P = 2.42e-35) (Figure 12F).

Figure 12. Depletion of YTHDF2 inhibits GBM cell proliferation and migration. (A) The expression of YTHDF2 in normal human astrocytes cells (NHA) and GBM cell lines (U87, A172, and U251). (B) The YTHDF2 knockdown efficiency in U251 was verified by qRT-PCR assay. (C) YTHDF2 knockdown significantly inhibited U251 cell proliferation examined by growth curve assay. (D) YTHDF2 knockdown significantly inhibited U251 cell migration examined by transwell assay. (E) YTHDF2 knockdown significantly inhibited U251 cell migration examined by wound healing assay. (F) Pearson correlation analysis the correlation between YTHDF2 and MKI67, VIM expression in TCGA LGG.