Expression Status and Prognostic Roles of m6A-related Genes in Multiple Myeloma: YTHDF2 as the independent prognostic factor

Abstract

Background: N6-methyladenosine is the most abundant RNA modification, which plays a prominent role in multiple biology processes, including tumorigenesis. Multiple myeloma (MM) is the second most frequent hematological cancer. However, the expression status and value of m6A-related genes in multiple myeloma remain elusive.

Methods: m6A-related gene expression data and clinical information of MM patients were obtained from the Gene Expression Omnibus (GEO) database. Based on differentially expressed analysis, protein-protein interaction (PPI) analysis, Pearson correlation analysis, Kaplan-Meier survival analysis and Cox regression analysis, the prognostic m6A-associated genes were identified. The receiver operation characteristic (ROC) curves were used to verify the prognostic and diagnostic efficiency. The molecular mechanisms were investigated by differentially expressed mRNA and microRNA analyses with the help of online database ENCORI and POSTAR.

Results: Among 22 m6A-related genes, HNRNPC, RBM15B, RBM15, YTHDF3, YTHDF2, HNRNPA2B1 and IGF2BP2 were significantly upregulated, while ZC3H13, FTO, IGF2BP3, ALKBH5 and YTHDC1 were significantly downregulated in MM patients. The high expression of HNRNPC, HNRNPA2B1, YTHDF2 and the low expression of ZC3H13 were associated with adverse survival. Furthermore, the expression level of YTHDF2 was the independent prognostic factor of MM. ROC curves suggested the great prediction performance for MM patients. Differentially expressed mRNA and microRNA analyses indicated the probable involvement of miR-205/YTHDF2/EGR1 axis.

Conclusions: Our study first systematically analyzed the expression status of m6A-related genes in multiple myeloma, identified HNRNPC, HNRNPA2B1, YTHDF2 and ZC3H13 that could be the potential prognostic biomarkers, especially YTHDF2, which may be implicated in the miR-205/YTHDF2/EGR1 axis.

1. Background

Multiple myeloma is a neoplastic hematological disease characterized by clonal proliferation of malignant plasm cells in bone marrow (BM). This furious proliferation results in the expansion and accumulation of monoclonal proteins to cause the typical CRAB clinical features: hypercalcemia, renal insufficiency/failure, anemia and bone lesions. MM has been the second most frequent hematological malignancy and accounts for 1.79% of all new cancer cases and 2.11% of all cancer deaths worldwide[1]. Proteasome inhibitor (PI)-based combination paradigm and the great introduction of novel targeted drugs and immune therapies, such as monoclonal antibodies (mAbs) and chimeric antigen receptor-engineered T cells (CAR-T), have become established as a cornerstone of therapy throughout the myeloma treatment algorithm[2]. However, multiple myeloma is still an incurable cancer. The two leading obstacles are: the high recurrence rate and highly heterogeneous cell population whose subclone content evolves over time, especially in the high-risk MM patients. Thus, it is urgent to explore the underlying mechanisms of cancer progression and identify potential biomarkers and therapeutic targets for clinical appliances.

N6-methyladenosine (m6A) was first discovered in eukaryotic mRNA in the 1970s, and it has been recognized as the most abundant RNA modifications among more than 160 types of distinct chemically posttranscriptional modifications[3-5]. However, it was not until 2012 that the first whole transcriptome high-throughput sequencing of m6A modification was finished, initiating the detailed investigation of m6A in various diseases, especially in cancers[6, 7]. Similar to DNA and histone methylation, m6A modification is dynamically reversible with the involvement of m6A “writers”, “erasers” and “readers”[8]. The installation of m6A is catalyzed by writers/methyltransferase complex (MTX), including methyltransferase-like protein 3 (METTL3), METTL14, METTL16, WT1-associated protein (WTAP), RNA-binding motif protein 15 (RBM15), RBM15B, KIAA1429 (or named VIRMA) and zinc finger CCCH-type containing 13 (ZC3H13) [9, 10]. Eraser is a class of demethylases, including fat mass and obesity-associated protein (FTO), α-ketoglutarate-dependent dioxygenase homolog 5 (ALKBH5) and ALKBH3 [11-13]. Read proteins can recognize specific m6A sites and regulate the biofunction of target RNA. YT521-B homology domain family proteins 1/2/3 (YTHDF1/2/3), and YTH domain containing 1/2 (YTHDC1/2) are associated with RNA splicing, translation efficiency and degradation[14, 15]. Insulin-like growth factor 2 mRNA-binding proteins 1/2/3 (IGF2BP1/2/3) can promote mRNA stability and translation[16]. Heterogeneous nuclear ribonucleoprotein A2B1 (HNRNPA2B1) plays an important role in microRNA maturation. HNRNPG and HNRNPG, which can mediate RNA abundance and splicing, are considered as “indirect” readers, given that their tendency to preferentially bind to an RNA structure switch induced by m6A modifications [17, 18]. Emerging studies have demonstrated that aberrant m6A modifications of coding and non-coding RNAs play critical roles in tumorigenesis, metastasis and resistance of various cancers, including glioblastoma (GBM), acute myeloid leukemia (AML), hepatocellular carcinoma (HCC) and breast cancer [19-22]. Several m6A genes are involved in both normal and malignant hematopoiesis, such as METTL3, YTHDF2 and FTO [23]. However, the function of m6A in multiple myeloma is still extremely elusive, and there are only a few published articles focusing on this field up to now[24-26]. In addition, previous articles have only taken into consideration one or few m6A-related regulators, no research have been able to conduct a systematic study about the field.

In the present study, we first used the RNA-seq datasets and clinical data of multiple myeloma from the GEO database to explore the differential expression m6A-related genes. The interactions among these regulators were assessed via PPI network and Pearson correlation analysis. Furthermore, Kaplan-Meier survival analysis and Cox regression model were established to identify prognostic predictors, which were further verified with the use of ROC curves. Finally, YTHDF2 was identified as the independent prognostic factor of MM, and the underlying molecular mechanism was investigated by differential mRNA/microRNA expression analyses. Our study first dissected the systematically abnormal expression status of m6A-related genes in MM and may provide novel and promising biomarkers and therapeutic targets.

2. Methods

2.1 Data preparation and processing

