Prognostic Value of m6A RNA Methylation Modulators and Potential Clinical Application in Pancreatic Cancer

m6A is the most prevalent and abundant form of mRNA modication and plays a dual role in cancer development. The high incidence and mortality of pancreatic cancer are critical obstacles worldwide. In this study, we investigated the function of m6A RNA methylation modulators in pancreatic cancer. Expression of 13 m6A RNA methylation modulators and clinical data from patients with pancreatic adenocarcinoma were obtained from TCGA database. Differences in the expression of 13 m6A RNA methylation modulators between tumour (n = 178) and healthy (n = 4) samples were compared by Wilcoxon test. LASSO Cox regression was used to select m6A RNA methylation modulators for analysis of the relationship between expression and clinical characteristics by univariate and multivariate regression. The pathways of the m6A RNA methylation modulators were examined by gene set enrichment analysis (GSEA) and we found enrichment in chemokine, ribosome, and mTOR signalling pathways.

Abstract Background m6A is the most prevalent and abundant form of mRNA modi cation and plays a dual role in cancer development. The high incidence and mortality of pancreatic cancer are critical obstacles worldwide. In this study, we investigated the function of m6A RNA methylation modulators in pancreatic cancer.

Methods
Expression of 13 m6A RNA methylation modulators and clinical data from patients with pancreatic adenocarcinoma were obtained from TCGA database. Differences in the expression of 13 m6A RNA methylation modulators between tumour (n = 178) and healthy (n = 4) samples were compared by Wilcoxon test. LASSO Cox regression was used to select m6A RNA methylation modulators for analysis of the relationship between expression and clinical characteristics by univariate and multivariate regression. The pathways of the m6A RNA methylation modulators were examined by gene set enrichment analysis (GSEA) and we found enrichment in chemokine, ribosome, and mTOR signalling pathways.

Results
WTAP had a low expression in tumour samples compared with healthy samples. Furthermore, our analyses revealed that the m6A RNA methylation modulators YTHDF1, ALKBH5, METTL3, METTL14, and KIAA1429 correlated with high-risk patients, resulting in an elevated risk score and a lower overall survival. High-risk score correlated with clinical characteristic and was an independent prognostic indicator for pancreatic adenocarcinoma. The pathways involved were identi ed by GSEA to explore the potential mechanism of action.

Conclusion
Modulators involved in m6A RNA methylation were associated with the development of pancreatic cancer. A risk score based on the expression of YTHDF1, ALKBH5, METTL3, METTL14, and KIAA1429 may be an independent prognostic indicator.

