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


 Background: N6-methyladenosine (m6A) 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 m6A regulators and their relationship with immune microenvironment infiltration have not been illuminated.Methods: The expression of m6A 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 m6A-sorted clusters were revealed. The m6A 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 m6A low/high-risk DLBCL were studied. Gene set enrichment analysis (GSEA) was adopted to unveil the underlying molecular mechanism in m6A low/high-risk DLBCL.Results: Differentially expressed m6A regulators were able to molecularly discriminate DLBCL as two clusters based on consensus clustering and principal component analysis. A six m6A regulators-based risk prediction signature was established and validated as an independent predictor, which separated patients into m6A low- and high-risk groups. High-risk m6A 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 m6A patients. Up-regulated (PDCD1 and KIR3DL1) and down-regulated (TIGIT, DO1, and BTLA) immune checkpoint genes in high-risk m6A patients were identified. GSEA analysis unveiled high-risk m6A related to tumor proliferation associated process, while low-risk m6A related to defense response associated process. Conclusions: This study provided a comprehensive analysis for the clinical significance of m6A regulators and their association with immune microenvironment infiltration. An m6A 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 m6A risk, which may be part of the mechanism of m6A regulators in DLBCL.

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 rst-line therapy for DLBCL patients based on the results of several phase III clinical trials [2][3][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 bene t from the therapy on genomic level. For example, according to gene expression pro le, 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 7 . 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 8 . 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 14 , and AKT inhibitor MK-2206 15 , 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 satis ed yet. As a consequence, the research on exploring the DLBCL eld 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 6 -Methyladenosine (m 6 A) is one of the common and abundant post-transcriptional modi cations in mRNAs, which exerts a pivotal regulatory role in tumorigenesis and cancer progression.
The process of m 6 A is mainly achieved by three types of m 6 A regulators, that is methyltransferases ("writers"), demethylases ("erasers"), and binding proteins ("readers"). m 6 A modi cations are assembled by "writers", removed by "erasers", and deciphered by "readers". Growing studies have revealed important roles of the m 6 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 16 . 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 17 . 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 18 . The combination of m 6 A regulators for constructing prognostic models has recently been a novel signature for predicting the outcome of tumor patients. For instance, m 6 A highrisk has been clari ed as an unfavorable indicator in clear cell renal cell carcinoma 19 , pancreatic cancer 20 , non-small cell lung cancer 21 , head and neck squamous cell carcinoma 22 , gastric cancer 23 , hepatocellular carcinoma 24 , colorectal carcinoma 25 , ovarian cancer 26 , and cervical squamous cell carcinoma 27 , which indicates worse survival of patients. However, in DLBCL, no study has reported the combined prognostic value of m 6 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 rst to comprehensively investigate the prognostic value and construct a prognostic signature of m 6 A regulators in DLBCL. More importantly, we will provide new perspectives for the relationship between m 6 A and immune cell in ltration 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.
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 6 A regulators
To identify signi cantly expressed m 6 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 6 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 6 A regulators in DLBCL, a Spearman correlation analysis was conducted for the differentially expressed m 6 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 signi cant.

Exploration of m 6 A-sorted DLBCL clusters and their clinical signi cance
We wondered whether the differentially expressed m 6 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 33 . The uncovered DLBCL clusters by m 6 A regulators was validated using a Principal Component Analysis (PCA) 34 . The correlation between the m 6 A-classed clusters and the clinical parameters of DLBCL was investigated. The Chisquared test was applied when estimating the relationship between different m 6 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 6 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 6 A regulators. A p value<0.05 was used as a cut-off to screen statistically signi cant prognostic m 6 A regulators for further analysis. The Least Absolute Shrinkage and Selection Operator (Lasso) regression algorithm, a popular method for regression analysis with highdimensional features 35 , was further adopted for variable determination and regularization so as to raise the prediction accuracy and interpretability of the prognostic model. The m 6 A risk strati cation based on a risk score, which was computed as follows: risk score= . In the formula, the X i represents the expression of each m 6 A regulator in the model, while Coef i represents the coe cient of each m 6 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 classi ed as m 6 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 rst two years after diagnosis 36 , 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 6 A low/high-risk DLBCL was evaluated using Chi-squared test.
To determine whether m 6 A prognostic signature is an independent prognostic factor in DLBCL, univariate and multivariate Cox regression analyses were performed for m 6 A signature along with age, gender, stage, extranodal site, ECOG performance status, LDH ratio, and molecular subtype in GSE10846. Hazard ratio (HR), 95% con dence interval (CI), and p-value were calculated.

Exploration of m 6 A signature and immune cell in ltration
We wonder whether m 6 A risk strati cation 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 in ltration difference between m 6 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) 39 in m 6 A low-risk and high-risk groups were also investigated to uncover the potentially possible relationship between m 6 A and immunotherapy. A p-value under 0.05 was deemed to be statistically signi cant.

