The m6A-Related mRNA Signature Predicts the Prognosis of Pancreatic Cancer Patients.

N6-methyladenosine (m6A) has an important epitranscriptomic modification that controls cancer self-renewal and cell fate. The addition of m6A to mRNA is a reversible modification. The deposition of m6A is encoded by a methyltransferase complex involving three homologous factors, jargonized as "writers," "erasers," and "readers." However, their roles in pancreatic adenocarcinoma (PAAD) are underexploited. With the use of The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases, we provided an mRNA signature that may improve the prognostic prediction of PAAD patients based on the genetic status of m6A regulators. PAAD patients with genetic alteration of m6A regulators had worse disease-free and overall survival. After comparing PAAD groups with/without genetic alteration of m6A regulators, we identified 196 differentially expressed genes (DEGs). Then, we generated a 16-mRNA signature score system through least absolute shrinkage and selection operator (LASSO) Cox regression analysis. Multivariate cox regression analysis demonstrated that a high-risk score significantly correlates with poor prognosis. Moreover, time-dependent receiver operating characteristic (ROC) curves revealed it was effective in predicting the overall survival in both training and validation sets. PAH, ZPLD1, PPFIA3, and TNNT1 from our signature also exhibited an independent prognostic value. Collectively, these findings can improve the understanding of m6A modifications in PAAD and potentially guide therapies in PAAD patients.

Introduction

Pancreatic adenocarcinoma (PAAD) remains a worldwide lethal disease. Despite recent advances, the 5-year survival rate remains low. Although surgery offers the best long-term survival, most PAAD patients miss the chance of tumor resection due to atypical symptoms at an early stage. Moreover, current adjuvant therapies (such as chemotherapy and radiotherapy), which are typically guided by the tumor, node, metastasis (TNM) staging system cannot effectively improve patient prognosis, suggesting biological heterogeneities exist among PAAD patients, and the stratification by the TNM staging system alone may be inadequate. Recently, an increasing number of studies have begun to subgroup PAAD patients in different ways, such as stratifying PAAD patients by their genetic/epigenetic features.2, 3, 4, 5

Modifications of RNAs, especially mRNA, are vital in the post-transcriptional regulation of gene expression. The N6-methyladenosine (m6A) modification regulates different stages of mRNA metabolism, including maturation, folding, translation, export, and decay, and thus, consequently, drives numerous biological processes. m6A in mRNA has emerged as an important epitranscriptomic modification that controls cancer self-renewal and cell fate. As one of the most abundant post-transcriptional modifications present in mammalian mRNA, several studies suggest that changes in m6A modification patterns are implicated in tumorigenesis, leading to various cancers, such as breast, lung, glioblastoma, and many more. m6A methylations are usually enriched around stop codons in the 5′ and 3′ untranslated regions and within long internal exons expressing the consensus sequence RRACH, where R = purine, A = m6A, and H = A, C, or U. Such modifications can help modulate biological processes, including RNA splicing, mRNA translocation, degradation and translation, cell proliferation, differentiation, and survival. In general, m6A modifications are manipulated by m6A regulators. There are three types of m6A regulators: writers, readers, and erasers.

Writers consist of RNA methyltransferases, which install the modifications. Examples of writers include the complex METTL3/METTL14 and METTL16 and their cofactors. Erasers (fat mass and obesity-associated protein [FTO] and alkB homolog 5 [ALKBH5]) are RNA demethylases reversing m6A methylation to balance mRNA modification. Readers (YTHDC1–2, YTHDF1–3, and IGF2BP1–3) bind to m6A-modified mRNA and exert biological functions, such as mRNA translocation, degradation, and translation., Whereas a previous study found that m6A writer METTL3 could promote chemo- and radioresistance in PAAD cells, one recent study discovered that several hypomethylated genes correlated with poor overall survival of PAAD patients. Erasers (FTO and ALKBH5) are RNA demethylases reversing m6A methylation to balance mRNA modification. Readers are proteins exerting regulatory effects on mRNA metabolism by selectively binding to m6A and exerting biological functions, such as mRNA translocation, degradation, and translation., There are several families of m6A-binding proteins. One such family is YTHDC (i.e., YTHDC1–2, YTHDF1–3). YTHDC1 proteins are found in the nucleus, directing mRNA splicing, whereas YTHDC2 and YTHDF proteins are predominantly cytoplasmic, mediating translational efficiency and decay of m6A-modified mRNAs. Other readers include IGF2BP1/2/3, all found in the nucleus. Although m6A regulators have been well studied in other diseases,11, 12, 13 like leukemia,, their roles in PAAD are still underexploited.

