Online ISSN: 1945-4589
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 14, Issue 15 pp 6269—6298

The comprehensive expression and functional analysis of m6A modification “readers” in hepatocellular carcinoma

Sha Qin1,2, , Gaoming Liu3, , Haoer Jin1,2, , Xue Chen7, , Jiang He5, , Juxiong Xiao4, , Yan Qin4, , Yitao Mao4,6, , Luqing Zhao1,2,6, ,

  • 1 Department of Pathology, Xiangya Hospital, Central South University, Changsha, Hunan, China
  • 2 Department of Pathology, School of Basic Medical Science, Xiangya School of Medicine, Central South University, Changsha, Hunan, China
  • 3 Xiangya School of Medicine, Central South University, Changsha, Hunan, China
  • 4 Department of Radiology, Xiangya Hospital, Central South University, Changsha, Hunan, China
  • 5 Center for Molecular Medicine, Xiangya Hospital, Central South University, Changsha, Hunan, China
  • 6 National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, Hunan, China
  • 7 Early Clinical Trial Center, Hunan Cancer Hospital and The Affiliated Cancer Hospital of Xiangya School of Medicine, Central South University, Changsha, Hunan, China

Received: April 1, 2022       Accepted: July 21, 2022       Published: August 12, 2022      

https://doi.org/10.18632/aging.204217
How to Cite

Copyright: © 2022 Qin et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

N6-methyladenosine (m6A) modification regulators are essential for the diagnosis and treatment of various cancers. However, the comprehensive analysis about roles of m6A “readers” in hepatocellular carcinoma (HCC) remains unclear. UALCAN, GEPIA2, HPA, Kaplan Meier plotter, cBioPortal, STRING WebGestalt, Metascape and TIMER 2.0 database and Cytoscape software were used to comprehensively analyze the bioinformatic data. We found that m6A “readers” were upregulated at the mRNA level and protein level in HCC patients. Highly expressed YTHDF1, IGF2BP3 and NKAP were positively correlated with advanced HCC stage and had a poor prognosis in OS and PFS. The gene alterations of m6A “readers” happened frequently, and YTHDF3 had the highest mutation rate. The function of m6A “readers” on HCC may be closely correlated with splicing related proteins (including HNRNP family, SNRP family, and SR family), metabolic process, protein binding and RNA splicing related signaling pathways. Moreover, although the correlation of YTHDF3 and CD8+ T cell infiltration, and the correlation of IGF2BP3 and infiltration of mast cells and CAF are negative, most m6A “readers” had a positive correlation with immune cells (including CD8+ T cell, CD4+ T cell, Tregs, B cell, neutrophil, monocyte, macrophage, myeloid dendritic cell, nature killer cell, mast cell, and CAF). Macrophages, CD4+ T cell, Treg, B cell, monocyte, and myeloid dendritic cell had a positively strong correlation (Rho>0.4) with most m6A “readers” (such as YTHDC1, YTHDC2, YTHDF1, IGF2BP3, HNRNPA2B1 and HNRNPC). In conclusion, by comprehensive analysis of m6A “readers”, we found that they were involved in the prognosis of HCC, and m6A “readers” might regulate the development and progression of HCC by participating in some metabolism-related and RNA splicing-related signaling pathways as well as immune cell infiltration.

Introduction

Primary liver cancer was reported as the sixth most commonly diagnosed cancer as well as the third leading cause of cancer death worldwide in 2020. Statistical data shows that it was about 906,000 new cases and 830,000 deaths. And liver hepatocellular carcinoma (HCC) accounts for 75% to 85% of primary liver cancer [1]. With the progress of medical technology, the treatment of HCC has been greatly developed. The widespread use of semi-annual surveillance, increasing the proportion of tumors eligible for treatment, and systemic therapies, including emerging immunotherapies, can improve overall survival for HCC patients [2]. Currently, small molecule inhibitors and immunotherapy are emerging as promising treatment options of HCC [3, 4]. With the development of immunotherapy, therapeutic agents are used to target immune cells such as CD4+ CD25+ Foxp3+ regulatory T cells, which led to a huge achievement in HCC treatment [5]. And the biomarker-driven therapies may significantly improve the survival of patients at advanced stages of HCC [6]. Therefore, the identification and combined use of more therapeutic targets and biomarkers may bring new ideas for the treatment of HCC.

N6-methyladenosine (m6A) modification has been a research hot spot in recent years. m6A modification is the most abundant modification in eukaryotic RNA, with roles including maintenance of RNA stability, mRNA precursor splicing, translation initiation, nuclear export, degradation, etc. [7]. It was found that m6A modifications are usually enriched in the stop codon and 3′UTR regions [8]. m6A modifications include m6A “writer”, “eraser” and “reader”, whose functions are to add m6A modifications, remove m6A modifications, and recognize m6A modifications on RNA, respectively. And they were reported to provide some possibilities for the early diagnosis as well as treatment of various cancers [9, 10]. m6A “readers”, which mainly including YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, NKAP, IGF2BP1, IGF2BP2, and IGF2BP3, are a kind of protein that could recognize the m6A modification in RNA, thus playing roles in RNA metabolism, tumorigenesis, hematopoiesis, viral replication, immune response, and adipogenesis [11]. Nowadays, an increasing number of researches have demonstrated that some m6A “readers” have a potential to predict the prognosis of HCC, such as YTHDF1 [12, 13], and YTHDC2 [14]. However, the comprehensive analysis of m6A “readers” to predict prognosis and diagnosis value of HCC has not been reported.

In this paper, we sort out a comprehensive bioinformatics analysis of the m6A “readers” (including YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP) in HCC. In addition, we discussed their values to be potential therapeutic targets and prognostic biomarkers through various public bioinformatics databases. Furthermore, we selected some key genes from the co-expression genes of m6A “readers”. And these key genes as well as some immune cells infiltration may provide may provide new ideas for the subsequent study of the mechanism of m6A modification in HCC.

Results

The up-regulation of the m6A “readers” expression in HCC patients

We used UALCAN database to check the mRNA expression level of m6A “readers” between HCC and normal liver tissue. The expression level of m6A “readers” such as YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were explored. To our surprise, all of these “readers” expressed significantly higher in HCC than that in normal tissues (Figure 1). Immunohistochemistry could provide us with more information at the tissue level. To determine the differential expression profile of m6A “readers” between normal tissues and HCC., we collected the typical images from the HPA database and analyzed the immunohistochemical staining results by comparing their integrated optical density (IOD) values. (Figure 2). From the Human Protein Atlas (HPA) database, YTHDC1, YTHDF2, IGF2BP1, HNRNPA2B1, HNRNPC and NKAP were detected in live cell nuclear, while YTHDC2, IGF2BP2, IGF2BP3 were expressed in liver cell cytoplasm and membrane. However, the protein expression data of YTHDF1 and YTHDF3 were missed in HPA database. In addition, these stain results showed that the protein expression level of YTHDC1, YTHDC2, YTHDF2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were higher in HCC tissue than in normal liver tissue. These results were consistent with the differential expression profile of m6A “readers” between normal tissues and HCC (as we showed in Figure 1).

The mRNA expression profile of m6A “readers” in HCC and normal liver tissues (UALCAN database). m6A “readers” YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were expressed significantly higher in HCC than that in normal tissues.

Figure 1. The mRNA expression profile of m6A “readers” in HCC and normal liver tissues (UALCAN database). m6A “readers” YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were expressed significantly higher in HCC than that in normal tissues.

The protein expression of the m6A “readers” in HCC and normal liver tissues (the Human Protein Atlas database). The protein expression level of YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were higher in HCC tissue than in normal liver tissue. In addition, the integrated optical density (IOD) was used to analyze the expression difference between HCC tissues and normal liver tissues.

Figure 2. The protein expression of the m6A “readers” in HCC and normal liver tissues (the Human Protein Atlas database). The protein expression level of YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were higher in HCC tissue than in normal liver tissue. In addition, the integrated optical density (IOD) was used to analyze the expression difference between HCC tissues and normal liver tissues.

UALCAN was used to evaluate the association between the mRNA expression level of the m6A “readers” and the pathological stages of HCC patients (Figure 3). we used the WebGestalt database YTHDF2, YTHDF3, IGF2BP2, IGF2BP3, HNRNPA2B1 was positively correlated with tumor stage (stage1, stage2, and stage3).Additionally, the mRNA expression levels of IGF2BP1, HNRNPC and NKAP were positively correlated with tumor stage (stage1, stage2, stage3, and stage4). To further verify the correlation between m6A “readers” and tumor stage, we used GEPIA2 database to explore this correlation. And the results showed that YTHDC1, YTHDF1, YTHDF2, IGF2BP3, HNRNPA2B1, and NKAP were involved in the stage of HCC (Figure 4). Similarly, we also focused on the relationship between mRNA expression level of m6A “readers” and HCC lymph node metastasis. Propensity score matching was used to eliminate confounding factors between groups. And we found that the mRNA expression level of m6A “readers” was high in HCC patients with lymph node metastasis. But, only YTHDC1, YTHDF2 and HNRNPC showed significant correlations with HCC metastasis (Figure 5). In addition, there is a lack of data between m6A “reader” and multiple lymph node metastases in HCC, which may be due to insufficient sample size or other reasons.

