Analysis and Validation of circRNA-miRNA Network in Regulating m6A RNA Methylation Modulators Reveals CircMAP2K4/miR-139-5p/YTHDF1 Axis Involving the Proliferation of Hepatocellular Carcinoma

The m6A RNA methylation modulators play a crucial role in regulating hepatocellular carcinoma (HCC) progression. The circular RNA (circRNA) regulatory network in regulating m6A RNA methylation modulators in HCC remains largely unknown. In this study, 5 prognostic m6A RNA methylation modulators in HCC were identified from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) projects. The differentially expressed microRNAs (DEmiRNAs) and circRNAs (DEcircRNAs) between paired tumor and normal tissues were screened out from TCGA and or Gene Expression Omnibus (GEO) database to construct the circRNA-miRNA- m6A RNA methylation modulator regulatory network, which included three m6A RNA methylation modulators (HNRNPC, YTHDF1, and YTHDF2), 11 DEmiRNAs, and eight DEcircRNAs. Among the network, hsa-miR-139-5p expression was negatively correlated with YTHDF1. Hsa-miR-139-5p low or YTHDF1 high expression was correlated with high pathological grade, advanced stage and poor survival of HCC. Additionally, cell cycle, base excision repair, and homologous recombination were enriched in YTHDF1 high expression group by GSEA. A hub circRNA regulatory network was constructed based on hsa-miR-139-5p/YTHDF1 axis. Furthermore, hsa_circ_0007456(circMAP2K4) was validated to promote HCC cell proliferation by binding with hsa-miR-139-5p to promote YTHDF1 expression. Taken together, we identified certain circRNA regulatory network related to m6A RNA methylation modulators and provided clues for mechanism study and therapeutic targets for HCC.


INTRODUCTION
Hepatocellular carcinoma (HCC) is one of the most common malignant tumors worldwide, accounting for 75%-85% of primary liver cancer (1). Despite advances in diagnosis and therapeutic strategies in recent years, the prognosis of HCC is still not ideal. Therefore, there exists an urgent need to identify sensitive and specific biomarkers and therapeutic targets for the early diagnosis and treatment of HCC (2).
Among the chemical modification of RNA, N6methyladenosine (m 6 A) methylated at the N6 position of adenosine is viewed as the most common, abundant and conservative internal transcriptional modification for various kinds of RNA. The m 6 A modifications are involved in RNA processing, transporting, translation and metabolism (3). Based on the different functions of m 6 A RNA methylation modulators, they are usually classified into "writers", "erasers", and "readers" (4). The "writers" catalyze the formation of m 6 A, including Methyltransferase-like 3 (METTL3) (5), METTL14 (6), Wilms tumor 1-associated protein (WTAP) (7), RNA Binding Motif Protein 15 (RBM15) (8) and KIAA1429 (9). The "erasers", removing m 6 A modification from RNA, compose of fat mass and obesity-associated protein (FTO) (10) and alkB homologue 5 (ALKBH5) (11). The m 6 A readers YT521-B homology (YTH) domain-containing proteins (YTHDF1/2 and YTHDC1/2) function as m 6 A binding proteins that recognize m 6 A methylation and generate a functional signal (12). Accumulating evidence has demonstrated that m 6 A modifications participate in the progression of cancers, such as glioma, breast cancers and hepatocellular carcinoma (HCC) (13).
Circular RNA (circRNA) is a class of covalently closed singlestranded circular RNA molecules formed by back-splicing. CircRNA is considered as the RNA with tissue-, developmental stage-and disease-specificity (14). Notably, it is well documented that circRNAs play crucial roles in cancer by acting as microRNA (miRNA) sponge to modulate the miRNA-mRNA regulatory axis, thereby affecting the initiation and progression of cancer (15). However, whether circRNAs can serve as miRNA sponge to affect the miRNA in the regulation of m 6 A RNA methylation modulators in HCC has not been reported yet.
The flowchart of this study design was shown in Figure S1. Briefly, we identified prognostic m 6 A RNA methylation modulators in HCC patients from The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) projects. According to the differentially expressed (DE) circRNAs and miRNAs between paired normal and tumor tissues, the circRNA-miRNA-mRNA regulatory network was constructed based on the m 6 A RNA methylation modulators. Among the 3 co-expressed miRNA-m 6 A RNA methylation modulators pairs, hsa-miR-139-5p low or YTHDF1 high expression was significantly correlated with high pathological grade, advanced stage and poor survival of HCC. Therefore, a hub circRNA regulatory network was constructed based on hsa-miR-139-5p/YTHDF1 axis. Among this hub network, circMAP2K4 was validated to promote HCC cell proliferation by binding with hsa-miR-139-5p to promote YTHDF1 expression. These findings indicate certain circRNA regulatory network is involved in the regulation of m 6 A RNA methylation modulators and provide clues for mechanism study and therapeutic strategy development for HCC.