With the creation of next-generation sequencing technologies, it is now feasible to obtain a clearer picture of the mutational and transcriptional landscape of most tumors. This has been elucidated through many large-scale studies, such as The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) databases. These analyses will identify many of the core genetic pathways activated in PAAD and help enable identification of distinct molecular subtypes associated with differences in therapy response. Our current study aims to assess the landscape of m6A regulators in a PAAD cohort from TCGA database. We compared the clinical and molecular features between groups with distinct m6A regulator alteration status and produced and validated a prognostic mRNA gene signature by least absolute shrinkage and selection operator (LASSO) Cox regression analysis. Ultimately, we hope to provide a stable prognostic tool, as well as potential therapeutic targets, for PAAD treatment.

Results

Genetic Alterations of m6A Regulators in TCGA PAAD Patients

From 149 TCGA PAAD cases analyzed, the m6A writer gene VIRMA (6.7%) demonstrated the highest frequency of genetic alteration (mutations and/or copy number variations [CNVs]), followed by HAKAI (4%), YTHDF3 (3.4%), IGF2BP1 (3.4%), and ALKBH5 (2.7%) (Figure 1A; Table 1). Among CNVs of m6A regulatory genes, amplifications (number of events = 33) were the most dominant event compared to deep deletions (number of events = 4) (Figure 1B; Table 1). 75.17% of PAAD cases did not have genetic alteration of the m6A regulator (Figure 1C). Subsequently, although no difference regarding disease-free survival (DFS) (p = 0.300) was found in the PAAD training set, PAAD cases with the alteration of m6A regulators exhibited significantly worse overall survival (OS) (p = 0.0341) compared with those without the alteration of m6A regulators (Figures 1D and 1E).

Figure 1

Genetic Alteration of the m6A Regulator in PAAD

(A) Percentage of PAAD samples with genetic alteration (mutation and/or CNV) from TCGA database. (B) Events of copy number amplification (gain) or deep deletion (loss) in PAAD samples. (C) Genetic alteration patterns of m6A regulators in PAAD samples from TCGA database. (D and E) Disease-free (D) and overall survival (E) analysis of PAAD patients with/without genetic alteration of m6A regulators, realized through online website cBioPortal. p < 0.05 was considered as significant difference.

Table 1

Different Genetic Alteration Patterns of m6A Regulators in PAAD Samples (n = 149)

		No Alterations	Genetic Alteration				Altered/Profiled
			Missense Mutation	Truncating Mutation	Amplification	Deep Deletion
Writers	METTL3	146	2	0	1	0	2.0%
	METTL14	148	0	1	0	0	0.7%
	METTL16	148	0	0	0	1	0.7%
	WTAP	146	2	1	0	0	2.0%
	RBM15	148	0	1	0	0	0.7%
	RBM15B	148	1	0	0	0	0.7%
	HAKAI	143	2	0	4	0	4.0%
	VIRMA	139	0	0	10	0	6.7%
	ZC3H13	147	2	0	0	0	1.3%
Erasers	FTO	147	1	0	1	0	1.3%
	ALKBH5	145	1	0	3	0	2.7%
Readers	YTHDC1	147	1	0	0	1	1.3%
	YTHDC2	148	1	0	0	0	0.7%
	YTHDF1	147	1	0	1	0	1.3%
	YTHDF2	148	0	0	0	1	0.7%
	YTHDF3	144	0	0	5	0	3.4%
	IGF2BP1	144	1	0	3	1	3.4%
	IGF2BP2	146	1	0	2	0	2.0%
	IGF2BP3	146	0	0	3	0	2.0%

Clinicopathological Features of PAAD Groups with Different m6A Regulator Status

Next, we explored the relationship between genetic alterations of m6A regulators and clinicopathological features of PAAD patients. The results demonstrated that PAAD patients with/without genetic alterations of m6A regulators had no significantly distinct features, including age, gender, pathological stage, histological grade, or TNM stage (Table 2).

Table 2

Clinicopathological Features of PAAD Patients with or without Mutation/CNV of m6A Regulators

		Without Mutation and/or CNVa	With Mutation and/or CNVa	p
Age	<65	39	28	0.524
	≥65	44	39	–