The relationship between m6A “readers” mRNA expression and pathological stage of HCC patients (UALCAN database). The mRNA expression level of YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP2, IGF2BP3, and HNRNPA2B1 were positively correlated with tumor stage (stage1, stage2, and stage3). The mRNA expression levels of IGF2BP1, HNRNPC and NKAP were positively correlated with tumor stage (stage1, stage2, stage3, and stage4).

Figure 3. The relationship between m6A “readers” mRNA expression and pathological stage of HCC patients (UALCAN database). The mRNA expression level of YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP2, IGF2BP3, and HNRNPA2B1 were positively correlated with tumor stage (stage1, stage2, and stage3). The mRNA expression levels of IGF2BP1, HNRNPC and NKAP were positively correlated with tumor stage (stage1, stage2, stage3, and stage4).

The relationship between m6A “readers” mRNA expression and pathological stage of HCC patients (GEPIA2 database). The mRNA expression level of YTHDC1, YTHDF1, YTHDF2, IGF2BP3, HNRNPA2B1 and NKAP were correlated with tumor stage.

Figure 4. The relationship between m6A “readers” mRNA expression and pathological stage of HCC patients (GEPIA2 database). The mRNA expression level of YTHDC1, YTHDF1, YTHDF2, IGF2BP3, HNRNPA2B1 and NKAP were correlated with tumor stage.

The relationship between m6A “readers” mRNA expression and lymph node metastasis of HCC patients (UALCAN database). YTHDF2 and HNRNPC showed significant difference between normal and N1. YTHDC1 showed significant difference between N0 and N1, while YTHDC2, YTHDF1, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, and NKAP showed significant difference between normal and N0.

Figure 5. The relationship between m6A “readers” mRNA expression and lymph node metastasis of HCC patients (UALCAN database). YTHDF2 and HNRNPC showed significant difference between normal and N1. YTHDC1 showed significant difference between N0 and N1, while YTHDC2, YTHDF1, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, and NKAP showed significant difference between normal and N0.

Prognostic value of m6A “readers” in HCC patients

Next, we analyzed the prognostic values of mRNA expression level of m6A “readers” in HCC patients through Kaplan-Meier plotter. It was obvious that HCC patients with higher mRNA levels of YTHDF1, YTHDF2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPC, NKAP displayed shorter overall survival (OS) time (Figure 6). The mRNA expression level of YTHDC1, YTHDF1, IGF2BP2, IGF2BP3, HNRNPA2B1, and NKAP were negatively correlated with progression-free survival (PFS) time in HCC patients. However, high expression of YTHDF3 seems to be a protective factor for PFS in HCC patients (Figure 7). Next, we combined the OS and PFS results of m6A “reader” in HCC. And we found that YTHDF1, IGF2BP2, IGF2BP3 and NKAP had significant value to predict both OS and PFS time in HCC patients.

The overall survival curve of m6A “readers” in HCC patients (Kaplan-Meier plotter database). The high expression level of YTHDF1, YTHDF2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPC, NKAP in HCC patients displayed shorter OS time.

Figure 6. The overall survival curve of m6A “readers” in HCC patients (Kaplan-Meier plotter database). The high expression level of YTHDF1, YTHDF2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPC, NKAP in HCC patients displayed shorter OS time.

The progression free survival curve of m6A “readers” in HCC patients (Kaplan-Meier plotter database). YTHDC1, YTHDF1, IGF2BP2, IGF2BP3, HNRNPA2B1, and NKAP were negatively related with PFS time in HCC patients. YTHDF3 was positively related with PFS time in HCC patients.

Figure 7. The progression free survival curve of m6A “readers” in HCC patients (Kaplan-Meier plotter database). YTHDC1, YTHDF1, IGF2BP2, IGF2BP3, HNRNPA2B1, and NKAP were negatively related with PFS time in HCC patients. YTHDF3 was positively related with PFS time in HCC patients.

Genetic alteration of the m6A “readers” in HCC patients

To further analyze the functions of m6A “readers” in HCC patients, we focused on the alteration profiles of m6A “readers”. cBioPortal database was used to evaluate the genetic alterations of the m6A “readers” (Figure 8). According to the results in Figure 8, all of these m6A “readers” had some genetic alterations, which including “inframe mutation”, “missense mutation”, “splice mutation”, “truncating mutation”, “amplification”, “deep deletion”, “mRNA high” and “mRNA low”. And among these alterations, “mRNA high”, “mRNA low”, and “amplification” were the most common types. YTHDF3 had the highest mutation rate, with a mutation rate of up to 25%. In addition, the mutation rates of YTHDF1 (18%), YTHDF2 (11%), IGF2BP1 (10%), IGF2BP3 (10%), and NKAP (12%) were more than 10%. And the mutation rates of YTHDC1 (7%), YTHDC2 (8%), IGF2BP2 (9%), HNRNPA2B1 (8%), and HNRNPC (9%) were more than 5%.

Genetic alteration of the m6A “readers” in HCC patients (cBioPortal database). All of these m6A “readers” had some genetic alteration, including “inframe mutation”, “missense mutation”, “splice mutation”, “truncating mutation”, “amplification”, and “deep deletion”.

Figure 8. Genetic alteration of the m6A “readers” in HCC patients (cBioPortal database). All of these m6A “readers” had some genetic alteration, including “inframe mutation”, “missense mutation”, “splice mutation”, “truncating mutation”, “amplification”, and “deep deletion”.

Interaction analyses of the m6A “readers” in HCC patients

To know the regulatory mechanisms of m6A “readers” in HCC, we picked out the top one co-expression genes of m6A “readers” (Figure 9). As is showed in Figure 9, ELF2, DMXL1, GMEB2, PPP1R8, VCPIP1, DNMT3A, LRRC1, MYBL2, SRSF1, SNRPD1 and RNF113A were picked out. And all of these genes present a positive correlation to m6A “readers”, with the spearman correlation coefficient greater than at least 0.4. Gene co-expression networks can help us study those genes with unknown functions [15]. So, we further explored the top 20 co-expression genes of each m6A “reader” through cBioPortal database (Supplementary Table 1). All of these co-expression genes had a high spearman’s correlation value (>0.4) with m6A “readers”. In addition, we found the spearman’s correlation value among most of these co-expression genes and m6A “readers” are more than 0.6. More importantly, the spearman’s correlation between VCPIP1 and YTHDF3 was more than 0.85. To learn the correlations among these co-expression genes, we put them into STING database. Then, we constructed a protein-protein interactions (PPIs) with these co-expression genes (Figure 10A), and found out the top 10 hub genes as well as the shortest paths (Figure 10B). As it was shown in Figure 10B, the top 10 hub genes were HNRNPA1, HNRNPK, HNRNPL, HNRNPH1, SNRPG, SNRPE, SNRPF, SNRPD1, SRSF1, SRSF3. To our surprise, these hub genes originated from three families (HNRNP family, SNRP family and SR family), which were all found to be the regulatory proteins for selective splicing.

The top one co-expression genes of m6A “readers” in HCC (cBioPortal database). All of these m6A “readers” had a positive correlation with top one co-expression genes (including ELF2, DMXL1, GMEB2, PPP1R8, VCPIP1, DNMT3A, LRRC1, MYBL2, SRSF1, SNRPD1, and RNF113A).

Figure 9. The top one co-expression genes of m6A “readers” in HCC (cBioPortal database). All of these m6A “readers” had a positive correlation with top one co-expression genes (including ELF2, DMXL1, GMEB2, PPP1R8, VCPIP1, DNMT3A, LRRC1, MYBL2, SRSF1, SNRPD1, and RNF113A).

Co-expression network analysis of the m6A “readers” related co-expression genes in HCC. (A) The PPI network based on top 20 co-expression genes of each m6A “reader” from cBioPortal database. This network was edited by STRING database and Cytoscape software. The spearman's correlation value between m6A co-expression genes was depicted in different color (purple for more than 0.4, pink for more than 0.6, blue for more than 0.5, orange for more than 0.7 and red for more than 0.8). (B) The top 10 hub genes (little circle in the core) and their shortest paths of these co-expression genes. These data were analyzed by Cytoscape database.

