Systematic Analyses of the Role of the Reader Protein of N⁶-Methyladenosine RNA Methylation, YTH Domain Family 2, in Liver Hepatocellular Carcinoma

Abstract

Background:

YTH domain family (YTHDF) 2 acts as a “reader” protein for RNA methylation, which is important in tumor regulation. However, the effect of YTHDF2 in liver hepatocellular carcinoma (LIHC) has yet to be elucidated.

Methods:

We explored the role of YTHDF2 in LIHC based on publicly available datasets [The Cancer Genome Atlas (TCGA), International Cancer Genome Consortium (ICGC), and Gene Expression Omnibus (GEO)]. A bioinformatics approach was employed to analyze YTHDF2. Logistic regression analyses were applied to analyze the correlation between YTHDF2 expression and clinical characteristics. To evaluate the effect of YTHDF2 on the prognosis of LIHC patients, we used Kaplan–Meier (K–M) curves. Gene set enrichment analysis (GSEA) was undertaken using TCGA dataset. Univariate and multivariate Cox analyses were used to ascertain the correlations between YTHDF2 expression and clinicopathologic characteristics with survival. Genes co-expressed with YTHDF2 were identified and detected using publicly available datasets [LinkedOmics, University of California, Santa Cruz (UCSC), Gene Expression Profiling Interactive Analysis (GEPIA), and GEO]. Correlations between YTHDF2 and infiltration of immune cells were investigated by Tumor Immune Estimation Resource (TIMER) and GEPIA.

Results:

mRNA and protein expression of YTHDF2 was significantly higher in LIHC tissues than in non-cancerous tissues. High YTHDF2 expression in LIHC was associated with poor prognostic clinical factors (high stage, grade, and T classification). K–M analyses indicated that high YTHDF2 expression was correlated with an unfavorable prognosis. Univariate and multivariate Cox analyses revealed that YTHDF2 was an independent factor for a poor prognosis in LIHC patients. GSEA revealed that the high-expression phenotype of YTHDF2 was consistent with the molecular pathways implicated in LIHC carcinogenesis. Analyses of receiver operating characteristic curves showed that YTHDF2 might have a diagnostic value in LIHC patients. YTHDF2 expression was associated positively with SF3A3 expression, which implied that they may cooperate in LIHC progression. YTHDF2 expression was associated with infiltration of immune cells and their marker genes. YTHDF2 had the potential to regulate polarization of tumor-associated macrophages, induce T-cell exhaustion, and activate T-regulatory cells.

Conclusion:

YTHDF2 may be a promising biomarker for the diagnosis and prognosis of LIHC and may provide new directions and strategies for LIHC treatment.

Introduction

Liver cancer is one of the deadliest types of cancer. About 90% of primary liver cancers are liver hepatocellular carcinoma (LIHC). The latter results in >700,000 deaths and 850,000 new cases annually globally and is the fifth most prevalent cancer type worldwide (Tomita et al., 2011; Torre et al., 2015; Llovet et al., 2016). LIHC remains a major challenge of treatment due to its high morbidity and mortality, lack of effective diagnosis and treatment, and extremely poor prognosis (Mittal and El-Serag, 2013).

Efficacious treatments for LIHC are surgery, chemoradiotherapy, neoadjuvant therapy, local ablation, liver transplantation, or immunotherapy. However, the treatment results are not good, and mortality rates are high (Hsiao et al., 2017; Lee et al., 2017; Shen et al., 2017). There is a lack of effective biomarkers to predict the prognosis. Thus, discovering reliable prognostic biomarkers for LIHC is very important.

The initiation and progression of LIHC are regulated by several factors (El-Serag, 2011; Laursen, 2014). N⁶-Methyladenosine (m⁶A) modification was identified first in the 1970s and is the predominant form of mRNA methylation in eukaryotes (Du et al., 2019). m⁶A modification plays an important part in mRNA splicing; translation of protein biosynthesis; export, localization, and stability of mRNA; virus infection; and stem-cell differentiation (Fu et al., 2014; Meyer and Jaffrey, 2014; Chandola et al., 2015; Mendel et al., 2018). Recent studies have shown that methylation of m⁶A RNA is related to the initiation and progression of diseases, embryonic development (Wang Y. et al., 2014), and cancer (Cui et al., 2017; Kwok et al., 2017; Deng et al., 2018). m⁶A methylation is regulated by a group of proteins, including methyltransferases [“writers” (METTL3, METTL14, and WTAP)], binding proteins {“readers” [YTH domain family (YTHDF) 1, YTHDF2, YTHDF3, YTHDC1, YTHDC2]}, and demethylases [“erasers” (ALKBH5 and FTO)] (Li et al., 2017; Yang et al., 2018).

The outcome of m⁶A-modified mRNAs is based on the YTH (YT521-B homology) domain family protein (Meyer and Jaffrey, 2017). Among them, YTHDF2 was the first m⁶A reader protein to be identified and studied most widely and affects mRNA stability (Wang X. et al., 2014). YTHDF2 combines specifically with m⁶A-containing mRNA, regulates the stability of target RNA, and participates in regulation of a series of physiologic or pathologic processes (Wang X. et al., 2014; Li L. J. et al., 2018). Some studies have reported that YTHDF2 has an important role in acute myeloid leukemia (Paris et al., 2019). Inhibition of YTHDF2 expression advances expansion of hematopoietic stem cells (Li Z. et al., 2018). In pancreatic cancer, knockdown of YTHDF2 expression inhibits cell proliferation and promotes the migration and invasion of cells (Chen et al., 2017).

