M6A regulator expression patterns predict the immune microenvironment and prognosis of non-small cell lung cancer

The m6A methylation modification is one of the most common mRNA modifications, and involved in a variety of biological processes, such as cell death, cancer stem cell formation and tumorigenesis. Increasing evidences have demonstrated that the expression patterns of m6A regulators are significantly correlated with PD-L1 level some solid tumors, but few study has explored the function of m6A regulators in the immune microenvironment and prognosis in non-small cell lung cancer (NSCLC). Survival analysis was independently conducted for 20 m6A regulators to explore their prognostic value in NSCLC, and then the prognostic risk model based on m6A regulator expression profiles is built to stratify NSCLC patients. Also, the correlation analysis between immune infiltrating cells and m6A regulators is used to reveal the impact of m6A on immune microenvironment of NSCLC. Furthermore, to explore the function of m6A as biomarker of anit-PD-L1 therapeutic effect, we explored the associations of tumor mutation burden (TMB) and PD-L1 levels to 20 m6A regulator expression patterns in NSCLC. First, the expressions of 20 m6A regulators in NSCLC tissues were significantly increased compared to normal tissues. Survival analysis revealed that three genes, METTL3, HNRNPC and VIRMA, were markedly correlated to the prognosis of NSCLC patients. In particular, cox regression analysis verified that METTL3 could be used as an independent prognostic factor to predict the survival rate of NSCLC patients. Second, the risk prognostic model built on seven m6A regulators can effectively stratify NSCLC patients, and the low-risk subgroup had better prognosis compared to high-risk group. Finally, a few m6A regulators showed significant associations with immune microenvironment, as well as TMB and PD-L1 level, suggesting that the m6A RNA methylation is indicative of therapeutic effect of anti-PD-L1 treatment. Our study identified some m6A regulatory factors as independent risk factors for the prognosis of NSCLC, and the expression patterns of m6A regulators are also correlated to the immune infiltration, as well as TMB and PD-L1 level in NSCLC. The m6A regulators could be used as biomarkers indicative of immunotherapy to NSCLC patients.