Figure 10. Co-expression network analysis of the m6A “readers” related co-expression genes in HCC. (A) The PPI network based on top 20 co-expression genes of each m6A “reader” from cBioPortal database. This network was edited by STRING database and Cytoscape software. The spearman's correlation value between m6A co-expression genes was depicted in different color (purple for more than 0.4, pink for more than 0.6, blue for more than 0.5, orange for more than 0.7 and red for more than 0.8). (B) The top 10 hub genes (little circle in the core) and their shortest paths of these co-expression genes. These data were analyzed by Cytoscape database.

Functional enrichment analysis of the m6A “readers”

Functional enrichment analysis allows hundreds of genes to be allocated in different pathways thus contributing to the understanding of an unknown gene. In this study, we used the WebGestalt database to evaluate the biological functions of the m6A “readers” from pan-cancer analysis.

From the GO analysis, it was clearly that the most highly enriched biological process (BP) category was metabolic process, followed by biological regulation, cellular component organization and response to stimulus. In the cellular component (CC) categories, nucleus, membrane-enclosed lumen, protein-containing complex, membrane, and cytosol were highly enriched. As for the molecular function (MF) category, the co-expression genes of m6A “readers” were mainly enriched in protein binding, nucleic acid binding, and ion binding (Figure 11A).

Functional enrichment analysis of the m6A “readers” in HCC. (A) GO enrichment analysis (molecular functions, biological processes and cell components) of the co-expression genes of m6A “readers”. These data were collected from WebGestalt database. (B) Bar graph of KEGG enrichment using WebGestalt website. The co-expressed gene of m6A “readers” were enriched in sister chromatid segregation, nuclear chromosome segregation, RNA splicing via transesterification reactions with bulged adenosine, mRNA splicing via spliceosome, RNA splicing via transesterification reactions and mRNA processing. (C) Bar graph of KEGG enrichment using Metascape. The co-expressed gene of m6A “readers” were enriched in sister chromatid segregation and mRNA processing pathways.

Figure 11. Functional enrichment analysis of the m6A “readers” in HCC. (A) GO enrichment analysis (molecular functions, biological processes and cell components) of the co-expression genes of m6A “readers”. These data were collected from WebGestalt database. (B) Bar graph of KEGG enrichment using WebGestalt website. The co-expressed gene of m6A “readers” were enriched in sister chromatid segregation, nuclear chromosome segregation, RNA splicing via transesterification reactions with bulged adenosine, mRNA splicing via spliceosome, RNA splicing via transesterification reactions and mRNA processing. (C) Bar graph of KEGG enrichment using Metascape. The co-expressed gene of m6A “readers” were enriched in sister chromatid segregation and mRNA processing pathways.

In addition, from the KEGG pathway results (Figure 11B), it was clearly that co-expression genes were enrichment in sister chromatid segregation, nuclear chromosome segregation, RNA splicing via transesterification reactions with bulged adenosine, mRNA splicing via spliceosome, RNA splicing via transesterification reactions and mRNA processing. We found that these co-expressed genes were enriched in several splicing-related pathways. This result echoed the functions of hub genes in PPI network. Besides, we used Metascape database to confirm this result. As is shown in Figure 11C, sister chromatid segregation and mRNA processing pathways were highly enriched. The detailed information was in Table 1.

Table 1. Representative top 20 clusters enriched for co-expressed genes of m6A “readers”.

GOCategoryDescription
GO:0000819GO Biological Processessister chromatid segregation
GO:0006397GO Biological ProcessesmRNA processing
GO:0008608GO Biological Processesattachment of spindle microtubules to kinetochore
GO:0006913GO Biological Processesnucleocytoplasmic transport
R-HSA-194315Reactome Gene SetsSignaling by Rho GTPases
GO:0022613GO Biological Processesribonucleoprotein complex biogenesis
GO:0007052GO Biological Processesmitotic spindle organization
GO:0044773GO Biological Processesmitotic DNA damage checkpoint
GO:0071824GO Biological Processesprotein-DNA complex subunit organization
CORUM:5615CORUMEmerin complex 52
GO:0016569GO Biological Processescovalent chromatin modification
GO:1903312GO Biological Processesnegative regulation of mRNA metabolic process
GO:0034502GO Biological Processesprotein localization to chromosome
GO:0006281GO Biological ProcessesDNA repair
GO:0046822GO Biological Processesregulation of nucleocytoplasmic transport
GO:0044839GO Biological Processescell cycle G2/M phase transition
CORUM:3082CORUMDGCR8 multiprotein complex
CORUM:1183CORUMCDC5L complex
R-HSA-162599Reactome Gene SetsLate Phase of HIV Life Cycle
R-HSA-2990846Reactome Gene SetsSUMOylation

In addition, we integrated data from Kaplan-Meier plotter database, UALCAN database and GEPIA2 database, and it was not difficult to find that YTHDF1, IGF2BP3 and NKAP were involved in OS, PFS, and stage in HCC (Supplementary Figure 1A). So, we analyzed their functional enrichment by Metascape based on the top 20 co-expression genes of them (Supplementary Figure 1B–1D). The results showed that YTHDF1 and IGF2BP3 were both related to cell cycle pathway, and NKAP was involved in protein location process. To further understand the roles of YTHDF1, IGF2BP3 and NKAP in other cancers, we collected their expression levels from pan-cancer analysis (Supplementary Figure 2A–2C). The data was obviously showed that YTHDF1, IGF2BP3 and NKAP were highly expressed in most cancer types.

Immune cell infiltration of the m6A “readers” in HCC

To explore the correlation between the m6A “readers” and immune cell infiltration in HCC, we searched some information from the TIMER 2.0 database (Figure 12 and Supplementary Figures 3, 4). As is shown, we discussed 11 types of immune cells including CD8+ T cell, CD4+ T cell, Tregs, B cell, neutrophil, monocyte, macrophage, myeloid DC, NK cell, mast cell, and CAF. Tumor purity is a major confounding factor in immune cell infiltration analysis. Therefore, all of our results had a purity adjustment. And in general, most of these 11 immune cells had a positive correlation with m6A “readers”, such as YTHDC1 (Figure 12A), YTHDF1 (Figure 12C), YTHDF2 (Figure 12D), IGF2BP1 (Supplementary Figure 3A), IGF2BP2 (Supplementary Figure 3B), HNRNPA2B1 (Supplementary Figure 4A), HNRNPC (Supplementary Figure 4B) and NKAP (Supplementary Figure 4C). However, the CD4+ T cell was negatively related to YTHDC2 (Figure 12B); the CD8+ T cell was negatively related to YTHDF3 (Figure 12E); the mast cell was negatively related to IGF2BP3 (Supplementary Figure 3C). In addition, there were no significant difference between IGF2BP1 and neutrophil. When we explored the relationship between immune cells and m6A “readers” such as IGF2BP3 and NKAP, the TIMER 2.0 database mixed macrophages and monocytes together for analysis. These may because some subpopulations of macrophages originated from adult monocytes [16]. In addition, we found that macrophages had a positively strong correlation (Rho>0.4) with most m6A “readers”, such as YTHDC1, YTHDC2, YTHDF1, YTHDF2, IGF2BP3, HNRNPA2B1, and NKAP. Besides, CD4+ T cell had a positively strong correlation with YTHDC2, IGF2BP2, IGF2BP3, and HNRNPC; Treg had a positively strong correlation with YTHDC1, IGF2BP2, HNRNPA2B1, and HNRNPC; B cell had a positively strong correlation with YTHDC1, HNRNPA2B1, and HNRNPC, monocyte had a positively strong correlation with YTHDF1, IGF2BP3, HNRNPA2B1, and HNRNPC, and DC had a positively strong correlation with YTHDC1, YTHDF1, IGF2BP3, HNRNPA2B1, and HNRNPC; Neutrophil had a positively strong correlation with YTHDC2; CAF had a positively strong correlation with YTHDC1 and YTHDF1; while CD8+ T cell, NK, and mast cell existed a much weaker correlation with all of these m6A “readers”. HNRNPC, HNRNPA2B1 and YTHDC1 had positive correlation with more than four immune cells, including Tregs, B cell, macrophage, DC and so on. However, YTHDF3, IGF2BP1, and NKAP showed weaker correlation with all of the immune cells we researched. More importantly, we obtained the most related immune cells of m6A “readers”. For example, the most related immune cell of HNRNPC was myeloid dendritic cell (Rho = 0.535). More detailed information was shown in Table 2 (the correlation coefficients greater than 0.4 have been bolded, and the negative correlation coefficients have been bolded and italicized).