The gene expression profiles containing clinical information of multiple myeloma patients were downloaded from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). The raw count data of mRNA expression profiles were preprocessed with normalization and log2 transformation. The overall design and workflow of this study are given in Fig. 1.

GSE47552 upon on the GPL6244 platform (Affymetrix Human Gene 1.0 ST Array transcription gene version) was selected to investigate m6A-related differentially expressed genes (DEGs) and the expression level of m6A-related genes in different stages of multiple myeloma. GSE47882 contains 5 normal plasma cell (NPC) samples, 20 patients with monoclonal gammopathy of unknown significance (MGUS), 33 with high-risk smoldering multiple myeloma (SMM) and 41 with MM. All samples are newly diagnosed without any treatment;

GSE24080 based on GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) was selected to conduct the survival analysis and the analysis of the relationship between expression level and clinical characteristics. GSE24080 includes totally 559 MM patients, and to ensure the homogeneity of this analysis, 313 IgG type MM were finally included. The detailed information of this dataset is shown as below: a. gender:195 (62.3%) males and 119 (37.7%) females; b. age: the median age is 57.75 (29.70 to 76.50); c. race: 270 (86.3%) white people and 43 other races (13.7%);

Four datasets from GEO database were employed to identify DEGs grouped by the YTHDF2 expression level (1) GSE47552; (2) GSE24080; (3) GSE19784 based on GPL570 includes 320 samples obtained from newly diagnosed MM patient; (4) GSE16558 based on GPL6244 platform contains 60 MM patients.

The GSE17498 dataset from GEO database was selected to identify differentially expressed microRNAs (DEMi) in MM; it upon GPL8227 platform (Agilent-019118 Human miRNA Microarray 2.0 G4470B, miRNA ID version) includes 38 multiple myeloma samples, 2 plasma cell leukemia (PCL) samples and 3 normal donors.

2.2 Identification of differentially expressed m6A genes

