NEFM DNA methylation correlates with immune infiltration and survival in breast cancer

Background This study aims to determine whether NEFM (neurofilament medium) DNA methylation correlates with immune infiltration and prognosis in breast cancer (BRCA) and to explore NEFM-connected immune gene signature. Methods NEFM transcriptional expression was analyzed in BRCA and normal breast tissues using Oncomine and Tumor Immune Estimation Resource (TIMER) databases. The relationship between NEFM DNA methylation and NEFM transcriptional expression was investigated in TCGA. Potential influence of NEFM DNA methylation/expression on clinical outcome was evaluated using TCGA BRCA, The Human Protein Atlas and Kaplan–Meier plotter databases. Association of NEFM transcriptional expression/DNA methylation with cancer immune infiltration was investigated using TIMER and TISIDB databases. Results High expression of NEFM correlated with better overall survival (OS) and recurrence-free survival (RFS) in TCGA BRCA and Kaplan–Meier plotter, whereas NEFM DNA methylation with worse OS in TCGA BRCA. NEFM transcriptional expression negatively correlated with DNA methylation. NEFM DNA methylation significantly negatively correlated with infiltrating levels of B, CD8+ T/CD4+ T cells, macrophages, neutrophils and dendritic cells in TIMER and TISIDB. NEFM expression positively correlated with macrophage infiltration in TIMER and TISIDB. After adjusted with tumor purity, NEFM expression weekly negatively correlated with infiltration level of B cells, whereas positively correlated with CD8+ T cell infiltration in TIMER gene modules. NEFM expression/DNA methylation correlated with diverse immune markers in TCGA and TISIDB. Conclusions NEFM low-expression/DNA methylation correlates with poor prognosis. NEFM expression positively correlates with macrophage infiltration. NEFM DNA methylation strongly negatively correlates with immune infiltration in BRCA. Our study highlights novel potential functions of NEFM expression/DNA methylation in regulation of tumor immune microenvironment. Supplementary Information The online version contains supplementary material available at 10.1186/s13148-021-01096-4.


Introduction
Breast cancer (BRCA) is the most common malignancy among females worldwide. Clinical outcome has been improved over the past two decades with currently available modalities, including surgery, chemotherapy, endocrine therapy, radiotherapy, and targeted therapy, but BRCA treatment remains challenging because of high heterogeneity [1][2][3]. Immunotherapy is emerging as new therapeutics in BRCA. Several immunotherapeutic agents have been explored in various tumors, including adoptive cell therapies, vaccines, oncolytic viruses, and most notably immune check point blockade (ICB). Agents of ICB such as inhibitors of cytotoxic T-lymphocyte-associated antigen (CTLA-4), programmed cell death receptor1 (PD-1), and programmed cell death1 ligand1 (PD-L1) have been widely used in solid tumors, refractory cancers harboring microsatellite instability and classical Hodgkin lymphoma. Notably, anti-PD-L1 antibody atezolizumab in combination with nab-paclitaxel has been approved for the treatment of metastatic triple-negative breast cancer (TNBC) [1][2][3][4]. Expression of PD-L1 in infiltrating immune cells is required for response to atezolizumab plus nab-paclitaxel in IMpas-sion130 trial [5].
Tumor-infiltrating lymphocytes (TILs) comprise a mixture of cytotoxic T cells, helper T cells, B cells, macrophages, natural killer cells, and dendritic cells, which have been observed in many solid tumors, including BRCA. TILs may provide prognostic and predictive clues in BRCA and other cancers. To date, robust predictive biomarkers for immunotherapy have not been established in BRCA [1,6]. TILs are more commonly observed at higher levels in TNBC and HER2-positive BRCA compared with estrogen receptor (ER)-positive and HER2negative BRCA [1,7,8]. TILs may be associated with improved prognosis and better response rates to neoadjuvant therapy [7].
In ovine amniotic epithelium (oAECs) isolated from late amnia, NEFM mRNA levels were significantly increased, while immunomodulatory effect of inhibiting lymphocyte proliferation was lost, and global DNA methylation was enhanced. Myelin oligodendrocyte glycoprotein induced incomplete tolerance of CD4 (+) T cells specific for myelin and neuronal self-antigen NEFM in mice [24,25]. These studies suggest that NEFM is related to immune response. However, the relationship of NEFM with TILs in tumor progression or immunotherapy remains unclear.
In this study, association between NEFM expression and prognosis of BRCA was explored using TCGA (The Cancer Genome Atlas), The Human Protein Atlas, Oncomine and Kaplan-Meier plotter. In addition, association between NEFM DNA methylation and NEFM transcriptional expression was analyzed using BRCA samples in TCGA. Moreover, the relationship of NEFM transcriptional expression and NEFM DNA methylation with tumor-infiltrating immune cells was investigated in TCGA BRCA based on Tumor Immune Estimation Resource (TIMER) and TISIDB (tumor-immune system interactions).