Implementation of gene set enrichment analysis (GSEA)
GSEA is a knowledge-based computational approach for interpreting a priori de ned set of genome-wide expression pro les, which helps determine statistically signi cant 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 signi cant gene ontology terms (GO), which includes GO biological process, GO cellular component, and GO molecular function, in m 6 A low/high-risk groups. GSEA version 4.10 software and Molecular Signatures Database v7.2 were applied for the current analysis.

Results
Differentially expressed m6A regulators in DLBCL The expression of the 22 m 6 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 6 A regulators. And the violin plot in Figure 2B visualized the expression difference of each m 6 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 6 A regulators. We can observe that the correlations among the m 6 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 6 A-sorted DLBCL clusters and clinical implication
As we can observe in Figure 3A-3C and Figure S1, when clustering the samples using a speci ed cluster count (k) from 2 to 9, only k=2 showed the best performance to distinguish the DLBCL samples by the differentially expressed m 6 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 signi cant between cluster 1 and cluster 2, which indicated that m 6 A regulators were capable to differentiate the DLBCL samples into molecularly distinguishable subtypes. The nding 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 6 A-sorted DLBCL clusters are reliable. The clinical signi cance of the two m 6 Asorted DLBCL clusters was also investigated and demonstrated in Figure 4. Figure 4A displayed the expression heatmap of the differentially expressed m 6 A regulators in m 6 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). Six out of 19 m 6 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 6 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 6 A low-risk have higher overall survival than those in m 6 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 6 A signature in the training set.
The m 6 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 signi cance 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 6 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 signi cant difference was observed in the K-M survival curve analysis for the combined samples (p<0.001). Patients with m 6 A low-risk showed better survival than those in m 6 A high-risk group ( Figure   6E). Figure 6F displayed the predictive e cacy of the m 6 A signature, which yielded an AUC for 1-, 2-, and 5-year of 0.605, 0.640, and 0.652, respectively.

Estimation of immune cell in ltration
The immune cell in ltration in m 6 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 6 A low/high-risk groups as a heatmap. The immune cell in ltration 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 6 A low-risk and high-risk groups. The B cell naïve, NK cells resting, NK cells activated, and Mast cells activated were signi cantly upregulated in m 6 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 signi cantly down-regulated in m 6 A highrisk 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 signi cantly increased in m 6 A high-risk group, while TIGIT (p<0.001), IDO1 (p=0.016), and BTLA (p=0.008) was signi cantly decreased. Regards LAG3, CTLA4, VISTA, and HAVCR2, no statistical signi cance was observed.
Enriched pathways in m6A low/high-risk DLBCL To reveal the potential biological function difference between the m 6 A low-risk and high-risk groups, GSEA analysis was carried out. The top 10 signi cantly enriched terms in the m 6 A low-risk group and high-risk group were selected for displaying, respectively. As we can see in Figure 10

