Skip to main content
Download PDF

Quantification of m6A RNA methylation modulators pattern was a potential biomarker for prognosis and associated with tumor immune microenvironment of pancreatic adenocarcinoma

BMC Cancer volume 21, Article number: 876 (2021) Cite this article

Abstract

Background

m6A is the most prevalent and abundant form of mRNA modifications and is closely related to tumor proliferation, differentiation, and tumorigenesis. In this study, we try to conduct an effective prediction model to investigated the function of m6A RNA methylation modulators in pancreatic adenocarcinoma and estimated the potential association between m6A RNA methylation modulators and tumor microenvironment infiltration for optimization of treatment.

Methods

Expression of 28 m6A RNA methylation modulators and clinical data of patients with pancreatic adenocarcinoma and normal samples were obtained from TCGA and GTEx database. Differences in the expression of 28 m6A RNA methylation modulators between tumour (n = 40) and healthy (n = 167) samples were compared by Wilcoxon test. LASSO Cox regression was used to select m6A RNA methylation modulators to analyze the relationship between expression and clinical characteristics by univariate and multivariate regression. A risk score prognosis model was conducted based on the expression of select m6A RNA methylation modulators. Bioinformatics analysis was used to explore the association between the m6Ascore and the composition of infiltrating immune cells between high and low m6Ascore group by CIBERSORT algorithm. Evaluation of m6Ascore for immunotherapy was analyzed via the IPS and three immunotherapy cohort. Besides, the biological signaling pathways of the m6A RNA methylation modulators were examined by gene set enrichment analysis (GSEA).

Results

Expression of 28 m6A RNA methylation modulators were upregulated in patients with PAAD except for MTEEL3. An m6Ascore prognosis model was established, including KIAA1429, IGF2BP2, IGF2BP3, METTL3, EIF3H and LRPPRC was used to predict the prognosis of patients with PAAD, the high risk score was an independent prognostic indicator for pancreatic adenocarcinoma, and a high risk score presented a lower overall survival. In addition, m6Ascore was related with the immune cell infiltration of PAAD. Patients with a high m6Ascore had lower infiltration of Tregs and CD8+T cells but a higher resting CD4+ T infiltration. Patients with a low m6Ascore displayed a low abundance of PD-1, CTLA-4 and TIGIT, however, the IPS showed no difference between the two groups. The m6Ascore applied in three immunotherapy cohort (GSE78220, TCGA-SKCM, and IMvigor210) did not exhibit a good prediction for estimating the patients’ response to immunotherapy, so it may need more researches to figure out whether the m6A modulator prognosis model would benefit the prediction of pancreatic patients’ response to immunotherapy.

Conclusion

Modulators involved in m6A RNA methylation were associated with the development of pancreatic cancer. An m6Ascore based on the expression of IGF2BP2, IGF2BP3, KIAA1429, METTL3, EIF3H and LRPPRC is proposed as an indicator of TME status and is instrumental in predicting the prognosis of pancreatic cancer patients.

Peer Review reports

Background

Pancreatic cancer is an aggressive cancer of the digestive system with high morbidity and mortality. It is the seventh leading cause of cancer death in both males and females worldwide because of its poor prognosis. The high death rate of pancreatic cancer has been a continuous challenge to the medical field [1]. 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 [2]. Studies have predicted that it could be the second leading cause of cancer death by 2030. The estimated incidence of pancreatic cancer in the United States in 2020 is around 57,600 with a projected 47,050 deaths, indicating 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 85% of pancreatic cancer [3]. While the precise cause and underlying mechanisms of PAAD are still unclear, genetics is an essential factor related to tumour development and mutations in DNA or RNA may also drive the initiation of PAAD [4].

RNA modification is a crucial part of epigenetics, which together with gene and protein modification, plays an important role in regulating cellular processes. According to the data from the MODOMICS update in 2017, about 163 different modifications have been identified in RNA molecules so far, including N1-methyladenosine, N7-methyladenosine, 5-methylcytosine, pseudouridine, N6,2′-O-dimethyladenosine (m6A) and 2′-O-methylation. In eukaryotes, m6A is the most prevalent and abundant form of mRNA modification, which was identified in the 1970s [5, 6]. m6A is a dynamic and reversible process of mRNA modification and has been widely studied, although the exact mechanism of m6A modification 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 reflection of m6A modification. These proteins are classified into three categories: writers, erasers, and readers, which function as modulators in the process of m6A modification. Writer proteins mainly include METTL3, METTL14, WTAP, CBLL1 and METTL16 [7,8,9,10]. RBM15, ZC3H13, and KIAA1429 have been newly identified 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 members that recognise m6A modification, including the family of YT521-B homology (YTH) domain-containing proteins (YTHDF1, YTHDF2, YTHDF3, YTHDC1, and YTHDC2) and the IGF-2 mRNA-binding proteins (IGF2BP1, IGF2BP2, and IGF2BP3) [17, 18], heterogeneous nuclear ribonucleoproteins hnRNPA2B1 and hnRNPC also have been identified as m6A RNA readers [19].

m6A plays a vital role in the splicing, translation, and stability of RNA and regulates chromatin state and transcription [20]. This has significant consequences in various human diseases where m6A alterations may enhance heart failure [21] influence 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, and play a dual role in cancer development. On the one hand, m6A regulates the expression of oncogenes or tumor suppressor genes, thereby affecting cancer progression. On the other hand, it can regulate the expression and activity of the m6A enzyme, thereby affecting the role of m6A in cancer. In bladder cancer, m6A modulator METTL3 promotes oncogenes’ expression [25]. Furthermore, YTHDF1 has been linked with EIF3C, leading to the occurrence and metastasis of ovarian cancer [26]. Studies have also found that m6A RNA methylation modulators cause genetic alterations such as mutations and copy number variations across diverse cancer types [27]. Moreover, m6A RNA methylation modulators have been shown to correlate with clinical parameters and have been utilized as tumour biomarkers to determine the diagnosis and prognosis of several cancer types and monitor tumor 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]. 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 through the Wnt signaling pathway [33]. However, therapies targeting signaling pathways in pancreatic cancer lack efficacy due to the heterogeneity of tumor microenvironment (TME) in patients, TME prevents antitumor drug resistance, immune escape and tumor cell’s infiltration during cancer progression. Increasing evidence demonstrated that m6A modification is associated with TME of cancer, and the polarization of tumor-associated macrophages-M2 enabled the oxaliplatin resistance via the elevation of METTL3-mediated m6A modification [34]. Dali et al. found a risen level of CD8+ cytotoxic T cells and natural killer (NK) cells while a reduced infiltration of myeloid-derived suppressor cells (MDSC) from Ythdf1−/− tumor mice compared to WT mice, besides, the therapeutic efficacy of PD-L1 checkpoint blockade is enhanced in Ythdf1−/− mice [35].

Since the tumor is a result of the coordinated action of multiple factors and genes, we integrated the gene expression profiles information from The Cancer Genome Atlas (TCGA) and GTEx database, extracted expression of 28 widely reported m6A RNA methylation modulators from 140 patients with PAAD, assessed the influence of m6A RNA methylation modulators in the clinical and pathological characteristics of PAAD comprehensively, explored its relationship with the TME cell-infiltrating characteristics and predicted the efficacy of immune checkpoint inhibitors.