The genes and pathways connected with NEFM transcriptional expression/DNA methylation
Differentially expressed genes associated with NEFM expression or NEFM DNA methylation were profiled through comparison between NEFM/NEFMmet high and low groups in TCGA BRCA cohort. Totally, 164 up-regulated and 546 down-regulated genes were significantly associated with NEFM expression, while 103 up-regulated and 641 down-regulated genes were significantly associated with NEFM DNA methylation (with absolute value of log2foldchange > 1, and adjust p value < 0.05; Fig. 4a, d; Additional files 1, 2: Tables S1-S2). The top 50 differentially expressed genes were presented as expression heatmaps (Fig. 4b, e). Critical signal transduction pathways involved in NEFM expression included neuroactive ligand-receptor interaction, protein digestion and absorption, chemical carcinogenesis, cAMP signaling pathway, IL-17 signaling pathway, and cytokine-cytokine receptor interaction by KEGG enrichment analysis (Fig. 4c). Cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, primary immunodeficiency, hematopoietic cell lineage, chemokine signaling pathway, neuroactive ligand-receptor interaction, T cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, IL-17 signaling pathway, and NF-kappa B signaling pathway were the top 10 pathways closely associated with NEFM DNA methylation based on KEGG enrichment analysis. Notably, some pathways involved in immune response such as Th17 cell differentiation, graft-versus-host disease, intestinal immune network for IgA production, Th1 and Th2 cell differentiation, as well as PD-L1 expression and PD-1 checkpoint pathway in cancer, were significantly associated with NEFM methylation (Fig. 4f ).

Correlation of NEFM transcription/DNA methylation with immune infiltration in breast cancer
Relationship of NEFM transcriptional expression/DNA methylation with immune infiltration in breast cancer was assessed using correlation analysis and TISIDB databases. NEFM transcriptional expression was weakly (R < 2) to moderately (2 < R < 3) positively associated with infiltration levels of macrophages and neutrophils using correlation analysis and TISIDB database (Fig. 5a, c). NEFM transcriptional expression was weakly positively associated with infiltration levels of CD8 + T cells, CD4 + T cells by correlation analysis, whereas weakly to moderately negatively associated with infiltration levels of activated CD8 + T cells, activated CD4 + T cells in TISIDB database (Fig. 5a, c). Since the different results from correlation analysis and TISIDB databases, TIMER gene modules were applied to evaluate the relationship of NEFM transcriptional expression with immune infiltration in breast cancer. In TIMER gene modules, NEFM transcriptional expression positively correlated with infiltration levels of CD8 + T cell, macrophage, neutrophil, and dendritic cell, and negatively correlated with infiltration level of B cell and tumor purity, whereas not with infiltration level of CD4 + T cell. After adjusted with tumor purity, NEFM expression weekly negatively correlated with infiltration level of B cell and positively correlated with macrophage and CD8 + T cell. NEFM DNA methylation was moderately to strongly (R > 3) negatively associated with infiltration levels of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells using correlation analysis and TISIDB database (n = 785) (Fig. 5b, d). Interestingly, NEFM transcriptional expression weakly negatively correlated to infiltration levels of M2 macrophage, while NEFM DNA methylation weakly negatively correlated to infiltration levels of M1 macrophage and positively correlated to infiltration levels of M2 macrophage with correlation analysis (Fig. 5e). Collectively, NEFM expression positively correlated with macrophage infiltration in TIMER and TISIDB; after adjusted with tumor purity, NEFM expression also weekly negatively correlated with infiltration level of B cell and positively correlated with CD8 + T cell in TIMER gene modules. However, NEFM DNA methylation was significantly negatively associated with immune infiltration in breast cancer. NEFM expression/ DNA methylation might play a specific role in immune infiltration in BRCA.

