DLK1 Distinguishes Subsets of NF1-Associated Malignant Peripheral Nerve Sheath Tumors with Divergent Molecular Signatures

Abstract

Purpose: Malignant peripheral nerve sheath tumor (MPNST) is the leading cause of premature death among individuals with neurofibromatosis type 1 (NF1), and the transcriptional aberrations that precede malignant transformation and contribute to MPNST tumorigenesis remain poorly defined. Alterations involving CDKN2A and components of PRC2 have been implicated as early drivers of peripheral nerve sheath tumor (PNST) evolution, but these events do not occur in all MPNST. Accordingly, emerging data have begun to highlight the importance of molecular-based stratification to improve outcomes in patients with NF1-PNST.

Experimental design: In this study, we perform an integrated analysis of multiple, independent datasets obtained from human patients with NF1 to gain critical insights into PNST evolution and MPNST heterogeneity.

Results: We show that delta-like noncanonical Notch ligand 1 (DLK1) is significantly increased in MPNST and provide evidence that DLK1 overexpression may precede histologic changes consistent with malignancy. In complementary analyses, we find that serum levels of DLK1 are significantly higher in both mice and humans harboring MPNST compared with those without malignancy. Importantly, although DLK1 expression is increased in MPNST overall, through the integration of multiple, independent datasets, we demonstrate that divergent levels of DLK1 expression distinguish MPNST subsets characterized by unique molecular programs and potential therapeutic vulnerabilities. Specifically, we show that overexpression of DLK1 is associated with the reactivation of embryonic signatures, an immunosuppressive microenvironment, and a worse overall survival in patients with NF1-MPNST.

Conclusions: Collectively, our findings provide critical insights into MPNST tumorigenesis and support prospective studies evaluating the utility of DLK1 tissue and serum levels in augmenting diagnosis, risk assessment, and therapeutic stratification in the setting of NF1-PNST.

Figure 1.. DLK1 overexpression is associated with a worse overall survival in NF1-MPNST.

(A) Representative photomicrographs of human tumor sections across the PNST continuum immunohistochemically stained for DLK1. Magnification is denoted by 100 μm scale bars with inset high-power magnification as shown. The HALO cytonuclear mask used to quantify DLK1 staining is shown in the bottom panel. Blue cells indicate negative staining for DLK1, yellow, weak positive (1+), orange, moderate positive (2+) and red, strong positive (3+) for DLK1 staining. (B) Box and whisker plot depicting DLK1 positive cells as a percentage of total cells per field generated in R studio. Dots represent individual regions of interest (ROI). Error bars represent the 95% confidence interval. The center line represents the median. The box spans the 25^th to 75^th percentiles. Data beyond the whiskers are outliers and are plotted as individual points. PNF (n=9, ROI=48), ANNUBP (n=7, ROI=38) and MPNST (n=18, ROI=92) were analyzed by one-way ANOVA using Tukey’s multiple comparisons tests between groups. Tumors with greater than 50% of ROIs identified as statistical outliers by ROUT method (Q=0.1%) were considered outliers (PNF: n=2 of 9, ANF: n=1 of 7). Graph represents single experiment with outliers included. (C) Subject IU-26, a 17-year-old female with NF1, presented with a nodular, contrast-enhancing, FDG-PET avid, left pelvic mass (top). Histopathological consensus could not be reached, with the lesion (IU-26-A) being classified as either ANNUBP or LGMPNST following independent, expert reviews. Immunohistochemical staining of the resected lesion revealed diffuse positivity for DLK1 (bottom-right). Box and whisker plot generated in R studio comparing DLK1 positive cells as a percentage of total cells per field of IU-26-A to the other lesions classified as ANNUBP from our institutional cohort (bottom-left). Dots represent individual regions of interest (ROI). Error bars represent the 95% confidence interval. The center line represents the median. The box spans the 25^th to 75^th percentiles. Data beyond the whiskers are outliers and are plotted as individual points. P-value represents unpaired, two-tailed t-tests between groups. (D) Surface plots depicting DLK1 expression within the spatially profiled ANNUBP and the contiguous MPNST. Yellow corresponds to highest expression and dark purple to decreased expression as indicated. Grey represents no expression. (E) Kaplan Meier curve depicting overall survival of patients from our institutional cohort whose MPNST were characterized by either greater than (Positive, n=7) or less than (Negative, n=8) 1% of cells staining positively for DLK1 on IHC (347 vs. 1739 days, p=0.0443). P-value represents log-rank Mantel-Cox test. (F) Box and whisker plot generated in R studio depicting concentration (pg/mL) of DLK1 in the serum of Nf1/Cdkn2a Cre+ mice with MPNST (MPNST, n=19) and without MPNST (No Tumor, n=8). Dots represent individual mice. Error bars represent the 95% confidence interval. The center line represents the median. The box spans the 25^th to 75^th percentiles. Data beyond the whiskers are outliers and are plotted as individual points. P-value represents unpaired, two-tailed t-tests between groups. Graph depicts single experiment. Results were validated in a repeat experiment with a subset of samples. (G) Box and whisker plot generated in R studio depicting DLK1 concentration (ng/mL) in serum obtained from human NF1 patients with PNF (n=19) or MPNST (n=11). A single statistical outlier was identified by Grubbs’ method (α=0.0001) and was excluded from further analysis. Graph depicts single experiment with one outlier excluded (p=0.0452). Dots represent individual patient samples. Error bars represent the 95% confidence interval. The center line represents the median. The box spans the 25^th to 75^th percentiles. Data beyond the whiskers are outliers and are plotted as individual points. P-value represents unpaired, two-tailed t-tests between groups. Results were validated in a repeat experiment conducted in duplicate. ****= p-value ≤0.0001, ***=p-value ≤0.001 **= p-value≤0.01 *=p-value ≤ 0.05, ns=not significant.

