Gene Signature and Identification of Clinical Trait-Related m6 A Regulators in Pancreatic Cancer

Pancreatic cancer (PC) has a very poor prognosis and is usually diagnosed only at an advanced stage. The discovery of new biomarkers for PC will help in early diagnosis and a better prognosis for patients. Recently, N6-methyladenosine (m6A) RNA modifications and their regulators have been implicated in the development of many cancers. To investigate the functions and mechanisms of m6A modifications in the development of PC, 19 m6A regulators, including m6A-methyltransferases (ZC3H13, RBM15/15B, WTAP, KIAA1429, and METTL3/14), demethylases (FTO and ALKBH5), and binding proteins (YTHDF1/2/3, YTHDC1/2, IGF2BP1/2/3, HNRNPC, and HNRNPA2B1) were analyzed in 178 PC tissues from the cancer genome atlas (TCGA) database. The results were verified in PC cell lines Mia-PaCa-2, BXPC-3, and the control cell line HDE-CT. The m6A regulators-based sample clusters were significantly related to overall survival (OS). Further, lasso regression identified a six-m6A-regulator-signature prognostic model (KIAA1429, HNRNPC, METTL3, YTHDF1, IGF2BP2, and IGF2BP3). Model-based high-risk and low-risk groups were significantly correlated with OS and clinical traits (pathologic M, N, and clinical stages and vital status). The risk signature was verified as an independent prognostic marker for patients with PC. Finally, gene set enrichment analysis revealed m6A regulators (KIAA1429, HNRNPC, and IGF2BP2) were related to multiple biological behaviors in PC, including adipocytokine signaling, the well vs. poorly differentiated tumor pathway, tumor metastasis pathway, epithelial mesenchymal transition pathway, gemcitabine resistance pathway, and stemness pathway. In summary, the m6A regulatory factors which related to clinical characteristics can be involved in the malignant progression of PC, and the constructed risk markers may be a promising prognostic biomarker that can guide the individualized treatment of PC patients.


