Experiments involving animals were conducted under protocols approved by the University of Wisconsin–Madison Institutional Animal Care and Use Committee in accordance with published National Institutes of Health and United States Department of Agriculture guidelines. We have previously published on the creation and genomic validation of our three distinct NF1 porcine models, including the one utilized in this study (24). The present study used tissue samples (two tumors and an adjacent normal region) collected from a 3.83-year-old male NF1 pig harboring an NF1 exon a31 excision producing an alternatively spliced NF1 allele, whose regulation is associated with disease severity (24, 26). Development of the masses was first noted when the animal was about 11 months of age; the masses increased in size with age until the animal was euthanized for necropsy and tissue collection. The pig was bred and housed at the UW Swine Research and Teaching Center, a closed-herd, specific-pathogen-free (SPF) facility. The pig was transported to an on-campus facility on the day of the necropsy. The pig was initially sedated with intramuscular administration of a cocktail of Telazol^® (3.3 mg/kg; Zoetis Inc., Kalamazoo, MI, USA), xylazine (1.6 mg/kg; Bimeda, Inc., Oakbrook Terrace, IL, USA), and ketamine (1.6 mg/kg; Zoetis Inc.). An additional dose of ketamine (2 mg/kg) was given intramuscularly before euthanasia with intracardiac administration of Fatal-Plus^® Euthanasia Solution (Vortech Pharmaceutical Ltd, Dearborn, MI, USA) according to the manufacturer’s recommendations. The necropsy was led by two board-certified veterinary pathologists.

Caudodorsal subcutaneous masses were identified on the right and left sides of the thorax and resected, providing two samples we have called Mass 1 and Mass 2 ( Figure 1 ). A non-tumor, normal adjacent sample was collected to serve as the control. Upon gross examination, the masses were primarily located in the subcutaneous tissue and extended into the overlying dermis. Like human neurofibromas, both masses appeared light-yellow to tan in color, semi-firm, and variably demarcated from surrounding tissue. Mass 1 was located on the right caudodorsal thorax, 88 cm caudal to the base of the right pinna, and 12 cm to the right of midline, while Mass 2 was located left caudodorsal thorax, at approximately the same level as the right caudodorsal mass, 10.5 cm left of midline. Once incised around the masses, Mass 1 was approximately 6 x 5 x 3.2 cm, while Mass 2 was 11.7 x 11.7 x 3.2 cm. Each tumor tissue sample was taken from the center of the mass and divided into two segments for histological and immunohistochemical validation and the scRNA-seq analysis.