Using various publicly available databases, we investigated YTHDF2 expression and its relationship with the prognosis of cancer patients. Also, the correlation between YTHDF2 expression and clinical characteristics was investigated, and gene set enrichment analysis (GSEA) and co-expression analyses were undertaken. Finally, the correlation between YTHDF2 and infiltration of immune cells was studied in the tumor environment via the Tumor Immune Estimation Resource (TIMER) database.

We used databases to analyze associations between YTHDF2 expression and the clinical prognosis of patients with LIHC. In this way, we provided suggestions and evidence for further study of the underlying mechanism between YTHDF2 and the tumorigenesis and progression of LIHC. Our study may aid in the use of a promising biomarker for the diagnosis and treatment of LIHC.

Materials and Methods

Dataset Acquisition and Bioinformatics Analyses

RNA-seq transcriptome data and corresponding clinicopathologic information were obtained from publicly available datasets. Datasets showing LIHC expression (GSE14520, GSE63898, and GSE64041) were downloaded from the Gene Expression Omnibus (GEO)¹ database. Datasets showing the profile of mRNA expression of normal liver tissue were obtained from the Genotype-Tissue Expression (GTEx)² database. RNA-seq data from The Cancer Genome Atlas (TCGA)-LIHC and International Cancer Genome Consortium (ICGC)-LIRI-JP cohorts were obtained from TCGA³ database and ICGC⁴ database, respectively. In TCGA-LIHC and ICGC-LIRI-JP cohorts, data on level 3 expression (RNA-seq transcriptome data) of LIHC patients and clinical-feature data of patients were arranged and analyzed further (Supplementary Table 1). The clinical-feature data of patients were extracted retrospectively from medical records (age, sex, clinical stage, histology grade, and Tumor–Node–Metastasis (TMN) classification). All RNA-seq data underwent normalization. Analyses of protein expression in clinical specimens were done using the Human Protein Atlas (HPA)⁵ database. We analyzed YTHDF2 expression in different cancer types in the TIMER⁶ database. Analyses of YTHDF2 expression in TCGA-LICH and ICGC-LIRI-JP cohorts were done using box plots to visualize differences in expression for discrete variables.

Analyses of the Association Between YTHDF2 Expression and Clinical Characteristics, and Survival Analyses

Logistic regression analyses and the Wilcoxon signed-rank test were applied to identify the correlation between YTHDF2 expression and the clinical characteristics of LIHC. Correlation between variables of clinicopathologic characteristics and overall survival (OS) was determined by univariate and multivariate Cox regression analyses. Disease-free survival (DFS) and OS in LIHC patients with high and low expression of YTHDF2 were compared by Kaplan–Meier analyses with log-rank tests using the survival and survminer package within R (R Project for Statistical Computing, Vienna, Austria). To reduce the “noise” of disease-irrelevant deaths, in the ICGC-LIRI-JP cohort, survival duration >3 years was truncated to “3 years,” and the status of the corresponding patient was set to “alive.” In TCGA-LIHC cohort, survival duration >5 years was truncated to 5 years, and the status of the corresponding patient was set to “alive.”

Gene Set Enrichment Analysis

GSEA is a computational method that determines whether a set of genes defined a priori shows significant concordant differences between two biological states (groups of high and low YTHDF2 expression in the present study) (Mootha et al., 2003; Subramanian et al., 2005). In the present study, a Hallmark gene set (h.all.v6.2.symbols.gmt), Biocarta gene set (c2.cp.biocarta.v6.2.symbols.gmt), and Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set (c2.cp.kegg.v6.2.symbols.gmt) were applied to explore the potential mechanisms underlying YTHDF2 expression in LIHC, and GSEA was used to normalize data using the GSEA-3.0 tool (Liberzon et al., 2011). In general, gene sets with a false discovery rate (FDR) <0.25, absolute value of the normalized enrichment score (NES) ≥ 1.0, and normalized P < 0.05 were considered to be enriched significantly.

Co-expression Genes Analysis

Analyses of Co-expressed Genes

The LinkedOmics database⁷ was used to screen genes that were co-expressed with YTHDF2 (absolute value of Pearson’s correlation coefficient > 0.4) and they were shown as heatmaps or volcano plots (Vasaikar et al., 2018). We used “clusterProfiler” within R package to undertake functional annotations for co-expressed genes. The P-adjust and q-value had to be < 0.05 to be considered significant categories in terms used in the Gene Ontology (GO) database and KEGG database. To evaluate the potential protein–protein interaction (PPI) of co-expressed genes, a PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database⁸. PPI pairs were extracted with a minimum required interaction score of 0.4, and the PPI network was visualized via Cytoscape 3.7.1⁹. The Reactome FI plugin within Cytoscape was employed to analyze networks of gene interactions. The top 10 core genes of the gene interaction network and PPI network were identified by the CytoHubba plugin within Cytoscape according to the degree score of each gene node. The Gene Expression Profiling Interactive Analysis (GEPIA) database¹⁰ was used to analyze the correlation between YTHDF2 and SF3A3, and we verified this correlation in other publicly available datasets [GEO; University of California, Santa Cruz (UCSC) Xena¹¹ ].

Infiltration of Immune Cells

TIMER is a comprehensive resource for systematic analyses of infiltration of immune cells across diverse cancer types and includes 10,897 samples across 32 cancer types from TCGA. We analyzed the relationship between YTHDF2 expression and infiltration of immune cells [B cells, cluster of differentiation (CD)4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells (DCs)] in LIHC using the “gene” module of TIMER, as well as the tumor purity. Furthermore, correlations between LAYN expression and marker genes of different immune cells were analyzed via correlation modules.

Statistical Analyses

