Expression of N6-Methyladenosine (m6A) Regulators Correlates With Immune Microenvironment Infiltration and Predicts Prognosis in Diffuse Large Cell Lymphoma (DLBCL)

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

Download PDF

Primary research

Expression of N6-Methyladenosine (m6A) Regulators Correlates With Immune Microenvironment Infiltration and Predicts Prognosis in Diffuse Large Cell Lymphoma (DLBCL)

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

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: N⁶-methyladenosine (m⁶A) and immune microenvironment infiltration have been widely reported to play important roles in various cancers. However, in diffuse large cell lymphoma (DLBCL), the clinical significance of m⁶A regulators and their relationship with immune microenvironment infiltration have not been illuminated.

Methods: The expression of m⁶A regulators in DLBCL was investigated using The Cancer Genome Atlas and Genotype-Tissue Expression databases. The capacity of m6A regulators in dividing molecular clusters of DLBCL was determined using Consensus Clustering algorithm and validated via principal component analysis. The clinical traits and prognosis difference in m⁶A-sorted clusters were revealed. The m⁶A prognostic signature was established and validated based on Gene Expression Omnibus dataset using univariate Cox and LASSO regression analysis. The immune cell infiltration and the expression of immune checkpoint genes in m⁶A low/high-risk DLBCL were studied. Gene set enrichment analysis (GSEA) was adopted to unveil the underlying molecular mechanism in m⁶A low/high-risk DLBCL.

Results: Differentially expressed m⁶A regulators were able to molecularly discriminate DLBCL as two clusters based on consensus clustering and principal component analysis. A six m⁶A regulators-based risk prediction signature was established and validated as an independent predictor, which separated patients into m⁶A low- and high-risk groups. High-risk m⁶A indicates worse survival, for which the predictive AUC at 1-, 2-, and 5-year achieved 0.605, 0.640, and 0.652, respectively. Immune cell infiltration analysis revealed the B cells naïve and T cells gamma delta were the top increased and decreased immune cells in high-risk m⁶A patients. Up-regulated (PDCD1 and KIR3DL1) and down-regulated (TIGIT, DO1, and BTLA) immune checkpoint genes in high-risk m⁶A patients were identified. GSEA analysis unveiled high-risk m⁶A related to tumor proliferation associated process, while low-risk m⁶A related to defense response associated process.

Conclusions: This study provided a comprehensive analysis for the clinical significance of m⁶A regulators and their association with immune microenvironment infiltration. An m⁶A regulators-based risk signature may be applied for the risk stratification of DLBCL patients, thus may facilitate the clinical management of DLBCL. Immune microenvironment was found to be closely related to m⁶A risk, which may be part of the mechanism of m⁶A regulators in DLBCL.

Cancer Biology

m6A regulator

immune microenvironment

prognostic signature

DLBCL

Latest cancer statistics reported that there were estimated 81,560 new cases and 20,720 deaths of non-Hodgkin lymphomas (NHL) in the United State, 2021. NHL is becoming the 7^th and 6^th in morbidity for male and female, respectively. The mortality ranked 9^th for both male and female, seriously threatening human health ¹. Diffuse large cell lymphoma (DLBCL), accounting for 25-30% of all NHL, is a greatly heterogeneous B cell lymphoid neoplasm with substantial variations in genome and genetic alterations, which cause diverse clinical phenotype and response to therapy. Over the past two decades, the combination of rituximab, cyclophosphamide, doxorubicin, vincristine, and prednisone (R-CHOP) has been established as standard first-line therapy for DLBCL patients based on the results of several phase III clinical trials ^2-4. About 50-70% of DLBCL patients can be clinically cured using the R-CHOP regimen, while the remaining patients turn out to be either refractory or relapsed using this approach. Worse still, only about 10% of the refractory or relapsed DLBCL patients have the fortune to be cured using intensive salvage immunochemotherapy followed by autologous stem cell transplantation, whereas the 90% rest of the patients suffer from dismal survival ^5,6. The urgent need for effective therapeutic strategies is unmet and requires unremitting efforts.

With the advent of the high-throughput genome sequencing technique, the construction of genetic landscape of DLBCL has become a reality. It is realizable to decipher which patient is likely or unlikely to benefit from the therapy on genomic level. For example, according to gene expression profile, the DLBCL could be divided into two molecularly distinct subtypes: germinal center cell (GCB)-like and activated B cell (ABC)-like DLBCL, which was a milestone for interpreting why DLBCL patients have a different response to the same R-CHOP therapy ⁷. Moreover, increasing researches have progressively unveiled the pivotal driver genes and pathways in DLBCL, such as TP53, MYC, BCL6, BCL2, MYD88, BCR pathway, NF-κB pathway, PI3K-AKT-mTOR pathway, and JAK-STAT pathway, etc., which help better understand the biological and pathological processes in DLBCL ⁸. What even more inspiring is growing novel targeted agents such as BTK inhibitor Ibrutinib ^9-11, BCL2 inhibitor Venetoclax ^12,13, PI3K inhibitor CUDC-907 ¹⁴, and AKT inhibitor MK-2206 ¹⁵, etc. have shown promising therapeutic potential in relapsed/refractory DLBCL. Therefore, providing a deeper insight into the mechanism in DLBCL is no doubt the essential foundation of novel targeted therapy development, which apparently is far from being satisfied yet. As a consequence, the research on exploring the DLBCL field from molecular mechanism to clinical application is urgently needed.

