Correlative study of epigenetic regulation of tumor microenvironment in spindle cell melanomas and cutaneous malignant peripheral nerve sheath tumors

The tumor microenvironment (TME) plays critical roles in tumor growth and progression, however key regulators of gene expression in the TME of cutaneous malignant peripheral nerve sheath tumor (C-MPNST) and spindle cell melanoma (SCM) have not been well elucidated. Herein, we investigate the epigenetic regulation of promoters and gene bodies and their effect on the TME composition of C-MPNSTs and SCMs. A cohort of 30 patients was analyzed using differential gene expression (DGE) and gene set enrichment analysis (GSEA) using the Nanostring platform. Methylation analysis was carried out utilizing an Infinium Methylation EPIC array targeting 866,562 methylation site (CpG) islands. DGE revealed overexpression of genes related to mast cells in the TME of SCMs, and a predominance of exhausted CD8+ T cells and macrophages in the TME of C-MPNSTs. Interestingly, we further observed promoter hypermethylation in key overexpressed genes and corresponding gene body hypomethylation. Analysis using ENCODE ChIP-sequencing data identified CTCF as the common transcription factor at the site of the hypomethylated probe. These findings support that the TME composition of C-MPNSTs and SCMs is at least partially independent on promoter methylation status, suggesting a possible relationship between gene body enhancers and expression of key TME genes in both entities.

Scientific RepoRtS | (2020) 10:12996 | https://doi.org/10.1038/s41598-020-69787-1 www.nature.com/scientificreports/ treatment is complete surgical excision with wide margins. SCMs are variants of melanoma that display fascicles of spindle cells arranged in a sheet-like manner with < 90% desmoplastic stroma in chronically sun-damaged skin. While conventional BRAF and NRAS mutations are not highly prevalent in SCMs, more recently, alterations in genes involved in aberrant activation of the mitogen-activated protein kinase (MAPK) and phosphoinositide 3-kinase (PI3K) signaling pathways have been identified with possibility for targeted therapies with small molecule inhibitors 8 . The tumor microenvironment (TME) establishes an intricate niche whereby tumor cells interact with neighboring stromal and extracellular matrix elements to create a collaborative setting that promotes tumor progression and an aptitude for metastasis 9,10 . Exploration of these entities has provided dynamic insights into the molecular classification and methylome signatures that aid in distinguishing between C-MPNSTs and SCMs, while further laying the foundation for potential future targeted therapies 11,12 . Recently, our group identified pivotal immune pathway perturbations in the TME detected in both C-MPNSTs and SCMs that included the janus kinase-signal transducer and activator of transcription (JAK/STAT) pathway, nuclear factor-κB (NF-κB), and CXCL12-CXCR4, in addition to CD274 (PD-L1) and CTLA4 overexpression 11 . Furthermore, subsequent interrogation of genomic methylation profiles identified a distinct methylome signature implicating the promoter region of the Branched-Chain Aminotransferase 1 (BCAT1) gene. Hypomethylation of the BCAT1 promoter region was seen in the C-MPNST group but not in SCMs, resulting in transcriptional activation and significant enrichment in gastric and breast cancer-related genes as emphasized on functional genomic pathway analysis 12 . Yet to date, the epigenetic regulation of the TME and its key players have not been well investigated in C-MPNSTs and SCMs.
Herein, we seek to investigate the epigenetic mechanisms governing the TME composition of C-MPNSTs and SCMs in an attempt to identify specific epigenetic events that would explain, at least partially, the difference in clinical behavior and provide actionable targets that warrant investigation in future studies. This was pursued implementing a tri-dimensional approach integrating differential gene expression (DGE), tandem genomewide methylation with differential methylation analysis, and gene set enrichment analysis (GSEA). Finally, we employed a methylation driven approach to estimate the TME composition in both entities (MethylCIBERSORT) and correlate its findings with the DGE and GSEA results.