Figure 2.. Overexpression of DLK1 is associated with the reactivation of embryonic gene signatures.

(A) Bar plot comparing DLK1 log₂ normalized counts between TCGA-MPNST (n=6) and GTEx normal nerve samples (n=278). Dots represent individual samples. Red dots indicate TCGA-MPNST with high DLK1 expression (DLK1_Hi), whereas blue dots indicate those with low (DLK1_Lo) expression. Error bars reflect the standard error of the mean (SEM). P-value represents unpaired, two-tailed t-tests between groups. (B) Principal component analysis (PCA) demonstrating global variation in gene expression between DLK1_Hi (red) and DLK1_Lo (blue) TCGA-MPNST samples based on principal components 1 and 2 (PC1 and PC2). (C) Lollipop plot generated in ShinyGO 0.80 depicting GO Biological Process 2023 enrichment of developmental and differentiation pathways by genes upregulated (adj. p-value≤0.05, log₂FC≥1) in DLK1_Hi compared to DLK1_Lo TCGA-MPNST. (D) Bar plot of DLK1 log₂ (CPM+4) normalized counts in MPNST obtained from patients cared for at our institution with high (DLK1_Hi, red) or low (DLK1_Lo, blue) expression of DLK1. Error bars reflect the SEM. P-value represents unpaired, two-tailed t-tests between groups. (E) Principal component analysis (PCA) demonstrating global variation in gene expression between DLK1_Hi (red) and DLK1_Lo (blue) MPNST samples from our institutional cohort based on principal components 1 and 2 (PC1 and PC2). (F) Heatmap depicting the expression pattern of the top 1000 variable genes sorted by 2-way, unsupervised hierarchical clustering. Columns represent individual genes, and their respective log₂ transformed, z-score normalized expression values across each sample in the data, with red corresponding to increased expression and blue to decreased expression according to the figure legend. DLK1 expression level is denoted by red and blue bars as indicated. Rows represent individual samples. (G) Lollipop plot generated in ShinyGO 0.80 depicting Go Biological Process 2023 enrichment of protein-coding genes (adj. p-value≤0.05, log₂FC≥1) upregulated in DLK1_Hi compared DLK1_Lo MPNST samples (n=10) from our institutional cohort.****= p-value ≤0.0001, ***=p-value ≤0.001 **= p-value≤0.01 *=p-value ≤ 0.05, ns=not significant.

