Expression patterns and prognostic value of m6A RNA methylation regulators in adrenocortical carcinoma

Abstract Adrenocortical carcinoma (ACC) is considered a rare cancer with poor prognosis. We used public datasets from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases to assess the relationships between N6-methyladenosine (m6A)-related genes and ACC. We used the Wilcoxon signed-rank test to compare m6A-related gene expression in ACC tissues with that in normal tissues. Then, ACC patients were grouped based on a cluster analysis of m6A-related gene expression. m6A-related genes that were significantly associated with survival were incorporated into a risk signature, and 2 groups were divided according to median risk score. Fisher exact tests were utilized to analyze differences in clinical variables between groups. We compared the overall survival (OS) rates of the groups by means of Kaplan–Meier curves and Cox regression analyses. We found that RBM15, ZC3H3, YTDHF1, YTDHF2, and ALBH5 were overexpressed in ACC and that KIAA1429, YTHDC1, HNRNPC, WTAP, METTL3, and FTO were down regulated in ACC. In addition, membership in cluster 2 or the high-risk group was associated with advanced clinical factors and poor prognosis. The univariable and multivariable Cox regression analyses showed that risk score can be considered an independent prognostic factor for ACC. We found that the expression of m6A-related genes could be used as an independent prognostic factor in ACC. However, the current study has some limitations, and further studies of m6A-related genes in ACC are needed.


Introduction
Adrenocortical carcinoma (ACC) is considered a rare cancer, having an estimated annual incidence of 2 per one million people. [1] However, the invasiveness of ACC is high, being second only to anaplastic thyroid cancer in endocrine carcinoma, and ACC patients with metastasis have poor prognosis. [1][2][3] The median overall survival (OS) and 5-year survival rate of ACC patients are 3.21 years and 15% to 44%, respectively. [4] ACC exhibits a bimodal age distribution, with the first peak occurring at 5 years old and the second peak occurring at 40 to 60 years old. [5] Radical resection is considered the first-line treatment for ACC. In recent years, many hub genes in ACC have been identified and confirmed to be associated with prognosis, suggesting their potential as therapeutic targets. [6][7][8] Much research has been conducted on N6-methyladenosine (m 6 A), the most common RNA modification in eukaryotic mRNAs and lncRNAs. [9] M 6 A methylation is associated with multiple RNA biofunctions and is dynamically regulated via a variety of genes known as "writers", "erasers," and "readers". [10] The "writers" (methyltransferases) catalyze m 6 A modification and upregulate the levels of m 6 A methylation. These genes include Wilms tumor 1-associated protein (WTAP), methyltransferase-like 3 (METTL3), vir like m 6 A methyltransferase associated (KIAA1429/VIRMA), zinc finger CCCH-type containing 13 (ZC3H13), RNA binding motif protein 15 (RBM15) and methyltransferase-like 14 (METTL14). [11] "Eras-Editor: Ning Zhang.
This study did not require ethical approval. We will disseminate our findings by publishing results in a peer-reviewed journal.
ers" (demethylases) include fat mass and obesity-associated protein (FTO) and alkB homolog 5 (ALKBH5), which downregulate the level of m 6 A methylation. [12] "Readers" recognize methylation sites, bind to the RNA and carry out biological functions, [13] they include YTH N6-methyladenosine RNA binding proteins 1 and 2 (YTHDF1 and YTHDF2), YTH domain containing 1 and 2 (YTHDC1 and YTHDC2), and heterogeneous nuclear ribonucleoprotein C (HNRNPC). [14] An increasing number of studies have suggested that the m 6 A modification plays a vital role in various malignancies. For example, the upregulation of METTL3 has been found to inhibit the proliferation, migration and invasion of colorectal cancer cells by p38/ERK pathways [15] and predict poor prognosis in hepatocellular carcinoma. [16] In cervical squamous cell carcinoma, FTO expression was found to be upregulated in tumor tissues and to regulate chemoradiotherapy resistance by targeting betacatenin. [17] However, the relationships between m 6 A-related genes and ACC remain unknown. Hence, we assessed m 6 Arelated gene expression in ACC and evaluated the correlations between m 6 A modifications and ACC prognosis using public datasets. Gene expression data and clinical data for ACC were downloaded from The Cancer Genome Atlas (TCGA), and gene expression data for normal adrenal gland tissue were obtained from the Genotype-Tissue Expression (GTEx) databases.

