RNA-sequencing highlights differential regulated pathways involved in cell cycle and inflammation in orbitofacial neurofibromas.

Although most commonly benign, neurofibromas (NFs) can have devastating functional and cosmetic effects in addition to the possibility of malignant transformation. Orbitofacial NFs, in particular, may cause progressive, disfiguring tumors of the lid, brow, temple, face, and orbit, and clinical evidence suggests that they may have increased local aggressiveness compared to NFs developing at other sites. The purpose of this study was to identify biological differences between orbitofacial NFs and those occurring at other anatomic sites. We performed RNA-sequencing in orbitofacial (n = 10) and non-orbitofacial (n = 9) NFs. Differential gene expression analysis demonstrated that a variety of gene sets including genes involved in cell proliferation, interferon, and immune-related pathways were enriched in orbitofacial NF. Comparisons with publicly available databases of various Schwann cell tumors and malignant peripheral nerve sheath tumor (MPNST) revealed a significant overlap of differentially expressed genes between orbitofacial versus non-orbitofacial NF and plexiform NF versus MPNST. In summary, we identified gene expression differences between orbitofacial NF and NFs occurring at other locations. Further investigation may be warranted, given that orbitofacial NF are notoriously difficult to treat and associated with disproportionate morbidity.

INTRODUCTION

Neurofibromatosis type 1 (NF1) is a neurocutaneous disorder and tumor predisposition syndrome with heterogeneous clinicopathologic manifestations affecting approximately 1 in 3000 people worldwide (1). Cutaneous neurofibroma (NF) is the most prevalent tumor in these patients. Plexiform NF, defined as a NF involving multiple peripheral nerve fascicles, is less frequent but more restricted to NF1 patients. Malignant transformation develops in 5–10% of NF1 patients, typically in a plexiform NF subtype (2).

Orbitofacial neurofibromatosis type 1 has been historically studied as a clinical variant of NF1 (3). Affecting between 1% and 22% of patients, orbitofacial NFs may cause progressive, disfiguring tumors of orbital and periorbital structures (4, 5, 6). In several clinical studies, large orbitofacial NFs are felt to be comparatively more aggressive and infiltrative compared to NFs arising in other anatomical body sites (7), although the precise reason for this remains unclear given difficulties to adjust for when comparing anatomic sites. Recurrence after excision is relatively high, but clinical aggressiveness may improve as the individual ages (7, 8, 9). In the study by Chai et al., the 1‐year recurrence rate after surgical reconstruction was 83% with a median time to recurrence of 7 months (10). An orbitofacial NF designation is typically limited to patients with infiltrative NFs involving the orbit and ocular adnexa. In our published diagnostic pathology experience, NF is the predominant nerve sheath tumor type involving the orbit and ocular adnexa. In our study, NFs (n = 63) were predominantly plexiform (38%), diffuse (25%), or with mixed patterns (5%), and the majority developed in patients with known NF1 (89%) (6). Even though orbitofacial NFs may be locally aggressive, malignant transformation appears to be paradoxically rare, with only one case of malignant peripheral nerve sheath tumor (MPNST) reviewed among 90 nerve sheath tumors from 67 patients in this anatomic site (6). The patient, a 63‐year‐old woman, developed clinical progression of a childhood orbitofacial NF over several decades, with eventual malignant transformation, and death with craniospinal dissemination (11).

Recent advances in our understanding of cancer facilitated by high throughput platforms have clarified an important role for genetic and epigenetic events. We know that genetic alterations in plexiform NFs are essentially limited to NF1 mutations (12). Chai and colleagues reported six novel NF1 mutations associated with orbitofacial NF in Chinese families and observed more severe ocular phenotypes in successive generations (10). We have also previously studied epigenetic alterations in orbitofacial NFs using global methylation profiling and found promoter hypomethylation and increased expression of several HOX family gene members (13). We hypothesized that global transcriptomic analysis would be helpful in further clarifying the biological differences between orbitofacial NFs and those developing in other body sites. We have excluded localized cutaneous NFs from the study group and controls given that they are small and have a different biology than orbitofacial NF as we understand it.