The Wilcoxon signed-rank test was applied to assess differences in expression of YTHDF2 mRNA between the LIHC group and normal group. Logistic regression analyses were done to evaluate the association between YTHDF2 expression and the clinical characteristics of patients with LIHC. The relationship between clinicopathologic variables and OS was analyzed by Cox regression analyses. In logistic regression, we classified high and low groups according to the median expression of YTHDF2 and SF3A3. In Cox regression analyses, expression of YTHDF2 and expression of SF3A3 were used as continuous variables. To evaluate the diagnostic efficacy of YTHDF2 expression, we drew the receiver operating characteristic (ROC) curve via the pROC package. The correlation between relative gene expression was undertaken by linear regression and Pearson’s correlation analyses. Data analyses were completed by R 3.6.0 and ActivePerl^® 5.26. P < 0.05 was considered significant.

Results

Clinical Data of Patients

Clinical-feature data and mRNA-expression data of LIHC patients were captured from TCGA and ICGC databases on June 2019 (Supplementary Table 1). In the TCGA-LIHC cohort, the clinical data of 348 LIHC patients were extracted retrospectively from medical records and comprised median age at the diagnosis (61 years), clinical stage (I = 173, II = 84, III = 86, and IV = 5), histology grade (G1 = 45, G2 = 171, G3 = 119, and G4 = 13), pathology classification for M (M0 = 267, M1 = 4, and MX = 77), pathology classification for T (T1 = 175, T2 = 86, T3 = 77, and T4 = 10), and pathology classification for N (N0 = 256, N1 = 4, and NX = 88). In the ICGC-LIRI-JP cohort, only the median age at the diagnosis (69 years), clinical stage (I = 40, II = 117, III = 80, and IV = 23), survival duration, and survival status of 260 LIHC patients were recorded.

High YTHDF2 Expression in Liver Hepatocellular Carcinoma

Initially, we used an online tool (TIMER) to investigate YTHDF2 expression in different cancer types. In multiple cancer types (Figure 1A), YTHDF2 expression in cancer tissues was significantly higher than that in adjacent normal tissues. These cancer types were bladder urothelial carcinoma, cholangiocarcinoma, colon adenocarcinoma, esophageal carcinoma, head and neck cancer, chromophobe renal cell carcinoma, renal clear cell carcinoma, renal papillary cell carcinoma, LIHC, lung adenocarcinoma, adenocarcinoma of the prostate gland, gastric adenocarcinoma, carcinoma of the thyroid gland, and uterine corpus endometrial carcinoma. We focused mainly on the role of YTHDF2 in LIHC. To demonstrate the relationship between YTHDF2 expression and clinical relevance, we further analyzed YTHDF2 expression in publicly available datasets [HPA, TCGA-LIHC cohort, and ICGC-LIRI-JP (validation cohort)] (Figure 1). LIHC patients (n = 374), normal samples (n = 50), and 50 LIHC tissues and corresponding adjacent non-tumor tissues were included in TCGA-LIHC cohort. The ICGC-LIRI-JP cohort comprised 260 LIHC patients, 202 normal samples, and 199 LIHC tissues, and the corresponding adjacent non-tumor tissues were included in TCGA-LIHC cohort. YTHDF2 expression was significantly higher in LIHC tissues than that in non-tumor tissues (P < 0.05) (Figures 1B,C) or paired non-tumor tissues (P < 0.05) (Figures 1D,E). In the HPA, we analyzed expression of YTHDF2 protein. We discovered that compared with YTHDF2 staining in normal liver tissue, YTHDF2 staining was stronger in LIHC tissue (Figure 1F). Next, in other exposed datasets (GEO and GTEx), we assessed YTHDF2 expression in LIHC tissues and normal tissues (Table 1). The results were consistent with those in TCGA-LIHC and ICGC-LIRI cohorts. These data showed that YTHDF2 expression was upregulated in LIHC.

FIGURE 1

High Expression of YTHDF2 Is Associated With Clinical Characteristics

There were a total of 375 and 232 LIHC patients with YTHDF2 expression and clinical information in TCGA-LIHC and ICGC-LIRI, respectively. Then, we evaluated the association between YTHDF2 expression and the clinical characteristics of LIHC patients. Univariate logistic regression analyses showed that YTHDF2 expression was correlated substantially with poor prognostic clinical factors (grade, T classification, and stage) in LIHC patients (Table 2). In TCGA-LIHC-cohort, high expression of YTHDF2 was connected significantly with high clinical stage [odds ratio (OR) = 2.56 for III–IV vs. I; 2.16 for II vs. I], histology grade (1.66 for poor vs. good), and T classification (2.51 for T3–T4 vs. T1) (P < 0.05 for all). In the ICGC-LIRI cohort, a high clinical stage (IV vs. I) was correlated significantly with high YTHDF2 expression. High expression of YTHDF2 was related significantly to T classification, stage, and histology grade (P < 0.05 for all) (Figures 2A–D). These results demonstrated that in LIHC patients, the high-expression group of YTHDF2 tended to progress to a more advanced grade and stage than the low-expression group of YTHDF2.

FIGURE 2

YTHDF2 Is an Independent Predictor of Poor Survival (Disease-Free Survival and Overall Survival) in Liver Hepatocellular Carcinoma