Figure 3.. DLK1 expression is increased in human MPNST cells characterized by a stem-like phenotype.

(A) Bar plot of Notch1 transcript expression (RPKM) at indicated days of embryonic (E13.5, E17.5) and postnatal (P1, P5, P14, P24, and P60) development obtained from the Sciatic Nerve ATlas (SNAT). (B) Box and whisker plot of SOX11 mRNA expression comparing DLK1_Hi (orange) and DLK1_Lo (teal) MPNST from our institutional cohort generated in R studio. Dots represent individual samples. Error bars represent the 95% confidence interval. The center line represents the median. The box spans the 25^th to 75^th percentiles. Data beyond the whiskers are outliers and are plotted as individual points. P-value represents unpaired, two-tailed t-tests between groups. (C) Pairwise scatter plot for the co-expression of DLK1 and SOX11 across MPNST from our institutional cohort. Dots indicate individual samples with orange indicating DLK1_Hi and teal, DLK1_Lo, samples. The straight line corresponds to the linear regression fit. Pearson r= 0.80. (D) Western blot of DLK1 expression in human MPNST cell lines. GAPDH serves as a loading control. (E,F) Western blots of DLK1 expression in human MPNST cell lines S462 (E) and ST-8814 (F) grown in standard, adherent (Adh) or low-adhesion, spheroid-promoting (LA) conditions. Vinculin serves as the loading control. (G) Bar graph depicting DLK1 expression in human MPNST cell lines grown under spheroid-promoting (LA) conditions (LA) plates normalized to control (Adherent, Adh). Dots represent independent experiments. Error bars represent the standard error of the mean (SEM). P-value reflects one sample Wilcoxon test (theoretical mean=1). Graph represents pooled results of three independent cell lines. (H,I) Western blots of SOX11 expression in human MPNST cell lines S462 (H) and ST-8814 (I) grown in standard, adherent (Adh) or low-adhesion, spheroid promoting (LA) conditions. GAPDH and vinculin serve as loading controls. (J) Bar plot depicting the percent viability of a human MPNST cell line (S462) following siRNA-mediated depletion of DLK1 normalized to control. Dots represent independent experiments. Error bars represent the SEM. P-value reflects one sample Wilcoxon test (theoretical mean=100). Graph reflects pooled results from three independent experiments. (K) Representative images depicting a decrease in cell number of a human MPNST cell line (S462) following siRNA-mediated depletion of DLK1 compared to control. (L) Western blot depicting the expression of SOX11, DLK1, cleaved NOTCH1, and phospho-ERK (pERK) in a human MPNST cell line (S462) following siRNA-mediated depletion of DLK1 compared to control. Vinculin and GAPDH serve as loading controls. (M) Bar graph depicting DLK1 expression in a human MPNST cell line following siRNA-mediated depletion of DLK1 compared to control. Dots represent independent experiments. Error bars represent the SEM. P-value reflects one sample Wilcoxon test (theoretical mean=1). Graph represents pooled results of six independent experiments. (N) Bar graph depicting cleaved NOTCH1 expression in a human MPNST cell line following siRNA-mediated depletion of DLK1 compared to control. Dots represent independent experiments. Error bars represent the SEM. P-value reflects one sample Wilcoxon test (theoretical mean=1). Graph represents pooled results of three independent experiments. (O) Bar graph depicting SOX11 expression in human MPNST cell lines following siRNA-mediated depletion of DLK1 compared to control. Dots represent independent experiments. Error bars represent the SEM. P-value reflects one sample Wilcoxon test (theoretical mean=1). Graph represents pooled results of three independent cell lines. (P) Bar graph depicting phospho-ERK (pERK) expression in a human MPNST cell line following siRNA-mediated depletion of DLK1 compared to control. Dots represent independent experiments. Error bars represent the SEM. P-value reflects one sample Wilcoxon test (theoretical mean=1). Graph represents pooled results of three independent experiments. ****= p-value ≤0.0001, ***=p-value ≤0.001 **= p-value≤0.01 *=p-value ≤ 0.05, ns=not significant.

Figure 4.. DLK1_Hi MPNST are characterized by an immunosuppressive microenvironment.