METHODS

Tumor samples and controls

Snap frozen tissue from orbitofacial NF (n = 10) obtained from middle eastern patients and non‐orbitofacial NF (n = 4 plexiform, n = 5 localized intraneural) of various ethnicities were studied. All orbitofacial NF were large tumors with plexiform components. No localized cutaneous/dermal NFs were studied. All patients studied satisfied clinical criteria for NF1. Clinicopathologic features of these patients are illustrated in Figure 1 and summarized in Table 1 (orbitofacial NFs) and Table S1 (non‐orbitofacial). Informed consent (or appropriate waiver of consent) was obtained on all patients and the study was performed with local institutional review board approval of the participating institutions. A human Schwann cell line was obtained from ScienCell (ScienCell Research Laboratories Inc., Carlsbad, CA).

FIGURE 1

Clinicopathologic features of orbitofacial neurofibromas. Representative clinical images of orbitofacial neurofibromas including patient 10 (right with associated cutaneous pigmentation) (A), patient 9 (left, with associated skull deformity), and patient 6 (right upper lid neurofibroma)

TABLE 1

Clinicopathologic features of orbitofacial neurofibromas

Case #	Age	Gender	Age at presentation	Family History	Duration of complain	Clinical presentation	Systemic manifestation	Management	Imaging
1	32	Male	24	Negative	5 years	Right upper lid mass and madarosis	No systemic manifestation	Surgical excision	None
2	22	Female	16	Negative	Few years	Right upper lid mass	No systemic manifestation	Surgical debulking	None
3	34	Female	18	Negative	Multiple years	Left upper lid mass	No systemic manifestation	Surgical debulking	None
4	27	Male	20	Negative	6 years	Right eyebrow swelling	No systemic manifestation	Surgical debulking	None
5	44	Male	38	Negative	Multiple years	Right upper lid mass	Pituitary lesion/hamartoma in the brain	Multiple surgical debulking	MRI &CT plexiform neurofibromatosis at the distal zygomatic and temporal branches of the right facial nerve extending into the right extraconal space
6	9	Male	5	Positive	Few years	Left eyelid mechanical ptosis and proptosis/Globe dystopia	Mental retardation	Surgical debulking with ptosis repair	CT scan Left upper lid plexiform neurofibromatosis with sphenoid wing dysplasia
7	36	Female	33	Negative	8 years	Localized right eyelid swelling/facial asymmetry	None	Surgical debulking with ptosis repair and blepharoplasty	MRI + CT scan asymmetrical soft tissue swelling at the right preseptal region in the lower, upper, and cheek
8	36	Female	15	Negative	Since childhood	Right eye blindness/left eye ptosis with lid lesion	Astrocytoma of the right optic nerve	Craniotomy, surgical debulking with multiple ptosis repair	MRI + CT scan craniofacial neurofibromatosis on the left side and right side gliomatosis changes from previous craniotomy
9	21	Male	14	Negative	Since birth	Right eyelid mass	Sickle cell trait and Mitral valve regurge	Craniotomy, Multiple surgical debulking with lateral canthoplasty and ptosis repair	MRI + CT scan craniofacial neurofibromatosis of the right side & orbit extending to the cavernous sinus, sphenoid wing dysplasia
10	10	Female	3	Negative	Since birth	Glaucoma/right facial mass	None	Surgical debulking with multiple ptosis repair	MRI + CT scan craniofacial neurofibromatosis of the right side & orbit extending to the cavernous sinus, sphenoid wing dysplasia

Immunohistochemistry