Discussion
In this study, NEFM transcriptional expression was downregulated and negatively correlated with DNA methylation in breast cancer. Enhanced DNA methylation on six loci within NEFM located on promoter CpG island or shore was associated with poor survival. Besides, NEFM transcriptional expression correlated with better prognosis and correlated with increased macrophage; after normalized with tumor purity, NEFM expression correlated with increased CD8 + T cell, whereas decreased B cell infiltration in BRCA. NEFM DNA methylation correlated with decreased infiltration levels of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells and diverse immune markers. Therefore, our study provides new evidence to support a role of NEFM transcriptional expression/NEFM DNA methylation in BRCA. NEFM polypeptide is one of the four subunits comprising neurofilaments, the most abundant intermediate filaments in nervous system. In addition, NEFM, NEFL and  [2,33,34]. Therefore, the identification of novel tumor-suppressive genes targeted by promoter methylation can reveal tumor-suppressive pathways in breast carcinogenesis and explore alternative approaches for diagnostic and therapeutic evaluation. We investigated the relationship between 24 loci within NEFM gene and prognosis in BRCA and identified enhanced NEFM DNA methylation of six out of 16 loci located on promoter CpG island related to poor survival. Elevated levels of anti-NEFM antibodies were detected in various neurological diseases, including autoimmune diseases, non-immune-mediated conditions, and even in individuals being considered normal or with disorders unrelated to intrathecal space, such as multiple sclerosis, schizophrenia, spondylogenic headache or neurastenia. Therefore, anti-NEFM antibodies may be regarded as natural circulating auto-antibodies [35,36]. Poly-specific T cells targeting distinct self-antigens have been identified in healthy individuals as well as in the context of autoimmunity. T cell recognizes NEFM protein, with implications for aggravation and perpetuation of central nervous system autoimmunity [37]. However, whether NEFM is involved in regulating antitumor immunity with clinical significance in breast cancer remains unknown. In this study, positive association of NEFM expression with infiltration level of macrophage was replicated by correlation analysis, TIMER2.0 gene module and TISIDB; however, relationship of NEFM expression with infiltration levels of other TILs varied, possibly due to different numbers of available TCGA samples used for batch correction and differences in calculation methods. TIMER gene module (version 2.0) provides abundance of immune infiltration estimated by multiple immune deconvolution methods and adjusted by tumor purity, which is a major confounding factor in this analysis, and thus, the results are more accurate, with more reliable biological significance. Relationship of NEFM transcriptional expression with immunomodulators has implicated its involvement in regulating tumor immunology in BRCA. Firstly, macrophage markers IL6, CSF1R, CXCL12 were weak-to-strong positively associated with NEFM expression, which could reveal a potential role of NEFM transcriptional expression in regulating polarization of tumor-associated macrophage (TAM). In addition, NEFM transcriptional expression was positively associated with levels of T-cell exhaustion markers, specifically ADORA2A, VISTA and CCR4 [38,39]. Moreover, NEFM was positively associated with ectonucleotidases CD39 and CD73, novel checkpoint inhibitors that interfere with anti-tumor immune responses [40]. NEFM negatively correlated with PVR (CD155), an immune checkpoint on tumor cells and interacting with CD96, CD226, and TIGIT (T cell immune receptor with immunoglobulin and ITIM domains) on TILs to modulate immune function in tumor microenvironment [41]. In addition, NEFM significantly positively correlated with TGFB1, TGFBR1. Depending on the presence of other secreted factors and cell surface co-receptors, TGFB can either suppress adaptive immune responses (through induction and stabilization of Tregs and directly suppressing Th1 cell, Th2 cell and CD8 + T cell) or enhance adaptive responses (through induction of Th17 cell, Th9 cell and CD4 + CTL-like effector cell) [42]. The relationship of NEFM with TILs in BRCA may partially rely on   chemokines and chemokine receptors. More and more studies have shown that chemokines and chemokine receptors are closely related to the immunity of breast cancer. NEFM transcriptional expression was negatively associated with CCL7, CCL8 and CCL18, which would recruit monocytes to differentiate into tumor-associated macrophages (TAMs) at the tumor site, indicating that NEFM may cause decreased M2 macrophage infiltration [43]. Furthermore, NEFM transcriptional expression was positively associated with CCL14, CCL21, CXCL12, CXCL14, CCL28, CCR10 and CX3CR1. Notably, CX3CR1 promotes macrophage recruitment during mammary tumor formation. Macrophages are attracted to tumor sites expressing chemotactic factors such as CCL7, CCL8 and CXCL12 [43,44]. Additionally, CXCL12 promotes neutrophil infiltration to tumors. Moreover, CXCL12 is a potent attractant of dendritic cells (DCs); CCL21 recruits DCs and regulatory T cells (Tregs) [45]. CCL14 participates in the infiltration of the tumor by anti-cancer TILs. CXCL14 is responsible for immune cell recruitment and maturation and is critical to upregulating major histocompatibility complex class I expression on tumor cells. CCL28 activates CCR10 and causes B cell and T cell migration [46][47][48]. Epigenetic mechanisms, including DNA methylation, histone posttranslational modifications and chromatin structure regulation, are critical for tumor microenvironment (TME) (including immune cells) interaction. Emerging evidence supports that tumors commonly hijack various methylation mechanisms to escape immune surveillance. Recent studies have identified a strong connection between epigenetics and cytokine production in tumorigenesis [49,50]. Methyltransferases regulate production of interferons, cytokines and chemokines [51]. KMT3A (SETD2), a methyltransferase, is required for interferon pathway by catalyzing the methylation of STAT1, a key transcription factor of interferon response [52]. DNMT suppresses MHC-I expression on tumor cells [53]. The upregulation of PD-L1 on tumor cells likely results from selection pressure exerted by T cell immune response. Epigenetic mechanisms certainly contribute to upregulation of PDL1 [54]. Methylation regulators KMT6A (EZH2), MBD2, TET2 and demethylase KDM5B have been implicated in lymphocyte development [55]. DNMT3A controls fate decision of early effector CD8 + T cell. Loss of DNMT3A leads to ineffective repression of genes that are supposed to be silenced in effector cells, thus generating fewer effector cells [56]. These studies have revealed special relationship between TME-infiltrating immune cells and DNA methylation modification, beyond RNA degradation. NEFM DNA methylation significantly negatively correlated with various immunomodulators and most chemokines and receptors listed in TISIDB, which also contributed to decreased immune cells infiltration in BRCA. M1 macrophages are involved in normal Th1 immune responses, whereas M2 macrophages support survival and dissemination of cancer cells via secretion of various factors, including cytokines, chemokines and enzymes, which recruit Tregs intratumorally to suppress antitumor cytotoxicity [57]. We demonstrated that NEFM transcriptional expression was significantly associated with IL-17 signaling pathway and cytokine-cytokine receptor interaction by KEGG enrichment analysis. In breast cancer, we infer that NEFM transcriptional expression may affect survival partially through the decreased M2 macrophage infiltration and anti-tumor cytotoxicity of cytokine-cytokine receptor interaction and IL-17 signaling pathway. However, molecular mechanisms should be further investigated. One potential mechanism by which NEFM methylation associated with poor survival may be NEFM methylation inducing tumor immunosuppression depending on decreased TILs. Other mechanisms underlying relationship of NEFM DNA methylation with immune infiltration and poor prognosis in BRCA may include cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway, natural killer cellmediated cytotoxicity, primary immunodeficiency, T cell receptor signaling pathway, IL-17 signaling pathway and PD-L1 expression and PD-1 checkpoint pathway in cancer, which are significantly associated with NEFM DNA methylation by KEGG enrichment analysis.
It has been reported that interrogation of site-specific CpG sites may be another option for assessing immune infiltration in tumors and may possibly predict response to checkpoint inhibitors [58]. Since NEFM DNA methylation significantly negatively correlated with TILs and many immune pathways (especially PD-L1 expression and PD-1 checkpoint pathway in cancer) in BRCA, the six CpG sites within NEFM promoter associated with poor prognosis may serve as biomarkers for predicting immune infiltration in BRCA. Our results also suggest that demethylation of NEFM might be a strategy to improve the efficiency of immunotherapy. Thus, it is required to further explore the detailed mechanism and function of transcriptional expression/DNA methylation in regulating tumor microenvironment.