Gender	female	38	29	0.76
	male	45	38	–
Pathological stage	I	5	11	0.153
	II	74	51	–
	III	2	1	–
	IV	1	2	–
	N/Ab	1	2	–
Histological grade	G1	12	13	0.204
	G2	51	29	–
	G3	19	23	–
	G4	1	1	–
	Gxb	0	1	–
T stage	T1	2	3	0.787
	T2	10	10	–
	T3	69	51	–
	T4	2	1	–
	Txb	0	1	–
	N/Ab	0	1	–
N stage	N0	18	21	0.125
	N1	64	42	–
	Nxb	1	3	–
	N/Ab	0	1	–
M stage	M0	39	25	0.341
	M1	1	2	–
	Mxb	43	40	–

N/A, not applicable.

With mutation and/or CNV: TCGA PAAD patients with mutant or CNV or mutant + CNV; without mutant and/or CNV: TCGA PAAD patients with neither mutant nor CNV.

Ambiguous variables were excluded from chi-square test or Fisher exact test.

Gene Ontology (GO) Enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis of Differentially Expressed Genes (DEGs) and Immune Cell Fraction Estimation

To identify distinct biological activities between PAAD groups with a different m6A regulator status, we detected a total of 196 DEGs (Figure 2A; Table S1) and subsequently performed GO enrichment and KEGG pathway analysis (Figures 2B–2D). Notably, the top 3-enriched terms were “NABA matrisome associated,” “pancreatic secretion,” and “digestion.”

Figure 2

The m6A-related Differences

(A) Heatmap showed 196 DEGs between the group with genetic alteration of the m6A regulator and the group without genetic alteration of the m6A regulator in TCGA database. (B–D) Gene Ontology (GO) enrichment and KEGG pathway analysis of the DEGs, realized through online website Metascape. Relevant annotations were on the right side. (E) Violin plot for comparison of the immune cell fraction difference between the group with genetic alteration of the m6A regulator (red) and the group without genetic alteration of the m6A regulator (green). p < 0.05 was considered as significant difference.

To explore one potential mechanism, we estimated the immune cell fraction in PAAD groups with a different m6A regulator status, summarized in Figure 2E. In both groups, uncommitted macrophages (M0), alternatively activated macrophages (M2), and cluster of differentiation (CD)4+ memory resting T cells accounted for the three largest fractions from all of the immune cells. Interestingly, the fraction of monocytes was significantly higher in the PAAD group without alteration of the m6A regulator. Other types of immune cells did not exhibit a significant intergroup difference.

Development of Prognostic mRNA Signature from the Training Set

The aforementioned PAAD groups with different m6A regulator status were treated as the training set to produce a prognostic mRNA signature on OS. As demonstrated in Table 2, there existed no significant difference in age, gender, pathological stage, histological grade, or TNM stage. We also conducted multivariate analysis using the Cox proportional hazards model to determine prognostic factors. After univariate COX regression analysis, 21 genes were selected with p < 0.1 from the DEGs for further analysis (Table S2). Then, LASSO coefficient profiles of DEGs were analyzed, as previously described, and a coefficient profile plot was produced after the log2 transformation of the lambda (λ) value (Figure 3A). A vertical line was drawn at the value selected by 10-fold cross-validation (Figure 3B). The resistance coefficient method (k method) resulted in 16 optimal coefficients. Among the 16 optimal coefficients, six mRNAs were shown to be upregulated, whereas ten mRNAs were downregulated in the PAAD group with m6A regulator alteration compared to the PAAD group without m6A regulator alteration. With the use of the LASSO Cox regression model, we generated a 16-mRNA signature to calculate the risk score for every PAAD patient, based on the expression level of these 16 mRNAs weighted by their regression coefficients: risk score = (−0.07744 × expression of CDHR3) + (−0.09939 × expression of CELSR3) + (0.29421 × expression of epidermal growth factor [EGF]) + (0.29019 × expression of fibroblast growth factor 10 [FGF10]) + (−0.60946 × expression of GAD1) + (0.08812 × expression of MT1H) + (0.09277 × expression of NMUR2) + (−0.48445 × expression of PAH) + (0.05361 × expression of PGC) + (−0.23567 × expression of PGM5) + (−0.33405 × expression of POPDC2) + (0.06505 × expression of PPFIA3) + (0.05930 × expression of SERPINA4) + (−0.37688 × expression of TMEM145) + (0.20989 × expression of TNNT1) + (0.38926 × expression of ZPLD1).

