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

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 outcome of m 6 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 6 A reader protein to be identified and studied most widely and affects mRNA stability (Wang X. et al., 2014). YTHDF2 combines specifically with m 6 Acontaining 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 . In pancreatic cancer, knockdown of YTHDF2 expression inhibits cell proliferation and promotes the migration and invasion of cells .
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.

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) 1 database. Datasets showing the profile of mRNA expression of normal liver tissue were obtained from the Genotype-Tissue Expression (GTEx) 2 database. RNAseq data from The Cancer Genome Atlas (TCGA)-LIHC and International Cancer Genome Consortium (ICGC)-LIRI-JP cohorts were obtained from TCGA 3 database and ICGC 4 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) 5 database. We analyzed YTHDF2 expression in different cancer types in the TIMER 6 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.

Analyses of Co-expressed Genes
The LinkedOmics database 7 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 coexpressed 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 8 . PPI pairs were extracted with a minimum required interaction score of 0.4, and the PPI network was visualized via Cytoscape 3.7.1 9 . 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 10 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 11 ].

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 R 5.26. P < 0.05 was considered significant.

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

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.

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 (  (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.

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 (    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).  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 highexpression 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 coexpressed 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).

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.

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.

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).

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)  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.

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 6 A is a predominant internal mRNA modification (Li and Mason, 2014;Roundtree et al., 2017). Studies have shown that methylation of m 6 A RNA is related to the initiation and progression of diseases, including cancer. YTHDF2 is a m 6 A reader. As a key enzyme participating in the regulation of RNA methylation, it can combine specifically with m 6 Acontaining 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 6 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 . 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 6 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 cellcycle 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.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.