Univariate Cox analyses indicated that the prognosis was related to the clinical stage, T classification, M classification, and YTHDF2 expression (Table 3). High expression of YTHDF2 was correlated significantly with a poor prognosis [TCGA-LIHC cohort: OS, hazard ratio (HR) = 2.05, 95% confidence interval (CI): 1.39–3.04, P < 0.001; DFS, 1.72, 1.23–2.40, P < 0.001; ICGC-LIRI cohort: OS, 1.06, 1.01–1.11, P < 0.001]. Multivariate analyses demonstrated that YTHDF2 expression was an independent predictor for survival in LIHC patients (TCGA-LIHC cohort: OS, HR = 1.80, 95% CI: 1.20–2.80; P < 0.05; DFS, 1.50, 1.07–2.10, P < 0.05. ICGC-LIRI cohort: OS, 1.34, 1.01–1.79, P < 0.05) (Figures 3A–C). OS and DFS of high- and low-expression groups of YTHDF2 were compared by Kaplan–Meier analyses with log-rank tests in LIHC patients. Kaplan–Meier curves revealed that LIHC patients with high YTHDF2 expression were correlated with a poor prognosis (P < 0.05) (Figures 3D–F).

FIGURE 3

Functional Annotation of YTHDF1

Gene Set Enrichment Analysis and Potential Function of YTHDF2 Expression

We wished to investigate the potential function of YTHDF2 expression in LIHC. We undertook GSEA and used the MSigDB Collection (c2.cp.biocarta.v6.2.symbols.gmt, h.all.v6.2.symbols.gmt, c2.cp.kegg.v6.2.symbols.gmt) to examine the enriched gene sets at different expression of YTHDF2. We selected the most significantly enriched signaling pathways based on the NES (Supplementary Table 2). GSEA suggested that mitotic spindle, protein secretion, integrin pathway, Wnt pathway, meiosis of oocytes, regulation of autophagy, spliceosome, and vascular endothelial growth factor, as well as expression of mammalian target of rapamycin (mTOR), phosphoinositide 3-kinase/protein kinase B/mTOR, and p53 pathways were enriched differentially in the YTHDF2 high-expression group. These results were consistent with the molecular pathways implicated in LIHC carcinogenesis (Couri and Pillai, 2019). Due to the limited space, only six pathways of high expression are listed in Supplementary Figure 1, but all indicated the potential role of YTHDF2 in LIHC development.

Co-expression of YTHDF2 mRNA in Liver Hepatocellular Carcinoma

We wished to investigate further the underlying regulation of YTHDF2 expression in LICH. Thus, we analyzed the co-expression profiles with YTHDF2 in LIHC using the LinkedOmics database. Seventy-eight genes were co-expressed significantly with YTHDF2 mRNA (absolute value of Pearson’s correlation ≥ 0.4, FDR > 0.05) (Supplementary Table 3). The volcano plot (Figure 4A) demonstrated that the global genes were correlated with YTHDF2 as identified by Pearson’s test. Heatmaps (Figures 4B,C) revealed the top 50 genes negatively and positively associated with YTHDF2 in LIHC. Analyses using the GO database showed that the genes co-expressed with YTHDF2 were located mainly in nuclear speckles, nuclear body, nucleoplasm, and nuclear lumen (Figure 4E), where they participated primarily in binding of nucleic acids and RNA (Figure 4F), thereby suggesting effects on RNA processing, RNA splicing, and metabolic processes (Figure 4D). Pathway analyses using the KEGG database revealed that these genes were enriched mainly in spliceosome pathways (Figure 4G).

FIGURE 4

Construction of a Network of YTHDF2 Co-expressed Genes and Identification of Potential “Hub” Genes

A PPI network of 78 co-expressed genes was constructed using the STRING database and visualized via Cytoscape 3.7.1 (Figure 5A). We used the Reactome FI plugin within Cytoscape for analyses of gene interaction networks (Figure 5B). The top 10 core genes (Table 4 and Figures 5C,D) were identified by the CytoHubba plugin within Cytoscape according to the degree score of each gene node. We explored the biological processes of five duplicate hub genes using the BINGO plugin of Cytoscape. These hub genes participated mainly in the cellular metabolic processes of nitrogen compounds and nucleic acid metabolism, suggesting that they may affect protein biosynthesis in cells (Figure 5E). Among these core genes, the highest degree score was for SF3A3, so it may be the most important key core gene.

FIGURE 5

Co-expression of YTHDF2 and SF3A3

Correlation analyses using the GEPIA database showed that expression of YTHDF2 and SF3A3 had a high correlation coefficient (R = 0.8, P < 0.05) (Figure 6A). In other publicly available datasets (GEO and cBioPortal), a positive correlation between YTHDF2 and SF3A3 transcripts was validated (Figures 6B–F). Finally, we confirmed this phenomenon further in the UCSC Xena database (Figure 6G). These results suggested that YTHDF2 expression and SF3A3 expression had a strong positive correlation, suggesting that they may be “functional partners” in LIHC.

FIGURE 6

SF3A3 Expression and Prognostic Value in Liver Hepatocellular Carcinoma

We wished to analyze gene expression of SF3A3 in LIHC, so we used publicly available datasets. Dataset analyses indicated that SF3A3 expression was significantly higher in LIHC tissues than in non-tumor tissues (Figures 7A–F). In the HPA, we analyzed expression of SF3A3 protein: compared with SF3A3 staining in normal liver tissue, SF3A3 staining was stronger in LIHC tissue (Figure 7G). Subsequently, we investigated the prognostic value of SF3A3 in TCGA-LIHC and ICGC-LIRI-JP cohorts, and we confirmed that high expression of SF3A3 was associated significantly with a decrease in OS and DFS in LIHC (Figures 7H–J).

FIGURE 7

Diagnostic Value of High YTHDF2 Expression in Liver Hepatocellular Carcinoma

