Dysregulation of N6-methyladenosine regulators predicts poor patient survival in mantle cell lymphoma.

N6-methyladenosine (m6A) is the most abundant eukaryote mRNA modification, modulated by regulators known as epigenetic writers, erasers and readers, which are known to serve crucial roles in mRNA metabolism. However, the role of m6A during B-cell development and B-cell tumorigenesis remains poorly understood. By analyzing the gene expression profile of 123 mantle cell lymphoma cases from the Gene Expression Omnibus database, the present study demonstrated that one-half of the m6A regulators were able to predict patient survival in mantle cell lymphoma, notably the m6A.index. The expression levels of the m6A regulators were regarded as good classifiers in mantle cell lymphoma. The m6A.index-low mantle cell lymphoma type exhibited a poor patient survival and lower mRNA levels from the total transcriptome. The m6A regulators may be associated with the cell division and the RNA metabolic pathways, which may result in poor survival of patients with mantle cell lymphoma.

Introduction

N6-methyladenosine (m6A) is the most abundant eukaryote messenger RNA modification, modulated by regulators known as writers, erasers and readers (1). It is known to serve crucial roles in mRNA metabolism and is primarily involved in RNA stability, mRNA splicing and protein translation (2–9).

The proteins that are involved in m6A modifications consist of ‘writers’, ‘erasers’ and ‘readers’. The m6A writers are the m6A methyltransferase enzymes, a 70-kDa complex consisting of 3 components: Methyltransferase like (METTL) 3, METTL14 and WT1 associated protein (WTAP), which methylate the adenosine motif (A) at the N6 position (10–12). The m6A erasers comprise the m6A demethyltransferase enzymes, including alpha-ketoglutarate-dependent dioxygenase FTO (FTO) and RNA demethylase ALKBH5 (ALKBH5), which are the first and second m6A demethyltransferase enzymes (6,13–15). The m6A readers consist of effectors (m6A RNA binding protein) that decode the m6A methylation code. The YTH domain family [YTH domain-containing family protein (YTHDF) 1, YTHDF2 and YTHDF3] and ELAV-like protein 1 (ELAVL1) are known as m6A readers (2,3,16,17). The m6A modification is associated with cancer progression. The m6A demethylase ALKBH5 sustains forkhead box protein M1 expression and cell proliferation and maintains the tumorigenesis of glioblastoma stem-like cells (18). m6A RNA methylation regulates the tumorigenesis of glioblastoma and the self-renewal of glioblastoma stem cells (19).

METTL3 is an m6A writer associated with the formation of undifferentiated myeloid cells in acute myeloid leukemia (AML), as well as with chemo- and radio-resistance of pancreatic cancer cells (20,21). In addition, METTL3 serves an important role in the growth, survival and invasion of human lung cancer cells (22). Protein virilizer homolog (KIAA1429) was defined as a writer of m6A in 2014 (23). KIAA1429 is a unique type of m6A writer: i) Mammalian KIAA1429 (202 kDa) is the largest known component within the m6A methyltransferase complex; ii) among all the components examined, the depletion of KIAA1429 resulted in the largest decrease in m6A levels [KIAA1429 depletion led to a 60% decrease in m6A levels, while METTL3, METTL14, and WTAP depletion resulted in 30, 40, and 50% decreases in m6A levels, respectively (12)]; and iii) biochemical studies have indicated that KIAA1429 recruits METTL3/METTL14/WTAP, the catalytic core components (24–26). Therefore, KIAA1429 may serve as a scaffold molecule of the methyltransferase complex, and serve a unique role that is different from those of the catalytic core components METTL3, METTL14 and WTAP.

FTO is a member of the m6A eraser family of proteins and has been demonstrated to promote AML and lung squamous cell carcinoma (LUSC) progression (27,28).

The expression levels of the m6A protein readers, including YTHDF1 and YTHDF2, were identified to be markedly associated with malignancy and poor prognosis of hepatocellular carcinoma (29,30). YTHDF2 and YTHDF3 were defined as readers of m6A in 2012 (16). YTHDF2 and YTHDF3 regulate messenger RNA stability (2,31). YTHDF1 was defined as a reader of m6A in 2014. YTHDF1 facilitates messenger RNA translation initiation, while ELAVL1 inhibits protein translation (3,32). It appears that YTHDF2 and YTHDF3 are more likely to regulate mRNA stability, while YTHDF1 and ELAVL1 are more likely to regulate mRNA translation.

