Clinical and prognostic pan-cancer analysis of m6A RNA methylation regulators in four types of endocrine system tumors

Results

The expression of m6A regulatory genes is related to endocrine system tumors

In view of the important role of m6A regulatory genes plays in the development of tumors, we systematically studied the relationship between 23 m6A regulatory genes and each endocrine system tumor. The expression level of each m6A regulatory gene in endocrine system tumors was presented in the form of a heat map. The results show that, except for HNRNPC gene, the expression level of other m6A regulatory genes in ACC (Figure 1A) was distinct from that in normal tissues. Except for the IGF2BP3 gene, the expression of m6A regulatory genes in PCPG (Figure 1B) generally differed from normal tissue. The expression of 23 m6A regulatory genes in THCA (Figure 1A) was generally different from that in normal tissue. Except for IGF2BP2 and FTO gene, the expression of other m6A regulatory genes in THYM (Figure 1B) was generally distinct from normal tissue.

Figure 1. The distribution of m6A regulatory genes in endocrine system tumors. (A–C) Heatmaps defined by the expression level of 23 m6A regulatory genes (red for up-regulated, green for down-regulated and the tree on the left shows the clustering results of different genes in different samples) and vioplots visualized the differential m6A regulatory genes in endocrine system tumors (assuming blue is normal tissue and red is tumors); (D) Correlation analysis of expression of 23 m6A regulatory genes in endocrine system tumors.

Then, due to different clinical variables, we divided endocrine system tumors into “Total 1” (ACC and THCA datasets) and “Total 2” (PCPG and THYM datasets). The expression level of m6A regulatory genes of each “Total” was analyzed. In the dataset of “Total 1” (Figure 1C), 23 genes were found to be different from normal tissue. In the dataset of “Total 2” (Figure 1C), except for ZNF217, the expression of m6A regulatory genes is generally different from that of normal tissue. These results indicate that most of the expression of m6A regulatory genes are associated with endocrine tumors.

We then analyzed the correlation of 23 gene expressions in tumors. Our study found that most of expression level of the genes are negatively correlated with METTL3 expression, but there are some exceptions. With the expression level of WTAP is elevated, METTL3 expression is in high level. This can be observed in each dataset (Figure 1D). These datasets illustrate the interaction of gene expression level between m6A regulatory genes.

Consensus clustering method for m6A regulators identified two clusters of endocrine system tumors

In the view of the result that most of the expression of m6A regulatory genes are associated with endocrine tumors, we further investigated the relationship between m6A regulatory genes and the clinicopathological characteristics of the tumors. Clustering classification was performed based on the expression correlation of the m6A RNA methylation regulators, and the clustering stability was increased from k = 2 to 10. The clinicopathological features of the two subgroups clustered were cluster 1 and cluster 2 according to k=2 in ACC datasets (Figure 2A). We further used principal component analysis (PCA) to compare the expression curves between cluster 1 and cluster 2 subgroups, and found that there were significant differences between subgroups (Figure 2A). It is showed that the cluster 1 subgroup (Figure 2D) was significantly correlated with better survival status (P < 0.01), higher who grade (P < 0.01), and higher T status (P < 0.01). The cluster 2 subgroup was associated with worse survival status, lower who grade (P < 0.01), and lower T degree (P < 0.01). In addition, we observed that the overall survival rate (OS) of the cluster 1 subgroup was significantly shorter than that of the cluster 2 subgroup (Figure 2D). We then clustered the other tumors in the same way, observed some similar results and other differences. According to the THCA dataset, we observed no significant difference in OS (Figure 2D) between the cluster 1 subgroup and the cluster 2 subgroup, and it may be necessary to observe the survival in patients with the disease for a longer time. Besides, if only from the aspects of living condition and gender analysis PCPG and THYM datasets (Figure 2B, E), we found no significant difference of clinical features.

Figure 2. Identification of consensus clusters by m6A RNA methylation regulators in endocrine system tumors. (A–C) Consensus clustering matrix for k = 2, and consensus clustering cumulative distribution function (CDF) for k = 2 to 9, and relative change in area under CDF curve for k = 2–9, and principal component analysis (PCA) of total RNA expression profiles in tumors data from the cancer genome atlas (TCGA) dataset; (D, E) Heatmaps, and Kaplan–Meier overall survival (OS) curve and clinicopathologic features of two clusters defined by consistent expression of the m6A regulatory genes (clusters1/2).