Immune cell infiltration of the m6A “readers” (YTHDC1/2, and YTHDF1/2/3) in HCC. (A) All of these 11 immune cells had a positive correlation with YTHDC1. (B) All these immune cells we researched had positive correlations with YTHDC2, except for CD4+ T cell. (The CD4+ T cell was negatively related to YTHDC2). (C) All of these 11 immune cells had a positive correlation with YTHDF1. (D) All of these 11 immune cells had a positive correlation with YTHDF2. (E) All these immune cells we researched had positive correlations with YTHDF3, except for CD8+ T cell. (The CD8+ T cell was negatively related to YTHDF3). These data were collected from TIMER 2.0 database.

Figure 12. Immune cell infiltration of the m6A “readers” (YTHDC1/2, and YTHDF1/2/3) in HCC. (A) All of these 11 immune cells had a positive correlation with YTHDC1. (B) All these immune cells we researched had positive correlations with YTHDC2, except for CD4+ T cell. (The CD4+ T cell was negatively related to YTHDC2). (C) All of these 11 immune cells had a positive correlation with YTHDF1. (D) All of these 11 immune cells had a positive correlation with YTHDF2. (E) All these immune cells we researched had positive correlations with YTHDF3, except for CD8+ T cell. (The CD8+ T cell was negatively related to YTHDF3). These data were collected from TIMER 2.0 database.

Table 2. The spearman's correlation value between the m6A “readers” and immune cells.

T cell CD8+T cell CD4+TregsB cellNeutrophilMonocyteMacrophageDCNK cellMast cellCAF
YTHDC10.2460.3270.4060.4200.3410.3480.5050.4390.3590.2830.456
YTHDC20.2200.4590.3420.3290.4550.3460.4200.3180.2160.1570.364
YTHDF10.2860.2960.3800.3980.3200.4140.4290.4770.3300.2870.416
YTHDF20.1850.2740.3280.3200.3860.3870.4490.3770.2220.2370.352
YTHDF3−0.2360.3490.2200.2920.2630.2450.3840.2350.2490.2780.309
IGF2BP10.1720.2630.2640.334−0.2550.2720.2800.2250.2520.142
IGF2BP20.2550.4180.4080.3810.1270.3460.3460.3480.2370.2230.269
IGF2BP30.2730.4640.3750.3040.1100.4160.4160.4010.281−0.165−0.236
HNRNPA2B10.2970.3810.4120.4390.3190.4400.4640.4760.3910.3120.374
HNRNPC0.3400.4530.4140.4730.3160.4000.4170.5350.3670.3700.382
NKAP0.2270.3160.2190.3110.2310.3550.3550.3720.2980.1430.263

Validation the expression level and functions of IGF2BP3 in HCC cell lines

To verify the reliability of the results, we carried out further studies using IGF2BP3 as a representative gene (because of its association with OS, PFS, and stage in HCC). First, we used immunofluorescence results from HPA database and demonstrated that IGF2BP3 was mainly localized in the cytoplasm (Figure 13A). IGF2BP3 was stained in green, the nucleus was stained in blue, and the microtubules were stained in red. This immunofluorescence information was based on A-431 (Figure 13A left) and U-251 (Figure 13A right) cell lines. In addition, CCLE database was used to display the mRNA expression level of IGF2BP3 in multiple HCC cell lines (Figure 13B).

Validation the expression of IGF2BP3 in cell line models. (A) The immunofluorescence results of IGF2BP3. IGF2BP3 was localized in the cytoplasm. (IGF2BP3 was stained in green, the nucleus was stained in blue, and the microtubules were stained in red.) The data were derived from the HPA database. (B) The mRNA expression level of IGF2BP3 in HCC cell lines. The data were derived from CCLE database.

Figure 13. Validation the expression of IGF2BP3 in cell line models. (A) The immunofluorescence results of IGF2BP3. IGF2BP3 was localized in the cytoplasm. (IGF2BP3 was stained in green, the nucleus was stained in blue, and the microtubules were stained in red.) The data were derived from the HPA database. (B) The mRNA expression level of IGF2BP3 in HCC cell lines. The data were derived from CCLE database.

Based on the expression level of IGF2BP3 from CCLE database, we stably knocked down the IGF2BP3 mRNA level in Hep3B cell, and the knocked-down (KD) efficiency were checked by RT-PCR (Figure 14A) and Western Blot (Figure 14B). In addition, the migration rate of IGF2BP3 stably KD cells were significantly reduced compared to the non-targeted control (NC) (Figure 14C, 14D). Besides, we performed the sphere formation assay to verify stemness characteristics of IGF2BP3 in HCC. As it was showed in Figure 14E, 14F, IGF2BP3 stably KD cells showed significantly reduced sphere formation capacity.

Validation the functions of IGF2BP3 in HCC cell line. (A) The knock-down efficiency of IGF2BP3 mRNA expression in Hep3B cell. (B) The knock-down efficiency of IGF2BP3 protein expression in Hep3B cell. (C) The representative migration images of IGF2BP3 stably KD and NC Hep3B cells. (D) The statistical analysis of migration assay for Hep3B cells with IGF2BP3 stably KD and NC. *P ***P ****P E) The representative sphere formation assay images of IGF2BP3 stably KD and NC Hep3B cells. (F) The statistical analysis of sphere formation assay for Hep3B cells with IGF2BP3 stably KD and NC. *P ***P ****P

Figure 14. Validation the functions of IGF2BP3 in HCC cell line. (A) The knock-down efficiency of IGF2BP3 mRNA expression in Hep3B cell. (B) The knock-down efficiency of IGF2BP3 protein expression in Hep3B cell. (C) The representative migration images of IGF2BP3 stably KD and NC Hep3B cells. (D) The statistical analysis of migration assay for Hep3B cells with IGF2BP3 stably KD and NC. *P < 0.05; ***P < 0.001; and ****P < 0.0001. (E) The representative sphere formation assay images of IGF2BP3 stably KD and NC Hep3B cells. (F) The statistical analysis of sphere formation assay for Hep3B cells with IGF2BP3 stably KD and NC. *P < 0.05; ***P < 0.001; and ****P < 0.0001.

Discussion

The biological importance of m6A modifications mainly depends on the roles of m6A “readers” which control RNA fate and function. Studies showed that some m6A “readers” were involved in the development of HCC. For instance, Tanabe A and his colleague proved that through activation of c-Jun and ATF-2, YTHDC2 played a significant role in promoting HCC cell growth [17]. In addition, Liu et al. implicated that YTHDC2 might be a target to improve the therapeutic efficacy for HCC patients [14]. Zhao et al. believed that YTHDF1 could serve as independent prognostic factors for HCC through bioinformation analysis [18]. Liu et al. revealed that YTHDF1 accelerated the translational output of FZD5 mRNA and promoted the progression of HCC through the WNT/β-catenin pathway [19]. Zhong et al. found that YTHDF2 could be a tumor suppressor to inhibit cell proliferation and growth in HCC through promoting the degradation of EGFR mRNA [20]. In addition, YTHDF2 has also been reported to promote stem cell phenotype and metastasis of HCC through regulating the m6A methylation level of OCT4 mRNA [21]. Pu et al. indicated that IGF2BP2 could maintain FEN1 expression in an m6A-IGF2BP2-dependent method, and it may be a potential biomarker for the prognostic prediction in HCC [22]. However, no comprehensive analysis based on multiple m6A “readers” in HCC has been integrated so far.

In this study, we conducted a comprehensive study on the expression level and clinical prognostic value of m6A “readers” in HCC. We showed the different expression of m6A “readers” between HCC and normal tissues in mRNA and protein level. Our results showed that the mRNA expression level of the common m6A “readers” including YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP were highly expressed in HCC. And the protein expression level of most of m6A “readers” was the same as mRNA expression. In addition, YTHDC1, YTHDF1, YTHDF2, IGF2BP3, HNRNPA2B1 and NKAP were confirmed by two databases to be positively correlated with tumor stage, which deserved more attention in the treatment of HCC. Besides, our results implied that YTHDC1, YTHDF1, YTHDF2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, and NKAP had a prognosis value in either OS or PFS. But only YTHDF1, IGF2BP2, IGF2BP3 and NKAP had significant value to predict both OS and PFS time in HCC patients. This seems a little inconsistent with the previous results [14]. Due to the researches on YTHDC2 in HCC were really rare, we could not interpret this paradoxical result. More mechanism and clinical researches were needed.

