The Prognostic Values of m6A RNA Methylation Regulators in Thyroid Carcinoma

Background: N6-methyladenosine (m6A) is the most common RNA modication and regulates RNA splicing, translation, translocation, and stability. Aberrant expression of m6A has been reported in various types of human cancers. m6A RNA modication is dynamically and reversibly mediated by different regulators, including methyltransferase, demethylases, and m6A binding proteins. However, the role of m6A RNA methylation regulators in thyroid cancer remains unknown. The aim of this study is to investigate the effect of the 13 main m6A RNA modication regulators in thyroid carcinoma. Methods: The gene expression prole of m6A RNA modication regulations and clinical information of patients with thyroid carcinoma were obtained from The Cancer Genome Atlas database. Consensus clustering was applied to identify two clusters of thyroid carcinomas with different clinical outcome. LASSO Cox regression analysis was used to construct gene-based prognostic signature based on the expression of m6A RNA methylation regulators. Kyoto Encyclopedia of Genes and Genomes, Gene Ontology and gene set enrichment analyses were performed to explore differential cellular processes and signaling pathways between the two groups based on risk signature. Results: We found that most of the m6A RNA modication regulators are down-regulated in 450 patients with thyroid carcinoma. We identied two clusters based on the gene expression proles of 13 m6A RNA modication regulators using consensus clustering. The cluster 2 subgroup correlates with an unfavorable outcome compared with the cluster 1 subgroup. In addition, we derived a three m6A RNA modication regulator genes-based risk signature (FTO, RBM15 and KIAA1429), that is an independent prognostic biomarker in patients with thyroid carcinoma. There were signicantly different signaling pathways between high and low risk group by Kyoto Encyclopedia of Genes and Genomes, Gene Ontology and gene set enrichment analyses. Moreover, we found that this risk signature could better predict outcome in male than female. Conclusion: Our study revealed the prognostic value of m6A RNA methylation regulators in patients with thyroid carcinoma. Via Cox univariate analysis, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was employed to construct a prognostic gene signature with three m6A methylation regulatory genes. The differences in the biological pathways were also explored by Kyoto Encyclopedia Genes and Genomes (KEGG), Gene Ontology (GO) and gene set enrichment analysis (GSEA) between subgroups with different prognosis. Our results emphasized the signicance of m6A methylation-related genes in thyroid cancer and established a prognostic gene signature for predicting the prognosis in patients with thyroid cancer.


Introduction
Thyroid carcinoma is one of the most prevalent human malignancies, and its incidence has signi cantly increased over the past ten years [1]. According to the histological characteristics, thyroid carcinoma can be classi ed into four types, including anaplastic, follicular, medullary, and papillary thyroid carcinoma [2]. Although thyroid cancer has been considered a disease with favorable outcome, the patients with advanced thyroid carcinoma have a poor 5-year survival rate [3][4]. Therefore, it is important to identify effective diagnostic markers, prognostic indicator and therapeutic targets for thyroid carcinoma treatment.
N6-methyladenosine (m6A), a modi cation occurs at the N6 position of RNA adenine (A), is the most abundant modi cation on RNA molecules. m6A was rst discovered in 1974, occurring in various types of RNAs, including mRNA, tRNA, snRNA, and long noncoding RNA [5][6]. The m6A methylation is the most prevalent and well recognized modi cation among all RNA modi cations. Similar to DNA methylation, m6A methylation regulates the post-transcriptional expression of target RNAs by affecting RNA alternative splicing, stability, nuclear export, mRNA decay, and translation [7][8][9][10]. Its regulatory effects are mediated by the dynamic and reversible interactions among methyltransferases ("writers"), demethylase ("erasers") and RNA-binding proteins ("readers"), which can add, remove, and recognize m6A-modi ed sites and generate different functions [11]. Methyltransferases, known as the "writers", are composed of METTL3/14/, RBM15, WTAP, KIAA1429, and ZC3H13, regulating the process of methylated modi cation of RNA. Demethylase, which are considered m6A "erasers", consist of FTO and ALKBH5, mediating the demethylation process of RNA. RNA-binding proteins, known as m6A "readers", recognize RNA methylated information and participate in the translation and degradation of downstream RNA, including YTHDC1, YTHDC2, YTHDF1, YTHDF2, and HNRNPC. Increasing studies have shown that abnormal m6A modi cation plays an important role during tumorigenesis and progression in a variety types of human cancers, including oral squamous cell carcinoma, acute myeloid leukemia (AML), glioblastoma, hepatocellular carcinoma, pancreatic cancer and breast cancer [12][13][14][15][16][17][18]. However, the association between m6A RNA modi cation regulators and thyroid carcinoma remains unknown.
In present study, we identi ed that 13 m6A RNA methylation regulators were abnormally expressed in thyroid carcinoma tissues compared to those in normal tissues from The Cancer Genome Atlas (TCGA) thyroid carcinoma database. Based on the signature of 13 m6A genes, we divided patients into two clusters by consensus clustering method and compared the clinical outcomes of two clusters. The Via Cox univariate analysis, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was employed to construct a prognostic gene signature with three m6A methylation regulatory genes. The differences in the biological pathways were also explored by Kyoto Encyclopedia Genes and Genomes (KEGG), Gene Ontology (GO) and gene set enrichment analysis (GSEA) between subgroups with different prognosis. Our results emphasized the signi cance of m6A methylation-related genes in thyroid cancer and established a prognostic gene signature for predicting the prognosis in patients with thyroid cancer.