Immunohistochemistry was performed using a tissue microarray containing orbitofacial and non‐orbitofacial NFs as previously characterized (13). Immunohistochemistry was performed using antibodies directed against WISP2 (Mouse monoclonal, LifeSpan Biosciences Cat# LS‐B6392; 1:250), pERK (Cell Signaling #4370S, 1:250), and MYC (Rabbit monoclonal, Epitomics, catalog #1472‐1, 1:250). A semiquantitative four‐tiered scale was used for interpretation, ranging from negative (0) to strong positivity in a majority of cells (3+) by a neuropathologist.

RNA extraction

Total RNA was extracted using TRIzol Reagent (Life Technologies, cat. #15596026) and RNeasy Mini Kit (Qiagen, cat. #74104) or RNeasy Micro Kit (Qiagen, cat. #74004). Briefly, 10–25 mg frozen tissue was placed into 50‐ml conical tubes and 1 ml TRIzol added. The mixture was placed in a PowerGen 125 homogenizer (Fisher Scientific) for 20–40 s, or a Kontes Pestle Cordless Motor (Fisher Scientific) for smaller samples. Homogenate was left for 5 min at room temperature and then 200 µl chloroform added to the tube and vortexed for 15 s. After Leaving the homogenate for 2 min at room temperature, the sample was transferred into a Phase Lock Gel Heavy tube (5 Prime, yellow‐color tube, spun 20 s prior to use). Samples were centrifuged at 12,000 g for 10 min at 4°C to separate phases, and subsequently, the upper aqueous phase was transferred to a new tube. One volume of 70% ethanol was added and mixed by pipetting. Samples were subsequently pipetted into an RNeasy mini‐column and the RNeasy Mini Kit procedure for RNA isolation was performed per manufacturer recommendations.

Sample preparation and analysis

For the low‐input RNA library preparation workflow, the quality of total RNA was measured by the Agilent Bioanlayzer to determine RNA integrity. Starting material was between 500 pg and 100 ng of total RNA and construction of the whole transcriptome library was prepared as directed in the Nugen Ovation RNA‐Seq System V2 Preparation Guide. Libraries were run on a high sensitivity chip using the Agilent Bioanalyzer to assess size distribution and overall quality of the amplified library. Quantification of the libraries was performed by qPCR or by the Agilent Bioanalyzer and equimolar concentrations of each library were pooled together, clustered and sequencing performed on an Illumina HiseqX for a 151 bp × 151 bp, paired‐end sequencing.

Illumina's CASAVA 1.8.4 was used to convert BCL files to FASTQ files.

Preprocessing and quantification

Quality checks and preprocessing were performed to ensure the quality of the libraries. Gene expression quantification was performed using salmon (v1.3.0) (14) and transcript models from the FANTOM Cage‐Associated‐Transcriptome (FANTOM‐CAT) robust set (15), which we previously used to analyze coding and noncoding gene expression in normal and tumor tissues (16).

Differential gene expression analysis

To detect differentially expressed genes, we normalized the expression values using the trimmed mean of M‐values (TMM) (17) and fitted a generalized linear model (GLM) approach coupled with empirical Bayes moderation of standard errors and voom precision weights (18, 19). The p‐values were adjusted to control for multiple hypothesis testing using the Benjamini‐Hochberg method and genes with false discovery rate (FDR) equal or less than 0.05 were reported (20).

Gene set enrichment analysis

The results from the differential gene expression analysis were ranked by t‐statistics. The ranked lists were tested for gene set enrichment using a Monte Carlo adaptive multilevel splitting approach, implemented in the fgsea (https://doi.org/10.1101/060012) package. A collection of gene sets (Hallmarks, REACTOME, and GO Biological Processes) were obtained from the Broad Institute MSigDB database. Gene sets with less than 15 and more than 1500 genes were removed from the analysis, except for the GO biological processes whose max size was set to 300 in order to avoid overly generic gene sets. The enriched pathways were collapsed to maintain only independent ones using the function collapsePathways from fgsea.

Comparison with publicly available data

To identify consensus of differentially expressed mRNA between orbitofacial NF and benign and malignant Schwann cell neoplasms in other sites, we analyzed differential mRNA expression using publicly available microarray data sets from the NCBI Gene Expression Omnibus database (GSE14038 and GSE32029) (21). mRNA expression profiles that separate different tumor groups were studied in the datasets. The profiles were ranked by t‐statistic and the probability of similarly ranked genes across the contrasts was computed based on a hypergeometric distribution, as previously described (22). Correspondence at the Top (CAT) plots implemented in the R/Bioconductor package matchbox (23, 24) was used to visualize the results.

RESULTS

Orbitofacial and non‐orbitofacial NFs have distinct gene expression profiles

A total of 804 genes (adjusted p‐value < 0.05) were differentially expressed between orbitofacial and non‐orbitofacial NF (333 overexpressed, 471 underexpressed; Table S2); 101 genes were differentially expressed between orbitofacial NF and non‐orbitofacial plexiform NF (23 overexpressed, 78 underexpressed; Table S3); and 573 genes were differentially expressed between orbitofacial NF and non‐orbitofacial localized intraneural NF (187 overexpressed, 386 underexpressed; Table S4). The top genes relatively overexpressed in orbitofacial NF included PRRX1, CXCL14, CCDC80, MMP27, CCN5 (WISP2), EDNRA, XG, MME, FAT4, OSR1, and AR. The top genes relatively underexpressed in orbitofacial NF included GAP43, MEST, SLCO2B1, SDC3, CIITA, SERPINA3, SH3TC1, ITGB8, SLC11A1, HOXB7, and CSF3R. Importantly, the differential expression of these genes remained when comparing to the non‐orbitofacial plexiform and localized intraneural groups individually. Furthermore, the expression of a set of these genes (CXCL14, MMP27, CCN5, XG, OSR1) was particularly high (logFC > 5) in orbitofacial NF compared to a non‐neoplastic human Schwann cell line (Table S5). Hierarchical clustering analysis outlined expression differences between the three different groups (Figure 2).

FIGURE 2

Hierarchical clustering of the top differentially expressed genes between orbitofacial NF and non‐orbitofacial NF types. The color palette and intensity reflect centered and scaled (z‐score) counts‐per‐million

Gene pathways enriched in orbitofacial NFs

The enrichment analysis revealed a variety of gene sets from the REACTOME, HALLMARK, and GO BP collection to be enriched in orbitofacial NF and non‐orbitofacial NF groups. Orbitofacial NF showed pathways involved in cell cycle and proliferation (e.g., MYC targets and G2/M Checkpoint) upregulated in contrast with plexiform and localized intraneural NF. Conversely, immune‐related pathways were found downregulated in orbitofacial NF (e.g., Interferon and cytokine signaling, neutrophil degranulation). Plexiform NF showed an increased expression of genes involved in immune‐related pathways in contrast with localized NF (Figure 3). Next, to further evaluate the significance of these altered genes and pathways, we performed immunohistochemistry, evaluating protein levels (WISP2, MYC) and MAPK pathway activation (pERK) (Figure 4). Of note, a prior study highlighted a role for WISP2 expression and MAPK (pERK) activation in ovarian cancer cells (25). Furthermore, ERK is an important mediator of proliferation resulting from neurofibromin loss (26). We detected high (3+) immunopositivity for pERK in 22/43 (51%) orbital/periocular NFs compared to 7/28 (28%) NFs developing in other anatomic sites, a difference that was statistically significant (p = 0.04, Fisher exact test). Conversely, high expressors of WISP2 were similar between the two groups (20/41, 49% orbital/periocular versus 17/31, 55% nonorbital/periocular NFs, p = 0.64). Results were similar when comparisons included only tumors that were clearly plexiform and/or diffuse (results not shown). MYC immunoreactivity was absent in both groups, suggesting that the MYC target genes were upregulated by mechanisms other than MYC overexpression.

FIGURE 3

Top enriched gene sets enriched across multiple contrasts. Heatmap of mean‐centered log2 signed p‐values (normalized enriched score multiplied by log10 of adjusted p‐value) showing the top 10 enriched gene sets of each collection (ranked by signed p‐value)

FIGURE 4

WISP2 and pERK immunoreactivity in orbitofacial and non‐orbitofacial neurofibromas. Immunoreactivity ranged from 3+ to 1+ based on staining intensity in the majority of the cells

Overlap with differentially expressed genes in the Schwann cell tumors

By comparing publicly available datasets encompassing various Schwann cell primary tumors/cell lines and MPNST, we could observe that the differences between orbitofacial NF and non‐orbitofacial NF were significantly similar to the differences between NF Schwann cells and MPNST cell lines, suggesting that the orbitofacial NF may be more biologically aggressive than the non‐orbitofacial counterparts despite lacking a malignant phenotype (Figure 5).

FIGURE 5

Correspondence‐at‐the‐top (CAT) plot between the differential expression of non‐orbitofacial versus orbitofacial and multiple contrasts. The figure shows the proportion of genes that are similarly differentially expressed in non‐orbitofacial versus orbitofacial in comparison with other contrasts. Lines represent the proportion of overlap with an increasing list size (up to the top 500 differentially expressed genes). Black‐to‐light gray shades represent the decreasing probability of agreeing by chance based on the hypergeometric distribution, with intervals ranging from 0.999999 (light gray) to 0.95 (dark gray). Lines outside this range represent agreement in different cohorts with a higher agreement than expected by chance

DISCUSSION

The wide range of orbitofacial NF clinical manifestations and the difficulties to predict the onset or the severity of new features, consequences, or complications make its management a significant challenge. In addition, orbitofacial NF seems to behave in a more locally aggressive fashion than their counterparts at other body sites, including relatively high recurrence rates. Paradoxically, malignant transformation of these tumors is exquisitely rare, with only one convincing case of transformation to MPNST in our prior studies (11).

The biological basis for this clinical behavior remains understudied, but data are starting to emerge explaining some aspects. Our group recently identified hypomethylation with associated increased expression of individual HOX genes in orbitofacial NF compared with NFs developing at other sites (13). The HOX genes have a better‐understood role in embryogenesis, but there is evidence that they are also involved in stem cell biology (27) and cancer predisposition (28). Another important observation in Chinese families with orbitofacial NF is the phenomenon of anticipation, with more severe phenotypes developing in subsequent generations (10). In the latter study, there was a female predominance (F = 17 and M = 9) and a mean age of 26.9 (range 6–50 years). We observed a similar age in our study set with a mean age of 27 (range 9–44) while there was a 1:1 ratio of females and males.

In the current study, we performed RNA sequencing to uncover further differences between orbitofacial NF and NFs developing at other body sites. We excluded cutaneous/dermal NFs as the pathogenesis of these tumors is distinct and they may have alternative cells of origin (29). In our study, several genes of interest were differentially expressed in orbitofacial NF compared with the other groups. WISP2 (WNT1 inducible signaling pathway protein 2), which is also known as CCN5, was overexpressed in orbitofacial NF compared to the non‐orbitofacial NF (logFC = 6 compared to non‐orbitofacial plexiform NF and logFC = 5 compared to localized intraneural NF) and a non‐neoplastic Schwann cell line (logFC = 8). WISP2 has been found to promote proliferation in ovarian cancer cells probably through ERK signaling (25). At the protein level we detected high WISP2 immunopositivity in approximately half of orbitofacial and non‐orbitofacial NFs, although MAPK signaling (pERK immunoreactivity) was higher in the orbitofacial group. This does not exclude that there are biologically relevant differences of WISP2 at the protein level in orbitofacial NF, given the known limitations of quantifying immunohistochemistry. However, it is intriguing that pERK levels were increased in the orbitofacial group, which is consistent with the gene enrichment analyses documenting upregulation of pathways involved in cell cycle and proliferation, as well as our prior detection of increased ki67 proliferative indices in orbitofacial NF (13).

Orbitofacial NF has also clinically demonstrated increased invasiveness, and histologically contains plexiform and diffuse components. The tumor microenvironment, as in other tumors, is likely to play a role, in the form of complex processes involving matrix metalloproteinases and chemokines. Matrix metalloproteinase 27 (MMP27) and C–X–C motif chemokine ligand 14 (CXCL14) were also overexpressed in the orbitofacial NF group. A variety of metalloproteinases play distinct roles in tumor–microenvironment interactions, including MMP27. CXCL14 is a chemokine that may be secreted by the tumor microenvironment and demonstrates protumor properties in glioblastoma, prostate, and breast carcinoma (30). The androgen receptor (AR) was also relatively overexpressed in the orbitofacial NF compared to the non‐orbitofacial NF groups (but not compared to nonneoplastic Schwann cells). In one study, AR signaling supported tumor growth and angiogenesis in NF (22).

We previously documented increased ki67 proliferative indices in orbitofacial NF compared to NF developing in other sites (13), and there is emerging evidence that inflammation contributes to NF formation (31). Our gene set enrichment analyses also supported these findings. We found that genes involved in cell cycle and proliferation (i.e., MYC targets, G2M checkpoint) pathways were overexpressed in orbitofacial NF. Conversely, genes associated with immune and inflammatory response (i.e., interferon and cytokine signaling and inflammatory response) were found downregulated in orbitofacial NF. We also demonstrated that the biological differences between orbitofacial and non‐orbitofacial NF were significantly similar to the differences observed between Schwann cell primary tumors/cell lines and MPNST by comparing our results with other publicly available datasets, which supports the findings of orbitofacial NF being more aggressive than its non‐orbitofacial counterparts.

There are several caveats in interpreting the findings in our study. First, the orbitofacial NFs were obtained from middle eastern patients while the non‐orbitofacial controls were obtained from various ethnic backgrounds. Second, it is still unclear from the available studies and our data whether the reported local aggressiveness of orbitofacial NF is related to intrinsic biologic properties or to the practical hurdles involved in resecting tumors in a complex anatomical region containing critical structures. Adjusting for these factors will require prospective, multi‐institutional well‐controlled clinical studies in the future. Third, despite our attempt to use various controls to specify the transcriptome of orbitofacial NF, it is at the moment unclear whether the differences in phenotype are related to anatomic site or other intrinsic biologic properties. It would be of great interest in the future to perform functional experiments validating expression data using cell lines and xenograft models. At the moment, models of these disease are lacking, even in well catalogued NF1 repositories (32), further highlighting the difficulties in studying this disease. Despite all these limitations, our data and that reported in the limited clinical literature available suggest that indeed orbitofacial NFs have a distinct biology, and may be associated with increased growth and recurrence rates compared to their non‐orbitofacial counterparts while demonstrating a very low incidence of malignant progression.

In summary, we have identified key biologic differences between orbitofacial NF and NFs developing at other sites, with altered gene levels involved with the tumor microenvironment, proliferation, and inflammatory pathways. Further functional studies and correlation with other genetic factors will be useful in establishing the precise role that each of these genes play in the biology of NF developing in special anatomic sites and hopefully lead to novel, much‐needed therapeutics.

COMPLIANCE WITH ETHICAL STANDARDS/CONFLICT OF INTEREST

The authors declare that they have no conflict of interest. Raw data files will be available in the Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra).

AUTHOR CONTRIBUTIONS

Eddie Luidy Imada generated and analyzed data, conceived the project, wrote portions of the manuscript. Diego Strianese, Deepak P. Edward, Rawan alThaqib, and Hailah Al‐Hussain provided data, reviewed, and approved the manuscript. Antionette Price and Antje Arnold performed experiments, provided data, reviewed, and approved the manuscript. Luigi Marchionni and Fausto J. Rodriguez analyzed data, supervised, conceived the project, and wrote portions of the manuscript.

DISCLOSURE

None.

Table S1 Clinical features of non‐orbitofacial neurofibromas group

Click here for additional data file.

Table S2 Differential gene expression of non‐orbitofacial (plexiform + localized interaneural) versus orbitofacial neurofibromas. Genes with positive logFC are relatively underexpressed in the orbitofacial group

Click here for additional data file.

Table S3 Differential gene expression of non‐orbitofacial plexiform versus orbitofacial neurofibromas. Genes with positive logFC are relatively underexpressed in the orbitofacial group

Click here for additional data file.

Table S4 Differential gene expression of non‐orbitofacial localized intraneural versus orbitofacial neurofibromas. Genes with positive logFC are relatively underexpressed in the orbitofacial group

Click here for additional data file.

Table S5 Differential gene expression of orbitofacial versus non‐neoplastic Schwann cell line. Genes with positive logFC are relatively underexpressed in Schwann cells

Click here for additional data file.