Then we conducted a cluster analysis on “Total”. Due to incomplete information of age and tumor metastasis in ACC clinical variables, these two variables will not be analyzed. Both the clinicopathological characteristics of cluster 1 and cluster 2 subgroups clustered by k=2 (Figure 2C), principal component analysis (PCA) was further used to compare the expression curves between cluster 1 and cluster 2 subgroups (Figure 2C). we found that the cluster 1 subgroup was associated with lower survival status (P < 0.01), higher who grade (P < 0.01), greater T (P < 0.01), and more lymph node metastasis (P < 0.01) in “Total 1” (Figure 2D). We observed cluster 1 subgroup of OS is not different from cluster 2 subgroup, only from the aspects of living condition and gender analysis, found the cluster 1 subgroup and gender (P < 0.01) in “Total 2” (Figure 2E). In general, the results showed that the expressions of m6A regulatory genes were closely related to the clinicopathological characteristics of the tumors.

The prognostic role of m6A regulators in endocrine system tumors

In order to evaluate the prognostic effect of the m6A regulatory genes in tumors, we performed univariate Cox regression analysis on the expression level of m6A regulatory genes in the dataset. The results showed that IGF2BP1, HNRNPA2B1, METTL14 and IGF2BP3 were significantly correlated with OS of ACC patients (Figure 3A). METTL14 was the protective gene for HR = 1, and IGF2BP1, HNRNPA2B1 and IGF2BP3 were the risk genes for HR = 1. For genes with prognostic values in PCPG, RBMX, HNRNPC, IGF2BP3 and HNRNPA2B1 were significantly correlated with the OS of PCPG patients (Figure 3B), and were all risk genes with HR = 1. For genes with prognostic value in THCA, YTHDF3 and FTO were significantly correlated with OS of THCA patients (Figure 3A), and were all risk genes with HR = 1. For genes with prognostic values in THYM, CBLL1, RBM15B, KIAA1429, WTAP were significantly correlated with OS of THYM patients (Figure 3B), all of which were HR = 1 protective genes.

Figure 3. Risk signatures with three m6A RNA methylation regulators in endocrine system tumors. (A–C) The coefficients calculated by multivariate Cox regression using LASSO are shown, and Kaplan–Meier overall survival (OS) curve of patients was divided into high-risk and low-risk groups according to the risk score.

We also analyzed the prognostic value of each m6A RNA methylation regulator in “Total 1” (Figure 3C). The results showed that among the 23 tested genes, METTL14, WTAP, RBM15B, METTL16, FTO and IGF2BP1 were significantly correlated with OS (P <0.05). Among the six genes, IGF2BP1 is the risk gene of HR = 1, METTL14, WTAP, RBM15B, METTL16 and FTO are the protective genes of HR = 1. In “Total 2” (Figure 3C), WTAP and METTL16 were significantly correlated with OS (P<0.05), and they were both protective genes of HR = 1. In general, WTAP, METTL16 and IGF2BP3 had a significant influence in the occurrence and development of endocrine system tumors. The reduction expression of WTAP and METTL16 can cause poorer survival in endocrine system tumors. IGF2BP3 is danger genes.

The prognostic risk score is closely related to the clinical characteristics of endocrine system tumors

Based on the above data analysis results that m6A regulated gene expression was closely related to the clinical characteristics and prognosis of tumors, we further studied the relationship between prognostic risk score and endocrine tumors. To better determine the clinical outcome of endocrine system tumors using m6A regulatory genes, we applied the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm to 23 prognostic genes in the TCGA dataset, obtained genes significantly associated with prognosis to establish risk characteristics based on the minimum criteria, and obtained the coefficients from the LASSO algorithm to calculate the risk scores of the validation dataset. To investigate the effect of genetic risk on prognosis, patients with endocrine system tumors in the dataset were separated into the low-risk group and the high-risk group according to the median risk scores. The significant differences in OS between the two groups were observed (P<0.05).

In our analysis of the tumors, we found that there were significant differences in survival status, T status and who grade between the high-risk group of ACC (Figure 4A) and the low-risk group, and significant differences in survival status between the high-risk group and the low-risk group of PCPG (Figure 4A), THCA (Figure 4A) and THYM (Figure 4A). Besides, we observed that the OS of high-risk patients in each disease dataset was much shorter than that of low-score patients.

Figure 4. Relationship between risk score, m6A regulatory genes and clinicopathologic characteristics. (A, B) Heatmaps showed the expression levels of different m6A regulatory genes in low-risk and high-risk patients. The distribution of clinicopathological features in low-risk and high-risk groups was compared. (C) ROC curve shows the signature of risk prediction efficiency in endocrine system tumors.

We also analyzed the expression of six selected m6A RNA methylation regulators in the “Total 1” dataset in high-risk and low-risk patients. We observed significant differences in “Total 1” (Figure 4B) between the high-risk group and the low-risk group, and the high-risk group associated with lower survival status (P < 0.01), higher who grade (P < 0.01), greater T (P < 0.01), and more lymph node metastasis (P < 0.01). We then analyzed the “Total 2” dataset (Figure 4B) and failed to find the difference between the high-risk group and the low-risk group in the clinical information of gender and survival status. But OS of the high-risk patients in both “Total” datasets was much shorter than the OS of the low-score patients.