Dataset
All data were obtained from the TCGA database (https://portal.gdc.cancer.gov/), containing RNA-seq transcriptome data and the corresponding clinicopathological information on 450 thyroid carcinoma patients (450 cases of tumor tissues and 58 cases of normal tissues; Table S1).
Differential expression analysis of m6A RNA methylation regulators Limma package and wilcox-test were used for screening the differentially expressed genes (DEGs) between tumor tissues and normal tissues. We extracted the expression matrix of the 13 known m6A RNA methylation regulators, including METTL3, YTHDF1, KIAA1429, YTHDC2, ALKBH5, YTHDF2, YTHDC1, ZC3H13, METTL14, FTO, WTAP, RBM15 and HNRNPC. Heatmap and vioplot were used to display the different expression of the 13 m6A RNA methylation regulators in normal tissues and tumor tissues. Pearson correlation analysis was used to evaluate the association among different m6A RNA methylation regulators with "corrplot" R package.

Consensus clustering analysis
The TCGA thyroid carcinoma cohort was divided into different clusters by consensus expression of m6A RNA methylation regulators with "ConsensusClusterPlus" R package [19]. Principal component analysis (PCA) was performed by limma and ggplot2 package. Kaplan-Meier survival curves and the log-rank test were used to evaluate patients of thyroid carcinoma with different clusters. Chi-square test was used to compare the distribution of age, gender and stage between different clusters.

Prognostic signatures construction and prediction
Univariate Cox analysis was performed to assess the prognostic value of m6A RNA methylation regulators. Regulators with a hazard ratio (HR) < 1 or > 1 were considered as protective or risk genes, respectively. The LASSO Cox regression algorithm was used to construct the optimal prognostic model out of the 13 m6A regulators. A risk score for each patient was calculated as the sum of each gene's score, which was obtained by multiplying the expression of each gene and its coe cient. We divided the patients from the TCGA thyroid carcinoma cohort into high-risk and low-risk groups based on their risk scores being below or above the median value. The survival difference in overall survival between highrisk group and low-risk group was analyzed by the Kaplan-Meier method with a two-sided log-rank test.
Receiver operating characteristic (ROC) curve was constructed to evaluate the accuracy of the model prediction. Chi-square test was performed to compare the differences of clinicopathological features between high-risk group and low-risk group. Heatmaps were generated using the pheatmap R package and used to visualize the differences. Univariate and multivariate Cox regression analyses were used to identify the independent prognostic factors for the TCGA thyroid carcinoma cohort. The survival difference between high-risk group and low-risk group subgrouped by age, gender, tumor size and stage was further evaluated.
Signaling pathways and cellular processes affected by m6A regulators GSEA was used to explore the biological pathways underlying the high-risk and low-risk groups de ned by the prognostic signature in the thyroid carcinoma patients. KEGG gene sets (v7. 0) and phenotype label (HighRisk vs. LowRisk) les were generated and loaded into the GSEA software (v4.0.3; Broad Institute, Cambridge, MA). The permutation test run 1, 000 times. The GO and KEGG pathway enrichment analyses were applied for exploring the relevant signaling pathways enriched by DEGs between the highrisk group and the low-risk group using clusterPro ler R package.

Statistical analysis
Data analysis was performed with the R v3.6.1. All statistical tests were two-sided. A P value of less than 0.05 was considered statistically signi cant.

Patient demographics and tumor characteristics
The thyroid carcinoma cohort was downloaded from TCGA and the clinicopathological characteristics of patients were listed in Table S1. Only cases with clinicopathological and prognostic information were included in this study (n = 450).
The landscape of m6A RNA methylation regulators in thyroid carcinoma To investigate the role of m6A RNA methylation regulators in thyroid carcinoma development and progression, we rst determined the expression levels of 13 m6A RNA methylation regulators 450 thyroid carcinoma tissues and 58 normal thyroid tissues from TCGA dataset by heatmap (Fig. 1A). We observed that the expression level of 11 m6A RNA methylation regulators (METTL3, METTL14, RBM15, WTAP, YTHDC1, YTHDC2, YTHDF1, KIAA1429, FTO, ZC3H13 and ALKBH5) was signi cantly decreased in thyroid carcinoma tissues compared to those in normal thyroid tissues, whereas the expression level of HNRNPC was increased in thyroid carcinoma tissues (Fig. 1B). The alternation of m6A RNA methylation regulators ratio is an intrinsic feature that can characterize individual differences [20]. The correlations between 13 m6A RNA methylation regulators were analyzed by Spearman correlation analysis. As shown in Fig. 1C, the relationship between the 13 m6A RNA methylation regulators was weakly to strongly positively correlated. Moreover, the METTL14 gene was most strongly correlated with KIAA1429 and YTHDC1 genes (r = 0.78), and the ZC3H13 gene and the KIAA1429 gene were also most relevant (r = 0.78). In addition, we investigated the association between each individual m6A RNA methylation regulator and the clinicopathological factors of patients with thyroid cancer, including age, gender, stage, tumor size, lymph node status and metastasis and found a robust relationship between m6A RNA methylation regulator and clinicopathological factors of patients with thyroid cancer (Fig. 1D).
Consensus clustering of m6A RNA methylation regulators identi ed two subgroups of thyroid carcinoma To provide insight into the clinical relevance of m6A RNA methylation regulators in thyroid carcinoma, we divided the patients into subgroups based on the expression of m6A RNA methylation regulators . "k" was used to represent the number of subgroups ( Fig. 2A and 2B). After consensus clustering analysis, k = 2 was regarded as the most appropriated selection to separate the patients with thyroid carcinoma into two subgroups, namely cluster 1 (n = 231) and cluster 2 (n = 219) (Fig. 2C). To verify our classi cation, we analyzed the two subgroups by PCA. The results indicated that both cluster 1 subgroup and cluster 2 subgroup could be gathered together (Fig. 2D). Furthermore, we observed that the cluster 1 subgroup had a signi cantly favorable outcome compared with that in the cluster 2 subgroup (P = 0.013; Fig. 2E). There were signi cant differences in clinicopathological characteristics between the two groups, including gender, tumor status and lymph node status (Fig. 2F).

Prognostic value of risk signature
To better understand the prognostic role of the m6A RNA methylation regulators in thyroid carcinoma, we performed a univariate Cox regression analysis on the expression level of the 13 m6A RNA methylation regulators genes. As shown in Fig. 3A, high expression of RBM15 (HR = 4.341, 95% CI = 1.328 -14.192, P = 0.015), FTO (HR = 1.623, 95% CI = 1.021 -2.580, P = 0.041) and KIAA1429 (HR = 1.925, 95% CI = 1.079 -3.434, P = 0.026) had a poor overall survival in patients with thyroid carcinoma. We next applied the LASSO Cox regression algorithm to the m6A RNA methylation regulator genes in the TCGA dataset and constructed an optimal prognostic gene signature. Three genes (RBM15, FTO and KIAA1429) were selected to construct the optimal prognostic gene signature based on the minimum criteria, and the coe cients obtained from the LASSO algorithm were used to calculate the risk score for TCGA dataset ( Fig. 3B and 3C). Based on the median risk score, we separated the thyroid carcinoma patients into highand low-risk subgroups. Although there was no signi cant difference in clinicopathological characteristics between the two subgroups (Fig. 3D), the high-risk subgroup had a worse overall survival in patients with thyroid carcinoma (Fig. 3E). The area under the ROC curve (AUC) was also calculated to evaluate the ability of the prediction model (Fig. 3F).

Investigation of biological pathways affected by m6ARNA methylation regulators
We next explored the biological signaling pathways underpinning the different outcome between the highrisk and low-risk subgroups by GSEA analysis. We observed an increased expression of genes involved in  Fig.4A). A total of 214 differentially expressed genes (log FC > 1 and P < 0.05) were identi ed between the two subgroups and the differentially expressed genes were then annotated using GO (Fig. 4B) and KEGG (Fig.  4C) analyses. Several malignancy-related biological processes or signaling pathways were noticeable, including Ras signaling pathway, phosphoinositide 3-kinase (PI3K)-AKT signaling, peroxisome proliferator-activated receptor (PPAR) signaling.
The prognostic risk score was an independent prognostic factor in patients with thyroid carcinoma Next, we performed univariate and multivariate Cox regression analyses to determine whether the risk score is an independent prognostic biomarker in patients with thyroid carcinoma. The univariate analysis showed that the age (P < 0.001, HR = 1.140, 95% CI = 1.081 -1.203), stage (P < 0.001, HR = 2.750, 95% CI = 1.594 -4.744), tumor size (P = 0.003, HR = 2.852, 95% CI = 1.415 -5.746), and risk score (P < 0.001, HR = 1.007, 95% CI = 1.003 -1.011) were signi cantly correlated with the overall survival (Fig. 5A). Multivariate Cox regression analysis was performed using variables signi cant in univariate analysis. As shown in Fig. 5B, the age (P < 0.001, HR = 1.127, 95% CI = 1.063 -1.195) and risk score (P = 0.001, HR = 1.008, 95% CI = 1.003 -1.012) were identi ed as the independent prognostic factors in patients with thyroid carcinoma. We also determined whether the risk score are able to predict clinical prognosis in thyroid carcinoma strati ed by age (Fig. 6A and 6B), gender ( Fig. 6C and 6D), stage ( Fig. 6E and 6F) and tumor size (Fig. 6G and 6H). We found that the risk score could predict outcome better in Male than Female (Fig. 6C and 6D). Discussion m6A, a member of RNA epigenetic modi cation families, is the most abundant mRNA modi cation in most eukaryotes and involves almost all steps of RNA metabolism, including splicing, stability, translocation and translation [21]. Although increasing evidence indicates that the m6A RNA methylation regulators act as both oncogene or tumor suppressor in many types of human malignant tumors by regulating the mRNA expression of tumor-related genes, the speci c regulatory role of m6A RNA methylation regulators in cancer development and progression remains to be clari ed. In the present study, we investigated the role of 13 m6A RNA methylation regulators in thyroid carcinoma and found that the expression of m6A RNA methylation regulators was closely associated with the outcome in patients with thyroid carcinoma.
Analyses of TCGA database revealed that most of m6A RNA methylated regulators except YTHDF2 were dysregulated in thyroid carcinoma, suggesting that m6A RNA modi cation plays a critical role in thyroid cancer development. Among them, only HNRNPC showed an increased expression in thyroid carcinoma tissues compared with that in normal tissues. HNRNPC, as the "readers" of m6A RNA modi cation, is well recognized for its regulatory roles in RNA splicing, sequence-unspeci c RNA exportation, stability and translation. Abnormal increased expression of HNRNPC was observed in multiple types of tumors, including breast cancer, lung cancer, glioblastoma, and hepatocellular carcinoma [22][23][24][25]. Consistent with their work in other cancers, our results suggest that HNRNPC functions as an oncogene in thyroid carcinoma.
We identi ed two subgroups of thyroid carcinoma by consensus clustering based on the expression of the m6A RNA methylation regulators, and the cluster 2 subgroup associated with a poorer outcome. In addition, we constructed a risk signature by using 3 m6A RNA methylation regulators (FTO, RBM15, and KIAA1429) and the risk score is an independent prognostic biomarker in patients with thyroid carcinoma. FTO is the rst mammalian RNA m6A demethylase, which belongs to the AlKB family [26][27]. Evidences have shown that FTO involves in the process of malignant tumors. FTO promotes AML carcinogenesis by enhancing the stability MYC and CEBPA mRNA [28]. In cervical cancer, FTO enhances the chemoradiotherapy-resistance by reducing the m6A modi cation of β-catenin [29]. KIAA1429, known as VIRMA, mainly regulates the m6A mRNA methylation in 3' UTR and near stop codon [30]. It is reported that KIAA1429 promotes breast cancer progression by regulating CDK1 [31]. RBM15 was originally described as a 5' translocation partner of the MAL gene in t (1; 22) (p13; q13) infant acute megakaryocytic leukemia. Thus, RBM15 plays an important role in hematopoietic development, and c-MYC is a potential target of RBM15 in the regulation of adult hematopoietic stem cell and megakaryocyte development [32][33]. In chronic myelogenous leukemia (CML), RBM15 promotes CML carcinogenesis through its effect on the Notch signaling [34]. The functional and mechanistic studies of these genes should be further performed to de nite its role in the carcinogenesis of thyroid carcinoma.
To explain the different prognosis of the two groups, we attempted to explore the mechanisms through GSEA and KEGG analyses. GSEA indicated that the signi cantly enriched gene sets in the high-risk subgroup with poor survival were some crucial signaling pathways and cellular processes, including ubiquitin mediated proteolysis, Wnt signaling pathway, mTOR signaling pathway, TGF-β signaling pathway, pathways in cancer, MAPK signaling pathway, RNA degradation, and cell cycle. KEGG analysis of differentially expressed genes between the two groups indicated that the m6A RNA methylation regulators were involved in several cancer-related signaling pathways, including PI3K-Akt signaling pathway and PPAR signaling pathway.
In conclusion, we determined the expression, potential function and prognostic value of the m6A RNA methylation regulators in thyroid carcinoma. The prognostic risk score was an independent prognostic factor in patients with thyroid carcinoma. Our study provides important evidence for further investigation of the role of m6A RNA methylation in thyroid carcinoma. Figure 1 The   Construction of risk signature with three m6A RNA methylation regulators. A, Univariate Cox analysis of the 13 m6A RNA methylation regulators to identify the genes that signi cantly correlated with overall survival. B and C, Coe cients of three m6A RNA methylation regulators calculated by LASSO Cox Page 15/18 regression algorithm. D, The relationship between the three m6A RNA methylation regulators expression and clinicopathologic characteristics in high-/low-risk subgroup of thyroid carcinoma patients. E, Kaplan-Meier analysis of overall survival in high-/low-risk subgroup of thyroid carcinoma patients. F, The prediction e ciency of risk signature analyzed by ROC curve.

Figure 4
Differential signaling pathways and cellular processes between the high-risk subgroup and the low-risk subgroup de ned by a 3-gene expression signature. A, The highly enriched gene sets in high-risk group analyzed by GSEA. B, Biological processes of differentially expressed genes between the high-risk subgroup and the low-risk subgroup. C, KEGG pathways in which differentially expressed genes were distributed.

Figure 5
Evaluation of prognostic values of risk scores in patients with thyroid carcinoma. A, Univariate Cox analyses of the association between clinicopathological characteristics and overall survival in patients with thyroid carcinoma. B, Multivariate Cox analyses of the association between clinicopathological characteristics and overall survival in patients with thyroid carcinoma.