Prognostic values of YTHDF1 regulated negatively by mir‐3436 in Glioma

Abstract M6A methylation is likely to be closely associated with the occurrence and development of tumours. In this study, we demonstrated that the transcription levels of the m6A RNA methylation regulators are closely related to the prognosis of glioma. Univariate Cox analysis was performed on the expression levels of methylation regulators and selected three hub genes in glioma. Next, we systematically compared the expression of these m6A RNA methylation regulators in gliomas with different clinicopathological features. The overall survival (OS) curve of the hub genes was initially established based on TCGA database information. YTHDF1 was selected from the hub genes following survival and prognosis analysis. A nomogram was developed to predict the survival probability. We further performed cell function and in vivo xenograft tumour experiments to further verify its role in tumour progression. Next, based on the miRanda and miRDB databases, we predicted one microRNA, hsa‐mir‐346, that might regulate and bind to 3'UTR of YTHDF1, which was confirmed by our fluorescent enzyme reporter gene experiment. In summary, m6A RNA methylation regulators play a potential role in the progression of gliomas. YTHDF1 may have an essential function in glioma diagnosis, treatment and prognosis.


| INTRODUC TI ON
Glioma is a common primary central nervous system tumour.
Approximately 200,000 new patients are diagnosed every year around the world. The most frequent type of malignant glioma (WHO III-IV) is glioblastoma multiforme (GBM), which accounts for 45.2% of all central nervous system malignancies, with a median overall survival of only 12-18 months. 1 The treatment of glioma has evolved from a single surgical treatment to a combined treatment such as surgery plus radiotherapy and chemotherapy, but the prognosis has not improved significantly. At present, the clinical diagnosis of glioma includes imaging and histopathology that also have limitations. Therefore, it is important to reveal the mechanisms of glioma treatment and prognosis to find new diagnostic targets.
N6-adenine (m6A) is the most abundant chemical modification in eukaryotes, and m6A methylation is an epigenetic modification of RNA molecules. In the 1970s, geneticists discovered this methylation modification on eukaryotic messenger RNA (message RNA, mRNA). 2 The function and role of m6A methylation in different biological processes have again attracted attention. A large number of studies have confirmed that with mainly three types of proteins are involved in m6A methylation. The first is m6A methyltransferase, whose coding gene is called Writer. [3][4][5] The second type of protein is m6A demethylase, whose encoding gene is Erasers, which can remove m6A methylation groups in RNA. 6 Their coding genes are called readers, including two subtypes of YTHDFs and YTHDCs. 7,8 YTHDF1, as a member of the m6A-modified RNA-binding protein family, inhibits its expression in normal lung epithelial cells to resist hypoxia-induced apoptosis and is highly expressed in non-small cell lung cancer tumour tissues and cell lines. 9 Nishizawa et al. confirmed that YTHDF1 can affect tumour progression by modifying m6A methylation levels of inhibitory genes. 10 Another analysis of clinical data showed that patients with m6A hypomethylation had significantly lower disease-free survival (DFS) and overall survival (OS) and a higher recurrence rate (P < 0.01). 11 However, there is a lack of literature on the function and prognostic value of regulators in the malignant progression of gliomas.
We systematically compared the expression of these m6A methylation genes in gliomas with different clinicopathological features of their TCGA datasets. YTHDF1 was selected from the hub genes based on survival and prognosis analysis. Multivariate cox regression analyses were performed, and a nomogram was built with potential risk factors based on a multivariate Cox analysis to predict survival probability. We further performed cell function and in vivo xenograft tumour experiments to further verify its role in tumour progression. Next, based on the miRanda and miRDB databases, we predicted one microRNA, hsa-mir-346, that might regulate and bind to 3'UTR of YTHDF1, which was confirmed by a fluorescent enzyme reporter gene experiment.

| Microarray data
The gene expression profiles and clinical information were downloaded from The Cancer Genome Atlas (TCGA). We selected the expression data of thirteen m6A RNA methylation regulatorrelated gene from the TCGA databases and then compared the associations between different clinical parameters and m6A RNA methylation regulator expression in gliomas. The differential expression of any gene of interest between tumour and adjacent normal tissues was studied in gliomas using R software (version 3.5.1). Analyses of the differential expression and the pairwise difference map of these m6A RNA methylation regulators were performed using R software.