Results
case selection. Our cohort is comprised of 30 patients with either C-MPNST (n = 15) or SCM (n = 15; mixed type) with representative photomicrographs shown in Fig. 1a-d. No significant difference in age (C-MPNST [mean = 56 years] vs. SCM [mean = 62 years]; p = 0.3, Student's t-test) was identified among the patient cohort, whereas a male predominance was noted in both groups (73% in C-MPNST vs. 80% in SCM). All cases were reviewed and confirmed by two different experienced dermatopathologists and soft tissue pathologists in order to reach a final diagnosis. The workup included ancillary immunohistochemical studies with SOX10 and S100 . Immunostains for SOX10 and S100 are shown for C-MPNST (e,f) and SCM (g,h), respectively (20× Objective Magnification). Note the diffuse pattern of both SOX10 (nuclear) and S100 (nuclear and cytoplasmic) immunoreactivity in SCM, and the patchy SOX10 (nuclear) and S100 (nuclear and cytoplasmic) immunoreactivity in C-MPNST. www.nature.com/scientificreports/ for both C-MPNST (Fig. 1e,f) and SCM (Fig. 1g,h), respectively, in addition to the identification of dermal and subcutaneous plexiform and/or diffuse neurofibromas for C-MPNST.

Differential gene expression in C-MPNST and SCM reveals distinct TME composition.
To investigate potential immunomodulatory differences in the TME composition of C-MPNSTs and SCMs we employed multiplex gene expression analysis using a customized 770-gene panel that encompasses 24 different immune cell types, common checkpoint inhibitors, cancer/testis antigens, and genes related to the adaptive and innate immune response. Select markers related to various cell types were scrutinized to detect changes in gene expression ( Table 1). Exploitation of 770 genes, including 60 validated marker genes defining 14 cell types, was applied for DGE and GSEA using the Nanostring platform. While both entities showed common activation of key pathways such as NF-κB and JAK/STAT pathways, important differences were noted at the gene level. DGE analysis of C-MPNST cases revealed overexpression of key genes corresponding to macrophages (CD68, CD163), exhausted CD8 + T cells (LAG3), and cytotoxic cells (KLRK1) (Fig. 2a). On the contrary, a predominance of mast cells was observed in the TME of SCMs as evidenced by overexpression of MS4A (MS4A1 and MS4A2) in the 60 marker genes (Fig. 2a) and STAT4 in the JAK/STAT pathway (Fig. 2b). Of note, gene transcripts associated with  www.nature.com/scientificreports/ activation of M1 macrophages (STAT1, STAT2), M2 macrophages (JAK1, JAK3) and CD274 were overexpressed in C-MPNSTs (Fig. 2b). These findings are concordant to the perturbed pathways identified in a previous study using network analysis 11 . Following identification of disparate immune signatures encountered in these two entities, we subsequently carried out GSEA to elucidate affected cancer immune pathways. Application of GSEA identified a significant upregulation of the "Rheumatoid arthritis" gene set in C-MPNST (Fig. 2c). Our analysis showed strong enrichment in specific genes related to pivotal growth factors (TGFB2, VEGFA), pro-inflammatory cytokines (CXCL5, CXCL6), lysosomal endopeptidases (CTSL), and specific human leukocyte antigen (HLA) alleles (HLA-DQA1, HLA-DQB1) in the C-MPNST group which have been reported to be intricately involved in the pathogenesis of rheumatoid arthritis. These pronounced differences in cellular composition provide evidence that the immune profiles of the TME in C-MPNSTs and SCMs are different.
Genome-wide methylation profiling reveals distinct genomic aberrations correlating gene expression and tMe composition. In order to investigate potential epigenetic events underlying the expression of the 60 above-mentioned genes, as well as those implicated in the JAK/STAT pathway, we applied an Infinium Methylation EPIC array targeting 866,562 methylation site (CpG) islands to assess genome-wide methylation levels. As expected, a positive correlation was seen between promoter hypomethylation and expression levels of genes catalogued among the 60 marker gene set and the JAK/STAT pathway. We further explored the methylation status of both promoter and body regions of the different genes of interest within this gene set. Surprisingly, key overexpressed genes from both the 60 above-mentioned gene set (LAG3 and CPA3) and from the JAK/STAT pathway (STAT4, JAK1, CD274) showed hypermethylation of the promoters and hypomethylation of gene bodies as highlighted by the color coded heatmap (Fig. 3a,b). We then assessed the methylation status of gene bodies and expression of genes involved in the JAK/STAT pathway and found a positive correlation (Pearson correlation coefficient: r = 0.37; p = 0.01). This observation raised the possibility whether expression of these genes is not regulated via the promoter regions but through epigenetic mechanisms pertaining to gene bodies.
To gain further mechanistic insight, we mapped the hypomethylated probe sequences to the Encyclopedia of DNA Elements (ENCODE) consortium data set to investigate protein-chromatin interactions from the transcription factor ChIP-sequencing (TF ChIP-seq) pipeline. Analysis of publically available data revealed binding sites for numerous transcription factors such as TAF1, POLR2A, ZNF263 and CTCF. Comparison of these gene lists using Venn diagrams identified CTCF as the common transcription factor at the site of the hypomethylated probe for all the above-mentioned overexpressed genes (Fig. 4a,b). These findings suggest an important role for CTCF in regulating expression of genes through gene bodies. www.nature.com/scientificreports/ Epigenetic deconvolution of methylation data corroborates TME differences identified on DGE analysis. Using MethylCIBERSORT we extended our analysis and identified marked alterations between the TME composition of C-MPNSTs and SCMs (Fig. 5a,b). In C-MPNSTs, MethylCIBERSORT detected a predominance of endothelial cells (p = 0.023, Kruskal-Wallis test), neutrophils (p = 0.031, Kruskal-Wallis test), and cancer-associated fibroblasts in the TME (p = 0.028, Kruskal-Wallis test), while a strong enrichment for CD8 + T cells (p = 0.0032, Kruskal-Wallis test) and CD14 + cells (p = 0.006, Kruskal-Wallis test) was seen in the SCM group set (Fig. 5c)

Discussion
In the current study, we pursued to explore the epigenetic events that could potentially explain the TME composition and differences in the distribution of various cell types among C-MPNSTs and SCMs to identify targets that would merit further investigation in future studies. Application of DGE analysis revealed distinct cell infiltration signatures within the TME of C-MPNSTs and SCMs. Results showed marked enrichment in exhausted CD8 + T cells and macrophages within the TME of C-MPNSTs as supported by LAG3, CD68 and CD163 overexpression, in contrast to a predominance of mast cells encountered in the TME of SCMs with overexpression of MS4A (MS4A1 and MS4A2), CPA3 and STAT4. Furthermore, GSEA identified a significant upregulation of the "Rheumatoid arthritis" gene set in C-MPNSTs. Interestingly, we observed promoter hypermethylation and gene body hypomethylation in a subset of overexpressed genes, including STAT4, LAG3, CPA3, CD274 and JAK1. These findings suggest a potential role implicating gene bodies as regulators for expression of these respective genes. DNA methylation plays pivotal roles in governing cellular processes such as regulation of gene transcription, chromatin structure remodeling, stability, tissue differentiation and X-chromosome inactivation [13][14][15][16] . The functionality of promoter hypermethylation has been well investigated in numerous studies across multiple decades, however the precise role of methylated gene bodies is poorly understood to date. Promoter hypermethylation results in the transcriptional silencing of genes, and aberrant methylation is now regarded as one of the hallmarks of cancer 13,17 . As such, investigators have pursued to mechanistically study whether a causal relationship exists between gene body methylation and gene expression [18][19][20][21] .
Chromatin structure and spatial positioning of the genome both have crucial implications on genome function 22 . CCCTC-Binding Factor (CTCF) is an essential DNA-binding protein involved in diverse functions such as mediating functional intra-and interchromosomal interactions to maintain genomic integrity, and restriction of the CTCF binding interface can result in alterations in gene expression 23,24 . Mechanistically, hypermethylation of genome CTCF binding sites or microdeletions present within these sites have been implicated in tumorigenesis and induction of oncogene expression 25,26 . CTCF has also been reported to pause RNA polymerase II, which in turn may have repercussions on splicing 27 . Our study represents an example of a possible link between gene body methylation status and gene expression, while additionally supporting the disparate nature of This was further supported by the fact that overlay with ChIP-seq analysis identified important binding sites on these hypomethylated probe sequences to pivotal transcription factor CTCF. Our analysis suggests that CTCF may represent a shared central player in the activation of the JAK/STAT pathway involved in tumorigenesis through differential expression of specific genes in both studied entities (C-MPNST and SCM). In fact, CTCF gene downregulation was shown to suppress cell proliferation, cell invasion and facilitate cell apoptosis in prostate cancer cell lines 28 . Whether inhibition of CTCF would have a similar effect in C-MPNST and SCM certainly warrants further investigation as it could represent an additional therapeutic modality to commonly used targeted therapies in C-MPNST and immunotherapy in SCM.
Interrogation of C-MPNSTs in our cohort revealed a preponderance of neutrophils, cancer-associated fibroblasts and endothelial cells. The intricate and dynamic disposition of the TME and signaling networks within this framework may play substantial roles in regulating responses to immunotherapeutic modalities. Neutrophils are highly abundant in the TME and studies have shown that the presence of a neutrophil-related gene signature has been associated with dismal prognosis across 39 types of cancer, whereas additional findings in mouse models have linked the presence of neutrophils to higher rates of metastatic tumor progression [29][30][31][32] . Moreover, cancer-associated fibroblasts are implicated in facilitating tumor progression and metastasis via cytokine secretion and epithelial-to-mesenchymal transition 33 . These lines of evidence suggest that the relative abundance of neutrophils and cancer-associated fibroblasts could potentially explain the more aggressive behavior and lack of immunotherapeutic response of C-MPNSTs in comparison to SCMs.
Infiltration of CD8 + T cells within tumor tissue is required to mediate cytotoxic effects against tumor cells. Tumor infiltrating lymphocytes can exert immunologic responses, whereas tumor regression has been shown in patients with metastatic melanoma following adoptive T cell transfer 34,35 . The predominance of CD8 + T cells encountered in SCMs, designating them as 'hot tumors' , could contribute to reported immunotherapeutic responses in these tumors notably with FDA approved therapies such as pembrolizumab or nivolumab.
In conclusion, our findings demonstrate that the TME in two morphologically overlapping entities is regulated by distinct epigenetic mechanisms which could explain the difference in clinical outcomes and prognosis between them. Furthermore, our deconvolution of the TME composition highlights specific cell populations such as neutrophils that warrants further investigation as potential targets for novel selective neutrophil exocytosis inhibitors (Nexinhibs) 36 in addition to currently available targeted therapies in the C-MPNST group. Finally, www.nature.com/scientificreports/ based on the identification of CTCF as a common underlying epigenetic modulator for different gene expression profiles in both entities, a centered therapeutic approach focusing on this particular transcription factor using selective CTCF inhibitors merits further investigation.

Methods case selection.
After institutional review board (IRB) approval at the 2 collaborating institutions, thirty cases including 15 cases of C-MPNSTs and 15 cases of SCMs were selected from patients diagnosed during the period from 2000 through 2015. All available slides for each case were reviewed by two pathologists certified in dermatopathology and bone and soft tissue pathology in order to confirm the diagnoses and select the most viable tumor areas for DNA and RNA extraction. All SCM cases were from resected tumors with an intra-epidermal in situ component with S100 (diffuse and strong nuclear and cytoplasmic), nestin, and SOX10 immunopositivity (all staining performed as part of the initial clinical workup), and moderate-to-severe solar elastosis. Criteria for diagnosing C-MPNST included the following: (1) lack of intra-epidermal in situ component and only dermal/ subcutaneous involvement; (2) lack of moderate-to-severe solar elastosis; (3) patchy weak-to-moderate nuclear and cytoplasmic S100 immunopositivity; and (4) the presence of a remnant neurofibroma from which the tumor arose (whenever applicable). In addition, any lesion exhibiting severe cytologic atypia with geographic necrosis and brisk mitotic activity was considered a C-MPNST. Using the following system, all cases had sufficient histologic and immunohistochemical findings that enabled an unequivocal diagnosis.
Tissue selection, nucleic acid isolation, DNA and RNA processing. For each specimen, ten 10 μm-thick formalin-fixed, paraffin-embedded (FFPE) sections were cut from a single representative block. Macrodissection was performed with a scalpel and focused on the areas of greatest tumor cell density to try to ensure that tumor cells accounted for at least 70% of the material analyzed for each case. DNA was recovered using the PicoPure DNA extraction kit (Thermo Fisher, Life Technologies), which enables recovery of genomic DNA from FFPE tissue 11 . DNA was then subjected to bead purification with the Sphere quality control kit (Thermo Fisher, Life Technologies Gene expression analysis. Multiplex gene expression analysis with 770 genes from 24 different immune cell types, common checkpoint inhibitors, cancer/testis antigens, and genes covering both the adaptive and innate immune response was performed using extracted RNA from patients' samples as previously described 11 . Differential gene expression (DGE) (log2 FC) and gene set enrichment analysis (GSEA) of 770 genes including 60 validated marker genes, defining 14 cell types in the tumor microenvironment, was performed using the Nanostring platform ( Table 1). The set of 60 genes was then validated against immunohistochemistry and flow cytometry in immune cells from tumor microenvironment across 23 types of TCGA solid tumors. DGE and GSEA, correlation between gene expression and methylation status of promoter and body regions was performed using R project with p < 0.05 and FDR < 0.01.

Methylation analysis.
The Infinium Methylation EPIC array (Illumina) was used to determine the DNA methylation status of 866,562 CpG sites, following the manufacturer's instructions. Minfi R package was used to process and analyze the methylation data 37 . Using minfi, the probes were quantile normalized and background adjusted. The resulting set of samples and probes was used for differentially methylated probes analysis. The differentially methylated probes finder function in minfi was used to identify differentially methylated probes for the 60 validated marker genes of interest in two comparisons: (1) SCM vs. C-MPNST, and (2) C-MPNST vs. SCM. Probes with statistical significance using Benjamini-Hochberg false discovery rate q cutoff of q < 0.05 were considered most significant, and corresponding heatmaps are shown. β-values for all 866,562 CpG sites tested were defined as the ratio of fluorescence intensity of the methylated probe over the overall intensity of probes. β < 0.2 indicated hypomethylation (blue); β > 0.8 indicated hypermethylation (red). For TME cell population assessment we used MethylCIBERSORT 38 -an R package for deconvolution that helps in estimating the tumor purity and cellular composition from DNA methylation data. Raw idats were processed and beta values were obtained. Signature genes were generated as described previously 38 . Beta matrix for tumor samples were loaded into CIBERSORT software along with the signature genes and the samples were deconvoluted to get the accurate estimate of the immune cell populations 38 . All statistical analysis and modeling was performed using the open-source software R (https ://www.r-proje ct.org/). All graphs and heatmaps were generated using the R package. Unsupervised hierarchical clustering was done with Euclidean measure for distance matrix and complete agglomeration method for clustering was used for unsupervised hierarchical clustering. ethics statement. Our research was conducted after approval by the ethics committee of the Office of Science and Research at NYU Langone Health (Study # i18-00637; Correspondence signed by director of the IRB: Helen Panageas) and followed the guidelines imposed by our institution. All methods listed in our research follow our institutional guidelines. A waiver of authorization for consent was approved by the review committee as there was no direct interaction with the patients and the retrospective nature of the study. The study involved