Discussion
Emerging evidence substantiated that m 6 A modi cation exerts an indispensable regulatory effect in neuronal disorders, osteoporosis, metabolic disease, viral infection, and various cancers 42,43 . Targeting m 6 A regulators for cancer therapy has been a closely focused eld 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 signi cantly suppress the progression of AML primary cells in xenotransplanted mice 44 . In recent years, the prominent roles of m 6 A modi cation in cancer immunotherapy has been gaining more and more attention. Han et al. recently discovered that m 6 A-binding protein YTHDF1 can control anti-tumor immunity through recognizing m 6 A-marked transcripts encoding lysosomal proteases to increase their translation in dendritic cells. Speci cally, the de ciency of YTHDF1 elevated the cross-presentation of tumor antigens and antigen-speci c CD8+ T cell antitumor response. More remarkably, loss of YTHDF1 enhanced the therapeutic e cacy of PD-L1 checkpoint blockade in vivo 45 . Despite the upsurging wave of m 6 A research in cancer, few research has reported the pathological role of m 6 A regulators and their clinical signi cance in DLBCL. Cheng et al. once reported that downregulated methyltransferase METTL3 functionally inhibited the DLBCL cell proliferation through reducing the m 6 A methylation and total mRNA level of pigment epithelium-derived factor 46 . Another study revealed that up-regulated piRNA-30473 was associated with aggressive phenotype and poor prognosis of DLBCL patients by virtue of m 6 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 6 A level 47 . Nevertheless, the above studies only focused on a single m 6 A regulator, lacking a comprehensive understanding of the potential clinical value of different m 6 A regulators as a whole. Given that, our current study primarily focused on evaluating the prognostic value of 22 known m 6 A regulators, constructing a prognosis predictive model, and investigating the potential molecular mechanism of m 6 A risk in DLBCL.
The majority of m 6 A regulators showed differential expression in DLBCL and their interactions could be complicated. The same m 6 A regulator could be positively correlated to a regulator while negatively correlated to another regulator. This indicated the m 6 A regulation network in DLBCL could be sophisticated, which deserves further analysis. More interestingly, we found that based on differentially expressed m 6 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 signi cant 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 pro le. 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 48 . This suggests more positive therapy and follow-up measures should be taken when managing ABC-DLBCL patients.
Similarly, different m 6 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 49 . Whether m 6 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 6 A-based prognostic model for predicting the overall survival of DLBCL patients. Patients with m 6 A high-risk was veri ed to yield poor outcomes. m 6 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. In the current study, the 1-year, 2-year, and 5-year predictive AUC of m 6 A signature for the 380 DLBCL patients were 0.605, 0.640, 0.652, respectively. Compared with previously published prognostic models, the m 6 A signature showed passable predictive performance, which is quali ed 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 broblasts 53 . Although in ltrated 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. 54 . 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 rst 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 55 . However, in relapsed/refractory DLBCL, PD-1 blockade therapy has been disappointing, yielding an objective response of merely 36% 56 . 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 in ltration is an essential step. As a result, our current study provided novel ndings on the relationship between m 6 A and immune cell in ltration, 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 6 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 6 A high-risk group. These indicated the complicated regulatory relationship between m 6 A and immune microenvironment. How m 6 A exactly in uences immune cell in ltration 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 6 A low-risk and high-risk DLBCL is a noteworthy aspect. Therefore, we explored the possible molecular mechanism in m 6 A low/high-risk DLBCL. We discovered that in m 6 A high-risk group, the enriched biological function mainly included cell cycle, DNA replication, transcription, post-transcriptional modi cation, and DNA repair relative pathways, which are manifestations of tumor malignancy features. And in m 6 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][58][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 6 A high-risk group suffered from worse survival than those in m 6 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 6 A-based prognostic signature and their correlation with tumor immune microenvironment in ltration, several limitations should be acknowledged in this study. First, even though we constructed and validated the m 6 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 6 A risk and immune microenvironment in ltration, 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 6 A regulators and tumor immune microenvironment in ltration.

Conclusions
This study comprehensively investigated the clinical signi cance of multiple m 6 A regulators and established a novel m 6 A risk scoring signature for predicting the survival of DLBCL patients for the rst time, which achieved satisfactory predictive performance. More importantly, we unveiled several highlights of immune cell in ltration, immune checkpoint genes and underlying molecular enriched pathways in m 6 A low/high-risk DLBCL pioneeringly, which shed light on the potential regulatory relationship between m 6 A and immune microenvironment in ltration. Nevertheless, the obtained results especially the exact mechanism on how m 6