YTHDF2 correlates with tumor immune infiltrates in lower-grade glioma

Immunotherapy is an effective treatment for many cancer types. However, YTHDF2 effects on the prognosis of different tumors and correlation with tumor immune infiltration are unclear. Here, we analyzed The Cancer Genome Atlas and Gene Expression Omnibus data obtained through various web-based platforms. The analyses showed that YTHDF2 expression and associated prognoses may depend on cancer type. High YTHDF2 expression was associated with poor overall survival in lower-grade glioma (LGG). In addition, YTHDF2 expression positively correlated with expression of several immune cell markers, including PD-1, TIM-3, and CTLA-4, as well as tumor-associated macrophage gene markers, and isocitrate dehydrogenase 1 in LGG. These findings suggest that YTHDF2 is a potential prognostic biomarker that correlates with LGG tumor-infiltrating immune cells.


INTRODUCTION
Despite therapeutic advances in recent years, cancer still ranks as a leading cause of death [1]. The Cancer Genome Atlas (TCGA) program and Gene Expression Omnibus (GEO) data, which provide important information for further understanding of tumor biology, are available to users via multiple web-based platforms ( [2][3][4][5][6][7][8][9][10][11]). This knowledge is essential and has already been incorporated into clinical practice, improving our ability to diagnose, treat, and prevent cancer.
In this study, we analyzed YTHDF2 expression and its correlation with the prognosis of cancer patients via a pan-cancer analysis using various web-based platforms. We also investigated the relationship between YTHDF2 expression and tumor-infiltrating immune cells (TIICs) in various cancers. Moreover, we analyzed the correlation of YTHDF2 with isocitrate dehydrogenase 1 (IDH1) in LGG. Finally, we performed the enrichment analysis of YTHDF2 in LGG. These results shed light on the important role of YTHDF2 in LGG and provide an underlying mechanism between YTHDF2 and tumor-immune interactions.