ROC curves were used to assess the diagnostic value of YTHDF2 in LIHC patients and healthy individuals from GTEx + TCGA-LIHC and ICGC-LIRI datasets. In the normal group vs. all tumor groups, the area under the ROC curve (AUC) of YTHDF2 was 0.65 in the GTEx + TCGA-LIHC, and ICGC-LIRI cohorts (Figures 8A,E). Furthermore, we assessed the diagnostic value of YTHDF2 in a subgroup of LIHC patients with regard to clinical stage. Analyses of ROC curves demonstrated that high YTHDF2 expression might have a diagnostic value for LIHC patients of different clinical stages: the AUC was 0.63 and 0.64 for stage I, 0.68 and 0.60 for stage II, and 0.65 and 0.72 for stage III + IV in the TCGA and ICGC cohorts, respectively (Figure 8). Combined diagnosis of YTHD2 and SF3A3 could improve the diagnostic efficiency significantly, with all AUC values > 0.85 (Figure 8). These results suggested that expression of YTHDF2 and SF3A3 might have a diagnostic value for LIHC patients.

FIGURE 8

YTHDF2 Expression Is Correlated With Infiltration of Immune Cells in Liver Hepatocellular Carcinoma

Tumor-infiltrating lymphocytes are independent prognostic elements of survival and can reflect metastasis to sentinel lymph nodes (Ohtani, 2007; Azimi et al., 2012). To further broaden the understanding of YTHDF2, we investigated if YTHDF2 expression was associated with the degree of infiltration of immune cells in LIHC using the TIMER database. YTHDF2 expression had a significant correlation with various “immune signatures,” which included infiltration levels of B cells (r = 0.19, P = 3.70 × 10^–4), CD8+ T cells (0.20, 1.83 × 10^–4), CD4 + T cells (0.37, 8.97 × 10^–13), macrophages (0.42, 7.23 × 10^–16), neutrophils (0.48, 3.02 × 10^–21), and DCs (0.34, 6.71 × 10^–4) in LIHC (Figure 9A). These results strongly indicated that YTHDF2 expression was a crucial participant in infiltration of immune cells in LIHC, particularly macrophages and neutrophils.

FIGURE 9

Correlation Analyses Between Immune Marker Genes and YTHDF2 Expression

Immune cells can be activated in different states (Fridlender et al., 2009; Murray, 2017). We focused on the relationship between YTHDF2 and marker genes of various types of immune cells [CD8+ T cells, monocytes, T cells (general), M1 and M2 macrophages, B cells, neutrophils, tumor-associated macrophages (TAMs), natural killer cells, AND DCs] in LIHC in GEPIA and TIMER databases. We also analyzed T cells with different functions, such as T-helper type 1 (Th1), Th17, Th2, T-regulatory (Treg), T follicular helper, and exhausted T cells (Table 5). After tumor purity-correlated adjustments, we found a significant correlation between YTHDF2 expression and most marker genes of various types of immune cells and different functional T cells in LIHC. YTHDF2 expression was strongly correlated with expression of monocytes, TAMs, DCs, and M1 and M2 macrophages (Table 5). Expression of C–C motif chemokine ligand-2, interleukin (IL) 10 and CD68 of TAMs, interferon regulatory factor-5, and cyclo-oxygenase (COX) 2 of the M1 phenotype, as well as MS4A4A, VSIG4, and CD163 of the M2 phenotype was correlated significantly with YTHDF2 expression in LIHC (P < 0.0001) (Figures 9B-E). In addition, we assessed expression of YTHDF2 and the immune markers stated above in GEPIA and ICGC-LIRI databases. The correlation of YTHDF2 expression with marker genes of monocytes and TAMs was consistent with that in the TIMER database (Table 6). These findings suggested that YTHDF2 may regulate macrophage polarization in LIHC.

High expression of YTHDF was associated with high infiltration of DCs. Expression of DC marker genes (e.g., HLA-DRA, HLA-DRA1, BDCA-4, and CD11c) was correlated significantly with YTHDF2 expression. Furthermore, for Tregs, expression of FOXPM3 and TGFB1 was correlated positively with YTHDF2 expression in LIHC. DCs promote tumor metastasis by increasing the number of Tregs and inhibiting the cytotoxicity of CD8 + T cells (Sawant et al., 2012). Whether YTHDF2 is a driving factor mediating DCs and tumor metastases merits further study.

We also observed that YTHDF2 expression had a significant correlation with the marker genes of Tregs and T-cell exhaustion (e.g., CCR8, FOXP3, STAT5B, and TIM-3) (Table 5). FOXP3 can inhibit cytotoxic T cells attacking cancer cells. TIM-3 is a key gene that regulates T-cell exhaustion and has a high positive correlation with YTHDF2 expression, indicating that high expression of YTHDF2 has a crucial role in TIM-3-mediated T-cell exhaustion. Hence, these results further demonstrate that YTHDF2 expression was correlated significantly with infiltration of immune cells in LIHC. Hence, in the tumor microenvironment, YTHDF2 has an important role in immune escape. YTHDF2 may also affect the occurrence and development of LIHC by regulating the infiltration of immune cells of different phenotypes.

Discussion

