Multiomics integrative analysis for gene signatures and prognostic values of m⁶A regulators in pancreatic adenocarcinoma: a retrospective study in The Cancer Genome Atlas project

Mutations and CNAs of m⁶A regulatory genes in TCGA_PAAD cohorts

M⁶A-related research has grown rapidly in recent years. To include all genes that have been experimentally confirmed to yield m⁶A regulating functions as much as possible, we first searched PubMed for all literatures relating to m⁶A. After carefully archiving, 2 reviews published in 2019 that comprehensively summarized all the m⁶A regulatory genes was adopted [15, 16]. Finally, twenty-four genes were included in subsequent studies by intersecting with the genome and transcriptome-sequencing gene set on TCGA, including 9 “writers”, 13 “readers” and 2 “erasers”. Protein-protein interaction (PPI) network analysis was then applied using “STRING” database (https://string-db.org/) and CytoScape software (National Resource for Network Biology, USA) to display the interaction between these m⁶A regulators. As shown in Figure 1, three “writer” genes, METTL3, METTL14 and KIAA1429 (also known as VIRMA) seemed to be the hub genes in this network. This result brought some preliminary hints for our next research.

To analyze somatic mutations of m⁶A regulatory genes, three cohort from different countries with complete genomic sequencing data and clinical information were included, they were TCGA_PAAD cohort, ICGC_AU_PAAD cohort and ICGC_CA_PAAD cohort, respectively (ICGC: international cancer genome consortium; TCGA: from America; AU: from Australia; CA: from Canada). All included genomic sequencing data were got from the same sequencing platform (Illumina HiSeq 2000, SAN DIEGO, USA). Among all 753 patients with complete sequencing data and clinical information, mutations of 24 m⁶A regulatory genes were found in 53 (7.04%) independent samples (Figure 2, Table 1, and Supplementary Table 1). On an individual level, most of these 53 cases had only one single nucleotide mutation of a certain m⁶A regulatory gene (Table 1 and Supplementary Table 1). The only one exception was from TCGA_PAAD cohort (ID: TCGA.IB.7651.01) who suffered a huge m⁶A mutation burden, with 21 non-synonymous mutations occurred in 11 m⁶A regulatory genes. From a genetic point of view, none of the included m⁶A genes had a mutation frequency of more than 1% in TCGA and ICGC database; compared with the mutation frequency of control genes KRAS and TP53. Thus, we had enough reason to believe that mutations of m⁶A regulatory genes were not a common event in PAAD patients (Figure 2). More specifically, a total of 6 patients had 6 mutation events in WTAP, made it the gene with the highest incidence of mutations in all 753 included patients; another 4 and 3 patients had 10 mutations in IGF2BP1 and YTHDC1, respectively. IGF2BP1 and YTHDC1 became the genes with the highest number of mutation events.

Given that no copy number alteration (CNA) related data was provided in two datasets from ICGC database, TCGA_PAAD cohort became our only remaining data source for the next CNA analysis. In this cohort, 125 of these 149 cases (83.89%) had CNA events of m⁶A regulatory genes, which meant that CNAs occurred very frequently in PAAD patients. The “writer” gene METTL16 had the highest frequency of CNA events (57.05%, 85/149), followed by WTAP and ALKBH5, which were “writer” and “eraser” genes, respectively (Figure 3A).

Then we evaluated CNA patterns in TCGA_PAAD cohort. The total number of CNA events was 1084 and most events led to the copy number loss of m⁶A genes (60.15%, 652/1084) (Figure 3B, 3C), which consisted of homozygous and hemizygous deletion, indicating the loss of two alleles, or one (or less than one) allele for a certain gene from genome in a diploid organism. In detail, hemizygous deletion of METTL16 (82 cases) was the most frequent alteration in all kinds of CNAs of m⁶A regulatory genes (Table 2), followed by the hemizygous deletions of WTAP (79 cases) and ALKBH5 (69 cases). The second most frequent CNA event was copy number gain (35.98%, 390/1084), which happened to VIRMA and HNRNPA2B1 in 40 patients, tied for the first place, followed by the “reader” gene IGF2BP3 (38 cases). To verify the reliability of this cohort, we also detected KRAS (all 149 cases with mutations and 39 cases with CNAs) and TP53 (94 cases with mutations and 72 with CNAs) as control, which were also widely recognized by the published literatures.

Alterations of m⁶A regulatory genes were associated with clinicopathological characteristics and survival of PAAD patients

The relationship between alterations (CNAs and/or mutations) of m⁶A regulatory genes and the clinicopathological characteristics of PAAD patients were analyzed. Results showed that alterations of m⁶A genes were significantly associated with the patients’ gender, tumor grade and the maximum diameter of resected tumor. More concretely, male patients had a higher frequency of mutations or CNAs than female patients; patients with larger tumor volumes and higher histopathological grade also tended to show increased alterations of m⁶A genes (Table 3).

The above results made us interested in the impact of these alterations on patients’ prognosis. Thus, we analyzed the effects of overall alterations of m⁶A regulatory genes and some certain high-frequency CNAs on patients’ overall survival (OS) and disease-free survival (DFS). However, results did not seem to be ideal. Individuals with or without mutations or CNAs of m⁶A regulatory genes didn’t have any significant correlation with OS and DFS (Figure 4A, 4B). Furthermore, separate analysis of the included 24 genes revealed that most single CNA events didn’t affect OS and DFS (Supplementary Figure 1, 2). The only exception is IGF2BP2, a m⁶A “reader” gene. Patients with copy number gain of IGF2BP2 showed significantly poorer OS and DFS than other patients without this alteration (Figure 4C, 4D).

Then we divided the 24 m⁶A regulatory genes into 3 classes (writers, readers and erasers) based on their functions in m⁶A modification and detected whether the loss or gain of overall copy number for a class of m⁶A genes could influence patients’ OS and DFS. We found that copy number loss mainly happened to “writers” and “erasers”, yet “readers” mainly gained copy number. Therefore, patients were divided by whether “writer loss”, “eraser loss” or “reader gain” happened. It turned out that patients with “reader gain” combined with “writer loss”, and those with “reader gain” and “eraser loss” had worse OS (Figure 4E, 4F). Besides, no significant difference was found about DFS of these groups (Supplementary Figure 3A–3D). These results highlighted the potential importance of m⁶A reader genes on the prognosis of PAAD patients.

CNAs of m⁶A regulatory genes were associated with mRNA and protein expression levels

We wondered whether CNAs of m⁶A regulatory genes at the genomic level were related to their mRNA expression status. The results of correlation analysis showed that mRNA expression level was significantly associated with most of the CNA patterns in PAAD samples. In 21 of 24 detected genes, copy number gains and high-level amplifications were related to higher mRNA expression; while homozygous and hemizygous deletions resulted in a declined level of mRNA (Figure 5A–5R and Supplementary Figure 4A–4C). No significant relation was found between the CNAs of FTO, RBM15 and ZC3H13 and their mRNA expression level (Supplementary Figure 4D–4F).

Then we detected whether mRNAs of these 24 m⁶A regulatory genes are differentially expressed between PAAD patients and normal individuals. RNA-sequencing data from TCGA_PAAD, ICGC_AU_PAAD and ICGC_CA_PAAD dataset, and normal pancreas transcriptome data from Genotype-Tissue Expression (GTEx) database were enrolled, including 279 tumor samples and 328 normal pancreas samples. Unfortunately, none of the m⁶A regulatory genes met the conventional criterion of differential expression (log2FoldChange > 1.5, P < 0.05) (Supplementary Table 2). There were only 4 genes which showed statistically significant differential expression (P <0.05), among which the most differentially expressed one was Insulin-like growth factor 2 mRNA binding protein 2 (IGF2BP2) (log2FoldChange = 0.341, P = 0.048). Heatmap showed that the expression of m⁶A -regulated genes cannot fully distinguish patients from normal individuals through clustering analysis, either (Supplementary Figure 5). However, some interesting hints drew our attention. First, the mRNA expression level of m⁶A regulatory genes were consistent with their copy number status; those genes that mostly showed copy number loss, like METTL3, ALKBH5 and YTHDF family, all showed a tendency of lower expression compared with normal tissues; On the contrary, IGF2BP family, which mainly suffered copy number gain in PAAD patients, showed a tendency of mRNA over-expression in PAAD compared with normal pancreas. Besides, it was not difficult to notice from the heatmap that the most capable genes to distinguish tumors from normal tissues in Cluster Analysis were from IGF2BP protein family as well.

Based on the above results, we had a strong interest in the expression status and mechanism of IGF2BP family in PAAD pathogenesis. IGF2BP family contained three members. As mentioned in our previous results, IGF2BP1 had the highest mutation frequency, while IGF2BP3 showed the highest CNA frequency; however, only IGF2BP2 was associated with patient prognosis, and most likely there was a significant increase in its mRNA expression. Thus, IGF2BP2 was chosen as the candidate gene for further experimental detection. Immunohistochemistry staining for IGF2BP2 expression was firstly performed in 3 pairs of PAAD tissues and their adjacent normal tissues. We found that IGF2BP2 protein was highly expressed in all three pairs of PAAD tissues than that of normal (Figure 6A and Supplementary Figure 6). Then, 20 pairs of PAAD clinical specimens with their adjacent normal pancreatic tissue were extracted to detect the mRNA expression level of IGF2BP2. It turned out that IGF2BP2 was significantly overexpressed in PAAD cells from clinical samples than normal pancreatic cells (Figure 6B). We further utilized 5 cell lines (four PAAD cell lines: Panc-1, CFPAC-1, BXPC-3, SW1990; one normal pancreatic ductal epithelial cell line: HPDE) to verify IGF2BP2 expression at mRNA and protein level in vitro. qRT-PCR and Western Blotting results further confirmed the presence of up-regulated mRNA and protein expression of IGF2BP2 in PAAD (Figure 6C, 6D).

Apart from the above research, survival analysis was also performed by dividing patients of TCGA_PAAD cohort into two groups according to the median mRNA expression level of every single gene. However, none of the mRNA expression of the 24 genes was found to significantly influence the OS and DFS of PAAD patients (Supplementary Figure 7).

Predictive potential of IGF2BP2 for the prognosis of PAAD patients and downstream Enrichment Analysis for IGF2BP2 overexpression

The results so far made us very interested in the m⁶A “reader” gene IGF2BP2. To preliminarily explore the clinical and biological function of IGF2BP2, univariate and multivariate Univariate and multivariate Cox regression analysis was applied to further confirm the relationship between copy number gain of IGF2BP2 and prognosis of PAAD patients. Results suggested that copy number gain of IGF2BP2 alone was an independent risk factor for worse OS (for multivariate Cox regression, HR = 2.392, 95%CI: 1.392-4.112, p<0.001) and DFS (for multivariate Cox regression, HR = 2.400, 95%CI: 1.236-4.659, p=0.010) of PAAD patients (Table 4). In summary, we believed that IGF2BP2 was an important indicator leading to higher severity and poorer prognosis of PAAD patients.

We then determined to explore the role of genetic alteration of IGF2BP2 in the pathogenesis of PAAD. Gene Set Enrichment Analysis (GSEA) suggested that high expression of IGF2BP2 was related to some critical biological processes including 1) regulation of centrosome cycle; 2) telomere organization; 3) β-catenin destruction complex disassembly; 4) cytokinesis (Figure 7A–7D and Table 5). All these 4 processes were highly conserved and related to the malignant behavior of tumors, including cell proliferation, cell immortalization and tumor immunization. To preliminarily explore how many m⁶A sites IGF2BP2 could bind to as a m⁶A “reader” protein, RMBase V2.0 (http://rna.sysu.edu.cn/rmbase/index.php) based on epitranscriptome sequencing data was then explored. Result of the data mining showed that IGF2BP2 protein was capable of binding to a total of 60,935 m⁶A sites on tens of thousands of mRNAs. These results might provide some new ideas for subsequent researchers who would like to study the molecular mechanism of IGF2BP2 in pancreatic cancer.

Ethics statement

All the CNAs, mutations, mRNA expression data and clinical information were retrieved from TCGA database by cBioportal platform [30] through R package “cgdsr” and UCSC-Xena program [31] which are open sourced to the public under certain guidelines. The ICGC database was further applied to prove TCGA results. Thus, it is confirmed that all written informed consent was achieved. For patients’ tissue array from the third Xiangya hospital, written informed consent was obtained from all patients and the study was approved by the ethics committee of the third Xiangya Hospital.

Data processing

We identified 3 datasets (TCGA_PAAD, ICGC_AU_PAAD and ICGC_CA_PAAD) with complete clinical information, of which all three datasets included somatic mutation and RNA-seq raw data, while CNA data was only available in TCGA_PAAD cohort. For somatic mutation research, raw file in “MAF” format were got from R package “TCGAmutations” and ICGC data portal. The R packages “Maftools” and “Genvisr” were then applied to analyze mutation frequency and draw “waterfall plot”. For CNA research, the loss and gain of copy number changes was identified using segmentation analysis and GISTIC algorithm. When clinical significance of the status of CNA and/or mutation was discussed, this PAAD cohort was usually divided into two subgroups, “with CNAs and/or mutation of all m⁶A regulatory genes (or a certain m⁶A regulatory gene)” and “without CNA and mutation”. The mRNA expression data were calculated from RNA-Seq V2 RSEM release and log2 fold change was applied before analyzing the relationship between CNA and mRNA expression. The R package “edgeR” [32] was utilized to analyze mRNA differential expression.

Statistical analysis

All statistical data and figures were analyzed by R (Lucent Technologies, New Jersey, USA) or GraphPad Prism 8.0 (GraphPad Software, La Jolla, CA, USA) to ensure the aesthetics and editability of the figures and tables. The association between CNAs of m⁶A regulatory genes and clinicopathological characteristics were analyzed using chi-square test or Mann-Whitney U test. Kaplan-Meier curve and log-rank test were used to evaluate the prognosis value of all kinds of CNAs that the frequency is greater than 15% of total sample number. Cox proportional hazard regression model was performed using R package “survival” [33] and “survminer” [34]. All statistical results with a p-value <0.05 were considered significant.

Clinical samples

A total of 20 pancreatic cancer tissues and adjacent normal tissues were collected from Department of Hepatobiliary and Pancreatic Surgery, The Third Xiangya Hospital. All tissues were frozen immediately in liquid nitrogen after surgical excision and stored at −80°C. Clinicopathological information was retrieved from the hospital database. The ethic approval was mentioned above.

Quantitative real-time PCR (qRT-PCR), Western Blotting (WB) and Immunohistochemistry (IHC)

For qRT-PCR, Total RNAs were isolated from PC tissues and cell lines by TRIzol reagent (Invitrogen) and reverse transcription was performed with a Prime Script RT reagent kit (Takara Biotechnology, China). Qrt-PCR reactions were carried out on Applied Biosystems 7500 Real-time PCR Systems (Thermo Fisher Scientific). GAPDH was used as control for IGF2BP2 mRNA expression. The relative expression level of the target RNA was calculated by 2^−ΔΔCt. In our study, all involved primers were produced by Sangon Biotech (Shanghai, China).

For WB analysis, Cells were lysed in 1x RIPA buffer with 1mM PMSF and 1x complete protease inhibitor cocktail (Roche, Basel, Switzerland). Western blottings were conducted under standard procedures. Briefly, Protein lysates were separated by SDS-PAGE, and electrophoretically transferred to PVDF membrane from MilliporeTM with nominal pore size of 0.22 um. After blocking, the membrane was incubated at 4°C overnight with anti-IGF2BP2 (Proteintech) or β-Actin (Santa Cruz Biotechnology) diluted in non-fat dry-milk (Santa Cruz Biotechnology), followed by HRP-labeled secondary (Santa Cruz Biotechnology) antibody incubation. Chemiluminescence signal was developed by ECL Plus Western Blotting Detection Reagents (GE Healthcare Life Sciences, Munich, Germany).

For IHC test, IHC staining of IGF2BP2 was conducted on tissue microarray consisting of 12 PDAC tissues and 12 paired normal tissues. IHC staining was carried out using Histostain-Plus kit according to manufacturer’s protocol. Briefly, antigen retrieval was conducted by heating the sections in boiling sodium citrate buffer for 20 mins. After 3% hydrogen peroxide and BSA blocking, the tissues were incubated with anti-IGF2BP2 (Proteintech, Chicago, IL) at 4°C overnight. After washing, the tissues were incubated with HRP-conjugated secondary antibody (Santa Cruz Biotechnology, Santa Cruz, CA). IHC signal was developed by DAB substrate, and counter-stained by haematoxylin. Three random fields at 200x magnification were captured per sections for evaluation. Scoring of the IHC was based on the percentage of positive cells and staining intensity under a light microscope. Four categories (0, 1+, 2+ and 3+) were denoted as negative, weak, moderate and strong respectively, and final score was calculated by averaging the score of three fields.

Downstream pathway and target analysis for IGF2BP2

Gene set enrichment analysis (GSEA) provided by the JAVA program (Version 4.0.3) with MSigDB v6.1 was applied to explore the downstream biological processes affected by differential expression of IGF2BP2. Patients in the cohort was divided into two groups according to the median expression level of IGF2BP2. Finally, 25880 genes were enrolled into the GSEA process. Hallmark gene set “c5.bp.v7.0.symbols.gmt” was used in this study [35]. Gene sets which got highest Enrichment Score (ES) with normalized p-value <0.05, and the false discovery rate (FDR) <0.25 were considered significantly enriched.