(A) Gene set enrichment plot depicting enrichment of a stemness gene signature by MPNST from our institutional cohort with high (DLK1_Hi) and low (DLK1_Lo) DLK1 expression. Black vertical bars indicate the rank of genes comprising each signature. The black curve corresponds to the running enrichment score for the gene set. Normalized enrichment score (NES) and q-value for DLK1_Hi compared to DLK1_Lo samples are as shown. (B) Bar plot generated using EnrichR and depicting Go Biological Process 2023 enrichment of downregulated protein-coding genes (adj. p-value ≤ 0.05, log2FC≤−1) acquired through iDEP.96. (C) Cell type mapping panel generated in Omics Playground representing inferred cell types in DLK1_Hi and DLK1_Lo samples from our institutional cohort. ImmunoStates served as the reference dataset. (D) Heatmap showing suppression of macrophage associated genes in clusters (bold) characterized by elevated levels of DLK1 expression. Rows represent individual genes, and their respective log₂ transformed, z-score normalized expression values across each cluster with red corresponding to increased expression and blue to decreased expression according to the figure legend. (E-H) Surface plots depicting the expression of macrophage/monocyte markers CD68 (G) and CD14 (H). Surface plot depicting cluster designations is shown in (E). DLK1 expression as shown previously in Figure 1D serves as a reference (F). (I,J) Geneset enrichment plot depicting enrichment of a tumor-associated antigen (TAA) (I) and anti-PD1 response (J) gene signature by MPNST from our institutional cohort. Black vertical bars indicate the rank of genes comprising each signature. The black curve corresponds to the running enrichment score for the gene set. NES and q-value for DLK1_Hi compared to DLK1_Lo samples are as shown.

Figure 5.. Spatial transcriptomic profiling reveals aberrant expression of DLK1 preceding histopathological evidence of malignancy.