Recently, the research on epitranscriptome in cancers has progressed in leaps and bounds owing to the rapid development of high-through sequencing such as chromatin immunoprecipitation sequencing, methylated RNA immunoprecipitation sequencing, and Assay for Transposase-Accessible Chromatin using sequencing. N⁶-Methyladenosine (m⁶A) is one of the common and abundant post-transcriptional modifications in mRNAs, which exerts a pivotal regulatory role in tumorigenesis and cancer progression. The process of m⁶A is mainly achieved by three types of m⁶A regulators, that is methyltransferases (“writers”), demethylases (“erasers”), and binding proteins (“readers”). m⁶A modifications are assembled by “writers”, removed by “erasers”, and deciphered by “readers”. Growing studies have revealed important roles of the m⁶A regulators in cancers. For example, methyltransferases METTL3 has been reported as an oncogene involving in the cell proliferation, differentiation, invasion, migration, and apoptosis of various tumors including acute myeloid leukemia, breast cancer, liver cancer, gastric cancer, bladder cancer, prostate cancer, lung cancer, and pancreatic cancer ¹⁶. Demethylase ALKBH5 also exerts effects in tumor proliferation, invasion, and metastasis in lung cancer, gastric cancer, pancreatic cancer, colon cancer, glioblastoma, osteosarcoma, and ovarian cancer as an oncogene or tumor suppressor ¹⁷. And the “readers” YTH domain-containing proteins, including YTHDF1-3 and YTHDC1-2, have also been proved to be closely in connection with poor prognosis of colorectal cancer, hepatocellular carcinoma, breast cancer, and ovarian cancer ¹⁸. The combination of m⁶A regulators for constructing prognostic models has recently been a novel signature for predicting the outcome of tumor patients. For instance, m⁶A high-risk has been clarified as an unfavorable indicator in clear cell renal cell carcinoma ¹⁹, pancreatic cancer ²⁰, non-small cell lung cancer ²¹, head and neck squamous cell carcinoma ²², gastric cancer ²³, hepatocellular carcinoma ²⁴, colorectal carcinoma ²⁵, ovarian cancer ²⁶, and cervical squamous cell carcinoma ²⁷, which indicates worse survival of patients. However, in DLBCL, no study has reported the combined prognostic value of m⁶A regulators. As a result, it is imperative to investigate the prognostic value of combined m6A regulators in DLBCL, which is promising for novel targeted therapy development and clinical management of DLBCL patients.

In the current study, we will be the first to comprehensively investigate the prognostic value and construct a prognostic signature of m⁶A regulators in DLBCL. More importantly, we will provide new perspectives for the relationship between m⁶A and immune cell infiltration in DLBCL (Figure 1). We hope the work achieved here will throw light on the interactions between epitranscriptome and immune microenvironment, as well as their clinical potential in DLBCL.

Collection of m⁶A regulators and DLBCL data

Based on currently published literature ^28-32, a total of 22 common m6A regulators, which included 8 “Writers” (METTL3, METTL14, METTL5, WTAP, ZC3H13, RBM15, RBM15B, KIAA1429), 2 “Erasers” (ALKBH5, FTO), and 12 “Readers” (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, EIF3A, FMR1) were collected. The 22 m⁶A regulators were used for the subsequent analysis.

The Cancer Genome Atlas (TCGA), Genotype-Tissue Expression (GTEx), and Gene Expression Omnibus (GEO) databases were used for obtaining the transcriptome expression and clinical information of DLBCL samples and normal controls. The RNA-seq data of 48 DLBCL from TCGA and 337 normal controls from GTEx were downloaded, removed batch effect, merged, and normalized with log2(FPKM+1). An mRNA microarray GSE10846, which contains 410 DLBCL samples and corresponding clinical parameters such as age, gender, stage, extranodal site, ECOG performance status, molecular subtype, LDH ratio, overall survival, and survival status was included in the current study. After removing the patients with an overall survival under 90 days for a better analysis quality, a total of 380 DLBCL samples from GSE10846 were eventually included. The mRNAs expression in GSE10846 were extracted and normalized with log2 for subsequent analysis.

Differential expression and correlation analysis for m⁶A regulators

To identify significantly expressed m⁶A regulators in DLBCL, the Wilcoxon signed-rank test was adopted and achieved in R version 4.02. A fold-change of 2 or 0.5 and a p-value below 0.05 were utilized for screening statistically up-regulated or down-regulated m⁶A regulators in DLBCL. A violin plot and a heatmap plot achieved with vioplot and pheatmap packages in R version 4.02 were used to visualize the differentially expressed m6A regulators.

In order to understand the correlations among the m⁶A regulators in DLBCL, a Spearman correlation analysis was conducted for the differentially expressed m⁶A regulators in DLBCL samples from the dataset GSE10846. A corrplot R package was utilized for correlation analysis and visualization. A p-value under 0.05 was regarded as statistically significant.

Exploration of m⁶A-sorted DLBCL clusters and their clinical significance

We wondered whether the differentially expressed m⁶A regulators can contribute to distinguishing different molecular subtypes of DLBCL. Herein, an unsupervised class discovery approach named Consensus Clustering was applied via a ConsensusClusterPlus R package ³³. The uncovered DLBCL clusters by m⁶A regulators was validated using a Principal Component Analysis (PCA) ³⁴. The correlation between the m⁶A-classed clusters and the clinical parameters of DLBCL was investigated. The Chi-squared test was applied when estimating the relationship between different m⁶A-sorted clusters and clinical characteristics such as age, gender, stage, extranodal site, ECOG performance status, LDH ratio, and molecular subtype. Meanwhile, the Kaplan-Meier survival analysis was utilized to compare the survival of DLBCL patients in different clusters. A statistical p value<0.05 was used to evaluate the difference.

Construction and validation of m⁶A prognostic signature in DLBCL

