Research Paper Advance Articles pp 23931—23944
Clinical and prognostic pan-cancer analysis of m6A RNA methylation regulators in four types of endocrine system tumors
- 1 Guangdong Key Laboratory for Research and Development of Natural Drugs, The Marine Biomedical Research Institute, Guangdong Medical University, Zhanjiang 524023, Guangdong, China
- 2 The Marine Biomedical Research Institute of Guangdong Zhanjiang, Zhanjiang 524023, Guangdong, China
- 3 Southern Marine Science and Engineering Guangdong Laboratory Zhanjiang, Zhanjiang 524023, Guangdong, China
- 4 Cancer Center, Affiliated Hospital, Guangdong Medical University, Zhanjiang 24023, Guangdong, China
Received: July 24, 2020 Accepted: August 28, 2020 Published: November 20, 2020https://doi.org/10.18632/aging.104064
How to Cite
Copyright: © 2020 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
N6-methyladenosine (m6A), internal modification of mRNA, has recently been reported to be an important regulatory mechanism affecting tumor proliferation. However, its role in endocrine system tumors is poorly understood. We obtained datasets for four types tumors from the TCGA database, analyzed the GTEx database as a supplement to the control group, and used “Perl” and “R” software to analyze the datasets. Then we differentiated the expression level, used it to cluster consensus. Besides, we established lasso regression model to screen variables, used univariate and multivariate cox analyses to explore the independent risk factors associated with cancer prognosis. The results indicated that except for WTAP, the expression level of METTL3 was negatively correlated with other genes. The expression level of WTAP and METTL16 was positively correlated with overall survival (OS). Moreover, we found that different clinical subtypes of adrenal cortical carcinoma had significant differences in survival status, histologic grading, pathological T grade, and OS. Furthermore, different clinical subtypes of thyroid carcinoma had significant differences in histologic grading and pathological T grade. The differential expression of m6A regulatory genes is closely associated with the presence of endocrine-system-related tumors, and risk scores can be used to assess prognosis.
According to the recent data, the incidence of endocrine system tumors is increasing year by year, although the proportion is not large which cannot be ignored . Despite advances in the diagnosis and treatment of endocrine tumors (adrenal cortical carcinoma (ACC), pheochromocytoma and paraganglioma (PCPG), thyroid carcinoma (THCA) and thymoma (THYM)) in recent decades, some patients still have a poor prognosis, especially ACC . Recently, a number of studies have found that m6A regulatory genes play important roles in the occurrence and development of tumors , which also include the endocrine system tumors .
N6-methyladenosine (m6A), also known as RNA methylation modification, occurs on the sixth nitrogen atom (N) of adenine (A) and is enriched in mRNA . It is mainly distributed in the coding sequence and 3'UTR region of mRNA. The modification process is dynamic and reversible, involving three parts: writers, erasers and readers [6, 7]. Writers (METTL3, METTL14, METTL16, WTAP, RBM15, RBM15, KIAA1429, CBLL1, ZC3H13) are involved in the methylation process, and erasers (FTO, ALKBH5) are involved in the demethylation process. Readers (YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, ZNF217, RBMX) identify RNA methylation modification of information, and participate in the process of translation, degradation of RNA downstream. Regulatory gene interactions affect m6A regulation . For example, CircNSUN2, which is formed by mRNA modification, can bind to intracellular m6A reader YTHDC1 and be regulated by YTHDC1 to determine the nuclear localization. When removed to the cytoplasm, CircNSUN2 bind to reader IGF2BP2, combine with downstream HMGA2 mRNA to improve the stability of HMGA2 mRNA, and ultimately promote liver metastasis of colorectal cancer tumors .
M6A is widely involved in cell proliferation and differentiation, immunity, tumor generation and metastasis, and other life processes. Its poor regulation affects the normal life process, resulting in decreased cell proliferation, immune changes and cell carcinogenesis [10, 11]. Previous studies have found that m6A regulates genes and plays a key role in blood development . Peng et al. 's study found that inhibition of the FTO-FOXO1 pathway resulted in weight loss and decreased blood sugar, which could be considered for the treatment of the metabolic syndrome . M6A gene regulation is closely related to cell proliferation and malignant transformation, especially the METTL3 which recruited downstream translation initiation factors by identifying GGAC sequences, rather than by affecting methylation pathways, to affect translation, thus increasing the expression of oncogenes and affecting the occurrence of cancer [14, 15]. A study has shown that METTL3 is highly expressed in NSCLC tissues, and the expression level of METTL3 is positively correlated with the expression of MiR-33a in NSCLC tissues . METTL3 promotes the development of colorectal cancer , gastric cancer , bladder cancer , breast cancer , renal cell carcinoma , and pancreatic cancer cells , and the occurrence of non-small cell lung cancer is related to METTL3 . The decrease of the METTL3 level can promote the apoptosis of cancer cells through MiR-33a. Previous studies have reported that changing the expression of the METTL3 gene to increase the m6A regulatory genes level in U251 cells affects the occurrence and development of glioma . Moreover, many studies have indicated that m6A regulatory gene greatly affects tumor proliferation. Kwok et al. found that mutations in m6A regulatory gene are closely related to the presence of TP53 mutations in acute myeloid leukemia patients . A study of Ji et al. showed that patients with prostate cancer can be found to have high expression of reader protein and methyltransferase complex and poor survival benefit . Tang et al. suggested that m6A RNA methylation regulator plays a critical role in the potential malignant progression and prognostic value of uveal melanoma . The study of Meng et al. indicated that overall survival was worse in pancreatic cancer patients with m6A regulatory gene change . However, little is known about the relationship between m6A-related genes and endocrine tumors. We investigated and evaluated different clinical features of m6A regulating tumor-related genes, and found the relationship between the expression changes and the clinical pathology of endocrine system tumors, including survival, provide a reference for subsequent research guidance value.
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.
|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.|
M6A regulatory genes have been found to be associated with malignant transformation and development of tumors , but its role in endocrine tumors remains unclear. Due to the heterogeneity of each tumor, we found that in four tumors, only a few m6A regulatory genes showed no difference in the expression of this gene, while most m6A regulatory genes showed differences in the expression of endocrine system tumors compared with normal tissues (Figure 1). All above prove that the expression of regulatory factors of m6A regulatory genes is closely related to endocrine system tumors.
Studies have shown that the core proteins of the m6A methyltransferase do not exist in isolation, they can function by forming complexes . For example, WTAP plays an important role in recruiting METTL3 and METTL14 . A study of Ries et al. showed that m6A modification can promote the co-occurrence and separation of mRNA and binding protein YTHDF protein, and the separation from YTHDF protein can in turn regulate the localization and translation efficiency of mRNA . We then investigated the expression correlation of the m6A regulatory genes. It was found that there was a positive correlation between WTAP and METTL3 expression in each dataset. When the WTAP gene is up-regulated, the METTL3 gene may also be up-regulated. There was an article to explain the effect of METTL3 on WTAP, and the upregulation of METTL3 can lead to an increase of the WTAP protein level in leukemia, causing the carcinogenic effect of leukemia . In addition to the positive correlation between ALKBH5 and FTO expression in THYM (Figure 1B), ALKBH5 showed a negative correlation with other genes. Some articles have demonstrated the role of ALKBH5 in various cancers, such as pancreatic cancer , breast cancer , and gastric cancer . However, ALKBH5 in colorectal cancer  plays an anti-cancer effect. These data also illustrate the interaction of gene expression level between m6A regulatory genes.
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, 2E), we found no significant difference of clinical features.
Then we conducted cluster analysis based on the correlation of expression of m6A regulatory genes in endocrine system tumors, and analyzed ACC (Figure 2D) findings that cluster 1 subgroup was significantly correlated with poor survival status (P < 0.01), higher who grade (P < 0.01), higher T degree (P < 0.01) and poor prognosis. The cluster 1 subgroup of THCA (Figure 2D) was found to be correlated with higher who grade (P < 0.01) and higher T degree (P < 0.01). It may still take a long follow-up to find the difference in OS in THCA dataset. We observed that 'Total 1” (dataset of ACC and THCA) found that its cluster 1 subgroup (Figure 2D) was associated with poorer survival status (P < 0.01), higher who grade (P < 0.01), greater T degree (P < 0.01), more lymph node metastasis (P < 0.01) and poorer prognosis. The above conclusions are sufficient to indicate that m6A regulatory genes are closely related to the clinicopathological characteristics of endocrine system tumors.
A study found that the m6A regulatory gene not only causes tumor progression, but also leads to tumor suppression . METTL14 is an anti-metastasis factor in primary Hepatocellular carcinoma (HCC), and down-regulation of METTL14 suggests poor prognosis in HCC patients and leads to HCC progression and metastasis . However, METTL14 plays a carcinogenic role in endometrial cancer by regulating AKT activation through m6A modifications . According to the gene expression and its clinical factors, we observed that METTL14 is the protective gene in ACC dataset (Figure 4A). The latest literature shows that IGF2BP plays an important role in the development of tumors . We found that IGF2BP3 is closely related to the prognosis of PCPG (Figure 4A), and is risk gene. It has been reported that YTHDF3 plays an important role in m(6) A-modified transcripts , and FTO is also extremely important for tumorigenesis . It is showed that YTHDF3 and FTO are closely related to the prognosis of THCA (Figure 4A), and they are all risk genes. We found that WTAP is protective gene which is closely related to the prognosis of THYM (Figure 4A), “Total 1” data (Figure 4B), and “Total 2” (Figure 4B). Current literature has shown that WTAP plays an important role in inhibiting cancer in gastric cancer , ovarian cancer , and bladder cancer , which is different from the results of this paper, possibly due to tumor heterogeneity. We obtained m6A regulatory genes closely related to OS through the analysis of “Total 1” data (Figure 4B), including METTL14, WTAP, RBM15B, METTL16, FTO and IGF2BP1. Among them, IGF2BP1 is a risk gene. In “Total 2” (Figure 4B), WTAP and METTL16 have a significant correlation with OS, and they are all protective genes. In general, it can be considered that WTAP, METTL16 and IGF2BP3 play important roles in the development of endocrine system tumors, WTAP, METTL16 are protective genes, IGF2BP3 is danger genes.
We also conducted risk assessment, univariate and multivariate Cox regression analysis, and T status (P <0.01) and risk scores (P =0.01) were significantly correlated with OS of ACC (Table 1). The risk scores (P =0.002) remained significantly correlated with OS of PCGA (Table 1). The risk scores remained significantly correlated with OS of THYM (Table 1). Age (P <0.01) was significantly correlated with the OS of THCA (Table 1). In “Total 1” (Table 1), the risk scores (P <0.001), grade (P =0.002), and T status (P =0.012) remained significantly correlated with OS. In “Total 2” (Table 1), the risk scores (P =0.013) remained significantly correlated with OS. In general, we found that gender was not associated with the prognosis of endocrine tumors, and risk scores could be used to greatly predict the prognosis of patients with endocrine system tumors.
In conclusion, our findings systematically demonstrate that m6A regulatory genes expression is closely related to endocrine system tumors, though the expressions of m6A regulatory genes are heterogeneous in the endocrine system tumors. We not only analyzed the correlation between tumors and m6A regulatory genes, but also studied the common characteristics between tumors and m6A regulatory genes in endocrine system tumors, analyzed the clinical prognostic factors of m6A regulatory genes based on its interactivity and expression characteristics. We observed that WTAP, METTL16 and IGF2BP3 play important roles in the development of endocrine system tumors, WTAP, METTL16 are protective genes, IGF2BP3 is danger genes. We also found that gender was not relative with the prognosis of endocrine tumors, and risk score is a prognosis prediction of patients with endocrine system tumors. In summary, WTAP, METTL16 and IGF2BP3 can be regarded as a new promising biomarker, and risk score can be used for the prognosis of UM and the formulation of treatment strategies.
Materials and Methods
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  from the TCGA database and compared the relationship between the expression of 23 m6A-related genes and the clinicopathological variables of endocrine system tumors.
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.
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 . Risk score calculation is as follows:
where Coefi is the coefficient, and xi is the expression value of each selected gene. This formula was used to calculate the risk scores for each patient.
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.
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.
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.
- 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. m6A 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 m6A 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. N6-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 m6A 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. N6-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 m6A 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. m6A 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 m6A-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 m6A 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 m6A-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 m6A 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 N6-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. m6A 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. m6A 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 m6A-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 m6A 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 N6-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. m6A 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 N6-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 m6A 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. m6A-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 m6A-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 m6A 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 m6A 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]