Mantle cell lymphoma is a type of non-Hodgkin B cell lymphoma with a median age of diagnosis at 60 years (33,34). Mantle cell lymphoma has an aggressive phenotype and a rapid rate of progression, with a short median survival of 5–7 years (35). The investigation of the molecular mechanisms that contribute to the aggressive phenotype of mantle cell lymphoma may provide novel treatment strategies. Recently, previous studies demonstrated that m6A mRNA methylation controls T cell homeostasis and modulates the characteristics of hematopoietic stem and progenitor cells (36,37). However, the role of m6A during B-cell development and B-cell tumorigenesis remains poorly understood (38). The present study demonstrated that the imbalanced expression levels of m6A regulators may be used for the prediction of poor survival in patients with mantle cell lymphoma.

Materials and methods

Data source

The Affymetrix Human Genome U133 Plus 2.0 array of 123 mantle cell lymphoma samples was retrieved from the NCBI Gene Expression Omnibus (GEO) database (GSE93291 dataset) (39,40). The detailed patient demographic data and disease characteristics for this dataset were published previously (39). The Affymetrix Human Genome U133 Plus 2.0 Array that contained 64 mantle cell lymphoma samples was retrieved from the NCBI GEO database (GSE21452) (41). GSE21452 was the first phase of the GSE93291 dataset.

Gene expression analysis

The probeset measures of all the arrays were calculated by robust multiarray averaging. The relative RNA expression values were log-transformed (log2). The data were analyzed with an unpaired Student's t-test and are presented as the mean ± standard error of the mean (SEM). P<0.05 was considered to indicate a statistically significant difference. Only genes with a fold change (log2) >1 or <-1 were defined as differentially expressed genes.

Definition of m6A.index for survival prediction

A comprehensive m6A.index was defined to predict the survival in patients with mantle cell lymphoma. The m6A.index was calculated using a previously described method (1), as follows:

where m6A.indexj represents the index of m6A of the jth sample for survival prediction. Fj represents the product of favorable gene expression of the jth sample. A total of 7 out of 10 m6A genes exhibited a hazard ratio <1 and were defined as ‘favorable genes’, which were favorable for the survival of mantle cell lymphoma. Hj represents the product of ‘harmful gene’ expression of the jth sample. A total of 3 out of 10 m6A genes (YTHDF1, KIAA1429 and ELAVL1) exhibited a hazard ratio >1 and were defined as ‘harmful genes’, which were harmful for the survival of mantle cell lymphoma.

The median of the m6A.index value from a cohort of mantle cell lymphoma, for example 123 patients with mantle cell lymphoma, was defined as the cut-off value for the m6A.index-low and m6A.index-high groups consisting of 62 and 61 samples, respectively. The correlation of the m6A.index with the gene expression levels of the marker of proliferation Ki-67 (Ki-67) was assessed using the Spearman's correlation test and the pairwise colored scatterplot was drawn based on a Kernel Density Estimation using the LSD package (version 4; cran.r-project.org/web/packages/LSD/index.html) in R.

Gene Ontology (GO) analysis

The Database for Annotation, Visualization and Integrated Discovery tool with default parameters was used for GO analysis (42). All enriched GO terms identified in the present study were manually prepared so that only selected, non-redundant GO terms in the ‘Biological Process’ category were identified.

Statistical analysis

The R software v3.1.3 (ggplot2 package) was used for the statistical analysis. Kaplan-Meier curves were used to plot survival curves of YTHDF3, METTL14, ALKBH5, ELAVL1 and KIAA1429 genes. For the YTHDF3, METTL14, ALKBH5 and ELAVL1 genes, the median of the gene expression value from a cohort of mantle cell lymphoma (123 patients with mantle cell lymphoma) was defined as the cut off value for the low and high expression groups. For KIAA1429 genes, the maximally selected rank statistics algorithm (survminer package) was used to define the low and high expression groups. Survival analysis of those genes was performed using the log-rank test. The heatmap depicted the cosine correlation similarity between 10 m6A regulators. Unpaired Student's t-tests were used for the statistical analysis of quantitative variables. The data are expressed as the mean ± SEM in scatter plots. P<0.05 was considered to indicate a statistically significant difference.