Introduction
Non-small cell lung cancer (NSCLC) is the leading cause of lung cancer-related mortality around the world (Gupta et al. 2020). Most NSCLC patients suffer from poor prognosis Xue Liu and Changsheng Ma contributed equally to this study.
* Zhiqiang Sun jungler@sina.cn * Judong Luo judongluo@163.com with low 5-year survival rate, and its incidence and mortality rate is increasing in recent years (Chen et al. 2014;Bray et al. 2018). NSCLC includes two major histological subtypes: lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD) (Coudray et al. 2018). Although surgery is the major clinical treatment, it is restricted to early stage lung cancer. Almost all patients with advanced and unresectable NSCLC have limited treatments other than chemotherapy and radiotherapy, which often lead to serious side effects (Hsu et al. 2020;Jin et al. 2019). In fact, chemotherapy has only 50% remission rate to NSCLC (Nagasaka and Gadgeel 2018). Therefore, identifying novel biomarkers and therapeutic targets for NSCLC patients is an urgent but challenge task. In recent years, immune checkpoint inhibitors (ICI) has showed strong clinical effects on the therapy of a variety of tumors (Wang et al. 2016). Its clinical advantage includes but not limited to continuous anti-tumor immune response with relatively weak side effect, low recurrence rate, and even complete remission for some advanced cancers. However, the fraction of NSCLC patients that benefits from ICI treatment remains very limited (Herbst and Boshoffc 2018;Hugo et al. 2016). Tumor mutation burden (TMB), microsatellite instability (MSI) and PD-L1 level have relative Fig. 1 Differential expression analysis of m6A regulators in NSCLC and normal tissues, as well as difference in LUSC and LUAD subtypes. a Heatmap of M6A regulators expression profiles in both NSCLC and normal tissues. b Volcano Plot of differentially expressed genes in NSCLC relative to normal tissues, in which blue indicates down-regulated genes and red indicates up-regulated ones. c Boxplots of expression levels of 20 m6A-related genes in two subtypes and normal tissues, respectively value as predictive biomarkers of anti-PD-L1 treatment in advanced/metastatic NSCLC (Brody et al. 2017;Larsen et al. 2019;Camidge et al. 2019). However, more than 50% of tumors with high PD-L1 level are resistant to ICI, and the underlying mechanism of drug resistance remains unclear (Anichini et al. 2020;Zhu and Lang 2017).
N6 Methyladenosine (m6A) is a type of prevalent posttranscriptional modification at the sixth nitrogen atom of adenosine on eukaryote mRNA. The m6A methylation mediates more than 60% of RNA methylation , and involves in the development of tumors by regulating the mRNA expression level of related oncogenes or suppressor genes (Huang et al. 2020). The m6A methylation level mainly depends on the expression of m6A methylation regulators (Ma et al. 2019), including the methyltransferases ("writers") and demethylases ("eraser"). Its status is also governed by the activities of m6A readers such as YTH N6-methyladenosine RNA binding proteins (YTHDF) (Wang and He 2014). The role of m6A regulatory factors has attracted much attention. For example, it is reported that IGF2BP2 is related to the occurrence, development and prognosis in colorectal cancer . The overexpression of METTL3 promotes the proliferation and migration of hepatocellular carcinoma cells , while RBM15B is proven to play tumor inhibitory role in renal cell carcinoma (Fang et al. 2020). More and more evidence indicated that m6A plays an important role in NSCLC. Also, previous studies have shown that yes associated protein (YAP) is the main mediator protein of Hippo pathway, which has been proved to promote tumor progression, drug resistance and metastasis of NSCLC (Hsu et al. 2020). M6A methylation can promote YAP translation and indirectly induce NSCLC resistance and metastasis by regulating MALAT1-miR-1914-3p-YAP axis, thus enhancing the activity of YAP ). In addition, the expression of YAP protein is closely relevant to PD-L1 in tumor microenvironment, and the mechanism of YAP regulating PD-L1 expression in human NSCLC has been clarified (Hsu et al. 2019;Miao et al. 2017). To sum up, m6A and PD-L1 are linked through YAP protein and participate in the resistance to ICI treatment. Nevertheless, the interactions between M6A methylation and ICI therapeutic effect remain unknown in NSCLC to date.
In this paper, we aim to systematically evaluate the impact of m6A methylation regulators on the prognosis and anti-PD-L1 therapy of NSCLC patients. In particular, we explore m6A methylation prognostic value in the lung squamous cell carcinoma and lung adenocarcinoma, respectively. We established a risk model and identified independent prognostic factors, and further researched the effect of m6A methylation regulatory factors on immune infiltration in NSCLC. Moreover, we also run the analysis of the NSCLC TMB to

M6A regulators differentially expressed in NSCLC
Differential expression analysis was conducted for the 20 m6A regulators in NSCLC compared to normal tissues. As shown in Fig. 1a, the expression levels of m6A regulatory genes were significantly different between cancer and normal groups. Figure 1b shows the up-regulation and downregulation of 3337 genes (Supplementary Table 1). Taking the overlap with 20 m6A regulatory genes, we found that the expressions of IGF2BP1, IGF2BP2 and IGF2BP3 were significantly up-regulated in NSCLC tissues. In contrast, there was little difference in expression patterns of 20 m6A regulators between LUSC and LUAD subtypes (Fig. 1c). These results suggest that m6A regulators may play an important biological role in the occurrence and development of NSCLC, whereas they did not discriminate LUSC and LUAD subtypes.

Correlation among m6A regulators in NSCLC
As shown in Fig. 2a m6A regulator network analysis was conducted to reveal the correlation among 20 m6A regulators in NSCLC. We found that m6A regulators belonging to same functional category showed similar expression patterns. Spearman correlations were used to describe the correlation between m6A regulators (Fig. 2b), and we observed significant correlation among writers, readers and erasers. For example, a few gene pairs, including RBM15B and ZC3H13, RBM15 and HNRNPA2B1, were strongly correlated. It is suggested that m6A regulatory factors are interrelated and function to form complex regulatory components.

