Spatial Gene-Expression Profiling Unveils Immuno-oncogenic Programs of NF1-Associated Peripheral Nerve Sheath Tumor Progression

Neurofibromatosis type 1 (NF1) is a cancer predisposition syndrome that affects approximately 1 in 3,000 individuals and results from mutations in the NF1 tumor-suppressor gene, which encodes neurofibromin, a GTPase activating protein for p21Ras (1, 2). Plexiform neurofibromas (PNF) are benign peripheral nerve sheath tumors (PNST) and are a hallmark of NF1, occurring in 50% of patients with NF1 (3, 4). PNF are heterogeneous neoplasms that despite similar histologic appearances can exhibit disparate biological behavior and growth kinetics across the lifespan of an individual (4, 5).

The lifetime risk of developing malignant peripheral nerve sheath tumors (MPNST) is approximately 8% to 16% in individuals with NF1 (6, 7). MPNST is a devastating form of sarcoma that is recalcitrant to chemotherapy and radiation and represents the leading cause of premature death in persons with NF1. Cure is only possible through wide marginal excision, which is often not feasible due to entanglement with vital structures. For patients presenting with unresectable disease, the median overall survival is approximately 5 months (8). Progression of PNF to MPNST often occurs with little to no early warning, and clinical indicators such as pain, change in consistency on palpation, nodularity, increased FDG-PET avidity, and rapid growth may only manifest once malignant transformation has already taken place. Thus, identifying key molecular alterations and biomarkers that precede malignant transformation is of critical importance for early detection and to inform risk-adapted care.

Malignant transformation of PNF often occurs through the development of atypical neurofibroma (ANF), which typically appear as distinct nodular lesions (DNL; ref. 9) and exhibit characteristic histopathologic features including atypia, increased cellularity, loss of neurofibroma architecture, and/or mitotic index >1/50 high-power fields (HPF) and <3/10 HPF (10). Recent consensus criteria have categorized these lesions as atypical neurofibromatous neoplasms of uncertain biological potential (ANNUBP) when two or more of these criteria are fulfilled, whereas neurofibromas with isolated atypia in the absence of other features are termed neurofibroma with atypia (NF with atypia; ref. 10). It remains unknown whether ANNUBP may exhibit other distinguishing characteristics from NF with atypia on a molecular level. Loss of the 9p21 locus harboring the CDKN2A tumor-suppressor gene has been identified as a key “second hit” in the development of ANF/ANNUBP and is present in the majority of ANF/ANNUBP and MPNST (11–16). Pathogenic variants in PRC2 complex genes and other diverse genetic events have been identified in MPNST (12, 13, 16, 17), yet are characteristically absent in PNF and ANF/ANNUBP (15). Loss of Cdkn2a and/or its alternate reading frame (Arf) in genetically engineered mice with conditional knockout of Nf1 is sufficient to engender the development of both ANF/ANNUBP and MPNST (18, 19), but the penetrance of malignant transformation is only 50% with loss of a single Arf allele (18). Similarly, not all patients with histopathologically defined ANF or DNL on imaging will develop MPNST, and the timing of malignant transformation remains unpredictable. Therefore, we reasoned that additional factors intrinsic to neoplastic Schwann cell precursors or within the tumor microenvironment may critically govern the malignant potential of ANF/ANNUBP precursor lesions.

Although genomic studies have advanced our knowledge of key driver events promoting neurofibroma progression, our understanding of transcriptional and microenvironmental dynamics orchestrating the progression and malignant transformation of PNF and ANF/ANNUBP precursor lesions remains limited. PNST are highly heterogeneous neoplasms, and evolution across multiple stages of the neurofibroma-to-MPNST continuum can be observed even within a single specimen. To overcome this challenge, we utilized a multiplexed, hybridization-based approach to quantify gene expression in spatially defined tissue regions. Through multiregional and longitudinal profiling of human PNST samples (n = 35), we identify alterations in signal transduction, immune response, and mitotic fidelity associated with disease progression and trace their dynamics across multiple stages of the neurofibroma-to-MPNST continuum within individual patients, providing insight into the nascent events associated with the evolution of MPNST precursor lesions. We further demonstrate that these core molecular signatures are retained in genetically engineered mice harboring conditional loss of Nf1 and Cdkn2a/Arf in neural crest-derived Schwann cell precursors that recapitulate the development of MPNST from existing PNF and ANNUBP. Finally, gene signature and pathway enrichment analyses comparing isolated lesions at discrete stages of PNST development as well as neurofibromas contiguous with adjacent MPNST provide preliminary insight into potential biomarkers of disease progression, that if validated prospectively in future studies, may help guide risk-adapted care.

Archived samples, associated clinical data, and imaging were collected under approval by the Indiana University Institutional Review Board (IRB Protocol No. 17332). Animal studies were conducted according to the guidelines established by the Indiana University Institutional Animal Care and Use Committee (IACUC Protocol No. 21009).

Samples were selected retrospectively by querying the Indiana University pathology archives with the search terms “neurofibromatosis type 1,” “ANNUBP,” “atypical neurofibroma,” “neurofibroma with atypia,” “plexiform neurofibroma,” and “MPNST.” Priority was given to cases of NF with atypia or ANNUBP that had additional samples from the same patient classified as PNF and/or MPNST. Several isolated PNST cases were also included in the analysis.

Hematoxylin and eosin (H&E)–stained tissue sections from each tumor were scanned using an Aperio ScanScope CS digital slide scanner at 20× magnification. A board-certified pathologist, with expertise in diagnosing NF1-associated nerve sheath tumors, digitally annotated tissue areas meeting consensus diagnostic criteria (10). Areas of hemorrhage and necrosis were excluded from the analysis.