ROC curves showed that the risk score was a perfect predictor of 5-year survival in ACC patients (AUC = 80.1%) (Figure 4C), PCPG (AUC = 88.8%) (Figure 4C), THCA (AUC =72.8%) (Figure 4C) and THYM (AUC = 81.4%) (Figure 4C), “Total 1” (AUC =68.3%) (Figure 4C), and “Total 1” (AUC =70.9%) (Figure 4C).

We then performed univariate and multivariate Cox regression analyses for the dataset to determine whether the risk characteristics were independent predictors. Univariate regression analysis of ACC data showed that risk scores, who rating and T status were all correlated with OS (P <0.01, Table 1). When these factors were included in multivariate Cox regression, T status (P <0.01) and risk scores (P =0.01) remained significantly correlated with OS (Table 1). Univariate regression analysis of PCPG data showed that the risk score was correlated with OS (P =0.004, Table 1). When these factors were included in the multivariate Cox regression, the risk scores (P =0.002) remained significantly correlated with OS (Table 1). Univariate regression analysis of THCA data showed that risk scores (P =0.014), grade (P =0.003), T status (P =0.028), and age (P <0.01) were all associated with OS (Table 1). When these factors were included in multivariate Cox regression, age (P <0.01) remained significantly correlated with OS (Table 1). A univariate regression analysis of THYM data shows that the risk score (P <0.01) is OS related. In “Total 1”, risk scores (P <0.001), who rating (P <0.001), and T status (P <0.001) were all associated with OS by univariate analysis (Table 1). When these factors were included in the multivariate Cox regression, the risk scores (P <0.001), who grade (P =0.002), and T status (P =0.012) remained significantly correlated with OS (Table 1). In “Total 2”, the risk scores (P =0.012) was correlated with OS by univariate analysis (Table 1). When these factors were included in the multivariate Cox regression, the risk scores (P =0.013) remained significantly correlated with OS (Table 1). In general, risk scores could be used to greatly predict the prognosis of patients with endocrine system tumors. Gender was not associated with the prognosis of endocrine tumors.

Table 1. The univariate and multivariate COX regression analysis of m6A regulatory genes for patients' overall survival (OS) in endocrine system tumors.

Cancer types		Univariate		Multivariate
		HR (95%CI)	P	HR	P
ACC	gender	1.100591489	0.814544864	1.038052354	0.930365584
	stage	2.760634449	<0.01	0.732881514	0.571016476
	T	3.237383176	<0.01	3.747329892	0.005749074
	N	1.829492414	0.269274473	2.918855029	0.137599231
	riskScore	1.193940656	<0.01	1.13395091	0.010312723
PCPG	gender	3.131482468	0.198106512	10.54784495	0.0623891
	riskScore	1.001445171	0.004298715	1.002063944	0.001882145
THCA	age	1.155357049	<0.01	1.166359177	0.000826538
	gender	0.868700107	0.860752257	0.731485506	0.748013828
	stage	2.729834716	0.002868951	0.796474584	0.823787006
	T	2.4428418	0.028302365	2.563634273	0.295472549
	M	2.560644338	0.379583475	3.906592763	0.292240082
	N	2.086807866	0.298770285	0.388357974	0.367125281
	riskScore	1.020547508	0.014232509	1.015493346	0.147944117
THYM	gender	0.629466902	0.49190441	0.66426715	0.550430644
	riskScore	<0.01	0.000837672	>1	0.000987518
Total1	gender	1.596802381	0.15281085	1.342885697	0.38062215
	stage	2.73244358	<0.01	2.124943581	0.002348671
	T	3.351611782	<0.01	1.668245363	0.051414123
	N	0.616791825	0.152388789	0.400896738	0.011567273
	riskScore	164.1049178	<0.01	47.82932365	<0.01
Total2	gender	1.271577618	0.642975417	1.180127876	0.749370968
	riskScore	414.944273	0.012117844	407.8618477	0.01288766
Note: Ambiguous variables (Nx, Mx, N/A, discrepancy and Gx) were excluded. ACC, Adrenal cortical carcinoma;
PCPG, Pheochromocytoma and Paraganglioma; THCA, Thyroid carcinoma; THYM, Thymoma;
Total1, Adrenal cortical carcinoma and Thyroid carcinoma; Total2, Pheochromocytoma and Paraganglioma and Thymoma.

Materials and Methods

Data sources