Figure 3

The mRNA Signature Associated with m6A Regulator Status

(A and B) LASSO coefficient profiles of the reaming 21 DEGs after univariate Cox regression analysis. (A) The profile of the LASSO coefficient. (B) The partial likelihood deviance is shown against log(Lambda). At the value fitting 10-fold cross-validation, a vertical line was drawn. (C and D) Time-dependent ROC curves at 1, 2, and 3 years (C) and Kaplan-Meier survival analysis between TCGA PAAD patients with low and high risk based on our mRNA signature (D). (E and F) Time-dependent ROC curves at 1, 2, and 3 years (E) and Kaplan-Meier survival analysis between ICGC PAAD patients with low and high risk based on our mRNA signature (F).p < 0.05 was considered as significant difference.

Association between Risk Score and Survival of PAAD Patients

Next, to explore the prognostic value of the mRNA signatures in PAAD patients, we analyzed the effects of risk score on the OS and DFS. As expected, both OS and DFS in univariate/multivariate regression models in the training set were correlated with risk score: for OS, univariate: hazard ratio (HR): 1.741 (95% confidence interval [CI]: 1.476–2.053, p = 0.000); multivariate: HR: 1.659 (95% CI: 1.230–2.237, p = 0.001), whereas for DFS, univariate: HR: 2.087 (95% CI: 1.631–2.670, p = 0.000); multivariate: HR: 1.944 (95% CI: 1.290–2.930, p = 0.001) (Table 3). Furthermore, we conducted a time-dependent receiver operating characteristic (ROC) curve analysis at 1, 2, and 3 years in the training set to assess the prognostic accuracy of the prognostic mRNA signature. We observed the area under the curve (AUC) of 1-, 2-, and 3-year survival to be 0.752, 0.586, and 0.628, respectively (Figure 3C).

Table 3

Cox Regression Analysis of m6A Regulatory Genes for Overall Survival (OS) and Disease-Free Survival (DFS) of PAAD Patients

Variablea	OS					DFS
	Univariate		Multivariate		Univariate			Multivariate
	HR (95% CI)	p	HR (95% CI)	p	HR (95% CI)		p	HR (95% CI)	p
Age	1.421 (0.906–2.229)	0.126	–	–	1.624 (0.897–2.941)		0.109	–	–
Gender	0.770 (0.498–1.193)	0.242	–	–	0.685 (0.386–1.217)		0.197	–	–
Stage	1.131 (0.355–3.597)	0.835	–	–	1.541 (0.212–11.236)		0.670	–	–
T	2.330 (1.159–4.683)	0.018b	1.078 (0.363–3.205)	0.892	1.964 (0.868–4.442)		0.105	1.795 (0.526–6.135)	0.350
N	2.097 (1.192–3.688)	0.010b	1.246 (0.525–2.958)	0.618	2.154 (1.063–4.362)		0.033b	1.268 (0.470–3.416)	0.639
M	1.239 (0.294–5.211)	0.770	1.135 (0.227–5.682)	0.877	1.059 (0.142–7.904)		0.955	1.143 (0.151–8.670)	0.897
Grade	1.455 (0.919–2.302)	0.110	–	–	1.505 (0.821–2.760)		0.186	–	–
Risk score	1.741 (1.476–2.053)	0.000b	1.659 (1.230–2.237)	0.001b	2.087 (1.631–2.670)		0.000b	1.944 (1.290–2.930)	0.001b

Ambiguous variables were excluded (–).

Variables were grouped and compared as follows: age (≥65 versus <65); gender (male versus female); stage (III–IV versus I–II); T (T3–T4 versus T1–T2); N (N1 versus N0); M (M1 versus M0); grade (G3–G4 versus G1–G2).

P < 0.05.

Subsequently, PAAD patients in the training set were divided into a low-risk group and a high-risk group by the median risk score as the cutoff value. As expected, the high-risk group of PAAD patients had worse OS (p < 0.001; Figure 3D). In the validation set (ICGC cohort), we report similar findings: the AUC of time-dependent ROC analysis was 0.535, 0.636, and 0.602 in 1-, 2-, and 3-year survival, respectively (Figure 3E). Furthermore, the high-risk group had a poorer prognosis (p = 0.043) (Figure 3F).

