Molecular characterization and prognostic relevance of m6A regulators in pancreatic adenocarcinoma

doi:10.21203/rs.3.rs-4948418/v1

Download PDF

Article

Molecular characterization and prognostic relevance of m6A regulators in pancreatic adenocarcinoma

https://doi.org/10.21203/rs.3.rs-4948418/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

RNA N6-methyladenine (m6A) modification played an essential role in the occurrence and development of malignant tumors. m6A modification patterns in immune response and tumor microenvironment (TME) remains an enigma.

Methods

25 m6A regulators were collected, the molecular alterations and clinical relevance of which were explored. The mutation landscape of the pancreatic adenocarcinoma (PAAD) patients was explored by using TCGA data. The expression difference of the m6A regulators was identified by TCGA and HPA data. The prognosis value of the m6A regulators was measured by TCGA and ICGC data. Consensus clustering analysis was used for different m6A modification patterns identification. CIBERSORT and ESTIMATE algorithms were used to explore the landscape of TME cell infiltration. DEG analysis was used for m6A-related gene identification. m6A-score signature was established by using univariate Cox regression analysis and PCA.

Results

CNV amplification of m6A regulators led to up-regulated of them in tumor tissues in comparison with normal tissues. 13 of the 25 regulators showed oncogenic features. Two distinct m6A modification patterns were defined. PAAD patients in m6Acluster A occupied better survival compared to m6Acluster B. The relationships between the two m6A patterns and different types of immune infiltrating cells were further identified. A consolidated scoring system to quantify the m6A modification pattern of individual patients was established. Patients in low m6A-score group had better OS compared with these in high m6A-score group. Subsequent analysis proved that m6A methylation modification patterns was associated with response to anti-PD-L1 immunotherapy.

Conclusions

The molecular alterations and prognostic implications of m6A regulators were analyzed. The distinct m6A modification patterns are crucial for understanding the heterogeneity and complexity of individual tumor microenvironments (TMEs). A comprehensive assessment of m6A modification in tumors enhances our understanding of TME infiltration characteristics and facilitates more effective immunotherapy strategies.

Health sciences/Biomarkers/Predictive markers

Health sciences/Biomarkers/Prognostic markers

Health sciences/Gastroenterology/Gastrointestinal diseases

Health sciences/Gastroenterology/Gastrointestinal system

m6A modification pattern

tumor microenvironment

pancreatic adenocarcinoma

prognosis

In recent years, with the development of ribonucleic acid (RNA) deep sequencing technology and bioinformatics methods, the regulatory role of biological epigenetic modification in the occurrence and development of malignant tumors is set to gain more and more attention(Nombela et al., 2021). Among these modifications, the N6-methyladenine (m6A) RNA epigenetic modification has gradually become a research hotspot in the field of biological science(Imam et al., 2020; Nombela et al., 2021). m6A modification is a methyltransferase catalyzed methylation of adenine at 6N position, which is the most common epigenetic modification of mRNA in all higher eukaryotes(Liu et al., 2014). It is mainly regulated by m6A methyltransferase, demethylase and RNA binding protein(Meyer and Jaffrey, 2014). At the level of epigenetic modification, m6A modification impacts on the occurrence and development of tumors through regulating the expressions of proto oncogene and anti-oncogene(Meyer and Jaffrey, 2014). The whole process is dynamic and reversible, which involves the fine regulation of many complex cell processes(Liu et al., 2018). By regulating the proliferation, invasion and metastasis, drug resistance and stem cell stem of tumor cells, m6A modification could play an essential role in tumor process and progression(Liu et al., 2018). Some m6A RNA methylation regulators have been reported to be significantly associated with progression and prognosis of malignancies(Tang et al., 2020). Methyltransferases METTL3 and METTL14, demethylases FTO and ALKBH5 have been proved to abnormally expressed in digestive tract tumors, which mainly enhance the metastasis and invasion ability of tumor cells by affecting the stability of mRNA, regulating the proliferation of cancer cells, affecting the extracellular matrix or changing the EMT ability of cells(Tang et al., 2020). In fact, there is growing evidence that expression disorder and genetic variation of m6A regulators paly indispensable roles in a variety of biological process, such as tumor malignant progression and immunomodulatory abnormality(Fang and Declerck, 2013; Binnewies et al., 2018; Han et al., 2019).

Immunotherapy is a type of cancer treatment, which helps immune system fight cancer(Marabelle et al., 2017). At present, more and more studies have indicated that tumors could be treated by immunotherapy effectively and safely(Martini et al., 2019; Sprooten et al., 2019; Xie et al., 2019). Pancreatic adenocarcinoma (PAAD) mainly refers to a malignant tumor derived from pancreatic duct epithelial cells and follicular cells, which is a common malignant tumor of digestive system(Bengtsson et al., 2020). According to the data of the National Comprehensive Cancer Network Center in 2020, about 540000 people were diagnosed with PAAD in the United States in 2019, and the number of deaths due to PAAD was about 430000(Siegel et al., 2020). The incidence rate of PAAD is the top ten among the leading cancers in the United States, ranking third in the cancer-related death causes, as the American Cancer Society reported(Siegel et al., 2020). Although the therapy of PAAD has made progress with the continuous efforts of scientists, its curative effect is unsatisfactory. Increasing evidence indicated that tumor microenvironment (TME) also non-negligibly impacted on tumor progression, which included cancer cells, stromal cells, distant recruited cells, bone marrow-derived cells, and secreted factors(Pinello et al., 2018; Tong et al., 2018; Chokkalla et al., 2020). Thus, by comprehensively exploring the TME landscape heterogeneity and complexity, we might improve the guiding and prediction potential of immunotherapeutic responsiveness by further screening out different tumor immune subtypes. Recently, a variety of studies have focused on exploring the relationship between TME infiltrating immune cells and m6A modification. However, the role of m6A methylation in PAAD has been only preliminarily recognized. Thus, a comprehensive understanding of the invasion characteristics of TME cell infiltration mediated by multiple m6A regulators will help us to deepen our understanding of the immune regulation of TME.