Conclusion
NEFM transcriptional expression positively correlates with favorable prognosis and increased levels of macrophage infiltration in BRCA. After adjusted by tumor purity, NEFM expression correlates with increased infiltration of CD8 + T cell, whereas decreased infiltration of B cell. NEFM DNA methylation correlates with poor prognosis and decreased immune infiltration of B cells, CD8 + and CD4 + T cells, macrophages, neutrophils and dendritic cells in BRCA. Moreover, NEFM expression/DNA methylation correlates with diverse immune markers and pathways in BRCA. Therefore, our study highlights potential clinical significance of NEFM transcriptional expression/DNA methylation in breast cancer and provides insight into a novel role of NEFM expression/DNA methylation in tumor immune infiltration.

Expression levels of NEFM in various types of cancer
The expression profiling of NEFM in various types of cancer was identified in the Oncomine database (https:// www. oncom ine. org/ resou rce/ login. html) [26]. The threshold was determined according to the following values: p-value of 0.001, fold change of 1.5, and gene ranking of all. We also analyzed NEFM expression in different types of cancer in TIMER database.

Prognosis assessment
The Kaplan-Meier plotter and univariable Cox proportional hazards regression models were used to estimate association between NEFM transcriptional expression or DNA methylation and median survival time. Multivariable Cox proportional hazards regression models were used to evaluate impacts of NEFM expression on OS in the presence of other known risk factors. NEFM associated with OS and RFS of BRCA patients was validated in Kaplan-Meier Plotter database (http:// kmplot. com/ analy sis/) [27] among 5,143 BRCA patients. Potential effects of NEFM expression on OS were evaluated in Pan-cancer RNA-seq in Kaplan-Meier plotter database. The association of NEFM protein expression with BRCA OS was analyzed by The Human Protein Atlas database (http:// www. prote inatl as. org/).

Immune infiltration
The data of immune infiltration in BRCA were downloaded from TIMER database. Relationship of NEFM transcriptional expression or DNA methylation with the abundance of immune infiltration was evaluated, including B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, M1 macrophages, M2 macrophages and dendritic cells by R package. We also analyzed the relationship of NEFM expression with the abundance of immune infiltration using gene modules in TIMER. Gene expression level normalized with tumor purity was displayed on the leftmost panel.
There were 28 TIL elements in TISIDB database. Relationship of NEFM transcriptional expression/ DNA methylation with abundance of TILs (including activated CD8 + T cell, activated CD4 + T cell, activated dendritic cell, activated B cell, macrophage, and neutrophil) was examined in TISIDB database.
In addition, the relationship of NEFM transcriptional expression or DNA methylation with immune gene markers was explored via Spearman's correlation. These immune gene markers included immunomodulators collected from Charoentong's study [28], chemokines and receptors based on TISIDB database [29]. The correlation scatter plots between NEFM transcriptional expression/DNA methylation and immune infiltration levels of immune cells, in BRCA, together with Spearman's correlation and estimated statistical significance, were described. The log2RSEM value of NEFM expression and log2β value of NEFM DNA methylation were used for x-axis, whereas related immune infiltration levels of immune cells for y-axis. Specific levels of gene markers were displayed with log2RSEM. TISIDB (http:// cis. hku. hk/ TISIDB/ index. php) is a web portal for tumor and immune system interaction, which integrates multiple heterogeneous data sets. The relative abundance of TILs as demonstrated by 28-gene immune-related signature from Charoentong's study was estimated by using gene set variation analysis (GSVA) based on gene expression profile in TISIDB database (Additional file 3) [29].
TIMER [30] is a comprehensive resource for systematic analysis of immune infiltrates across diverse cancer types (https:// cistr ome. shiny apps. io/ timer/). TIMER applies a deconvolution, a previously published statistical method, to infer relative abundance of tumorinfiltrating immune cells from gene expression profiles. TIMER database includes 10,897 samples across 32 cancer types from The Cancer Genome Atlas (TCGA) to estimate relative abundance of immune infiltration.

Enrichment analysis
Differentially expressed genes associated with NEFM expression and levels of NEFM DNA methylation were analyzed with DESeq2 R package. Volcano plots and heatmaps were presented. KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were performed with R package to identify pathways related to NEFM transcriptional expression/NEFM DNA methylation. A p-value of < 0.05 was considered as statistically.

Statistical analyses
Overall survival (OS) was calculated from the date of diagnosis to death due to any causes or to last followup. Recurrence-free survival (RFS) was calculated from the date of diagnosis to local relapse/recurrence or regional relapse/recurrence or death (all causes) whichever occurs first. The Kaplan-Meier method and log-rank test were used to estimate the relationship of NEFM transcriptional expression/DNA methylation with OS and RFS. The Fisher exact and Wilcoxon rank-sum tests were used, respectively, for categorical and continuous variables, to assess the relationship of NEFM expression levels and clinical or molecular characteristics. Multivariable Cox proportional hazards regression models were used to evaluate potential impact of NEFM expression on OS in the presence of other known risk factors. Student's t-test and multiple hypothesis correction (false discovery rate, FDR) were used to identify differences in genome-wide genes, methylation profiles between NEFM high and NEFM low groups. Spearman correlation analysis was performed to evaluate the relationship of NEFM methylation with transcriptional expression or other genes. A p value of less than 0.05 was considered statistically significant. All analyses were performed using R 3.6.1 software packages.