Results

Specific m6A regulators predict patient survival in mantle cell lymphoma

To investigate the association between the m6A regulators [METTL3, METTL14, WTAP and KIAA1429 (writers); FTO and ALKBH5 (erasers); and YTHDF1, YTHDF3, YTHDF2 and ELAVL1 (readers)] and the survival of the patients with mantle cell lymphoma, the expression profiles of 123 mantle cell lymphoma samples from the GSE21452 dataset were analyzed. A total of 5 out of 10 m6A regulators revealed expression levels that were significantly associated with the survival of patients with mantle cell lymphoma (P<0.05, log-rank test). The 10 m6A regulators were classified according to the hazard ratio values. A total of 7 out of 10 m6A genes had a hazard ratio value <1 and were defined as ‘favorable genes’. YTHDF3 exhibited the highest statistical significance of all the ‘favorable genes’ with a hazard ratio of 0.51 [95% confidence interval (CI), 0.27–0.96]. A total of 3 out of 10 m6A genes (YTHDF1, KIAA1429 and ELAVL1) exhibited a hazard ratio value of >1 and were defined as ‘harmful genes’ (Table I; Fig. 1). ELAVL1 exhibited the highest significance of all the ‘harmful genes’ with a hazard ratio of 2.46 (95% CI, 1.35–4.48). Kaplan-Meier curves were conducted to assess the overall survival of the 123 patients with mantle cell lymphoma in association with 5 m6A regulators. The comparisons were performed using a log-rank test (Fig. 2). In addition, the 123 patients with mantle cell lymphoma were divided into KIAA1429-high (n=24) and KIAA1429-low (n=99) groups using maximally selected rank statistics algorithm (survminer package), and the survival curves of the two groups were compared using a log-rank test (Fig. S1; P<0.0001). The m6A regulators of both ‘favorable’ and ‘harmful’ genes were deemed significant predictors for mantle cell lymphoma patient survival.

Table I.

Survival analysis of 10 m6A regulators.

Gene	Hazard ratio	95% confidence interval	P-value	m6A role
ELAVL1	2.46	1.35–4.48	0.0030	Reader
KIAA1429	2.09	1.21–3.62	0.0086	Writer
YTHDF3	0.51	0.27–0.96	0.0356	Reader
METTL14	0.68	0.47–0.99	0.0435	Writer
ALKBH5	0.70	0.49–1.00	0.0493	Eraser
WTAP	0.81	0.65–1.02	0.0683	Writer
METTL3	0.59	0.32–1.07	0.0840	Writer
FTO	0.61	0.31–1.20	0.1519	Eraser
YTHDF2	0.77	0.51–1.16	0.2113	Reader
YTHDF1	1.29	0.51–3.30	0.5881	Reader

ELAVL1, ELAV-like protein 1; KIAA1429, Protein virilizer homolog; YTHDF, YTH domain-containing family protein; METTL, methyltransferase like; ALKBH5, RNA demethylase ALKBH5; WTAP, WT1 associated protein; FTO, alpha-ketoglutarate-dependent dioxygenase FTO.

Figure 1.

Forest plots of 10 N6-methyladenosine regulators associated with survival. The black lines indicate lower and upper 95% confidence of the hazard ratios.

Figure 2.

Kaplan-Meier curves measuring associations between overall survival and 4 N6-methyladenosine regulators in 123 patients with mantle cell lymphoma. YTHDF3 (P=0.0356), METTL14 (P=0.0435), ALKBH5 (P=0.0493) and ELAVL1 (P=0.003). The log-rank test was used to compare the Kaplan-Meier curves. YTHDF3, YTH domain-containing family protein 3; METTL14, methyltransferase like 14; ALKBH5, RNA demethylase ALKBH5; ELAVL1, ELAV-like protein 1.

Expression patterns of m6A regulators predict poor or favorable patient survival in patients with mantle cell lymphoma