The 380 DLBCL samples in GSE10846 were randomly divided into two groups, which is the training set and testing set, equally. In the training set, 190 DLBCL patients were included for the univariate Cox regression analysis to select prognostic m⁶A regulators. A p value<0.05 was used as a cut-off to screen statistically significant prognostic m⁶A regulators for further analysis. The Least Absolute Shrinkage and Selection Operator (Lasso) regression algorithm, a popular method for regression analysis with high‐dimensional features ³⁵, was further adopted for variable determination and regularization so as to raise the prediction accuracy and interpretability of the prognostic model. The m⁶A risk stratification based on a risk score, which was computed as follows: risk score= . In the formula, the X_i represents the expression of each m⁶A regulator in the model, while Coef_irepresents the coefficient of each m⁶A regulator in the Lasso regression analysis. In the testing set and the whole combined samples set, the same formula was utilized to calculate risk scores and construct validation models. The DLBCL patients were classified as m⁶A low-risk and high-risk groups according to a median value of the m6A risk score. A survminer package and a survivalROC package in R were used to depict the K-M survival curve and Receiver Operating Characteristic curve (ROC) curve, respectively. Since most of the DLBCL adverse events occur in the first two years after diagnosis ³⁶, the area under the curve (AUC) of the time-dependent ROC, which is calculated for evaluating the predictive capacity of the m6A signature, was assessed at time points of 1, 2, and 5 years. Moreover, the statistical difference of DLBCL clinical traits (age, gender, stage, extranodal site, ECOG performance status, LDH ratio, and molecular subtype) in m⁶A low/high-risk DLBCL was evaluated using Chi-squared test.

To determine whether m⁶A prognostic signature is an independent prognostic factor in DLBCL, univariate and multivariate Cox regression analyses were performed for m⁶A signature along with age, gender, stage, extranodal site, ECOG performance status, LDH ratio, and molecular subtype in GSE10846. Hazard ratio (HR), 95% confidence interval (CI), and p-value were calculated.

Exploration of m⁶A signature and immune cell infiltration

We wonder whether m⁶A risk stratification is correlated to tumor immune microenvironment. As a result, we adopted CIBERSORT, a powerful analytical tool to estimate the abundances of immune cells using gene expression data ^37,38, to identify 22 different kinds of immune cells in the DLBCL samples from GSE10846. Then, we compared the immune cell infiltration difference between m⁶A low-risk and high-risk groups using Wilcoxon signed-rank test. The Spearman correlation analysis was also used to construe the correlation among the differentially distributed immune cells. Moreover, the expression of nine previously reported immune checkpoint genes (CTLA4, LAG3, TIGIT, HAVCR2, PDCD1, IDO1, VISTA, KIR3DL1, BTLA) ³⁹ in m⁶A low-risk and high-risk groups were also investigated to uncover the potentially possible relationship between m⁶A and immunotherapy. A p-value under 0.05 was deemed to be statistically significant.

Implementation of gene set enrichment analysis (GSEA)

GSEA is a knowledge-based computational approach for interpreting a priori defined set of genome-wide expression profiles, which helps determine statistically significant concordant differences of biological function between two phenotypes ^40,41. As a result, for the purpose of illuminating the underlying molecular pathways related to m6A risk, GSEA was performed to select the most significant gene ontology terms (GO), which includes GO biological process, GO cellular component, and GO molecular function, in m⁶A low/high-risk groups. GSEA version 4.10 software and Molecular Signatures Database v7.2 were applied for the current analysis.

Differentially expressed m6A regulators in DLBCL

The expression of the 22 m⁶A regulators in 48 DLBCL and 337 normal controls were compared. As shown in Figure 2A, the heatmap displayed the expression distribution of the 22 m⁶A regulators. And the violin plot in Figure 2B visualized the expression difference of each m⁶A regulator. In summary, 19 out of 22 regulators were differentially expressed, among which, 10 were up-regulated (EIF3A, HNRNPA2B1, YTHDC2, ZC3H13, ALKBH5, RBM15, METTL5, YTHDF1, RBM15B, YTHDF2) and 9 were down-regulated (METTL3, FTO, FMR1, YTHDC1, HNRNPC, WTAP, YTHDF3, IGF2BP2, KIAA1429) in DLBCL. In Figure 2C, the Spearman correlation analysis was achieved for the 19 differentially expressed m⁶A regulators. We can observe that the correlations among the m⁶A regulators were complicated. For example, the “Writers” can be positively correlated such as METTL3 and RMB15B (coef=0.46), or negatively correlated such as METTL3 and KIAA1429 (coef=-0.32). The two “Erases” (ALKBH5 and FTO) were positively associated with each other (coef=0.58). And the “Readers” can also be positively related like YTHDF3 and YTHDC1 (coef=0.44), or negatively related such as YTHDF3 and HNRNPA2B1 (coef=-0.51). Similarly, a “Writer” can be positively related to a “Reader” such as METTL3 and HNRNPA2B1 (coef=0.68), or negatively related to another “Reader” such as METTL3 and YTHDF3 (coef=-0.61).

m⁶A-sorted DLBCL clusters and clinical implication

As we can observe in Figure 3A-3C and Figure S1, when clustering the samples using a specified cluster count (k) from 2 to 9, only k=2 showed the best performance to distinguish the DLBCL samples by the differentially expressed m⁶A regulators. In Figure 3D-3E, the heatmap showed that the samples in cluster 1 were highly correlated, so were those in cluster 2. There was barely any correlation significant between cluster 1 and cluster 2, which indicated that m⁶A regulators were capable to differentiate the DLBCL samples into molecularly distinguishable subtypes. The finding was validated using a PCA analysis, which was displayed in Figure 3F. The samples in cluster 1 or cluster 2 were highly clustered, which further proved that m⁶A-sorted DLBCL clusters are reliable. The clinical significance of the two m⁶A-sorted DLBCL clusters was also investigated and demonstrated in Figure 4. Figure 4A displayed the expression heatmap of the differentially expressed m⁶A regulators in m⁶A sorted clusters and their correlation with DLBCL clinical parameters. The gender, age, ECOG performance status, number of extranodal sites, and survival status showed statistical differences in cluster 1 and cluster 2 (p<0.05). And in Figure 4B, the survival of DLBCL patients in cluster 1 and cluster 2 also showed statistical differences. The overall survival rate of DLBCL patients in cluster 2 outperformed than those in cluster 1 (p<0.001).