INTRODUCTION
Pancreatic cancer (PC) is one of the most lethal malignant neoplasms and has become one of the leading causes of cancerrelated deaths in developed countries (Ilic and Ilic, 2016). About 85% of PCs are adenocarcinoma, and less than 5% are pancreatic endocrine tumors (Wolfgang et al., 2013). There are usually no noticeable symptoms in the early stages of PC. When symptoms are specific enough to suggest PC, the disease might have reached an advanced stage. By the time of diagnosis, PC has often spread or metastasized to other parts of the body (Mohammed et al., 2014). With the development of medical techniques, PC can be diagnosed by ultrasound or computed tomography combined with blood tests and examination of tissue samples (biopsies). However, screening the general population for the early stage of the disease is not effective (Welinsky and Lucas, 2017). Pancreatic cancer can be treated with surgery, chemotherapy, targeted therapy, radiotherapy, palliative care, immunotherapy, or a combination of these based on the cancer stage (Mcguigan et al., 2018). With the current treatment methods pancreatic adenocarcinoma has a very poor prognosis: only 25% of patients with PC survive one year and 5 year-overall survival (OS) is lower than 5% (Ansari et al., 2015). Current treatment and diagnostic methods are not enough for the management of PC. Therefore, key goals for PC research are to develop novel prognostic markers, improve the early diagnostic rate, and find new targets for molecular targeted therapy. N6-methyladenosine (m 6 A), a potential biomarker, is a chemical modification present in multiple RNA species, which take part in various biological processes in cancer (Liu et al., 2018).
N6-methyladenosine regulators are involved in more than 60% of all RNA [messenger (mRNA), transport RNA (tRNA), and ribosomal RNA (rRNA)] modifications, which is an intense area of research for post-transcriptional regulation including translation, mRNA splicing, and mRNA stability (Yu et al., 2018). The level of modification of transcripts with m 6 A is regulated by methyltransferases, binding proteins and demethylases (Koh et al., 2019). The methyltransferases (including ZC3H13, RBM15, RBM15B, KIAA1429, METTL3/14, and WTAP), act as "writers, " and add the methyl group to the nitrogen on the sixth carbon of the aromatic ring of an adenosine residue (Meyer and Jaffrey, 2017). The cellular m 6 A status is reverted by demethylases (FTO, and ALKBH5; called "erasers"), and is recognized by m 6 A-binding proteins (HNRNPC, YTHDF1/2/3, YTHDC1/2, IGF2BP1/2/3, and HNRNPA2B1; called "readers") (Zaccara et al., 2019). N6-methyladenosine, a potential biomarker, is a chemical modification present in multiple RNA species, which take part in various biological processes in cancer (Liu et al., 2018). The dysregulation of m 6 A regulators is involved in the occurrence and development of multiple cancers, including bladder cancer, prostate cancer, head and neck squamous cell carcinoma, gastric cancer, breast cancer, hepatocellular carcinoma, and colorectal cancer (Hong, 2018). For example, METTL14 which suppresses colorectal cancer progression via regulating m 6 Adependent miR-375/yes-associated protein 1 (YAP1) pathway, is downregulated in colorectal cancer tissues and cell lines (Chen et al., 2020). FTO, a key m 6 A demethylase, is up-regulated in human breast cancer and is significantly associated with poor survival rates (Niu et al., 2019). FTO mediates m 6 A demethylation in the 3'UTR of BNIP3 mRNA and induces its degradation via an YTHDF2 independent mechanism, which indicates that FTO can serve as a novel potential therapeutic target for breast cancer (Niu et al., 2019). It has also been reported that IGF2BP2 regulates lncRNA DANCR through m6A modification, and IGF2BP2 and DANCR jointly promote the stemness-like characteristics of cancer and the pathogenesis of PC . Although more and more studies have shown that m6A regulatory factors play a crucial role in the pathogenesis and development of cancer, the fundamental relationship between m6A regulatory factors and PC remains unclear (Xia et al., 2019). The construction of prognostic signal based on m6A regulators that predicting the prognosis of PC will be helpful for prediction, prevention and personalized treatment.
This study used ConsensusClusterPlus to find that m 6 A regulators were closely related to PC OS rates in different clusters. Furthermore, lasso regression was used to identify a six-gene signature model (KIAA1429, HNRNPC, METTL3, YTHDF1, IGF2BP2, and IGF2BP3). Most of the genes identified were consistent with previous data (Taketo et al., 2018). For example, the m 6 A eraser ALKBH5, which was indicated as a potential therapeutic target for PC, was downregulated in PC cells and immortalized human pancreatic duct epithelial (HPDE6-C7) cells (He et al., 2018). Immunohistochemistry (IHC), western blots, and RT-qPCR were used to detect the expression of METTL3 in PC, and the results showed that METTL3 protein and mRNA levels were significantly higher in tumor samples than in paracancer samples. Down-regulation of METTL3 reduced the proliferation, invasion and migration of PC cell lines (Xia et al., 2019). While it is known that m 6 A plays important roles in different types of cancers, the available clinical traitrelated m 6 A regulator studies in PC are insufficient. Single-gene analysis are used to predict prognosis and to guide therapy in cancer. However, RNA-Seq is helpful for the construction of a prediction model using multiple genes. Here, we analyzed the gene signatures in different PC cell lines and identified clinical trait-related m 6 A regulators in PC. Additionally, potential related enrichment pathways of m 6 A regulators might be useful to further study their mechanisms of action.
We systematically searched for PC gene expression datasets that were publicly available and reported full clinical annotations. we download data from GSE28735 "Microarray gene-expression profiles of 45 matching pairs of pancreatic tumor and adjacent non-tumor tissues from 45 patients with pancreatic ductal adenocarcinoma" to validate the reliability of the built model. The raw data from the microarray datasets generated by Affymetrix and Illumina were downloaded from the Gene Expression Omnibus 2 .

Protein-Protein Interactions Network Construction and Correlation Analysis
The STRING database 3 was used for analyzing the proteinprotein interactions (PPI) among m 6 A regulators. The association among different m 6 A regulators was revealed by Spearman correlation coefficient with R package.