| Kaplan-Meier analysis
We compared the levels of genes expression at the clinical stage using R software. The overall survival analyses of the hub genes were performed trough TCGA databases. Kaplan-Meier plots and mutation of hub genes were analysed by the open-source web tool cBioPortal (http://www.cbiop ortal.org/index.do).

| Prognostic nomogram
Multivariate analyses were conducted to determine the prognostic value of the hub genes and clinical characteristics. A nomogram and risk score were built with potential risk factors (P < 0.05) based on a multivariate Cox analysis using the R software package. The predictive performance of the nomogram was then measured by receiver operating characteristic (ROC) curves and analysed by Kaplan survival probability.

| Cell culture and transfection
Glioma cell lines were cultured in DMEM medium containing 10% foetal bovine serum under the conditions of 37°C, 5% CO 2 and 100% saturated humidity. Passive culture was performed when the cells confluent to 70% to 80%. The cells were cultured in 6-well plates and synchronized with serum-free medium for 24 hours.
Further, the cells were divided into a Si-NC group and a Si-YTHDF1 group. Next, we transfected the blank plasmid and YTHDF1 interfering plasmid following the instructions of the Lipofectamine 2000 transfection kit (Guangzhou Xiangbo Biotechnology Co., Ltd., Guangzhou, China).

| Quantitative real-time polymerase chain reaction and Western blotting
Qqt-PCR was used to detect the changes in the YTHDF1 expression.
The aforementioned reverse transcription products were tested by Takara's SYBRPremixExTaqTM on an ABI7900 instrument using qPCR, and GAPDH expression was utilized as an internal reference. The following primers were used: YTHDF1-5′ACCTGTCCAGCTATTACCCG, TGGTGAGGTATGGAATCGGAG; GAPDH-5′ACCCACTCCT CCACC TTTGA, 3′CTGTTGCTGTAGCCAAATTCGT. Cells from each group were collected, and total cell protein was extracted using cell lysate.
Then, the protein concentration was determined by a BCA kit.
According to the results of protein concentration detection, 25 μL of protein samples was subjected to Western blot detection.

| Colony formation assay
We diluted the cell suspension of a single-layer culture cell in logarithmic growth phase by multiple gradients and inoculated the culture dish with an appropriate cell density (according to the proliferation ability). Next, we discarded the supernatant and carefully dipped it twice with PBS, added 5 mL of pure methanol or 1:3 acetic acid/methanol and fixed it for 15 minutes. Then, we removed the fixative solution, added an appropriate amount of Giemsa stain and applied the staining solution for 10-30 minutes. Finally, we inverted the plate, superimposed a grid of transparencies and counted the clones directly with the naked eye.

| In vivo xenograft tumour growth
Twelve nude mice were divided into an experimental and a control group (n = 6). Nude mice were kept in an environment of about 25°C with a suitable humidity and were provided with sufficient food and water during the administration. Experimental grouping: About 5 × 106 cells transfected with pc-YTHDF1 in SHG-44 cell (si-YTHDF1 in U87 cell) were suspended in 0.1 mL of PBS and injected subcutaneously into nude mice. The control group was given the same amount of normal saline. The volume and weight of the xenograft tumours in the nude mice were then measured every 10 days. Thirty days later, the nude mice were killed, and the transplanted tumours were removed.
The morphological changes in the tumour tissues were observed by H & E staining, and the expression of Ki67 was also detected.

| miRNA database analysis
The potential miRNAs targeting YTHDF1 were downloaded from miRanda and miRDB databases. We performed survival analysis of hsa-miR-346 in brain glioma, downloaded from TCGA. We then conducted survival analysis of hsa-miR-346 in glioma using a TCGA project. Coexpression analysis for hsa-miR-346 and YTHDF1 was performed via ENCORI, which mainly focuses on miRNA-target interactions and is an open-source platform for studying RNA-RNA interactome data (http://starb ase.sysu.edu.cn/index.php).

| Luciferase reporter assay
The YTHDF1 gene promoter was cloned by RT-q PCR and DNA fragments from the 3′-UTR of YTHDF1 inserted into the luciferase reporter vector pGL3. Then, the pGL3-YTHDF1-3′UTR reporter plasmids (100 ng) plus 5 ng of pRL-TK renilla plasmid and increasing levels of NC, Hsa-mir-346 or mutant Hsa-mir-346 mimics were cotransfected into the U87 cell. Luciferase activity analysis was next performed to calculate the luciferase activity ratio of the reporter plasmid and the internal reference.

| Statistical analysis
Statistical analysis was conducted using SPSS (IBM, Chicago, IL, USA) and R software (version 3.5.1). Factors were identified as significant at P < 0.1 in the univariate analysis. Comparisons between groups were performed using independent sample t tests, and multiple group comparisons were done using single-factor variance.
Pairwise comparisons were performed using LSD Lt test. P < 0.05 was considered statistically significant.

| Clinical significance of YTHDF1
We examined the microarray data of all 605 glioma cases from the TCGA database. Heat maps and violin plot showed different gene expression between the normal and the tumour group ( Figure 1A were decreased in glioma ( Figure 1B). Next, we performed univariate cox regression analysis on the expression levels in the TCGA dataset ( Figure 1C). As can be seen in Figure 1D, YTHDF1, RBM15 and METTL14 are associated with clinicopathological features (P < 0.05; Figure 1D). Of these genes, YTHDF1 is the risky gene with HR > 1, whereas RBM15 and METTL14 are protective genes.
To evaluate the prognosis of the three tested genes in glioma, the overall survival (OS) curve was performed first on the Kaplan-Meier Plotter Database. The high expression of YTHDF1 in glioma was associated with worse OS (P < 0.05; Figure 2A). OS of RBM15 and METTL14 were not statistically significant ( Figure 2B,C).
Interestingly, the OS and DFS of the YTHDF1 mutations in glioma were performed via the cBioPortal. The results we obtained showed that glioma cases with related genes mutations had better OS and DFS (P < 0.05; Figure 2D,E). Next, we investigated the association between the YTHDF1 gene expression and the clinical stage using the TCGA database. The results indicated that high levels of YTHDF1 are correlated with advanced stages. The above results indicated that YTHDF1 may contribute to glioma progression ( Figure 2F,G).

| Prognosis of YTHDF1 in glioma patients
In the present study, a total number of 605 glioma patients were included from TCGA. The clinicopathological characteristics of the patients are listed in Table 1. The results of the multivariate analysis revealed that the dimension, subtype, age, grade, stage and YTHDF1 were statistically significant factors for glioma progression ( Figure 2F). The prognostic nomogram that integrated all significant independent factors is presented in Figure 3A. The ROC curve area evaluating the prognostic nomogram for OS prediction was 0.750 (95% CI, 0.729 to 0.759; Figure 3B). The calibration plot for the probability of survival at three or five years revealed an optimal agreement with the prediction by the nomogram, which appeared to be good in the model ( Figure 3C,D).

| YTHDF1 promotes glioma tumour growth in vivo
To further study the biological significance of YTHDF1 in gliomas, we injected U87-siYTHDF1 and SHG-44-pcYTHDF1 cells and their  Figure 6A,B). In addition, the low YTHDF1 levels in glioma cells were associated with low tumour growth rates and smaller tumour weight, whereas the high YTHDF1 levels in glioma cells were associated with higher tumour growth rates and tumour weight ( Figure 6C). Immunohistochemical experiments established that the Ki67 positive rate in the high-YTHDF1 group was significantly higher than that in the control group ( Figure 6D). Taken together, these data suggest that YTHDF1 promotes the growth of glioma tumours in vivo.

| Hsa-mir-346 targets YTHDF1
Finally, we selected two sets of miRNAs related to YTHDF1 from mi-Randa and miRDB databases and obtained miR-486-5p as a possible candidate by drawing Veen map ( Figure 7A). Next, we performed survival analysis of hsa-miR-346 in brain glioma via ENCORI. We found that glioma patients with high expression of Hsa-mir-346 had better overall survival than those with low expression (P < 0.05; Figure 7B).

Moreover, the miRNA-target coexpression between hsa-mir-346
and YTHDF1 via ENCORI showed that hsa-mir-346 expression was negatively correlated with the YTHDF1 expression ( Figure 7C). Importantly, we established that hsa-mir-346 may be combined with 3′UTR of YTHDF1 ( Figure 7D). Therefore, we speculated  Figure 7G). Moreover, we analysed the expression levels of YTHDF1 in u87 cells transfected with hsa-mir-346 and the negative control (NC) by qRT-PCR and Western blotting and found that the expression of YTHDF1 with hsa-mir-346 was significantly lower than that of the negative control group ( Figure 7H). Therefore, these results indicated that YTHDF1 is regulated by hsa-mir-346 and may have an impact on the prognosis of glioma patients.

| D ISCUSS I ON
Previous studies have reported that the pathological grade, age at onset, surgical methods and post-operative adjuvant treatment of glioma patients are directly related to their prognosis. 13,14 F I G U R E 3 A, The prognostic nomogram that integrated all significant independent factors from the multivariate analysis for OS in the training cohort. B, The ROC curve area evaluating the prognostic nomogram for OS prediction was 0.750 (95% CI, 0.729 to 0.759). C, The calibration plot for the probability of survival at 3 or 5 years after surgery showed an optimal agreement between the prediction by the nomogram, respectively. *P < 0.05 versus control However, many factors affect the prognosis of glioma patients, most of which are uncertain. 15 Furthermore, much of the literature on the relationship between M6A methylation and tumours has been mostly focused on liver and cervical cancer. For example, Dominissini confirmed that the knockout of YTHDF1 in cell lines induced apoptosis of liver cancer cells in in vitro experiments and speculated that its mechanism may be through the activation of the p53 signalling pathway. 16 JmzhaoMa found that the expression of YTHDF1in liver cancer was reduced, which was negatively correlated with the prognosis of liver cancer patients. 17 Patients with m6A hypomethylated cervical cancer had significantly lower disease-free survival (DFS) and overall survival (OS), whereas hypomethylation led to a higher recurrence rate (P < 0.01). 11 In this study, we demonstrated that the expression of regulators (RNA m6A RNA methylation) is also closely associated with the progress and prognosis of gliomas. Among them, YTHDF1 is associated with OS and clinicopathological features (P < 0.05). Importantly, YTHDF1 is a risky gene with HR > 1 (P < 0.05).
Nishizawa et al found that the YTHDF1 gene was highly expressed in colorectal cancer patients and was related to the tumour diameter, clinical stage, but its specific regulatory mechanism was not studied. 19 Additionally, Hao et al conducted research on liver cancer data in TCGA and found that YTHDF1 was associated with tumour progression and may be an indicator factor for poor prognosis. 20 The results of these studies are consistent with our findings that the high expression of YTHDF1 in glioma is associated with worse OS and advanced stages.
Furthermore, glioma patients with YTHDF1 mutations have a better overall survival and disease-free survival than those of the wild-type.  It is noteworthy that miR-346 was found to be up-regulated in human cervical cancer tissues and promoted malignant progression. 28 Here, using ENCORI, we found that glioma patients with high expression of hsa-mir-346 had better overall survival than those with low expression (P < 0.05). Importantly, the miRNA-target coexpression between hsa-mir-346 and YTHDF1 established by ENCORI showed that the hsa-mir-346 expression was negatively correlated with the YTHDF1 expression. Next, we confirmed that hsa-mir-346 mimic significantly reduced the luciferase activity of YTHDF1-3'UTR by luciferin experiments. There is evidence that the m6A reader YTHDF1 can enhance the translation of EIF3C and its expression, promoting the malignant progression of ovarian cancer. 29 Therefore, we speculated that hsa-miR-346 is a potential key upstream negative regulator of YTHDF1 that may be related to cancer treatment and prognosis. However, this notion needs to be further explored.
In summary, our study showed that YTHDF1 is associated with glioma progression, and high YTHDF1 expression also predicts a poor prognosis in glioma patients. We also developed a nomogram to predict the overall survival (OS) of glioma patients. Finally, we identified hsa-mir-346 as an upstream regulator of YTHDF1 expression, which may be involved in the development of a new treatment reducing the development of glioma.

ACK N OWLED G EM ENT
The authors gratefully acknowledge contributions from the TCGA Network. This study was generously supported by the Science and Technology Research Project of Henan Province (Grant no. 182102311180).

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interest.

AUTH O R CO NTR I B UTI O N
SL and CX conceived and designed the study, which were proofed by SL. CX analysed the data and wrote the manuscript. BY and TH carried out the concepts, design and manuscript preparation. BY, TH and BD provided assistance for data acquisition. All authors have read and approved the content of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.