RNA-seq transcriptome data and the corresponding clinicopathological and prognostic information were obtained for 79 ACC patients,183 PCPG patients, 510 THCA patients, 119 THYM patients from the TCGA database (https://cancergenome.nih.gov/). And TCGA normal group data were obtained for 0 ACC patients, 3 PCPG patients, 58 THCA patients, 2 THYM patients. Since the number of TCGA controls was insufficient, we combined the Genotype-Tissue Expression (GTEx) (http://www.bio-info-trainee.com/3705.html) database with a small number of TCGA normal group data to form the control group. We added GTEx_adrenal gland to the ACC control group (n=128), added GTEx_adrenal gland to the PCPG control group (n=128), added GTEx_Thyroid to the THCA control group (n=279), and added GTEx_Heart (Atrial Appendage, Left Ventricle) to the THYM control group (n=377) (Table 2).

Table 2. RNA-seq dataset of tumor group and control group in this study.

Cancer types	Tumor count	Normal count	Amount
ACC	TCGA_ACC 79	GTEx_adrenal gland 128	207
PCPG	TCGA_PCPG 183	TCGA_PCPG 3 + GTEx_adrenal gland 128 = 131	314
THCA	TCGA_THCA 510	TCGA_THCA 58 + GTEx_Thyroid 279 = 337	847
THYM	TCGA_THYM 119	TCGA_THYM 2 + GTEx_Heart (Atrial Appendage, Left Ventricle) 377 = 379	498
Total 1	TCGA_ACC 79+ TCGA_THCA 510=589	GTEx_adrenal gland 128+ TCGA_THCA 58 + GTEx_Thyroid 279 = 465	1054
Total 2	TCGA_PCPG 183+ TCGA_THYM 119=302	TCGA_PCPG 3 + GTEx_adrenal gland 128+ TCGA_THYM 2 + GTEx_Heart (Atrial Appendage, Left Ventricle) 377 =510	812
Note: ACC, Adrenal cortical carcinoma; PCPG, Pheochromocytoma and Paraganglioma; THCA, Thyroid carcinoma; THYM, Thymoma; Total 1, Adrenal cortical carcinoma and Thyroid carcinoma; Total 2, Pheochromocytoma and Paraganglioma and Thymoma.

In order to study the common characteristics of the expression value of m6A regulatory genes between endocrine system tumors, we divided the four groups of tumors for different clinical variables into two groups. “Total 1” is the dataset of ACC and THCA, mainly analyzing the clinical information of their common survival time, survival status, gender, grade, tumors range and lymph node metastasis. “Total 2” is the dataset of PCPG and THYM, mainly analyzing the only clinical information of survival time, survival status and gender (Figure 5).

Figure 5. Clinical data processing and subsequent studies in the pan-oncogenic analysis of m6A regulatory genes in endocrine system tumors. Total 1 is the dataset of ACC and THCA, and Total 2 is the dataset of PCPG and THYM. The figure shows from TCGA and GTEx after downloading data in the database of datasets, Exclude patients with incomplete information after the rest of the sample size and the process of analysis of m6A regulatory genes.

Selection of m6A regulatory genes

We compiled a list of 23 tumor-related m6A regulatory genes [46–49] based on published literature, and then restricted the list to the genes with available RNA expression data in TCGA data. We obtained mRNA expression data of 23 m6A-related genes [50] from the TCGA database and compared the relationship between the expression of 23 m6A-related genes and the clinicopathological variables of endocrine system tumors.

Bioinformatic analysis

We used the “R” package (R v3.6.2) for analysis (Figure 5). And the Practical Extraction and Report Language (Perl) was used to accurately handle text format that require R package analysis. Firstly, we used the “limma” package to conduct gene difference analysis, used the “pheatmap” and the “vioplot” packages to visualize the expression of 23 genes, and then used the “corrplot” package to analyze the correlation of each gene expression in the tumors. Due to the expression similarity, the “ConsensusClusterPlus” package was used to cluster endocrine system tumors into different groups, and the “ggplot2” package was used to analyze the gene expression of each cluster group. Moreover, we removed patients with incomplete clinical information and used the “survival” package to determine the survival prognosis, and used the “pheatmap” package to visualize the difference expression of 23 genes between cluster 2 and cluster 1.

Next, we used the “forestplot” package to construct univariate COX proportional regression model, and evaluated different m6A expression values and used the “glmnet” package to obtain m6A regulatory genes. Moreover, we use the “survival” package and the “survivalROC” package to analyze the survival of the cluster, and used the “pheatmap” package to visualize the relationship between the expression of relative genes and the different clinical factors. Finally, we used the “forestplot” and the “survival” packages to screen out m6A regulatory genes and clinical factors with a significant influence on prognosis.

Statistical analysis

In order to determine m6A RNA methylation regulating factors in the prognostic role of endocrine tumor patients, we use the Cox univariate analysis data from the TCGA and GTEx database. From this, we selected value genes significantly associated with survival (P < 0.05), which we chose for further analysis of the function and developed a potential risk signature with the LASSO Cox regression algorithm [51]. Risk score calculation is as follows:

$\text{Risk\hspace{0.17em}score\hspace{0.17em}}={\displaystyle \sum }_{i=1}^n{\text{Coef}}_i\times {\text{\hspace{0.17em}x}}_i$

where Coef_i is the coefficient, and x_i is the expression value of each selected gene. This formula was used to calculate the risk scores for each patient.

Abbreviations

ACC: Adrenal cortical carcinoma; PCPG: Pheochromocytoma and paraganglioma; THCA: Thyroid carcinoma; THYM: Thymoma; HCC: Hepatocellular carcinoma; METTL3: Methyltransferase-Like 3; METTL14: Methyltransferase-Like 14; METTL16: Methyltransferase-like 16; WTAP: Wilms Tumor 1- Associated Protein; RBM15: RNA binding motif protein 15; RBM15B: RNA binding motif protein 15B; ZC3H13: domain-containing protein 13; CBLL1: cbl proto-oncogene like 1; ZC3H13: zinc finger CCCH-type containing 13; KIAA1429: vir-Like m6A methyltransferase associated; FTO: obesity-associated protein; ALKBH5: αketoglutarate-dependent dioxygenase alkB homolog 5; YTHDF1/2/3: YTH domain family 1/2/3; YTHDC1/2: YTH domain containing1/2; IGF2BP1/2/3: insulin-like growth factor 2 mRNA-binding proteins 1/2/3; HNRNPA2B1: heterogeneous nuclear ribonucleoprotein A2/B1; HNRNPC: heterogeneous nuclear ribonucleoprotein C; ZNF217: Zinc finger gene217; RBMX: RNA binding motif protein X-linked; TCGA: The Cancer Genome Atlas; GTEx: Genotype-Tissue Expression; LASSO: least absolute shrinkage and selection operator; ROC: receiver operating characteristic; OS: overall survival; m6A: N6-methyladenosine.

Author Contributions

X.Z. and H.L. designed the study. X.Z. prepared the TCGA and GTEx data. K.L. and X.Z. performed bioinformatics and biostatistical analyses. K.L. wrote the manuscript. Hai.L. interpreted the clinical data and edited the manuscript. X.Z. checked the statistical and bioinformatic accuracy as an expert in statistics and bioinformatics. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

Funding

This work was supported partly by National Natural Science Foundation of China (81541153); Guangdong Science and Technology Program (2016A050503046, 2015A050502048, 2016B030309002 and 2019B090905011); Guangdong Province Natural Science Foundation (2018A030310155); the Fund of Southern Marine Science and Engineering Guangdong Laboratory (Zhanjiang) (ZJW-2019-007); The Public Service Platform of South China Sea for R&D Marine Biomedicine Resources (GDMUK201808); and the Science and Technology Program of Zhanjiang (2017A06012). The funders had no role in the design of the study; the collection, analysis, and interpretation of the data; the writing of the manuscript; and the decision to submit the manuscript for publication.

References

• 1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020; 70:7–30. https://doi.org/10.3322/caac.21590 [PubMed]
• 2. Vaidya A, Nehs M, Kilbridge K. Treatment of adrenocortical carcinoma. Surg Pathol Clin. 2019; 12:997–1006. https://doi.org/10.1016/j.path.2019.08.010 [PubMed]
• 3. Huang H, Weng H, Chen J. m⁶A modification in coding and non-coding RNAs: roles and therapeutic implications in cancer. Cancer Cell. 2020; 37:270–88. https://doi.org/10.1016/j.ccell.2020.02.004 [PubMed]
• 4. Li Y, Xiao J, Bai J, Tian Y, Qu Y, Chen X, Wang Q, Li X, Zhang Y, Xu J. Molecular characterization and clinical relevance of m⁶A regulators across 33 cancer types. Mol Cancer. 2019; 18:137. https://doi.org/10.1186/s12943-019-1066-3 [PubMed]
• 5. Zhu W, Wang JZ, Xu Z, Cao M, Hu Q, Pan C, Guo M, Wei JF, Yang H. Detection of N6-methyladenosine modification residues (review). Int J Mol Med. 2019; 43:2267–78. https://doi.org/10.3892/ijmm.2019.4169 [PubMed]
• 6. Liu L, Wang Y, Wu J, Liu J, Qin Z, Fan H. N⁶-methyladenosine: a potential breakthrough for human cancer. Mol Ther Nucleic Acids. 2020; 19:804–13. https://doi.org/10.1016/j.omtn.2019.12.013 [PubMed]
• 7. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol. 2019; 20:608–24. https://doi.org/10.1038/s41580-019-0168-5 [PubMed]
• 8. Pinello N, Sun S, Wong JJ. Aberrant expression of enzymes regulating m⁶A mRNA methylation: implication in cancer. Cancer Biol Med. 2018; 15:323–34. https://doi.org/10.20892/j.issn.2095-3941.2018.0365 [PubMed]
• 9. Chen RX, Chen X, Xia LP, Zhang JX, Pan ZZ, Ma XD, Han K, Chen JW, Judde JG, Deas O, Wang F, Ma NF, Guan X, et al. N⁶-methyladenosine modification of circNSUN2 facilitates cytoplasmic export and stabilizes HMGA2 to promote colorectal liver metastasis. Nat Commun. 2019; 10:4695. https://doi.org/10.1038/s41467-019-12651-2 [PubMed]
• 10. Chen XY, Zhang J, Zhu JS. The role of m⁶A RNA methylation in human cancer. Mol Cancer. 2019; 18:103. https://doi.org/10.1186/s12943-019-1033-z [PubMed]
• 11. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019; 112:108613. https://doi.org/10.1016/j.biopha.2019.108613 [PubMed]
• 12. Zhang C, Chen Y, Sun B, Wang L, Yang Y, Ma D, Lv J, Heng J, Ding Y, Xue Y, Lu X, Xiao W, Yang YG, Liu F. m⁶A modulates haematopoietic stem and progenitor cell specification. Nature. 2017; 549:273–76. https://doi.org/10.1038/nature23883 [PubMed]
• 13. Peng S, Xiao W, Ju D, Sun B, Hou N, Liu Q, Wang Y, Zhao H, Gao C, Zhang S, Cao R, Li P, Huang H, et al. Identification of entacapone as a chemical inhibitor of FTO mediating metabolic regulation through FOXO1. Sci Transl Med. 2019; 11:eaau7116. https://doi.org/10.1126/scitranslmed.aau7116 [PubMed]
• 14. Lin S, Choe J, Du P, Triboulet R, Gregory RI. The m(6)A methyltransferase METTL3 promotes translation in human cancer cells. Mol Cell. 2016; 62:335–45. https://doi.org/10.1016/j.molcel.2016.03.021 [PubMed]
• 15. Schumann U, Shafik A, Preiss T. METTL3 Gains R/W Access to the Epitranscriptome. Mol Cell. 2016; 62:323–324. https://doi.org/10.1016/j.molcel.2016.04.024 [PubMed]
• 16. Du M, Zhang Y, Mao Y, Mou J, Zhao J, Xue Q, Wang D, Huang J, Gao S, Gao Y. MiR-33a suppresses proliferation of NSCLC cells via targeting METTL3 mRNA. Biochem Biophys Res Commun. 2017; 482:582–89. https://doi.org/10.1016/j.bbrc.2016.11.077 [PubMed]
• 17. Li T, Hu PS, Zuo Z, Lin JF, Li X, Wu QN, Chen ZH, Zeng ZL, Wang F, Zheng J, Chen D, Li B, Kang TB, et al. METTL3 facilitates tumor progression via an m⁶A-IGF2BP2-dependent mechanism in colorectal carcinoma. Mol Cancer. 2019; 18:112. https://doi.org/10.1186/s12943-019-1038-7 [PubMed]
• 18. Yue B, Song C, Yang L, Cui R, Cheng X, Zhang Z, Zhao G. METTL3-mediated N6-methyladenosine modification is critical for epithelial-mesenchymal transition and metastasis of gastric cancer. Mol Cancer. 2019; 18:142. https://doi.org/10.1186/s12943-019-1065-4 [PubMed]
• 19. Han J, Wang JZ, Yang X, Yu H, Zhou R, Lu HC, Yuan WB, Lu JC, Zhou ZJ, Lu Q, Wei JF, Yang H. METTL3 promote tumor proliferation of bladder cancer by accelerating pri-miR221/222 maturation in m6A-dependent manner. Mol Cancer. 2019; 18:110. https://doi.org/10.1186/s12943-019-1036-9 [PubMed]
• 20. Cai X, Wang X, Cao C, Gao Y, Zhang S, Yang Z, Liu Y, Zhang X, Zhang W, Ye L. HBXIP-elevated methyltransferase METTL3 promotes the progression of breast cancer via inhibiting tumor suppressor let-7g. Cancer Lett. 2018; 415:11–19. https://doi.org/10.1016/j.canlet.2017.11.018 [PubMed]
• 21. Li X, Tang J, Huang W, Wang F, Li P, Qin C, Qin Z, Zou Q, Wei J, Hua L, Yang H, Wang Z. The M6A methyltransferase METTL3: acting as a tumor suppressor in renal cell carcinoma. Oncotarget. 2017; 8:96103–16. https://doi.org/10.18632/oncotarget.21726 [PubMed]
• 22. Taketo K, Konno M, Asai A, Koseki J, Toratani M, Satoh T, Doki Y, Mori M, Ishii H, Ogawa K. The epitranscriptome m6A writer METTL3 promotes chemo- and radioresistance in pancreatic cancer cells. Int J Oncol. 2018; 52:621–29. https://doi.org/10.3892/ijo.2017.4219 [PubMed]
• 23. Chen WW, Qi JW, Hang Y, Wu JX, Zhou XX, Chen JZ, Wang J, Wang HH. Simvastatin is beneficial to lung cancer progression by inducing METTL3-induced m6A modification on EZH2 mRNA. Eur Rev Med Pharmacol Sci. 2020; 24:4263–4270. https://doi.org/10.26355/eurrev_202004_21006 [PubMed]
• 24. Li F, Zhang C, Zhang G. m6A RNA methylation controls proliferation of human glioma cells by influencing cell apoptosis. Cytogenet Genome Res. 2019; 159:119–25. https://doi.org/10.1159/000499062 [PubMed]
• 25. Kwok CT, Marshall AD, Rasko JE, Wong JJ. Genetic alterations of m⁶A regulators predict poorer survival in acute myeloid leukemia. J Hematol Oncol. 2017; 10:39. https://doi.org/10.1186/s13045-017-0410-6 [PubMed]
• 26. Ji G, Huang C, He S, Gong Y, Song G, Li X, Zhou L. Comprehensive analysis of m6A regulators prognostic value in prostate cancer. Aging (Albany NY). 2020; 12:14863–84. https://doi.org/10.18632/aging.103549 [PubMed]
• 27. Tang J, Wan Q, Lu J. The prognostic values of m6A RNA methylation regulators in uveal melanoma. BMC Cancer. 2020; 20:674. https://doi.org/10.1186/s12885-020-07159-8 [PubMed]
• 28. Meng Z, Yuan Q, Zhao J, Wang B, Li S, Offringa R, Jin X, Wu H. The m⁶A-related mRNA signature predicts the prognosis of pancreatic cancer patients. Mol Ther Oncolytics. 2020; 17:460–70. https://doi.org/10.1016/j.omto.2020.04.011 [PubMed]
• 29. Jaffrey SR, Kharas MG. Emerging links between m⁶A and misregulated mRNA methylation in cancer. Genome Med. 2017; 9:2. https://doi.org/10.1186/s13073-016-0395-8 [PubMed]
• 30. Liu J, Yue Y, Han D, Wang X, Fu Y, Zhang L, Jia G, Yu M, Lu Z, Deng X, Dai Q, Chen W, He C. A METTL3-METTL14 complex mediates mammalian nuclear RNA N⁶-adenosine methylation. Nat Chem Biol. 2014; 10:93–95. https://doi.org/10.1038/nchembio.1432 [PubMed]
• 31. Ries RJ, Zaccara S, Klein P, Olarerin-George A, Namkoong S, Pickering BF, Patil DP, Kwak H, Lee JH, Jaffrey SR. m⁶A enhances the phase separation potential of mRNA. Nature. 2019; 571:424–28. https://doi.org/10.1038/s41586-019-1374-1 [PubMed]
• 32. Sorci M, Ianniello Z, Cruciani S, Larivera S, Ginistrelli LC, Capuano E, Marchioni M, Fazi F, Fatica A. METTL3 regulates WTAP protein homeostasis. Cell Death Dis. 2018; 9:796. https://doi.org/10.1038/s41419-018-0843-z [PubMed]
• 33. Tang B, Yang Y, Kang M, Wang Y, Wang Y, Bi Y, He S, Shimamoto F. m⁶A demethylase ALKBH5 inhibits pancreatic cancer tumorigenesis by decreasing WIF-1 RNA methylation and mediating Wnt signaling. Mol Cancer. 2020; 19:3. https://doi.org/10.1186/s12943-019-1128-6 [PubMed]
• 34. Zhang C, Samanta D, Lu H, Bullen JW, Zhang H, Chen I, He X, Semenza GL. Hypoxia induces the breast cancer stem cell phenotype by HIF-dependent and ALKBH5-mediated m⁶A-demethylation of NANOG mRNA. Proc Natl Acad Sci USA. 2016; 113:E2047–56. https://doi.org/10.1073/pnas.1602883113 [PubMed]
• 35. Zhang J, Guo S, Piao HY, Wang Y, Wu Y, Meng XY, Yang D, Zheng ZC, Zhao Y. ALKBH5 promotes invasion and metastasis of gastric cancer by decreasing methylation of the lncRNA NEAT1. J Physiol Biochem. 2019; 75:379–89. https://doi.org/10.1007/s13105-019-00690-8 [PubMed]
• 36. Yang P, Wang Q, Liu A, Zhu J, Feng J. ALKBH5 holds prognostic values and inhibits the metastasis of colon cancer. Pathol Oncol Res. 2020; 26:1615–23. https://doi.org/10.1007/s12253-019-00737-7 [PubMed]
• 37. Wang S, Chai P, Jia R, Jia R. Novel insights on m⁶A RNA methylation in tumorigenesis: a double-edged sword. Mol Cancer. 2018; 17:101. https://doi.org/10.1186/s12943-018-0847-4 [PubMed]
• 38. Ma JZ, Yang F, Zhou CC, Liu F, Yuan JH, Wang F, Wang TT, Xu QG, Zhou WP, Sun SH. METTL14 suppresses the metastatic potential of hepatocellular carcinoma by modulating N⁶-methyladenosine-dependent primary MicroRNA processing. Hepatology. 2017; 65:529–43. https://doi.org/10.1002/hep.28885 [PubMed]
• 39. Liu J, Eckert MA, Harada BT, Liu SM, Lu Z, Yu K, Tienda SM, Chryplewicz A, Zhu AC, Yang Y, Huang JT, Chen SM, Xu ZG, et al. m⁶A mRNA methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer. Nat Cell Biol. 2018; 20:1074–83. https://doi.org/10.1038/s41556-018-0174-4 [PubMed]
• 40. Huang H, Weng H, Sun W, Qin X, Shi H, Wu H, Zhao BS, Mesquita A, Liu C, Yuan CL, Hu YC, Hüttelmaier S, Skibbe JR, et al. Recognition of RNA N6-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nat Cell Biol. 2018; 20:285–295. https://doi.org/10.1038/s41556-018-0045-z [PubMed]
• 41. Shi H, Wang X, Lu Z, Zhao BS, Ma H, Hsu PJ, Liu C, He C. YTHDF3 facilitates translation and decay of N⁶-methyladenosine-modified RNA. Cell Res. 2017; 27:315–28. https://doi.org/10.1038/cr.2017.15 [PubMed]
• 42. Akbari ME, Gholamalizadeh M, Doaei S, Mirsafa F. FTO gene affects obesity and breast cancer through similar mechanisms: a new insight into the molecular therapeutic targets. Nutr Cancer. 2018; 70:30–36. https://doi.org/10.1080/01635581.2018.1397709 [PubMed]
• 43. Li H, Su Q, Li B, Lan L, Wang C, Li W, Wang G, Chen W, He Y, Zhang C. High expression of WTAP leads to poor prognosis of gastric cancer by influencing tumour-associated T lymphocyte infiltration. J Cell Mol Med. 2020; 24:4452–65. https://doi.org/10.1111/jcmm.15104 [PubMed]
• 44. Yu HL, Ma XD, Tong JF, Li JQ, Guan XJ, Yang JH. WTAP is a prognostic marker of high-grade serous ovarian cancer and regulates the progression of ovarian cancer cells. Onco Targets Ther. 2019; 12:6191–201. https://doi.org/10.2147/OTT.S205730 [PubMed]
• 45. Chen L, Wang X. Relationship between the genetic expression of WTAP and bladder cancer and patient prognosis. Oncol Lett. 2018; 16:6966–70. https://doi.org/10.3892/ol.2018.9554 [PubMed]
• 46. Liu S, Li Q, Chen K, Zhang Q, Li G, Zhuo L, Zhai B, Sui X, Hu X, Xie T. The emerging molecular mechanism of m⁶A modulators in tumorigenesis and cancer progression. Biomed Pharmacother. 2020; 127:110098. https://doi.org/10.1016/j.biopha.2020.110098 [PubMed]
• 47. Zhao Y, Shi Y, Shen H, Xie W. m⁶A-binding proteins: the emerging crucial performers in epigenetics. J Hematol Oncol. 2020; 13:35. https://doi.org/10.1186/s13045-020-00872-8 [PubMed]
• 48. Heck AM, Russo J, Wilusz J, Nishimura EO, Wilusz CJ. YTHDF2 destabilizes m⁶A-modified neural-specific RNAs to restrain differentiation in induced pluripotent stem cells. RNA. 2020; 26:739–55. https://doi.org/10.1261/rna.073502.119 [PubMed]
• 49. Han SH, Choe J. Diverse molecular functions of m⁶A mRNA modification in cancer. Exp Mol Med. 2020; 52:738–49. https://doi.org/10.1038/s12276-020-0432-y [PubMed]
• 50. Zhao W, Qi X, Liu L, Ma S, Liu J, Wu J. Epigenetic regulation of m⁶A modifications in human cancer. Mol Ther Nucleic Acids. 2020; 19:405–12. https://doi.org/10.1016/j.omtn.2019.11.022 [PubMed]
• 51. Sauerbrei W, Royston P, Binder H. Selection of important variables and determination of functional form for continuous predictors in multivariable model building. Stat Med. 2007; 26:5512–28. https://doi.org/10.1002/sim.3148 [PubMed]