RNA methylation is the major modification of RNA. m⁶A is a predominant internal mRNA modification (Li and Mason, 2014; Roundtree et al., 2017). Studies have shown that methylation of m⁶A RNA is related to the initiation and progression of diseases, including cancer. YTHDF2 is a m⁶A reader. As a key enzyme participating in the regulation of RNA methylation, it can combine specifically with m⁶A-containing mRNA and can regulate the stability of target RNA (Wang X. et al., 2014; Li L. J. et al., 2018). The reader proteins of m⁶A methylation are YTHDF1, YTHDF2, YTHDF3, YTHDC1, and YTHDC2. We found that only YTHDF1 and YTHDF2 were related to the survival of LIHC patients. The prognostic role of YTHDF1 in LIHC has been reported (Zhao et al., 2018). Studies on the potential prognostic role of YTHDF2 in LIHC are lacking: the present study was conducted to fill this knowledge gap. We comprehensively analyzed the expression and role of YTHDF2 in a large cohort of LIHC patients. We employed bioinformatics to analyze R-seq data on TCGA-LIHC (training group) and ICGC-LIRI-JP (validation group) cohorts. We demonstrated, for the first time, that high expression of YTHDF2 in LIHC was correlated significantly with clinicopathologic variables (high clinical stage, histology grade, and T classification) and a poor prognosis.

Differences in expression of YTHDF2 mRNA between the LIHC group and normal group were analyzed in multiple LIHC datasets. Then, Kaplan–Meier survival analyses, univariate and multivariate Cox regression analyses, and ROC curve analyses were used to ascertain the value of YTHDF2. Logistic regression analyses were applied to identify the relationship between YTHDF2 expression and clinical characteristics. In addition, GSEA, as well as analyses of co-expression, and infiltration of immune cells were done to study how YTHDF2 participates in the evolution and progression of LIHC.

Studies have shown that knockdown of YTHDF2 expression inhibits cell proliferation and promotes the migration and invasion of pancreatic cancer cells (Chen et al., 2017). Zhong et al. (2019) showed that YTHDF2 can inhibit the proliferation and growth of cells by destroying the stability of epidermal growth factor receptor mRNA in LIHC. Chen et al. (2018) demonstrated that RNA m⁶A-METTL3 promotes LIHC progression via YTHDF2-dependent post-transcriptional silencing of SOCS2 expression. The conclusions of the two studies mentioned above are inconsistent, which suggests the dual role of YTHDF2 in tumor cells.

Here, we revealed that YTHDF2 expression was significantly higher in LIHC than in normal liver tissue and was correlated significantly with clinicopathologic variables (high clinical stage, histology grade, and T classification). Analyses of Kaplan–Meier survival curves showed that high YTHDF2 expression was associated with a poor prognosis (OS and DFS) in LIHC patients. Univariate and multivariate Cox regression analyses demonstrated that YTHDF2 was an independent factor for a poor prognosis in LIHC patients. Hence, YTHDF2 may be a valuable prognostic biomarker for patients with LIHC. The GSEA results of the high-expression phenotype of YTHDF2 are consistent with the molecular pathways implicated in LIHC pathogenesis. ROC curve analyses showed that YTHDF2 might have a diagnostic value in LIHC patients. Taken together, our data suggest the potential role of YTHDF2 in LIHC development.

Co-expression analyses demonstrated that YTHDF2 expression was associated positively with SF3A3 expression. SF3A3 is obtained from purified spliceosomes and is an important component of the SF3A RNA splicing complex (Chiara et al., 1994). SF3A3 participates in the replication and repair of DNA and transcriptional regulation (Yang et al., 2012). SF3A3 can also inhibit the activity of the tumor-suppressor gene p53. Silencing SF3A3 expression can increase p53 expression, thereby inducing the cell-cycle arrest and death of tumor cells (Siebring-van Olst et al., 2017). In addition, the interaction between SF3A3 and CSR1 is essential for cell death (Zuo et al., 2017). Therefore, we analyzed SF3A3 expression in LIHC tissues and non-tumor tissues: SF3A3 expression in LIHC tissues was significantly higher than that in non-tumor tissues. Survival results showed that high expression of SF3A3 was associated with a poor prognosis in LIHC. We speculate that YTHF2 and SF3A3 may participate jointly in the development and progression of LIHC.

Another important aspect of this research was the correlation between YTHDE2 expression and immune signatures in LIHC. YTHDF2 expression was associated with infiltration of immune cells and marker genes of immune cells (Tables 5, 6). The correlation between YTHDF2 expression and markers of M2 macrophages was slightly stronger than that of markers of M1 macrophages and was related to TAM markers, which suggested the potential regulatory role of YTHDF2 in TAM polarization. Furthermore, YTHDF2 expression was correlated positively with markers of Tregs (TGFB1, CCR8, FOXP3, and STAT5B), which indicated that YTHDF2 could activate Tregs. YTHDF2 expression was also related to TIM-3 (a key gene for T-cell exhaustion) expression. Thus, we speculate that YTHDF2 has the potency to induce T-cell exhaustion. Taken together, these results show that YTHDF2 has an important role in the recruitment and regulation of infiltrating immune cells in LIHC.

More in-depth experiments and clinical trials are needed to verify the value of YTHDF2 in LIHC. One limitation of our study was that we did not consider potential confounders for the association between YTHDF2 expression and survival (e.g., type of radiotherapy/chemotherapy and location of patient enrollment).

Conclusion

High expression of YTHDF2 was associated with a poor prognosis and, together with increased infiltration of immune cells, could be an independent prognostic indicator for patients with LIHC. GSEA results were consistent with the molecular pathways implicated in LIHC pathogenesis in a high-expression phenotype of YTHDF2. Co-expression analyses confirmed the upregulation of expression of YTHDF2 and SF3A3 in LIHC, both of which were associated with a poor prognosis of LIHC. In addition, YTHDF2 was involved in regulation of TAMs and T_regs. Thus, YTHDF2 may have an important role in infiltration of immune cells and help to improve immunomodulatory therapy in LIHC patients.

We believe that YTHDF2 may be a promising biomarker for the diagnosis and prognosis of LIHC patients and may provide new directions and strategies for their treatment. However, the specific regulatory mechanisms merit further exploration.

References