Somatic mutation and clonal selection were two major reasons for cancer development [23, 24]. There were large amounts of evidences showed that genetic alterations play a significant role in tumor occurrence [25, 26]. Researches showed that there are two categories of clinically important gene mutations. The first category is single nucleotide mutations (including synonymous mutations, missense mutations, nonsense mutations and splice mutation), which involves mutations that alter the coding region of a gene, resulting in the loss of function or abnormal activation of the protein product of a gene. The second type of mutation is structural mutations (including deletions, insertions, duplications, inframe mutation, amplifications, and inversions), which may involve the non-coding sequences that make up 98% of the genome [27]. In this study, all of the m6A “readers” shows frequent gene alterations in HCC. The major genetic alteration is the change of mRNA expression level. In addition, some clinically important gene mutations such as nonsense mutations, splice mutation, inframe mutation, amplifications, and deletions, are also accounted for a certain percentage. These results partly explained the genetic alteration of the m6A “readers” probably played a critical role in HCC.

The co-expression genes may give some new thoughts to the development and progression of HCC. For instance, they may be potential targets of m6A “readers”, or be regulators of m6A “readers”. To determine the functions of these m6A “readers” in HCC, we depicted a PPI network based on the co-expression genes of m6A “readers”. The top 10 hub genes seem from three families (HNRNP family, SNRP family, and SR family). In addition, HNRNP, SNRP, and SR proteins were all discovered as regulators of alternative splicing. HNRNP family was a large family of RNA-binding proteins (RBPs) that participated in multiple aspects of nucleic acid metabolism, which including alternative splicing, mRNA stabilization, transcriptional regulation and translational regulation [28]. Small nuclear ribonucleoprotein polypeptide (SNRP) family is precursors of spliceosome [29, 30]. And they have attracted large attention because of their implicated roles in tumorigenesis and tumor development [31–33]. Serine and arginine-rich (SR) family were known as constitutive and alternative splicing regulators. Shuttling SR family had a potential to be regulators for translation in the cytoplasm [34]. Besides, SR and HNRNP family were discovered as alternative splicing regulators. SR and HNRNP family significantly affect the responses to chemotherapy through acting as modulators of drug-induced apoptosis [35]. In conclusion, the function of m6A “readers” on HCC may link with these splicing related proteins, and these families may have a potential to be biomarkers or therapeutic targets.

The GO and KEGG enrichment analysis helped us to understand the functional meanings of genes at a molecular level [36, 37]. And according to our enrichment analysis, the metabolic process and RNA splicing associated pathway were enriched. These results, to a certain extent, reconfirmed the critical role of splicing related mechanisms of m6A “readers” in HCC.

In recent years, immunotherapy becomes a popular research field. It is noteworthy that clinically detected cancers must evade the anti-tumor immune response before they can grow gradually. Tumor immune monitoring is essential for the body to eliminate tumor cells. Studies have shown that tumors can resist immune attack through immune system suppression pathways and immune system exclusion or neglect [38]. The tumor progression and recurrence may be associated with immune cell infiltration, and comprehending the interaction between the tumor and host immune system is essential for the discovery of prognostic biomarkers, the reduction of drug resistance and the development of new therapeutic approaches [39]. Furthermore, in previous researches, many immune cells played a key role in the development of tumors. For instance, dendritic cell immunotherapy has been shown to be effective against tumors beyond the central nervous system [40]. CD8+ T cell can destroy tumor cells by cytotoxic molecules such as granzyme and perforin [41]. The goal of cancer immunotherapy is to promote the activity of cytotoxic T lymphocytes (CTL) within the tumor and to establish effective and durable anti-tumor immunity. Specific dendritic cells will help signal transmission from T cells CD4+ to T cells CD8+ to optimize the size and quality of the CTL response [42, 43]. Furthermore, studies have shown that B cells and B cell receptor-associated kinases, such as BTK, function in the microenvironment of squamous cell carcinoma and pancreatic cancer, so targeting B cells or B cell receptor-associated kinases may have more potent anticancer activity than B cell malignancies [44]. In addition, m6A methylation also was reported to modulate immunity and predict patient outcomes [45, 46]. Our study exhibited the relationship between m6A “readers” and immune cell infiltration in HCC.

We found that the expression of the m6A “readers” is significantly positively linked to the infiltration of CD8+ T cell, CD4+ T cell, Tregs, B cell, neutrophil, monocyte, macrophage, DC, NK cell, and mast cell. It was worthy to pay attention to the correlation between HNRNPC and DC. They had a high correlation value 0.535. In addition, among these immune cells, macrophage might play a core role in m6A “readers”, because it showed positively strong correlation with most m6A “readers”. Macrophage was reported as potent immune effector cells, and it had a functional plasticity which contributed to both antitumor and protumor function in different background [47]. In this study, macrophage acted as a protumor factor and the detailed mechanisms were worthy of digging. Besides, the distinct changes in the microenvironment may induce a series of reversible metabolic programs which might form the basis of the activation of macrophage [48]. Thus, m6A “readers” may exert their functions in the microenvironment by interacting with macrophages. However, it needs more researches to clarify the mechanisms. In addition, m6A “readers” also interact with various immune cells. For instance, more than six immune cells showed positively strong correlation with HNRNPC. These results suggested that m6A “readers” played a key role in promoting tumor progression through modulating tumor immune microenvironment [49].

In further analysis and validation of IGF2BP3 functions in cell model, we found that IGF2BP3 showed high expression in HCC cell lines according to the CCLE database, and knockdown of IGF2BP3 inhibited the migration and sphere formation ability of HCC cells. These results suggested that IGF2BP3 might act as an oncogene and was related to the poor prognosis of HCC patients by enhancing cell migration and stemness ability.

Conclusions and limitations

In this study, we systematically collated the expression data of the common m6A “readers” in HCC, and analyzed their clinical value. In addition, we speculated their possible functions and discussed their correlation with immune cell infiltration. It was worthy to note that YTHDF1, IGF2BP3 and NKAP may be potential therapeutic targets and biomarkers. Because they were expressed differently, and were negatively correlated to OS, PFS, and tumor stage in HCC. More importantly, YTHDF1, IGF2BP3 and NKAP were also highly expressed in pan-cancer, which provided the promise of YTHDF1, IGF2BP3 and NKAP as potential therapeutic targets. From the results of the functional enrichment analysis, the main functions of these three potential therapeutic targets are involved in the regulation of cell cycle or protein localization. In addition, the functions of m6A “readers” in HCC may be related to splicing proteins such as HNRNP, SNRP and SR. And the immune cells especially the macrophage, CD4+ T cell, Treg, B cell, monocyte, and myeloid dendritic cell may play a significant role in the progression of HCC. Besides, immune microenvironment of m6A “readers” was also worthy to be explored. In conclusion, our findings in this study could give some new idea to the development of multiple molecular diagnostic biomarkers and treatment targets for HCC patients to improve their prognosis.

Protein expression patterns in cancers contain crucial diagnostic and prognostic information. However, from the HPA database, we could not obtain the protein expression information of YTHDF1 and YTHDF3 in HCC. In addition, there was no data to present the correlations between m6A “readers” and metastasis of multiple lymph node in HCC, and it may be caused by insufficient sample size or other reasons. Besides, we found that most of the m6A readers were associated with the prognosis of HCC patients, and the m6A “readers” were found to be closely related to the immune microenvironment. So, in the near future, we will compose a prognostic model incorporating the m6A “readers” (especially YTHDF1, IGF2BP3 and NKAP) and immune microenvironment components (especially the macrophage, CD4+ T cell, Treg, B cell, monocyte, and myeloid dendritic cell) to better complement this study. Furthermore, our research was mainly based on public databases, and the results of this study should be further validated by more laboratory experiments.

Materials and Methods

Multiple bioinformatics database resources were accessed to explore the expression levels of m6A “readers” in HCC samples. The website tools (UALCAN, GEPIA2, and cBioPortal databases) were based on the TCGA database, which contains 371 HCC samples and 50 normal liver samples. And TPM (Transcripts per million) data was applied for the mRNA expression level analysis.

UALCAN database

The UALCAN database is an interactive portal for in-depth analysis of gene expression data from the TCGA database. It compares the tumor and normal tissue samples based on pathological stage, lymph nodes metastasis and other clinicopathological characteristics [50]. UALCAN database was available from http://ualcan.path.uab.edu/index.html. We used UALCAN database to evaluate the expression difference between HCC and normal tissue. And the differences in transcriptional expression were compared by Student’s t-test, we considered P < 0.05 as statistically significant. In addition, the figures of pathological stage and lymph nodes metastasis were also derived from this database.

The human protein atlas (HPA) database

The Human Protein Atlas database is an immunohistochemistry-based map of protein expression profiles in normal tissues, cancers and cell lines. Proteome analysis was based on 26941 antibodies targeting 17165 unique proteins. It allows complex queries (including expression profiles, protein classes and chromosome location), so as to provide a resource for pathology-based research [51]. And the Human Protein Atlas was available from http://www.proteinatlas.org. In this study, we used the HPA database to show the different expression of m6A “readers” between HCC and normal tissues.