Taken together, our results from both the training and validation sets suggest that the alteration of m6A regulators may predict poor survival. Lastly, we also explored the effect of each gene from our mRNA signature on OS and found that PAH (HR = 0.66; log-rank p = 0.048), ZPLD1 (HR = 1.8; log-rank p = 0.008), PPFIA3 (HR = 0.58; log-rank p = 0.0095), and TNNT1 (HR = 1.7; log-rank p = 0.0084) may act as independent OS indicators (Figures 4A–4D). Moreover, unlike ZPLD1 and PPFIA3, PAH and TNNT1 were differentially expressed between PAAD cancer tissues (red box) from TCGA and normal pancreatic tissue (gray box) from TCGA and Genotype-Tissue Expression (GTEx) (Figures 4E–4H).

Figure 4

Prognostic Genes from the mRNA Signature

(A–D) The expression level of PAH (A), ZPLD1 (B), PPFIA3 (C), and TNNT1 (D) successfully predicted the overall survival of PAAD patients. (E–H) Different expression levels between PAAD tissue (red box) and normal tissue (gray box) groups for PAH (E), ZPLD1 (F), PPFIA3 (G), and TNNT1 (H). p < 0.05 was considered as significant difference.

Discussion

The first aim of our study is to provide an mRNA signature to improve the prognostic accuracy of PAAD patients. Consortium efforts, such as those of TCGA and the ICGC, have been a tremendous asset for this purpose. Whereas pancreatic cancer survival rates have been improving from decade to decade, improving survival outcomes of PAAD patients still challenges medical decisions, including surgical resection and/or adjuvant therapies. Currently, surgery and/or adjuvant therapies can bring about unavoidable risks, such as potential relapse. Furthermore, because the definitions of borderline resectable and locally advanced pancreatic cancer vary among institutions and countries, it is impossible to compare survival rates according to clinical stage in pancreatic cancer patients. Moreover, the data published by most institutions do not include patients with metastatic or locally advanced pancreatic cancer. Consequently, an accurate prognosis for PAAD patients is essential for proper individualized therapy. Genes involved in PAAD development allow for advances in risk stratification, which potentially can outperform the current pathological staging system.17, 18, 19 Moreover, the identification of specific genetic changes in PAAD can result in a better understanding of the molecular mechanisms in the development of PAAD and help determine effective therapeutic strategies. Thus, in this study, we depicted the landscape of PAAD, stratified by genetic alteration status, and derived an mRNA prognostic signature.

m6A modification indicates new directions for the treatment of various cancers. Regulators or inhibitors of m6A modifications may provide the potential therapeutic strategies for cancers.

m6A is the most common and abundant methylation modification in mRNA. The methylation modification of m6A has been proven to be reversible through the regulation of methyltransferase (writers), demethylase (erasers), and proteins that recognize m6A modification (readers). “Writers” catalyze the formation of m6A; “erasers,” which include FTO and ALKBH5, selectively remove the methyl code from target mRNAs; and “readers” are capable of decoding m6A methylation and generating a functional signal, including YTH domain-containing protein, eukaryotic initiation factor (eIF) 3, IGF2BP families, and heterogeneous nuclear ribonucleoprotein (HNRNP) protein families. The YTH domain can recognize m6A through a conserved aromatic cage, and another two proteins—FMR1 and LRPPRC—“read” this modification. Contrary to the conventional writer–eraser–reader’ paradigm, few studies reveal METTL3/16 as an m6A writer or reader. m6A RNA modification is a dynamic and reversible process that was corroborated by the discovery of eraser.

Previous studies have indicated that genetic alterations, like the mutations in the m6A regulator genes, may cause various functional alterations and thus, influence physiological and/or pathological processes, including axon guidance, viral replication, and tumor progression.20, 21, 22, 23, 24 In this study, based on distinct genetic alteration status (including mutations and/or CNVs) of m6A regulators (i.e., m6A writers, erasers, and readers), we divided a PAAD cohort from TCGA database into two groups: a group with m6A regulator alteration and a group without m6A regulator alteration. After stratification, survival analysis—specifically, overall survival—demonstrated intergroup heterogeneity (Figures 1D and 1E), which suggests different biological activities between two PAAD groups.