• 1

  AzimiF.ScolyerR. A.RumchevaP.MoncrieffM.MuraliR.McCarthyS. W.et al (2012). Tumor-infiltrating lymphocyte grade is an independent predictor of sentinel lymph node status and survival in patients with cutaneous melanoma.J. Clin. Oncol.302678–2683. 10.1200/jco.2011.37.8539

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  ChandolaU.DasR.PandaB. (2015). Role of the N6-methyladenosine RNA mark in gene regulation and its implications on development and disease.Brief. Funct. Genomics14169–179. 10.1093/bfgp/elu039

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  ChenJ.SunY.XuX.WangD.HeJ.ZhouH.et al (2017). YTH domain family 2 orchestrates epithelial-mesenchymal transition/proliferation dichotomy in pancreatic cancer cells.Cell Cycle162259–2271. 10.1080/15384101.2017.1380125

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  ChenM.WeiL.LawC. T.TsangF. H.ShenJ.ChengC. L.et al (2018). RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2.Hepatology672254–2270. 10.1002/hep.29683

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  ChiaraM. D.Champion-ArnaudP.BuvoliM.Nadal-GinardB.ReedR. (1994). Specific protein-protein interactions between the essential mammalian spliceosome-associated proteins SAP 61 and SAP 114.Proc. Natl. Acad. Sci. U.S.A.916403–6407. 10.1073/pnas.91.14.6403

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  CouriT.PillaiA. (2019). Goals and targets for personalized therapy for HCC.Hepatol. Int.13125–137. 10.1007/s12072-018-9919-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  CuiQ.ShiH.YeP.LiL.QuQ.SunG.et al (2017). m(6)A RNA methylation regulates the self-renewal and tumorigenesis of glioblastoma stem cells.Cell. Rep.182622–2634. 10.1016/j.celrep.2017.02.059

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  DengX.SuR.WengH.HuangH.LiZ.ChenJ. (2018). RNA N(6)-methyladenosine modification in cancers: current status and perspectives.Cell Res.28507–517. 10.1038/s41422-018-0034-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  DuK.ZhangL.LeeT.SunT. (2019). m(6)A RNA methylation controls neural development and is involved in human diseases.Mol. Neurobiol.561596–1606. 10.1007/s12035-018-1138-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  El-SeragH. B. (2011). Hepatocellular carcinoma.N. Engl. J. Med.3651118–1127. 10.1056/NEJMra1001683

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  FridlenderZ. G.SunJ.KimS.KapoorV.ChengG.LingL.et al (2009). Polarization of tumor-associated neutrophil phenotype by TGF-beta: “N1” versus “N2” TAN.Cancer Cell16183–194. 10.1016/j.ccr.2009.06.017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  FuY.DominissiniD.RechaviG.HeC. (2014). Gene expression regulation mediated through reversible m6A RNA methylation.Nat. Rev. Genet.15293–306. 10.1038/nrg3724

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  HsiaoJ. H.TsaiC. C.LiangT. J.ChiangC. L.LiangH. L.ChenI. S.et al (2017). Adjuvant hepatic arterial infusion chemotherapy is beneficial for selective patients with Hepatocellular carcinoma undergoing surgical treatment.Int. J. Surg.4535–41. 10.1016/j.ijsu.2017.07.071

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  KwokC. T.MarshallA. D.RaskoJ. E.WongJ. J. (2017). Genetic alterations of m(6)A regulators predict poorer survival in acute myeloid leukemia.J. Hematol. Oncol.10:39. 10.1186/s13045-017-0410-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  LaursenL. (2014). A preventable cancer.Nature516S2–S3. 10.1038/516S2a

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  LeeJ.LeeJ. H.YoonH.LeeH. J.JeonH.NamJ. (2017). Extraordinary radiation super-sensitivity accompanying with sorafenib combination therapy: what lies beneath?Radiat. Oncol. J.35185–188. 10.3857/roj.2017.00262

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  LiA.ChenY. S.PingX. L.YangX.XiaoW.YangY.et al (2017). Cytoplasmic m(6)A reader YTHDF3 promotes mRNA translation.Cell. Res.27444–447. 10.1038/cr.2017.10

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  LiL. J.FanY. G.LengR. X.PanH. F.YeD. Q. (2018). Potential link between m(6)A modification and systemic lupus erythematosus.Mol. Immunol.9355–63. 10.1016/j.molimm.2017.11.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  LiS.MasonC. E. (2014). The pivotal regulatory landscape of RNA modifications.Annu. Rev. Genomics Hum. Genet.15127–150. 10.1146/annurev-genom-090413-025405

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  LiZ.QianP.ShaoW.ShiH.HeX. C.GogolM.et al (2018). Suppression of m(6)A reader Ythdf2 promotes hematopoietic stem cell expansion.Cell. Res.28904–917. 10.1038/s41422-018-0072-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  LiberzonA.SubramanianA.PinchbackR.ThorvaldsdóttirH.TamayoP.MesirovJ. P. (2011). Molecular signatures database (MSigDB) 3.0.Bioinformatics271739–1740. 10.1093/bioinformatics/btr260

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  LlovetJ. M.Zucman-RossiJ.PikarskyE.SangroB.SchwartzM.ShermanM.et al (2016). Hepatocellular carcinoma.Nat. Rev. Dis. Primers2:16018. 10.1038/nrdp.2016.18

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  MendelM.ChenK. M.HomolkaD.GosP.PandeyR. R.McCarthyA. A.et al (2018). Methylation of structured RNA by the m(6)a writer METTL16 is essential for mouse embryonic development.Mol. Cell.71986.e11–1000.e11. 10.1016/j.molcel.2018.08.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  MeyerK. D.JaffreyS. R. (2014). The dynamic epitranscriptome: N6-methyladenosine and gene expression control.Nat. Rev. Mol. Cell. Biol.15313–326. 10.1038/nrm3785

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  MeyerK. D.JaffreyS. R. (2017). Rethinking m(6)a readers, writers, and erasers.Annu. Rev. Cell. Dev. Biol.33319–342. 10.1146/annurev-cellbio-100616-060758

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  MittalS.El-SeragH. B. (2013). Epidemiology of hepatocellular carcinoma: consider the population.J. Clin. Gastroenterol.47(Suppl.(0)), S2–S6. 10.1097/MCG.0b013e3182872f29

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  MoothaV. K.LindgrenC. M.ErikssonK. F.SubramanianA.SihagS.LeharJ.et al (2003). PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes.Nat. Genet.34267–273. 10.1038/ng1180

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  MurrayP. J. (2017). Macrophage polarization.Annu. Rev. Physiol.79541–566. 10.1146/annurev-physiol-022516-034339

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  OhtaniH. (2007). Focus on TILs: prognostic significance of tumor infiltrating lymphocytes in human colorectal cancer.Cancer Immun.7:4.

  • Google Scholar