An unsupervised clustering of the expression levels of the 10 m6A regulators was conducted in 123 patients with mantle cell lymphoma, and the cosine correlation similarity was presented as a heatmap (Fig. 3). Notably, it was identified that the 10 m6A regulators were classified into two groups; WTAP, ALKBH5 and METTL14 in one group, and the other 7 genes in a second group. ELAVL1 and KIAA1429, two of the most significantly ‘harmful’ genes, were clustered together in one group. The other 3 most significantly ‘favorable’ genes, YTHDF3, METTL3 and FTO, were also clustered together in one group. The clustering of ‘harmful’ and ‘favorable’ genes together in single groups suggests that the genes that are the most closely associated with survival cluster together. Furthermore, it was observed that the mantle cell lymphoma cancer type could be classified into two groups by fuzzy clustering (Fig. 4A). The ratio of the expression levels of ‘harmful genes’ to ‘favorable genes’ was estimated, which was termed as the m6A.index. The m6A.index was highly associated with the survival of patients with mantle cell lymphoma (Fig. 4B; P<0.05). The hazard ratio of the m6A.index was 0.39 (95% CI, 0.24–0.65). The m6A.index reflected the imbalanced expression between ‘harmful genes’ and ‘favorable genes’ of the m6A regulators. The m6A.index-low group was associated with a poor patient survival in mantle cell lymphoma, while the m6A.index-high group was associated with a favorable patient survival in mantle cell lymphoma.

Figure 3.

Unsupervised clustering of the expression levels of 10 N6-methyladenosine regulators in 123 patients with mantle cell lymphoma. The cosine correlation similarity was depicted by the heatmap. ELAVL1, ELAV-like protein 1; KIAA1429, Protein virilizer homolog; YTHDF, YTH domain-containing family protein; METTL, methyltransferase like; ALKBH5, RNA demethylase ALKBH5; WTAP, WT1 associated protein; FTO, alpha-ketoglutarate-dependent dioxygenase FTO.

Figure 4.

Groups of 10 m6A regulators were used as a classifier in the 123 patients with mantle cell lymphoma. (A) The fuzzy clustering of 123 patients with mantle cell lymphoma as determined by the expression levels of the 10 m6A regulators. PC1 and PC2 were the first and second components, respectively. Each point indicated a mantle cell lymphoma sample. Colors 1 (red) and 2 (green) denote two clusters. (B) Kaplan-Meier curves for assessment of the association between overall survival of 123 patients mantle cell lymphoma with m6A.index (P<0.001). The log-rank test was used to compare the Kaplan-Meier curves. m6A, N6-methyladenosine.