Cell Lines and Cell Culture
Two PC cell lines (Mia-PaCa-2 and BXPC-3) and one control cell line (HDE-CT) were purchased from China Center for Type Culture Collection (CCTCC, Shanghai, China). HDE-CT is a normal human pancreatic cell line and is cultured in DMEM medium (Corning, NY, United States) supplemented with 10% fetal bovine serum (GIBCO, South America, NY, United States). Mia-PaCa-2 with a KRAS mutation and BXPC-3 with wild type KRAS are human PC cell lines. Mia-PaCa-2 was cultured in DMEM medium with 10% fetal bovine serum, and BXPC-3 was cultured in RPMI-1640 medium with 10% fetal bovine serum. All cell lines were maintained in 5% CO 2 atmosphere at 37 • C.

RNA Extraction and qRT-PCR Verification
Total RNA of the four PC cell lines (Mia-PaCa-2 and BXPC-3) and HDE-CT were extracted with an RNA extraction kit (QIAGEN) according to the manufacturer's instructions. Briefly, 1 × 10 7 cells were collected and lysed for 10 min, genomic DNA was removed with an adsorption column, the samples were washed once with 75% ethyl alcohol and twice with wash buffer, and the samples were resuspended in RNA-grade enzymefree water. Total RNA was reversely transcribed into cDNA and used to perform quantitative real-time PCR (qRT-PCR) with SYBR Premix ExTaq (TaKaRa). GAPDH was used as a reference 2 https://www.ncbi.nlm.nih.gov/geo/ 3 http://string-db.org gene. Primers (Table 1) were synthesized by Sangon Biotech (Shanghai, China).

Primer name
Primer sequence (from 5 to 3 ) Frontiers in Genetics | www.frontiersin.org clustering methods use the proportion of ambiguous clustering (PAC) to infer optimal K (K-means) in order to identify and classify patients for further analysis (Lock and Dunson, 2013). Cluster analysis was performed using the ConsensusClusterPlus R package with cycle computation for 1000 times to ensure the stability and reliability of the classification (Wilkerson and Hayes, 2010). The Kaplan-Meier method was used for the OS analysis in different clusters.

Lasso Regression for PC Tissues
Lasso is a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the statistical model it produces. The best subset selection and the connections between lasso coefficient estimates can be identified to construct a prognostic model (Alhamzawi and Ali, 2018). Lasso regression was constructed to examine the relationship between gene signatures and PC risk. Further, clinical characteristics associated with OS were analyzed in patients with PC using Cox regression (including univariate and multivariate models) and the Kaplan-Meier method to evaluate the availability of the prognostic model. Pheatmap R package was used to correlate clinical data with the risk score (high or low).

Gene Set Enrichment Analysis for KIAA1429, HNRNPC, and IGF2BP2 in PC Tissues
Gene set enrichment analysis (GSEA) is widely used to analyze genome or proteome data, linking disease phenotypes with many different functional gene sets. The 178 patients with PC were divided into high expression groups and low expression groups according to the median expression values of KIAA1429, HNRNPC, and IGF2BP2. Two groups of TCGA data were analyzed by GSEA. Gene set enrichment analysis was also conducted in different sample risk groups based on the LASSO regression model. The 178 patients with PC were divided into high risk score group and low risk score group according to the median value of risk score.

Statistical Analysis
Gene expression data of FPKM form is used as input. WilcoxTest is used to get the p value for different expression between different clusters. The relationships between clusters or different risk score groups were analyzed using the Chi-square test. In all cases, p < 0.05 was considered statistically significant. Spearman correlation coefficient was calculated for the molecular pairing between m 6 A regulator genes. The student's t-test in SPSS 13.0 (SPSS Inc., Chicago, United States) was used to assess the expression differences between HDE-CT and PC cancer cells. Each experiment was repeated at least three times. Benjamini-Hochberg for multiple testing, and false discovery rate (FDR) were calculated to correct the p-value in GSEA.

Consensus Clustering for PC Tissues Based on the Expression of m6A Regulators
To determine whether the expression levels of m 6 A regulators were associated with PC prognosis, the TCGA PC cohort was clustered into different groups by consensus expression of m 6 A regulators with the ConsensusClusterPlus R package. Gene signatures of m 6 A regulators in PC are shown in Supplementary  Table S1. When the consensus matrix k value was equal to 2, there was no crossover between PC samples ( Figure 1A, Supplementary Figure S1 and Supplementary Table S2). The OS difference between different clusters was calculated by the Kaplan-Meier method and log-rank test ( Figure 1B and Supplementary Table S2). A heatmap was generated to visualize the expression pattern of m 6 A regulators between different clusters ( Figure 1C).