Data Collection
Regarding the expression data of m6A RNA methylation modulators, we obtained transcriptome data of TCGA-LIHC project from TCGA data portal (https://tcga-data.nci.nih.gov/ tcga/) and ICGC-LIRI-JP project from ICGC data portal (https:// dcc.icgc.org/), respectively. Regarding the miRNA data, in order to maintain the consistency of data sources, we downloaded miRNA-seq data from the TCGA-LIHC project for subsequent difference and co-expression analysis. For the verification of the prognostic value of miRNA, we searched the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/gds/) with "microRNA", "hepatocellular carcinoma", and "survival" as keywords. In order to ensure the reliability of the results, we only selected datasets with more than 100 cases for analysis, and finally included the GSE31384 into this study. Regarding the circRNA data, we searched the GEO database with "circular RNA", "hepatocellular carcinoma", and "microarray" as keywords, and finally included GSE94508, GSE97332 and GSE78520 into this study. Criteria for study inclusion were: 1) The disease was diagnosed as HCC. 2) HCC caused by different etiologies was acceptable. 3) The case had a complete expression profile. 4) The case had clinical information. Criteria for study exclusion were: 1) The survival data was unknown or survival time was less than 30 days. 2) The clinical staging and or pathological grade was unknown. The analysis flowchart of HCC cases with complete expression data was shown in Figure S2.

Identification of DERNAs
The mRNA expression level of 13 m 6 A RNA methylation modulators were compared between 50 or 199 paired tumor and non-tumor samples from TCGA and ICGC projects by Mann-Whitney-Wilcoxon Test, respectively. And the prognostic value of the m 6 A RNA methylation modulators in both TCGA and ICGC were further assessed by univariate Cox regression survival analyses. Finally, those prognostic m 6 A RNA methylation modulators in both TCGA and ICGC were identified for the following network construction. The DEmiRNAs were screened out from 49 paired tumor and nontumor samples from TCGA by using R package "Bioconductor Limma". The adjusted P value (false discovery rate, FDR) of each gene was calculated by Benjamini Hochberg method and the threshold for DEmiRNA selection was FDR <0.05 and | log2FC | > 1. Finally, DEmiRNA is visualized by volcano graph and fold change (FC) filtering. According to the significance scores <0.01 and | log2FC |> 2, the DEcircRNAs between tumor and nontumor cases from multiple studies was determined by a robust rank aggregation method (16). And the DEcircRNAs were visualized by heatmap.

Construction of CircRNA-miRNA Network Involved in Regulating m 6 A RNA Methylation Modulators
The miRNAs potentially targeting m 6 A RNA methylation regulators were predicted by microRNA Data Integration Portal (miRDIP), which integrated more than 20 miRNA related databases for miRNA target or miRNA prediction (17). Among the DEmiRNAs, the potential miRNAs targeting m 6 A RNA methylation regulators were selected with the very high score (top 1%) in miRDIP. Next, the DEcircRNAs targeted miRNAs were predicted by Cancer-Specific CircRNA Database (CSCD, https:// http://gb.whu.edu.cn/CSCD/). As the sponge of miRNA, the expression level of circRNA usually does not influence the expression of miRNA. In addition, some miRNAs may inhibit highly expressed mRNA in a compensatory elevated expression manner (18). Therefore, the selection of circRNA-miRNA or miRNA-mRNA pairs was not limited by their expression patterns that must be reversed. Finally, the circRNA-miRNA-mRNA regulatory network was constructed after taking the intersection of DEcircRNA-miRNA pairs and DEmiRNA-m 6 A RNA methylation regulator pairs. The regulatory network was visualized using Cytoscape 3.4.0 (http://cytoscape.org/).

Gene Set Enrichment Analysis (GSEA)
The HCC samples from TCGA or ICGC were divided into highand low-expression groups according to the expression level of YTHDF1, respectively. GSEA (http://software.broadinstitute. org/gsea/index.jsp) was carried out to compare the potential biological pathways between two groups. The annotated gene set list c2.cp.kegg. v5.2.symbols.gmt was utilized as the reference gene set. The cut-off criteria were defined as FDR < 0.25 and a nominal P < 0.01. The gene sets with top 5 normalized enrichment score (NES) in high-and low-expression groups were selected for visualization.

Cell Culture and Transfection
Human HCC cell lines Huh7, Hep3B, MHCC97H, HCCLM3, and normal LO2 cells were gained from Shanghai Advanced Research Institute, Chinese Academy of Sciences. Cells were cultured in Dulbecco's Modified Eagle's Medium (DMEM) (pH 7.4) supplemented with 10% (v/v) fetal bovine serum (Gibco). YTHDF1 siRNA, hsa-miR-139-5p mimics, circRNA overexpressing plasmid and their corresponding negative control were purchased from GenePharma (Shanghai, China). The cells in 24-well plates were transfected with 1ug plasmid, 50 nM mimics or siRNA by using Lipofectamine 3000 reagent (Invitrogen) according to the manufacturer's recommendation. The specific siRNA sequences for YTHDF1 were provided in the Table S1. Three independent experiments were carried out for cell transfection.

Quantitative Reverse Transcription Polymerase Reaction (qRT-PCR) Analysis
Total RNA were isolated from cells using TRIzol reagent (Invitrogen) and cDNA were synthesized by utilizing the Prime Script RT reagent kit (Takara Bio, Shiga, Japan). The SYBR ® Premix Ex Taq ™ (Takara) were used for qRT-PCR detection through real-time detection system (ABI7500, USA). The primer sequences for detection were provided in Table S2. GAPDH was used as an internal standard control. Gene expression level was quantified using 2 -△△Ct method. The results were obtained from three independent experiments.

Western Blot
Protein was extracted using RIPA (Beyotime, China) and separated on 10% SDS-PAGE gels and transferred onto polyvinylidene fluoride membranes (Millipore, USA). The primary antibodies of anti-YTHDF1 (#86463) and anti-GAPDH (#2118) were purchased from Cell Signaling Technology (Danvers, MA, USA). After incubating with the primary antibodies at 4°C overnight, the membranes were then subjected to HRP-conjugated secondary antibody (Cell Signaling Technology, USA) at room temperature for 1 h. The blots were visualized using an imaging system (Bio-Rad, USA).

Agarose Gel Electrophoresis Analysis
Total RNA isolation and cDNA synthesis were mentioned as above. Total DNA was extracted using the SteadyPure Universal Genomic DNA Extraction Kit (Accurate Biology Co. Ltd., Changsha, China). The specific divergent primers and convergent primers for circMAP2K4 were used for amplification. The primer sequences were listed in Table S2. Then the amplification products were detected in 2% agarose gel electrophoresis. And those products amplified by divergent primers were used for Sanger sequencing to detect the junction site of circMAP2K4. Three independent experiments were carried out for agarose gel electrophoresis analysis.

Cell Proliferation Assay
The transfected cells were seeded in 96-well plates at a density of 2,000 cells per well. Cell viability was accessed from 12 to 120 h by using the Cell Counting Kit-8 (CCK-8) according to the manufacturer's recommendation (Dojindo, Kumamoto, Japan). The optical density (OD) was recorded at 450 nm by an automatic microplate reader (Synergy4; BioTek, Winooski, VT, USA). The results were obtained from three independent experiments.

Luciferase Reporter Assay
The wild type or mutated circRNA sequence containing hsa-miR-139-5p binding site, the wild type or the mutated 3' untranslated region (UTR) of YTHDF1, were respectively synthesized and inserted into pmiR-RB-REPORT ™ vector (RIBOBIO, Guangzhou, China). The above vectors and hsa-miR-139-5p mimics or negative control were co-transfected into cells using Lipofectamine3000 (Invitrogen). 48 h after transfection, the cells were harvested for firefly and renilla luciferase activities detection by using the dual-luciferase reporter assay system (Promega, Massachusetts, USA). Renilla luciferase served as the internal control for luciferase activity. The results were obtained from three independent experiments.

Statistical Analysis
OS differences between high and low expression groups were evaluated by Kaplan-Meier survival analysis and log-rank test.
The differences of gene expression between each clinicopathological characteristics were evaluated by Mann-Whitney-Wilcoxon Test. Data differences between in-vitro experimental groups were analyzed by Student's t-test or one-way analysis of variance (ANOVA). All tests were analyzed using R software version 3.4.2 and P < 0.05 was considered statistically significant.

Most m 6 A RNA Methylation Modulators Are Up-Regulated and Correlated With the Prognosis of HCC
In TCGA project, all m 6 A RNA methylation modulators were up-regulated in HCC, but the expression difference of METL14 and ZC3H13 were not statistically significant ( Figure 1A). In ICGC project, most m 6 A RNA methylation modulators except ZC3H13 were up-regulated in HCC, but the up-regulation of METTL14 was not statistically different ( Figure 1B). Next, we analyzed the prognostic values of 11 commonly DE m 6 A RNA methylation modulators in both TCGA and ICGC. High expression of YTHDF2, YTHDF1, KIAA1429, HNRNPC, WTAP, METTL3, or RBM15 was correlated with the poor survival of HCC in TCGA project by univariate Cox regression survival analysis ( Figure 1C). While in ICGC project, high expression of METTL3, YTHDF2, HNRNPC, YTHDF1, YTHDC2, RBM15, or ALKBH5 was associated with poor survival of HCC ( Figure 1D). Thus, the commonly prognostic m 6 A RNA methylation modulators, namely METTL3, YTHDF2, HNRNPC, YTHDF1 and RBM15, in both TCGA and ICGC were used for the following study.

Construction of CircRNA-miRNA-m 6 A RNA Methylation Modulator Regulatory Network in HCC
By screening of miRNA-seq data from paired tumor and adjacent non-tumor tissues in TCGA HCC cases, a total of 121 DEmiRNAs (29 up and 92 down) were obtained (Figure 2A). From the circRNA microarray data of paired tumor and adjacent nontumor tissues in 3 GEO datasets, a total of 22 DEcircRNAs (eight up and 14 down) were identified ( Figure 2B). 209 miRNA-m 6 A RNA methylation modulators pairs and 1260 circRNA-miRNA pairs were predicted by miRDIP and CSCD, respectively. After taking the intersection of these RNA pairs, 3 m 6 A RNA methylation modulators (HNRNPC, YTHDF1, and YTHDF2), 11 DEmiRNAs, and eight DEcircRNAs were utilized to construct a circRNA-miRNA-m 6 A RNA methylation modulator regulatory network. This network contained 16 circRNA-miRNA pairs and 11 miRNA-mRNA pairs ( Figure 2C). In order to verify the expression stability of DEmiRNAs and DEmRNAs in the regulatory network, we further analyzed the expression of DEmiRNAs in HCC by using dbDEMC 2.0, a database of differentially expressed miRNAs in human cancers (19) and the expression of DEmRNAs in HCC by using HCCDB, a database of hepatocellular carcinoma expression atlas (20), respectively. The results showed that most DEmiRNAs had the same expression trend in multiple GEO datasets, which were consistent with the  results from TCGA and ICGC (Table S3). Especially, hsa-miR-139-5p was down-regulated in HCC tissues from even six datasets. Similarly, in accordance with the TCGA and ICGC results, the expression of YTHDF1, HNRNPC or YTHDF2 was up-regulated in HCC from seven, six, or four datasets, respectively (Table S4).

Clinicopathological Characteristics Correlation and GSEA of m 6 A RNA Methylation Modulators
In accordance with the negative correlation of hsa-miR-139-5p and YTHDF1 expression, the expression of YTHDF1 were higher in high pathological grade and advanced TNM stage in TCGA project ( Figure 4A). While the expression of HNRNPC were significant different between different pathological grades but not between early and advanced TNM stage in TCGA project ( Figure 4B). The relationships between YTHDF1 or HNRNPC expression and clinicopathological characteristics could not be fully investigated due to the lack of data about pathological grade in ICGC project. The expression of YTHDF1 or HNRNPC was significantly higher in the advanced TNM stage of HCC from ICGC project ( Figure S3A). Based on the clinical significance of YTHDF1 in both TCGA and ICGC, we further performed GSEA to explore whether biological pathways differ between high and low YTHDF1 expression groups. Among top five gene sets based  on NES, cell cycle, base excision repair and homologous recombination were enriched in YTHDF1 high expression group in both TCGA and ICGC ( Figure 4C, Figure S3B). While fatty acid metabolism, retinol metabolism, complement and coagulation cascades were enriched in YTHDF1 low expression group in both TCGA and ICGC ( Figure 4C, Figure S3B).

Prognostic Value of the hsa-miR-139-5p and YTHDF1 Signature and Hub circRNA Network Construction
Based on the co-expressed results of miRNA and m 6 A RNA methylation modulators, we further explored the prognostic value of 3 miRNA. Survival analysis revealed that only the hsa-miR-139-5p expression status was correlated with the OS for HCC patients ( Figure 5A). Similar results were observed in GSE31384 ( Figure  5B). According to the results that hsa-miR-139-5p high or YTHDF1 low expression was associated with the better OS of HCC, we further evaluated the prognosis of HCC patients with hsa-miR-139-5p high and YTHDF1 low expression. The results showed that HCC patients with hsa-miR-139-5p high and YTHDF1 low expression had longer OS time than those with contrast expression level ( Figure 5C). Based on the clinical significance of hsa-miR-139-5p and YTHDF1, a hub circRNA-miRNA-mRNA regulatory network was constructed finally. This hub network contained two regulatory axes, namely h s a _ c i r c _ 0 0 0 7 4 5 6 / h s a -m i R -1 3 9 -5 p / Y T H D F 1 a n d hsa_circ_0091570/hsa-miR-139-5p/YTHDF1 ( Figure 5D).

CircMAP2K4 Promotes HCC Proliferation by Modulating hsa-miR-139-5p/YTHDF1
Two potentially dysregulated circRNAs were predicted to involve in regulating hsa-miR-139-5p/YTHDF1 axis in this study. We further analyzed the potential interaction of circRNA and  miRNA through miRanda v3.3a, a microRNA target scanning algorithm. The result showed that hsa_circ_0007456 is predicted to has a score of 140 and energy of -19.98 kCal/Mol to interact with hsa-miR-139-5p, which is favorable for hsa_circ_0007456 serving as hsa-miR-139-5p sponge. But no predicting results were provided for the interaction of hsa_circ_0091570 and hsa-miR-139-5p. Thus, we selected the hsa_circ_0007456 for the following study. Hsa_circ_0007456 (circMAP2K4) derived from mitogen activated protein kinase kinase 4 (MAP2K4) gene and its position is located in chr17:11984672-12016677. The expression of circMAP2K4 or hsa-miR-139-5p was the highest in Huh7; moderate in Hep3B and MHCC97H; and lowest in HCCLM3 ( Figure S4A). In contrast, the mRNA and protein expression of YTHDF1 were the highest in HCCLM3; moderate in Hep3B and MHCC97H; and lowest in Huh7 ( Figure S4A). Agarose gel electrophoresis results showed that the amplified product of divergent primers for circMAP2K4 could be detected in cDNA but not in gDNA. (Figure S4B). Sanger sequencing also confirmed the junction site of circMAP2K4 provided by Circular RNA Interactome ( Figure S4C). These results indicated that circMAP2K4 was a covalently closed-loop RNA. In order to explore the function of YTHDF1, two siRNAs were designed to knockdown the expression of YTHDF1 in MHCC97H and HCCLM3 cells. The results showed that these two siRNAs could effectively decrease the mRNA and protein expression of YTHDF1 in MHCC97H and HCCLM3 cells ( Figure S4D). Transfection of YTHDF1 siRNAs significantly inhibited the proliferation of MHCC97H and HCCLM3 cells ( Figure 6A). Next, we explored the function of circMAP2K4 on regulating hsa-miR-139-5p/YTHDF1 axis. Transfection of hsa-miR-139-5p mimics or circMAP2K4 expressing plasmid could effectively increase the expression of hsa-miR-139-5p ( Figure S4E) or  circMAP2K4 ( Figure S4F), respectively. The luciferase assay showed that the luciferase activity was inhibited when cotransfection of hsa-miR-139-5p mimics and reporter plasmids containing circMAP2K4 wide type sequence with hsa-miR-139-5p binding site. While the luciferase activity had no obvious change when co-transfection of hsa-miR-139-5p mimics and reporter plasmids containing circMAP2K4 mutant sequence ( Figure S4G). Using the miRanda algorithm, we found that Similarly, the luciferase activity of YTHDF1 wide type reporter plasmids was significantly inhibited by transfection of hsa-miR-139-5p mimics ( Figure S4G). Transfection of hsa-miR-139-5p mimics significantly down-regulated the mRNA and protein expression of YTHDF1 in HCC cells, whereas these effects were reversed by circMAP2K4 overexpression ( Figure 6B). Moreover, enforced hsa-miR-139-5p expression significantly inhibited the proliferation of HCC cells. However, additive circMAP2K4 overexpression partly abrogated the inhibitory effect of hsa-miR-139-5p on cell proliferation ( Figure 6C). These results suggested that circMAPK4 acts as has-miR-139-5p sponge to regulate the expression and activity of YTHDF1.

DISCUSSION
The m 6 A RNA modification is a dynamic and reversible process, which is related to various diseases such as obesity, infertility and cancer (21). Numerous studies have confirmed that circRNAs act as miRNA sponges to modulate the pathogenesis of cancer, which facilitates them to serve as diagnostic and prognostic biomarkers and even therapeutic targets for tumors, including HCC (22). In the present study, we demonstrated that most m 6 A RNA methylation modulators are up-regulated and correlated with the prognosis of HCC. Based on the prognostic m 6 A RNA methylation modulators, a circRNA-miRNA-mRNA regulatory network was constructed. Among this network, hsa-miR-139- 5p/YTHDF1 axis was illustrated to be associated with clinicopathological characteristics and prognosis of HCC. Moreover, circMAP2K4 could serve as the hsa-miR-139-5p sponge to up-regulate YTHDF1 expression and promote HCC proliferation. In this study, METTL3, YTHDF2, HNRNPC, YTHDF1, and RBM15 were identified as the commonly prognostic m 6 A RNA methylation modulators for HCC in both TCGA and ICGC projects. Finally, 3 m 6 A RNA methylation modulators, namely HNRNPC, YTHDF1, and YTHDF2, were included for the construction of circRNA regulatory network. Our findings were consistent with previous reports that high expression of YTHDF1, HNRNPC, or METTL3 is related to the poor prognosis of HCC (23)(24)(25). While high level of METTL3 or YTHDF2 can be used as the poor prognostic factor for hepatoblastoma (26). In addition, the combination of YTHDF1 and METTL3 can reflect the malignant degree and evaluate the prognosis of HCC (27). Combining the biomarkers reported in the previous reports and the prognostic m 6 A RNA methylation modulators found in this study, we may try to construct a predictive signature related to the prognosis of HCC in the future study. This may be more beneficial to the evaluation of the prognosis of patients with HCC.
Accumulating evidence revealed that m 6 A writers, erasers, and readers participate in the development and progression of HCC by targeting various tumor-related genes. For example, overexpression of METTL3 promotes cell proliferation, migration, and clonal formation by inhibiting suppressor of cytokine signaling 2 mRNA expression in a m 6 A-YTHDF2dependent manner (28). RAD52 motif 1 (RDM1) binds to the tumor suppressor p53 and enhances its stability, while METTL3 overexpression can significantly reduce the expression of RDM1 mRNA through m 6 A modification (24). In addition, knockdown of METTL3 results in the down-regulation of Snail, a key transcription factor of epithelial-mesenchymal transition (EMT), thereby reducing the invasion and EMT of HCC cell lines (29). As an m 6 A reader, YTHDF is also involved in the occurrence of HCC. YTHDF1 promotes HCC progression by enhancing FZD5 mRNA translation or AKT/GSK-3b/b-catenin signaling activation (30,31). YTHDF2 is involved in the decay of IL11 and Serpine2 mRNA, which are important genes that regulate the normalization and inflammation of vessels (32). The expression of YTHDF2 in HCC is specifically induced by hypoxia, and overexpression of YTHDF2 inhibits cell proliferation, tumor growth and the activation of MEK and ERK. Mechanistically, YTHDF2 can bind to the 3′UTR m 6 A modification site of EGFR, and down-regulate the expression of EGFR mRNA in HCC cells (33). However, the specific mechanisms of HNRNPC in the development of HCC remain unclear. In this study, GSEA results showed that cell cycle, base excision repair and homologous recombination were enriched in YTHDF1 high expression group. Emerging studies had showed that cell cycle dysregulation and DNA damage repair are involved in HCC progression (34,35), which indicates that high expression of YTHDF1 may promote the progression of HCC by regulating cell cycle progression and DNA damage repair. We also confirmed that knockdown the expression of YTHDF1 could inhibit the proliferation of HCC. These findings explain to some extent the reasons for high YTHDF1 expression was associated with high pathological grade, advanced TNM stage and poor survival of HCC in both TCGA and ICGC project. Combining the reported results with our findings indicate that m 6 A RNA methylation modulators play important roles in the occurrence and development of HCC.
Previous studies have demonstrated that some miRNAs participate in the regulation of RNA methylation modulators. For instance, hsa-miR-145 and YTHDF2 mRNA levels were negatively correlated in HCC tissues and overexpression of hsa-miR-145 could down-regulate the expression of YTHDF2, thereby increasing the mA level in HCC cells (36). In hepatoblastoma, METTL3 is identified as a direct target of hsa-miR-186, and the hsa-miR-186/METTL3 axis participates in the progress of hepatoblastoma through the Wnt/b-catenin signaling pathway (26). In this study, YTHDF1 was predicted to be the target gene of hsa-miR-139-5p and hsa-miR-767-5p with a coexpression relationship, and HNRNPC was co-expressed with and predicted to be targeted by hsa-miR-335-5p. There are no reports about these miRNAs acting on these m 6 A regulators right now. However, all three miRNAs have been shown to be involved in the development of HCC. For example, high expression of hsa-miR-139-5p and hsa-miR-335-5p can inhibit the proliferation and invasion of HCC cells and induce tumor shrinkage (37,38). While HCC cells with hsa-miR-767-5p overexpression have significantly higher proliferation, migration and invasion potential (39). In addition, many studies have confirmed that high expression of hsa-miR-139-5p is associated with a better prognosis of HCC (37,40). Different from previous studies, we firstly found that hsa-miR-139-5p could act on m 6 A RNA methylation modulator to regulate the progress of HCC. The present study demonstrated that the expression of hsa-miR-139-5p is negatively correlated with YTHDF1. High expression of hsa-miR-139-5p is associated with a lower grade, an earlier clinical stage, and a better prognosis of HCC, while high expression of YTHDF1 shows the contrast relationship. Moreover, in-vitro experiments also demonstrated that overexpression of hsa-miR-139-5p could inhibit the proliferation of HCC by targeting YTHDF1. These findings suggest that hsa-miR-139-5p/YTHDF1 regulatory axis play an important role in the development of HCC.
According to the important role of hsa-miR-139-5p/YTHDF1 regulatory axis in the progression of HCC, we further evaluated the candidate circRNAs involving in regulating hsa-miR-139-5p and hsa_circ_0007456 and hsa_circ_0091570 were identified finally. As for hsa_circ_0091570, it serves as hsa-miR-1307 sponge and its inhibition promotes HCC cell proliferation, migration and tumor growth in the mouse xenograft model (41). However, to date, the function of hsa_circ_0007456 (circMAP2K4) has not been reported in HCC, which encourages us to test whether circMAP2K4 can also function as miRNA sponge like hsa_circ_0091570 in regulating the proliferation of HCC. We firstly demonstrated that circMAP2K4, hsa-miR-139-5p and YTHDF1 participate in regulating the proliferation of HCC. Importantly, we validated that circMAP2K4 has the bind site for hsa-miR-139-5p and could reverse the repression of hsa-miR-139-5p on YTHDF1, thus eliminating the inhibitory effect of hsa-miR-139-5p on cell proliferation. YTHDF1 high expression was correlated with high pathological grade and advanced stage, which indicates YTHDF1 may be involved in the migration and metastasis of HCC. Indeed, previous studies have confirmed that upregulation of YTHDF1 improve the migratory and invasive capabilities of HCC cells (30,31), which provide clues for the investigation of circMAP2K4/ miR-139-5p/YTHDF1 axis in the migration and metastasis of HCC in our future study.
There are several limitations in this study. First, the datasets included in this study were from different sources, in which the circRNA data were obtained from the microarray data of GEO, while the data of miRNA and m 6 A RNA methylation modulators were obtained from the sequencing data of TCGA or ICGC. Different data sources may affect the reliability of the conclusions to a certain extent. There was no survival information related to the circRNA microarray of HCC in GEO database and due to our lack of HCC tissues, we were unable to evaluate the prognostic value of circRNAs for HCC. Second, due to the limited datasets included in this study, it may lead to that these identified DERNAs were not the most representative although we have verified the expression status of DERNAs through different databases. Moreover, with the advancement of technology and research, more and more m 6 A RNA methylation modulators are discovered. In this study, only 13 m 6 A RNA methylation modulators were included for analysis, and the threshold for circRNA selection was relatively strict, which may lead to the loss of some DEcircRNAs and m 6 A RNA methylation modulators.
In conclusion, we utilize the m 6 A RNA methylation modulators with prognostic value combined with DEcircRNAs and DEmiRNAs to construct a circRNA regulatory network in HCC. Among this network, the expression of hsa-miR-139-5p was negatively correlated with YTHDF1. hsa-miR-139-5p low or YTHDF1 high expression was illustrated to be associated with high grade, advanced stage and poor prognosis of HCC. The hub circRNA regulatory network was constructed based on hsa-miR-139-5p/YTHDF1 axis. In the hub network, circMAP2K4 could serve as the hsa-miR-139-5p sponge to up-regulate YTHDF1 expression and promote HCC proliferation. Our findings indicate that certain circRNA regulatory network is involved in the regulation of m 6 A RNA methylation modulators and provide a novel insight into mechanism study and therapeutic targets for HCC.

DATA AVAILABILITY STATEMENT
Publicly available data sets were analyzed in this study. These data can be found here: TCGA; ICGC; GSE31384; GSE94508; GSE97332; GSE78520.

AUTHOR CONTRIBUTIONS
YuC designed the study and conducted the experiment. FC