The correlation of the m6A.index with the gene expression levels of marker of proliferation Ki-67 (Ki-67) in 123 mantle cell lymphoma samples was analyzed (Fig. S2; Cor =−0.52; P<0.001; Spearman's correlation test). The m6A.index exhibited a negative correlation with Ki-67 and DNA polymerase genes (DNA polymerase alpha 1, catalytic subunit, DNA polymerase alpha 2, accessory subunit, DNA polymerase delta 1, catalytic subunit, DNA polymerase delta 2, accessory subunit and DNA polymerase theta). Although, the mutational status of tumor protein 53 (TP53), ATM serine/threonine kinase (ATM) and MYC proto-oncogene, BHLH transcription factor (MYC) could not be assessed using the data from the present study, the survival analysis of 123 patients with mantle cell lymphoma was analyzed in association with Ki-67, ATM, MYC and TP53 gene expression (Fig. S2). Increased expression levels of MKI67 and MYC corresponded with poorer survival in patients with mantle cell lymphoma (Fig. S3; P=2.9×10−11 and P=4.5×10−4, respectively). The high expression levels of ATM and TP53 demonstrated the trend for predicting favorable survival in patients with mantle cell lymphoma (Fig. S3; P=7.5×10−2 and P=1.1×10−1, respectively). ATM and Ki-67 were differentially expressed between the m6A.index-low and m6A.index-high groups (Fig. S3; P<0.05), while TP53 and MYC were not differentially expressed (P>0.05). Furthermore, the correlations between the m6A.index with the expression of 17 proliferation-associated genes in 123 mantle cell lymphoma samples were calculated. A total of 16 of 17 proliferation-associated genes were highly correlated with the m6A.index (Fig. S4).

Association between the upregulation of gene expression and the m6A.index-high group in mantle cell lymphoma

The m6A.index-high and the m6A.index-low groups were identified to be two different classes of mantle cell lymphoma. Therefore, the expression profiles of the m6A.index-high and the m6A.index-low groups in mantle cell lymphoma were compared (Fig. 5A). A total of 280 upregulated genes and 54 downregulated genes were identified between the m6A.index-high and the m6A.index-low groups in mantle cell lymphoma (Fig. 5B; P<0.05). The m6A.index-high mantle cell lymphoma group exhibited an increased number of upregulated genes compared with the m6A.index-low group, which suggested that the m6A.index-high mantle cell lymphoma type had a different RNA metabolism process from the m6A.index-low mantle cell lymphoma type. The cumulative distribution of the expression of RNA molecules corresponding to different genes with regard to the m6A.index-high and m6A.index-low mantle cell lymphoma types additionally demonstrated that the m6A.index-high mantle cell lymphoma type exhibited high RNA levels compared with the total transcript profile in the GSE93291 dataset (Fig. 5C; P<0.001). This result was also validated in the secondary GSE21452 dataset (n=64).

Figure 5.

Different expression of the genes in the m6A.index-high and the m6A.index-low groups of patients with mantle cell lymphoma. (A) The heatmap indicates the different expression of the genes in the m6A.index-high and m6A.index-low groups of patients with mantle cell lymphoma. Red, high expression; green, low expression; white, moderate expression. Only the top 12 upregulated and downregulated genes were noted. The two bar plots in the left of the heatmap refer to the fold-change (log2, left; green and red) difference and the P-value (-log10, right; blue), respectively. (B) The total number of the upregulated (280 genes) and the downregulated genes (54 genes) between the m6A.index-high and m6A.index-low mantle cell lymphoma types. (C) Cumulative distribution of RNA levels (log2-fold) of the differentially expressed genes between the m6A.index-high (red) and the m6A.index-low (green) mantle cell lymphoma groups. The left plot represents the GSE93291 dataset (n=123), and the right plot represents the GSE21452 dataset (n=64) (D) The different expression levels of the MALAT1 and the WHSC1 genes between the m6A.index-high and the m6A.index-low mantle cell lymphoma groups. The left plot corresponds to the GSE93291 dataset (n=123), and the right plot corresponds to the GSE21452 dataset (n= 64). A two sided unpaired Student's t-test was used. *P<0.05 and ****P<0.0001. m6A, N6-methyladenosine.

Cell division and RNA metabolism pathways are significantly enriched pathways of m6A in mantle cell lymphoma

A characteristic difference between the m6A.index-high and m6A.index-low mantle cell lymphoma types was evident. WHSC1, which encodes for a histone-lysine N-methyltransferase, was the top downregulated gene in the m6A.index-high group compared with the m6A.index-low group (Fig. 5D; P=4.81×10−6). Metastasis associated in lung adenocarcinoma transcript 1 (MALAT1) was the top upregulated gene in the m6A.index-high group compared with the m6A.index-low group (Fig. 5D; P=1.98×10−6). The differential expression of WHSC1 and MALAT1 was also validated in the GSE21452 dataset (n=64). A pathway analysis of the differentially expressed genes between the m6A.index-high and the m6A.index-low mantle cell lymphoma types was conducted. The cell division and RNA metabolism pathways were demonstrated to be the most significantly enriched pathways based on the differential expression of specific genes (Fig. 6A). The RB transcriptional corepressor 1, cell division cycle 25A and kinesin family member C1 genes were included in the cell division pathways group, and were demonstrated to be differentially expressed, out of a total of 9 genes (Fig. 6B). Therefore, the m6A regulators may regulate the cell division pathways, contributing to poor patient survival in mantle cell lymphoma.

Figure 6.

Significantly enriched pathways of differentially expressed genes between the m6A.index-high and the m6A.index-low groups of patients with mantle cell lymphoma. (A) The bar plot indicates the significantly enriched cell division (left) and the RNA metabolism (right) pathways. (B) The differentially expressed genes identified between the m6A.index-high and the m6A.index-low mantle cell lymphoma groups in the cell division pathway. The left plot corresponds to the GSE93291 dataset (n=123), while the right plot corresponds to the GSE21452 dataset (n=64). A two-sided unpaired Student's t-test was used. P<0.05 for all histograms. m6A, N6-methyladenosine.

Discussion

Recent studies have indicated that m6A mRNA methylation controls T-cell homeostasis and modulates hematopoietic stem and progenitor cell differentiation (36,37). However, the role of m6A during B-cell development and B-cell tumorigenesis remains poorly understood (38). The present study demonstrated that the imbalanced expression of m6A regulators may be used for the prediction of patient survival in mantle cell lymphoma. A previous study indicated that m6A regulator participates in the innate immunity via the RNA helicase probable ATP-dependent RNA helicase DDX46. Therefore, the m6A gene and the m6A regulators are involved in B-cell lymphoma development, T-cell homeostasis and innate immunity.

Mantle cell lymphoma is an aggressive type of B cell lymphoma with a short median patient survival time. The identification of novel biomarkers for the prediction of patient survival in mantle cell lymphoma is considered a challenging task (43). The mantle cell lymphoma international prognostic index (MIPI) score is currently the most common prognostic model for mantle cell lymphoma in clinical practice (44). MIPI includes the age, Eastern cooperative oncology group performance status, leukocyte count and lactate dehydrogenase activity (45). However, these models lack a component that incorporates gene expression analysis. The present study demonstrated that the expression levels of 5 m6A regulators were significantly associated with the survival of patients with mantle cell lymphoma. In addition, a comprehensive m6A.index was constructed to predict the survival of mantle cell lymphoma. The m6A.index was better compared with each individual m6A regulator for survival prediction with a hazard ratio of 0.39 (95% CI, 0.24–0.65).

It was surprising that half (5 of 10) of the m6A regulators were able to predict patient survival in mantle cell lymphoma (P<0.05). It appeared that the m6A regulators were commonly associated with mantle cell lymphoma patient survival. In addition, the m6A.index was a better survival predictor compared with the single m6A regulator, which suggested that the imbalanced expression of m6A regulators may predict poor patient survival in mantle cell lymphoma. Furthermore, the data from the present study provided evidence to support a marked association between m6A expression and mantle cell lymphoma incidence: i) The expression of m6A regulators was a good classifier in mantle cell lymphoma; ii) the samples of mantle cell lymphoma were divided into two groups according to the m6A.index (m6A.index-high and m6A.index-low groups), with the m6A.index-high group exhibiting high RNA levels compared with the total transcript profile; and iii) the differentially expressed genes were associated with the cell division and RNA metabolism pathways in the m6A.index-high group that may result in poor patient survival.