To our knowledge, studies on m6A modification in PAAD are limited, yet in other cancer types, studies have well illustrated the role of certain m6A regulators, including METTL3, VIRMA, ALKBH5, FTO, and YTHDF2, being associated with the tumor proliferation, differentiation, tumorigenesis, proliferation, invasion, and metastasis and functioning as oncogenes or anti-oncogenes in malignant tumors. With the m6A writer METTL3 as an example, in leukemia, promoter-bound METTL3 induces m6A modification within the coding region of the related mRNA transcript and enhances the translation. Furthermore, downregulation of METTL3 could lead to cell-cycle arrest, differentiation of leukemic cells, and a failure to establish leukemia mouse models. In lung cancer, METTL3 acts as an oncogene in lung cancer by increasing EGF receptor (EGFR) and Tafazzin (TAZ) expression and promoting cell growth, survival, and invasion. METTL3-eIF3-caused mRNA circularization promotes the translation and oncogenesis of lung adenocarcinoma. Besides, small ubiquitin-like modifier (SUMO)ylation of METTL3 is of importance for the promotion of tumor growth at lysine residues K177, K211, K212, and K215 in non-small cell lung carcinoma (NSCLC). In liver cancer, METTL3 is frequently upregulated and capable of reflecting malignancy, as well as acting as an independent poor prognostic factor., Likewise, METTL3 could promote the progression of bladder cancer, but METTL14 is an anti-metastatic factor and serves as a favorable factor in hepatocellular carcinoma (HCC) by regulating m6A-dependent microRNA (miRNA) processing. In conclusion, METTL3 upregulation or METTL14 downregulation predicts poor prognosis in patients with HCC and contributes to HCC progression and metastasis. One way to explain the pathological role of METTL3 refers to epithelial-mesenchymal transition (EMT), an important step for cancer cell metastasis. Deletion of METTL3 could impair the EMT of cancer cells both in vitro and in vivo. However, m6A modification is multifaceted in the cancer environment. Reduced expression of METTL3 may be responsible for reductions in m6A methylation in 70% of endometrial tumors, and such changes could lead to increased proliferation and tumorigenicity of endometrial cancer cells through the AKT pathway, in which expression of a negative AKT regulator, like PHLPP2 is decreased, and expression of a positive AKT regulator, like mammalian target of rapamycin complex 2 (mTORC2), is increased. Similar to the multifaceted role of METTL3, the expression of m6A writer VIRMA, the gene with greatest amplification in our study, also discriminated between seminomas and nonseminomatous tumors, according to a recent study.

As for m6A erasers, the silencing of ALKBH5 inhibits cancer growth and invasion by disturbing EMT and angiogenesis-related transcripts, including transforming growth factor (TGF)-β signaling pathway genes. Another m6A eraser, FTO, could enhance leukemic oncogene-mediated cell transformation and leukemogenesis and inhibit all-trans retinoic acid-induced acute myeloid leukemia (AML) cell differentiation through regulating expression of targets, such as Ankyrin Repeat And SOCS Box Containing 2 (ASB2) and Retinoic Acid Receptor Alpha (RARA), by reducing m6A levels in these mRNA transcripts. Moreover, FTO could also serve as a novel potential therapeutic target for breast cancer. Among several m6A readers, YTHDF2 promotes targeted mRNA decay. Mapping of m6A in RNAs from mouse hematopoietic stem and progenitor cells and human umbilical cord hematopoietic stem cells demonstrated m6A enrichment in mRNAs encoding transcription factors essential for stem-cell self-renewal. In both YTHDF2 knockout hematopoietic stem and progenitor cells and YTHDF2 knockdown human umbilical cord hematopoietic stem cells, these mRNAs were stable and facilitated hematopoietic stem-cell expansion while knocking down one of the key targets of YTHDF2: Tal1 mRNA, which partially rescued the phenotype. Collectively, these studies corroborate the functional importance of m6A modifications, such as METTL3, METTL14, FTO, and YTHDF2, and they provide profound insights into development and maintenance of AML and self-renewal of leukemia stem/initiation cells through the downstream MYC and Tal1 pathways.

The above information demonstrates that m6A modification can target multiple genes participating in various critical biological processes, such as transcription, cell proliferation, and cancer-related pathways. However, to our surprise, unlike previously reported studies,33, 34, 35 this study, with the focus on PAAD, did not reveal a very meaningful association between m6A status and immune-infiltration patterns, even though a significant fraction change of monocytes could be observed. Thus, the understanding of the role of m6A modification in PAAD is essential for understanding the potential consequences of therapeutic intervention, as well as making an accurate prognosis. In similar studies on glioma and clear cell renal cell carcinoma, m6A regulator alterations were correlated with a poorer prognosis, and an indication that novel therapeutic strategies for m6A RNA methylation should be further explored in the treatment of cancer.