GEPIA2 database

GEPIA2 was available to analyze gene expression quantification both in the gene level and the transcript level. In addition, it was able to analyze and compare different cancer subtype [52]. GEPIA2 analyzed RNA sequencing data from the GTEx and TCGA database, and it was publicly accessible at http://gepia2.cancer-pku.cn/. In this study, we used GEPIA2 database to show the correlation between m6A “readers” and HCC stage.

Kaplan-Meier plotter

The Kaplan Meier plotter could evaluate the effect of 54k genes (mRNA, miRNA, protein) on survival in 21 cancer types. The correlations between gene expression and survival were computed using the Cox proportional hazards regression. And the sources for the databases contain GEO, EGA, and TCGA [53]. Kaplan-Meier plotter was available from http://kmplot.com/analysis/. We used Kaplan-Meier plotter to analyze the prognostic value including overall survival (OS) and progression free survival (PFS) of the m6A “readers” in HCC patients. We divided 371 HCC specimens into high and low expression groups according to their mRNA expression by automatically selecting the best cut-off value. Kaplan-Meier analysis was evaluated by log-rank test, and we considered P < 0.05 as statistically significant difference.

cBioPortal

cBioPortal database provides a Web resource for exploring, visualizing, and analyzing cancer genomics data. It provides and simplifies the summaries of gene-level data from multiple platforms [54]. cBioPortal was available from https://www.cbioportal.org. In this study, we used cBioPortal to analyze the genome map of the m6A “readers” in HCC and obtain the mutation and mRNA expression data. Besides, we downloaded the co-expression genes of m6A “readers” from cBioPortal for further analysis.

STRING

The STRING database collects and integrates publicly available sources of protein-protein interaction information, and provides the computational predictions. And it aims to achieve a comprehensive and objective global network among proteins [55]. STRING was available from https://string-db.org/. In our study, STRING was applied to explore and analyze possible protein-protein interactions (PPIs).

Cytoscape

Cytoscape is a popular tool for analyzing and visualizing multiple proteins or genes and for integrating biomolecular interaction networks into a unified conceptual framework [56]. In our study, we used Cytoscape to conduct functional integration on 220 co-expression genes of m6A “readers” screened from cBioPortal. These 220 co-expression genes integrated from top 20 genes that had a high spearman’s correlation with m6A “readers”. According to the spearman’s correlation value between the interacting proteins, we made genes with high spearman’s correlation value (>0.6) pink, and those genes with low spearman’s correlation value (<0.6) blue. Besides, we used it to figure out top 10 hub genes of these co-expression genes.

WebGestalt

WebGestalt is a tool for the interpretation of gene lists. It supports three well-established and complementary methods for enrichment analysis, including over-representation analysis (ORA), gene set enrichment analysis (GSEA), and network topology-based analysis (NTA) [57]. WebGestalt was available from http://www.webgestalt.org/. In our study, we used WebGestalt to provide enrichment analysis results including Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of co-expression genes of m6A “readers”.

Metascape

Metascape, a gene annotation and analysis resource, integrated functional enrichment, interactome analysis, gene annotation, and membership together [58]. Metascape was available from https://metascape.org/gp. In this study, we used Metascape to analyze the functional enrichment of m6A “readers” and their co-expression genes.

TIMER2.0

TIMER is a comprehensive resource for the analysis of immune infiltrates among various cancer types. TIMER2.0 provides immune infiltrates’ abundances estimated by multiple immune deconvolution methods [59]. TIMER2.0 was available from http://timer.cistrome.org/. In this study, we used TIMER2.0 to depict the relationship between m6A “readers” and immune cells including CD8+ T cell, CD4+ T cell, T cell regulatory (Tregs), B cell, neutrophil, monocyte, macrophage, myeloid dendritic cell (DC), nature killer (NK) cell, mast cell, and cancer associated fibroblast (CAF). All of these data had a purity adjustment process. And log2 RSEM was used for displaying the gene expression levels.

CCLE