The Interaction and Correlation Among the m 6 A Regulators
The relationship between m 6 A regulators were further supported by the correlation analysis. Some highly correlated (|correlation coefficient| ≥ 0.5, p < 0.05) m 6 A regulator pairs were identified, including IGF2BP2 and IGF2BP3, IGF2BP2 and ALKBH5, YTHDC1 and YTHDC2, YTHDC1 and METTL14, YTHDC1 and ZC3H13, YTHDC2 and METTL14, YTHDC2 and ZC3H13, YTHDC2 and YTHDF3, METTL14 and FTO, METTL1 and ZC3H13, METTL14 and YTHDF3, FTO and ZC3H13 (Figure 2 and Supplementary Table S3). The interactions among the 19 m 6 A regulators are shown in Figure 3A. All m 6 A regulators have interactions in the same network. The results of the interaction network showed that IGF2BP1 and IGF2BP3, WTAP and KIAA1429, HNRNPC and HNRNPA2B1, WTAP and ZC3H13, METTL14 and METTL3, KIAA1429 and ZC3H13, METTL14

Gene Signature of m 6 A Regulators in PC Cell Lines
The expression of m 6 A regulators, including the m 6 A methyltransferases, the demethylases, and the m 6 A-binding proteins were analyzed by qRT-PCR in the PC cell lines, Mia-PaCa-2 and BXPC-3, and the control cell line HDE-CT. The results showed that some m 6 A regulators were differentially expressed in PC and control cell lines ( Figure 3B).

Lasso Regression Identified the Six-Gene Signature Prognostic Model
In order to determine the optimal prognostic model, lasso regression was performed using the glmnet R package. Lasso regression is a generalized linear model, and the adjustment degree of lasso regression complexity is controlled by lambda. The optimal six-gene signature prognostic model was identified when log (lambda) was between −2 and −3 (Supplementary  Figures S2A,B), where the coefficient of KIAA1429 was 0.28, the coefficient of HNRNPC was 0.34, the coefficient of METTL3 was −0.11, the coefficient of YTHDF1 was −0.37, the coefficient of IGF2BP2 was 0.28, and the coefficient of IGF2BP3 was 0.04. According to the median risk score, patients were divided into low-and high-risk groups (Supplementary Table S4). There was a significant difference in the OS rate between the two groups, and the OS rate of the high-risk group was significantly lower than that of the low-risk group (Figure 4A, p = 5.286e-04). A Receiver Operating Characteristiccurve (ROC) was used to evaluate the prediction efficiency of the prognostic signature. The prognostic signature model showed good prediction efficiency with the value of the area under the ROC curve (AUC) equal to 0.796 ( Figure 4B). Additionally, the KM plotter showed that the six selected m 6 A regulators were significantly related with OS according to OncoLnc 4 ( Figure 4C). Importantly, the heat map shows the expression of the six selected m 6 A regulators and clinicopathological variables in the high-and low-risk groups. Significant differences were found for the neoplasm histologic grade, pathologic M stage, pathologic N stage, pathologic stage, and vital status between high-and low-risk groups (Figure 5 and Supplementary Table S5).