m⁶A signature and clinical implication

Univariate Cox regression algorithm was carried out for the 19 differentially expressed m⁶A regulators. Six out of 19 m⁶A regulators (ALKBH5, FMR1, HNRNPC, RBM15B, YTHDC2, and YTHDF1) show statistically prognostic value in DLBCL patients (Figure 5A). Through further Lasso regression, all the six prognostic m⁶A regulators were determined for prognostic signature construction (Figure 5B-5C). The calculated risk score curve and survival information of the patients in the training set were displayed in Figure 5D. The K-M survival curve suggested the patients with m⁶A low-risk have higher overall survival than those in m⁶A high-risk group in the training set (Figure 5E). ROC curve revealed the AUC at 1-, 2-, and 5-year was 0.674, 0.699, and 0.691 which suggested a certain predictive capacity of the m⁶A signature in the training set.

The m⁶A prognostic signature in the training set was validated using a testing set and a combined data of the training set and testing set as a whole. In Figure 6A, the risk score curve and survival information of patients in the testing set was visualized. Although there was no statistical significance when conducting a K-M survival analysis for the samples in the testing set (p=0.081), it still showed an obvious trend that patients with m⁶A low-risk have better survival (Figure 6B-6C). The AUC of the 1-, 2-, and 5-year ROC curve achieved 0.552, 0.584, and 0.611 in the testing samples. Figure 6D displayed the calculated risk score and survival information of the combined samples of the training set and testing set. A statistically significant difference was observed in the K-M survival curve analysis for the combined samples (p<0.001). Patients with m⁶A low-risk showed better survival than those in m⁶A high-risk group (Figure 6E). Figure 6F displayed the predictive efficacy of the m⁶A signature, which yielded an AUC for 1-, 2-, and 5-year of 0.605, 0.640, and 0.652, respectively.

Whether m⁶A signature is an independent factor in DLBCL was investigated. Displaying in Figure 7A, the univariate Cox regression revealed that m⁶A signature (HR=1.994, p<0.001), age (HR=1.768, p=0.001), molecular subtype (HR=2.207, p<0.001), ECOG performance status (HR=2.660, p<0.001), stage (HR=1.890, p<0.001), and LDH ratio (HR=2.641, p<0.001) were all statistically significant prognostic predictors. And multivariate Cox regression analysis displayed in Figure 7B further verified that m⁶A signature (HR=1.994, p<0.001), age (HR=1.574, p=0.017), molecular subtype (HR=1.888, p=0.003), ECOG performance status (HR=1.981, p=0.001), stage (HR=1.719, p=0.007), and LDH ratio (HR=1.877, p=0.003) were all independent prognostic predictors in DLBCL.

Estimation of immune cell infiltration

The immune cell infiltration in m⁶A low/high-risk groups was explored. The bar plot in Figure 8A visualized the proportion of the 22 immune cells in each DLBCL sample on the whole. And Figure 8B also displayed the 22 immune cells’ distribution in m⁶A low/high-risk groups as a heatmap. The immune cell infiltration difference in m6A low-risk vs. high-risk groups was compared by statistics. As shown in Figure 9A, ten out of 22 immune cells showed a statistical difference between m⁶A low-risk and high-risk groups. The B cell naïve, NK cells resting, NK cells activated, and Mast cells activated were significantly up-regulated in m⁶A high-risk group. The Plasma cells, T cells CD4 memory activated, T cells gamma delta, Dendritic cells resting, Mast cells resting, and Eosinophils were significantly down-regulated in m⁶A high-risk group. The correlations among the 10 differentially distributed immune cells were depicted in Figure 9B. We can notice that T cells CD4 memory activated was positively correlated to T cells gamma delta (coef=0.41), while negatively correlated to B cell naïve (coef=-0.35). The T cells gamma delta was negatively related to B cell naïve (coef=-0.44) and NK cells resting (coef=-0.49).

Immune checkpoint is a critical mechanism of tumor immune escape, which makes it a prospective target for cancer immunotherapy. Therefore, the expression of nine well-studied immune checkpoint genes in low/high-risk groups was analyzed and displayed in Figure 9C. We found that PDCD1 (P=0.045) and KIR3DL1 (P=0.012) were significantly increased in m⁶A high-risk group, while TIGIT (p<0.001), IDO1 (p=0.016), and BTLA (p=0.008) was significantly decreased. Regards LAG3, CTLA4, VISTA, and HAVCR2, no statistical significance was observed.

Enriched pathways in m6A low/high-risk DLBCL

To reveal the potential biological function difference between the m⁶A low-risk and high-risk groups, GSEA analysis was carried out. The top 10 significantly enriched terms in the m⁶A low-risk group and high-risk group were selected for displaying, respectively. As we can see in Figure 10, the top 10 significantly enriched terms in m⁶A high-risk group included Telomerase holoenzyme complex, Mitotic G2-M transition checkpoint, Blastocyst formation, Regulation of transcription by RNA polymerase III, Saga type complex, Base excision repair, ATPase complex, DNA replication checkpoint, Histone H4 acetylation, INO80 type complex. And GDP binding, Double stranded RNA binding, Positive regulation of endothelial cell apoptotic process, Nitric oxide synthase biosynthetic process, Barbed end actin filament capping, Positive regulation of viral genome replication, Regulation of defense response to virus, Positive regulation of interferon Alpha production, Ubiquitin dependent ERAD pathway, and Endoplasmic reticulum tubular network organization were the top 10 significantly enriched items in m⁶A low-risk group (Figure 11).