Formalin-fixed paraffin-embedded (FFPE) tissues were microdissected from 5-μm-thick sections on freshly cut slides. Total RNA was extracted using an RNeasy FFPE Kit (Qiagen). RNA concentration and fragmentation index were quantified on an Agilent bioanalyzer. Gene-expression analysis was conducted on the nCounter Sprint Profiler (NanoString Technologies, Inc.) using the nCounter human PanCancer IO360 probe set (NanoString Technologies, Inc.). The fragmentation index (percent greater than 200 bp) was used to normalize the loading of 50-ng total input RNA in a final volume of 5 μL with the capture probe set. Probes and target transcripts were hybridized for 20 hours at 65°C per the manufacturer's recommendations. Hybridized samples were loaded into the nCounter cartridge and quantified on the nCounter Sprint Profiler. Quality control, normalization, differential expression, and pathway analysis were performed using nSolver Advanced Analysis Software (version 4.0, NanoString Technologies, Inc.) as further detailed below.

Adequate signal was observed across all samples in the data set as indicated by the heat map of raw counts (Supplementary Fig. S1; Supplementary Table S1). Forty probes below threshold with raw counts less than 25 (below the level of the background) were excluded from further analysis (Supplementary Table S2).

Data were normalized using both positive control probes and housekeeping genes. Prior to hybridization, nonnatural RNA sequence (positive control) probes were added to the samples at known concentrations to control for variations in pipetting and hybridization across samples. The read counts of each sample were then normalized by a correction factor computed by the geometric mean level of the control probes of across all samples divided by the levels of the control probes in individual samples.

To account for differences in RNA loading and quality, a set of housekeeping genes exhibiting low variance across all samples profiled was selected using the geNorm algorithm to form a normalization constant for each sample (Supplementary Tables S3 and S4). Housekeeping probes with significant sample-to-sample variability were excluded from the housekeeping normalization. Normalized and log₂-transformed normalized counts for each probe and sample used in the downstream differential gene-expression analysis are presented in Supplementary Tables S5 and S6, respectively.

Statistical analysis of differentially expressed genes and signature pathway scoring was performed using the Differential Expression and PathView modules of the nSolver Advanced Analysis Platform (version 4.0, NanoString Technologies, Inc.). Multiplicity-adjusted P values <0.05 were considered statistically significant. Results of differential gene expression (Supplementary Tables S7 and S8) and PathView analysis (Supplementary Table S9) used in downstream plotting and data visualization are presented in the Supplementary data. Sample metadata used for analysis are provided in Supplementary Table S10.

RNA-seq data from murine PNST arising in DhhCre mice harboring Nf1-floxed alleles and Ink4a/Arf heterozygosity, as described by Chaney and colleagues (19), were obtained from the Gene-Expression Omnibus (GEO) database (GSE148249). Raw count matrices with sample metadata deposited on GEO were reanalyzed using the Omics Playground v2.7.18 Platform. RNA-seq data from Nf1^flox/flox;PostnCre+ and Nf1^flox/flox;Arf ^flox/flox;PostnCre+ mice are available under GEO accession GSE213786 and GSE232707. Publicly available RNA sequencing data were also obtained from the TARGET TCGA GTEx database using the UCSC Xena Browser (RRID:SCR_018938). Differential gene-expression analysis was performed using the UCSC Xena Browser to compare gene-expression profiles between TCGA MPNST (n = 10) and GTEx Normal Nerve samples (n = 335). A volcano plot of the resulting data was generated using GraphPad Prism 9.5.1 software.

Principal component analysis (PCA) was conducted in IDEP.92 (20) using log₂-transformed, normalized counts as input under default settings.

The pathway heat map in Fig. 2C was generated using the nSolver Advanced Analysis Platform (version 4.0, NanoString Technologies, Inc). The gene-expression heat map in Supplementary Fig. S3 was generated using iDEP.92 (20). Log₂-transformed, z-score normalized counts from the top 100 most variable genes were clustered by correlation coefficient with average linkage. Heat maps in Fig. 5B and C were generated using Morpheus with log₂-transformed, z-score normalized counts from all genes clustered by correlation coefficient with average linkage. Heat maps in Supplementary Fig. S8B and Fig. 6C were generated using the Omics Playground v2.7.18 platform (21) and were generated using the ComplexHeatmap R/Bioconductor package (RRID:SCR_017270) on scaled log-expression values (z-score) using Euclidean distance and Ward linkage. The standard deviation was used to rank the genes for the reduced heat maps. Rows represent individual genes, and their respective expression log₂-normalized, z-score transformed expression values across each sample in the data. Columns represent individual samples, with annotation tracks below indicating sample IDs and above indicating the status of contiguity with MPNST.

Uniform manifold approximation and projection (UMAP) plots were generated by Omics Playground v2.7.18 (21). Clustering of features was computed using a normalized log-expression matrix (logCPM), with covariance as the distance metric. UMAP clustered gene sets were colored by standard deviation (sd.X). Labels reflect Hallmark gene set enrichment.