The Effect of m 6 A Regulators on PC Prognosis
To investigate the effect of m 6 A regulators on PC prognosis, we performed Cox univariate ( Figure 6A) and multivariate 4 http://www.oncolnc.org/ analysis ( Figure 6B). The six-gene signature was consistent with the single-factor analysis of genes using Cox regression.  Figure 6A). The multivariate analysis revealed that age at initial pathologic diagnosis [hazard ratio (HR): 1.033; 95% confidence interval   Figure 6B). The factor of risk score based on the optimal six-gene signature prognostic model was significant both at univariate and multivariate analyses.  Table S6). Some interesting pathways are shown in Figure 7. It's not containing GSEA analysis for METTL3, IGF2BP1, and IGF2BP3. We did it, but there were no significant results for METTL3, IGF2BP1, and IGF2BP3. We conducted GSEA analysis in different sample risk score groups based on the LASSO regression model. The GSEA result showed that it is significantly related to cancer survival, oncogenesis by met, gemcitabine resistance, response to UV, HOXC6 targets cancer, recurrent liver cancer, WTAP targets, FIGURE 5 | The heatmap of sample risk groups and related pancreatic cancer clinical characteristics. Age at initial pathologic diagnosis (patients were aged 35-88), alcohol history documented (yes or no), alcoholic exposure category (daily drinker, weekly drinker, occasional drinker, social drinker, and non-drinker), anatomic neoplasm subdivision (body of pancreas, head of pancreas, tail of pancreas, and other parts), family history of cancer (yes or no), gender (male and female), history of chronic pancreatitis (yes or no), history of diabetes (yes or no), count of lymph nodes examined (from 1 to 57), neoplasm histologic grade (G1, G2, and G3), pathologic M (M represents tumor metastasis, including M0, M1,and MX), pathologic N (N represents tumor lymph node metastasis, including N0, N1, N2, and NX), pathologic T (T represents tumor size, including T1, T2, T3,T4, and TX), pathologic stage (Stages I, II, III, and IV), vital status (alive or dead). *p < 0.05 and **p < 0.01. tumor differentiated well vs. poorly, epithelial mesenchymal transition, hypoxia pathway, TGFB1 targets, cancer meta signature, and so on (Supplementary Table S7). Some interesting pathways are shown in Supplementary Figure S3, and those pathways closely related with tumorigenesis and development.

The Independent Verification by GEO
The different expression of m6A regulators between cancer tissue and normal tissue, including the m6A methyltransferases, the demethylases, and the m6A-binding proteins were analyzed based on the independent verification by GEO (Supplementary Figure S4 and Supplementary Table S8). In view of some similarities of identified different genes in TCGA data and GEO data, it is believed that the prognostic m6A regulators might not just be due to chance. For example, the overlapping genes that are significant were including RBM15B, KIAA1429, ALKBH5, YTHDF1, IGF2BP 2/3, and HNRNPC. Furthermore, the testing dataset based on GEO showed the different expression of m6A regulators in PC and validate the reliability of the built model based on TCGA. The optimal six-gene signature prognostic model was validated. According to the median risk score, patients from GEO were divided into low-and high-risk score groups (Supplementary Table S9). There was a significant difference in the OS rate between the two groups, and the OS rate of the highrisk score group was significantly lower than that of the low-risk score group (Supplementary Figure S5, p = 0.0012).

Experimental Validation
The inhibition of KIAA1429, HNRNPC, and IGF2BP2, respectively, significantly suppressed the proliferation abilities of PC cells based on CCK8 (Figure 8A). The EdU assay further showed that KIAA1429, HNRNPC, and IGF2BP2 inhibitors reduced DNA replication in both Mia-PaCa-2 and BXPC-3 cells (Figure 8B).