Emerging evidence substantiated that m⁶A modification exerts an indispensable regulatory effect in neuronal disorders, osteoporosis, metabolic disease, viral infection, and various cancers ^42,43. Targeting m⁶A regulators for cancer therapy has been a closely focused field by scholars. For instance, Huang et al. recently developed two promising FTO inhibitors, named FB23 and FB23-2 using structure-based rational design. Encouragingly, FB23-2 was found to dramatically inhibit proliferation and enhance the differentiation/apoptosis of acute myeloid leukemia (AML) cell line in vitro. More importantly, in vivo experiment, FB23-2 could significantly suppress the progression of AML primary cells in xeno-transplanted mice ⁴⁴. In recent years, the prominent roles of m⁶A modification in cancer immunotherapy has been gaining more and more attention. Han et al. recently discovered that m⁶A-binding protein YTHDF1 can control anti-tumor immunity through recognizing m⁶A-marked transcripts encoding lysosomal proteases to increase their translation in dendritic cells. Specifically, the deficiency of YTHDF1 elevated the cross-presentation of tumor antigens and antigen-specific CD8+ T cell antitumor response. More remarkably, loss of YTHDF1 enhanced the therapeutic efficacy of PD-L1 checkpoint blockade in vivo ⁴⁵. Despite the upsurging wave of m⁶A research in cancer, few research has reported the pathological role of m⁶A regulators and their clinical significance in DLBCL. Cheng et al. once reported that down-regulated methyltransferase METTL3 functionally inhibited the DLBCL cell proliferation through reducing the m⁶A methylation and total mRNA level of pigment epithelium-derived factor ⁴⁶. Another study revealed that up-regulated piRNA-30473 was associated with aggressive phenotype and poor prognosis of DLBCL patients by virtue of m⁶A dependent regulatory manners. Mechanistically, piRNA-30473 exerts its oncogenic effect via increasing the expression of methylase WTAP and its critical target gene HK2, thus enhance the global m⁶A level ⁴⁷. Nevertheless, the above studies only focused on a single m⁶A regulator, lacking a comprehensive understanding of the potential clinical value of different m⁶A regulators as a whole. Given that, our current study primarily focused on evaluating the prognostic value of 22 known m⁶A regulators, constructing a prognosis predictive model, and investigating the potential molecular mechanism of m⁶A risk in DLBCL.

The majority of m⁶A regulators showed differential expression in DLBCL and their interactions could be complicated. The same m⁶A regulator could be positively correlated to a regulator while negatively correlated to another regulator. This indicated the m⁶A regulation network in DLBCL could be sophisticated, which deserves further analysis. More interestingly, we found that based on differentially expressed m⁶A regulators, the DLBCL patients could be divided into two clusters. More importantly, the gender, age, extranodal sites, ECOG performance status, and overall survival of the patients in the two clusters show statistically significant differences. The patients in cluster 1 had worse survival. As we all known, DLBCL can be mainly divided into two subtypes namely ABC-DLBCL and GCB-DLBCL according to gene expression profile. ABC-DLBCL is associated with more malignant biological properties and worse clinical outcome than GCB-DLBCL. Population-based studies showed evidence that the 5-year overall survival rate of ABC-DLBCL patients is 35%, while it is 60% for GBC-DLBCL patients ⁴⁸. This suggests more positive therapy and follow-up measures should be taken when managing ABC-DLBCL patients. Similarly, different m⁶A-sorted DLBCL clusters also have distinct overall survival, which is worthy of clinical attention. Since the patients in cluster 1 and ABC-DLBCL subtype both yield worse survival, their relationship is worth further investigation. Why ABC-DLBCL is more refractory was deemed to be associated with constitutive activation of the NF-κB and BCR signaling pathways ⁴⁹. Whether m⁶A regulators are implicated in the NF-κB and BCR signaling pathways, resulting in corresponding clinical phenotype features also bears thinking about.

Another highlight achieved in the current study is we constructed an m⁶A-based prognostic model for predicting the overall survival of DLBCL patients. Patients with m⁶A high-risk was verified to yield poor outcomes. m⁶A risk scoring model was validated as an independent prognostic predictor, which is the same as age, molecular subtype, ECOG performance status, and stage. Scholars have previously constructed prognostic signatures for DLBCL using other predictors. For example, Zhou et al constructed a prognostic immunoscore model using immune cell infiltration, yielding an AUC of 0.562 in DLBCL ⁵⁰. And Hu et al. constructed an integrated prognosis model of a pharmacogenomic gene signature for DLBCL with a predictive AUC of 0.67 ⁵¹. Moreover, Zhang et al. utilized a combined five types of alternative splicing events to construct prognostic predictors for DLBCL patients, which showed an AUC of 0.564 ⁵². In the current study, the 1-year, 2-year, and 5-year predictive AUC of m⁶A signature for the 380 DLBCL patients were 0.605, 0.640, 0.652, respectively. Compared with previously published prognostic models, the m⁶A signature showed passable predictive performance, which is qualified for potential clinical application.

The concept of tumor microenvironment (TME) has been come up with for years and its development has never been kept from moving with the times. It’s now clear that TME consists of tumor cells, immune cells, stromal cells, endothelial cells, and cancer-associated fibroblasts ⁵³. Although infiltrated by immune cells, the cancer cells can somehow escape from immune supervision and destruction through multiple tricky mechanisms. The tumor immune privilege mechanisms mainly include reduced expression of cancer antigens and major histocompatibility complex class I, elevated expression of immune checkpoints, as well as increased recruitment of immunosuppressive cells, such as T regulatory cells, tumor-associated macrophages, and myeloid-derived suppressor cells, etc. ⁵⁴. Immune blockades have been developing vigorously over the years, among which, PD-1 is living proof. PD-1 inhibitor, mainly putting brakes on unrestricted cytotoxic T effector function, was first approved by U.S. Food and Drug Administration for the treatment of unresectable/metastatic melanoma cancer and non-small-cell lung cancer second-line alternative supported by National Comprehensive Cancer Network guideline ⁵⁵. However, in relapsed/refractory DLBCL, PD-1 blockade therapy has been disappointing, yielding an objective response of merely 36% ⁵⁶. As a result, making immune blockades a success in treating DLBCL is still a far way, while taking a deeper insight into the immune cell infiltration is an essential step. As a result, our current study provided novel findings on the relationship between m⁶A and immune cell infiltration, which has not been reported previously. We discovered a bunch of up-regulated (B cell naïve, NK cells resting, NK cells activated, and Mast cells activated) and down-regulated (Plasma cells, T cells CD4 memory activated, T cells gamma delta, Dendritic cells resting, Mast cells resting, and Eosinophils) immune cells in m⁶A high-risk patients. And immune checkpoint related genes could be up-regulated (PDCD1 and KIR3DL1) or down-regulated (TIGIT, IDO1, and BTLA) in m⁶A high-risk group. These indicated the complicated regulatory relationship between m⁶A and immune microenvironment. How m⁶A exactly influences immune cell infiltration and the expression of immune checkpoint, thereby impacts the outcome of DLBCL patients, remains unclear, which still requires further experimental studies.