scRNA-seq data from four human MPNST samples with NF1 were downloaded from GEO (GSE179043; ref. 22). Data preprocessing and visualization were performed using Scanpy (RRID:SCR_018139) in Python. Low-quality cells with less than 200 genes, genes expressed in less than three cells, and outlier cells with less than the second percentile or greater than the 98th percentile of the number of genes detected, or greater than 20% mitochondrial counts were discarded. Filtered data from the four samples were concatenated. The top 4,000 highly variable genes were selected with the Seurat v3 method. An scVI model (23) was trained on the raw counts with sample and quality metrics (percentage of mitochondrial counts and total counts) as covariates. The latent representation and normalized expression values from the scVI model were extracted for downstream analysis to compute the neighborhood graph, UMAP coordinates, and Leiden clusters at a resolution of 0.5. After identifying cellular communities with similar transcriptional profiles, the decoupler package (24) was used to perform overrepresentation analysis (ORA) of cell type markers against canonical human cell type markers obtained from the PanglaoDB database (25). ORA estimates were converted to activation scores and ranked to identify the top five predicted cell types for each cluster. Final annotations were assigned manually following assessment of the enrichment results. To further distinguish malignant from nonmalignant cells, the infercnvpy package (https://github.com/icbi-lab/infercnvpy) was used to detect copy-number variations (CNV) from single-cell transcriptomic data. The average gene expression over genomic regions was compared between immune cell types (as reference cells) and query cells, using a window size of 25. CNV scores were computed for each cell based on the number and magnitude of CNV events detected and projected onto the UMAP coordinates for visual confirmation.

Five-micrometer-thick tissue sections were deparaffinized, hydrated, and transferred to 0.1M EDTA (pH 8.0) for antigen retrieval in a pressure cooker. Sections were then treated with 3% hydrogen peroxide for 10 minutes, rinsed and blocked with 5% normal goat serum in TBST (TBS buffer with 0.1% Tween-20), and incubated overnight at 4°C or 45 minutes at room temperature (CD31) with primary antibodies diluted in blocking buffer: CD31 (3528S, 1:100, Cell Signaling Technology), Survivin (MA5-15077, 1:200, Invitrogen), CENPF (PA5-84637, 1:200, Invitrogen), and H3K27me3 (9733S, 1:100, Cell Signaling Technology). Sections were then incubated with secondary antibodies for 1 hour at room temperature (goat anti-rabbit, ab205718, 1:1,000, Abcam or goat anti-mouse, ab205719 1:1,000, Abcam). VECTASTAIN DAB was applied for 10 minutes, and the reaction was terminated by rinsing in distilled water. Counterstaining was performed with modified Mayer's hematoxylin (Vector), and the sections were dehydrated, cleared, and coverslipped. Slide images were acquired on an Aperio ScanScope CS at 20× magnification. CD31⁺ staining blood vessels were counted manually on five randomly selected HPF per slide. Quantitative IHC analysis for BIRC5 and CENPF was conducted using the Cytonuclear IHC module of HALO Image Analysis software (version 2.0.5, Indica Labs). Cytonuclear analysis settings were optimized for each stain, and the intensity of nuclear staining was scored as negative (0, blue), weakly positive (1+, yellow), moderately positive (2+, orange), or strongly positive (3+, red). All positive cells (1–3+) were used for statistical analysis in GraphPad Prism software as described below.

Human MPNST cell lines JH-2-002 (26) and JH-2-103 (27) were obtained from the Johns Hopkins NF1 Biospecimen Repository. ST-8814 and S462 cells (28) were obtained from Dr. Andrew Tee (Cardiff University). Immortalized human normal (hTERT ipn02.3 2λ) and NF1-deficient neurofibroma Schwann cells (hTERT NF1 ipNF95.6; ref. 29) were obtained from Dr. Peggy Wallace. Cells were authenticated by STR analysis. Cells were cultured in either DMEM (NF90.8, ST88-14, ipn02.3 2 2λ and ipNF95.6) or DMEM/F12 (JH-2-002 and JH-2-103) media supplemented with 10% FBS (Harvest Midsci), 1% glutamine (Gibco), 1% penicillin/streptomycin (Lonza), and prophylactic 5 μg/mL Plasmocin (Invivogen). Trypsin-EDTA 0.05% (Gibco) was used to dissociate cells for passaging upon reaching confluence. Cultures were tested for mycoplasma and confirmed to be negative prior to experimentation.

Embryonic day 13.5 (E13.5) mice were extracted by cesarean section from freshly euthanized pregnant females. Dorsal root ganglia and spinal nerve roots were isolated under a stereomicroscope and treated with 20 mg/mL collagenase. After centrifugation and washing in PBS (No. 20012027, Gibco), the cell pellet was resuspended in DNSC complete media consisting of serum-free DMEM/F12 media (Gibco) supplemented with 0.2% heparin (StemCell), 30% glucose (Gibco), 7.5% sodium bicarbonate (Gibco), 1M HEPES (Gibco), 1% N2 supplement (Gibco), 1% glutamine (Gibco), 1% sodium pyruvate (Gibco), 1% penicillin/streptomycin (Lonza), 20 ng/mL epidermal growth factor (Sigma-Aldrich), 40 ng/mL basic fibroblast growth factor (PeproTech), 2% B27 (without vitamin A; Gibco), 40 μg/mL amphotericin B/Fungizone (Gibco), and prophylactic 5 μg/mL Plasmocin (Invivogen). Cells were cultured in low-adhesion plates for 7 to 10 days to form neurospheres, which were then transferred to fibronectin-coated plates for subsequent passaging.

Genotyping of Nf1^flox/flox; Arf ^flox/flox and Nf1^flox/flox;Cdkn2a^flox/flox DNSC was confirmed by PCR using the following primers:

P1: 5′-AATGTGAAATTGGTGTCGA GTAAGGTAACCAC-3′,

P2: 5′-TTAAGAGCATCTGCTGCTCTTAGAGGGAA-3′,

P3: 5′-TCAGACTGATTGTTGTACCTGAT GGTTGTACC-3′

P1: Forward (19243): 5′-ACT GCA GCC AGA CCA CTA GG-3′

Reverse (19244): 5′-AGC TCG GAG ATT GAG AAA GC-3′

P1: Forward (INK20): 5′-GTTTCCATTGCGAGGCTGCTCCGTAAGC-3′

Reverse (INK21): 5′-CTTTAGGGCGTTCCTTTCCCACTTCTGC-3′

P2: Forward (INK10): 5′-CCAAGTGTGCAAACCCAGGCTCC-3′

Reverse (INK11): 5′-TTGTTGGCCCAGGATGCCGACATC-3′

Nf1^flox/flox; Arf^flox/flox and Nf1^flox/flox;Cdkn2a^flox/flox DNSC were incubated with Cre recombinase expressing adenovirus (Ad5CMVCre-GFP High Titer, University of Iowa Viral Vector Core) at 1:1,000 dilution for 24 hours. The cells were then replenished with fresh DNSC complete media. Cre-mediated recombination of the floxed Nf1, Arf, and Cdkn2a alleles in the resulting Nf1^−/−Arf^−/− and Nf1^−/−Cdkn2a^−/− DNSC was confirmed by PCR and Western blot. Dissociation of DNSC for passaging was performed using 1X TrypLE express enzyme (Thermo Fisher Scientific). Prior to experimentation, cultures were tested for Mycoplasma and confirmed to be negative.

Human MPNST cell lines and murine DNSC were plated at a density of 5,000 cells per well in 96-well plates. Cells were allowed to adhere overnight and then treated with increasing concentrations of YM155 (MedChemExpress, HY-10194) or LQZ-7I (MedChemExpress, HY-136538). Cell viability was assessed using the CellTiter-Glo Assay (Promega) 48 to 72 hours after treatment. Endpoint luminescence was measured using a SynergyH4 plate reader with filters and settings as follows: 528/20 and hole filter sets, top read, 4-mm read height, gain 135, 0.5-second integration time.

Human MPNST cell lines were reverse transfected with 10 nmol/L siRNA against CENPF (Santa Cruz Biotechnology, sc-37563), BIRC5 (Santa Cruz Biotechnology, sc-29499), or scrambled control (Santa Cruz Biotechnology, sc-37007) using Lipofectamine 3000 transfection reagent (Invitrogen) according to the manufacturer's instructions. Cells were plated in the transfection complex at a density of 5,000 cells per well in 96-well plates. After 24 hours of transfection, cell lysates were collected to confirm protein knockdown by Western blot. Viability was assessed via CellTiter-Glo Assay 48 to 72 hours after transfection as described above. Luminescent values were normalized to the average of the siControl for each cell line.

Lysis buffer was prepared using cOmplete Mini Protease Inhibitor Cocktail (No. 11836153001, Roche), PhosSTOP Phosphatase Inhibitor Cocktail (No. 4906837001, Roche), and xTractor buffer (No. 635671, Takara Bio). Cell lysates were collected, and protein concentrations were determined using Pierce BCA Protein Assay Kit (No. 23227, Thermo Fisher Scientific). Isolated proteins were fractionated using NuPAGE 4% to 12% Bis–Tris Gels (Invitrogen, cat. No. NP0322BOX) and electro-transferred to PVDF membranes. Immunoblots were carried out using antibodies specific to CENPF (ab223847, 1:1,000 dilution, Abcam), Survivin (MA5-15077, 1:1,000, Invitrogen), and Vinculin (ab219649, 1:1,000 dilution, Abcam). After incubation with the primary antibody, an appropriate HRP-conjugated secondary antibody was used (anti-rabbit, No. NA934V, 1:2,000 dilution, GE Healthcare).

Statistical analyses were performed in R or using GraphPad Prism 9.5.1 software (GraphPad). One-way analysis of variance (ANOVA) with post hoc correction for multiple comparisons or unpaired t tests was used to evaluate statistically significant differences between groups as detailed in the accompanying figure legends. Adjusted P values ≤ 0.05 were considered statistically significant. For dose-response curves, the IC₅₀ values were calculated for each line using nonlinear regression in GraphPad Prism.

Raw RCC files, raw counts, and sample metadata pertaining to spatial gene-expression profiling of human PNST specimens described here are available from GEO under accession GSE239561. Bulk RNA-seq for comparing TCGA MPNST vs. GTEx normal nerve are in Supplementary Tables S11 and S12 and publicly available through the UCSC RNA-seq Toil Recompute Compendium (30) hosted through the UCSC Xena Browser. Bulk RNA-seq from Nf1^flox/flox;DhhCre and Ink4a/Arf ^+/−;Nf1^flox/flox;DhhCre murine PNST generated by Chaney and colleagues (19) are available from GEO under accession number GSE148249. Bulk RNA-seq data from Nf1^flox/flox;PostnCre and Nf1^flox/flox;Arf ^flox/flox;PostnCre + mice are available from GEO under accession numbers GSE213786 and GSE232707.

A cohort of 35 PNST samples was assembled by a retrospective review of the Indiana University Pathology archives. Thirty-two samples were obtained from subjects who met clinical diagnostic criteria for NF1, two additional samples were from a subject where the diagnosis of NF1 was strongly suspected pending genetic confirmation, and one case consisted of a sporadic neurofibroma in which the subject did not have NF1 (Fig. 1A). Some samples contained multiple PNST subtypes within a single lesion, and in select cases, multiple specimens were collected longitudinally from the same individual during clinically indicated procedures (Supplementary Table S13; Fig. 1B). NF with atypia were distinguished from ANNUBP based on published consensus diagnostic criteria (10). Two samples with a histopathologic diagnosis of neurofibroma (NF) and one with cellular NF were included in the hierarchical clustering and pathway analysis but excluded from statistical comparisons based on tumor histology due to their underrepresentation in the data set. Four out of six MPNST were negative for H3K27Me3 staining by IHC (Supplementary Fig. S2). Tissues were microdissected from selected regions of interest conforming to consensus diagnostic criteria (Fig. 1C) and subjected to RNA expression profiling across a high-content panel of 770 genes involved in tumor progression, microenvironment, and immune response (PanCancer IO360 Panel, NanoString Technologies, Inc.), thus allowing for spatially restricted, multiplexed quantification of gene expression from tumors across the neurofibroma to MPNST continuum (Fig. 1D).

Copy-number loss of the 9p21 chromosome cytoband encoding the CDKN2A tumor-suppressor locus has been identified in the majority of ANF/ANNUBP and MPNST (11–15). Seeking confirmation of these findings, we compared CDKN2A mRNA expression across the various PNST subtypes in our data set. CDKN2A exhibited variable expression in PNF, NF with atypia, and ANNUBP. Nonetheless, we observed reduced CDKN2A expression both in ANNUBP and MPNST relative to PNF (Fig. 2A). These findings are concordant with previously published data suggesting that increased Cdkn2a expression in PNF impedes the genesis of ANNUBP and malignant transformation by inducing senescence-mediated growth arrest in Schwann cells with biallelic Nf1 loss (18).

PNF, ANNUBP, and MPNST were largely distinguishable by PCA and unsupervised hierarchical clustering of the top variable features within the data set (Fig. 2B; Supplementary Fig. S3), although considerable transcriptional heterogeneity was observed. Notably, NF with atypia were highly heterogeneous and did not cluster discretely by unsupervised methods. Nonmalignant lesions from the same subject clustered more closely regardless of tissue diagnosis, whereas MPNST segregated distinctly. Taken together, these findings indicate that PNF, ANNUBP, and MPNST broadly possess distinctive global gene-expression signatures, but substantial transcriptional heterogeneity exists even within tumors of the same histopathologic subtype.

To further dissect the molecular heterogeneity across the PNST continuum, we implemented an unsupervised approach to computationally stratify samples agnostic of tumor histology based on pathway enrichment analysis of gene signatures related to cancer signal transduction, immune response, and the tumor microenvironment. Our analysis revealed three distinct subclusters predominated by PNF (cluster 1), ANNUBP (cluster 2), and MPNST (cluster 3), each characterized by unique transcriptional programs (Fig. 2C). As above, NF with atypia exhibited heterogeneous dispersion across all subclusters. CDKN2A mRNA expression was highest in cluster 1 (PNF predominant) and significantly reduced in both cluster 2 (ANNUBP predominant) and cluster 3 (MPNST predominant; Fig. 2D). Cluster 2 exhibited robust enrichment of pathways associated with immune system processes, antigen presentation, cytotoxicity, interferon signaling, hypoxia, and epigenetic regulation (Fig. 2E), while cluster 3 showed enhanced signatures of MAPK, PI3K/AKT, Notch, WNT, TGF-beta, JAK/STAT, Hedgehog, and NF kappa B pathway transcriptional output as well as enrichment of gene signatures associated with autophagy, angiogenesis, metabolic stress, matrix remodeling and metastasis, DNA damage repair, and cellular proliferation (Fig. 2F and G). Overall, these findings suggest a coordinated deregulation of both microenvironmental and tumor cell autonomous processes as neurofibromas progress across the PNST continuum.

To further examine the transcriptional alterations in signatures of antigen presentation and immune response across the neurofibroma-to-MPNST continuum, we performed differential gene-expression analysis, which revealed genes involved in antigen presentation, including HLA-DPA1 and HLA-DPB1, to be among the top upregulated genes in cluster 2 lesions (Fig. 3A–C). Conversely, multiple HLA-associated genes including HLA-DQA1 and HLA-DQB1 were found to be markedly downregulated in MPNST and cluster 3 lesions (Fig. 3D–F). Concordantly, we observed that ANNUBP are heavily infiltrated with CD4⁺FOXP3⁻ and CD8⁺FOXP3⁻ T cells, which were significantly diminished in MPNST and replaced by an increase in regulatory, FOXP3⁺ T cells (Fig. 3G and H). Signatures of cytotoxicity (Fig. 3I) in MPNST and cluster 3 lesions were also reduced and accompanied by increased expression of genes associated with immune exhaustion including CTLA4 (Fig. 3J).

Angiogenesis has been implicated in the malignant transformation of PNF in children with NF1 (31). In comparing differentially expressed genes in cluster 3 tumors and MPNST relative to cluster 1 lesions and PNF, respectively (Supplementary Fig. S4A and S4B), we identified genes associated with angiogenesis including MMP9 (32), FGF18, and VEGFA (33) to be among the top upregulated genes associated with malignant transformation (Supplementary Fig. S4C). IHC staining for CD31-positive blood vessels across PNF, ANNUBP, and MPNST in our data set was consistent with these findings, demonstrating a significant increase in tumor vascularity in MPNST (Supplementary Fig. S4D and S4E). Reanalysis of integrated, publicly available scRNA data sets (22) demonstrated that proangiogenic marker genes are highly expressed by a number of cell lineages in MPNST including both tumor cells, which highly express VEGFA, SERPINH1, HEY1, and FGFR1 as well as constituents of the tumor microenvironment including macrophages and dendritic cells, expressing MMP9 and fibroblasts, which abundantly express FGF18 (Supplementary Figs. S4F and S4G; S5). Collectively, these data provide additional insight regarding the angiogenic signatures characterizing MPNST identified in our data set.

Aberrant cellular proliferation and prosurvival signaling are hallmarks of malignant transformation (34). Genes associated with cellular proliferation and survival/apoptosis were also highly dysregulated in MPNST and cluster 3 tumors (Supplementary Fig. S6A–S6D). CENPF and BIRC5 were among the top differentially expressed genes comprising these signatures (Supplementary Fig. S6E; Fig. 4A and B). CENPF encodes centromere protein-F, a kinetochore-associated protein that safeguards chromosome segregation during mitosis (35). IHC staining confirmed elevated CENPF protein expression in ANNUBP and a marked increase in CENPF in MPNST as compared with both PNF and ANNUBP (Fig. 4C and D). BIRC5 encodes the Baculoviral Inhibitor of Apoptosis (IAP) Repeat Containing 5 protein, also known as survivin, which plays an essential role in preventing apoptosis via caspase inactivation. Survivin also participates in chromosomal alignment and segregation during mitosis as part of the chromosomal alignment and segregation during mitosis as part of the chromosome passenger complex (CPC; ref. 36). IHC analysis independently confirmed an increase in the percentage of BIRC5-positive cells in MPNST relative to PNF (Fig. 4E and F). We further validated these findings by comparing gene expression in normal nerve and MPNST from the Genotype-Tissue Expression (GTEx) and The Cancer Genome Atlas (TCGA) databases, respectively, which revealed BIRC5 and CENPF to be among the most highly upregulated genes in MPNST (Supplementary Fig. S7).

ANF/ANNUBP development and malignant transformation occur only in a small subset of individuals with NF1. Due to the limited availability of clinical samples, preclinical disease models that recapitulate key genetic driver events in specific cells of origin have played a critical role in understanding the biology and evolution of tumors along the neurofibroma-to-MPNST continuum. Conditional genetic inactivation of Nf1 and Arf in Schwann cell precursors driven by Periostin Cre (PostnCre) gives rise to tumors that recapitulate human ANNUBP and progress to MPNST with high penetrance, dependent on gene dosage (18). Similar results were observed in mice with conditional biallelic inactivation of Nf1 driven by Desert hedgehog Cre (DhhCre) and germline loss of Cdkn2a (Ink4a/Arf; ref. 19).

We reanalyzed publicly available RNA-seq data from tumors arising in Ink4a/Arf ^+/−;Nf1^flox/flox;DhhCre mice, as well as RNA-seq data generated from Nf1^flox/flox;Arf ^flox/flox;PostnCre mice developed in our laboratory. We evaluated the coherence of gene signature enrichments within PNST derived from these genetically engineered mouse models (GEMM; refs. 18, 19) and our patient-derived neurofibromas and MPNST described above. T-distributed stochastic neighbor embedding (t-SNE; Supplementary Fig. S8A) and unsupervised hierarchical clustering of top variable features (Supplementary Fig. S8B) demonstrated distinct transcriptional profiles of plexiform neurofibromas, paraspinal tumors (PST) reminiscent of human ANF, and GEM-PNST akin to human MPNST, arising in Ink4a/Arf ^+/−;Nf1^flox/flox;DhhCre mice, as reported by Chaney and colleagues (19). Enriched genes were queried against publicly available databases including Gene Ontology, KEGG, and Reactome. DNA repair emerged as the top enrichment in cluster S4 and was strongly upregulated in GEM-PNST and to a lesser extent in PST. This finding is consistent with results from clinical samples in our data set, which also identified enrichment of DNA-repair signatures in MPNST compared with PNF and ANF/ANNUBP (Supplementary Fig. S9A and S9B). Gene set enrichment analysis (GSEA) also demonstrated enhancement of cell-cycle/proliferation and apoptosis-related signatures in GEM-PNST (Supplementary Fig. S8C). We observed robust expression of Birc5 in GEM-PNST/ MPNST arising in both models (Supplementary Fig. S8D and S8E), consistent with our analysis of patient-derived samples. Cenpf expression demonstrated upregulation in MPNST compared with both PNF (P = 0.0647) and ANNUBP (P = 0.0783) in the PostnCre model, though not achieving significance. In the DhhCre model, Cenpf expression was significantly increased in GEM-PNST compared with PST (P = 0.0035), but was not significant compared with NF (P = 0.1126; Supplementary Fig. S8F and S8G).

Having observed concordant molecular features between human NF1-associated PNST samples and preclinical murine models, we aimed to investigate the biological significance of BIRC5 and CENPF as potential regulators of prosurvival and mitotic fidelity pathways in MPNST. We treated human MPNST cell lines (ST88–14, JH-2–002), primary murine MPNST precursors lacking Nf1 and either Arf (Nf1^−/−Arf ^−/−) or Cdkn2a (Nf1^−/−Cdkn2a^−/−), human wild-type (ipn02.3 2λ) and NF1-mutant Schwann cells (ipNF95.6) with pharmacologic inhibitors of BIRC5 (survivin), YM155 and LQZ-7I, and observed a significant reduction in cell viability with both drugs in human MPNST and murine MPNST precursor cell lines, with minimal effect on the viability of human wild-type and NF1-mutant Schwann cells (Fig. 4G; Supplementary Fig. S10), suggesting that BIRC5 inhibitors preferentially target MPNST cells. Likewise, siRNA-mediated depletion of CENPF (Fig. 4H and I) or BIRC5 (Supplementary Fig. S10B–S10E) impaired the viability of human MPNST cell lines. Together, these results demonstrate the utility of integrated, cross-species analysis to delineate putative targets with functional relevance to MPNST biology.

NF1-associated PNST exhibit substantial intratumoral heterogeneity. Lesions at multiple stages of the neurofibroma-to-MPNST continuum can coexist discretely within an individual or even within a single tumor itself. This heterogeneity poses unique challenges for diagnosis and treatment, as early identification of malignant transformation is crucial for preventing morbidity and mortality in persons with NF1. Needle biopsy specimens can be confounded by sampling bias, and furthermore, even neurofibromas that conform to uniform histopathologic criteria can exhibit diverse molecular features as shown above.

To further illustrate this concept, subject IU12-M, a 17-year-old female with NF1, presented with a large, FDG-avid, right upper extremity mass, which was confirmed to be high-grade MPNST on biopsy (Fig. 5A). Due to extensive involvement of the brachial plexus, complete surgical resection was not feasible. She received preoperative radiation followed by a partial resection and debulking procedure. Her tumor remained stable for 2 years following initial radiation and four cycles of chemotherapy with ifosfamide and doxorubicin. The patient subsequently returned to clinic with severe pain and MRI revealed significant tumor growth. A forequarter amputation was performed with resection of the shoulder and chest wall tumor. Multiregional tissue sampling revealed regions of tumor consistent with PNF, ANNUBP, and MPNST. The PNF component of the lesion (35_PNF) retained a molecular signature associated with predominantly benign lesions (cluster 1). However, the ANNUBP portion of the mass (34_ANNUBP) stratified molecularly in cluster 3 exhibiting enhanced expression of BIRC5, CENPF, MET, CDK6, and MKI67, which overlapped closely with the adjacent malignancy (35_MPNST; Fig. 5B).

Similarly, another individual, subject IU16-L, a 13-year-old female with NF1 presented to the clinic at age 13 with a large, FDG-PET avid, right, cervical, and occipital neck mass (Fig. 5C). Molecular profiling of multiple tumor regions from an initial needle biopsy and a subsequent surgical resection revealed components of the tumor spanning the neurofibroma-to-MPNST continuum. Notably, several histopathologically benign-appearing regions of tumor obtained from the initial needle biopsy (29_NF with atypia) and subsequent surgical resection (26_PNF) were characterized by a cluster 3-like molecular signature (Fig. 2A) and were transcriptionally similar to the adjacent malignancy (27_MPNST).

We conducted a detailed review of the clinical records from subjects in our data set to identify PNF and ANF/ANNUBP precursor lesions that were contiguous with adjacent MPNST. All microdissected samples were scrutinized by an expert clinical pathologist to confirm a lack of morphologic evidence of cross-contamination of malignant cells, thus allowing for the transcriptomic analysis of histopathologically benign lesions contiguous with MPNST. Interestingly, we found that none of the tumors in cluster 2 (ANNUBP predominant) were contiguous with MPNST. This may reflect a clinical bias to preemptively resect lesions histologically classified as ANNUBP due to their recognized potential for malignant transformation. Conversely, four out of the five (80%) histopathologically benign lesions computationally stratified in cluster 3 (Fig. 2A; 34_ANNUBP, 29_NF with atypia, 26_PNF, 13_PNF) were associated with the development of MPNST. The fifth lesion, 2_PNF, was obtained from subject IU13-B, an 8-year-old female who had innumerable NF1-associated neurofibromas involving her mediastinum and cervical–thoracic paraspinal region. Over a 4-year period, she underwent multiple biopsies and tumor resections, which identified lesions histopathologically consistent with PNF (2_PNF, 21_PNF) and ANNUBP (4_ANNUBP, 22_ANNUBP). 2_PNF corresponded to a DNL that was resected due to interval growth and clustered within the MPNST-predominant cluster 3 (Fig. 6A). Extensive internal and external histopathologic review revealed no evidence of malignant transformation within the lesion. Notably, after 5 years of follow-up, the patient developed a large posterior chest mass consistent with high-grade MPNST and ultimately succumbed to widespread metastatic disease following aggressive surgery, chemotherapy, and radiation.

Based on these findings, we conducted PCA and differential gene-expression profiling to compare precursor lesions contiguous with MPNST to discrete PNF and ANF/ANNUBPs not associated with malignant transformation. Lesions contiguous with MPNST were broadly distinguishable from isolated PNF, NF with atypia, and ANNUBP by PCA and unsupervised hierarchical clustering (Fig. 6B and C), indicating that they harbor distinct global gene-expression programs. Notably, the sole outlier, 2_PNF, was a nodular neurofibroma resected from subject IU13-B that clustered with MPNST-contiguous lesions (Fig. 6B and C, labeled), and this patient subsequently developed high-grade MPNST as described above.

Further differential gene-expression profiling (Fig. 6D) and GSEA (Fig. 6E–G) revealed upregulation of genes involved in TNF-alpha signaling, inflammatory responses, hypoxia, epithelial mesenchymal transition, and angiogenesis in neurofibromas contiguous with MPNST. Conversely, genes involved in progenitor cell differentiation, interferon-gamma signaling, and antigen presentation (HLA-DQA1, HLA-DQB1, HLA-DRB1) were downregulated (GO Biological Processes). IL6, CXCL8, and HLA-DQB1 emerged as the top three most significantly differentially expressed genes in neurofibromas contiguous with MPNST. To determine the ability of IL6, CXCL8, and HLA-DQB1 (Supplementary Fig. S11A–S11C) to distinguish neurofibromas contiguous with MPNST from discrete lesions not associated with malignant transformation, multiple logistic regression analysis was performed (Supplementary Fig. S11D) that revealed an area under the receiver-operating characteristic curve of 0.9778 (SE = 0.02246; 95% CI, 0.9338–1.000, P < 0.0001), suggesting that IL6, CXCL8, and HLA-DQB1 exhibited excellent discriminative power in distinguishing neurofibromas contiguous with MPNST in our data set (Tjur's R squared = 0.72, NPP = 95.00%, PPP = 88.89%, cutoff = 0.5; Supplementary Fig. S11E). Together, our findings suggest that neurofibromas contiguous with MPNST exhibit distinct alterations in core oncogenic and immune surveillance programs. Further investigation is needed to evaluate the prospective validity of these signatures as molecular diagnostic tools to complement conventional histopathologic diagnosis and guide risk-adapted care.

This study provides insights into the nascent transcriptional programs present in PNF and ANF/ANNUBP precursor lesions, both in isolation and concurrent with malignant transformation. Although PNF, ANNUBP, and MPNST exhibit distinct global gene-expression profiles, our findings demonstrate significant molecular heterogeneity even among tumors of a given histologic subtype. NF with atypia, in particular, were transcriptionally diverse and did not segregate discretely by PCA or unsupervised clustering methods. These observations suggest a potential role for molecular profiling to more accurately ascribe biological phenotypes within the context of consensus histopathologic classification schemes (10).

Differential gene expression and pathway-level enrichment analyses revealed deregulation of tumor cell autonomous factors and the tumor microenvironment at discrete stages of PNST development. ANNUBP exhibited gene signature alterations associated with antigen presentation and immune response, including dendritic cell function, immune cell chemotaxis, and T-cell activation. The ANNUBP predominant cluster, cluster 2, also exhibited robust enrichment of pathways associated with hypoxia and epigenetic regulation. Interestingly, previously published studies suggest that hypoxic tumor microenvironments drive epigenetic dysregulation and subsequent immunosuppression (37). Specifically, hypoxia impairs immune cell cytotoxicity and reduces interferon-gamma signaling, both of which are critical to the antitumor immune response and demonstrated enrichment in ANNUBP in this data set (37, 38). Conversely, compared with ANNUBPs, MPNSTs were characterized by further disruptions in apoptotic and angiogenic signature genes, as well as genes involved in mitotic processes and cell-cycle integrity.

Unsupervised, computationally based stratification identified alterations in core molecular signatures across the PNST continuum agnostic of tumor histology. These findings were cross-validated in two independent genetically engineered mouse models that spontaneously recapitulate the progression of MPNST from existing PNF and ANF/ANNUBP lesions driven by combined inactivation of Nf1 and Cdkn2a/Arf in Schwann cell precursors. Targeting key regulators of survival and mitotic fidelity pathways, including CENPF and BIRC5 (Survivin), potently disrupted the viability of human MPNST cell lines and primary murine Nf1-Ink4a/Arf-mutant Schwann cell precursors in vitro. These results are congruent with other published work implicating survivin as a potential therapeutic target in MPNST (39) and demonstrate the utility of integrated, cross-species omic analysis to identify functionally relevant therapeutic candidates in NF1-associated PNST and other orphan cancers where patient samples are limited.

Early detection of malignant transformation in patients with NF1 is crucial to prevent adverse outcomes. Intratumoral heterogeneity and sampling bias in needle biopsy specimens pose unique diagnostic challenges. Moreover, even neurofibromas that conform to uniform histopathologic criteria can exhibit diverse molecular phenotypes, growth kinetics, and clinical behavior. Neurofibromas contiguous with MPNST exhibited enrichment of genes involved in epithelial–mesenchymal transition, angiogenesis, proliferation, and hypoxic responses. Conversely, genes involved in antigen presentation were downregulated, suggesting a potential role of impaired antitumor immune responses in driving PNST evolution and transformation, which warrants further investigation. These results suggest that some of the molecular events driving the progression of neurofibroma to MPNST may be present before such changes manifest histopathologically.

Although this study represents a key step toward establishing molecular signatures that define PNST behavior, there are several key limitations. First, our study cohort was intentionally selected, focusing on subjects with ANF (NF with atypia and ANNUBP lesions) who also developed other tumors along the PNST continuum, including PNF and MPNST. Given the inherent selection bias, the incidence of MPNST in our study cohort is significantly higher than the 8% to 13% rate of malignant conversion in the general NF1 patient population. Furthermore, samples in this study were obtained retrospectively from subjects who were not followed as part of a formal natural history study. Therefore, it is possible that the prevalence of malignant transformation could be underestimated if subjects subsequently developed MPNST and sought care outside of our health care system. Although previous treatments including radiation, chemotherapy, and targeted therapies were documented based on a review of the electronic medical record (Supplementary Table S13), the sample size of this cohort is insufficiently powered to investigate the impact of these interventions on the tumor transcriptome. Genomic characterization of the tumor samples in our study was also limited primarily to a subset of MPNST that underwent sequencing for clinical indications, and information regarding germline pathogenetic variants and/or microdeletions of NF1 was not available. Finally, although the probe set used here included a high-content panel of genes with established relevance to cancer biology, tumor progression, and microenvironment, it does not provide comprehensive whole-transcriptome coverage. Thus, other key transcriptional alterations may not have been captured in these data.

In addition to transcriptomics, combining multiple technologies including but not limited to ATAC-seq, methylation profiling, and cell-free DNA (cfDNA) as predictive biomarkers may likely enhance our ability to accurately forecast the biological behavior and stratify risk of malignant transformation of neurofibroma precursor lesions. For instance, cutaneous neurofibroma exhibit site-specific methylation changes in genes associated with RAS-dependent signaling that may contribute to their distinct growth kinetics in comparison with PNF (40). Recently, global methylation profiling revealed that most ANF exhibit an epigenetic profile similar to benign PNST (41). However, similar to our study, a subset of NF with atypia and ANNUBP clustered distinctly with MPNST, suggesting that histopathologic diagnosis alone may be insufficient to ascribe the potential risk of malignant transformation associated with these lesions. Ultra-low-pass whole-gene sequencing of cfDNA has shown promise in distinguishing PNF from MPNST based on copy-number alterations and has demonstrated potential for monitoring treatment responses and early detection of relapsed disease (42). Recently, multiomic analysis conducted by the international Genomics of MPNST (GeM) consortium has further identified specific somatic copy-number aberrations in both MPNST samples and cfDNA that function as a surrogate for H3K27 trimethylation loss and predict prognosis (16).

Prospective validation of biomarkers predictive of malignant transformation is limited by the practice of preemptively resecting DNLs/ANNUBP at many NF centers including our own. DNLs/ANNUBP that are inoperable, however, offer a valuable opportunity to elucidate molecular features that prognosticate the biological behavior of these lesions. Our findings suggest that suppression of immune surveillance and antigen presentation, along with upregulation of transcriptional programs associated with proliferation, survival, and developmental regression may indicate the emergence of malignant transformation. If validated prospectively, targeted gene-expression analysis could augment conventional histopathologic criteria to improve risk stratification and clinical decision-making for surveillance and early intervention in the treatment and prevention of MPNST.

C.A. Pratilas reports grants from Neurofibromatosis Therapeutic Acceleration Program during the conduct of the study, grants from Novartis and Kura Oncology, and personal fees from Day One Therapeutics outside the submitted work; in addition, C.A. Pratilas has U.S. Patent No. 7,812,143 B2 issued and a patent for U.S. Provisional App. No. 63184422 pending. S.D. Rhodes reports grants from NIH/NINDS (K08-NS128266-02) to S.D. Rhodes, grants from Neurofibromatosis Therapeutic Acceleration Program (2004757180) to S.D. Rhodes, grants from NIH/NINDS (R01-NS128025-01) to D.W. Clapp, and grants from NIH/NCI (U54-CA196519-07) to D.W. Clapp during the conduct of the study. No disclosures were reported by the other authors.

D.K. Mitchell: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. B. Burgess: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. E.E. White: Data curation, formal analysis, validation, investigation, visualization, methodology. A.E. Smith: Data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. E.A. Sierra Potchanant: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–review and editing. H. Mang: Data curation, methodology. B.E. Hickey: Data curation, formal analysis, methodology, writing–review and editing. Q. Lu: Data curation, validation. S. Qian: Data curation, validation. W. Bessler: Data curation, visualization. X. Li: Data curation, validation. L. Jiang: Data curation, validation. K. Brewster: Investigation. C. Temm: Data curation, validation. A. Horvai: Data curation, formal analysis, validation, investigation, visualization, writing–review and editing. E.A. Albright: Resources, data curation, validation. M.L. Fishel: Resources, validation, methodology, writing–review and editing. C.A. Pratilas: Resources, visualization, writing–review and editing. S.P. Angus: Writing–review and editing. D.W. Clapp: Resources, supervision, funding acquisition, project administration, writing–review and editing. S.D. Rhodes: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration.

This work was supported by a K08 Mentored Clinical Scientist Research Career Development Award to S.D. Rhodes by the National Institute of Neurological Disorders and Stroke (NIH/NINDS, K08-NS128266-02), R01-NS128025-01 (to D.W. Clapp) and a Developmental and Hyperactive RAS Tumor SPORE funded by the National Cancer Institute (NIH/NCI, U54-CA196519-07) led by D.W. Clapp. S.D. Rhodes is also supported by the Francis S. Collins Scholars Program in Neurofibromatosis Clinical and Translational Research funded by the Neurofibromatosis Therapeutic Acceleration Program (2004757180) and startup funding from the Department of Pediatrics at the Indiana University School of Medicine, the Indiana University Simon Comprehensive Cancer Center, and the Indiana University Precision Health Initiative. We thank Katie Jackson and Heather Daniel for administrative support.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.