M6A regulators are indicative of NSCLC prognosis
Survival analysis was conducted with respect to single m6A regulator on the TCGA cohort, and the Kaplan-Meier curves were shown in Fig. 3. The results indicated that METTL3 (Writers), VIRMA (Writers) and HNRNPC (Readers) were statistically significant in association to overall survival.
Particularly, METTL3 high expression improved NSCLC prognosis, while HNRNPC and VIRMA high expression level means poor prognosis (Fig. 3a-c). We further conducted prognosis analysis regarding m6A regulators separately on LUAD and LUSC subcohorts. As shown in Fig. 3d-i, six m6A-related genes, including RBM15, HNRN-PA2B1, HNRNPC, IGF2BP1, IGF2BP2 and IGF2BP3, were negatively correlated to LUAD prognosis, while WTAP and FTO had significant impact on LUAC prognosis (Fig. 3j-k). Our results showed that the expression patterns of m6A regulators were indicative of the prognosis of both LUAD and LUSC patients. In particular, m6A readers have more spread impact on the prognosis than m6A erasers and writers.

Establishment of risk prognostic model
The univariate Cox regression analysis found that VIRMA, METTL3 and IGF2BP1 were significantly associated with the prognosis of patients (P < 0.05; Fig. 4a). Next, the multivariate cox regression analysis further verified that HNRN-PA2B1 and METTL3 were independent prognostic factors of NSCLC patients (P < 0.05, Fig. 4b). In particular, METTL3 low expression reflected strong statistical significance as a poor prognostic indicator. Cox regression analysis was also performed on LUAD and LUSC subtypes respectively. In LUAD, we found three independent factors, HNRNPA2B1, IGF2BP1 and HNRNPC, which were indicative of poor prognosis ( Supplementary Fig. 1a, b).
To build a comprehensive risk prognostic model, lasso regression analysis was conducted regarding the expression of 20 m6A regulators on TCGA cohort. Seven regulators, YTHDC2, HNRNPC, HNRNPA2B1, METTL3, VIRMA and IGF2BP1, were included as risk factors. The risk score of each patient was calculated as: r isk sco r e = (0.034 9) * IGF2 BP1 + (− 0.0 428) * Y THDC1 + ( 0.115) * VIRMA + (− 0.19) * METTL3 + ( 0.1 4 6) * HNRN PA2B1 + (0.041) * HN RNPC + (− 0. 050 7) * YTHDC2. According to the median risk score, patients were divided into high-risk group and low-risk group (Fig. 5a). As shown in Fig. 5b, the oncogenic m6A regulators were highly expressed in the high-risk group, while the protective m6A regulators were up-regulated in the low-risk group. The Kaplan-Meier curve in Fig. 5c shows that the low-risk group has significantly better overall survival than high-risk group. The prediction models for 1-, 3-and 5-year survival status were built based on the seven m6A regulators, and the ROC curves are shown in Fig. 5d, respectively. Besides, we established the risk prognostic model for LUAD and LUSC subtypes respectively. In LUAD, the risk model included 14 m6A regulators ( Supplementary Fig. 2), and verified that the patients in the low-risk group had significantly better prognosis than high-risk group (Fig. 5e). The ROC-AUCs for the  (Fig. 5f). However, we did not get statistically significant result on LUSC subtype cohort. These results suggest that the expression patterns of m6A regulators conferred stronger predictive capacity in LUAD subtype.

Immune infiltration differed in high-and low-risk group
To investigate the relationship between m6A and immune microenvironment in NSCLC, we computed the proportions of 22 type of immune cells using CIBER-SORT. We compared the immune infiltration between high-risk and lowrisk groups, and found that Macrophages M0 infiltration was the most obvious in the high-risk group, Plasmacells infiltration was most obvious in the low-risk group (Fig. 6a). Within LUAD subtype, only T cell CD4 naive infiltration was obvious (Fig. 6b). These results demonstrate different degree of immune infiltration among different subtypes. We also calculated the correlation between the prognostic genes and 22 immune cells in NSCLC (Fig. 6c). It can be seen that B cell naive is closely related to YTHDC1, METTL3 and HNRNPA2B1, CD4 T memory cells activated showed significant association with HNRNPC. Within LUAD subtype, we run same analysis and found that the genes related to immune infiltration increased in the high-risk group, especially METTL3, HNRNPA2B1 and ZC3H13 (Fig. 6d). Moreover, some immune infiltrating cells increased in the high-risk group, especially the T cells follicular helper and T cells CD4 memory activated.