(A) Heatmap showing suppression of genes involved in the inflammatory response (GO:0006954) in the region of elevated DLK1 expression (cluster #5) of the MPNST-contiguous ANNUBP initially presented in Figure 1D. Rows represent individual genes, and their respective log₂ transformed, z-score normalized expression values across each cluster with red corresponding to increased expression and blue to decreased expression according to the figure legend. (B) Heatmap showing upregulation of genes involved in the nervous system development (GO:0007399) in the region of elevated DLK1 expression (cluster #5) of an MPNST-contiguous ANNUBP. Rows represent individual genes, and their respective log₂ transformed, z-score normalized expression values across each cluster with red corresponding to increased expression and blue to decreased expression according to the figure legend. (C-E) Surface plots depicting the expression of Schwann cell specific markers, S100B (D) and NES (E) in the MPNST-contiguous ANNUBP. Plot of DLK1 initially shown in Figure 1D serves as a reference (C). Arrow represents trajectory used to generate trajectory plots in K-O.(F-O) Surface and corresponding trajectory plots depicting the expression of genes associated with stem- and progenitor cell phenotypes including POSTN (F,K), GAP43 (G,L), NGFR (H,M), PTN (I,N) and PTPRZ1 (J,O). Yellow corresponds to increased expression and dark purple to decreased expression according to the figure legend. Grey indicates no expression. (P-S) Surface plots depicting the expression of genes involved in the inflammatory response including, IL-33 (P), CCL2 (Q), CCL5 (R), CCL14 (S) in the MPNST-contiguous ANNUBP. Yellow corresponds to increased expression and dark blue to decreased expression according to the figure legend. Grey indicates no expression. (T-W) Surface plots depicting the expression of T cell, CD2 (T), CD3E (U), CD8A (V), and macrophage/monocyte, CD68 (W) specific markers the MPNST-contiguous ANNUBP. Yellow corresponds to increased expression and dark blue to decreased expression according to the figure legend. Grey indicates no expression.

Figure 6.. Overexpression of DLK1 occurs independently of PRC2 loss.

(A) Schematic depicting the imprinted DLK1-DIO3 locus. Red represents genes expressed on the maternal allele and blue represents genes expressed on the paternal allele. Gray represents genes that are repressed. Filled circles represent methylated DMRs, whereas unfilled circles are unmethylated DMRs. (B,C) Scatter plots comparing DLK1 expression in TCGA Sarcoma samples (n=255) based on indicated gistic2-thresholded copy number alterations of PRC2-components SUZ12 (B) and EED (C) obtained through the Xena Browser database. “1” represents heterozygous loss (−) or gain (+), whereas “2” represents homozygous loss (−) or gain (+) and “0” indicates no copy number alteration (Neutral). Pearson correlation analysis conducted in GraphPad Prism revealed no significant correlation between copy number alterations of SUZ12 (Pearson r=−0.0247, 95% CI [−0.1471 to 0.09841], p=0.6943) or EED (Pearson r=−0.0611, 95% CI [−0.1826 to 0.06220], p=0.3311) and expression of DLK1. (D) Bar plot comparing inferred PRC2 complex activation in DLK1_Hi and DLK1_Lo TCGA-MPNST. Dots represent individual samples. Error bars reflect standard error of the mean (SEM). P-value represents unpaired, two-tailed t-tests between groups. (E) Heatmap depicting the methylation (beta value, β) pattern of TCGA-MPNST across indicated CpG sites (n=560) within the DLK1-DIO3 locus. Rows represent individual samples, and their respective β across each CpG site. Red corresponds to a higher, and blue to lower a β value according to the legend. DLK1 expression level is denoted by red and blue bars as indicated. Missing values are depicted in black (F) Bar plot comparing average methylation (beta value, β) across CpG sites within the DLK1-DIO3 locus between DLK1_Hi and DLK1_Lo TCGA-MPNST. Dots represent individual samples. Error bars reflect the SEM. P-value represents unpaired, two-tailed t-test between groups. (G) Bar plot comparing DLK1 mRNA expression across PRC2-deficient human MPNST cell lines. Dots represent individual samples. Error bars reflect the SEM. P-values represent one-way ANOVA using Tukey’s multiple comparisons tests between groups. Only significant comparisons are depicted on the graph. (H) Western blots of DLK1 expression in PRC2-deficient human MPNST cell lines JH-2-002 and NF90.8 following re-expression of HA-SUZ12 and restoration of H3K27me3. Yellow fluorescent protein with a nuclear localization signal (NLS-YFP) was used as an ectopic expression control. GAPDH serve as the loading control. (I) Bar plot comparing Suz12, Eed and Dlk1 mRNA expression in murine Dorsal root ganglia (DRG) NeuroSphere Cells (DNSCs) with (Nf1Arf-Cre) and without (Nf1-Cre, Nf1-GFP) malignant potential. Dots represent individual samples. Error bars reflect the SEM. P-values represent two-way ANOVA using Šídák's multiple comparisons test. ****= p-value ≤0.0001, ***=p-value ≤0.001 **= p-value≤0.01 *=p-value ≤ 0.05, ns=not significant.

Figure 7.. DLK1_Hi tumors are characterized by homozygous loss of CDKN2A and MTAP.

(A) Table depicting copy number alterations of CDKN2A and MTAP in MPNST (n=10) from our institutional cohort. −2 represents homozygous deletion and 0 represents neutral status as determined based on CNVkit thresholds. (B) Schematic adapted from cBioPortal depicting genes adjacent to CDKN2A and MTAP within the 9p21.3 locus. (C) Representative photomicrographs of serial MPNST sections from the institutional cohort (n=13, DLK1_Lo regions of interest (ROI)=62, DLK1_Hi ROI=78) immunohistochemically stained for DLK1 (top) and MTAP (bottom). Magnification is denoted by 500 μm scale bars with inset high-power magnification as shown. The HALO MTAP Stain 1 Mask is shown in the third row of panels and the cytonuclear mask used to quantify MTAP staining is shown in the bottom panel. (D) Bar plot depicting MTAP strong (3+) and moderate (2+) positive cells as a percentage of total cells per field. Dots represent individual ROIs. Error bars reflect standard error of the mean (SEM). P-value reflects unpaired, two-tailed students t-test between groups. ****= p-value ≤0.0001, ***=p-value ≤0.001 **= p-value≤0.01 *=p-value ≤ 0.05, ns=not significant.