Materials and methods

Data collection and pre-processing

RNA-seq transcriptome data (FPKM value) and corresponding clinical information of pancreatic cancer patients were obtained from TCGA database (https://cancergenome.nih.gov/). The gene expression profile was measured experimentally using the Illumina HiSeq2000 RNA Sequencing platform by the University of North Carolina TCGA genome characterization center. We obtained expression profiles of 178 patients with PAAD; however, 34 cases should not have been included in the PAAD study according to the quality annotation from TCGA, and one patient was a repeated measure, furthermore, PAAD patients with incomplete data on age, gender, survival time (futime), survival status (fustat), pathological grade of the tumor, and pathologic stage were excluded. A total of 140 patients with PAAD were eventually enrolled in our study from TCGA database in the end. The excluded and enrolled patients were listed in Table S1. RNA-seq transcriptome data of 167 normal human pancreatic tissues was download from the Genotype-Tissue Expression (GTEx) project (https://xenabrowser.net/). As to GTEx, FPKM values were extracted, and log2(x + 0.001) transformed into FPKM values, then the FPKM values from TCGA and GTEx database were transformed into transcripts per kilobase million (TPM) values. Data of gene expression quantification was integrated from the TCGA database and GTEx project; 140 patients with PAAD and 167 healthy populations were incorporated into our study in the end.

We obtained the list of immune-related genes from import database (https://immport.niaid.nih.gov/), a long-term, sustainable data warehouse developed to promote re-use of immunological data generated by NIAID DAIT and DMID funded investigators. Besides, the copy number of the genes in 140 PAAD patients was download in the UCSC database (https://xenabrowser.net/). Expression of prognosis m6A modulators mRNA levels was analyzed to compare normal and tumor tissues in various cancer types using the Oncomine database (https://www.oncomine.org). The GSE62452 dataset from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) was used for validation.

Extraction of m6A RNA methylation modulator expression data

We selected m6A RNA methylation modulators which were widely reported in a variety of cancers. Twenty-eight m6A RNA methylation modulators including 10 writers, 16 readers, and 2 erasers were included in this study (Table 1). mRNA expression data of 28 m6A RNA methylation modulators from TCGA and GTEx projects were extracted for further bioinformatics analysis.

Table 1 the lists of the 28 m6A RNA methylation regulative factors
Full size table

Constructing of m6A RNA methylation modulators signature in pancreatic cancer

Twenty-eight m6A RNA methylation modulators were incorporated into univariate Cox proportional hazard regression, then modulators with a p-value < 0.1 were involved in the subsequent LASSO regression for further screening, and the selected modulators were applied to conducted prognosis model. Here m6Ascore = Σni = 1Coef × xi where Coef represents the coefficient and xi is the expression value of each selected modulator, and nomograph was portrayed to predict the risk of death and the survival rate.

Pathways and function enrichment analysis

PAAD samples were categorised into high-risk and low-risk groups according to the median of m6Ascore, “limma” package of R was used to analyse the differentially expressed genes (DEGs) with stated threshold values where false discovery rate (FDR) < 0.05 and |Log2 fold change (FC)| > =1.0, all the DEGs was performed in volcan plot by “ggplot2” package in R, intersection of DEGs and immune-related genes from import database was displayed in venn plot by “VennDiagram” package of R software and exhibited in heatmap partly. Then the function of differentially expressed immune-related gens between two group was performed by Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, visualised by “ggplot2” package in R. KEGG pathways of m6A RNA methylation modulator genes were identified by GESA-3.0.jar software. Gene sets (c2.cp.kegg. v7.4 symbols.gmt) with nominal p-value < 0.05 and false discovery rate (FDR) q-value < 0.25 were considered significantly enriched.

Estimation of tumor-infiltrating immune cell

Estimation proportions of the 22 tumor-infiltrating immune cells (B cells naive, B cells memory, Plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), T cells gamma delta, NK cells resting, NK cells activated, Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Mast cells activated, Eosinophils, and Neutrophils) from each sample were calculated by using the “CIBERSORT” (R package) according to the gene expression profiles. The tumor cellularity was inferred through the ESTIMATE algorithm using the “estimate” package in R with default parameters. We acquired ImmuneScore, StromalScore, and ESTIMATEScore of each sample in pancreatic cancer to predict the level of infiltrating stromal and immune cells based on specific gene expression signatures of stromal and immune cells.

Immune response analysis

The immunophenoscore (IPS) of 140 PAAD patients was obtained from The Cancer Immunome Database (TCIA) (https://www.tcia.at/home), which provides results of comprehensive immunogenomic analyses of next-generation sequencing data (NGS) data for 20 solid cancers from The Cancer Genome Atlas (TCGA) and other data sources. Furthermore, three immunotherapeutic cohorts, including IMvigor210 cohort, TCGA-SKCM cohort, and GSE78220 cohort downloaded from GEO were used to predict the m6Ascore for curative effect of immunotherapy, expression data and detailed clinical information of IMvigor210 cohort were obtained from http://research-pub.gene.com/imvigor210corebiologies, and normalized by the DEseq2R package in R and then the count value was transformed into the TPM value. As to GSE78220 and TCGA-SKCM cohort, the FPKM data of gene expression profiles were also converted to the more comparable TPM value.

Clinical tissue specimens

A total of six paired pancreatic cancer and adjacent non-tumorous tissue used for immunohistochemistry (IHC) staining were collected from the First Affiliated Hospital of Anhui Medical University. All patients have already provided written informed consent. This work was approved by the Academic Committee of The First Affiliated Hospital of Anhui Medical University and was conducted following the Declaration of Helsinki.

Immunohistochemistry staining analysis

Collected tissue specimens were formalin-fixed and embedded with paraffin. Tissue sections (4 μm thickness) were used in IHC staining analyses, and the slices were treated with methanol containing 3% hydrogen peroxide to inactivate the endogenous peroxidase and treated with citric acid buffer (pH = 6.0) to obtain optimal antigen recovery. The slices were incubated in 1% bovine serum albumin in phosphate buffer for 30 min to block non-specific binding. In addition, the slices were stained with primary antibody overnight at 4 °C. Then, these sections were subjected to three 5-min mild washing in phosphate buffer saline, followed by staining with secondary antibody (HRP polymer) at 1:200 for 50 min. Diaminobenzidine was applied to the slices before being counterstained with hematoxylin. Finally, the samples were sealed, observed, and photographed by a light microscope.

Statistical analysis

R 4.0.3 (https://www.r-project.org/) was used for our statistical analyses. Wilcoxon test was used to perform a difference comparison between two group, including differential expression of m6A RNA methylation modulators between cancerous and healthy tissues and the DEGs between high and low m6Ascore group. The relationship between m6Ascore and age, gender, pathological stage and histologic grade were analysed by Chi-square test. Association between m6Ascore and tumor-infiltrating immune cells was computed by Spearman’s correlation. Multiple regression analysis was used to analyse the association between m6A RNA methylation modulators and clinical features. Cox regression analysis was used to examine the prognosis of pancreatic cancer. Kaplan-Meier method was used to analyses the difference in overall survival (OS) between the high-risk and low-risk group. P-values less than 0.05 on both sides were statistically significant. Graphics were visualized by R and Graphad prism 8.

Results

Expression of m6A RNA methylation modulators was upregulated in patients with PAAD and associated with cancer progression

Twenty-eight m6A RNA methylation modulators (10 writers, 2 erases and 16 readers) were identified in our study. Expression of 28 m6A RNA methylation modulators was extracted from 167 healthy samples from GTEx database and 140 patients with PAAD integrated from TCGA database (Fig. 1A), All the m6A RNA methylation modulators showed a high expression in pancreatic carcinoma than normal pancreatic tissue except METTL3. The landscape expression of all the 28 m6A RNA methylation modulators is displayed in Fig. 1B. We further analysed the change in copy number of these modulators (Fig. 1C and D), and the alteration frequency of copy number was observed to be common in these modulators, copy number of PRRC2A, KIAA1429, EIF3H, and IGF2BP2 showed a high increment frequency, in contrast, ELAVL1, YTHDF2, and METTL3 displayed a loss of copy number, indicating an abnormal expression of m6A RNA methylation modulators. A positive correlation among m6A RNA methylation modulators was shown in Fig. 1E and F, revealing m6A RNA methylation modulators interacted with each other to influence pancreatic cancer progression. The network results displayed 12 modulators function as risk factors in PAAD patients, including LRPPRC, IGF2BP2, IGF2BP3, HNRNPC, FMR1, EIF3H, EIF3A, FTO, CBLL1, RBM15, YTHDF3, and KIAA1429.

Fig. 1
figure 1

Expression of m6A RNA methylation modulators in 140 patients with PAAD and 167 normal samples based on TCGA and GTEx database. The expression levels of 28 m6A RNA methylation modulators in patients with PAAD presented in A a heatmap (green means high expression while blue means low expression and blue means normal samples while pink means tumour samples) and B a boxplot (blue means normal samples while red means tumour samples). The association between the expression of the modulators and survival time with statistically significant C The CNV variation frequency of 28 m6A regulators in 140 patients with PAAD. D The location of CNV alteration of m6A regulators on 23 chromosomes. E Spearman’s correlation analysis of the 28 m6A RNA methylation modulators in PAAD patients. F Interaction associations among 28 m6A RNA methylation modulators are visualized in the network and the relationship between m6A regulators and prognosis. *p < 0.05, **p < 0.01, ***p < 0.001

Full size image

m6A RNA methylation modulators were associated with clinical characteristic of PAAD

We built a prognosis prediction model to discover the influence of the 28 m6A RNA methylation modulators in PAAD. First, we conducted a Cox univariate analysis, which showed that the expression of KIAA1429, IGF2BP2, IGF2BP3, METTL3, EIF3H and LRPPRC was correlated with the survival (Fig. 2A). Modulators less than 0.1 in cox univariate analysis, including IGF2BP2, IGF2BP3, LRPPRC, KIAA1429, EIF3H, and METTL3, were subsequently filtrated. Variable selection was performed by LASSO Cox regression and the parameter λ indicated that the most suitable model to predict survival included KIAA1429, IGF2BP2, IGF2BP3, METTL3, EIF3H and LRPPRC with coefficients of − 0.052, 0.007, 0.005, 0.02, 0.004, and − 0.08, respectively (Fig. 2B-C). The expressions of the selected modulators were incorporated to form a nomogram to estimate 1-, 2-, and 3-year OS of pancreatic cancer (Fig. 2D). We used calibration curves to examine the precision of the nomogram (Fig. 2E-G) and observed that all three calibration curves were close to the ideal curve, which indicated an excellent prediction value. To analyze the association between expression of the six selected modulators and clinical parameters, we combined the expression data and clinical parameters of 140 PAAD patients, and the baseline characteristics of these patients are listed in Table 2. The survival curve of the six selected modulators showed high expression of LRPPRC, IGF2BP2, IGF2BP3, EIF3H, and KIAA1429 was related to decreased survival while overexpression of METTL3 displayed a favourable prognosis (Figure S1). We analyzed the relationship between the expression of the six modulators and the clinical characteristic of PAAD, and found a high expression of EIF3H, IG2BP3, IG2BP2, and LRPPRC was related to the maximum tumor dimension, while reduced expression of METTL3 may be associated with the maximum tumor dimension, overexpression of IGF2BP2 and LRPPRC was relevant to a higher histologic grade (Table 3).

Fig. 2
figure 2

Prognostic value of risk score from m6A RNA methylation modulators expression in PAAD patients. A Cox univariate analysis of m6A RNA methylation modulators. B and C Selection of genes by LASSO Cox regression. D Nomogram predicting 1-, 2-, and 3-year survival of patients with PAAD based on the expression of selected m6A RNA methylation modulators. The top row shows the point value for each variable, rows 2–6 indicate the expression of selected m6A RNA methylation modulators, the sum of these values is located on the Total Points axis, and the line drawn downward to the survival axes is used to determine the likelihood of 1-, 2- or 3-year survival. Calibration of nomograms predicting. Calibration of nomograms predicting E 1-year, F 2-year, and G 3-year survival

Full size image
Table 2 Demographic and clinicopathological parameters of PAAD patients based on TCGA datasets
Full size table
Table 3 Association between the expression of eight m6A RNA methylation modulators and clinical characteristic in patients with PAAD based on the TCGA datasets
Full size table

m6Ascore correlated with clinical prognostic of patient with PAAD and was an independent prognostic indicator

The expression of the six m6A RNA methylation modulators was merged using the m6Ascore formula and defined as a risk score, patients were then divided into high- and low-risk groups (Fig. 3A). The result of survival curve predicted that PAAD patients with high m6Ascore had a 1,2 and 3-year OS rate of 61.7, 23.3 and 17.5%, respectively, in patients with low-risk, 1,2 and 3-year OS rate was 82.7, 44.0 and 32% respectively. The risk score model showed a good prediction of 1,2 and 3-year OS, survival curve showed worse OS in high-risk patients than patients in the low-risk group (Fig. 3B and C). Then we investigated the relationship between the risk score and the clinical characteristic of patients with PAAD. Risk score was associated with living status and maximum tumor dimension of patients with PAAD (Fig. 3D, Table 2).

Fig. 3
figure 3

Relationship between m6A risk score and OS of PAAD patients. A landscape of risk grouped by risk score. B Survival analyses for low and high m6Ascore patient groups using Kaplan-Meier curves. C Time-dependent ROC analysis of risk score in predicting prognoses. D The difference of clinical features between high and low risk group

Full size image

We subsequently hypothesised that risk score could be an independent prognostic indicator and then incorporated risk score and relevant clinical and pathological factors including age, gender, tumour grade and stage into univariate and multivariate Cox regression to test our hypothesis. Univariate Cox regression revealed that risk score was significantly linked with OS (Fig. 4A). Moreover, risk score was still associated with OS in a multivariate Cox regression with HR = 1.702(95%CI:1.382–2.096, P < 0.001) (Fig. 4B). These results indicated that m6Ascore could be an independent prognostic indicator of PAAD, the death rate of patients with high-risk score was nearly two times higher than those in the low-risk group, we formed a nomogram with age, gender, tumour grade, stage, and risk core to predict 1-, 2-, and 3-year survival of patients with PAAD (Fig. 4C), calibration curves showed a good prediction value of our risk model (Fig. 4D-F).

Fig. 4
figure 4

The risk score is an independent prognostic indicator of PAAD patients. A Univariate and B multivariate Cox regression analysis of the association between clinical and pathological factors, including the risk score and OS of PAAD patients. C Nomogram based on age, gender, tumour grade and stage, and risk score. Calibration of nomograms predicting D 1-year, E 2-year, and F 3-year survival

Full size image

The function analysis and the pathways of the genes that enriched in the group of PAAD patients with high-risk and low-risk

Further analysis was conducted to determine the different expressed genes between high and low m6Ascore patients, all the genes between the two groups were presented in volcano plot (Fig. 5A), around 840 genes showed differential expression from TCGA, 1793 immune related-genes were acquired from immport database, the venn plot showed 94 immune related-genes among the different genes, including 6 upregulated and 88 downregulated genes in patients with high m6Ascore (Fig. 5B). the top 50 different immune related-genes were displayed in heatmap according to the following filter: false discovery rate (FDR) < 0.05 and |Log2 fold change (FC)| > =1.0 (Fig. 5C), IL36B, IL36RN and CCL26 were the top 3 downregulated gene in low m6Ascore group with more remarkable Log2 fold change, whereas CD19, AGTR2 and TNFRSF13C were the top 3 overexpression in patients with high m6Ascore (Fig. 5C). The GO functional enrichment analysis was carried out to exhibit the function of 94 immune related-genes between high m6Ascore and low m6Ascore group, a total of 344 GO terms were enriched, including 295 biological process terms, 5 cellular component terms and 44 molecular function terms. Kegg pathway analysis was performed to present the related biological pathways among these DEGs, a total of 19 Kegg pathways were enriched among these DEGs, the top 10 GO terms of three functional groups and all enriched Kegg pathways were listed in our results with the most significant P values and bigger counts (Fig. 5D-E). DEGs take part in the modulation of T cell costimulation, activating cell surface receptor signaling, activating signal transduction, and join the plasma membrane signaling receptor complex, T cell receptor complex and external side of the plasma membrane influenced receptor-ligand, signaling receptor activator, and cytokine activity, on the other hand. Analysis of Kegg showed an enrichment in cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer progression.

Fig. 5
figure 5

Identifying differentially expressed immune-related genes and enrichment analysis. A Differentially expressed genes between high and low m6Ascore group and visualized in volcano plot. B Intersection between differentially expressed genes from TCGA database and immune-related genes from Immport database showed in venn diagram. C Heatmap of the top 50 intersection immune-related genes with higher logFC. D Gene ontology enrichment analysis of intersection immune-related genes, BP biology process; CC cellular component; MF molecular function. E KEGG pathway enrichment analysis of intersection immune-related genes

Full size image

To identify the potential biological function of the m6A RNA methylation modulators that predicted the OS of patients with PAAD, we used GSEA to search for the associated Kegg pathways of m6A RNA methylation modulator genes (Fig. 6A-G). The enrichment score (ES), FDR, and nominal p-value were shown for each gene in Table 4. Although high expression of EIF3H expression was linked to cell cycle, no gene sets were significantly enriched in patients with high expression of METTL13. Gene sets were significantly enriched in Kegg pancreatic cancer with an upregulation of LRPPRC and IGF2BP3, we also observed a high expression of KIAA1429 in the Kegg pathways in cancer.

Fig. 6
figure 6

The associated KEGG pathways of m6A RNA methylation genes by GSEA. A EIF3H. B IGF2BP2. C IGF2BP3. D LRPPRC. E KIAA1429

Full size image
Table 4 gene set enrichment analysis for six m6A modulators with KEGG pathway using GSEA
Full size table

m6Ascore was related with the immune cell infiltration of PAAD

The composition of infiltrating immune cells between high and low m6Ascore group wsa explored by CIBERSORT algorithm. As for the tumor immune microenvironment in PAAD, immune scores and estimate scores were higher in PAAD patients with low m6Ascore (Fig. 7A-C), furthermore, we found a low reverse relationship between m6Ascore and immune score with r = − 0.303 (Fig. 7D-F), which pointed out that tumor purity was higher in the high m6Ascore group, indicating that patients with an unfavourable prognosis in the high m6Ascore group associated with the variation in tumor immune microenvironment of PAAD. Then we provided an insight into the correlations of m6Ascore with the 22 type of immune cells (Fig. 8A), patients with high m6Ascore had lower infiltration of Tregs and CD8+T cells but a higher resting CD4+ T infiltration, other immune cells showed no difference between the high and low groups. Furthermore, low CD8+T cells were significantly related to a poor OS, indicating that exhaustion of CD8+T cells that may lead to a poor OS in the high-risk group (Fig. 8B-D).

Fig. 7
figure 7

The association between tumor microenvironment and m6Ascore. Discrepancy of stromalscore (A), immunescore (B) and estimatescore (C) in two groups. Spearman’s correlation analysis of m6Ascore with tumor microenvironment. D stromalscore. E immunescore. F estimatedscore

Full size image
Fig. 8
figure 8

The different tumor-infiltrating immune cells between high and low m6Ascore patients with PAAD. A Profiles of 22 type of tumor-infiltrating immune cells in two group. The relationship between OS and A resting memory CD4+ T cells, B CD8+ T cells, C Tregs cells

Full size image

To predict the pancreatic cancer patients’ response to immunotherapy, we further investigated the difference in the expression of immune checkpoints (PD-1, PD-L1, CTLA-4, LAG3, TIGIT and TIM-3) between patients with low and high m6Ascore (Fig. 9A-F). Our results showed that patients with a low m6Ascore displayed a high abundance of PD-1, CTLA-4 and LAG3. IPS of 140 PAAD patients got from TCIA database was considered to be a superior predictor of response to anti-cytotoxic T lymphocyte antigen-4 (CTLA-4) and anti-programmed cell death protein 1 (anti-PD-1) antibodies, but IPS showed no discrepancy in two groups (Fig. 9H-K). Furthermore, three immunotherapy cohorts, including the IMvigor210 cohort, GSE78220 cohort and TCGA-SKCM cohort, were also used to investigate whether m6Ascore could predict patients’ responses to anti-PD-L1 therapy. in both GSE78220 and TCGA-SKCM cohort, patients with a low m6Ascore showed a high proportion of response to anti-PD-L1, but survival rate showed no difference in patients with high or low m6Ascore among all three anti-PD-L1 immunotherapy cohorts. Therapeutic effects of immunotherapies in patients with PAAD need to be verified using further studies because of the discrepancy in prediction results (Figure S2).

Fig. 9
figure 9

IPS and immunotherapy gene expression analysis. The levels of immune checkpoints molecules including PD- L1 (A), PD-1 (B), CTLA-4 (C), TIM-3 (D), LAG-3 (E), and TIGIT (F) in low-risk and high-risk groups. The association between IPS and the m6Ascore in PAAD patients based on TCIA database, H CTLA4− PD1− I CTLA4− PD1+ J CTLA4+ PD1− CTLA4+PD1+

Full size image

External validation of six key prognostic m6A RNA modulators

External pancreatic cancer databases from oncomine and GSE62452 were used to validate the expression of the six m6A RNA modulators in our prognosis model. In the oncomine database, expression of KIAA1429, LRPPRC, IGF2BP2, IGF2BP3, and EIF3H was higher in pancreatic cancer than normal tissues (Figure S3 A-F). In the GSE62452 dataset, we found a high expression of IGF2BP2 and IGF2BP3 in pancreatic tumor tissues compared with adjacent non-tumor tissue (Figure S4 A-F). Then we applied the risk score algorithm into the validation cohort to further verify the predictive value of the risk signature. In GSE62452 dataset, we also found that patients with a high m6Ascore showed poor prognosis, and the result from ROC analysis displayed that risk score had a predictive accuracy of prognosis, the AUC of 1-year, 3-year and 5-year survival was 0.61, 0.71, and 0.74, respectively (Figure S4 G-H).

IHC staining of six paired pancreatic cancer and adjacent non-tumorous tissues was used to validate the protein expression of the six prognosis genes (Fig. 10). The Compared with adjacent non-tumorous tissues, we found a higher level of EIF3H, IGF2BP2, IGF2BP3, KIAA1429 in all six pancreatic tumor tissues, LRPPRC was overexpressed in five pancreatic tumor tissues, while the expression of METTL3 decreased in four pancreatic tumor tissues.

Fig. 10
figure 10

IHC analysis of six m6A modulators in pancreatic cancer tissues and adjacent normal tissues

Full size image

Discussion

The high death rate of pancreatic cancer has been a continuous challenge in the field, it can progress rapidly with only mild symptoms, and there may be a multitude of causes, 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 multiple m6A methylation modulators in pancreatic cancer, identified the role of m6A methylation modulators signature in the TME cell infiltration and estimated the therapeutic effects of immunotherapy in pancreatic cancer.

We analysed 28 widely reported m6A RNA methylation modulators in pancreatic cancer according to the TCGA and GTEx database, a higher expression of these modulators except METTL3 was found in patients with PAAD compared to the normal tissues, suggesting an alterant level of m6A RNA methylation may play a carcinogenic role in pancreatic cancer. The prognosis model was conducted based on the expression of the selected methylation modulators, including KIAA1429, IGF2BP2, IGF2BP3, LRPPRC, METTL3, and EIF3H in our study. A similar study published in 2020 conducted a prognosis model based on the expression of KIAA1429, IGF2BP2, IGF2BP3, HNPNPC, METTL3, and YTHDF1. The discrepancy could be due to the difference of enrolled samples, but the similar finding was that METTL3 acts as a favourable factor while KIAA1429, IGF2BP2 and IGF2BP3 were risk factors in both pieces of research [36]. An overexpression of METTL3 in the MIA-PaCa-2 and BxPC-3 pancreatic cancer cell lines was involved in pancreatic carcinogenesis and promoted proliferation and invasion of pancreatic cancer cell lines in Xia T’s study [37]. Taketo K et al. pointed out that METTL3-depleted cells showed higher sensitivity to chemo- and radioresistance in pancreatic cancer cells [38], which was adverse with our analysis. IGF2BP2 and IGF2BP3, members of the mammalian IGF2 mRNA-binding protein family, were found correlated with an overall poor prognosis and metastasis of various cancer [39]. We found a higher pathologic stage and histologic grade associated with the overexpression of IGF2BP2 and IGF2BP3 in our study. The results from GSEA showed that overexpression of IGF2BP2 and IGF2BP3 was related to the progression of cancer, IGF2BP3 may involve rearrangement of peripheral actin, which led to the formation of additional membrane protrusions, promoting cell invasiveness and tumor metastasis of pancreatic cancer [40], IGF2BP2 was confirmed to be an oncogenic factor; miR-141 was a downstream regulatory gene in the progression of pancreatic cancer in a previous study [41]. EIF3H showed a function on cell cycle from GSEA in our study; however, the specific mechanism in PAAD need to be confirmed. Further, LRPPRC was a multifunctional protein involved in energy metabolism, and we found high LRPPRC enriched in pancreatic cancer, and it may play an essential role in pancreatic tumorigenesis by connecting an oncogenic gene expression with energy production. KIAA1429 overexpression is seen in many types of cancer that drive tumor growth and liver cancer cells’ metastasis by targeting GATA3 [42].

Tumorigenesis is a process that relies on the coordination of multiple factors, results showed that risk score was an independent prognosis indicator, and the combined expression of the selected modulators had a good prediction of 1-, 2-, and 3-year OS of PAAD, patients with high-risk score related with an unfavourable OS. Previous reports described m6A regulator-mediated methylation modification patterns played a non-negligible role in tumor microenvironment infiltration in gastric cancer [43]. Comparing the immune-related genes between high and low m6Ascore group, we found IL36B and IL36RN exhibited a higher level in the high m6Ascore group where IL36B was a gene involved TH17 pathway and played a role in polarizing T-helper responses. IL36RN encoded an interleukin-36-receptor antagonist [44] in renal cell carcinoma cells, and its overexpression inhibited proliferation, migration, invasion, and colony formation of tumor cells by suppressing β-catenin [45]. CD19, a target of chimeric antigen receptor (CAR) T-cell immunotherapy, was downregulated in PAAD patients with low m6Ascore. The results of GO and KEGG pointed out that different immune-related genes participated in the activation and differentiation of T cell and other immune-related pathways, revealing a significant role of m6Ascore signature in tumor’s immune environment. Besides, patients with high m6Ascore had a lower immune score and estimated score than patients with low m6Ascore, while a recent study declared that gliomas patients with high m6Ascore had a more abundance of immune infiltration with higher immune score and stromal score, which was conversed in PAAD [46]. As to the discrepancy of immune cell infiltration, patients with low m6Ascore had an elevation of Tregs and CD8+T cells but a low level of resting memory CD4+T cells, indicating an association between m6A modification and tumor immunity of PAAD. Moreover, a high level of CD8+T cell presented a high survival rate, where exhaustion of CD8+T cell may lead to an unfavourable OS in patients with high m6Ascore. Similar results in Han D’s study exhibited a deficiency of YTHDF1 that promoted the antigen-specific activation of CD8+ T cells and attributed to the enhanced capability of DCs to present tumor neoantigens [35]. Increasing evidence had shown that m6A RNA methylation modification plays a role in tumor immunity. Previous research reported that conduction of m6A signature was beneficial to reflect tumor microenvironment features of lung adenocarcinoma in different cohorts [47].

Immunotherapy had been shown been as an effective therapeutic strategy in diverse cancers, however, some patients faced an unsatisfactory effect due to individual variation. TME is composed of tumor vessels and different immune microenvironments, patients have different therapeutic responses due to the discrepancy of their structure. Checkpoint inhibitors have been the most mature and sufficient clinical research and the most widely used of tumor immunotherapy. Immune checkpoint therapy is a therapeutic method to inhibit tumor cells by regulating the activity of T cells through a series of pathways such as co-inhibition or co-stimulation. Various immune checkpoints have been used to treat PAAD, however, PD-1/PD-L1 and CTLA-4 blockade observed limited effect in PAAD treatment [48, 49]. We further predicted the efficacy of PD-1/PD-L1, CTLA-4, LAG-3, TIM-3 and TIGIT, and found a higher expression of PD-1, CTLA-4, and LAG-3 in the low m6Ascore group. Nevertheless, there was no difference in IPS between the two groups, and the estimated results from three immunotherapy cohorts indicated no association between m6Ascore and the curative effect of immunotherapy. Immunotherapy in pancreatic cancer is conflicted, and it needs more clinical trials to confirm the curative effect. For example, a phase I trial included one pancreatic patient who failed to show any efficacy with a treatment of anti-PD-L1 (MPDL-3280A) therapy [50]. In contrast, another clinical trial that enrolled 32 patients with pancreatic cancer showed the best response to MK-3475 and disease-free progression for 20 weeks after anti-PD-1 therapy [51]. Therapy targets immune checkpoints in pancreatic cancer showed a poor effect, and more efforts are required to overcome resistance to PD-1/PD-L1 targeted immunotherapy.

In general, our research found a discrepancy of m6A RNA modification modulators in pancreatic cancer based on TCGA and GTEx databased, and found an association between the expression of m6A RNA modification modulators and clinical characteristic. We conducted an m6A risk model and observed that high m6Ascore presented an unfavourable prognosis in PAAD patients. Besides, patients with high m6Ascore had a lower level of infiltration of CD8+ T cells, suggesting that m6A RNA modification may related to tumor microenvironment. However, the inadequate difference in IPS and m6Ascore showed no function in evaluating the therapeutic effects of checkpoint inhibitors treatment in PAAD. Furthermore, the exact mechanism of m6A modulators needs to be elucidated with future research.

Conclusion

m6A RNA methylation modulators play an important role in the progression of cancer. In this study, we generated a prediction model using 28 m6A RNA methylation modulators in PAAD and identified their association with the prognosis, and the prediction model was related to immune cell infiltration of PAAD. However, the function of the m6Ascore to evaluate the curative effect of immunotherapy needs to be figured out with further studies. Our study highlighted the involvement of m6A in pancreatic cancer opening up potential new research avenues and it also can be a potential measurement to diagnostic pancreatic cancer and to steer therapeutic decision making.

Availability of data and materials

The datasets analyzed during the current study are available from TCGA database (https://cancergenome.nih.gov/), GTEx database (https://xenabrowser.net/), Immport database (https://immport.niaid.nih.gov/), TCIA database (https://www.tcia.at/home), GEO database (https://www.ncbi.nlm.nih.gov/geo/), Oncomine database (https://www.oncomine.org).

References

  1. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–24. https://doi.org/10.3322/caac.21492. 

  2. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020;70(1):7–30. https://doi.org/10.3322/caac.21590.

    Article  PubMed  Google Scholar 

  3. Saad AM, Turk T, Al-Husseini MJ, Abdel-Rahman O. Trends in pancreatic adenocarcinoma incidence and mortality in the United States in the last four decades; a SEER-based study. BMC Cancer. 2018;18(1):688. https://doi.org/10.1186/s12885-018-4610-4.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Makohon-Moore A, Iacobuzio-Donahue CA. Pancreatic cancer biology and genetics from an evolutionary perspective. Nat Rev Cancer. 2016;16(9):553–65. https://doi.org/10.1038/nrc.2016.66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proc Natl Acad Sci U S A. 1974;71(10):3971–5. https://doi.org/10.1073/pnas.71.10.3971.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Desrosiers RC, Friderici KH, Rottman FM. Characterization of Novikoff hepatoma mRNA methylation and heterogeneity in the methylated 5′ terminus. Biochemistry. 1975;14(20):4367–74. https://doi.org/10.1021/bi00691a004.

    Article  CAS  PubMed  Google Scholar 

  7. Liu J, Yue Y, Han D, Wang X, Fu Y, Zhang L, et al. A METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nat Chem Biol. 2014;10(2):93–5. https://doi.org/10.1038/nchembio.1432.

    Article  CAS  PubMed  Google Scholar 

  8. Ping XL, Sun BF, Wang L, Xiao W, Yang X, Wang WJ, et al. Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell Res. 2014;24(2):177–89. https://doi.org/10.1038/cr.2014.3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Schwartz S, Mumbach MR, Jovanovic M, Wang T, Maciag K, Bushkin GG, et al. Perturbation of m6A writers reveals two distinct classes of mRNA methylation at internal and 5′ sites. Cell Rep. 2014;8(1):284–96. https://doi.org/10.1016/j.celrep.2014.05.048.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Wang X, Feng J, Xue Y, Guan Z, Zhang D, Liu Z, et al. Corrigendum: structural basis of N6-adenosine methylation by the METTL3-METTL14 complex. Nature. 2017;542(7640):260. https://doi.org/10.1038/nature21073.

    Article  CAS  PubMed  Google Scholar 

  11. Moindrot B, Cerase A, Coker H, Masui O, Grijzenhout A, Pintacuda G, et al. A pooled shRNA screen identifies Rbm15, Spen, and Wtap as factors required for Xist RNA-mediated silencing. Cell Rep. 2015;12(4):562–72. https://doi.org/10.1016/j.celrep.2015.06.053.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Patil DP, Chen CK, Pickering BF, Chow A, Jackson C, Guttman M, et al. M (6) a RNA methylation promotes XIST-mediated transcriptional repression. Nature. 2016;537(7620):369–73. https://doi.org/10.1038/nature19342.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, et al. VIRMA mediates preferential m6A mRNA methylation in 3'UTR and near stop codon and associates with alternative polyadenylation. Cell Discov. 2018;4(1):10. https://doi.org/10.1038/s41421-018-0019-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Knuckles P, Lence T, Haussmann IU, Jacob D, Kreim N, Carl SH, et al. Zc3h13/Flacc is required for adenosine methylation by bridging the mRNA-binding factor Rbm15/Spenito to the m6A machinery component Wtap/Fl (2) d. Genes Dev. 2018;32(5–6):415–29. https://doi.org/10.1101/gad.309146.117.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011;7(12):885–7. https://doi.org/10.1038/nchembio.687.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Tang C, Klukovich R, Peng H, Wang Z, Yu T, Zhang Y, et al. ALKBH5-dependent m6A demethylation controls splicing and stability of long 3′-UTR mRNAs in male germ cells. Proc Natl Acad Sci U S A. 2018;115(2):E325–33. https://doi.org/10.1073/pnas.1717794115.

    Article  CAS  PubMed  Google Scholar 

  17. Du H, Zhao Y, He J, et al. YTHDF2 destabilizes m (6) A-containing RNA through direct recruitment of the CCR4-NOT deadenylase complex. Nat Commun. 2016;7(1):12626. https://doi.org/10.1038/ncomms12626.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB. Dynamic m (6) a mRNA methylation directs translational control of heat shock response. Nature. 2015;526(7574):591–4. https://doi.org/10.1038/nature15377.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Liu N, Dai Q, Zheng G, He C, Parisien M, Pan T. N (6)-methyladenosine-dependent RNA structural switches regulate RNA-protein interactions. Nature. 2015;518(7540):560–4. https://doi.org/10.1038/nature14234.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Liu J, Dou X, Chen C, Chen C, Liu C, Xu MM, et al. N6-methyladenosine of chromosome-associated regulatory RNA regulates chromatin state and transcription. Science. 2020;367(6477):580–6. https://doi.org/10.1126/science.aay6018.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Berulava T, Buchholz E, Elerdashvili V, Pena T, Islam MR, Lbik D, et al. Changes in m6A RNA methylation contribute to heart failure progression by modulating translation. Eur J Heart Fail. 2020;22(1):54–66. https://doi.org/10.1002/ejhf.1672.

    Article  CAS  PubMed  Google Scholar 

  22. Livneh I, Moshitch-Moshkovitz S, Amariglio N, Rechavi G, Dominissini D. The m6A epitranscriptome: transcriptome plasticity in brain development and function. Nat Rev Neurosci. 2020;21(1):36–51. https://doi.org/10.1038/s41583-019-0244-z.

    Article  CAS  PubMed  Google Scholar 

  23. Gokhale NS, McIntyre A, Mattocks MD, et al. Altered m6A Modification of Specific Cellular Transcripts Affects Flaviviridae Infection. Mol Cell. 2020;77(3):542–555.e8.

    Article  CAS  PubMed  Google Scholar 

  24. Tran B, Nguyen ND, Center JR, Eisman JA, Nguyen TV. Association between fat-mass-and-obesity-associated (FTO) gene and hip fracture susceptibility. Clin Endocrinol. 2014;81(2):210–7. https://doi.org/10.1111/cen.12335.

    Article  CAS  Google Scholar 

  25. Cheng M, Sheng L, Gao Q, Xiong Q, Zhang H, Wu M, et al. The m6A methyltransferase METTL3 promotes bladder cancer progression via AFF4/NF-κB/MYC signaling network. Oncogene. 2019;38(19):3667–80. https://doi.org/10.1038/s41388-019-0683-z.

    Article  CAS  PubMed  Google Scholar 

  26. Liu T, Wei Q, Jin J, et al. The m6A reader YTHDF1 promotes ovarian cancer progression via augmenting EIF3C translation. Nucleic Acids Res. 2020;48(7):3816–31. https://doi.org/10.1093/nar/gkaa048.

  27. Li Y, Xiao J, Bai J, Tian Y, Qu Y, Chen X, et al. Molecular characterization and clinical relevance of m6A regulators across 33 cancer types. Mol Cancer. 2019;18(1):137. https://doi.org/10.1186/s12943-019-1066-3.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Chen M, Wong CM. The emerging roles of N6-methyladenosine (m6A) deregulation in liver carcinogenesis. Mol Cancer. 2020;19(1):44. https://doi.org/10.1186/s12943-020-01172-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Chen Y, Wang J, Xu D, et al. m6A mRNA methylation regulates testosterone synthesis through modulating autophagy in Leydig cells. Autophagy. 2020:1–19. https://doi.org/10.1080/15548627.2020.1810918.

  30. Tang C, Xie Y, Yu T, Liu N, Wang Z, Woolsey RJ, et al. m6A-dependent biogenesis of circular RNAs in male germ cells. Cell Res. 2020;30(3):211–28. https://doi.org/10.1038/s41422-020-0279-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Ni W, Yao S, Zhou Y, Liu Y, Huang P, Zhou A, et al. Long noncoding RNA GAS5 inhibits progression of colorectal cancer by interacting with and triggering YAP phosphorylation and degradation and is negatively regulated by the m6A reader YTHDF3. Mol Cancer. 2019;18(1):143. https://doi.org/10.1186/s12943-019-1079-y.

    Article  PubMed  PubMed Central  Google Scholar 

  32. He Y, Hu H, Wang Y, Yuan H, Lu Z, Wu P, et al. ALKBH5 inhibits pancreatic Cancer motility by decreasing long non-coding RNA KCNK15-AS1 methylation. Cell Physiol Biochem. 2018;48(2):838–46. https://doi.org/10.1159/000491915.

    Article  CAS  PubMed  Google Scholar 

  33. Tang B, Yang Y, Kang M, Wang Y, Wang Y, Bi Y, et al. m6A demethylase ALKBH5 inhibits pancreatic cancer tumorigenesis by decreasing WIF-1 RNA methylation and mediating Wnt signaling. Mol Cancer. 2020;19(1):3. https://doi.org/10.1186/s12943-019-1128-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Lan H, Liu Y, Liu J, Wang X, Guan Z, du J, et al. Tumor-associated macrophages promote Oxaliplatin resistance via METTL3-mediated m6A of TRAF5 and necroptosis in colorectal Cancer. Mol Pharm. 2021;18(3):1026–37. https://doi.org/10.1021/acs.molpharmaceut.0c00961.

    Article  CAS  PubMed  Google Scholar 

  35. Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, et al. Anti-tumour immunity controlled through mRNA m6A methylation and YTHDF1 in dendritic cells. Nature. 2019;566(7743):270–4. https://doi.org/10.1038/s41586-019-0916-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Hou J, Wang Z, Li H, Zhang H, Luo L. Gene signature and identification of clinical trait-related m6 a regulators in pancreatic Cancer. Front Genet. 2020;11:522. https://doi.org/10.3389/fgene.2020.00522.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Xia T, Wu X, Cao M, Zhang P, Shi G, Zhang J, et al. The RNA m6A methyltransferase METTL3 promotes pancreatic cancer cell proliferation and invasion. Pathol Res Pract. 2019;215(11):152666. https://doi.org/10.1016/j.prp.2019.152666.

    Article  CAS  PubMed  Google Scholar 

  38. Taketo K, Konno M, Asai A, Koseki J, Toratani M, Satoh T, et al. The epitranscriptome m6A writer METTL3 promotes chemo- and radioresistance in pancreatic cancer cells. Int J Oncol. 2018;52(2):621–9. https://doi.org/10.3892/ijo.2017.4219.

    Article  CAS  PubMed  Google Scholar 

  39. Lederer M, Bley N, Schleifer C, Hüttelmaier S. The role of the oncofetal IGF2 mRNA-binding protein 3 (IGF2BP3) in cancer. Semin Cancer Biol. 2014;29:3–12. https://doi.org/10.1016/j.semcancer.2014.07.006.

    Article  CAS  PubMed  Google Scholar 

  40. Taniuchi K, Furihata M, Hanazaki K, Saito M, Saibara T. IGF2BP3-mediated translation in cell protrusions promotes cell invasiveness and metastasis of pancreatic cancer. Oncotarget. 2014;5(16):6832–45. https://doi.org/10.18632/oncotarget.2257.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Xu X, Yu Y, Zong K, Lv P, Gu Y. Up-regulation of IGF2BP2 by multiple mechanisms in pancreatic cancer promotes cancer proliferation by activating the PI3K/Akt signaling pathway. J Exp Clin Cancer Res. 2019;38(1):497. https://doi.org/10.1186/s13046-019-1470-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Lan T, Li H, Zhang D, Xu L, Liu H, Hao X, et al. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3. Mol Cancer. 2019;18(1):186. https://doi.org/10.1186/s12943-019-1106-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020;19(1):53. https://doi.org/10.1186/s12943-020-01170-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Gresnigt MS, van de Veerdonk FL. Biology of IL-36 cytokines and their role in disease. Semin Immunol. 2013;25(6):458–65. https://doi.org/10.1016/j.smim.2013.11.003.

    Article  CAS  PubMed  Google Scholar 

  45. Zhu K, Sun D, Zou X, Liu R, Wang Z. Interleukin-36 receptor antagonist is associated with the progression of renal cell carcinoma. Int Immunopharmacol. 2020;84:106474. https://doi.org/10.1016/j.intimp.2020.106474.

    Article  CAS  PubMed  Google Scholar 

  46. Xu S, Tang L, Dai G, Luo C, Liu Z. Expression of m6A regulators correlated with immune microenvironment predicts therapeutic efficacy and prognosis in gliomas. Front Cell Dev Biol. 2020;8:594112. https://doi.org/10.3389/fcell.2020.594112.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Li Y, Gu J, Xu F, et al. Molecular characterization, biological function, tumor microenvironment association and clinical significance of m6A regulators in lung adenocarcinoma. Brief Bioinform. 2020. pp. bbaa225. https://doi.org/10.1093/bib/bbaa225. Epub ahead of print.

  48. Royal RE, Levy C, Turner K, Mathur A, Hughes M, Kammula US, et al. Phase 2 trial of single agent Ipilimumab (anti-CTLA-4) for locally advanced or metastatic pancreatic adenocarcinoma. J Immunother. 2010;33(8):828–33. https://doi.org/10.1097/CJI.0b013e3181eec14c.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Brahmer JR, Tykodi SS, Chow LQ, et al. Safety and activity of anti-PD-L1 antibody in patients with advanced cancer. N Engl J Med. 2012;366(26):2455–65. https://doi.org/10.1056/NEJMoa1200694.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Herbst RS, Soria JC, Kowanetz M, Fine GD, Hamid O, Gordon MS, et al. Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature. 2014;515(7528):563–7. https://doi.org/10.1038/nature14011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Shimizu T, Seto T, Hirai F, Takenoyama M, Nosaki K, Tsurutani J, et al. Phase 1 study of pembrolizumab (MK-3475; anti-PD-1 monoclonal antibody) in Japanese patients with advanced solid tumors. Investig New Drugs. 2016;34(3):347–54. https://doi.org/10.1007/s10637-016-0347-6.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thankful for the assistance from Yuan Zhang, Weisheng Cheng, Sheng Fang and Zhen Wang in proceeding data.

Funding

This study was supported by the Anhui Natural Science Foundation (grant number:9021138205).

Author information

Authors and Affiliations

  1. Department of Clinical Laboratory, the First Affiliated Hospital of Anhui Medical University, No. 218 Jixi Road, Hefei, 230032, Anhui, China

    Lianzi Wang, Shubing Zhang, Huimin Li, Yang Xu, Tao Li & Yuanhong Xu

  2. Department of Pathology, the First Affiliated Hospital of Anhui Medical University, No. 218 Jixi Road, Hefei, 230032, Anhui, China

    Qiang Wu

  3. The Key Laboratory of Microbiology and Parasitology of Anhui Province, the Key Laboratory of Zoonoses of High Institutions in Anhui, Anhui Medical University, Hefei, China

    Jilong Shen

Authors
  1. Lianzi Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Shubing Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Huimin Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yang Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Qiang Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Jilong Shen
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Tao Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Yuanhong Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

LZW, TL and YHX involved in the study design and writing of the article; HML and YX involved in downloading the data; SBZ involved in the data analysis; QW provided clinical tissue samples; JLS provided experimental environment for IHC; YX involved in the evaluation and analysis the results of IHC. All authors reviewed and approved the final manuscript.

Corresponding authors

Correspondence to Tao Li or Yuanhong Xu.

Ethics declarations

Ethics approval and consent to participate

This work was approved by the Academic Committee of The First Affiliated Hospital of Anhui Medical University and was conducted following the Declaration of Helsinki. All patients provided written informed consent.

Consent for publication

Not applicable.

Competing interests

The authors report no conflicts of interest in this work. None of the authors have any personal or financial involvement with organizations that have financial interest in the content of this manuscript.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

The information of excluded and enrolled samples of PAAD from TCGA database.

Additional file 2: Figure S1.

The association between the expression of the modulators and survival time with statistically significant (A) METTL3. (B) IGF2BP3. (C) IGF2BP2. (D) KIAA1429. (E) EIF3H. (F) LRPPRC. Figure S2. Predictive value of m6Ascore in anti-PD-1/L1 immunotherapy based on three immunotherapeutic cohorts. (A) Survival analyses for patients with high or low m6Ascore in GSE78220 cohort. (B) The proportion of patients with response to anti-PD-1/L1 immunotherapy in patients with high or low m6Ascore in GSE78220 cohort. (C) Survival analyses for patients with high or low m6Ascore in TCGA-SKCM cohort (D) The proportion of patients with response to anti-PD-1/L1 immunotherapy in patients with high or low m6Ascore in TCGA-SKCM cohort. (E) Survival analyses for patients with high or low m6Ascore in IMvigor210 cohort (F) The proportion of patients with response to anti-PD-1/L1 immunotherapy in patients with high or low m6Ascore in IMvigor210 cohort. Figure S3. Expression of of six m6A modulators in pancreatic cancer tissues and adjacent normal tissues regarding to oncomine database (A) EIF3H. (B) IGF2BP2. (C) IGF2BP3. (D) KIAA1429. (E) METTL3. (F) LRPPRC. Figure S4. External Validation of six Key Prognostic m6A RNA modulators in GSE62452 dataset. Expression of of six m6A modulators in pancreatic cancer tissues and adjacent normal tissues based on GSE62452 dataset. (A) METTL3. (B) LRPPRC. (C) KIAA1429. (D) IGF2BP3. (E) IGF2BP2. (F) EIF3H. (G) Time-dependent ROC analysis of m6Ascore in predicting prognosis. (H) Survival analyses for low and high m6Ascore patient groups using Kaplan-Meier curves.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, L., Zhang, S., Li, H. et al. Quantification of m6A RNA methylation modulators pattern was a potential biomarker for prognosis and associated with tumor immune microenvironment of pancreatic adenocarcinoma. BMC Cancer 21, 876 (2021). https://doi.org/10.1186/s12885-021-08550-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-021-08550-9

Keywords

BMC Cancer

ISSN: 1471-2407

Contact us