Here, the candidates for an mRNA signature were selected based on m6A regulators. After the division of PAAD patients into two groups based on m6A regulator alteration status, we generated an mRNA signature from the DEGs. Specifically, our stable prognostic 16-mRNA signature was established by LASSO Cox regression analysis and contains both up- and downregulated genes. Among upregulated genes, according to previous reports, Cadherin EGF LAG Seven-Pass G-Type Receptor 3 (CELSR3) is epigenetically dysregulated in 84% of small intestinal neuroendocrine tumors (SINETs). SINET is the most common malignancy of the small intestine. Epigenetic changes in GAD1 expression facilitates cancer metastasis by altering glutamate metabolism in the microenvironment of metastatic brain cancer. Conversely, among the downregulated genes in our mRNA signature, FGF10 signaling was previously reported to play a role in the development of aggressiveness in pancreatic cancer cells through interactions with FGFR2 (FGF receptor 2); the abolishment of PGC-1α was associated with the resistance of pancreatic cancer stem cells against the anti-diabetic drug metformin; EGF, as a part of the tumorigenic EGF-EGFR system in PAAD, could induce translocation of RhoA from the cytosol to the membrane fraction and cause cell rounding in human pancreatic cancer cells, and such a process could be reversed by high-mobility group-coenzyme A (HMG-CoA) reductase inhibitor fluvastatin. Noticeably, the EGFR inhibitor erlotinib offers a potential therapeutic tool in PAAD treatment. The combination of the EGFR inhibitor erlotinib with gemcitabine outperforms gemcitabine alone in PAAD treatment. One mechanism demonstrated that gemcitabine induces mitogen-activated protein kinase (MAPK) signaling, which could be dramatically inhibited by erlotinib. In addition, EGF-TGF-beta interactions could increase pancreatic cell invasion, which could be blocked by erlotinib and SB505124, a type I TGF-beta receptor inhibitor.

Aside from the above-studied genes in our mRNA signature, we also identified some under-exploited genes with independent prognostic value in PAAD, such as TNNT1, PPFIA3, PAH, and ZPLD1 (Figure 4). According to studies, TNNT1 is overexpressed in breast cancer and leiomyosarcoma;, PPFIA3 is methylated in most gastric cancer samples, whereas barely methylated in normal samples,, suggesting m6A writers could outperform erasers during such carcinogenesis. Likewise, in our 16-mRNA signature, the number of downregulated genes exceeds the number of upregulated genes. Such a phenomenon may be attributed to the dominance of m6A writers over m6A erasers and consequential aberrant downstream activity (such as EMT) caused by m6A methylated mRNA.

Overall, our study has contributed to the field in the following aspects. First, we established and validated an mRNA signature from DEGs between groups with distinct m6A regulator status; this can potentially help prognosticate PAAD. Second, genes from our mRNA signature could participate in the development of PAAD and therefore, serve as potential therapeutic targets. Third, rather than make comparisons between normal and cancer groups and use a conventional pathological staging system to sub group patients, we have developed a novel method to stratify PAAD patients by their m6A regulator alteration status, offering a new perspective to identify heterogeneities among PAAD patients.

Limitations to the study include the lack of biological verification. Future molecular studies on interactions between m6A regulators and members from our mRNA signature can benefit the understanding of PAAD. Furthermore, future studies with prospective validation are still warranted to support our findings.

In conclusion, we report, for the first time, an mRNA signature to prognosticate and potentially guide therapies in PAAD patients, as well as a novel m6A regulator-based method for PAAD patient stratification. Novel therapeutic strategies for m6A RNA methylation should be further explored in the treatment of PAAD.

Materials and Methods

Datasets and Data Preprocessing