Background
Pancreatic cancer is an aggressive cancer of the digestive system with high morbidity and mortality. Data from the National Centre for Health Statistics (NCSH) show a reduction in the occurrence of breast cancer, colorectal cancer, and prostate cancer in the past decade, but the incidence of pancreatic cancer continues to increase and it has become the fourth leading cause of cancer death [1] . Studies have predicted that it could be the second leading cause of cancer death by 2030 [2] . The estimated incidence of pancreatic cancer in the United States in 2020 is around 57,600 with a projected 47,050 deaths, which indicates a mortality above 80% [1] . Furthermore, the 5-year survival rate of pancreatic cancer is 9% and is the lowest across all types of cancer. Pancreatic adenocarcinoma (PAAD) is the most common type of pancreatic cancer, which encompasses up to 85% of pancreatic cancer [2] . While the precise cause and underlying mechanisms of PAAD are still unclear, genetics is an important factor related to tumour development and mutations in DNA or RNA may also drive the initiation of PAAD [3] .
RNA modi cation is a crucial part of epigenetics, which together with gene and protein modi cation play an important role in the regulation of cellular processes. Interest in this eld has been rising in recent year and according to the data from the MODOMICS update in 2017, 163 different modi cations have been identi ed in RNA molecules so far including N1-methyladenosine, N7-methyladenosine, 5-methylcytosine, pseudouridine, N6,2'-O-dimethyladenosine (m6A) and 2'-Omethylation [4] . In eukaryotes, m6A is the most prevalent and abundant form of mRNA modi cation, which was identi ed in the 1970s [5,6] . m6A is a dynamic and reversible process of mRNA modi cation and has been widely studied, although the exact mechanism of m6A modi cation remains unelucidated due to the technical limitations. The formation of m6A is composed of three individual proteins and it is possible to detect these m6A proteins as a re ection of m6A modi cation. The m6A modulator proteins are classi ed into three categories: writers, erasers, and readers. Writer proteins mainly include METTL3, METTL14, and WTAP [7][8][9][10] . RBM15, ZC3H13, and KIAA1429 have been newly identi ed and play a role in mediating the formation of the writer protein complex [11][12][13][14] . These writer proteins transfer the methyl group from S-adenosyl methionine to the RNA nucleotides. Erasers on the other hand are demethylases and remove the methyl group from the RNA molecule. Members of this group include FTO and ALKBH5 [15,16] . Finally, reader proteins the members recognise m6A modi cation and include YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2 [17,18] , and the newly discovered HNRNPC [19] .
m6A plays a vital role in the splicing, translation, and stability of RNA and regulates chromatin state and transcription [20] .
This has signi cant consequences in various human diseases where m6A alterations may enhance heart failure [21] and in uence brain development and function [22] , immune response to viral infection [23] , and bone metabolism [24] . An increasing amount of evidence indicates that m6A RNA methylation modulators are associated with the development and progression of several malignant tumours. In bladder cancer, m6A modulator METTL3 promotes the expression of oncogenes including MYC and AF4 [25] . Furthermore, YTHDF1 has been linked with EIF3C that leads to the occurrence and metastasis of ovarian cancer [26] . Studies have also found that m6A RNA methylation modulators cause genetic alterations across diverse cancer types including mutations and copy number variations [27] . Moreover, m6A RNA methylation modulators have been shown to correlate with clinical parameters and have been utilised as tumour biomarkers to determine the diagnosis and prognosis of several cancer types and monitor the development [27] .
m6A and its modulators are widely studied in hepatocellular cancer [28] , colorectal cancer, and genital and nervous system tumours [26,[29][30][31] . However, research on m6A RNA methylation modulators in pancreatic cancer is limited. Recently, He et al. demonstrated that ALKBH5 inhibits pancreatic cancer motility by demethylating long non-coding RNA KCNK15-AS1 [32] . Other studies indicate that ALKBH5 decreases WIF-1 RNA methylation and inhibits pancreatic cancer tumorigenesis via the Wnt signalling pathway [33] . In this study, we obtained the expression pro les and clinical data from The Cancer Genome Atlas (TCGA) database of pancreatic cancer to explore the expression of 13 widely reported m6A RNA methylation modulators in PAAD and examine the in uence of m6A RNA methylation modulators in the clinical and pathological characteristics of PAAD.

Data collection
RNA-seq transcriptome data and corresponding clinical information of patients with PAAD were obtained from TCGA database (https://cancer genome.nih.gov/). We obtained the expression pro les of 178 patients with PAAD and 4 healthy samples. Clinical data of patients are presented in Table 1. Patients with incomplete data on age, gender, survival time (futime), whether they were censored (fustat), pathological grade of the tumour, and World Health Organization classi cation were excluded. In addition, samples with futime = 0 were also removed and we ultimately included 171 samples for analysis.

Bioinformatics analysis
ActivePerl 5.26 was used to extract the expression of the m6A methylation modulators and clinical data from TCGA database. R 3.6.3 (https://www.perl.org) was used for analyses of the data with the following packages: the limma package analysed the difference in expression between healthy and tumour samples, pheatmap and vioplot packages visually illustrated these differences, the corrplot package identi ed the correlations among the 13 modulators, ROC curve was drawn using the timeROC package, the forestplot package displayed the results of the univariate and multivariate cox analyses, glmnet and survival packages assisted with the LASSO regression, and the rms package portrayed a nomograph to predict the risk of death and the survival rate. STRING database (http://www.string-db.org/) was used to assist with the analysis of the interactions among the 13 m6A RNA methylation regulators.

Gene set enrichment analysis
PAAD samples were categorised into two groups according to the expression level of the genes of m6A RNA methylation modulators. KEGG pathways of m6A RNA methylation modulator genes were identi ed by gene set enrichment analysis (GSEA). Gene sets with nominal p-value < 0.05 and false discovery rate (FDR) q-value < 0.25 were considered signi cantly enriched.
2.5 Statistical analysis R 3.6.3 was used for our statistical analyses. The differential expression of m6A RNA methylation modulators between cancerous and healthy tissue was compared by the Wilcoxon test. Cox regression analysis was used to examine the prognosis of pancreatic cancer. Patients with PAAD in TCGA database were separated into a high-risk group and a lowrisk group based on the median risk score. Kaplan-Meier method was used to analyses the difference in overall survival (OS) between the two groups and Spearman's correlation was used to explore the correlations among the 13 m6A RNA methylation modulators.

Different expression of m6A RNA methylation modulators in PAAD and healthy tissues
The expression of m6A RNA methylation modulators in PAAD and healthy samples were analysed based on TCGA database. Transcriptome from TCGA showed that most of m6A RNA methylation modulators were not differentially expressed between the two groups except for WTAP (Figs. 1A and 1B), which was lowered in tumour samples (p = 0.012).
Further analyses were conducted to explore the correlations and interaction among the 13 modulators. We found the 13 m6A RNA methylation modulators were positively correlated with each other and that the coe cient of association between YTHDC1 and YTHDC2 was the highest. In addition, YTHDC1 expression was highly associated with METTL14 expression and was the only one that was positively correlated with all other modulators. Our analyses revealed that the proteins of m6A regulators showed strong levels of interaction (Figs. 1C and 1D).

Prognostic value of risk score based on the expression of m6A RNA methylation modulators in PAAD
To discover the in uence of m6A RNA methylation modulators in PAAD, we built a risk prediction model. First, we conducted a Cox univariate analysis, which showed that the expression of YTHDF1, ALKBH5, METTL3, KIAA1429, and METTL14 signi cantly correlated with survival ( Fig. 2A). High expression of KIAA1429 was related to lower survival whereas high expression of YTHDF1, ALKBH5, METTL3, and METTL14 decreased the risk of death and acted as protective factors. Variable selection was performed by LASSO Cox regression and the parameter λ indicated that the most suitable model to predict survival included YTHDF1, ALKBH5, METTL3, METTL14, and KIAA1429 with coe cients of -0.023, -0.021, -0.124, -0.178 and 0.233, respectively (Figs. 2B and 2C). The risk score was used to divide the patients into high-and low-risk groups and calculated using the following formula: risk score = Σ n i=1 Coef × x i where Coef represents the coe cient and x i is the expression value of each selected modulators. We found that the survival curve showed signi cant worse OS in high-risk tumour patients compared with patients in the low-risk group (Fig. 2D).
The expressions of YTHDF1, ALKBH5, METTL3, METTL14, and KIAA1429 were incorporated to form a nomogram to estimate 1-, 3-, and 5-year OS of pancreatic cancer (Fig. 3A). We used calibration curves to examine the precision of the nomogram (Figs. 3B -3D) and observed that all three calibration curves were close to the ideal curve, which indicated a good prediction value. The result of ROC curve predicted that PAAD patients with high-risk cancer had a 3-year OS rate of 21.8% compared with 44.9% in patients with low-risk cancer (Fig. 3E).

Risk score correlated with clinical characteristic of patient with PAAD
We investigated the relationship between the risk score from our prediction model and the clinical characteristic of patients with PAAD. Multiple linear regression showed that the risk score was signi cantly correlated with the grade and stage of the tumour ( Table 1), suggesting that a high-risk score corresponded to a high cancer grade and stage. We also analysed the relationship between the expression of m6A RNA methylation modulators and the clinical characteristics and multivariate regression indicated that the level of expression of all modulators except METTL3 was associated with clinical characteristics of patients with PAAD. The expression of YTHDF1 was inversely correlated with the tumour stage and the expression of ALKBH5 was negatively correlated with age, gender, and tumour stage of patients. Furthermore, METTL14 expression was negatively associated with tumour grade while KIAA1429 expression was correlated positively.
Importantly, high expression of ALKBH5 and METTL3 were related to increased survival whereas high expression of KIAA1429 was associated with a low survival.

Risk score was an independent prognostic indicator
We subsequently hypothesised that risk score could be an independent prognostic indicator and incorporated risk score and relevant clinical and pathological factors such as age, gender, and tumour grade and stage into univariate and multivariate Cox regression to test our hypothesis. Univariate Cox regression revealed that age, tumour grade and stage, and risk score were signi cantly linked with OS (Fig. 5A). Moreover, age and risk score were still signi cantly associated with OS in a multivariate Cox regression (Fig. 5B). These results indicated that age and risk score could be independent prognostic indicators of PAAD and we formed a nomogram with age, gender, tumour grade and stage, and risks core to predict 1-, 3-, and 5-year survival of patients with PAAD (Fig. 5C).

The associated KEGG pathways of m6A RNA methylation modulator genes
To identify the potential biological function of the genes that predicted the OS of patients with PAAD, we used GSEA to search for enriched KEGG pathways. The enrichment score (ES), FDR, and nominal p-value are shown for each gene in Table 2

Discussion
The high death rate of pancreatic cancer has been a continuous challenge in the eld. It can progress rapidly with only mild symptoms and there may be a multitude of causes for pancreatic cancer including a combination of lifestyle, environment, and genetic mutation. Surgery is the preferred treatment for resectable pancreatic cancer with a poor prognosis. Pancreatic cancer is generally insensitive to both chemotherapy and radiotherapy. Therefore, new treatment strategies focusing on the immunotherapy and targeted therapy are urgently needed for patients with pancreatic cancer. m6A is involved in the development of many cancers and some studies have shown its potential in the treatment of diverse cancers. In this study, we explored the function of m6A RNA methylation modulators in pancreatic cancer according to the transcriptome data from TCGA.
We analysed 13 widely reported m6A RNA methylation modulators in pancreatic cancer and our analyses revealed that WTAP was the only one with differential expression in PAAD patients. However, we did not nd a relationship between the expression of WTAP and the clinical prognosis. A high expression of WTAP has been reported in human pancreas cancer cell lines PANC-1, which enhances cell migration and invasion and suppressed its chemoresistance to gemcitabine in a Fak-dependent manner [34] . Further analysis and research is needed to identify the function of WTAP in PAAD.
We combined the expression of 13 m6A RNA methylation modulators with the clinical parameters of PAAD and Cox univariate analysis indicated that YTHDF1, ALKBH5, METTL3, KIAA1429, and METTL14 were associated with patient survival. Then we used LASSO Cox regression to create a survival prediction model including the ve m6A RNA methylation modulators. The resulting risk score had a good prediction of 1-, 3-, and 5-year OS of PAAD and a high-risk score increased the chance for a high grade and stage tumour. Therefore, even though the modulators were not differentially expressed between healthy controls and PAAD patients, their expression was related with OS. Additionally, we discovered that the expression of YTHDF1 showed an inverse relationship with the tumour stage. A previous report described that high expression of YTHDF1 is present in ovarian cancer and promotes cancer progression while knockdown of YTHDF1 inhibits tumorigenesis and peritoneal metastasis [26] . However, we did not observe similar ndings in our study. We also found a negative association between ALKBH5 and tumour stage. Previous studies have demonstrated that ALKBH5 suppresses pancreatic cancer motility by demethylating long non-coding RNA [32] and a recent study found that decreased ALKBH5 levels predict poor clinical outcome in pancreatic ductal adenocarcinoma, which is consistent with our data. Furthermore, the authors observed that ALKBH5 serves as a tumour suppressor by altering Wnt signalling via a decrease of Wnt inhibitory factor 1 RNA methylation, which in turn inhibits tumorigenesis. A high expression of ALKBH5 also increased sensitivity to chemotherapy [33] . Expression of METTL3 also affected the risk score and in uenced the OS of PAAD patients. However, we did not note any relationship with clinical characteristics.
Others have reported that the knockdown of METTL3 in pancreatic cancer cell lines decreases m6A RNA modi cation modulators and cripples cell proliferation, invasion, and migration. In addition, high METTL3 expression is associated with high pathological and N stage [35] . Additional studies are required to con rm the function of METTL3. Several publications describe that downregulation of METTL14 in pancreatic cancer tissue sensitises the tumour cells to cisplatin by enhancing apoptosis. Autophagy is also improved via an mTOR signalling-dependent pathway [36] . Following these results, we observed reduced expression of METTL14 enriched the mTOR signalling pathway in our GSEA analysis of METTL14. KIAA1429 has been reported to predict poor prognosis in hepatocellular carcinoma patients [37] . We found that KIAA1429 was also a promoter in PAAD and was involved in Wnt signalling, MAPK signalling, ERBB signalling, and pathways in cell cycle and cancer, indicating a central role in the development of PAAD.
A study published in 2019 revealed the regulation mechanism of pentose phosphate-related metabolism in tumour cells, which targets the adjustment of the tumour's cell cycle, affecting the energy requirement of tumour cells [38] . Other studies pointed out that ribosomes involved in epithelial-mesenchymal transition accelerate the metastasis of cancer cells [39] and that chemokines have a regulatory role in immune cell recruitment and tumour and stromal cell biology [40] . These pathways are important in the development of cancer and we also found these to be enriched in high-risk PAAD patients.
Taken together, this suggests that the m6A RNA modi cation modulators mediate these pathways to manage the survival of tumour cells. However, the exact mechanism of action remains to be elucidated with further research.

Conclusion
m6A RNA methylation modulators play an important role in the progression of cancer. In this study, we generated a prediction model using 13 m6A RNA methylation modulators in PAAD and identi ed their association with the prognosis.
Our results gain insight into the potential role of m6A RNA methylation in pancreatic cancer.