We systematically analyzed the expression of 22 m6A genes, including seven m6A writers (METTL3, METTL14, WTAP, KIAA1429, RBM15, RBM15B and ZC3H13), three m6A erasers (FTO, ALKBH5, ALKBH3) and twelve m6A readers (IGF2BP1/2/3, YTHDF1/2/3, YTHDC1/2, HNRNPA2B1, HNRNPC, RBMX) in GSE47552 by using the LIMMA package [27] (version 3.46.0; http://bioinf.wehi.edu.au/limma/) from the R Studio (version 1.4.1103; https://rstudio.com/). P-value < 0.05 was used as a threshold parameter to identify differentially expressed m6A genes. The heatmap was generated via PHEATMAP package (version 1.0.12; https://cran.rstudio.com/web/packages/pheatmap/index.html).

2.3 PPI network and Pearson correlation analysis

The PPI network among the differentially expressed m6A regulatory genes was constructed using the STRING online database (www.string-db.org). Cytoscape software (version 3.8.2) was used for visualization[28]. Pearson correlation analysis was performed to demonstrate the association among different m6A-related genes based on the Pearson’s correlation coefficient (r) value between -1 and +1. The CORRPLOT package (version 0.84; https://github.com/taiyun/corrplot) was used to visualize the correlation matrix.

2.4 Survival analysis

The expression level of each gene was divided into high expression group and low expression group according to the mean value. Kaplan-Meier curve analysis with log-rank test was used to analysis survival rates between different subgroups. Univariate and multivariate Cox regression analysis was performed to identify the independent prognostic factor. Briefly, variables with P-value<0.10 in univariate analysis were eligible for multivariate Cox regression analysis with forward likelihood ration (LR) method. R studio with SURVIVAL package (version 3.2-7; https://github.com/therneau/survival) was used for visualization.

2.5 Investigation of underlying molecular mechanisms of YTHDF2

Firstly, 60 MM patients in ISS-III stage selected from GSE24080 dataset, 41 MM patients from GSE47552, 60 MM patients from GSE16558 and 320 MM patients from GSE19784 were stratified into high and low YTHDF2 expression groups according to the mean value. Differential expression analysis was performed by LIMMA package, and an absolute log2 fold change (FC)>0.6 and P-value<0.05 was used as the cutoff for DEGs. The microRNA expression data of GSE17498 was used to observe DEMi between MM patients and normal samples. Then, the online Venn diagram website (https://bioinfogp.cnb.csic.es/tools/venny/index.html) was used to obtain the intersection of significantly downregulated genes from these four datasets. Next, the datasets downloaded from the POSTAR3 database were used to identify YTHDF2-targeting RNA among these common downregulated genes [29] (http://lulab.life.tsinghua.edu.cn/postar/). Finally, the information of DEMi and the dataset of their downstream genes from ENCORI database [30] (http://starbase.sysu.edu.cn/) were used to identify the upstream microRNA of YTHDF2. An absolute log2 fold change (FC)>2.0, and P-value<0.05 were used as the threshold for significantly DEMi.

2.6 Statistical analysis

Statistical analyses were conducted using SPSS 26.0, R Studio and GraphPad Prism software (version 8.0.2). Continuous variables were described as the mean ± standard deviation (SD); Mann-Whitney test and Kruskal-Wallis test were used to test the difference of subgroups as appropriate; Student’s t test and One-way ANOVA followed by Bonferroni post hoc comparison were employed to compare the difference in expression of m6A genes. ROC curves and the area under the ROC curve (AUC) were examined to determine the predictive value in MM diagnosis and prognosis. P<0.05 was considered statistically significant if not specified. All statistical tests were two-sided.

3. Results

3.1 Twelve of twenty-two m6A-related genes were differentially expressed in multiple myeloma

The transcriptome data from the GSE47552 dataset was employed to identify the differential expression of twenty-two m6A-related genes in MM (n=41) compared with NPC (n=5). The heatmap showed that the expression of twelve m6A methylation regulators (HNRNPC, ZC3H13, RBM15B, FTO, RBM15, YTHDF3, IGF2BP3, ALKBH5, YTHDF2, YTHDC1, HNRNPA2B1 and IGF2BP2) were differentially expressed in MM patients compared to normal volunteers (Fig. 2a). We observed that HNRNPC, RBM15B, YTHDF3, RBM15, YTHDF2, HNRNPA2B1 and IGF2BP2 were significantly upregulated, while ZC3H13, FTO, YTHDC1, ALKBH5 and IGF2BP3 were significantly downregulated in MM patients (Fig. 2b). In addition, the expression of METTL3, KIAA1429, YTHDF1, METTL14, METTL16, YTHDC2, ALKBH3, WTAP, RBMX and IGF2BP1 was similar between MM and NPC samples (Additional file 1: Table S1).

3.2 PPI network and Pearson correlation analysis among m6A RNA methylation regulators

To further investigate the relationship among these 12 m6A RNA methylation regulators, PPI network was established via STRING database and CYTOSCAPE software (Fig. 3a). Of note, the PPI network demonstrated that RBM15, ALKBH5, YTHDC1, YTHDF2 and YTHDF3 were the hub genes. In addition, IGF2BP3 had no protein-protein association with other regulators, which was hidden in the figure. Moreover, the correlation of mRNA expression was analyzed by the Pearson correlation analysis (Fig. 3b). The correlation between YTHDC1 and HNRNPA2B1 was most positively significant (r=0.72); and the relationship between YTHDC1 and IGF2BP2 was most significantly negative (r=-0.38).

3.3 Identification the prognostic roles of m6A-related genes

To evaluate the prognostic roles of the m6A RNA methylation regulators in MM, Kaplan-Meier survival analysis was conducted among these 12 genes. GSE24080 dataset was selected, which containing 313 IgG type MM patient samples with clinical and survival information. Patients were divided into high expression group and low expression group according to the average expression level of each gene. We found that the high expression levels of YTHDF2, HNRPNC and HNRNPA2B1 were associated with poorer over survival (OS), albeit the OS of the low ZC3H13 expression group was significantly prolonged than that of the high expression group (Fig. 4a-d). However, the expression levels of RBM15B, FTO, RBM15, YTHDF3, IGF2BP3, ALKBH5, YTHDC1 and IGF2BP2 were not correlated with OS (Additional file 2: Fig. S1).

Furthermore, univariate Cox regression analysis demonstrated that International Staging System (ISS) stage, cytogenetic abnormalities, bone lesions, proportion of bone marrow plasm cell (BMPC), ZC3H13, HNRNPA2B1, HNRNPC and YTHDF2 were significantly associated with the OS, which were consistent with the results of Kaplan-Meier survival analysis. Based on these results, multivariate Cox regression analysis showed that ISS stage, cytogenetic abnormalities, bone lesions and the expression of YTHDF2 were independent negative prognostic factors (Table 1). Moreover, the performance of this prognostic predictor was confirmed by ROC curves. As shown in Additional file 3: Fig. S2, the AUC were 0.5611, 0.5815 and 0.5904 in the 3-, 4- and 5-year survival, respectively, which indicated some predictive power of YTHDF2.

Subgroup survival analysis of YTHDF2 showed that higher YTHDF2 expression level was correlated with poorer OS especially in ISS-III stage patients and patients with cytogenetic abnormalities (Fig. 5a-e). The subgroup analysis results of HNRNPC, HNRNPA2B1 and ZC3H13 (Additional file 4: Fig. S3; Additional file 5: Fig. S4) showed that low expression of ZC3H13 was correlated with poorer OS in ISS-II stage patients (P=0.00229).

Table 1. Univariate and multivariate analysis of OS using the Cox proportional hazard regression model

^a 46% is the mean value of the percentage of bone marrow plasm cell of 313 MM patients.

HR: hazard ratio

3.4 Correlation of prognostic m6A-related genes and clinical characteristics

Next, we analyzed the correlation between the expression of these four prognostic risk signature genes and the clinicopathological features of MM patients in the GSE24080 cohort (Table 2). We found gender-related differences in the expression of HNRNPC and ZC3H13, while each of the four was age-independent genes. High expression of YTHDF2 and HNRNPA2B1 were associated with advanced ISS stage. In addition, the higher expression of HNRNPA2B1 was correlated with higher BMPC ratio. Multiple myeloma patients with cytogenetic abnormalities were observed to have higher HNRNPC expression and lower ZC3H13 expression.

HNRNPC gene is located at 14q11.2, thus, we used the GSE14080 dataset to investigate the correlation between aberrant expression of HNRNPC and chromosome 14 translocations. However, although the HNRNPC expression level was increased in patients with chromosome 14 translocations, we found no significant difference between normal FISH (fluorescence in situ hybridization) MM patients and patients with t(4;14)/(14;16)/(11;14) (p=0.576, Additional file 6: Fig. S5A). Furthermore, ZC3H13 and RB1 gene (retinoblastoma-1) are co-located at chromosome 13q14, Del (13q) is considered as the adverse predictor.[31]. We observed the expression level of ZC3H13 was significantly decreased in MM patients with RB deletion compared to normal FISH patients (p=0.008, Additional file 6: Fig. S5B).

Table 2. The correlation of prognostic m6A RNA methylation regulators and clinical characteristics

3.5 Identification the diagnostic roles of m6A-related genes

The development of MM has been widely identified as a multistage process starting with the premalignant condition MGUS and followed by the intermediate stage SMM. The underlying mechanism of this continuous disease spectrum is considered as the genetic instability [32]. To figure out the roles of m6A RNA methylation regulators in the development of MM, GSE47552 dataset, including 5 NPC patients, 20 MGUS patients and 33 SMM patients, was employed. The results showed that the expression of HNRNPA2B1 was significantly different between NPC and SMM, albeit no significant difference was observed between NPC and MM (Fig. 6a). As for HNRNPC, statistically significant increase was observed in the development of MM (Fig. 6b). The expression of YTHDF2 was significantly different in MM compared with MGUS and SMM, while no significant difference was found between MM and NPC (Fig. 6c). In addition, ZC3H13 expression level was significantly decreased with the disease development (Fig. 6d).

Furthermore, ROC curve analysis was conducted to verify the diagnostic value of these four genes (Fig. 6e-h). The AUC value of HNRNPA2B1 in MM for NPC, MGUS and SMM were 0.6293, 0.6402 and 0.8004, respectively (Fig. 6e) As fig. 6f shows, the AUC value of HNRNPC in MM for NPC, MGUS and SMM were 0.9902, 0.7561 and 0.5595, respectively. The AUC value of YTHDF2 were 0.6537, 0.7805, 0.7502 for NPC, MGUS and SMM, respectively (Fig. 6g). In addition, The AUC value of ZC3H13 were 1.0000, 0.7646, 0.7302 for NPC, MGUS and SMM, respectively (Fig. 6h). Our results demonstrated that the expression levels of HNRNPA2B1, HNRNPC, YTHDF2 and ZC3H13 have accurate discernibility ability of the continuous disease spectrum of MM.

3.6 Potential molecular mechanisms of YTHDF2 in MM

Firstly, 481 MM patients from four different datasets (60 in ISS-III stage from GSE24080 dataset, 41 from GSE47552, 60 from GSE16558 and 320 from GSE19784) were stratified into high and low YTHDF2 expression groups according to the mean value, subsequently, differential expression analyses were performed (Additional file 7: Table S2-S5).

The primary biofunction of YTHDF2 is identified as the regulation of m6A-containing RNA degradation[33]. Therefore, to explore the target RNA of YTHDF2, we only selected significantly downregulated genes for the next study. The Venn diagram showed the intersection of significantly downregulated genes from these four datasets (Fig. 7a), which indicated that eight RNAs (EGR1, RNASE2, CLC, S100A12, RNASE3, HBA2///HBA1, HBB and CXCL12) were significantly downregulated among these four cohorts. Next, YTHDF2-targeting RNA was obtained through intersecting eight common downregulated RNAs and the data from the POSTAR3 database (Fig. 7b), which showed that EGR1 is the probable target of YTHDF2 (Additional file 8: Table S6). As expected, the expression of EGR1 was significantly decreased in MGUS, SMM and MM compared to NPC (Fig. 7c).

Finally, GSE17498 dataset including 722 human microRNAs was used to identify DEMi between MM patients and normal samples. As Table 3 shows, 4 miRNAs were upregulated and 10 miRNAs were downregulated in MM patients compared to healthy individuals. Furthermore, to investigate the cause of abnormal YTHDF2 expression, we screened these DEMi and their target mRNA in the ENCORI database. The results indicated that YTHDF2 was the target of miR-205 (Fig. 7d), which was significantly downregulated in MM patients. Furthermore, unsurprisingly, Pearson correlation analyses in GSE24080 and GES19784 indicated that the expression level of YTHDF2 was negatively associated with miR-205 in MM (r=-0.175 and -0.351, respectively; Fig. 7e-f).

Taken together, we found the probable mechanism of miR-205/YTHDF2/EGR1 in the development of MM.

Table 3. The DEMi in GSE17498 dataset between MM patients and healthy dornors

AveExpr: average expression; adj.P.Val: adjusted P value.

S100A12 (S100 calcium binding protein A12), RNASE3 (ribonuclease A family member 3), HBA2///HBA1 (hemoglobin subunit alpha 2///1), HBB (hemoglobin subunit beta), CXCL12 (C-X-C motif chemokine ligand 12).

4. Discussion

Thanks to the introduction of triple-drug combinations and the approval of new agents against MM, the survival of MM patients has dramatically extended in the last decade. However, the medical needs of high-risk patients remain unmet, given the difficulty of eradicating genetically high heterogeneity of constantly evolving subclones [34]. The present risk stratification models are mostly based on tumor load and genetic abnormalities, however, ignore the epigenetic and epitranscriptomic alterations. New inhibitor targeting significant epigenetic regulators histone deacetylase (HDAC) has been applied for relapsed/refractory MM in clinical trials, which showed effective and tolerable anti-tumor ability [35]. Thus, faced with this scenario, it will become increasing urgent to reformulate these models with epigenetic/epitranscriptomic prognosis predictors to allow the currently available therapeutic arsenal to be more individualized [36]. m6A RNA methylation, the most remarkable epitranscriptomic modification, has been widely considered to be involved in tumorigenesis and resistance [8], while its biofunction in MM progression is little known yet. In this study, we first explored the comprehensive roles of m6A-associated genes in the regulation of MM.

The expression status of 22 m6A RNA methylation genes were analyzed using the GEO database. We found 12 out of them were significantly expressed in MM patients (HNRNPC, RBM15B, YTHDF3, RBM15, YTHDF2, IGF2BP2 and HNRNPA1B1 were upregulated, ZC3H13, FTO, YTHDC1, ALKBH5 and IGF2BP3 were downregulated). These findings implicated that m6A-related genes, acting as oncogenes or tumor-inhibiting factors, were dysregulated in MM. Furthermore, we elucidated that the high expression level of HNRNPA2B1, HNRNPC, YTHDF2 and low expression level of ZC3H13 were negatively associated with the OS of MM patients through Kaplan-Meier survival analysis. Univariate and multivariate Cox regression analysis indicated that the ISS stage, cytogenetic abnormality, bone lesion and the expression level of YTHDF2 were independent prognostic factors. Additionally, ROC curves verified the prognostic prediction value of YTHDF2. The initiation and progression of multiple myeloma is a multistep process including MGUS and SMM. Our results also showed m6A-related genes were involved in this development with significantly diagnostic capacity confirmed by ROC curves.

Consistent with the published articles, HNRNPC, as an oncogene, was overexpressed in various cancers and identified as the indispensable factor for tumorigenesis in breast cancer [37]. ZC3H13 was downregulated in colorectal cancer and found to suppress proliferation and invasion via inactivating Ras-ERK pathway [38]. The overexpression of YTHDF2 was observed in liver cancer, AML and GBM, which was implicated in the stemness maintenance of cancer stem cells [19, 39, 40]. HNRNPA2B1 was highly expressed in multiple types of cancer, and it promoted the proliferation of breast cancer cells through STAT3 pathway [41]. However, our results are inconsistent with a few studies. Our results showed that ALKBH5 was upregulated in MM patients; the expression level of FTO and METTL3 had no significant difference. Xu et al. revealed that FTO was upregulated in myeloma and extramedullary myeloma patients, and promoted cancer progression by regulating heat shock transcription factor 1 (HSF1) in a YTHDF2-dependent way [26]. In contrast, Bach et al. emphasized that the overexpression of METTL3 in MM, while no significant differences in the expression of FTO and ALKBH5 [24]. Divergences among these studies are probably associated with sample size, clinical features, tumor heterogeneity and detection methods. Therefore, what remains to be established is the systematic and detailed function of each member in the elaborate m6A methylation network in mediating the progression of multiple myeloma and relative mechanisms by which m6A modification alters gene expression posttranscriptionally.

Noticeably, we confirmed the remarkable relationship between some m6A-related genes and clinical features, especially cytogenetic abnormalities. Our results showed HNRNPC was increased while ZC3H13 was reduced in MM patients with cytogenetic abnormalities. The most common translocations of MM are those involving the IGH gene (14q32); t(14;16) and t(4;16) have been identified as the high-risk predictors by the International Myeloma Group (IMWG)[42]. However, although the HNRNPC expression level was increased in patients with chromosome 14 translocations, we found no significant difference between MM with and without chromosome 14 abnormalities. Additionally, the low expression of ZC3H13 was significantly associated with 13p deletion, which may be one of the mechanisms of the negative prognostic impact of del(13p). Del(13p) was observed in about 46% of MM patients, including monosomy 13 and interstitial deletion [43]. Although chromosome 13 deletion has consistently been associated with poor prognosis [44], a retrospective study observed an intriguingly opposite effect of monosomy 13 and partial deletion. The OS and progression-free survival (PFS) of patients with monosomy 13 were significantly shorten, while the presence of partial deletion was associated with significantly better OS [45]. The tumor-inhibiting factor RB1 is simultaneously lost with chromosome 13 abnormality. A recent research found that RB1 could inhibit programmed death ligand 1 (PD-L1) expression and promote tumor immunity, which may provide new insight to overcome the tumor immune evasion [46]. In contrast, “eraser” FTO was demonstrated to rectify PD-L1 expression and promote immune evasion [47, 48]. However, the roles of other m6A regulators in immune escape remain unclear, and the deepen understanding of m6A RNA methylation may explain the dissatisfactory effect of immune check point inhibitors in certain cancers.

YTHDF2, as a member of m6A “readers”, can selectively recognize and bind to m6A-containing RNA to initiate its degradation in a deadenylation-dependent manner by interacting directly with the SH domain of CNOT1, the scaffolding subunit of the CCR4-NOT deadenylase complex [33]. Accumulating evidences have demonstrated that YTHDF2 can act as an oncogene cooperating with dysregulated “writers” or “erasers” via degrading m6A-containg tumor-inhibiting factors [49]. For example, the overexpression of METTL3 and YTHDF2 in bladder cancer can promote cell proliferation and migration via the decay of tumor suppressors Kruppel-like factor 4 (KLF4) and SET domain containing 7 (SETD7) [50]. In addition, the downregulation of ALKBH5 abrogated the anti-tumor function of period circadian protein 1 (PER1) in a YTHDF2-dependent manner, which was involved in the reactivation of ATM-CHK2-P53/CDC25C signaling pathway [51].

Therefore, to investigate the regulatory mechanism of YTHDF2, we first conducted PPI network and Pearson correlation analyses. The results demonstrated that YTHDF2 can interplay with the m6A “writers” RBM15/RBM15B and “erasers” ALKBH5. The reduction of ALKBH5 cooperating with increased RBM15/RBM15B in MM may result in the elevated m6A level of YTHDF2-targeting mRNAs. Then, we divided MM patients into two groups according to YTHDF2 expression, and found that EGR1, a tumor suppressor significantly downregulated in MM, may be the downstream of YTHDF2. EGR1 plays a significant role in the regulation of cell growth, differentiation, apoptosis and has potent anti-tumor function [52] via participating the tumor suppressor network, including p53 [53] , phosphatase and tensin homolog (PTEN) [54] and AP-1 [55]. Low expression of EGR1 is correlated with high risk and poor survival of MM; and overexpression of EGR1 in MM cell lines leaded to profound proliferation inhibition and apoptosis with more sensitive to bortezomib [55]. Moreover, the combination of lenalidomide and dexamethasone synergistically induced cell cycle arrest and apoptosis via promoting the expression of a series of tumor suppressors including EGR1 [56].

MicroRNA is a class of non-coding RNA (ncRNA) which plays an important role in gene silencing ; and the aberrant microRNA expression is one of the key determinants of oncogenesis, metastasis and resistance. Not only do m6A modifications have regulatory effects on miRNAs, such as HNRNPA2B1, but also miRNAs can regulate m6A modifications by a sequence pairing mechanism[22]. The involvement of microRNA dysregulation in YTHDF2 overexpression has been recognized in cancers. miR-145 was implicated in the regulation of m6A levels of HCC by targeting YTHDF2 [57]. Furthermore, decreased expression of miR-495 and miR-493-3p leaded to upregulated YTHDF2 and enhanced proliferation and invasion in prostate cancer [58, 59]. Therefore, we finally performed differential microRNA expression analysis and observed that YTHDF2 was the target of miR-205, which was decreased in MM and negatively correlated with YTHDF2. Taken together, we observed the potential pathway miR-205/YTHDF2/EGR1 leading to the adverse survival of MM.

However, our study does have some limitations. First, the majority of patients included in our study were white, the prognostic value of m6A-related genes remains to be further verified in Asian and African cohorts. Second, given that our results are based on the RNA-seq data, the expression of these genes at protein level should be further clarified. Third, some important clinical information, including R-ISS stage, mSMART risk stratification, therapeutic strategies, was not available in the GEO database. Thus, experiment verification and clinical investigation with innovative design and orchestrated implement are needed to translate these descriptive results into clinical benefits.

5. Conclusions

In conclusion, our study identified m6A-associated genes ZC3H13, HNRNPC, HNRNPA2B1 and YTHDF2, especially YTHDF2, that could be the potential diagnostic and prognostic predictors in multiple myeloma. Furthermore, the potential pathway miR-205/YTHDF2/EGR1 was constructed to explain the oncogenic role of YTHDF2 in MM. Epitranscriptomics, as a novel subfield of epigenetics, may provide new targets for cancer treatment. Future research focusing on this field is pressing.

Supplementary Material Legends

Additional file 1: Table S1. The results of differential expression m6A-related gene analysis

Additional file 2: Figure S1. The prognostic roles of m6A-related genes. Kaplan-Meier survival curve for YTHDF3, YTHDC1, IGF2BP3, IGF2BP2, RBM15B, RBM15, FTO and ALKBH5.

Additional file 3 Figure S2. The verification of the prognostic role of YTHDF2 in multiple myeloma. ROC curves in GSE24080 dataset.

Additional file 4: Figure S3. Subgroup analysis of ZC3H13, HNRNPA2B1 and HNRNPC for overall survival in MM. Kaplan-Meier survival curves comparing the OS between high and low expression among ISS-I/II/III patients.

Additional file 5: Figure S4. Subgroup analysis of ZC3H13, HNRNPA2B1 and HNRNPC for overall survival in MM. Kaplan-Meier survival curves comparing the OS between high and low expression among patients with/without cytogenetic abnormalities.

Additional file 6: Figure S5. The expression level of HNRNPC and ZC3H13 in MM patients. (A) The expression level of HNRNPC in MM patients with normal FISH and with chromosome 14 translocations. (B) The expression level of ZC3H13 in MM patients with normal FISH and with RB deletion. ns means no significant difference.

Additional file 7: Table S2. The results of significantly downregulated genes of GSE24080; Table S3. The results of significantly downregulated genes of GSE47552; Table S4. The results of significantly downregulated genes of GSE19784; Table S5. The results of significantly downregulated genes of GSE16558.

Additional file 8: Table S6. The binding information of YTHDF2 and EGR1 from POSTAR3 database

Abbreviations

ALKBH5: α-ketoglutarate-dependent dioxygenase homolog 5; AML: acute myeloid leukemia; AUC: area under the ROC curve; BM: bone marrow; BMPC: bone marrow plasm cell; CAR-T: chimeric antigen receptor-engineered T cells; CLC: Charcot-Leyden crystal galectin; CXCL12: C-X-C motif chemokine ligand 12; DEGs: differentially expressed genes; DEMi: differentially expressed microRNA; EGR1: early growth response protein 1; FTO: fat mass and obesity-associated protein; GBM: glioblastoma; GEO: Gene Expression Omnibus; HBA2///HBA1: hemoglobin subunit alpha 2///1; HBB: hemoglobin subunit beta; HCC: hepatocellular carcinoma; HDAC: histone deacetylase;

HNRNPA2B1: Heterogeneous nuclear ribonucleoprotein A2B1; HSC: hematopoietic stem cell;

HSF1: heat shock transcription factor 1; IGF2BP1/2/3: Insulin-like growth factor 2 mRNA-binding proteins 1/2/3; IMWG: International Myeloma Group; KLF4: Kruppel-like factor 4; LSC: leukemia stem cell; m6A: N6-methyladenosine; mAbs: monoclonal antibodies; METTL3: methyltransferase-like protein 3; MGUS: monoclonal gammopathy of unknown significance; MM: Multiple myeloma; mSMART: Stratification for Myeloma and Risk-Adapted Therapy; MTX: methyltransferase complex;

ncRNA: non-coding RNA; NPC: normal plasma cell; OS: overall survival; PCL: plasma cell leukemia;

PER1: Period circadian protein 1; PFS: progression free survival; PI: Proteasome inhibitor; PPI: Protein-protein interaction; RBM15: RNA-binding motif protein 15; RNASE2: ribonuclease A family member 2; ROC: receiver operation characteristic; S100A12: S100 calcium binding protein A12; SD: standard deviation; SETD7: SET domain containing 7; SMM: smoldering multiple myeloma; WTAP: WT1-associated protein; YTHDC1: YT521-B homology domain containing 1; YTHDF1/2/3: YT521-B homology domain family proteins 1/2/3; ZC3H13: zinc finger CCCH-type containing 13.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent to publication

Not applicable.

Availability of data and materials

The datasets of this article were generated from the GEO database, STRING database, POSTAR database and ENCORI database.

Acknowledgments

None.

Competing interests

The authors declare that there is no conflict of interests.

Author Contribution

RL, DW and MZ processed data and performed the analyses. RL prepared all the tables and figures. RL, YS, XMW and ALH conceived, designed, and wrote the manuscript. YS, XMW, DW, MZ, JB and ALH revised the manuscript. All authors approved the final manuscript.

Funding:

This work was supported by Natural Science Foundation of Shaanxi Province (Grant/Award Number: No. S2020-JC-QN-1401).

References

1. Siegel RL, Miller KD, Jemal A: Cancer statistics, 2020. CA: a cancer journal for clinicians 2020, 70(1):7-30.
2. Rajkumar SV, Kumar S: Multiple myeloma current treatment algorithms. Blood cancer journal 2020, 10(9):94.
3. Adams JM, Cory S: Modified nucleosides and bizarre 5'-termini in mouse myeloma mRNA. Nature 1975, 255(5503):28-33.
4. Desrosiers R, Friderici K, Rottman F: Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proceedings of the National Academy of Sciences of the United States of America 1974, 71(10):3971-3975.
5. Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, de Crécy-Lagard V, Ross R, Limbach PA, Kotter A et al: MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic acids research 2018, 46(D1):D303-d307.
6. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, Cesarkas K, Jacob-Hirsch J, Amariglio N, Kupiec M et al: Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 2012, 485(7397):201-206.
7. 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(7):1635-1646.
8. Wang T, Kong S, Tao M, Ju S: The potential role of RNA N6-methyladenosine in Cancer progression. Molecular cancer 2020, 19(1):88.
9. Zeng C, Huang W, Li Y, Weng H: Roles of METTL3 in cancer: mechanisms and therapeutic targeting. Journal of hematology & oncology 2020, 13(1):117.
10. Coker H, Wei G, Moindrot B, Mohammed S, Nesterova T, Brockdorff N: The role of the Xist 5' m6A region and RBM15 in X chromosome inactivation. Wellcome open research 2020, 5:31.
11. Chen Z, Qi M, Shen B, Luo G, Wu Y, Li J, Lu Z, Zheng Z, Dai Q, Wang H: Transfer RNA demethylase ALKBH3 promotes cancer progression via induction of tRNA-derived small RNAs. Nucleic acids research 2019, 47(5):2533-2545.
12. Wang J, Wang J, Gu Q, Ma Y, Yang Y, Zhu J, Zhang Q: The biological function of m6A demethylase ALKBH5 and its role in human disease. Cancer cell international 2020, 20:347.
13. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, Yi C, Lindahl T, Pan T, Yang YG et al: N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nature chemical biology 2011, 7(12):885-887.
14. Li J, Xie H, Ying Y, Chen H, Yan H, He L, Xu M, Xu X, Liang Z, Liu B et al: YTHDF2 mediates the mRNA degradation of the tumor suppressors to induce AKT phosphorylation in N6-methyladenosine-dependent way in prostate cancer. Molecular cancer 2020, 19(1):152.
15. Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, Sun HY, Li A, Ping XL, Lai WY et al: Nuclear m(6)A Reader YTHDC1 Regulates mRNA Splicing. Molecular cell 2016, 61(4):507-519.
16. Huang H, Weng H, Sun W, Qin X, Shi H, Wu H, Zhao BS, Mesquita A, Liu C, Yuan CL et al: Recognition of RNA N(6)-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nature cell biology 2018, 20(3):285-295.
17. Alarcón CR, Goodarzi H, Lee H, Liu X, Tavazoie S, Tavazoie SF: HNRNPA2B1 Is a Mediator of m(6)A-Dependent Nuclear RNA Processing Events. Cell 2015, 162(6):1299-1308.
18. Fischl H, Neve J, Wang Z, Patel R, Louey A, Tian B, Furger A: hnRNPC regulates cancer-specific alternative cleavage and polyadenylation profiles. Nucleic acids research 2019, 47(14):7580-7591.
19. Paris J, Morgan M, Campos J, Spencer GJ, Shmakova A, Ivanova I, Mapperley C, Lawson H, Wotherspoon DA, Sepulveda C et al: Targeting the RNA m(6)A Reader YTHDF2 Selectively Compromises Cancer Stem Cells in Acute Myeloid Leukemia. Cell stem cell 2019, 25(1):137-148.e136.
20. Zhou T, Li S, Xiang D, Liu J, Sun W, Cui X, Ning B, Li X, Cheng Z, Jiang W et al: m6A RNA methylation-mediated HNF3γ reduction renders hepatocellular carcinoma dedifferentiation and sorafenib resistance. Signal transduction and targeted therapy 2020, 5(1):296.
21. Huang H, Weng H, Chen J: m(6)A Modification in Coding and Non-coding RNAs: Roles and Therapeutic Implications in Cancer. Cancer cell 2020, 37(3):270-288.
22. Ma S, Chen C, Ji X, Liu J, Zhou Q, Wang G, Yuan W, Kan Q, Sun Z: The interplay between m6A RNA methylation and noncoding RNA in cancer. Journal of hematology & oncology 2019, 12(1):121.
23. Vu LP, Cheng Y, Kharas MG: The Biology of m(6)A RNA Methylation in Normal and Malignant Hematopoiesis. Cancer discovery 2019, 9(1):25-33.
24. Bach C, Leffler M, Kroenke J, Saul D, Mougiakakos D, Mackensen A, Bruns H: Role of N6-methyladenosine (m6A) RNA modification in Multiple Myeloma. Oncology Research and Treatment 2019, 42:123-123.
25. Sobh A, Kaestner C, Riva A, Licht JD: Identification of Genetic Vulnerabilities and Synthetic-Lethal Targets in NSD2-High Multiple Myeloma. Blood 2019, 134.
26. Xu A, Zhang J, Zuo L, Xu J, Yin X, Cheng Q, Zhang Y, Chen L, Zhao F, Fan F et al: Fto Promotes Extramedullary Progression of Multiple Myeloma By Regulation of HSF1 through m6a RNA Methylation. Blood 2019, 134.
27. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic acids research 2015, 43(7):e47.
28. 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 research 2003, 13(11):2498-2504.
29. Zhu Y, Xu G, Yang YT, Xu Z, Chen X, Shi B, Xie D, Lu ZJ, Wang P: POSTAR2: deciphering the post-transcriptional regulatory logics. Nucleic acids research 2019, 47(D1):D203-d211.
30. Li JH, Liu S, Zhou H, Qu LH, Yang JH: starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic acids research 2014, 42(Database issue):D92-97.
31. Mikhael JR, Dingli D, Roy V, Reeder CB, Buadi FK, Hayman SR, Dispenzieri A, Fonseca R, Sher T, Kyle RA et al: Management of newly diagnosed symptomatic multiple myeloma: updated Mayo Stratification of Myeloma and Risk-Adapted Therapy (mSMART) consensus guidelines 2013. Mayo Clinic proceedings 2013, 88(4):360-376.
32. Dutta AK, Fink JL, Grady JP, Morgan GJ, Mullighan CG, To LB, Hewett DR, Zannettino ACW: Subclonal evolution in disease progression from MGUS/SMM to multiple myeloma is characterised by clonal stability. Leukemia 2019, 33(2):457-468.
33. Du H, Zhao Y, He JQ, Zhang Y, Xi HR, Liu MF, Ma JB, Wu LG: YTHDF2 destabilizes m(6)A-containing RNA through direct recruitment of the CCR4-NOT deadenylase complex. Nat Commun 2016, 7.
34. Corre J, Munshi NC, Avet-Loiseau H: Risk factors in multiple myeloma: is it time for a revision? Blood 2021, 137(1):16-19.
35. Gao X, Shen L, Li X, Liu J: Efficacy and toxicity of histone deacetylase inhibitors in relapsed/refractory multiple myeloma: Systematic review and meta-analysis of clinical trials. Experimental and therapeutic medicine 2019, 18(2):1057-1068.
36. Cardona-Benavides IJ, de Ramón C, Gutiérrez NC: Genetic Abnormalities in Multiple Myeloma: Prognostic and Therapeutic Implications. Cells 2021, 10(2).
37. Wu Y, Zhao W, Liu Y, Tan X, Li X, Zou Q, Xiao Z, Xu H, Wang Y, Yang X: Function of HNRNPC in breast cancer cells by controlling the dsRNA-induced interferon response. The EMBO journal 2018, 37(23).
38. Zhu D, Zhou J, Zhao J, Jiang G, Zhang X, Zhang Y, Dong M: ZC3H13 suppresses colorectal cancer proliferation and invasion via inactivating Ras-ERK signaling. Journal of cellular physiology 2019, 234(6):8899-8907.
39. Dixit D, Prager BC, Gimple RC, Poh HX, Wang Y, Wu Q, Qiu Z, Kidwell RL, Kim LJY, Xie Q et al: The RNA m6A Reader YTHDF2 Maintains Oncogene Expression and Is a Targetable Dependency in Glioblastoma Stem Cells. Cancer discovery 2021, 11(2):480-499.
40. 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(23):4507-4518.
41. Gao LB, Zhu XL, Shi JX, Yang L, Xu ZQ, Shi SL: HnRNPA2B1 promotes the proliferation of breast cancer MCF-7 cells via the STAT3 pathway. Journal of cellular biochemistry 2021, 122(3-4):472-484.
42. Palumbo A, Avet-Loiseau H, Oliva S, Lokhorst HM, Goldschmidt H, Rosinol L, Richardson P, Caltagirone S, Lahuerta JJ, Facon T et al: Revised International Staging System for Multiple Myeloma: A Report From International Myeloma Working Group. Journal of clinical oncology : official journal of the American Society of Clinical Oncology 2015, 33(26):2863-2869.
43. Shaughnessy J, Tian E, Sawyer J, Bumm K, Landes R, Badros A, Morris C, Tricot G, Epstein J, Barlogie B: High incidence of chromosome 13 deletion in multiple myeloma detected by multiprobe interphase FISH. Blood 2000, 96(4):1505-1511.
44. Zojer N, Königsberg R, Ackermann J, Fritz E, Dallinger S, Krömer E, Kaufmann H, Riedl L, Gisslinger H, Schreiber S et al: Deletion of 13q14 remains an independent adverse prognostic variable in multiple myeloma despite its frequent detection by interphase fluorescence in situ hybridization. Blood 2000, 95(6):1925-1930.
45. Binder M, Rajkumar SV, Ketterling RP, Greipp PT, Dispenzieri A, Lacy MQ, Gertz MA, Buadi FK, Hayman SR, Hwa YL et al: Prognostic implications of abnormalities of chromosome 13 and the presence of multiple cytogenetic high-risk abnormalities in newly diagnosed multiple myeloma. Blood cancer journal 2017, 7(9):e600.
46. Jin X, Ding D, Yan Y, Li H, Wang B, Ma L, Ye Z, Ma T, Wu Q, Rodrigues DN et al: Phosphorylated RB Promotes Cancer Immunity by Inhibiting NF-κB Activation and PD-L1 Expression. Molecular cell 2019, 73(1):22-35.e26.
47. Tsuruta N, Tsuchihashi K, Ohmura H, Yamaguchi K, Ito M, Ariyama H, Kusaba H, Akashi K, Baba E: RNA N6-methyladenosine demethylase FTO regulates PD-L1 expression in colon cancer cells. Biochemical and biophysical research communications 2020, 530(1):235-239.
48. Su R, Dong L, Li Y, Gao M, Han L, Wunderlich M, Deng X, Li H, Huang Y, Gao L et al: Targeting FTO Suppresses Cancer Stem Cell Maintenance and Immune Evasion. Cancer cell 2020, 38(1):79-96.e11.
49. Wang JY, Lu AQ: The biological function of m6A reader YTHDF2 and its role in human disease. Cancer cell international 2021, 21(1):109.
50. Xie H, Li J, Ying Y, Yan H, Jin K, Ma X, He L, Xu X, Liu B, Wang X et al: METTL3/YTHDF2 m(6) A axis promotes tumorigenesis by degrading SETD7 and KLF4 mRNAs in bladder cancer. Journal of cellular and molecular medicine 2020, 24(7):4092-4104.
51. Guo X, Li K, Jiang W, Hu Y, Xiao W, Huang Y, Feng Y, Pan Q, Wan R: RNA demethylase ALKBH5 prevents pancreatic cancer progression by posttranscriptional activation of PER1 in an m6A-YTHDF2-dependent manner. Molecular cancer 2020, 19(1):91.
52. Inoue K, Fry EA: Tumor suppression by the EGR1, DMP1, ARF, p53, and PTEN Network. Cancer investigation 2018, 36(9-10):520-536.
53. Krones-Herzig A, Mittal S, Yule K, Liang H, English C, Urcis R, Soni T, Adamson ED, Mercola D: Early growth response 1 acts as a tumor suppressor in vivo and in vitro via regulation of p53. Cancer research 2005, 65(12):5133-5143.
54. Virolle T, Adamson ED, Baron V, Birle D, Mercola D, Mustelin T, de Belle I: The Egr-1 transcription factor directly activates PTEN during irradiation-induced signalling. Nature cell biology 2001, 3(12):1124-1128.
55. Chen L, Wang S, Zhou Y, Wu X, Entin I, Epstein J, Yaccoby S, Xiong W, Barlogie B, Shaughnessy JD, Jr. et al: Identification of early growth response protein 1 (EGR-1) as a novel target for JUN-induced apoptosis in multiple myeloma. Blood 2010, 115(1):61-70.
56. Gandhi AK, Kang J, Capone L, Parton A, Wu L, Zhang LH, Mendy D, Lopez-Girona A, Tran T, Sapinoso L et al: Dexamethasone synergizes with lenalidomide to inhibit multiple myeloma tumor growth, but reduces lenalidomide-induced immunomodulation of T and NK cell function. Current cancer drug targets 2010, 10(2):155-167.
57. Yang Z, Li J, Feng G, Gao S, Wang Y, Zhang S, Liu Y, Ye L, Li Y, Zhang X: MicroRNA-145 Modulates N(6)-Methyladenosine Levels by Targeting the 3'-Untranslated mRNA Region of the N(6)-Methyladenosine Binding YTH Domain Family 2 Protein. The Journal of biological chemistry 2017, 292(9):3614-3623.
58. Du C, Lv C, Feng Y, Yu S: Activation of the KDM5A/miRNA-495/YTHDF2/m6A-MOB3B axis facilitates prostate cancer progression. Journal of experimental & clinical cancer research : CR 2020, 39(1):223.
59. Li J, Meng S, Xu M, Wang S, He L, Xu X, Wang X, Xie L: Downregulation of N(6)-methyladenosine binding YTHDF2 protein mediated by miR-493-3p suppresses prostate cancer by elevating N(6)-methyladenosine levels. Oncotarget 2018, 9(3):3752-3764.