Correlation between m6A and anti-PD-L1 biomarkers
The PD-L1, an biomarker to evaluate the therapeutic effect of anti-PD-L1 therapy, showed strong correlation to m6A regulators. The p value and correlation coefficient between PD-L1 (CD274) and 20 m6A regulatory factors were shown in Supplementary Table 2. YTHDF2 has strong positive correlation to PD-L1 level, while YTHDC2 has negative correlation to m6A regulators. Besides, we explore the correlation between TMB and m6A regulators (Supplementary Table 3). The results showed that there was certain relationship c and e K-M survival curves and in high-and low-risk groups of NSCLC and LUAD, d and f are ROC curves and AUC values for the prediction of 1-, 3-, and 5-year survival status using the risk scores of NSCLC and LUAD Fig. 6 Immune microenvironment between high-and low-risk groups. a and b showed the difference of immune infiltrating cells in high-and low-risk groups of NSCLC and LUAD subtype respec-tively. c and d showed the correlation of high risk genes and immune infiltrating cells in NSCLC and LUAD respectively between m6A and TMB, and according to the above genes with high correlation coefficient have a significant impact on prognosis. Also, m6A had higher correlation to TMB in LUAD (Supplementary Table 3), whereas no statistical significance in LUSC was found (P > 0.05).

Data sources
Transcriptomic profiles RNA-seq data and clinical information of NSCLC patients were obtained from the public database TCGA https:// portal. gdc. cancer. gov/. Patients with incomplete survival information were excluded from the study. Finally, a total of 1013 NSCLC samples (501 LUSC samples, 513 LUAD samples) and 108 non-tumor samples were collected for our analysis. The expression profiles were plotted using R package "ggplot2" and "pheatmap".

Survival and cox regression analysis
Kaplan-Meier survival analysis were performed using R package "survival" and "surviviner". Univariate and multivariate Cox regression were conducted to identify the independent risk factors of prognosis. The forest maps were constructed by R package "forestplot", which shown the P value and HR (95% CI) of each m6A related gene.

Establishment of risk model
Lasso regression analysis was carried out to establish prognosis risk model based on the expression levels of m6A regulators. Patients were divided into two high-and lowrisk groups according to the risk scores using median as threshold. The R package "glmnet" was used to conduct Lasso regression.

Immune microenvironment analysis
CIBERSORT http:// ciber sort. stanf ord. edu/ as used to calculate the composition of immune cells from the gene expression profiles. We calculated the infiltration levels of 22 immune cells in NSCLC and LUAD, as well as the difference of immune score between different groups.

Correlation between m6A regulators and PD-L1
To explore the potential indication of m6A methylation to ICI efficacy in NSCLC, Spearman correlations were computed to describe the associations between m6A regulators and conventional indicators of ICI therapy. The heatmap was plotted by R package "ggstatsplot" to present the correlation between PD-L1 and 20 m6A regulators. Also, the correlation between TMB scores and 20 m6A regulatory factors were presented. In all the above experiments, P value < 0.05 was considered as significance threshold in all statistical tests.

Statistical analysis
All pictures in this paper are based on TCGA database and drawn with the help of R language. Significance was considered at either P < 0.05 was considered statistically significant and highlighted an asterisk in the figures, while P < 0.01 were highlighted using two asterisks and P < 0.001 highlighted using three asterisks in the figures.