All clinical and sequencing data are available in public databases. It is acknowledged that the necessary consent has been achieved. PAAD transcript profiles were obtained from the Genomic Data Commons (GDC) Data Portal (TCGA GDC: https://portal.gdc.cancer.gov/), and somatic mutation profiles of TCGA PAAD were masked as the training group. For the validation group, ICGC Australia Pancreatic Cancer (PACA) data were obtained from the University of California, Santa Cruz (UCSC), Xena (UCSC Xena: https://xenabrowser.net/) and normalized by log2 transformation. cBioPortal (cBioPortal: https://www.cbioportal.org/) performed the survival analysis and provided the mutation information.

Clinical information (age, gender, TNM stage, tumor histological grade, survival time, and survival status) was also obtained from the aforementioned databases. All probes were mapped to the latest gene symbols, according to Ensembl identification (ID). The median value was adopted as expression level if a particular gene had multiple probes.

Identification of m6A-Related DEGs

Alterations (mutation and/or CNV) of m6A regulatory genes (writers: METTL3, METTL14, METTL16, WTAP, RBM15, RBM15B, HAKAI, KIAA1429, and ZC3H13; erasers: FTO and ALKBH5; and readers: YTHDC1–2, YTHDF1–3, and IGF2BP1–3) were screened in each TCGA PAAD sample. Subsequently, TCGA cohorts were divided into 2 groups: a group with m6A alterations and a group without m6A alterations. Linear models for microarray (LIMMA) data were adopted to generate the DEGs between 2 groups. The threshold for DEGs was p <0.05 and fold change ≥2, and R package “pheatmap” was used to draw the heatmap.

GO Enrichment, KEGG Pathway Analysis, and Immune Cell-Type Fractions Estimation of DEGs

GO enrichment and KEGG pathway analysis of the DEGs were realized through Metascape (Metascape: http://metascape.org). To explore potential the role of DEGs, the CIBERSORT method was adopted to reflect proportions of immune cells in PAAD samples. This algorithm was conducted, as described in previous research. Briefly, a total of 547 gene-expression values were set as references and considered to represent the minimum for each cell type. According to these values, immune cell-type proportions could be inferred from the transcriptome data of tumor samples through support vector regression. For cases with p <0.05, fractions of the immune cell population could be considered as accurate CIBERSORT results. Such a method was used to determine differences in the fractions of various immune cell types between the group with m6A alterations and the group without m6A alterations. R package “vioplot” was used to visualize corresponding results.

Construction of Prognostic Model (mRNA Signature) and Further Analysis

With the use of univariate COX regression, we identified a list of mRNAs with p < 0.01 from DEGs. The LASSO Cox regression model was used to construct a prognostic model from these mRNAs. OS was the endpoint. R package “glmnet” offered a sequence of λs and various prognostic models. Ten-fold cross-validation minimum criteria were used to select the λs minimum (min) with the minimum mean across validation error. Every included case would obtain a risk score. R package “survival” calculated survival rates and intergroup difference of survival curves and produced Kaplan-Meier (K-M) curves with the log-rank test. A standard formula was generated to calculate the risk score for each patient. In this study, the LASSO Cox method was utilized to reduce the dimensionality and to secure the most significantly overall survival-associated DEGs to build a prognostic model using the Cox regression method. The risk score for each patient was calculated by a standard formula, which combines the expression levels of the mRNAs and LASSO Cox regression coefficients (λs). Risk score is equal to the result of (λ1 × expression of A) + (λ2 × expression of B) + (λ3 × expression of C) + … + (λn × expression of N), in which “λ” represents the regression coefficient of each gene.

To determine the predictive accuracy of the risk score, R package “timeROC” generated time-dependent ROC curves with corresponding AUC. Univariable and multivariable Cox regression models were used to test the independent prognostic ability of risk score. Furthermore, ICGC data were used to validate the risk score model through ROC and K-M curves. Finally, we evaluated the prognostic ability of each gene from our mRNA signature on overall survival using the Gene Expression Profiling Interactive Analysis (GEPIA) database (GEPIA: http://gepia.cancer-pku.cn/).

Statistical Analysis

Except for the aforementioned analysis by R software (version 3.25.0; R: https://www.r-project.org), all other statistical data and figures were analyzed by SPSS 20.0 (IBM, Chicago, IL, USA) and GraphPad Prism 6.0 (GraphPad Software, La Jolla, CA, USA). The associations between m6A mutation/CNV and different clinicopathological features were analyzed with chi-square test or Fisher exact test. p <0.05 was considered to be statistically significant.

Author Contributions

H.W., X.J., and Z.M. designed the study. Q.Y., Z.M., and J.Z. collected and analyzed data. Q.Y., B.W., and S.L. wrote the manuscript. Z.M., Q.Y., R.O., and X.J. participated in the revision of the manuscript. All authors read and approved the manuscript.

Conflicts of Interest

The authors declare no competing interests.