DISCUSSION
Treatment for PC has improved considerably, for example surgery with high success and lower complication rate is better than ever before, novel drug combinations (chemotherapy, target therapy, and immunotherapy) have been shown to improve survival rate, and advances in radiation therapy have achieved less toxicity; however, many researchers are focused on early diagnosis and prompt treatment as PC is still one of the deadliest solid malignancies (Chu et al., 2017). The development of multiomics has given us a better understanding of the fundamental genetics of PC. These advancements provide hope, but the survival rate of patients with PC is still poor (Cid-Arregui and Juarez, 2015). Biological functions of m 6 A were not studied extensively until around 2012, when major progress was made in the transcriptome profiling of m 6 A through antibody-based immunoprecipitation and high-throughput sequencing (Gan et al., 2019). Moreover, m 6 A regulators were shown to be related with the development of cancer (He et al., 2019). The process of m 6 A modification is reversible through the regulation of m 6 A methyltransferases, demethylases, and binding proteins. A series of m 6 A regulators have been described (Dominissini et al., 2012), including ZC3H13, RBM15/15B, KIAA1429, METTL14, YTHDC1/2, WTAP, METTL3, FTO, ALKBH5, YTHDF1/2/3, HNRNPA2B1, IGF2BP1/2/3, and HNRNPC (Lee et al., 2014). Therefore, it is necessary to explore the influence of m 6 A regulators on PC.
Recent studies have found that the m 6 A modification, when the related enzyme is abnormal, plays various roles in a series of human diseases such as neurological disorders, cancer, and embryonic developmental retardation . Both coding RNAs and some non-coding RNAs, such as lncRNA, microRNA, tRNA, and rRNA and RNA splice body, were regulated by an m 6 A modification before and after transcription (Yen et al., 2019). N6-methyladenosine modification is closely related to the metabolic processes of RNAs, for example, RNA processing, RNA transfer from the nucleus to the cytoplasm, RNA translation, RNA decay, and the biogenesis of RNA (Liang et al., 2020). The dynamic modification of RNA as a way of regulating genetic information is a new field of study, so there is still a lot of work to be done to understand the underlying mechanisms. Recently, a number of studies have found that m 6 A modifications are associated with cancer, having functions such as helping tumor stem cells to self-renew, promoting the growth and proliferation of cancer cells, and resisting radiotherapy or chemotherapy (Mao et al., 2019). All this evidence indicates that m 6 A regulators may be a target for cancer treatment (Bi et al., 2019;Ianniello et al., 2019). The regulation of m 6 A modifications is a collaboration between methyltransferases, demethylases, and binding proteins. The functions of these proteins in stem cell differentiation, stomach cancer, lung cancer, osteosarcoma, liver cancer, colorectal cancer, leukemia, neuroblastoma, renal cell carcinoma, and breast cancer have been extensivelyreported (Feng et al., 2019;Jin et al., 2019). For example, YTHDF1-deficient mice show an elevated antigen-specific CD8 + T cell antitumor response compared with wild-type mice, which indicated that durable neoantigen-specific immunity is regulated by mRNA m 6 A methylation through the m 6 A-binding protein YTHDF1 (Han et al., 2019). It was also reported that some drugs with antitumor activity, such as R-2-hydroxyglutarate (R-2HG), inhibited proliferation/survival of FTO-high cancer cells via targeting FTO/m 6 A/MYC/CEBPA signaling (Su et al., 2018). METTL3 which is independently of METTL14, binds to chromatin, and locates the transcription initiation site of active genes. The promoter bounding METTL3 induces m6A modification in the coding region of the relevant mRNA transcription and enhances its translation by alleviating ribosomal stalling. The gene regulated by METTL3 in this way is necessary for acute myeloid leukemia, suggesting that METTL3 may be a therapeutic target for acute myeloid leukemia (Barbieri et al., 2017). The researchers also found that that m 6 A mRNA demethylation by FTO increases melanoma growth and decreases response to anti-PD-1 blockade immunotherapy . Knockdown of FTO increased the methylation of m6A in the intrinsic genes of key primary melanoma cells such as PD-1 (PDCD1), CXCR4, SOX10, and so on, leading to increased attenuation of RNA in m6A reader YTHDF2, suggesting that FTO inhibition combined with anti-PD-1 blocking may abate the resistance of melanoma immunotherapy .
TCGA, a landmark cancer genomics project, described more than 20,000 primary cancers at the molecular level and matched normal samples of 33 cancer types. TCGA generated more than 2.5 petabytes of genome, epigenome, transcriptome and proteome data. The data has already lead to improvements in our ability to diagnose, treat, and prevent cancer (Blum et al., 2018). N6-methyladenosine RNA methylation regulators can lead to malignant progression and impact the prognosis of many kinds of cancer based on the TCGA database. For example, the lasso Cox regression model was applied to identify three m 6 A regulators in bladder cancer. The risk signature was constructed as follows: 0.164FTO -(0.081YTHDC1 + 0.032WTAP), which indicated that the three m 6 A regulators identified might be promising prognostic biomarkers to guide personalized treatment for patients with bladder cancer . Another study has built up a robust m 6 A regulators-based molecular signature that predicts the prognosis of patients with head and neck squamous cell carcinoma with high accuracy, which might provide important guidance for therapeutic strategies. The results revealed that the expression levels of YTHDF1, METTL3, KIAA1429, YTHDF2, RBM15, METTL14, ALKBH5, FTO, WTAP, and HNRNPC were significantly upregulated in head and neck squamous cell carcinoma samples, while YTHDC2 was remarkably downregulated (Zhao and Cui, 2019). In addition, a study identified two subgroups of gastric cancer (cluster1 and 2) by applying consistency clustering to the m6A regulators. Compared with the cluster1 subgroup, the prognosis of the cluster2 subgroup was poorer, and most of the 13 major m6A regulators were highly expressed in cluster2. This finding provides clues to understand epigenetic modifications of RNA in gastric cancer (Su et al., 2019). However, the prognostic role of m 6 A regulators in PC is poorly understood. In the present study, we are the first to show, by applying consensus clustering to m 6 A regulators, that there are two subgroups of PC (cluster1 and 2). The cluster2 subgroup correlates with a poorer prognosis, which suggests that m 6 A regulators may be promising prognostic biomarkers for PC. Furthermore, the lasso regression analysis identified a six-gene signature prognostic model (KIAA1429, HNRNPC, METTL3, YTHDF1, IGF2BP2, and IGF2BP3). These results agree with the results of previous studies. The major function of IGF2BP2 is to regulate cell metabolism . However, our results suggest that lncRNA DANCR is a novel target for IGF2BP2 through m 6 A modification in PC, and that it promotes cancer stemness-like properties and PC pathogenesis. Mechanistically, IGF2BP2 serves as a reader for the m 6 A modified DANCR (at adenosine 664), and the definite interaction site provides a novel target for PC therapy .
We did GSEA for KIAA1429, HNRNPC, and IGF2BP2. Many enrichment pathways were significantly related to cancer pathogenesis. We focused on some important events, for example, pathways of oncogenesis by Met, EIF4 pathway, downregulation of SMAD2-SMAD3-SMAD4 transcriptional activity, EZH2 targets, stemness, well vs. poorly differentiated tumor, epithelial mesenchymal transition, UV response via ERCC3, and metastasis. The identified pathways were consistent with reported data. The importance of m 6 A in the response to ultraviolet DNA damage was demonstrated, and the findings support that m 6 A RNA serves as a beacon for the selective, rapid recruitment of DNA polymerase κ to damage sites to facilitate repair and cell survival (Xiang et al., 2017). Meanwhile, many studies show that m 6 A-related genes work on stemness regulation in tumor relapse. For example, METTL3 was identified as a regulator for terminating murine naïve pluripotency. METTL3 knockout preimplantation epiblasts lead to early embryonic lethality, because it is associated with stability of key naïve pluripotency-promoting transcripts (Geula et al., 2015). Epithelial mesenchymal transition (EMT), as an important cellular program during tumor migration, invasion and metastasis, is also regulated by m 6 A mRNA methylation. N6-methyladenosine-sequencing and functional studies confirm that YTHDF1 mediates m 6 A-increased translation of Snail mRNA (a key transcription factor of EMT) . Interestingly, the process of m 6 A mRNA methylation was also regulated by cytokines (Li et al., 2017). The TGFβ pathway plays roles in disease through the intracellular effectors SMAD2 and SMAD3. SMAD2/3 promotes binding of the m 6 A methyltransferase complex to a subset of transcripts involved in early cell fate decisions. These aspects of m 6 A methyltransferase signaling could have far-reaching implications in the treatment of many cancers (Bertero et al., 2018).
In conclusion, this study is the first to identify and profile the gene signatures of clinical trait-related m 6 A regulatory genes in PC. We also developed a six-gene signature prognostic model, which might play a crucial role in determining the clinical progression of PC. With the development of m 6 Asequencing and methylated RNA immunoprecipitation, m 6 A regulatory genes might serve as promising molecular biomarkers for monitoring many kinds of cancers and providing important guidance for selecting therapeutic strategies.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in The Cancer Genome Atlas (TCGA), https://cancergenome.nih.gov/.

AUTHOR CONTRIBUTIONS
JH and LL designed the experiments. ZW and HL performed the experiments, data collection, and analysis. HZ wrote the manuscript. LL revised the manuscript and provided financial support. All authors approved the final version of the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2020.00522/full#supplementary-material   FIGURE S4 | Expression differences of m 6 A RNA methylation regulators in pancreatic cancer based on GEO data. Red and green represent relatively high or low expression, respectively. * p < 0.05, * * p < 0.01, and * * * p < 0.001.