Tissue sections were processed and embedded using an adapted protocol (Abcam, Cambridge, UK) (27). Briefly, following fixation with 10% neutral buffered formalin, tissues placed in processing/embedding cassettes were dehydrated using increasing ethanol concentrations. This was followed by clearing with three changes of xylene for 20 minutes each. Xylene was then exchanged with two paraffin changes (60 minutes each) to allow adequate tissue infiltration of the wax and then embedded. The embedded tissues were cut into 5 µm sections for the conventional hematoxylin and eosin (H&E) and IHC staining. The IHC protocol was adapted from R&D Systems (Minneapolis, MN, USA) (28). Heat-induced antigen retrieval was performed using a citrate buffer at pH 7.4 (29). The following antibodies were used to classify the tumors: S100 (1:4000; Agilent Dako, Santa Clara, CA, USA; GA50461-2), CD34 (1:400; Bioss, Woburn, MA, USA; bs-8996R), Sox10 (1:100; Abcam, Cambridge, UK; ab227680), CD117 (1:100, Cell Signaling; 37805S), GFAP (1:200, Abcam; ab16997), AIF1 (1:500, LSBio; LS-B2645), and SMA (1:250, LSBio; LS-B3933). Briefly, the primary antibodies were incubated overnight at 4°C while the secondary antibodies, including anti-rabbit, anti-goat, and anti-mouse (1:250, Abcam; ab6728, ab6885, ab6721, respectively), were incubated on slides for 30 minutes at room temperature. The targets of interest were detected using a DAB chromogen substrate kit (Cell Signaling; #8059) and counterstained using Mayer’s Hematoxylin for 30 seconds (Abcam; ab220365).

The samples for the scRNA-seq analysis were promptly placed into 1.5 mL conical tubes on ice until enzymatically digestion could be started within the hour to obtain a single cell suspension (30). Briefly, each tissue sample was minced into <4 mm pieces and digested with a Tumor Dissociation Kit (Cat #130-095-929; Miltenyi Biotec, Gaithersburg, MD, USA) according to the manufacturer’s instructions. The resulting suspensions were further digested using collagenase II to obtain sufficient cell counts for scRNA-seq (>50,000 cells/sample) (31). Libraries were constructed according to the Chromium Single Cell 3’ Reagent Kit v3.1 user guide (10x Genomics, Pleasanton, CA) by the UW Gene Expression Center. Briefly, cell concentration and cell viability of the single-cell suspension were quantified on the Countess II (Thermo Fisher Scientific, Waltham, MA) using 0.4% Trypan Blue (Invitrogen, Carlsbad, CA). The appropriate volume of cells was loaded onto the Chromium Single Cell Chip G (10X) required for yielding a cell recovery of approximately 50,000 cells. After completing the Chromium run, the gel beads-in-emulsion (GEMs) were transferred to emulsion-safe strip tubes for GEM-RT using an Eppendorf MasterCycler Pro thermocycler (Eppendorf, Hamburg, Germany). Following RT, GEMs were broken, and the pooled single-cell cDNA was amplified. Post-cDNA amplified product was purified using SPRIselect (Beckman Coulter, Brea, CA) and quantified on a Bioanalyzer 2100 (Agilent, Santa Clara, CA) using the High Sensitivity DNA kit (Agilent, Santa Clara, CA). Full-length cDNA was fragmented and used to generate cDNA libraries according to the standard 10X Genomics workflow. These libraries were sequenced on a NovaSeq 6000 Sequencing System (Illumina Inc., San Diego, CA) with paired-end 150bp sequencing by the UW DNA Sequencing Facility.

The scRNA-Seq data were analyzed by the UW Bioinformatics Resource Center. Experiment data were demultiplexed using the Cell Ranger Single Cell Software Suite and aligned to Sscrofa11.1, Barcode counting, UMI counting, and gene expression estimation for each sample according to the 10x Genomics documentation (32). The gene expression estimates from each sample were then aggregated using Cellranger (cellranger aggr) to compare experimental groups with normalized sequencing depth and expression data.

Analysis was performed with R Statistical Software (v4.1.3) and Seurat (v4.2.0) (33). Each sample (Mass 1, Mass 2, and Normal Adjacent Tissue) was individually processed for quality control before integration. Seurat objects were created with thresholds: min.cells = 3, min.features = 200. We used ENSEMBL to identify human orthologs of porcine genes, and manually used VGNC to identify remaining orthologs. Next, data was subsetted based on a threshold for 18 mitochondrial genes set at <10% representation and <1% representation of Hemoglobin Beta ( Supplementary Table 1 ). The threshold for unique features per cell was >200 and <6000 to avoid low-quality cells and doublets. Seurat objects for each dataset were merged then split to create lists that can be transformed and normalized using SCTransform (34).

Initial clustering divided the cell population into 14 individual clusters, which were automatically identified with ScType ( Supplementary Table 2 ) using a modified marker list generated by converting human gene symbols to orthologous Sus scrofa gene symbols (35). Lists of cell types were curated and included markers for all possible cell types within these tissues, and subsequent cell type scores were generated using the full scRNA-seq dataset and differentially expressed genes (DEGs) (35). Automatically scored and identified cell types were then verified manually using commonly known cell type markers ( Supplementary Table 1 ). Violin plots for cell types are depicted using the raw RNA counts for genes of interest.

Analysis of gene expression was performed using Seurat v4.2.0 ( Supplementary Document ), and DEGs were identified using the normalized RNA counts, or the SCT assay. To prepare the Seurat object to run differential expression testing on the SCT assay, we used PrepSCTFindMarkers(). For bulk comparisons of tumor vs. normal, we used RenameIdent() to merge both Mass 1 and Mass 2. Populations were subsetted for comparison using Subset(). Pairwise comparisons between clusters was made using FindMarkers() and specific identification of clusters. DEGs in one cluster compared to all other clusters were done using FindAllMarkers(). The default Wilcoxon Rank Sum test was used for analysis. DEGs were defined as being expressed in >20% of the cluster, with |logFC| >0.25, and adjusted p-value <0.05 ( Supplementary Tables 3 , 4 ).

DEG lists were populated with Sus scrofa Ensembl IDs. First, we automatically converted Ensembl IDs using “biomaRt” datasets “sscrofa_gene_ensembl” and “hsapiens_gene_ensembl”. This method was only successful in converting ~70% of Sus scrofa Ensembl IDs to human gene symbols. The remaining IDs were then converted manually, and the list was uploaded to R to be used for GSEA. The Fast Gene Set Enrichment Analysis (fgsea) package was used to analyze enriched pathways from the GO: Gene Ontology gene sets (36). The top 15 enriched pathways were visualized by adjusted p-value initially, and only pathways with an adjusted p-value <0.05 were used for further analysis. Enrichment plots were visualized using plotEnrichment() of specified pathways ( Supplementary Table 5 ).

Extensive phenotyping of the focal study animal (and over 30 others from our three NF1 porcine models) will be published elsewhere. Here, we focus on the validation of the two resected masses. The tissues were first stained with hematoxylin and eosin (H&E) and Masson’s Trichrome to examine tumor architecture and collagen deposition. Both masses demonstrated typical histologic features of diffuse neurofibroma involving subcutaneous adipose tissue and showing an admixture of Schwann cells with wavy, slightly hyperchromatic nuclei, admixed with fibroblasts and fibroblast-like cells, small axons, interspersed immune cells, thick or wavy collagen bundles, and blood vessels ( Figure 2 ). The presence of lesional collagen was highlighted by Masson’s Trichrome staining ( Figures 2 , 3 ). Hypercellular Schwann cell-rich zones, Verocay bodies, and nuclear palisading indicative of schwannomas were not observed (37). The adjacent normal tissue was primarily composed of dense collagen bundles surrounding areas of adipose. Mononuclear, interspersed immune cells morphologically resembling mast cells were present in low abundance, but we could not objectively confirm this cell type due to the absence of a compatible cKIT antibody. Mass 2 showed similar findings of diffuse neurofibroma; but some areas also showed S100⁺ Schwann cell-rich nodules surrounded by concentric rings of collagen, immune cells, and other cell types ( Figure 3 ). Neither case showed convincing precursor plexiform neurofibroma or could be traced to the originating nerve. Therefore, the overall morphologic findings strongly favored diffuse neurofibromas, most likely of plexiform origin (diffuse plexiform neurofibromas, dpNFs) rather than of cutaneous origin.

Immunohistochemistry revealed that both tumors displayed strong cytoplasmic and nuclear positivity for two proteins used in neurofibroma diagnosis, S100 positivity in Schwann cells and CD34, which highlighted fibroblasts, fibroblast-like cells, and a subset of Schwann cells (38). The adjacent normal tissue showed little to no positivity for these markers ( Figure 2 ). Tumors also focally stained positive for neural cell origin markers GFAP and SOX10, particularly in areas of hypercellularity (data not shown) (38). Immunostaining for smooth muscle actin (SMA), a myofibroblastic and smooth muscle marker, produced little to no positivity except for a few larger vessels that can also be visualized by red Masson’s Trichrome staining ( Figure 2 ) (39, 40). To summarize, we observed the typical histologic components of neurofibromas, including mixed cellularity, spindle-shaped Schwann cells, neurofibroma-associated fibroblasts, coarse and wire-like collagen bundles, blood vessels and moderate immune infiltration, including mononuclear cells histologically similar to mast cells. The lesional spindle cells were S100⁺, SOX10⁺, GFAP⁺, and CD34⁺, in keeping with known human neurofibroma markers. With this data, we can conclude that the two masses (Mass 1 and Mass 2) accurately recapitulate the histology of neurofibromas from human NF1 patients by being positive for necessary neurofibroma markers and reproducing many of the essential features common to neurofibromas.

After single-cell dissociation, library preparation, sequencing, data integration and normalization, we were left with 31,211 individual cell transcriptomes with an average of 2,980 unique genes per cell. To then identify shared clusters and associated cell types, all samples were subject to dimensional reduction and non-supervised data clustering as described in Materials and Methods. The UMAP (Uniform Manifold Approximation and Projection) from this analysis initially outputted 14 distinct clusters of cells that were later combined into 12 clusters ( Figure 4A ). Using both automatic and manual cell typing, we identified Schwann cells (S100A1⁺, SOX10⁺, NCAM1⁺ ), fibroblasts (COL1A1⁺, MMP2⁺, DCN⁺ ), myeloid cells (FCGR3A⁺, CSF1R⁺, SIRPA⁺ ), vascular smooth muscle cells/pericytes (ACTA2⁺, TAGLN⁺, DES⁺ ), endothelial cells (vWF⁺, PECAM1⁺, CDH5⁺ ), T cells (GPR183⁺ , CD3D⁺, CD3E⁺ ), and neutrophils (TGM3⁺, CXCL8⁺, S100A12⁺ ) ( Figures 4B , 4C ). Histological analysis of human neurofibromas and those from mouse models have established that the cells of origin of neurofibromas are from the Schwann cell lineage and that supporting cells within the tumors include fibroblasts, vascular smooth muscle cells, immune cells, and endothelial cells (11, 22, 41, 42). Schwann cells were found almost exclusively in our tumor samples when compared to the normal adjacent sample and made up 4.37% and 0.14% of all cells, respectively ( Figure 4D ). Because this analysis and further analyses did not identify noticeable differences between Mass 1 and Mass 2 sample data, data from these samples were pooled for further analysis. We did not specifically identify perineurial or endoneurial fibroblasts, but we suspect these cell types are among one or more fibroblast clusters.

The tumor microenvironment was enriched for myeloid cells compared to normal tissue (9.30% vs. 3.84%) ( Figure 4D ); the 2,334 myeloid transcriptomes represented 7.48% of all cell types. This cluster was defined as myeloid lineage cells based on the pan-expression of SIRPA, CSF1R, and FCGR3A ( Figure 5C ). SIRPA and CSF1R were recently confirmed to be pan-myeloid lineage markers in the characterization of the porcine immune system (43). However, we observed that markers for more precise subtypes of myeloid cells were not uniform in expression, suggesting that more specific myeloid subtypes could be elucidated. To identify these specific cell types within the myeloid cluster, we subsetted this population from the full dataset and performed principal component analysis (PCA) to generate a new UMAP, identifying eight unique clusters ( Figures 5A, B ). We then generated lists of the top significant DEGs for manual identification and GSEA. Based on the overall higher expression of several MHC Class II molecules, such as SLA-DRA, and the lower relative expression of CD163, we identified Clusters 6, 7, and 8, as dendritic cells (Dendritic Cells A, Dendritic Cells B, and Dendritic Cells C, respectively) ( Figure 5D ). The remaining clusters had higher relative expression of CD14 and CD163 and were identified broadly as macrophages. However, we found that three out of the five macrophage clusters (Clusters 2, 3, and 4) were best classified as M2 macrophages (M2 Macrophage A, M2 Macrophage B, and M2 Macrophage C, respectively), characterized by a higher ratio of CD206 (M2) to CD86 (M1) and in some cases, higher expression of macrophage scavenger receptors such as MSR1 ( Figure 5E ). Cluster 1 was the only macrophage cluster with higher levels of CD86 and SLA-DRA and thus was classified as M1 macrophages ( Figures 5D, E ). Furthermore, M1 macrophages were enriched for GO terms “Response to Cytokine” and “Immune Response,” primarily through the higher expression of CXCL8, IL1B, CCL5, and CXCL16 ( Supplementary Table 4 ). Surprisingly, Cluster 5 was defined by high expression of COL1A1, COL1A2, and DCN while still expressing pan-myeloid markers ( Figure 5F ). While this cluster was the smallest within the myeloid cell population, it was nearly exclusively represented by tumor-derived macrophages (91.53%). Fibrotic macrophages have been characterized in patients with pulmonary fibrosis and correlated with increased mortality and fibrosis (44). Additionally, a recent study using a mouse model of fibrotic scar formation in the heart has shown that monocyte-derived macrophages directly contributed to COL1A1 deposition (45). This cluster might represent a rare population of macrophages with a profound impact on a tumor microenvironment centered around extracellular matrix (ECM) reconstruction.

M2 macrophages represent a critical population in the tumor microenvironment and primarily promote immunosuppressive and pro-fibrotic pathways often found in advanced cancers (46). M2 Macrophage A was enriched for CCL2, MSR1, CD209, CD206, and SCARB2 ( Supplementary Table 4 ). GSEA revealed that this cluster, representing 21% of the total myeloid population, was enriched for GO terms “Innate Immune Response” and “Endocytosis,” suggesting a phagocytic M2 macrophage population. M2 Macrophage B was enriched for various chemokines, including CCL2, CXCL12, CCL19, and CXCL14 ( Supplementary Table 4 ). In addition to these cytokines, GSEA revealed that M2 Macrophage B was negatively associated with “Immune response” and “Cell Activation,” further suggesting a strong immunosuppressive phenotype (47–51). M2 Macrophage C had a higher relative expression of complement system components, such as C1QA/B, and was also negatively associated with “Cellular Response to Biotic Stimulus.” Evidence suggests that complement proteins on macrophages can drive immunosuppression and local tissue remodeling and indicate a poor prognosis (52, 53). Dendritic Cells A is the second largest myeloid cluster and notably was overrepresented by cells derived from tumors (89.69% vs. 83.60% expected, χ^2 = 18.3, p <0.00005). Cells in this cluster expressed immune checkpoint modulators IDO1 and CD274 (PD-L1), thus suppressing T cell activation and providing critical immunosuppressive action in the microenvironment ( Figure 5G ) (54, 55). Furthermore, we observed that IDO1/PD-L1 expression in the cluster was nearly exclusive to tumor-derived cells ( Figures 5H, I ). Taken together, we provide evidence for M2 macrophage polarization, IDO1⁺ /PD-L1⁺ dendritic cells, and a novel population of macrophages that combine to provide a strong immunosuppressive tumor microenvironment, further supporting the advanced stage of these neurofibromas.

Collagen accounts for up to 70% of lipid-free dry-weight in human neurofibromas (56, 57). A recent study of the matrisome of a human cutaneous neurofibroma found that neurofibroma fibroblasts preferentially deposited pro-tumorigenic collagens rather than classical pro-fibrogenic collagens (58). Our current study found that fibroblast-like cells represented the majority (over 70%) of cells in all samples and were positive for canonical fibroblast markers COL1A1, MMP2, and DCN ( Figures 4B, C ). We identified 24,532 transcriptomes derived from fibroblasts and fibroblast-like cells in tumors and normal adjacent tissue, giving us detailed insight into this population. All six clusters were found in all three samples at varying levels of representation. Five clusters (Fibroblasts 1, 3, 4, 5, and 6) were primarily found in the tumors, while Fibroblasts 2 were found mainly in the normal adjacent tissue ( Figure 6A ). However, before analyzing each cluster individually, we separated this population by sample to examine differences between normal adjacent-derived fibroblasts and tumor-derived fibroblasts ( Figure 6B ). While all fibroblasts expressed COL1A1 at variable levels, fibroblasts in the adjacent normal tissue were significantly enriched for canonical fibroblast markers CTHRC1, CCN2, COL1A1, and COL1A2 ( Supplementary Table 3 ) (59–61). In contrast, fibroblasts derived from the neurofibromas were significantly enriched for CXCL2, IL6, CCL19, and CCL2 ( Figure 6C , Supplementary Table 3 ). Top upregulated markers in fibroblasts from the neurofibroma were mostly pro-inflammatory and corresponded to markers of CAFs, while those of normal tissue were associated with resident dermal fibroblasts ( Supplementary Table 3 ) (59, 62–64). GSEA of fibroblast-related DEGs revealed that tumor-derived fibroblasts were significantly enriched in processes “Cytoplasmic Translation,” “Peptide Biosynthetic Process,” and “Defense Response” ( Supplementary Table 5 ). Neurofibroma-associated fibroblasts were negatively associated with “Collagen Fibril Organization,” “External Encapsulating Structure Organization,” and “Collagen Metabolic Process” ( Supplementary Table 5 ). The matrisome from the porcine neurofibromas recapitulates that of human neurofibromas described by Brosseau and colleagues (58).

After bulk comparisons of fibroblasts in the tumors and normal adjacent tissue, we explored each fibroblast cluster individually to elucidate fibroblast diversity. While this large group of cells can be placed under the umbrella of fibroblasts, initial data integration and analysis split the group into six distinct clusters ( Figure 7A ). As previously noted, we found that Fibroblast 2 was the only cluster represented by a majority of normal-derived fibroblasts, while Fibroblasts 1, 3, 4, 5, and 6 were represented by a majority of tumor-derived fibroblasts. Functionally, normal-associated fibroblast markers were enriched in Fibroblast 2, while cancer-associated markers were enriched in the remaining five clusters ( Figures 6D, E ). Because collagens are critical components of the neurofibroma matrisome, we examined the expression levels of all detected collagens in each fibroblast cluster. Fibroblasts 2 had significantly higher expression of most collagens analyzed in this study, highlighting the importance of that cluster in ECM construction ( Figure 7B ). Fibroblasts 2 was enriched for PCOLCE, a pro-collagen that enhances reconstruction of the ECM, as well as some of the most abundant pro-fibrotic collagens found in the body, COL1A1 and COL1A2 ( Figure 7B ) (65). However, additional noteworthy collagens were higher in the other clusters, perhaps highlighting a specialized role for each group. COL5A3, which promotes Schwann cell adhesion and neurite growth (66, 67), was significantly higher in Fibroblasts 6. COL4A2 was also enriched in Fibroblasts 6; Type IV collagens are localized to basement membranes and surround Schwann cells to promote axonal growth (66). Finally, COL15A1, a critical component in nerve regeneration (68, 69), was enriched in Fibroblasts 4 ( Figure 7B ). With this data, we hypothesize that Fibroblasts 2 is the primary driver of ECM deposition, while Fibroblasts 4 and 6 may play a more specialized role in supporting processes of nerve regeneration. Notably, 87.8% of Fibroblasts 4 and 74.7% of Fibroblasts 6 were found in the neurofibromas ( Figure 6B ).

We next examined the expression of several major ECM glycoproteins and ECM regulatory proteins (58). With a few exceptions, we found that Fibroblasts 2 was enriched for most major ECM glycoproteins, including MGP, CCN2, and SPARC ( Figure 7C ) (70, 71). Notable exceptions include FBLN5 and LAMA4, most highly expressed in Fibroblasts 4 ( Figure 7C ). FBLN5 is a critical factor in promoting Schwann cell proliferation and myelination, while LAMA4 works in combination with COL15A1 to promote nerve regeneration (72). This further suggests that Fibroblasts 4 may represent a specialized group of fibroblasts within the neurofibroma microenvironment that promotes nerve regeneration and maintenance processes. Fibroblasts 2 was enriched for most of the queried ECM regulatory protein genes, while Fibroblasts 5 was enriched for several cathepsins ( Figure 7D ). Cathepsins are proteases involved in the maintenance of ECM by promoting the turnover of extracellular proteins (73). Taken together, we hypothesize that Fibroblasts 2 and 5 primarily promote ECM reconstruction and maintenance, while Fibroblasts 4 and 6 are more specialized for nerve regeneration and maintenance.

GSEA was performed on each fibroblast cluster to examine critical pathways potentially driving tumor progression or tumor microenvironments. To elucidate roles unique to each fibroblast cluster in the context of the microenvironment, we compared each to all other cells in the dataset, including non-fibroblast cells. All significant upregulated and downregulated DEGs with a minimum log-fold change of 0.25 were then analyzed using the fgsea package against GO-Biological Processes (GOBP) gene sets. This analysis revealed distinct enriched pathways for each cluster of fibroblasts, further highlighting each cluster’s unique role in the tumor microenvironment beyond matrisome function.

Fibroblasts 1, a cluster composed primarily of tumor fibroblasts (86.65%), was highly enriched for ribosomal L and S mRNAs ( Supplementary Table 3 ) and other genes associated with translation and ribosome biogenesis ( Figure 7F ) (74). The significance of increased protein translation is difficult to pinpoint, but heightened translation is a critical component of tumorigenic and differentiating cells, suggesting this is a highly active population (75, 76). Dysregulation of mRNA translation is a major feature of neoplasia and stem cell maintenance (77–79). Plasticity is often a consequence of cells adapting to the tumor microenvironment, especially hypoxia, nutrient limitation, or cancer therapeutics, and may reflect an important role of change in translation in phenotypic switching of cancer cells (75, 80). Additionally, ribosomal proteins can have functions outside of translation, including induction of cyclins and NFKβ (81–83). Notable overexpressed genes for Fibroblasts 1 include IL-6, CSF-1, CXCL12, ITIH4, and CXCL14, all factors that promote inflammatory responses and commonly found in CAFs ( Figure 7E , Supplementary Table 3 ) (84–89). CSF1 is one of the ligands for CSF1R on macrophages, and upon binding, CSF1 promotes M2 polarization of macrophages. This suggests that Fibroblasts 1 may be partly responsible for macrophage class switching and immune evasion.

Fibroblasts 2, the only cluster predominantly composed of fibroblasts from adjacent normal tissue, was highly enriched for dozens of collagens and ECM regulatory proteins, as noted in the previous section. Consequently, the top three significant non-redundant GOBP processes enriched for this cluster were “Collagen Fibril Organization,” “External Encapsulating Structure Organization,” and “Animal Organ Morphogenesis” ( Figure 7F ). Many of the top overexpressed genes include the previously discussed collagens, but also cell-cell communication molecules CCN2, SPARC, and MGP ( Figure 7E ) (71, 90, 91). Interestingly, this cluster also overexpressed THY-1 and SPON2. While THY-1 and SPON2 can be expressed in normal fibroblasts and smooth muscle cells, they are also associated with synaptogenesis, neuron outgrowth, and tumorigenesis. They are often found in other nervous system cell types (92–96).

Spatially, Fibroblasts 3 was surrounded by the other five fibroblast clusters and was generally non-descriptive in enriched processes and genes. The top 3 enriched pathways relate to negative regulation of transcription, while other significant non-redundant processes include the “Rhythmic Process” and “Pattern Specification Process” ( Figure 7F ). A decreased rate of transcriptional activity is associated with senescence, differentiation, apoptosis, angiogenesis, and neoplastic transformation (97, 98). Although chromatin exposure and high transcription rates are found in embryonic stem cells, upon initiation of differentiation, cells acquire markers for transcriptional repression. A recent meta-analysis of human lung cancer samples in The Cancer Genome Atlas found that genes represented in the “Rhythmic Process” were associated with immune activation, circadian rhythm, and carcinogenic pathways (99). With this information in combination with enriched processes of pattern specification, regionalization, and various organ development processes, we hypothesize that this population represents a dedifferentiated group of fibroblast-like cells that are more involved in the tumor microenvironment than the GOBP terms suggest.

Fibroblasts 4 represented the highest proportion of tumor-derived cells within a single cluster at 87.75% ( Figure 6B ). This cluster was unique because many top upregulated genes are associated with neural development or nerve regeneration. These genes included NEGR1, MYOC, EPHA3, COL15A1, PTN, FBLN5, PDGFRβ, and THY-1 ( Figure 7E , Supplementary Table 3 ) (72, 94, 100–105). In addition to these genes, this cluster was enriched for both “Cell-Matrix Adhesion” and “Regulation of Cell Junction Assembly” ( Figure 7F ). These two processes are essential in the development of signaling networks and neural regeneration. Specifically, EPHA3 and NEGR1 are primarily expressed in the adult brain and promote neuronal synaptogenesis and neurite outgrowth (106, 107). Loss of NEGR1 in development is associated with decreased numbers of synapses and dendritic length, resulting in anxiety and depression-like behaviors in mice (100). Furthermore, FBLN5 and PDGFRβ are found to be upregulated immediately following sciatic nerve crushing in rats. It has been recently reported that FBLN5 promotes the proliferation of Schwann cells in vitro (72, 105). We hypothesize that this cluster is profoundly important in promoting a hospitable environment for neurofibroma development.

Like Fibroblasts 4, Fibroblasts 5 also expressed genes associated with the peripheral and central nervous systems. These markers include CRLF1, SERPINE2, NEGR1, GRIA2, NRXN1, SOX9, and LGALS3 ( Figure 7E , Supplementary Table 3 ) (100, 108–113). While both these clusters express neural markers, the processes they influence are distinct. Specifically, Fibroblasts 4 expresses markers consistent with promoting cell-cell interactions and processes necessary for synapse formation, but Fibroblasts 5 is enriched for several processes of early neuronal development ( Figure 7F , Supplementary Table 5 ). Perhaps the most notable of these processes is the enrichment of “Tube Development.” While many of the included genes in this set pertain to processes associated with angiogenesis and uretic bud formation, it should be noted that many of these markers have dual roles reliant on context. Specifically, the expression of mRNAs such as SFRP2 and GRIA2 in this cluster may favor smaller sub-processes related to the neural tube or neural development processes, such as “Gliogenesis” (114–116).

Finally, Fibroblasts 6, the smallest group among the six clusters and similar to Fibroblasts 3, has few notable features or upregulated processes. The top processes enriched were primarily driven by overexpression of mitochondrial-derived genes such as MT-ND5, despite that we excluded cells with over 10% representation of mitochondrial genes ( Figure 7E ). Mitochondrially-derived upregulated factors in Fibroblasts 6 favored the enrichment of several redundant processes related to cellular respiration and metabolism ( Figure 7F ). Excluding mitochondrial mRNAs, this cluster was also significantly enriched for processes affecting RNA splicing and mRNA processing ( Figure 7F ). This cluster is notably post-translationally and metabolically active, features common to cells undergoing differentiation (117, 118). Taken together, we provide evidence that each group of fibroblasts is specialized for a unique role in the tumor microenvironment. Despite high expression of neuronal-associated markers in some fibroblast clusters, all were negative for neural crest lineage marker SOX10 ( Figure 4C ) (119–121). The absence of this lineage marker suggests these classes did not arise from neural crest stem cells but may be resident or recruited fibroblasts reprogrammed for the nerve regeneration microenvironment by tumorigenic Schwann cells. Further, this supports the possibility that cell signaling pathways associated with nerve regeneration are drivers of neurofibroma formation (122).

We characterized Schwann cells using the same analysis methods described for the aforementioned cell types ( Figure 8A ). Because of the distinct lack of Schwann cells in the adjacent tissue, comparisons between tumorigenic and non-tumorigenic Schwann cells were not possible from this dataset. Instead, we compared Schwann cells universally and used DEGs to identify enriched pathways. The top three processes enriched in Schwann cells were all associated with processes of development and repair: “Cell Morphogenesis Involved in Differentiation,” “Gliogenesis,” and “Axon Development.” Beyond these top three pathways, the top ten enriched pathways were all related to cell morphogenesis or neural development processes ( Figure 8I ). Within these top ten pathways, there were 59 unique genes upregulated, which accounted for 20.4% of all significantly enriched Schwann cell genes. To confirm that identified pathways were not biased by comparing Schwann cells to non-Schwann cells, which would inherently highlight neuronal pathways, we queried genes representing enriched processes in the Sciatic Nerve Atlas (123). Within the top three enriched gene sets, we found that 52% of the genes displayed maximum expression in immature Schwann cells, which we defined as peak expression occurring in Schwann cells harvested no later than E13.5 or E17.5. Enriched genes representative of immature Schwann cells includes APOA1, CHL1, NRXN1, NGFR, ERBB3, and PTPRZ1 ( Figures 8C–H ). Interestingly, each one of these factors is also associated with various stages of Schwann cell activity during nerve injury (69, 111, 124–129). Other notable enriched mRNAs indicative of a nerve regeneration/repair phenotype include SPP1 and CLU, two secreted factors that promote motor and sensory neuron regeneration, respectively (130).

We also found that roughly 60% of Schwann cells in the microenvironment had significantly higher expression of IL-34, recently discovered to be a second ligand for the macrophage receptor, CSF1R; IL-34 binding to CSF1R potently promotes M2 macrophage polarization ( Figure 8B ) (131, 132). Given that we observed an overall enrichment of M2 macrophages in the tumor microenvironment, this data suggests that neurofibroma-associated Schwann cells may directly promote this polarization and, therefore, directly contribute to reprogramming the microenvironment to be immunosuppressive. Schwann cells are unique in the ability to de-differentiate from a fully differentiated mature state to an immature state following nerve injury (133–135). During this process, an inflammatory response is involved; Schwann cells proliferate, immune cells and other fibroblasts are recruited, and the ECM composition is modified (136). Once this process is complete, Schwann cells stop proliferating and become mature myelinating cells again (133). Considering processes occurring in other cell types, including ECM modifying fibroblasts and immunosuppressive functions of myeloid cells, we conclude that this neurofibroma is under a state of dysregulated nerve repair. The presence of immature Schwann cells, immunosuppressive myeloid cells, multiple types of fibroblasts, and modified ECM deposition supports this notion.