The m6A.index was correlated with enhanced proliferation. Ki-67 is clinically important for risk stratification and clinical management of mantle cell lymphoma (45–47). The m6A.index demonstrated a highly negative correlation with gene expressions of Ki-67 and DNA polymerases. A total of 17 proliferation-associated genes were defined as a proliferation ‘signature’ and associated with overall survival in mantle cell lymphoma (39). The present study further correlated the m6A.index with the expression of 17 proliferation-associated genes in 123 mantle cell lymphoma samples; 16 of these 17 proliferation-associated genes were highly correlated with m6A.index, which suggested that the m6A.index was associated with proliferation in mantle cell lymphoma.

Mammalian KIAA1429 is the largest known component within the m6A methyltransferase complex and serves as a scaffold of the methyltransferase complex, while METTL3/METTL14/WTAP serve as the catalytic core components (24). In the results from the present study, the majority of the m6A writer molecules, including METTL3, METTL14 and WTAP, were classified as ‘favorable genes’, while the remaining m6A writer (KIAA1429) was categorized into the group of ‘harmful genes’. Therefore, within the group of 4 m6A writers, KIAA1429 appeared to be the most significant in terms of its biological function and clinical implications.

In summary, the results from the present study demonstrated that the expression levels of m6A regulators were associated with the survival of patients with mantle cell lymphoma and may serve as potential biomarkers for prognosis.