Whether the underlying biological function and molecular pathways are different between m⁶A low-risk and high-risk DLBCL is a noteworthy aspect. Therefore, we explored the possible molecular mechanism in m⁶A low/high-risk DLBCL. We discovered that in m⁶A high-risk group, the enriched biological function mainly included cell cycle, DNA replication, transcription, post-transcriptional modification, and DNA repair relative pathways, which are manifestations of tumor malignancy features. And in m⁶A low-risk group, the enriched biological functions included several interesting items such as Positive regulation of endothelial cell apoptotic process, Nitric oxide synthase biosynthetic process, Regulation of defense response to virus, and Positive regulation of interferon Alpha production, which mostly correlate to defense process. In detail, we know that anti-angiogenesis is an important perspective for cancer therapy, in which inducing the apoptosis of vascular endothelial cells is a pivotal process ^57-59. And Nitric oxide has been found to be involved in immune response, in which its important synthase NOS2 could regulate macrophages, T cells, B cells, and myeloid derived suppressor cells. Interferon Alpha has been approving for the treatment of more than fourteen types of cancer, including hairy cell leukemia, melanoma, and renal cell carcinoma as an immune-based oncologic drug for years ^60,61. Collectively, these might help explain why patients in m⁶A high-risk group suffered from worse survival than those in m⁶A low-risk group on the molecular mechanism level. Still, experimental experiments are worthy and necessary to be carried out for further validation.

Despite of the highlights of m⁶A-based prognostic signature and their correlation with tumor immune microenvironment infiltration, several limitations should be acknowledged in this study. First, even though we constructed and validated the m⁶A risk signature through a cross-validation of 380 DLBCL cohorts, further validation with large-sample external datasets is imperative. Second, although we have investigated the correlation between m⁶A risk and immune microenvironment infiltration, their relationship to immunotherapy response remained unknown due to lack of corresponding clinical information, which required further studies. Finally, the current study design and results were based on bioinformatics analysis. Future experimental studies are required to validate the results and elaborate the exact molecular regulatory mechanism between m⁶A regulators and tumor immune microenvironment infiltration.

This study comprehensively investigated the clinical significance of multiple m⁶A regulators and established a novel m⁶A risk scoring signature for predicting the survival of DLBCL patients for the first time, which achieved satisfactory predictive performance. More importantly, we unveiled several highlights of immune cell infiltration, immune checkpoint genes and underlying molecular enriched pathways in m⁶A low/high-risk DLBCL pioneeringly, which shed light on the potential regulatory relationship between m⁶A and immune microenvironment infiltration. Nevertheless, the obtained results especially the exact mechanism on how m⁶A affects immune cell infiltration and immune checkpoint genes still need further experimental research. Anyhow, the findings in this study provided whole-new perspectives for the m⁶A epitranscriptome and immunomics in DLBCL, which are promising for the development of individualized and comprehensive management of DLBCL patients.

m⁶A: N⁶-methyladenosine; DLBCL: Diffuse large cell lymphoma; NHL: non-Hodgkin lymphomas; R-CHOP: Rituximab, Cyclophosphamide, Doxorubicin, Vincristine, and Prednisone; TCGA: The Cancer Genome Atlas (TCGA); GTEx: Genotype-Tissue Expression; GEO: Gene Expression Omnibus; PCA: Principal Component Analysis; LASSO: The Least Absolute Shrinkage and Selection ROC: Operator; Receiver Operating Characteristic curve; AUC: The area under the curve; HR: Hazard ratio; CI: Confidence interval; GSEA: Gene Set Enrichment Analysis; AML: acute myeloid leukemia; TME: Tumor microenvironment; GCB-DLBCL: Germinal center cell-like DLBCL; ABC-DLBCL: activated B cell-like DLBCL.

Acknowledgments

Not applicable

Authors’ contributions

ZGP designed and revised the manuscript. QYX and ZDH contributed to the data download and preprocessing, while MWL and HYH analyzed the data. ZCX re-checked all the data and wrote the entire manuscript. All authors have approved the final version and publication.

Funding

This study was funded by the International Scientific Exchange Foundation of China Project (Z2020LGX001), the China Anti-cancer Association Project (CORP-117), and the Guangxi Medical and Health Appropriate Technology Development and Application Project (S2020037). The funders had no participation in the article design, data analysis, decision for publication, or preparation of the manuscript.

Availability of data and materials