• 30

  ParisJ.MorganM.CamposJ.SpencerG. J.ShmakovaA.IvanovaI.et al (2019). Targeting the RNA m(6)a reader YTHDF2 selectively compromises cancer stem cells in acute myeloid leukemia.Cell Stem Cell25137.e6–148.e6. 10.1016/j.stem.2019.03.021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  RoundtreeI. A.EvansM. E.PanT.HeC. (2017). Dynamic RNA modifications in gene expression regulation.Cell1691187–1200. 10.1016/j.cell.2017.05.045

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  SawantA.HenselJ. A.ChandaD.HarrisB. A.SiegalG. P.MaheshwariA.et al (2012). Depletion of plasmacytoid dendritic cells inhibits tumor growth and prevents bone metastasis of breast cancer cells.J. Immunol.1894258–4265. 10.4049/jimmunol.1101855

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  ShenJ.HeL.LiC.WenT.ChenW.LuC.et al (2017). Nomograms to predict the individual survival of patients with solitary hepatocellular carcinoma after hepatectomy.Gut Liver11684–692. 10.5009/gnl16465

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  Siebring-van OlstE.BlijlevensM.de MenezesR. X.van der Meulen-MuilemanI. H.SmitE. F.van BeusechemV. W. (2017). A genome-wide siRNA screen for regulators of tumor suppressor p53 activity in human non-small cell lung cancer cells identifies components of the RNA splicing machinery as targets for anticancer treatment.Mol. Oncol.11534–551. 10.1002/1878-0261.12052

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  SubramanianA.TamayoP.MoothaV. K.MukherjeeS.EbertB. L.GilletteM. A.et al (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.Proc. Natl. Acad. Sci. U.S.A.10215545–15550. 10.1073/pnas.0506580102

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  TomitaM.ShimizuT.AyabeT.YoneiA.OnitsukaT. (2011). Preoperative neutrophil to lymphocyte ratio as a prognostic predictor after curative resection for non-small cell lung cancer.Anticancer Res.312995–2998.

  • Google Scholar

• 37

  TorreL. A.BrayF.SiegelR. L.FerlayJ.Lortet-TieulentJ.JemalA. (2015). Global cancer statistics, 2012.CA Cancer J. Clin.6587–108. 10.3322/caac.21262

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  VasaikarS. V.StraubP.WangJ.ZhangB. (2018). LinkedOmics: analyzing multi-omics data within and across 32 cancer types.Nucleic Acids Res.46D956–D963. 10.1093/nar/gkx1090

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  WangX.LuZ.GomezA.HonG. C.YueY.HanD.et al (2014). N6-methyladenosine-dependent regulation of messenger RNA stability.Nature505117–120. 10.1038/nature12730

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  WangY.LiY.TothJ. I.PetroskiM. D.ZhangZ.ZhaoJ. C. (2014). N6-methyladenosine modification destabilizes developmental regulators in embryonic stem cells.Nat. Cell. Biol.16191–198. 10.1038/ncb2902

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  YangY.HsuP. J.ChenY. S.YangY. G. (2018). Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism.Cell. Res.28616–624. 10.1038/s41422-018-0040-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  YangZ. M.ChenW. W.WangY. F. (2012). [Study on gene differential expressions of substance and energy metabolism in chronic superficial gastritis patients of Pi deficiency syndrome and of pi-wei hygropyrexia syndrome].Zhongguo Zhong Xi Yi Jie He Za Zhi321180–1187.

  • Google Scholar

• 43

  ZhaoX.ChenY.MaoQ.JiangX.JiangW.ChenJ.et al (2018). Overexpression of YTHDF1 is associated with poor prognosis in patients with hepatocellular carcinoma.Cancer Biomark.21859–868. 10.3233/cbm-170791

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  ZhongL.LiaoD.ZhangM.ZengC.LiX.ZhangR.et al (2019). YTHDF2 suppresses cell proliferation and growth via destabilizing the EGFR mRNA in hepatocellular carcinoma.Cancer Lett.442252–261. 10.1016/j.canlet.2018.11.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  ZuoZ. H.YuY. P.MartinA.LuoJ. H. (2017). Cellular stress response 1 down-regulates the expression of epidermal growth factor receptor and platelet-derived growth factor receptor through inactivation of splicing factor 3A3.Mol. Carcinog.56315–324. doi: 10.1002/mc.22494

  • Pubmed Abstract
  • CrossRef
  • Google Scholar