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
variables
|
Number (%) all=360
|
univariate analysis
|
multivariate analysis
|
HR (95%CI)
|
P value
|
HR (95%CI)
|
P value
|
Age, years
|
|
|
|
|
|
<65
|
236
|
|
|
|
|
>=65
|
77
|
0.953(0.593-1.531)
|
0.842
|
|
|
Gender
|
|
|
|
|
|
female
|
118
|
|
|
|
|
male
|
195
|
1.103(0.733-1.660)
|
0.638
|
|
|
Race
|
|
|
|
|
|
white
|
270
|
|
|
|
|
others
|
43
|
1.145(0.639-2.052)
|
0.649
|
|
|
ISS stage
|
|
|
|
|
|
I
|
194
|
|
|
|
|
II+III
|
119
|
2.531(1.705-3.759)
|
<0.001
|
2.323(1.553-3.476)
|
<0.001
|
BMPC a
|
|
|
|
|
|
<46%
|
157
|
|
|
|
|
>=46%
|
156
|
2.118(1.404-3.197)
|
<0.001
|
|
0.054
|
Cytogenetic abnormalities
|
|
|
|
|
|
yes
|
117
|
|
|
|
|
no
|
196
|
2.662(1.788-3.961)
|
<0.001
|
2.178(1.450-3.271)
|
<0.001
|
Bone lesions
|
|
|
|
|
|
yes
|
210
|
|
|
|
|
no
|
103
|
2.246(1.374-3.670)
|
0.001
|
2.372(1.447-3.890)
|
0.001
|
ZC3H13
|
|
|
|
|
|
low
|
169
|
|
|
|
|
high
|
144
|
0.602(0.405-0.893)
|
0.012
|
|
0.245
|
HNRNPA2B1
|
|
|
|
|
|
low
|
157
|
|
|
|
|
high
|
156
|
1.562(1.047-2.331)
|
0.029
|
|
0.091
|
HNRNPC
|
|
|
|
|
|
low
|
144
|
|
|
|
|
high
|
169
|
1.594(1.062-2.394)
|
0.025
|
|
0.600
|
YTHDF2
|
|
|
|
|
|
low
|
155
|
|
|
|
|
high
|
158
|
1.849(1.233-2.773)
|
0.003
|
1.684(1.115-2.541)
|
0.013
|
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
variables
|
YTHDF2, mean (SD)
|
P value
|
HNRNPC, mean (SD)
|
P value
|
HNRNPA2B1, mean (SD)
|
P value
|
ZC3H13, mean (SD)
|
P value
|
|
|
Age, years
|
|
|
0.220
|
|
|
0.146
|
|
|
0.242
|
|
|
0.348
|
|
<65
|
11.6718(0.31443)
|
|
10.2541(0.36895)
|
|
10.5382(0.39030)
|
|
9.4182(0.42882)
|
|
|
>=65
|
11.7052(0.38974)
|
|
10.3141(0.37951)
|
|
10.4734(0.35631)
|
|
9.4729(0.39573)
|
|
|
Gender
|
|
|
0.700
|
|
|
0.019
|
|
|
0.309
|
|
|
0.014
|
|
female
|
11.7026(0.38356)
|
|
10.2045(0.38197)
|
|
10.4957(0.38109)
|
|
9.3573(0.42531)
|
|
|
male
|
11.6937(0.36648)
|
|
10.3079(0.36105)
|
|
10.5384(0.38372)
|
|
9.4766(0.41293)
|
|
|
Race
|
|
|
0.564
|
|
|
0.675
|
|
|
0.694
|
|
|
0.352
|
|
white
|
11.6925(0.37966)
|
|
10.2711(0.37270)
|
|
10.5208(0.39060)
|
|
9.4431(0.41324)
|
|
|
others
|
11.7251(0.32592)
|
|
10.2548(0.37061)
|
|
10.5319(0.33259)
|
|
9.3598(0.46536)
|
|
|
ISS stage
|
|
|
0.045
|
|
|
0.266
|
|
|
0.044
|
|
|
0.106
|
|
I
|
11.6661(0.33939)
|
|
10.2704(0.34648)
|
|
10.4892(0.36052)
|
|
9.4680(0.41989)
|
|
|
II
|
11.7066(0.40567)
|
|
10.2152(0.42793)
|
|
10.5373(0.43813)
|
|
9.3743(0.44538)
|
|
|
III
|
11.7875(0.42817)
|
|
10.3167(0.39157)
|
|
10.6143(0.38487)
|
|
9.3701(0.39227)
|
|
|
BMPC
|
|
|
0.361
|
|
|
0.518
|
|
|
0.022
|
|
|
0.344
|
|
<46%
|
11.6864(0.34063)
|
|
10.2912(0.34700)
|
|
10.4804(0.35906)
|
|
9.4553(0.41615)
|
|
|
>=46%
|
11.7077(0.40273)
|
|
10.2464(0.39516)
|
|
10.5644(0.40183)
|
|
9.4077(0.42575)
|
|
|
cytogenetic abnormalities
|
|
|
0.116
|
|
|
0.009
|
|
|
0.210
|
|
|
0.035
|
|
yes
|
11.7491(0.38145)
|
|
10.3312(0.36061)
|
|
10.5672(0.39261)
|
|
9.3598(0.45651)
|
|
|
no
|
11.6659(0.36437)
|
|
10.2316(0.37440)
|
|
10.4954(0.37507)
|
|
9.4744(0.39323)
|
|
|
bone lesions
|
|
|
0.760
|
|
|
0.231
|
|
|
0.163
|
|
|
0.429
|
|
yes
|
11.7017(0.38267)
|
|
10.2908(0.37122)
|
|
10.5039(0.37900)
|
|
9.4185(0.42481)
|
|
|
no
|
11.6874(0.35224)
|
|
10.2241(0.37095)
|
|
10.5597(0.38924)
|
|
9.4583(0.41375)
|
|
|
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
Up-regulated miRNAs
|
|
logFC
|
AveExpr
|
t
|
P.Value
|
adj.P.Val
|
B
|
hsa-miR-148a
|
4.353993
|
9.578942
|
12.33122
|
8.27E-16
|
2.99E-13
|
25.82094
|
hsa-miR-365
|
2.539799
|
5.981127
|
4.097674
|
0.000178
|
0.007091
|
0.27561
|
hsa-miR-193b
|
2.329963
|
6.846663
|
3.487717
|
0.001125
|
0.03125
|
-1.47862
|
hsa-miR-451
|
2.312865
|
6.549557
|
3.042909
|
0.00396
|
0.084087
|
-2.65702
|
Down-regulated miRNAs
|
hsa-miR-363
|
-3.31986
|
2.266472
|
-12.6704
|
3.26E-16
|
2.36E-13
|
26.7257
|
hsa-miR-205
|
-3.08017
|
1.481864
|
-9.31884
|
6.14E-12
|
1.48E-09
|
17.09942
|
hsa-miR-28-5p
|
-3.13096
|
3.011141
|
-9.2066
|
8.75E-12
|
1.58E-09
|
16.75138
|
hsa-miR-203
|
-3.0607
|
1.811558
|
-7.5428
|
1.95E-09
|
2.35E-07
|
11.43158
|
hsa-miR-200c
|
-2.14786
|
2.10884
|
-6.22333
|
1.64E-07
|
1.58E-05
|
7.076365
|
hsa-miR-342-3p
|
-3.82554
|
3.944036
|
-6.20442
|
1.75E-07
|
1.58E-05
|
7.013807
|
hsa-miR-92a
|
-2.67398
|
6.138967
|
-5.43078
|
2.35E-06
|
0.000155
|
4.468974
|
hsa-miR-151-5p
|
-3.00543
|
3.450716
|
-4.26654
|
0.000105
|
0.005056
|
0.784062
|
hsa-miR-155
|
-4.53579
|
4.23315
|
-3.95619
|
0.000276
|
0.009972
|
-0.14346
|
hsa-miR-21
|
-2.41606
|
7.872264
|
-2.83685
|
0.006886
|
0.130841
|
-3.16799
|
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).