Discussion
M6A has recently been implicated in the occurrence of many kinds of human cancer (Chen and Wong 2020), such as breast cancer (Anita et al. 2020) and prostate cancer (Yuan et al. 2020). This post transcriptional modification is dynamically and reversibly mediated by different regulatory factors, including methyltransferase, demethylase and m6A binding protein (Shi et al. 2020). Up to date, few study paid attention to the association of m6A to tumor immune microenvironment and anti-PD-L1 therapeutic effect.
In this study, we confirmed the expression pattern and prognostic value of 20 m6A regulatory factors in NSCLC, and revealed their relationship with immune infiltration and TMB. In particular, the m6A readers of IGF2BP protein family are most significantly up-regulated in NSCLC compared with the other m6A genes. Our downstream study suggested that IGF2BPs were associated with patient prognosis. According to previous studies, IGF2BP protein family is known to be an active regulator of tumor progression and metastasis of NSCLC (Bell et al. 2013). Its functional mechanism could be that circNDUFB2 participates in the degradation of IGF2BPs and activation of anti-tumor immunity during NSCLC progression via the modulation of protein ubiquitination and degradation, as well as cellular immune responses . Also, we found that YTHDC1 has similar expression pattern to HNRNPA2B1, which are both m6A nuclear readers and regulate mRNA splicing (Roundtree et al. 2017;Alarćon et al. 2015). Our further study found that HNRNPA2B1 was an independent risk factor for prognosis in has not been reported yet. Both YTHDC1 and HNRNPA2B1. It is suggested that YTHDC1 may be indirectly related to the prognosis of NSCLC patients.
As prognostic risk factors of NSCLC, our univariate and multivariate Cox regression analysis verified that METTL3 could be used as an independent risk factor affecting the prognosis in NSCLC. We also found that m6A regulators have greater effect on the prognosis of LUAD than LUSC. Three regulators, IGF2BP1, HNRN-PA2B1 and HNRNPC, could be used as independent risk factors affecting the prognosis in lung adenocarcinoma. Within LUSC subtype, METTL3, VIRMA and HNRNPC genes were significantly correlated with the prognosis of patients. This is consistent with the current research results (Xu et al. 2021;Yan et al. 2019;Wang et al. 2020). However, the regression analysis yield no statistically significant results in LUSC.
It has been reported that m6A is associated with tumor immune infiltration. YTHDC1 is significantly related to immune infiltration in endometrial carcinoma (Ma and Yang 2021), and YTHDC2 is closely related to immune infiltration in head and neck squamous cell carcinoma (Li et al. 2020b). We also found significant correlation between m6A regulators and immune infiltrating cells in NSCLC. Three genes, YTHDC1, YTHDC2 and HNRNPC, are closely related to immune infiltration in NSCLC. Significant difference in immune infiltration between different NSCLC subtypes can be observed. In LUAD, we found that METTL3, HNRNPA2B1 had a great correlation with T cell follicular helper infiltration compared to LUSC. As high concentration of Treg immune cells can lead to immune escape and promote tumor progression in ESCC (Fakhrjou et al. 2014;Guo et al. 2021). The expression of METTL14 and ZC3H13 is negatively correlated with Treg cells in breast cancer (Gong et al. 2020). ALKHB5 inhibits immune cell aggregation in the tumor microenvironment to regulate the response to anti-PD-1 therapy (Li et al. 2020a). We speculate that there is a potential relationship between m6A regulatory factors and the efficacy of immunotherapy for NSCLC, which is worthy of further study.
Finally, the correlation between m6A regulators and biomarkers of anti-PD-L1 therapy in NSCLC was explored. We found YTHDC2 is most closely related to PD-L1, which is consistent with previous result in gastric adenocarcinoma (Mo et al. 2020). A recent study found that IGF2BP3 inhibits CD8 + T cell response by promoting deubiquitination of PD-L1 in non-small cell lung cancer to promote tumor immune escape (Liu et al. 2021). This tells us that m6A regulator can affect the prognosis of patients by affecting PD-L1 in NSCLC. In addition, the Spearman correlation between TMB and m6A regulators NSCLC showed that the HNRNPC highly related to TMB. Consistent to previous study, HNRNPC low expression means better prognosis. Within LUAD sub-type, the high expression of HNRNPC and IGF2BP3 significantly improve the prognosis. This suggests that they can be used as potential prognostic biomarkers.

Conclusions
In conclusion, this paper proposed a systematic evaluation of the prognostic value for m6A regulatory factors in NSCLC, and their prognostic ability in LUAD and LUSC subtypes was also conducted, respectively. We found that METTL3 could be an independent risk factor for the prognosis of NSCLC, while IGF2BP1, HNRNPA2B1 and HNRNPC could play as independent risk factors in LUAD subtype. The evaluation between anti-PD-L1 biomarkers and expression patterns of m6A regulators not only further the understanding of anti-PD-1/L1 therapy, but also suggested effective immunotherapy strategy.