In the present study, the genomic information of PAAD samples was retrieved from (The Cancer Genome Atlas) TCGA database, which was further `used to evaluate the m6A modification patterns in PAAD. The correlation between the m6A modification pattern and the invasion characteristics of TME cells was immediately analyzed. We revealed two different m6A modification patterns in PAAD, and unexpectedly found that the characteristics of TME in these two patterns were highly associated with the relative abundance of TME-infiltrating cells and prognosis of PAAD patients, indicating that m6A modification might play an essential role in the formation of PAAD tumor microenvironment. Finally, to quantify the m6A modification pattern of individual patients, we constructed a scoring system.

Collection of m6A regulators

Figure 1 showed the workflow of the present study. We firstly extracted m6A RNA methylation regulators by conducting a m6A-related literature review(Batista, 2017; Liu et al., 2017; Dai et al., 2018; Chai et al., 2019; Zhang et al., 2020). In total, 25 m6A regulators including 15 readers (EIF3A, ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, LRPPRC, RBMX, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3), 8 writers (CBLL1, METTL14, METTL3, RBM15, RBM15B, VIRMA, WTAP, ZC3H13, ALKBH5), and 2 erasers (ALKBH5, FTO) were included in this study. The basic information of the included m6A regulators were shown in Table S1.

Collection and preprocessing of pancreatic adenocarcinoma genome-related data

PAAD microarray data which displayed as Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values were retrieved from TCGA database (https://genomecancer.ucsc.edu/). We further transformed FPKM values into TPM values and removed tumor samples without complete clinical information from subsequent analysis. In total, 179 PAAD samples were included in this study. The adjusted normal samples (n = 171) were also obtained from this database and another database (Genome Tissue Expression (GTEX)). Expression differences of m6A regulators in PAAD samples and normal samples were measured by student t test or Wilcoxon rank-sum test. The single-nucleotide variants (SNV) data of PAAD which containing 158 samples were also obtained from TCGA database. The mutation landscape in PAAD patients was presented by using R package “maftools”(Mayakonda et al., 2018). Furthermore, we downloaded Copy Number Variation (CNV) data from cBioPortal (http://www.cbioportal.org/)(Cerami et al., 2012). GISTIC was applied to calculate the CNV variation type (gain or loss) and frequency of m6A regulators in TCGA-PAAD cohort. R package “likert” was used for the visualization of CNV analysis result. We also downloaded the PAAD samples and related survival information from the International Cancer Genome Consortium (ICGC) database (https://daco.icgc.org/) for prognostic value validation of m6A regulators. Furthermore, the protein expression levels of m6A regulators were obtained by IHC retrieving from The Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). 6 to 12 PAAD samples from the 25 m6A regulators were included and analyzed. IHC staining was indicated as four levels (high-, medium-, low- staining or not detected) automatically by this database.

Function enrichment analysis and cross talks among m6A regulators

To explore the potential functions of m6A regulators in PAAD, we firstly obtained function enrichment analysis including GO (Gene ontology) and KEGG (Kyoto encyclopedia of genes and genomes) based on package “clusterProfiler”(Yu et al., 2012) in R software. Also, in order to explore the interactions among the m6A regulators, we calculated the Pearson correlation coefficient among m6A regulators by using their own expression in TCGA-PAAD data. R package “corrplot” was used for visualization. Furthermore, a protein-protein interaction (PPI) network was constructed for the 25 m6A regulators by using the Search Tool for the Retrieval of Interacting Genes (STRING)(Szklarczyk et al., 2019). The degree of genes was further calculated by an appliance in Cytoscape(Shannon et al., 2003) software named network analyzer to quantify the relationships among m6A regulators.

m6A regulator prognostic relevance

To further explore the association between m6A regulators and PAAD patient survival, we performed survival analysis for each m6A regulator by using two independent PAAD cohort (TCGA-PAAD data, ICGC data). In each dataset, we split PAAD samples into high-/low- expression groups by setting the centermost element of each m6A regulator expression as the cut-off criteria. The P values were calculated by using R package “survival”(Therneau, 2015) and survival curves were drawn by using R package “forestplot”(Aut and Aut, 2016).

Consensus clustering for m6A regulators

Then based on the expression data profile of the 25 m6A regulators from TCGA-PAAD data, we performed consensus clustering analysis to explore the m6A modification patterns. To guarantee the stability of classification, we conducted 1000 iterations and kept 80% resampling rate. R package “ConsensuClusterPlus”(Wilkerson and Hayes, 2010) was chosen to achieve the above procedures. After choosing an appropriate number for clustering (k) by R package “Nbclust”, patients were further classified into m6A-related groups for subsequent analysis (m6A cluster A, B, C, etc). Principal component analysis (PCA) was performed to explore the difference between different modification patterns by using R package “factoextra”. Survival analysis between different m6A-related groups was also obtained by R package “survival”(Therneau, 2015; Aut and Aut, 2016). We also quantified the expression difference of m6A regulators between different modification patterns. For statistical significance, Wilcoxon rank-sum test was recommended for two groups and Kruskal Wallis test was chosen for multiple independent sets.

Landscape of tumor microenvironment (TME) cell infiltration between different m6A distinct phenotypes

The relative abundance of 22 kinds of TME-infiltrating cells in each PAAD tissue from TCGA-PAAD were quantified by using CIBERORT(Newman et al., 2015) algorithm. As a deconvolution algorithm, CIBERSORT infers cell type proportions by using support vector regression with a set of reference gene-expression values (a signature with 547 genes). We further explored the relationships between different TME-infiltrating cell types and measured the difference of relative abundance of TME cells between different m6A modification patterns. Furthermore, ESTIMATE scores, immune scores and stromal scores of these PAAD samples were firstly calculated by applying the ESTIMATE algorithm based on R package “estimate”(Yoshihara et al., 2013). We immediately explored the association between these scores and different m6A distinct phenotypes. Survival analysis for subgroup PAAD patients stratified by both m6A modification patterns and scores (immune/stromal) was also conducted as we previous did.

m6A-related gene identification

After classifying PAAD patients into different m6A distinct phenotypes, we further screened out differentially expressed genes (DEGs) between different m6A modification patterns by using R package “edgeR”(Robinson et al.). In the present study, adjusted P value < 0.001 and |log2FC| ≥ 1.5 were set as the standards to identify DEGs, which were regarded as m6A-related genes for subsequent analysis.

m6A gene signature generation

The overlapped m6A-related genes identified from different m6A modification patterns were firstly extracted. Then we also classified PAAD patients into several gene clusters by applying consensus clustering analysis as we previous did. Furthermore, univariate Cox regression analysis was conducted for each overlapped m6A-related gene by using R package “survival”. In the present study, we thought a m6A-related gene with P value < 0.05 was significantly associated with prognosis of patients with PAAD. Relying on the identified prognosis-related genes, principal component analysis (PCA) was further obtained by using R package “facroextra”, principal component 1 (PC1) was selected for signature score calculation. PCA is a classic multivariate statistical analysis method, which has obvious advantage in the integrated multi-variable problem. According to the HR (results from univariate Cox regression analysis) of each prognosis-related gene, a method similar to GGI was used for m6A-score calculation as follows: m6A-score = ∑PC1_i × Exp_i - ∑PC1_j × Exp_j. In which Exp represents the expression value of each prognosis-related gene, i represents the signature score of clusters whose HR is less than 1, j represents the signature score of clusters whose HR is more than 1. After defining the m6A-score, we compared the m6A-score difference between different m6Aclusters, different gene clusters. An alluvial diagram was drawn by using R package “ggplot2” to show the change of m6Aclusters, gene clusters, survival status, and m6A-score.

Relationship between m6A modification patterns and anti-PD-1/L1 immunotherapy

Two immunotherapeutic cohorts were collected and used in the present study. The related clinical information was also retrieved. One of the cohorts was IMvigor210 cohort(Mariathasan et al., 2018) which contained advanced urothelial cancer with intervention of atezolizumab. Atezolzumab acted as an anti-PD-L1 antibody. The expression data of IMvogor210 cohort displayed as count number was immediately retrieved from http://research-pub.Gene.com/imvigor210corebiologies. We then transformed the count value into the TPM value by using R package “DEseq2”. After that, totally 298 samples were included for subsequent analysis. Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) was a public functional genomics data repository including gene expression profiles. We further downloaded GSE78220 cohort(Hugo et al., 2016) including metastatic melanoma samples with pembrolizumab, an anti-PD-1 antibody from this database. We used R package “limma”(Ritchie et al., 2015) for the normalization, the FPKM value was also transformed into the TPM value. In total, 27 samples were used in the subsequent analysis. Also, the related clinical information was collected. The survival curves for prognostic analysis were generated by using R package “survival”. And receiver operating characteristic (ROC) curves were also drawn and the area under the curve (AUC) were further calculated by using R package “pROC”(Robin et al., 2011).

Cross talk of genetic variation of m6A regulators in PAAD

After reviewing the previous literature, 25 m6A regulators were screened out in total (Fig. 2A). These m6A regulators could been classified in three types, which contained 15 readers, 8 writers, and 2 erasers (Fig. 2A). The CNV alteration frequency for the 25 m6A regulators was investigated. CNV alterations of m6A regulators were prevalent. Several genes including VIRMA, HNRNPA2B1, IGF2BP3, YTHDF3, YTHDF1 showed strong frequency of amplification in copy number (Fig. 2B). On the contrary, WTAP, YTHDF2, and RBM15B occupied a widespread frequency of CNV deletion (Fig. 2B). We immediately identified frequently mutated genes and further summarized the incidence of somatic mutations of 25 m6A regulators in PAAD. The top 30 high-frequency mutated genes were shown in Fig. 2C. These genes altered in 124 (78.48%) of 158 samples. Among them, the top 5 high-frequency mutated genes were TP53 (54%), KRAS (54%), CDKN2A (17%), SMAD4 (17%), and TTN (11%). Among 158 samples, only 8 (5.06%) experienced mutations of m6A regulators (Fig. 2D). In general, for m6A regulators, the overall average mutation frequency was low, scoping from 0.00–1.00% (Fig. 2D).

Expression difference of m6A regulators between PAAD tissues and normal tissues

To explore the effect of genetic alterations in m6A regulators expression, the expression difference of m6A regulators between PAAD tissues and normal tissues was subsequently screened out. The overall average expression level was higher, compared to normal tissues. In detail, the mRNA expressions of 8 readers including ELAVL1, HNRNPA2B1, HNRNPC, IGF2BP2, IGF2BP3 were significantly increased in PAAD tissues compared to normal tissues (P < 0.05, Fig. 3A). Moreover, 4 writers including METTL14, RBM15, VIRMA, and ZC3H13 were significantly high expressed in PAAD samples compared to normal samples (P < 0.05, Fig. 3B). The erasers mRNA (ALKBH5 and FTO) were also over expressed in tumor samples (P < 0.05, Fig. 3C). Interestingly, we found that CNV amplification of m6A regulators led to up-regulated of them in tumor tissues in comparison with normal tissues (such as VIRMA, HNRNPA2B1, IGF2BP3, YTHDF3), which suggested that CNV alterations might be one of the factors influencing the mRNA expression level of m6A regulators. IHC data containing 6 to 12 independent PAAD samples were further collected from HPA database (Fig. 3D). YTHDF1, YTHDF3, and METTL3 were not evaluated. Overall, m6A readers including HNRNPA2B1, HNRNPC, RBMX, and YTHDC1 showed high staining. But we could only observe expression of IGF2BP1, IGF2BP3, YTHDF2 in 5/12, 6/12, 2/8 samples, respectively. Almost all the writers including METTL14, RBM15, RBM15B, VIRMA, and WTAP were not detected in several samples (4/12, 2/11, 10/11, 3/11, 3/11, separately). We observed that 12/12 and 11/12 samples expressed ALKBH5, and FTO erasers m6A regulators, separately (Fig. 3D).

Cross talks among m6A regulators

We further explored m6A regulators related biological process (BPs) and KEGG signaling pathways, for better comprehending the roles of m6A regulators in PAAD. The results demonstrated that m6A regulators significantly enriched in regulation of mRNA metabolic process, mRNA processing, RNA catabolic process, RNA splicing, regulation of mRNA stability, regulation of RNA stability, regulation of mRNA catabolic process, mRNA catabolic process, RNA modification, and RNA methylation (Table S2, Fig. 4A). And the 25 m6A regulators were only regulated in two KEGG signaling pathways including spliceosome and RNA transport (Fig. 4B). These results demonstrated that m6A regulators played nonnegligible roles in m6A modification. We immediately explored the co-occurrence of genetic alterations and expression correlation among m6A regulators. Surprisingly, the result demonstrated that there was also a high association between different types of m6A regulators (Fig. 4C). For example, IGF2BP2 reader was significantly associated with eraser ALKBH5 (R = -0.62, P < 0.001, Fig. 4C). Then we constructed a PPI network of m6A regulators (Fig. 4D). Though classified into different types of m6A regulators, readers, writers, and erasers interacted with each other closely in the PPI network. We further calculated the degree of connectivity of each m6A regulators. Among all the readers, HNRNPC (degree = 17) and HNRNPA2B1 (degree = 17) showed the highest degree (Fig. 3D). VIRMA (degree = 17), WTAP (degree = 17), and METTL14 (degree = 17) showed strongest association with each other among all the writers (Fig. 4D). The two erasers (ALKBH5: degree = 15; FTO: degree = 11) also frequently associated with other m6A regulators (Fig. 4D). Summary above, we thought these m6A regulators might work with each other in PAAD.

Survival landscape of m6A regulators in PAAD

To investigate the association between the m6A regulator expression and prognosis of PAAD, survival analysis was further conducted. Overall, half of the m6A regulators (13/25) showed significant correlation with PAAD patient survival by using TCGA data (Fig. 3E). In detail, all the 13 genes showed oncogenic features, PAAD patients in high expression group of these genes occupied obviously worse survival compared to these in low expression group (Fig. 4E). Seven m6A regulators including 5 readers (FMR1 (TCGA: P = 0.029; ICGC: P = 0.003), IGF2BP2 (TCGA: P < 0.001; ICGC: P = 0.028), IGF2BP3 (TCGA: P < 0.001; ICGC: P = 0.005), RBMX (TCGA: P = 0.016; ICGC: P = 0.028), YTHDC1 (TCGA: P = 0.013; ICGC: P < 0.001)) and 2 writers (METTL14 (TCGA: P = 0.028; ICGC: P = 0.017), VIRMA (TCGA: P = 0.001; ICGC: P = 0.016)) were validated to show the similar correlation with PAAD patient survival by using ICGC data (Fig. 4F). Unfortunately, none of erasers showed significant association with PAAD patient survival (Fig. 4E-F), which pointed out that types of m6A regulators might affect their roles in prognosis of PAAD patients.

m6A modification patterns defined via 25 m6A regulators

Based on R package “ConsensuClusterPlus”, unsupervised hierarchical clustering of the 179 tumors with matched m6A regulator expression profiles from TCGA data was obtained. Finally, as shown in Fig. 5A, 179 tumors were classified into two distinct modification patterns (80 samples in m6Acluster A, 99 samples in m6Acluster B). PCA demonstrated that there was significant distinction existed between different m6A modification patterns (Fig. 5B). Further analysis suggested that there was prognosis difference between m6A cluster A and m6Acluster B. PAAD patients in m6Acluster A occupied better survival compared to m6Acluster B (P = 0.02, Fig. 5C). We immediately explored the expression difference of m6A regulators between different m6A clusters, the results were shown as a heatmap (Fig. 5D). Both the two erasers were higher expressed in m6AclusterA compared to m6Acluster B (ALKBH5: P < 0.0001; FTO: P < 0.05). Several readers including IGF2BP1 (P < 0.0001), IGF2BP2 (P < 0.0001), and IGF2BP3 (P < 0.0001) were significantly up-regulated in m6Acluster B while RBMX (P < 0.05), YTHDC1 (P < 0.05), YTHDC2 (P < 0.05), and YTHDF1 (P < 0.05) were obviously down-regulated in m6Acluster B. The expressions of writers were consistent between different m6A modification patterns. In detail (Fig. 5D), METTL14, METTL3, RBM15B, and WTAP occupied a significant down-regulation in m6Acluster B across different clusters.

Association between m6A modification patterns and tumor microenvironment

As shown in Fig. 6A, after measuring the relative abundance of 22 TME-infiltrating cells, a stacked bar plot was used to show the relative proportion of these cell types in each PAAD sample. The results demonstrated that the relative proportion of the 22 cells were quite different across the PAAD samples. Obviously, some TME-infiltrating cells showed strong correlation with each other (Fig. 6B). For example, B cells memory showed negatively association with T cells CD4 naïve. T cells CD8 had strong correlation with Macrophages M0, positively. More relationships between different types of TME-infiltrating cells could be explored in Fig. 6B. What surprised us was that the immune infiltration levels of several immune cells were quite different between different m6A clusters (Fig. 6C). Macrophages M0 (P < 0.05), and Macrophages M1 (P < 0.05) were significantly enriched in m6Acluster B group, compared to m6Acluster A group (Fig. 6C). In addition, Monocytes, and T cells CD8 were enriched in m6Acluster A group when compared with m6Acluster B group (Fig. 6C). We further calculated immune scores and stromal scores of PAAD samples. Unfortunately, we could not explore significant relationship between these scores and m6A modification patterns. Not only immune scores (Fig. 6D) but also stromal scores (Fig. 6E) did not show significant expression difference between different m6A clusters. We further explored the influence of immune/stromal score on prognosis prediction by m6A modification patterns. In PAAD patients of low immune score, we found PAAD patients in m6Acluster A occupied better survival compared with these in m6Acluster B (Fig. 6F, P = 0.038). Furthermore, there was no survival difference between patients in different m6A modification patterns in PAAD patients of high immune score, which indicated that immune score could affect the survival possibility of PAAD patients in different m6A modification patterns. Unfortunately, subsequent analysis demonstrated that stromal score could not affect the survival probability of patients in different m6A modification patterns (Fig. 6G, P = 0.087).

Construction of the m6A gene signature

We immediately tried to construct a m6A gene signature which significantly associated with prognosis of patients with PAAD. As shown in Fig. 7A, 146m6A phenotype-related DEGs which differently expressed between m6Acluster A and m6Acluster B were firstly screened out in total, via the standards we set before (adjusted P value < 0.001 and |log2FC | ≥ 1.5). The detail information for these DEGs were shown in Table S3. Among the 146 m6A-related genes, 48 genes were significantly high expressed and 98 genes were obviously low expressed in m6Acluster A compared to m6Acluster B (Fig. 7B). Then univariate Cox regression analysis was used to identify prognosis-related gene from the 146 DEGs. 75 prognosis-related genes including 25 favorable factors and 50 risk factors were further desired. As the alluvial diagram showed, all the 25 favorable factors of prognosis belonged to the 48 down-regulated DEGs in m6Acluster B meanwhile 50 risk factors of prognosis were included in the 98 up-regulated DEGs in m6Acluster B (Fig. 7C). Furthermore, as shown in Fig. 7D, based on the 75 prognosis-related factors, unsupervised clustering algorithm was conducted and the 179 PAAD patients were further classified into two different genomic subtypes (m6A gene cluster A: 90 patients; m6A gene cluster B: 89 patients). Survival analysis indicated that PAAD patients in m6A gene cluster A subtype had worse overall survival (OS) compared with patients in m6A gene cluster B (P = 0.005, Fig. 7E). All of these outcomes indicated that m6A methylation modification acted as an indispensable role in the prognosis of patients with PAAD. Nonetheless, the above results were only in view of the patient population, which meant those might not play significant roles in m6A methylation modification pattern predication for individual patients. Therefore, based on the 75 prognosis-related factors, we further set up a consolidated scoring system for quantifying the m6A modification pattern of individual PAAD patients. By combining HR of each prognosis-related factor (Table S4) and PCA, we calculated m6A-score for each PAAD patient with the formula we mentioned above. The value of m6A-score in predicting outcome of PAAD patients was immediately explored. Based on the median value of m6A-score, patients were split into high m6A-score group (n = 89) and low m6A-score group (n = 90). Survival analysis demonstrated that PAAD patients in low m6A-score group had better OS compared with these in high m6A-score group (P < 0.001, Fig. 7F). We further compared the m6A-score level between different m6A modification patterns and genomic subtypes. As Fig. 7G showed, the m6A-score levels of patients in m6Acluster A group were higher than these of patients in m6Acluster B group (P < 0.001, Fig. 7G). Similarly, difference of m6A-score levels also existed between patients in different m6A gene related patterns (P < 0.001, Fig. 7H). Furthermore, an alluvial diagram was plotted to show the changes of m6Acluster, survival status, m6A gene cluster and m6A-score (Fig. 7I).

Role of m6A modification patterns in immunotherapeutic benefit prediction

Immunotherapy has been proved to be an effective and safe cancer treatment. Among these therapeutic methods, PD-L1 and PD-1 blockade represents a quantum leap for medical science. In IMvigor210, PAAD patients in high m6A-score group (n = 94) occupied significantly shorter survival (P < 0.001, Fig. 8A), compared with patients in low m6A-score group (n = 204). The predictive value of the m6A-score to anti-PD-L1 immunotherapy was further confirmed by using IMvigor210 (Fig. 8B-D). Patients in higher m6A-score group were more likely to benefit from anti-PD-L1 treatment (Fig. 8B-C), which was validated by Kruskal-Wallis test (P = 0.023, Fig. 8D). Next ROC analysis (based on IMvigor210) was used to evaluate the efficacy of anti-PD-L1 treatment by m6A-score, which indicated that m6A-score was a predictive biomarker to immunotherapeutic benefits (AUC: 0.782, Fig. 8E). In addition, we also tried to find some relationship between m6A-score and anti-PD-1 treatment via GSE78220 cohort. But unfortunately, m6A-score was not significantly related to OS (P = 0.830, Fig. 8F) and response to anti-PD-1 immunotherapy (Fig. 8G-I, Kruskal-Wallis test: P = 0.510). Further ROC analysis also demonstrated that m6A-score might not be an appropriate predictive tool to anti-PD-1 treatment benefits (AUC: 0.583, Fig. 8J). Talked together, the present study obviously demonstrated that m6A methylation modification patterns was associated with response to anti-PD-L1 immunotherapy. Thus, the constructed m6A modification signature might effectively give play to the prediction of response to anti-PD-L1 treatment.

More and more evidences indicated that m6A modification was an extremely important factor in inflammation, innate immunity and anti-tumor by interacting with multiple m6A regulators. Nonetheless, most of the previous researches only focused on one kind of TME cells or single m6A related gene. The global features of TME infiltration defined by the integration of various m6A regulators have not been fully appreciated. Recognizing the influence of different m6A modification patterns in the invasion of TME cells could assist better understanding the anti-tumor immune response of TME and guide more effective immunotherapy strategies. In the present study, we firstly collected 25 m6A regulators including 15 readers, 8 writers, and 2 erasers from previous studies. In general, in PAAD, the cross talk of genetic variation of m6A regulators was immediately explored, which indicated that CNV alteration of m6A regulators were prevalent. Furthermore, the whole level of mutation frequency of m6A regulators was under-represented. We also attempted to explore the expression difference of these regulators. As the result suggested, most of them including 8 readers, 4 writers and 2 erasers were significantly higher expressed in PAAD tissues compared to normal tissues. What interested us was that m6A regulators with CNV amplification led to higher expression, which demonstrated that CNV alterations was one of the features influencing the mRNA expression level of m6A regulators. GO and KEGG enrichment analyses indicated that m6A regulators played nonnegligible roles in m6A modification. These m6A regulators might work with each other suggested by correlation analysis and PPI network. We further explored the clinical relevance of m6A regulators in PAAD. Higher expression of 13 m6A regulators were significantly correlated to worse survival, which indicated that these regulators worked as oncogenic features in PAAD.

Then based on the 25 m6A regulators, we immediately identified two distinct m6A methylation modification patterns. PAAD patients in m6Acluster B group occupied poorer survival compared to patients in m6Acluster A group, which demonstrated that m6A modification patterns might play an essential role in prognosis of PAAD patients. Furthermore, we also screened out a series of prognosis-related genes, which significantly different expressed between different m6A modification patterns. What was consistent with the clustering results of the m6A modification phenotypes was that we also assigned two genomic subtypes. Similarly, PAAD patients included in m6A gene cluster B had better overall survival compared to these in m6A gene cluster A, which indicated that m6A modification was of great significance in prognosis of PAAD.

Cancer immunotherapy is an important tumor therapy potion for prevention and treatment of tumors and has attracted tremendous interests(Wesch et al., 2020). Nowadays, more and more researchers focus on screening out novel immune-related prognostic biomarkers and therapeutic target for immunotherapy. However, similar studies regarding PAAD remains scarce. We further developed the m6A-risk scoring system to quantify the m6A modification patterns of individual tumor, which could solve the individual heterogeneity of m6A modification. Further integrated analyses indicated that m6A-score might act as a potential prognostic biomarker for PAAD patients. Higher m6A-score was significantly associated with worse overall survival in PAAD patients. We then explored the role of the m6A modification patterns in anti-PD-1/L1 immunotherapy. ROC analysis demonstrated that the m6A-score signature was associated with response to anti-PD-L1 immune therapy, which might help improve the predictive strategy for anti-PD-L1 immunotherapy.

In brief, the m6A-score might act as a useful tool to comprehensively quantify the m6A gene modification pattern and predict the prognosis of a single patient, so as to further determine the tumor immunophenotype and increase the clinical application. We also demonstrated that m6A-score could be used to evaluate the clinicopathological features of patients, including the response to immunotherapy. Similarly, the m6A-score had been demonstrated to be an independent tool which could effectively and accurately predict the prognosis of patients with PAAD. Moreover, by using the m6A-score, the curative effect of adjuvant chemotherapy and the clinical response to anti-PD-L1 immune therapy for individual patient could be significantly forecasted. We can also predict the efficacy of adjuvant chemotherapy and the clinical response of patients to anti-PD-1 / PD-L1 immunotherapy by m6score. What’s more, the present research provided some fresh ideas for tumor immunotherapy. In detail, this study focused on m6A regulators or genes related to m6A phenotype, thus changing the m6A modification patterns, which might contribute to the development of new drug combination strategies or new immunotherapy in the years to come. The work in this study put forward some new insights to increase the clinical response to immune therapy for individual patient, explore various cancer immunophenotypes, and develop personalized precision immunotherapy strategy over the next few decades.

All told, the present study proved the m6a modification pattern was an independent prognostic biomarker in PAAD. And the difference of the m6A modification patterns was an indispensable feature for heterogeneity and complexity of the individual TME. The comprehensive evaluation of the modification mode of m6A in a single tumor would help us to improve our understanding of the infiltration characteristics of tumor microenvironment and guide more effective immunotherapy strategies.

Ethics approval and consent to participate

The experimental protocol was established according to the ethical guidelines of the Helsinki Declaration and was approved by the Human Ethics Committee of West China Hospital, Sichuan University. Written informed consent was obtained from individual participants or their guardian.

Consent for publication

Not applicable.

Availability of data and material

The data that support the findings of this study are openly available in The Human Protein Atlas (HPA) database (https://www.proteinatlas.org/), The Cancer Genome Atlas (TCGA) database at https://genomecancer.ucsc.edu/, Gene Expression Omnibus (GEO) database at (https://www.ncbi.nlm.nih.gov/geo/) and International Cancer Genome Consortium (ICGC) database at https://daco.icgc.org/.

Competing interests

There is no conflict of interest in the present study.

Funding

This study was supported by grants from the National Natural Science Foundation of China (No. 82070625, No. 82100650), and Technological Supports Project of Sichuan Province (No.2022YFS0257).

Authors' contributions

ML and CWC collected the literatures and drafted the initial manuscript. TFW, ML and CWC revised the manuscript and edited the language. TFW conceptualized and guaranteed the review. XYZ, CL, JYS and PW designed the figures and tables. XZ, CWC, YXZ, PW and CL formatted the references and whole manuscript. All authors contributed to the article and approved the submitted version. ML and CWC contributed to this paper equally.

Acknowledgements

None.

Aut, M.G., and Aut, T.L. (2016). forestplot: Advanced Forest Plot Using 'grid' Graphics. doi.
Batista, P.J. (2017). The RNA Modification N(6)-methyladenosine and Its Implications in Human Disease. Genomics Proteomics Bioinformatics 15, 154–163. doi: 10.1016/j.gpb.2017.03.002.
Bengtsson, A., Andersson, R., and Ansari, D. (2020). The actual 5-year survivors of pancreatic ductal adenocarcinoma based on real-world data. Sci Rep 10, 16425. doi: 10.1038/s41598-020-73525-y.
Binnewies, M., Roberts, E.W., Kersten, K., Chan, V., Fearon, D.F., Merad, M., Coussens, L.M., Gabrilovich, D.I., Ostrand-Rosenberg, S., Hedrick, C.C., Vonderheide, R.H., Pittet, M.J., Jain, R.K., Zou, W., Howcroft, T.K., Woodhouse, E.C., Weinberg, R.A., and Krummel, M.F. (2018). Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat. Med. 24, 541–550. doi: 10.1038/s41591-018-0014-x.
Cerami, E., Gao, J., Dogrusoz, U., Gross, B.E., Sumer, S.O., Aksoy, B.A., Jacobsen, A., Byrne, C.J., Heuer, M.L., and Larsson, E. (2012). The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2, 401–404. doi.
Chai, R.C., Wu, F., Wang, Q.X., Zhang, S., Zhang, K.N., Liu, Y.Q., Zhao, Z., Jiang, T., Wang, Y.Z., and Kang, C.S. (2019). m(6)A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gliomas. Aging (Albany NY) 11, 1204–1225. doi: 10.18632/aging.101829.
Chokkalla, A.K., Mehta, S.L., and Vemuganti, R. (2020). Epitranscriptomic regulation by m(6)A RNA methylation in brain development and diseases. J Cereb Blood Flow Metab 40, 2331–2349. doi: 10.1177/0271678X20960033.
Dai, D., Wang, H., Zhu, L., Jin, H., and Wang, X. (2018). N6-methyladenosine links RNA metabolism to cancer progression. Cell Death Dis. 9, 124. doi: 10.1038/s41419-017-0129-x.
Fang, H., and Declerck, Y.A. (2013). Targeting the tumor microenvironment: from understanding pathways to effective clinical trials. Cancer Res. 73, 4965-77. doi: 10.1158/0008-5472.CAN-13-0661.
Han, D., Liu, J., Chen, C., Dong, L., Liu, Y., Chang, R., Huang, X., Liu, Y., Wang, J., Dougherty, U., Bissonnette, M.B., Shen, B., Weichselbaum, R.R., Xu, M.M., and He, C. (2019). Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells. Nature 566, 270–274. doi: 10.1038/s41586-019-0916-x.
Hugo, W., Zaretsky, J.M., Sun, L., Song, C., Moreno, B.H., Hu-Lieskovan, S., Berent-Maoz, B., Pang, J., Chmielowski, B., Cherry, G., Seja, E., Lomeli, S., Kong, X., Kelley, M.C., Sosman, J.A., Johnson, D.B., Ribas, A., and Lo, R.S. (2016). Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 165, 35–44. doi: 10.1016/j.cell.2016.02.065.
Imam, H., Kim, G.W., and Siddiqui, A. (2020). Epitranscriptomic(N6-methyladenosine) Modification of Viral RNA and Virus-Host Interactions. Front Cell Infect Microbiol 10, 584283. doi: 10.3389/fcimb.2020.584283.
Liu, J., Yue, Y., Han, D., Wang, X., Fu, Y., Zhang, L., Jia, G., Yu, M., Lu, Z., Deng, X., Dai, Q., Chen, W., and He, C. (2014). A METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nat. Chem. Biol. 10, 93 − 5. doi: 10.1038/nchembio.1432.
Liu, N., Zhou, K.I., Parisien, M., Dai, Q., Diatchenko, L., and Pan, T. (2017). N6-methyladenosine alters RNA structure to regulate binding of a low-complexity protein. Nucleic Acids Res. 45, 6051–6063. doi: 10.1093/nar/gkx141.
Liu, Z.X., Li, L.M., Sun, H.L., and Liu, S.M. (2018). Link Between m6A Modification and Cancers. Front Bioeng Biotechnol 6, 89. doi: 10.3389/fbioe.2018.00089.
Marabelle, A., Tselikas, L., de Baere, T., and Houot, R. (2017). Intratumoral immunotherapy: using the tumor as the remedy. Ann. Oncol. 28, xii33-xii43. doi: 10.1093/annonc/mdx683.
Mariathasan, S., Turley, S.J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., Kadel, E.I., Koeppen, H., Astarita, J.L., Cubas, R., Jhunjhunwala, S., Banchereau, R., Yang, Y., Guan, Y., Chalouni, C., Ziai, J., Şenbabaoğlu, Y., Santoro, S., Sheinson, D., Hung, J., Giltnane, J.M., Pierce, A.A., Mesh, K., Lianoglou, S., Riegler, J., Carano, R., Eriksson, P., Höglund, M., Somarriba, L., Halligan, D.L., van der Heijden, M.S., Loriot, Y., Rosenberg, J.E., Fong, L., Mellman, I., Chen, D.S., Green, M., Derleth, C., Fine, G.D., Hegde, P.S., Bourgon, R., and Powles, T. (2018). TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554, 544–548. doi: 10.1038/nature25501.
Martini, D.J., Kline, M.R., Liu, Y., Shabto, J.M., Williams, M.A., Khan, A.I., Lewis, C., Collins, H., Akce, M., Kissick, H.T., Carthon, B.C., Shaib, W.L., Alese, O.B., Pillai, R.N., Steuer, C.E., Wu, C.S., Lawson, D.H., Kudchadkar, R.R., El-Rayes, B.F., Ramalingam, S.S., Owonikoko, T.K., Harvey, R.D., Master, V.A., and Bilen, M.A. (2019). Adiposity may predict survival in patients with advanced stage cancer treated with immunotherapy in phase 1 clinical trials. Cancer-Am. Cancer Soc.. doi: 10.1002/cncr.32576.
Mayakonda, A., Lin, D.C., Assenov, Y., Plass, C., and Koeffler, H.P. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28, 1747–1756. doi: 10.1101/gr.239244.118.
Meyer, K.D., and Jaffrey, S.R. (2014). The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol 15, 313 − 26. doi: 10.1038/nrm3785.
Newman, A.M., Liu, C.L., Green, M.R., Gentles, A.J., Feng, W., Xu, Y., Hoang, C.D., Diehn, M., and Alizadeh, A.A. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453-7. doi: 10.1038/nmeth.3337.
Nombela, P., Miguel-López, B., and Blanco, S. (2021). The role of m(6)A, m(5)C and Ψ RNA modifications in cancer: Novel therapeutic opportunities. Mol. Cancer 20, 18. doi: 10.1186/s12943-020-01263-w.
Pinello, N., Sun, S., and Wong, J.J. (2018). Aberrant expression of enzymes regulating m(6)A mRNA methylation: implication in cancer. Cancer Biol Med 15, 323–334. doi: 10.20892/j.issn.2095-3941.2018.0365.
Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., Shi, W., and Smyth, G.K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47. doi: 10.1093/nar/gkv007.
Robin, X., Turck, N., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J.C., and Muller, M. (2011). pROC: an open-source package for R and S + to analyze and compare ROC curves. BMC Bioinformatics 12, 77. doi: 10.1186/1471-2105-12-77.
Robinson, M.D., Mccarthy, D.J., and Smyth, G.K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. doi.
Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B., and Ideker, T. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498 − 504. doi: 10.1101/gr.1239303.
Siegel, R.L., Miller, K.D., and Jemal, A. (2020). Cancer statistics, 2020. CA: A Cancer Journal for Clinicians 70, 7–30. doi: 10.3322/caac.21590.
Sprooten, J., Ceusters, J., Coosemans, A., Agostinis, P., De Vleeschouwer, S., Zitvogel, L., Kroemer, G., Galluzzi, L., and Garg, A.D. (2019). Trial watch: dendritic cell vaccination for cancer immunotherapy. Oncoimmunology 8, e1638212. doi: 10.1080/2162402X.2019.1638212.
Szklarczyk, D., Gable, A.L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., Simonovic, M., Doncheva, N.T., Morris, J.H., Bork, P., Jensen, L.J., and Mering, C.V. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607-D613. doi: 10.1093/nar/gky1131.
Tang, B., Yang, Y., Kang, M., Wang, Y., Wang, Y., Bi, Y., He, S., and Shimamoto, F. (2020). m(6)A demethylase ALKBH5 inhibits pancreatic cancer tumorigenesis by decreasing WIF-1 RNA methylation and mediating Wnt signaling. Mol. Cancer 19, 3. doi: 10.1186/s12943-019-1128-6.
Therneau, T.M. (2015). survival: Survival Analysis. Technometrics 46, 111–112. doi.
Tong, J., Cao, G., Zhang, T., Sefik, E., Amezcua, V.M., Broughton, J.P., Zhu, S., Li, H., Li, B., Chen, L., Chang, H.Y., Su, B., Flavell, R.A., and Li, H.B. (2018). m(6)A mRNA methylation sustains Treg suppressive functions. Cell Res. 28, 253–256. doi: 10.1038/cr.2018.7.
Wesch, D., Kabelitz, D., and Oberg, H.H. (2020). Tumor resistance mechanisms and their consequences on γδ T cell activation. Immunol. Rev.. doi: 10.1111/imr.12925.
Wilkerson, M.D., and Hayes, D.N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572-3. doi: 10.1093/bioinformatics/btq170.
Xie, F., Xu, M., Lu, J., Mao, L., and Wang, S. (2019). The role of exosomal PD-L1 in tumor progression and immunotherapy. Mol. Cancer 18, 146. doi: 10.1186/s12943-019-1074-3.
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., Treviño, V., Shen, H., Laird, P.W., Levine, D.A., Carter, S.L., Getz, G., Stemke-Hale, K., Mills, G.B., and Verhaak, R.G. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi: 10.1038/ncomms3612.
Yu, G., Wang, L.G., Han, Y., and He, Q.Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi.
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y.L. (2020). m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol. Cancer 19, 53. doi: 10.1186/s12943-020-01170-0.

No competing interests reported.

Download PDF

Editor invited by journal
27 Aug, 2024
Submission checks completed at journal
26 Aug, 2024
First submitted to journal
20 Aug, 2024

You are reading this latest preprint version

Molecular characterization and prognostic relevance of m6A regulators in pancreatic adenocarcinoma

Status:

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

INTRODUCTION

MATERIALS AND METHODS

Collection of m6A regulators

Collection and preprocessing of pancreatic adenocarcinoma genome-related data

Function enrichment analysis and cross talks among m6A regulators

m6A regulator prognostic relevance

Consensus clustering for m6A regulators

Landscape of tumor microenvironment (TME) cell infiltration between different m6A distinct phenotypes

m6A-related gene identification

m6A gene signature generation

Relationship between m6A modification patterns and anti-PD-1/L1 immunotherapy

RESULTS

Cross talk of genetic variation of m6A regulators in PAAD

Expression difference of m6A regulators between PAAD tissues and normal tissues

Cross talks among m6A regulators

Survival landscape of m6A regulators in PAAD

m6A modification patterns defined via 25 m6A regulators

Association between m6A modification patterns and tumor microenvironment

Construction of the m6A gene signature

Role of m6A modification patterns in immunotherapeutic benefit prediction

DISCUSSION

CONCLUSIONS

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1