YTHDF2 expression in cancer
We used the Tumor Immune Estimation Resource (TIMER) database to study differences in YTHDF2 expression in tumor tissues and adjacent normal tissues. Figure 1A shows that YTHDF2 expression was substantially higher in BLCA (bladder urothelial carcinoma), breast invasive carcinoma, colon adenocarcinoma, esophageal carcinoma, LUAD (lung adenocarcinoma), stomach adenocarcinoma, prostate adenocarcinoma, and UCEC (uterine corpus endometrial carcinoma) tissues than in adjacent normal tissues. However, YTHDF2 expression was lower in head and neck squamous cell carcinoma, KICH (kidney chromophobe), KIRC (kidney renal clear cell carcinoma), kidney renal papillary cell carcinoma, and LIHC (liver hepatocellular carcinoma) tissues than in adjacent normal tissues. YTHDF2 expression was not expressed substantially between cholangiocarcinoma, lung squamous cell carcinoma, READ (rectum adenocarcinoma), and thyroid carcinoma tissues and adjacent normal tissues. Unfortunately, no data were available on YTHDF2 expression in adjacent normal tissues for the following cancers: adrenocortical carcinoma, DLBC (lymphoid neoplasm diffuse large B-cell lymphoma), GBM (glioblastoma multiforme), LAML (acute myeloid leukemia), LGG (lower-grade glioma), mesothelioma, OV (ovarian serous cystadenocarcinoma), PAAD (pancreatic adenocarcinoma), pheochromocytoma and paragonglioma, SARC (sarcoma), skin cutaneous melanoma, testicular germ cell tumor, thymoma, uterine carcinosarcoma, and uveal melanoma.
To provide a more comprehensive evaluation of YTHDF2 expression in cancers, we used the online database Gene Expression Profiling Interactive Analysis (GEPIA) to compare YTHDF2 expression across 33 TCGA cancer types and in TCGA and GTEx normal tissues. Figure 1B shows that YTHDF2 expression was elevated in many cancers, especially DLBC, GBM, PAAD, and THYM.
We then used the ONCOMINE database to compare YTHDF2 expression in human cancer and corresponding normal samples ( Figure 1C and Supplementary

Prognostic value of YTHDF2 in cancer
We investigated the impact of YTHDF2 expression on survival rates by using the PrognoScan database. The relationships between YTHDF2 expression and prognosis in different cancers are shown in Supplementary Table 2. YTHDF2 expression substantially impacted the prognosis of four cancer types, including brain, breast, colorectal, and soft tissue. However, the impact of YTHDF2 on survival was conflicting in two independent breast cancer cohorts.

AGING
We then used the "survival" TIMER module to confirm the prognostic value of YTHDF2 expression in LGG, LIHC, and SARC (Table 1). We explored the clinical impact of YTHDF2 and corrected for potential confounding factors with a multivariable Cox proportional hazard model. In the univariate analysis, YTHDF2, patient age, and all TIICs (B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils. and DCs) had a significant impact on OS in LGG.
YTHDF2, macrophages, and neutrophils had a significant impact on OS in LIHC, whereas YTHDF2, patient age, and CD4+ T cells had a significant impact on OS in SARC. In the multivariate analysis, we observed significant associations of YTHDF2, patient age, and macrophages with OS in LGG. However, only YTHDF2 was associated with OS in LIHC. In addition, associations between YTHDF2, patient age, CD4+ T cells, and OS were observed in SARC. By using the  LGG, Brain Lower Grade Glioma; LIHC, Liver hepatocellular carcinoma; SARC, Sarcoma; YTHDF2, YTH N6-methyladenosine RNA binding protein 2; Cor, R value of Spearman's correlation; None, correlation without adjustment. Purity, correlation adjusted by purity.P-value Significant Codes: 0 ≤ *** < 0.001 ≤ ** < 0.01 ≤ * < 0.05.
UALCAN database, higher YTHDF2 expression was associated with poor OS in LGG, LIHC, and SARC. YTHDF2 expression also impacted the OS in LGG, LIHC, and SARC with different clinicopathological parameters, such as gender and tumor grade (Supplementary Figure 1 and Supplementary Table 4). Although YTHDF2 expression was not significantly higher in LGG compared with normal samples (Supplementary Figure 2A), we found that YTHDF2 expression was higher in astrocytoma than in oligoastrocytoma and oligodendroglioma. YTHDF2 expression was higher in grade 3 LGG than in grade 2. In addition, higher YTHDF2 expression was associated with poor OS in all LGG and LGG with astrocytoma, but not oligoastrocytoma and oligodendroglioma (Supplementary Figure 2 and Supplementary Table 4).

YTHDF2 expression is correlated with the immune infiltration level in LGG
As stated previously, some TIICs were independent predictors of survival in cancers (Table 1). Therefore, we investigated the correlation of YTHDF2 expression with immune infiltration levels in 32 cancer types from the TIMER database. The analysis showed that YTHDF2 expression was associated with tumor purity AGING in 14 cancer types and B cell infiltration levels in 10 cancer types. In addition, YTHDF2 expression was associated with CD8+ T cell levels in 12 cancer types, CD4+T cell levels in 14 cancer types, macrophage levels in 14 cancer types, neutrophil levels in 12 cancer types, and DC levels in 12 cancer types (Supplementary Table 5).

Correlation analysis between YTHDF2 expression and immune markers
To better understand the relationship between YTHDF2 and various infiltrating immune cells, we analyzed the correlations between YTHDF2 expression and the marker genes of different immune cells and functional T cells in LGG, LIHC, and SARC with the TIMER database. Table 2 shows YTHDF2 expression was associated with most marker genes of the various immune and T cells in LGG. However, YTHDF2 expression was associated with only 14 markers in LIHC and 13 markers in SARC (Table 2).

YTHDF2 expression is correlated with IDH1 level in LGG
IDH1 mutations often occur in gliomas [47,48] and AML [49,50]. In addition, mutant IDH is highly associated with the regulation of the immune microenvironment in LGG [51]. Moreover, YTHDF2 is related to cancer stem cells (CSCs) in AML [52]. We attempted to find the relationship between YTHDF2 and IDH1 expression. We also analyzed the impact of the IDH1 mutation on immune infiltration levels in LGG. Interestingly, data from the GEPIA database showed that high IDH1 expression was associated with poor OS in LGG (HR = 1.7, P = 0.0061) ( Figure 4A).
LGG patients with IDH1 mutations had a superior OS according to the cBioPortal for Cancer Genomics analysis ( Figure 4B). Chinese Glioma Cooperative Group (CGGA) data also indicated that the IDH1 mutation led to a superior OS in glioma ( Figure 4C). However, the IDH1 mutation had no impact on OS in AML ( Figure 4D). In addition, YTHDF2 expression has a moderate positive relationship with IDH1 in LGG ( Figure 4E) and a weak positive relationship with IDH1 in AML ( Figure 4F). YTHDF2 expression was weakly related to TAM-related genes and markers in AML (Supplementary Table 7). More importantly, the levels of infiltration B cells, CD8+ T cells, macrophages, neutrophils, and DCs were higher in IDH1-wild-type LGG than IDH1-mutant LGG ( Figure 4G). These results suggest that YTHDF2 may play an important role in immune infiltration in LGG, especially IDH1wild-type LGG, but not in AML.

Enrichment analysis of YTHDF2 functional networks in LGG
We used the LinkedOmics database to analyze YTHDF2 mRNA sequencing data from 27 LGG patients. The volcano plot in Figure 5A shows that YTHDF2 was positively correlated with 241 genes (dark-red dots) but negatively correlated with 195 genes (dark-green dots) (FDR< 0.05). The 50 significant gene sets positively and negatively associated with YTHDF2 are shown in the heat map ( Figure 5B and 5C). The LinkedOmics GESA tool was used to perform the Gene Ontology and pathway enrichment analyses (Supplementary Table 8 and Figure  5D-5G). Supplementary Table 8 shows that in general the genes correlated with YTHDF2 were enriched in biological processes (double-strand break repair, DNA replication, cell cycle checkpoint, and mitotic cell cycle phase transition), cellular components (DNA packaging complex, protein-DNA complex, nuclear speck, replication fork, and chromosomal region), and molecular function (RNA polymerase II transcription factor binding, repressing transcription factor binding, NF-kappaB binding, nucleosome binding, and alcohol binding). Our results, demonstrating enrichment analyses for the KEGG, Panther, Reactome, and Wiki pathways, show the genes correlated with YTHDF2 were more enriched in cell cycle, TCA cycle, DNA replication, and the FAS signaling pathway.

DISCUSSION
In the present study, we first performed a pan-cancer analysis to analyze YTHDF2 expression and prognostic value. Comprehensive analysis suggested that the differences in YTHDF2 expression and prognostic values in different types of cancer may reflect underlying mechanisms associated with different biological characteristics. Importantly, multivariate analysis confirmed that high YTHDF2 expression was an independent prognostic factor in patients with LGG, LIHC, or SARC. We found that YTHDF2 expression AGING was higher in LGG compared with normal samples, although the difference was not significant.
LGG are a diverse group of primary brain tumors, which mainly include astrocytoma, oligoastrocytoma, and oligodendroglioma. Previous studies have shown that astrocytic tumor type (vs. oligodendroglioma or oligodominant) was a poor prognostic indicator in patients with LGG [53][54][55]. We also found that YTHDF2 expression was higher in astrocytoma than in the other tumor types (oligoastrocytoma and oligodendroglioma). Moreover, high YTHDF2 expression was a prognostic factor in LGG with astrocytoma but not with oligoastrocytoma and oligodendroglioma. Similarly, the expression of YTHDF2 was higher in grade 3 LGG than in grade 2, and high YTHDF2 expression was a prognostic factor in LGG with different tumor grades. These results implied that YTHDF2 was a prognostic factor in LGG, especially with the more malignant subtype or higher tumor grade. However, more research is needed to verify the findings.
A second important finding from this study is that YTHDF2 expression positively correlated with the levels of infiltrating B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and DCs in LGG.
Notably, an association was found between YTHDF2 expression and TAM markers, such as CCL2, CSF1, CSF1R, EGF, STAT3, STAT6, IL-6, IL-10, TLR4, TGFβ (TGFB1), LOX, PD-L1 (CD274), PD-L2 (PDCD1LG2), CD80, CD86, and MFGE8. TAMs play a special role in regulating different steps of tumor progression and metastasis [56]. In glioma, CSCs can induce M2 macrophages, which secrete many cytokines, including TGF-β1 and IL-10, and facilitate immunesuppression [57]. Secretion of IL-10 and TGF-β was shown to facilitate an immunosuppressive microenvironment by inhibiting T cell proliferation in oral squamous cell carcinoma [58]. Interestingly, colonystimulating factor-1 (CSF1) secreted from tumor cells was shown to induce macrophages to produce epidermal growth factor (EGF), which in turn promoted the migration of cancer cells [59]. In addition, inhibition of colony-stimulating factor-1 receptor (CSF1R) in TAMs suppressed the metastasis of pancreatic tumors [60]. The role of TAMs in immunosuppression has been widely studied. For instance, activation of the PD-1/PD-L1/PD-L2 and CTLA4/CD80/CD86 pathways leads to inhibition of TCR signal and T cell cytotoxic functions [61,62]. Previously, it has been suggested that TAMs are attractive therapeutic targets, based on their important role in the tumor immunosuppressive microenvironment in cancer patients [56]. Another interesting finding is the association between YTHDF2 expression and DCs, Treg cells, and T cell exhaustion markers, such as HLA-DPB1, HLA-DQB1, HLA-DRA, HLA-DPA1, TGFβ, and TIM-3. Notably, TIM-3 is a crucial T cell exhaustion regulator [63]. DCs can promote tumor metastasis by increasing Treg cells and reducing CD8+T cell cytotoxicity [64]. In addition, some markers (tumor mutational burden [TMB], PD-1, and PD-L1) have been identified as the effectors of immunotherapy. TMB can be used as a biomarker to identify pediatric glioblastoma (GBM) patients who may benefit from immunotherapy [65]. However, another study found that high TMB is only found in 3.5% of GBM patients, and that IDH1-mutant gliomas are not enriched for high TMB [66]. PD-1 (PDCD1) promoter methylation is a prognostic factor in patients with LGG harboring IDH mutations [20]. A previous study found that PD-L2 expression upregulated in higher grade glioma and IDH-wild-type glioma. High PD-L2 expression was associated with poor survival in GBM [67]. Importantly, several immunotherapies have been evaluated in patients with glioma, including peptide vaccines, DC vaccines, oncolytic viruses, CAR-T cells, and checkpoint inhibitor therapy [68][69][70]. However, a previous study reported the response rates were low in refractory high-grade gliomas treated with PD-1 inhibitors [71]. TIGIT and PD-1 dual checkpoint blockade enhances antitumor immunity and survival in a murine GBM model [72]. Blocking PD-1/PD-L1 interactions together with MLN4924 therapy is a potential strategy for glioma treatment [73]. Gliomas treated with DC vaccination ± murine anti-PD-1 monoclonal antibody blockade or a colony-stimulating factor 1 receptor inhibitor (PLX3397) had prolonged survival in vivo [74]. Previous studies indicate that combination therapy with immune checkpoint blockade is effective for the treatment of malignant tumors, including GBM [75,76].
Our third important finding is that YTHDF2 expression correlated with IDH1 expression in LGG. The analysis showed that high IDH1 expression was associated with poor OS in LGG. IDH1 mutations were associated with a superior OS. This is consistent with previous studies showing that IDH1 mutation is an independent favorable prognostic marker in glioma [47,48]. In addition, the immune infiltration levels were higher in IDH1-wild-type LGG than in IDH1-mutant LGG. We showed that significant infiltration of immune cells, such as B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and DCs, was linked to poor prognosis in LGG. In a previous study, IDH-wild-type gliomas exhibit a more prominent tumor infiltrating lymphocytes than IDH-mutant cases [77]. IDH1 mutations in gliomas caused leukocyte chemotaxis downregulation, resulting in suppression of the tumorassociated immune system [78]. As previously noted, IDH-mutant gliomas can escape the immune surveillance of natural killer cells [79]. More importantly, YTHDF2 expression has a positive AGING relationship with IDH1 level. These results indicate that the role of YTHDF2 in immune infiltration in LGG may depend on IDH1 status. However, further investigations are needed to verify our findings.
Pathway enrichment analysis of YTHDF2 in LGG by GESA found that the genes correlated with YTHDF2 were more significantly enriched in cell cycle, TCA cycle, DNA replication, and the FAS signaling pathway. Interestingly, the most significant gene positively associated with YTHDF2, FAF1, can regulate antiviral immunity ( [80,81]). Moreover, notch family genes (the pathway found in the enrichment analysis) were prognostic biomarkers and correlated with immune infiltrates in gastric cancer ( [82]). Because bioinformatics analysis was performed based on TCGA or GEO datasets, further biological experiments are needed to validate future results.
In summary, our data provide a comprehensive bioinformatics analysis of YTHDF2 expression and prognostic value in human cancers. High YTHDF2 expression correlates with poor prognosis and increased immune infiltration levels (including infiltration of B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils and DCs) in LGG. YTHDF2 expression positively correlated with expression of several immune cell markers, including exhausted T cell markers, PD-1, TIM-3, and CTLA-4 in LGG. In addition, YTHDF2 expression positively correlated with TAM gene markers in LGG. Interestingly, YTHDF2 expression positively correlated with IDH1 expression in LGG. These findings suggest that YTHDF2 is a potential prognostic biomarker and correlates with tumor immune cells infiltration in LGG.

GEPIA database analysis
GEPIA (http://gepia.cancer-pku.cn/index.html) [2] is an interactive web server for analyzing the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from the TCGA and the GTEx projects using a standard processing pipeline. GEPIA was used to analyze YTHDF2 expression and associated survival values (including OS and DFS) of YTHDF2 in 33 different cancer types. Using the Spearman method, correlation between YTHDF2 and IDH1 was determined. YTHDF2 values were represented on the xaxis, and IDH1 values were represented on the y-axis.

TIMER database analysis
The TIMER database (https://cistrome.shinyapps.io/ timer/) [3], which includes 10,897 samples across 32 cancer types from TCGA, is a comprehensive resource for estimating the abundance of six types of infiltrating immune cells, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and DCs. We analyzed YTHDF2 expression in different cancer types via different expression modules and the correlation of YTHDF2 expression with the abundance of immune infiltrates via the gene module. Partial correlations between variables, when considering tumor purity, are shown on the left-most panel of the figure or table [83]. In addition, relationships between YTHDF2 expression and publicly available gene markers of TIICs were explored via correlation modules [84]. The Spearman method was used to determine the correlation coefficient.
The search filters were set as the following: differential analysis (cancer vs normal), cancer type (breast cancer), sample type (clinical specimen), data type (mRNA), and gene (YTHDF2). Thresholds were set as gene rank, 10%; fold change, 2; and P-value, 0.05.

UALCAN database
UALCAN (http://ualcan.path.uab.edu/index.html) [5] is a portal for facilitating tumor subgroup gene expression and survival analyses. It was used to evaluate the mRNA levels and prognostic value of YTHDF2 in LGG patient and normal samples. A P value less than 0.05 was considered significant.

PrognoScan database analysis
The PrognoScan database (http://www.abren.net/ PrognoScan/) [6] was used to analyze the relationships between YTHDF2 expression and patient prognosis, such as OS and DFS, across publicly available cancer microarray datasets.

Kaplan-Meier plotter database analysis
The Kaplan-Meier plotter (http://kmplot.com/analysis/) [7] is capable of assessing the effect of 54,675 genes on survival in 21 cancer types. The correlation between YTHDF2 expression and survival was analyzed by the pan-cancer module of the Kaplan-Meier plotter. The HR with 95% CI and the log-rank P-value were determined.

OncoLnc database analysis
OncoLnc (http://www.oncolnc.org/) [8] is an interactive tool for exploring survival correlations, and for downloading clinical data coupled to expression data AGING for mRNAs, miRNAs, and long noncoding RNAs. The correlation between YTHDF2 expression and survival was analyzed by OncoLnc. The Cox correlation coefficient and P-value were calculated.

CGGA database analysis
A total of 118 glioma samples (82 samples with IDH1 mutation and 37 with wild-type IDH1) from CGGA were analyzed to determine the association of IDH1 with survival [9]. GraphPad Prism software was used to generate a survival curve, and the log-rank test was used to assess the statistical significance.

cBioportal for Cancer Genomics database analysis
The cBioportal Cancer Genomics database (https:// www.cbioportal.org) [10], which was originally developed at Memorial Sloan Kettering Cancer Center, enables users to visualize, analyze, and download largescale cancer genomics datasets. The survival associated with IDH1 alterations in LGG was analyzed, and the log-rank test P-value was calculated. Determination of the correlation between YTHDF2 and IDH1 was performed using the Spearman and Pearson methods.

LinkedOmics dataset
LinkedOmics (http://www.linkedomics.org/login.php) [11] is a publicly available portal that includes multiomics data from all 32 TCGA cancer types. It provides a unique platform for biologists and clinicians to access, analyze, and compare cancer multi-omics data within and across tumor types.