Patient information
We downloaded gene expression and clinical data of ACC patients from TCGA (https://portal.gdc.cancer.gov/). RNA-seq transcriptome data for normal adrenal tissue were obtained from GTEx (https://www.gtexportal.org/). We normalized and merged the transcriptome data from 79 ACC tissues and 127 normal adrenal tissues from the 2 databases into a single dataset. We used the single dataset to screen differentially expressed m 6 A-related genes. Among the 79 ACC tissues, 77 patients had both gene expression and clinical information, which were included for further prognosis analysis. Information on the 77 ACC patients is shown in Table 1. Then, 13 m 6 A-related regulatory genes described in the published literature were selected for analysis  (METTL3, METTL14, WTAP, KIAA1429, RBM15, ZC3H13,  YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC, FTO,  and ALKBH5). [10][11][12] Additionally, copy number variation (CNV) data of 90 ACC tissues and 90 normal tissues were also obtained in TCGA database.

Bioinformatic analysis
Based on the expression levels of m 6 A-related genes, cluster analysis was performed to cluster the patients of ACC into different groups using the Consensus Cluster Plus package (http:// www.bioconductor.org/packages/ConsensusClusterPlus) of R software. [18] We then utilized principal component analysis to observe the patterns of m 6 A-related gene expression in the groups.
To analyze the relationships between m 6 A-related genes and ACC prognosis, we performed univariable Cox regression and identified the genes that were closely related to survival (P < .05). Then, the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was performed to determine coefficients of each selected gene. All Cox regressions were completed via the survival package (https://cran.r-project.org/ web/packages/survial) of R software. We used the coefficients to calculate risk scores according to the following formula to build a potential risk signature: Where xi is the relative expression of genes closely related to survival. [19] The receiver operating characteristic (ROC) curve was utilized to validate the accuracy of risk score using the survival ROC package (https://cran.r-project.org/web/packages/ survivalROC) of R software.

Statistical analysis
To compare the expression of m 6 A-related genes in ACC tissues with that in normal tissues, we conducted Wilcoxon signed-rank tests using the limma package (http://www.bioconductor.org/ packages/limma) of R. Correlations between m 6 A-related genes were determined using the corrplot package (https://cran.rproject.org/web/packages/corrplot) of R. Then, ACC patients were grouped based on the cluster analysis or risk score (with the median value of risk score considered the cut-off value). Chisquare tests were used to evaluate the differences in CNV between tumor tissues and normal tissues. Fisher exact tests were used to evaluate the differences in clinical variables between clusters and between risk-score groups.
To compare OS between cluster groups or risk groups, Kaplan-Meier survival curves were generated and compared with the log-rank test via the survival package (https://cran.rproject.org/web/packages/survial) of R. Next, we performed univariable and multivariable Cox regression analyses to evaluate whether risk score is an independent predictor of poor OS in ACC patients. Before assessing for multivariable Cox regression, proportional hazard assumption for each covariate was tested via in minus in survival curves. By observing the in minus in survival curves of each covariate in 2 states (Age 48 vs Table 1 Characteristics of the included ACC patients obtained from the TCGA database.

m 6 A-related gene expression in ACC and normal tissues
To compare m 6 A-related gene expression in the 79 ACC tissues with that in the 127 normal tissues, we used the Wilcoxon signedrank test. RBM15, ZC3H3, YTDHF1, YTDHF2 and ALBH5 were overexpressed in ACC tissues relative to normal tissues (all P values <.001). The expression levels of KIAA1429, YTHDC1, HNRNPC, WTAP, METTL3, and FTO were downregulated in ACC tissues (all P values <.001). All the results were shown by heatmap, violin plot and volcano plot. (Fig. 1A-C). We then explored the correlations between m 6 A-related genes. Among all pairwise combinations of m 6 A-related genes, HNRNPC and METTL3 had the highest positive correlation of expression, and RBM15 and METTL3 had the highest negative correlation (Fig. 1D). The results of CNV analysis showed that compared with normal adrenal gland tissues, YTHDF2 was single deleted on chromosome 1, YTHDC1 and METTL14 were single gained on chromosome 4, YTHDC2 was single gained on chromosome 5, KIAA1429 was single gained on chromosome 8, FTO was single gained on chromosome 16, ALKBH5 was single deleted on chromosome 17, and YTHDF1 was single gained on chromosome 20 (Fig. 1E).

Grouping based on the cluster analysis was significantly correlated with prognosis
We used cluster analysis to cluster the ACC patients into different groups based on m 6 A-related gene expression ( Fig. 2A-D), ultimately dividing them into 2 groups: cluster 1 and cluster 2.
Principal component analysis was utilized to observe the patterns  (Fig. 2E), which were unambiguously distributed between cluster 1 and cluster 2.
Fisher exact tests showed that cluster 2 was significantly associated with advanced clinical stage (Fig. 2F) (Table 2). Additionally, the Kaplan-Meier survival curves showed that the prognosis of patients with ACC in cluster 2 was worse than that of patients in cluster 1 (P = .004) (Fig. 2G).

Grouping based on risk score was strongly associated with prognosis
We used univariable Cox regression to analyze the relationship between the OS of ACC patients and the expression of m 6 Arelated genes (Fig. 3A). The expression levels of 4 genes, namely, HNRNPC, RBM15, METTL14, and FTO, were found to be significantly associated with prognosis. We then applied LASSO Cox regression to establish gene coefficients. The results showed that the risk signature could be constructed using all 4 genes, and risk scores were calculated based on the coefficients from the LASSO Cox regression (Fig. 3B-C). Two groups of patients were established according to median risk score. The ROC curve results showed that the OS of ACC patients was perfectly predicted by risk score (AUC = 0.865) (Fig. 3D). The Fisher exact tests revealed that membership in the high-risk group of ACC patients was significantly associated with advanced clinical stage, high M classification and high T classification (Fig. 3E) ( Table 3). The Kaplan-Meier survival curves revealed that membership in the high-risk group of ACC patients implied poor prognosis (P < .001) (Fig. 3F). Univariable Cox regression showed that high M classification, high T classification, advanced clinical stage and high risk score were strongly correlated with poor prognosis (Fig. 3G). We then performed multivariable Cox regression, which identified risk score and T classification were associated with prognosis (Fig. 3H). This result indicates that a high risk  score is an independent predictor of poor OS in ACC. Before assessing for multivariable Cox regression, the 2 survival curves of each covariate were parallel, indicating that the proportional hazard assumption was true (Supplementary Fig. 1A-G, http:// links.lww.com/MD/F841).

Discussion
With the development of epigenetics, m 6 A methylation has become a promising research direction and has been demonstrated to be closely associated with the occurrence and development of tumors. [20][21][22][23] The overexpression of a "writer" of m 6 A, METTL3, was correlated with the poor prognosis of bladder carcinoma by positively modulating pri-miR221/222. [24] In pancreatic cancer cells, the knockdown of METTL3 was found to decrease chemo-and radio-resistance, [25] and increased expression of METTL3 was found to promote the invasion and migration of melanoma cells via the actions of matrix metallopeptidase 2. [26] A "reader" of m 6 A, YTHDF1, was found to be overexpressed in colorectal cancer and to stimulate the invasion of colorectal cancer cells by activating the Wnt/beta-catenin pathway. [27] In acute myeloid leukemia, the upregulation of YTHDF2 was found to decrease the half-life of diverse m 6 A transcripts, and YTHDF2 knockdown was shown to inhibit cancer stem cell development. [28] Furthermore, an "eraser" of m 6 A, FTO, was shown to be suppressed in clear cell renal cell carcinoma (ccRCC) tissues, and overexpression of FTO was found to activate oxidative stress and ROS production by increasing the expression of PGC-1a. [29] Although an increasing number of studies have confirmed that m 6 A modification is significantly associated with tumors, the relationships between m 6 A modification and ACC have to date been unknown.
Here, we compared the expression of m 6 A-related genes between ACC tissues and normal adrenal gland tissues. In ACC relative to normal tissue, RBM15, ZC3H3, YTDHF1, YTDHF2, and ALBH5 were overexpressed, and the expression levels of KIAA1429, YTHDC1, HNRNPC, WTAP, METTL3, and FTO in ACC were low. Based on the cluster analysis results, the ACC patients were grouped into clusters 1 and 2. The prognosis of cluster 2, which was associated with advanced clinical stage and high T classification, was worse than that of cluster 1. We constructed a risk signature through LASSO Cox regression based on 4 m 6 A genes identified as associated with prognosis by univariable Cox regression and calculated risk scores. The prognosis of the high-risk group, which was significantly related to advanced clinical stage, higher T classification and M classification, was worse than that of the low-risk group. The univariable and multivariable Cox regression results showed that the risk score of m 6 A modification can be used as an independent prognostic factor in ACC. Taken together, our results demonstrate that the expression of m 6 A-related genes is closely related to clinical variables and prognosis in ACC. These results indicate that expression of m 6 A-related genes can serve as an independent prognostic factor in ACC and that such genes offer potential new drug targets in ACC, as targeting m 6 A may be an effective treatment for tumors. [30] The biological functions of m 6 A related genes should also be focused. M 6 A methyltransferase complex ("writers") could mediate m 6 A methylation of RNAs, a modification that plays a role in the efficiency of mRNA splicing and RNA processing. Previous studies have found that METTL3 and METTL14 formed stable heterodimers, which interacted with WTP to constitute WMM complex (WTP-METTL3-METTL14), which promoted the deposition of m 6 A. [31] Recent in-depth studies indicated that KIAA1429 led to the preferential mRNA methylation in 3 0 UTR and near the stop codon, then recruited WMM as the core to guide m 6 A methylation at specific sites. [32] As an RNA binding protein, RBM15 was involved in the regulation of hematopoietic cell homeostasis and mRNA alternative splicing and recruited WMM complex to specific RNA sites during the m 6 A methylation. [33] ZC3H13 was considered to be a component of WMM complex, which could improve the catalytic function of WMM complex by interacting with WTAP. [34,35] As demethylase, the "erasers"(FTO and ALKBH5) used ferrous as cofactor, and a-ketoglutarate was used as cosubstrate to remove the methylation of m 6 A. [36] FTO promotes adipocyte cycle progression and adipogenesis by reducing the m 6 A level of cyclin A2 and cyclin dependent kinase 2, and controls mRNA splicing by inhibiting serine and arginine rich splicing factor 2 binding at splicing sites. [37,38] Demethylation induced by ALKBH5 was related to splicing and the production of long 3 '-UTR mRNA. [39,40] As "readers" of m 6 A methylation (YTHDC1, YTHDC2, YTHDF1, YTHDF2 and HNRNPC), the function was to recognize and bind to the m 6 A site. [41,42] YTHDF1 could promote the synthesis of protein and translation of mRNA. [43] By binding to m 6 A-modified mRNA, YTHDF2 activated transcription product degradation. [44] YTHDC1 was significantly associated with RNA splicing and export. [45,46] YTHDC2 enhanced the target RNAs translation efficiency, but decreased their abundance. [47] HNRNPC was involved in regulating mRNA splicing. [48] The current study has some limitations. The number of ACC tissues was low, and paired normal tissue controls were not available in the TCGA database; thus, we had to pool the data from 2 databases. In addition, the protein levels of m 6 A-related genes in ACC could not be reliably estimated; thus, we focused on the relationship between m 6 A-related gene expression and ACC prognosis. Furthermore, the biological functions and mechanisms of the m 6 A-related genes were not explored. Moreover, previous studies have shown that other genes with abnormal expression in ACC can affect prognosis in ACC. For example, the upregulation of urothelial carcinoma-associated 1 (UCA1) promoted the proliferation and inhibited the apoptosis of ACC cells though the miR-298-CDK6 axis. [49] Overexpression of HOX transcript antisense RNA (HOTAIR), by regulating the cell cycle, led to the poor OS in ACC. [50] High expression of topoisomerase alpha 2 (TOP2A) was demonstrated to activate ACC cell proliferation, with TOP2A representing a potential therapeutic target. [51] We tried to explore the relationships between m 6 A-related genes and these ACC-driver genes. But the relationships could not be clearly determined due to a lack of studies. Moreover, because the mechanisms of m 6 A regulating the development of tumors were very complex, we admitted that the researches on the related mechanisms were not deep enough, and no basic experiment were performed, which also became a limitation of our article. Therefore, further studies of m 6 A-related genes in ACC are needed.

Conclusions
In the current study, we compared m 6 A-related gene expression between ACC tissues and normal adrenal gland tissues. We found that the expression of m 6 A-related genes could be used as an independent prognostic factor in ACC. However, the current study has some limitations, and further studies of m 6 A-related genes in ACC are needed.  Differences in the characteristics of ACC patients between the high risk and low risk.

Basic information
Low risk High risk P value