Cancer Cell Line Encyclopedia (CCLE) database can be used to explore gene variants, candidate targets, small molecules and biotherapeutics, and characterize cancer cell lines [60]. In this study, the mRNA expression level of HCC cell lines was obtained from the CCLE dataset (https://portals.broadinstitute.org/ccle).

Cell culture

Hep3B cells were purchased from Zhong Qiao Xin Zhou Biotechnology Co., Ltd. (Shanghai, China) and were cultured in MEM medium (Gibco, USA) containing 10% fetal bovine serum (FBS) (Gibco, USA), 100 μg/mL streptomycin and 100 IU/mL penicillin (Gibco, USA). They were cultured at 37°C, 5% CO2.

Lentiviral transduction

Stable KD of IGF2BP3 was obtained by a lentiviral-based shRNA system. Target specific shRNAs or nontarget control (NC) were subcloned into GV493(hU6-MCS-CBh-gcGFP-IRES-puromycin) vector. To isolate infected Hep3B cells, Puromycin selection was performed at 2 ug/ml for 7 days. The KD efficiencies were verified at both mRNA and protein levels. The shRNA sequences used in this study are as follows. shIGF2BP3#1, 5′-ATAGGTTACATTTACAACTGC-3′, shIGF2BP3#2, 5′-TAATCCAGGAATTAAATGTGC-3′, shNC, 5′-TTCTCCGAACGTGTCACGT-3′. 1 μg/ml puromycin for selection was added to the medium of cells transduced with lentivirus.

RT-PCR

Total RNA was extracted from Hep3B cells which transfected with shNC and shIGF2BP3 lentivirus, according to the RNA Isolater (Vazyme, China). The concentration and purity of total RNA were measured. Through SureScript First-strand cDNA synthesis kit (GeneCopoeia, USA), total RNA was used to obtain cDNA by reverse transcription reaction. The RT-PCR assays were carried out by a protocol from Power SYBR Green (Takara, Hangzhou, Zhejiang, China). The relative expressions of genes were calculated and normalized using the 2−ΔΔCt methods relative to GAPDH. Specific primer sequences were as follows: IGF2BP3: 5′-ACGAAATATCCCGCCTCATTTAC-3′ (forward), 5′-GCAGTTTCCGAGTCAGTGTTCA-3′ (reverse); GAPDH: 5′-CTGGGCTACACTGAGCACC-3′ (forward), 5′-AAGTGGTCGTTGAGGGCAATG-3′ (reverse).

Western blot

Whole cell lysates were collected in RIPA lysis buffer (Beyotime, China). The antibodies used in this study are IGF2BP3 (#ab179807, Abcam), and β-actin (#20536-1-AP, Proteintech) antibodies. HRP-conjugated anti-rabbit secondary antibody was used for the ECL detection.

Migration assay

800 μl of conditioned medium containing 20% fetal bovine serum was used as a chemotactic agent for the bottom chamber of the 24-well plate. Serum-free medium of cells transfected with lentivirus shNC and shIGF2BP3 were prepared as suspension at 5 × 105 cells/ml and 200 μl of the suspension was inoculated in transwell chambers. After 48 hrs of incubation, the cells were fixed in methanol for 30 mins and stained with crystal violet for 30 mins.

Sphere formation assay

Lentiviral transfected shNC and shIGF2BP3 cells were digested with trypsin into single cell suspensions and seeded into special 24-well plates (with 1% polyHEMA (Sigma)-coated) at a density of 3000 cells per well. Cells were grown for 10–14 days in special medium (which including serum-free DMEM/F12 medium (Gibco, USA) supplemented with 4mg/mL insulin (Sigma-Aldrich, St. Louis, MO, USA), B27 (Gemini, USA), 20 ng/mL EGF (MCE, USA), and 20 ng/mL FGF (MCE, USA). The numbers of tumor spheres larger than 50 μm in diameter were calculated.

Statistical analysis

The differences between two groups were analyzed by Student’s t-tests, while differences among more than two groups were analyzed by ANOVA. The two patient cohorts are compared by a Kaplan-Meier survival plot, and the hazard ratio with 95% confidence intervals and log-rank P value are calculated. Pearson and Spearman correlation analysis were used to analyze correlations of the data. The P < 0.05 was considered as statistical significance.

All methods were performed in accordance with the relevant guidelines and regulations.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.

Abbreviations

BP: biological process; CAF: cancer associated fibroblast; CC: cellular component; CTL: cytotoxic T lymphocytes; DC: dendritic cell; GO: Gene ontology; GSEA: gene set enrichment analysis; HPA: The Human Protein Atlas; KEGG: Kyoto Encyclopedia of Genes and Genomes; HCC: hepatocellular carcinoma; m6A: N6-methyladenosine; MF: molecular function; NK: nature killer; NTA: network topology-based analysis; ORA: over-representation analysis; OS: overall survival; PFS: progression free survival; PPI: protein-protein interactions; RBPs: RNA-binding proteins; SNRP: small nuclear ribonucleoprotein polypeptide; SR: serine and arginine-rich; Tregs: T cell regulatory.

Author Contributions

LZ designed the study. SQ drafted the manuscript. SQ, GL, HJ, XC, JH, JX, YQ, YM and LZ collected the data and conducted the picture processing. LZ and YM revised the manuscript. All authors have read and approved the final version of manuscript.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

This study was supported by the National Natural Science Foundation of China (81602575, 81701847, 61572200), Natural Science Foundation of Hunan Province (2022JJ30794), Xiang Jiang Scholars Program/Hong Kong Scholars Program (XJ2016054), Basic Research Project of Changsha Science and Technology Plan (kq2004126), Changsha Municipal Natural Science Foundation (kq2202126), Graduate Education and Teaching Reform Project of Central South University (2022YJSKS018), the Student Innovation Project of Central South University (2022ZZTS0267), University Student Innovation and Entrepreneurship Training Program of Hunan Province (S2021105330529).

References

  • 1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021; 71:209–49. https://doi.org/10.3322/caac.21660 [PubMed]
  • 2. Bucci L, Garuti F, Lenzi B, Pecorelli A, Farinati F, Giannini EG, Granito A, Ciccarese F, Rapaccini GL, Di Marco M, Caturelli E, Zoli M, Borzio F, et al, and Italian Liver Cancer (ITA.LI.CA) group. The evolutionary scenario of hepatocellular carcinoma in Italy: an update. Liver Int. 2017; 37:259–70. https://doi.org/10.1111/liv.13204 [PubMed]
  • 3. Gryziak M, Woźniak K, Kraj L, Stec R. Milestones in the treatment of hepatocellular carcinoma: A systematic review. Crit Rev Oncol Hematol. 2021; 157:103179. https://doi.org/10.1016/j.critrevonc.2020.103179 [PubMed]
  • 4. Pinter M, Jain RK, Duda DG. The Current Landscape of Immune Checkpoint Blockade in Hepatocellular Carcinoma: A Review. JAMA Oncol. 2021; 7:113–23. https://doi.org/10.1001/jamaoncol.2020.3381 [PubMed]
  • 5. Granito A, Muratori L, Lalanne C, Quarneti C, Ferri S, Guidi M, Lenzi M, Muratori P. Hepatocellular carcinoma in viral and autoimmune liver diseases: Role of CD4+ CD25+ Foxp3+ regulatory T cells in the immune microenvironment. World J Gastroenterol. 2021; 27:2994–3009. https://doi.org/10.3748/wjg.v27.i22.2994 [PubMed]
  • 6. Nault JC, Villanueva A. Biomarkers for Hepatobiliary Cancers. Hepatology. 2021 (Suppl 1); 73:115–27. https://doi.org/10.1002/hep.31175 [PubMed]
  • 7. Qin S, Mao Y, Chen X, Xiao J, Qin Y, Zhao L. The functional roles, cross-talk and clinical implications of m6A modification and circRNA in hepatocellular carcinoma. Int J Biol Sci. 2021; 17:3059–79. https://doi.org/10.7150/ijbs.62767 [PubMed]
  • 8. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell. 2012; 149:1635–46. https://doi.org/10.1016/j.cell.2012.05.003 [PubMed]
  • 9. Liu ZX, Li LM, Sun HL, Liu SM. Link Between m6A Modification and Cancers. Front Bioeng Biotechnol. 2018; 6:89. https://doi.org/10.3389/fbioe.2018.00089 [PubMed]
  • 10. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019; 112:108613. https://doi.org/10.1016/j.biopha.2019.108613 [PubMed]
  • 11. Chen XY, Zhang J, Zhu JS. The role of m6A RNA methylation in human cancer. Mol Cancer. 2019; 18:103. https://doi.org/10.1186/s12943-019-1033-z [PubMed]
  • 12. Luo X, Cao M, Gao F, He X. YTHDF1 promotes hepatocellular carcinoma progression via activating PI3K/AKT/mTOR signaling pathway and inducing epithelial-mesenchymal transition. Exp Hematol Oncol. 2021; 10:35. https://doi.org/10.1186/s40164-021-00227-0 [PubMed]
  • 13. Zhao X, Chen Y, Mao Q, Jiang X, Jiang W, Chen J, Xu W, Zhong L, Sun X. Overexpression of YTHDF1 is associated with poor prognosis in patients with hepatocellular carcinoma. Cancer Biomark. 2018; 21:859–68. https://doi.org/10.3233/CBM-170791 [PubMed]
  • 14. Liu J, Wang D, Zhou J, Wang L, Zhang N, Zhou L, Zeng J, Liu J, Yang M. N6-methyladenosine reader YTHDC2 and eraser FTO may determine hepatocellular carcinoma prognoses after transarterial chemoembolization. Arch Toxicol. 2021; 95:1621–9. https://doi.org/10.1007/s00204-021-03021-3 [PubMed]
  • 15. van Dam S, Võsa U, van der Graaf A, Franke L, de Magalhães JP. Gene co-expression analysis for functional classification and gene-disease predictions. Brief Bioinform. 2018; 19:575–92. https://doi.org/10.1093/bib/bbw139 [PubMed]
  • 16. Kurotaki D, Sasaki H, Tamura T. Transcriptional control of monocyte and macrophage development. Int Immunol. 2017; 29:97–107. https://doi.org/10.1093/intimm/dxx016 [PubMed]
  • 17. Tanabe A, Konno J, Tanikawa K, Sahara H. Transcriptional machinery of TNF-α-inducible YTH domain containing 2 (YTHDC2) gene. Gene. 2014; 535:24–32. https://doi.org/10.1016/j.gene.2013.11.005 [PubMed]
  • 18. Zhao Z, Yang L, Fang S, Zheng L, Wu F, Chen W, Song J, Chen M, Ji J. The Effect of m6A Methylation Regulatory Factors on the Malignant Progression and Clinical Prognosis of Hepatocellular Carcinoma. Front Oncol. 2020; 10:1435. https://doi.org/10.3389/fonc.2020.01435 [PubMed]
  • 19. Liu X, Qin J, Gao T, Li C, He B, Pan B, Xu X, Chen X, Zeng K, Xu M, Zhu C, Pan Y, Sun H, et al. YTHDF1 Facilitates the Progression of Hepatocellular Carcinoma by Promoting FZD5 mRNA Translation in an m6A-Dependent Manner. Mol Ther Nucleic Acids. 2020; 22:750–65. https://doi.org/10.1016/j.omtn.2020.09.036 [PubMed]
  • 20. Zhong L, Liao D, Zhang M, Zeng C, Li X, Zhang R, Ma H, Kang T. YTHDF2 suppresses cell proliferation and growth via destabilizing the EGFR mRNA in hepatocellular carcinoma. Cancer Lett. 2019; 442:252–61. https://doi.org/10.1016/j.canlet.2018.11.006 [PubMed]
  • 21. Zhang C, Huang S, Zhuang H, Ruan S, Zhou Z, Huang K, Ji F, Ma Z, Hou B, He X. YTHDF2 promotes the liver cancer stem cell phenotype and cancer metastasis by regulating OCT4 expression via m6A RNA methylation. Oncogene. 2020; 39:4507–18. https://doi.org/10.1038/s41388-020-1303-7 [PubMed]
  • 22. Pu J, Wang J, Qin Z, Wang A, Zhang Y, Wu X, Wu Y, Li W, Xu Z, Lu Y, Tang Q, Wei H. IGF2BP2 Promotes Liver Cancer Growth Through an m6A-FEN1-Dependent Mechanism. Front Oncol. 2020; 10:578816. https://doi.org/10.3389/fonc.2020.578816 [PubMed]
  • 23. Martincorena I, Raine KM, Gerstung M, Dawson KJ, Haase K, Van Loo P, Davies H, Stratton MR, Campbell PJ. Universal Patterns of Selection in Cancer and Somatic Tissues. Cell. 2017; 171:1029–41.e21. https://doi.org/10.1016/j.cell.2017.09.042 [PubMed]
  • 24. Sondka Z, Bamford S, Cole CG, Ward SA, Dunham I, Forbes SA. The COSMIC Cancer Gene Census: describing genetic dysfunction across all human cancers. Nat Rev Cancer. 2018; 18:696–705. https://doi.org/10.1038/s41568-018-0060-1 [PubMed]
  • 25. Sharma A, Jiang C, De S. Dissecting the sources of gene expression variation in a pan-cancer analysis identifies novel regulatory mutations. Nucleic Acids Res. 2018; 46:4370–81. https://doi.org/10.1093/nar/gky271 [PubMed]
  • 26. Supek F, Miñana B, Valcárcel J, Gabaldón T, Lehner B. Synonymous mutations frequently act as driver mutations in human cancers. Cell. 2014; 156:1324–35. https://doi.org/10.1016/j.cell.2014.01.051 [PubMed]
  • 27. Manolio TA, Rowley R, Williams MS, Roden D, Ginsburg GS, Bult C, Chisholm RL, Deverka PA, McLeod HL, Mensah GA, Relling MV, Rodriguez LL, Tamburro C, Green ED. Opportunities, resources, and techniques for implementing genomics in clinical care. Lancet. 2019; 394:511–20. https://doi.org/10.1016/S0140-6736(19)31140-7 [PubMed]
  • 28. Geuens T, Bouhy D, Timmerman V. The hnRNP family: insights into their role in health and disease. Hum Genet. 2016; 135:851–67. https://doi.org/10.1007/s00439-016-1683-5 [PubMed]
  • 29. Mabonga L, Kappo AP. The oncogenic potential of small nuclear ribonucleoprotein polypeptide G: a comprehensive and perspective view. Am J Transl Res. 2019; 11:6702–16. [PubMed]
  • 30. Chen T, Zhang B, Ziegenhals T, Prusty AB, Fröhler S, Grimm C, Hu Y, Schaefke B, Fang L, Zhang M, Kraemer N, Kaindl AM, Fischer U, Chen W. A missense mutation in SNRPE linked to non-syndromal microcephaly interferes with U snRNP assembly and pre-mRNA splicing. PLoS Genet. 2019; 15:e1008460. https://doi.org/10.1371/journal.pgen.1008460 [PubMed]
  • 31. Lan Y, Lou J, Hu J, Yu Z, Lyu W, Zhang B. Downregulation of SNRPG induces cell cycle arrest and sensitizes human glioblastoma cells to temozolomide by targeting Myc through a p53-dependent signaling pathway. Cancer Biol Med. 2020; 17:112–31. https://doi.org/10.20892/j.issn.2095-3941.2019.0164 [PubMed]
  • 32. Anchi T, Tamura K, Furihata M, Satake H, Sakoda H, Kawada C, Kamei M, Shimamoto T, Fukuhara H, Fukata S, Ashida S, Karashima T, Yamasaki I, et al. SNRPE is involved in cell proliferation and progression of high-grade prostate cancer through the regulation of androgen receptor expression. Oncol Lett. 2012; 3:264–8. https://doi.org/10.3892/ol.2011.505 [PubMed]
  • 33. Dai X, Yu L, Chen X, Zhang J. SNRPD1 confers diagnostic and therapeutic values on breast cancers through cell cycle regulation. Cancer Cell Int. 2021; 21:229. https://doi.org/10.1186/s12935-021-01932-w [PubMed]
  • 34. Jeong S. SR Proteins: Binders, Regulators, and Connectors of RNA. Mol Cells. 2017; 40:1–9. https://doi.org/10.14348/molcells.2017.2319 [PubMed]
  • 35. Kędzierska H, Piekiełko-Witkowska A. Splicing factors of SR and hnRNP families as regulators of apoptosis in cancer. Cancer Lett. 2017; 396:53–65. https://doi.org/10.1016/j.canlet.2017.03.013 [PubMed]
  • 36. The Gene Ontology Consortium. The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res. 2019; 47:D330–8. https://doi.org/10.1093/nar/gky1055 [PubMed]
  • 37. Kanehisa M, Sato Y, Furumichi M, Morishima K, Tanabe M. New approach for understanding genome variations in KEGG. Nucleic Acids Res. 2019; 47:D590–5. https://doi.org/10.1093/nar/gky962 [PubMed]
  • 38. Gajewski TF, Schreiber H, Fu YX. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol. 2013; 14:1014–22. https://doi.org/10.1038/ni.2703 [PubMed]
  • 39. Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, Jiang P, Shen H, Aster JC, Rodig S, Signoretti S, Liu JS, Liu XS. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 2016; 17:174. https://doi.org/10.1186/s13059-016-1028-7 [PubMed]
  • 40. Sokratous G, Polyzoidis S, Ashkan K. Immune infiltration of tumor microenvironment following immunotherapy for glioblastoma multiforme. Hum Vaccin Immunother. 2017; 13:2575–82. https://doi.org/10.1080/21645515.2017.1303582 [PubMed]
  • 41. Tsukumo SI, Yasutomo K. Regulation of CD8+ T Cells and Antitumor Immunity by Notch Signaling. Front Immunol. 2018; 9:101. https://doi.org/10.3389/fimmu.2018.00101 [PubMed]
  • 42. Borst J, Ahrends T, Bąbała N, Melief CJM, Kastenmüller W. CD4+ T cell help in cancer immunology and immunotherapy. Nat Rev Immunol. 2018; 18:635–47. https://doi.org/10.1038/s41577-018-0044-0 [PubMed]
  • 43. Lee YS, Radford KJ. The role of dendritic cells in cancer. Int Rev Cell Mol Biol. 2019; 348:123–78. https://doi.org/10.1016/bs.ircmb.2019.07.006 [PubMed]
  • 44. Burger JA, Wiestner A. Targeting B cell receptor signalling in cancer: preclinical and clinical advances. Nat Rev Cancer. 2018; 18:148–67. https://doi.org/10.1038/nrc.2017.121 [PubMed]
  • 45. Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, Huang X, Liu Y, Wang J, Dougherty U, Bissonnette MB, Shen B, Weichselbaum RR, et al. Anti-tumour immunity controlled through mRNA m6A methylation and YTHDF1 in dendritic cells. Nature. 2019; 566:270–4. https://doi.org/10.1038/s41586-019-0916-x [PubMed]
  • 46. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020; 19:53. https://doi.org/10.1186/s12943-020-01170-0 [PubMed]
  • 47. Anderson NR, Minutolo NG, Gill S, Klichinsky M. Macrophage-Based Approaches for Cancer Immunotherapy. Cancer Res. 2021; 81:1201–8. https://doi.org/10.1158/0008-5472.CAN-20-2990 [PubMed]
  • 48. El Kasmi KC, Stenmark KR. Contribution of metabolic reprogramming to macrophage plasticity and function. Semin Immunol. 2015; 27:267–75. https://doi.org/10.1016/j.smim.2015.09.001 [PubMed]
  • 49. Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, Li J, Li F, Tan HB. Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett. 2020; 470:126–33. https://doi.org/10.1016/j.canlet.2019.11.009 [PubMed]
  • 50. Chandrashekar DS, Bashel B, Balasubramanya SAH, Creighton CJ, Ponce-Rodriguez I, Chakravarthi BVS, Varambally S. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia. 2017; 19:649–58. https://doi.org/10.1016/j.neo.2017.05.002 [PubMed]
  • 51. Pontén F, Jirström K, Uhlen M. The Human Protein Atlas--a tool for pathology. J Pathol. 2008; 216:387–93. https://doi.org/10.1002/path.2440 [PubMed]
  • 52. Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019; 47:W556–60. https://doi.org/10.1093/nar/gkz430 [PubMed]
  • 53. Nagy Á, Munkácsy G, Győrffy B. Pancancer survival analysis of cancer hallmark genes. Sci Rep. 2021; 11:6047. https://doi.org/10.1038/s41598-021-84787-5 [PubMed]
  • 54. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E, Cerami E, Sander C, Schultz N. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013; 6:pl1. https://doi.org/10.1126/scisignal.2004088 [PubMed]
  • 55. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 47:D607–13. https://doi.org/10.1093/nar/gky1131 [PubMed]
  • 56. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
  • 57. Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 2019; 47:W199–205. https://doi.org/10.1093/nar/gkz401 [PubMed]
  • 58. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019; 10:1523. https://doi.org/10.1038/s41467-019-09234-6 [PubMed]
  • 59. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020; 48:W509–14. https://doi.org/10.1093/nar/gkaa407 [PubMed]
  • 60. Ghandi M, Huang FW, Jané-Valbuena J, Kryukov GV, Lo CC, McDonald ER 3rd, Barretina J, Gelfand ET, Bielski CM, Li H, Hu K, Andreev-Drakhlin AY, Kim J, et al. Next-generation characterization of the Cancer Cell Line Encyclopedia. Nature. 2019; 569:503–8. https://doi.org/10.1038/s41586-019-1186-3 [PubMed]
Download PDF
Corresponding Authors

Luqing Zhao
luqingzhao@csu.edu.cn

Yitao Mao
maoyt@csu.edu.cn

Keywords
N6-methyladenosine (m6A) modification hepatocellular carcinoma (HCC) functional analysis expression profiles immune cell infiltration
Paper Sections