The clinical information, RNA-seq, and microarray expression data were obtained from the GEO (https://www.ncbi.nlm.nih.gov/geo/), GTEx (https://gtexportal.org/home/), and TCGA (https://portal.gdc.cancer.gov/) databases, which are all publicly available.

Ethics approval and consent to participate

Informed consent for all patients was obtained through the GEO, GTEx, and TCGA projects.

Consent for publication

Not applicable

Competing interests

The authors declare no conflict of interest involved.

Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA: a cancer journal for clinicians. 2021; 71(1): 7-33.
Feugier P, Van Hoof A, Sebban C, et al. Long-term results of the R-CHOP study in the treatment of elderly patients with diffuse large B-cell lymphoma: a study by the Groupe d'Etude des Lymphomes de l'Adulte. J Clin Oncol. 2005; 23(18): 4117-26.
Pfreundschuh M, Schubert J, Ziepert M, et al. Six versus eight cycles of bi-weekly CHOP-14 with or without rituximab in elderly patients with aggressive CD20+ B-cell lymphomas: a randomised controlled trial (RICOVER-60). Lancet Oncol. 2008; 9(2): 105-16.
Pfreundschuh M, Trümper L, Osterborg A, et al. CHOP-like chemotherapy plus rituximab versus CHOP-like chemotherapy alone in young patients with good-prognosis diffuse large-B-cell lymphoma: a randomised controlled trial by the MabThera International Trial (MInT) Group. Lancet Oncol. 2006; 7(5): 379-91.
Gisselbrecht C, Glass B, Mounier N, et al. Salvage regimens with autologous transplantation for relapsed large B-cell lymphoma in the rituximab era. J Clin Oncol. 2010; 28(27): 4184-90.
Friedberg JW. Relapsed/refractory diffuse large B-cell lymphoma. Hematology Am Soc Hematol Educ Program. 2011; 2011: 498-505.
Alizadeh AA, Eisen MB, Davis RE, et al. Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature. 2000; 403(6769): 503-11.
Miao Y, Medeiros LJ, Li Y, Li J, Young KH. Genetic alterations and their clinical implications in DLBCL. Nature reviews Clinical oncology. 2019; 16(10): 634-52.
Brower V. Ibrutinib promising in subtype of DLBCL. Lancet Oncol. 2015; 16(9): e428.
Lionakis MS, Dunleavy K, Roschewski M, et al. Inhibition of B Cell Receptor Signaling by Ibrutinib in Primary CNS Lymphoma. Cancer cell. 2017; 31(6): 833-843.
Wilson WH, Young RM, Schmitz R, et al. Targeting B cell receptor signaling with ibrutinib in diffuse large B cell lymphoma. Nature medicine. 2015; 21(8): 922-6.
de Vos S, Swinnen LJ, Wang D, et al. Venetoclax, bendamustine, and rituximab in patients with relapsed or refractory NHL: a phase Ib dose-finding study. Ann Oncol. 2018; 29(9): 1932-8.
Davids MS, Roberts AW, Seymour JF, et al. Phase I First-in-Human Study of Venetoclax in Patients With Relapsed or Refractory Non-Hodgkin Lymphoma. J Clin Oncol. 2017; 35(8): 826-33.
Younes A, Berdeja JG, Patel MR, et al. Safety, tolerability, and preliminary activity of CUDC-907, a first-in-class, oral, dual inhibitor of HDAC and PI3K, in patients with relapsed or refractory lymphoma or multiple myeloma: an open-label, dose-escalation, phase 1 trial. Lancet Oncol. 2016; 17(5): 622-31.
Petrich AM, Leshchenko V, Kuo P-Y, et al. Akt inhibitors MK-2206 and nelfinavir overcome mTOR inhibitor resistance in diffuse large B-cell lymphoma. Clin Cancer Res. 2012; 18(9): 2534-44.
Zeng C, Huang W, Li Y, Weng H. Roles of METTL3 in cancer: mechanisms and therapeutic targeting. J Hematol Oncol. 2020; 13(1): 117.
Wang J, Wang J, Gu Q, et al. The biological function of m6A demethylase ALKBH5 and its role in human disease. Cancer Cell Int. 2020; 20: 347.
Liu S, Li G, Li Q, et al. The roles and mechanisms of YTH domain-containing proteins in cancer development and progression. American journal of cancer research. 2020; 10(4): 1068-84.
Zhou J, Wang J, Hong B, et al. Gene signatures and prognostic values of m6A regulators in clear cell renal cell carcinoma - a retrospective study using TCGA database. Aging. 2019; 11(6): 1633-47.
Geng Y, Guan R, Hong W, et al. Identification of m6A-related genes and m6A RNA methylation regulators in pancreatic cancer and their association with survival. Ann Transl Med. 2020; 8(6): 387.
Shi H, Zhao J, Han L, et al. Retrospective study of gene signatures and prognostic value of m6A regulatory factor in non-small cell lung cancer using TCGA database and the verification of FTO. Aging. 2020; 12(17): 17022-37.
Zhao X, Cui L. Development and validation of a mA RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. American journal of cancer research. 2019; 9(10): 2156-69.
Su Y, Huang J, Hu J. mA RNA Methylation Regulators Contribute to Malignant Progression and Have Clinical Prognostic Impact in Gastric Cancer. Front Oncol. 2019; 9: 1038.
Huang H, Bai Y, Lu X, Xu Y, Zhao H, Sang X. N6-methyladenosine associated prognostic model in hepatocellular carcinoma. Ann Transl Med. 2020; 8(10): 633.
Ji L, Chen S, Gu L, Zhang X. Exploration of Potential Roles of m6A Regulators in Colorectal Cancer Prognosis. Front Oncol. 2020; 10: 768.
Han X, Liu J, Cheng G, Cui S. Gene Signatures and Prognostic Values of m6A RNA Methylation Regulators in Ovarian Cancer. Cancer Control. 2020; 27(1): 1073274820960460.
Pan J, Xu L, Pan H. Development and Validation of an m6A RNA Methylation Regulator-Based Signature for Prognostic Prediction in Cervical Squamous Cell Carcinoma. Front Oncol. 2020; 10: 1444.
Han SH, Choe J. Diverse molecular functions of mA mRNA modification in cancer. Exp Mol Med. 2020; 52(5): 738-49.
Zhao Y, Shi Y, Shen H, Xie W. mA-binding proteins: the emerging crucial performers in epigenetics. J Hematol Oncol. 2020; 13(1): 35.
Huang H, Weng H, Chen J. mA Modification in Coding and Non-coding RNAs: Roles and Therapeutic Implications in Cancer. Cancer cell. 2020; 37(3): 270-88.
Li Y, Xiao J, Bai J, et al. Molecular characterization and clinical relevance of mA regulators across 33 cancer types. Mol Cancer. 2019; 18(1): 137.
Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. mA regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. 2020; 19(1): 53.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26(12): 1572-3.
Abid A, Zhang MJ, Bagaria VK, Zou J. Exploring patterns enriched in a dataset with contrastive principal component analysis. Nat Commun. 2018; 9(1): 2134.
Sauerbrei W, Royston P, Binder H. Selection of important variables and determination of functional form for continuous predictors in multivariable model building. Stat Med. 2007; 26(30): 5512-28.
Maurer MJ, Habermann TM, Shi Q, et al. Progression-free survival at 24 months (PFS24) and subsequent outcome for patients with diffuse large B-cell lymphoma (DLBCL) enrolled on randomized clinical trials. Ann Oncol. 2018; 29(8): 1822-7.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol. 2018; 1711: 243-59.
Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12(5): 453-7.
Toor SM, Sasidharan Nair V, Decock J, Elkord E. Immune checkpoints in the tumor microenvironment. Seminars in cancer biology. 2020; 65: 1-12.
Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005; 102(43): 15545-50.
Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015; 1(6): 417-25.
Yang C, Hu Y, Zhou B, et al. The role of mA modification in physiology and disease. Cell death & disease. 2020; 11(11): 960.
Chen Y, Lin Y, Shu Y, He J, Gao W. Interaction between N-methyladenosine (mA) modification and noncoding RNAs in cancer. Mol Cancer. 2020; 19(1): 94.
Huang Y, Su R, Sheng Y, et al. Small-Molecule Targeting of Oncogenic FTO Demethylase in Acute Myeloid Leukemia. Cancer cell. 2019; 35(4): 677-691.
Han D, Liu J, Chen C, et al. Anti-tumour immunity controlled through mRNA mA methylation and YTHDF1 in dendritic cells. Nature. 2019; 566(7743): 270-4.
Cheng Y, Fu Y, Wang Y, Wang J. The m6A Methyltransferase METTL3 Is Functionally Implicated in DLBCL Development by Regulating m6A Modification in PEDF. Front Genet. 2020; 11: 955.
Han H, Fan G, Song S, et al. piRNA-30473 contributes to tumorigenesis and poor prognosis by regulating m6A RNA methylation in DLBCL. Blood. 2020; 2019003764.
Nowakowski GS, Feldman T, Rimsza LM, Westin JR, Witzig TE, Zinzani PL. Integrating precision medicine through evaluation of cell of origin in treatment planning for diffuse large B-cell lymphoma. Blood Cancer J. 2019; 9(6): 48.
Pasqualucci L, Dalla-Favera R. The genetic landscape of diffuse large B-cell lymphoma. Semin Hematol. 2015; 52(2): 67-76.
Zhou H, Zheng C, Huang D-S. A prognostic gene model of immune cell infiltration in diffuse large B-cell lymphoma. PeerJ. 2020; 8: e9658.
Hu J, Xu J, Yu M, et al. An integrated prognosis model of pharmacogenomic gene signature and clinical information for diffuse large B-cell lymphoma patients following CHOP-like chemotherapy. Journal of translational medicine. 2020; 18(1): 144.
Zhang R, Lin P, Yang X, et al. Survival associated alternative splicing events in diffuse large B-cell lymphoma. American journal of translational research. 2018; 10(8): 2636-47.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011; 144(5): 646-74.
Saleh R, Elkord E. Treg-mediated acquired resistance to immune checkpoint inhibitors. Cancer Lett. 2019; 457: 168-79.
Vaddepally RK, Kharel P, Pandey R, Garje R, Chandra AB. Review of Indications of FDA-Approved Immune Checkpoint Inhibitors per NCCN Guidelines with the Level of Evidence. Cancers (Basel). 2020; 12(3): 738.
Lesokhin AM, Ansell SM, Armand P, et al. Nivolumab in Patients With Relapsed or Refractory Hematologic Malignancy: Preliminary Results of a Phase Ib Study. J Clin Oncol. 2016; 34(23): 2698-704.
Qin S, Li A, Yi M, Yu S, Zhang M, Wu K. Recent advances on anti-angiogenesis receptor tyrosine kinase inhibitors in cancer therapy. J Hematol Oncol. 2019; 12(1): 27.
Shojaei F. Anti-angiogenesis therapy in cancer: current challenges and future perspectives. Cancer Lett. 2012; 320(2): 130-7.
Wang J, Zhu C. Anticoagulation in combination with antiangiogenesis and chemotherapy for cancer patients: evidence and hypothesis. OncoTargets and therapy. 2016; 9: 4737-46.
Belardelli F, Ferrantini M, Proietti E, Kirkwood JM. Interferon-alpha in tumor immunity and immunotherapy. Cytokine Growth Factor Rev. 2002; 13(2): 119-34.
Vidal P. Interferon α in cancer immunoediting: From elimination to escape. Scand J Immunol. 2020; 91(5): e12863.

FigureS1.jpg
Additional file 1: Figure S1. Heatmap and dendrogram of the consensus matrix for k from 3 to 9. The cluster memberships are marked by colored rectangles.

Download PDF

Version 1

posted

You are reading this latest preprint version

Expression of N6-Methyladenosine (m6A) Regulators Correlates With Immune Microenvironment Infiltration and Predicts Prognosis in Diffuse Large Cell Lymphoma (DLBCL)

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1