Somatic Landscape of Intracranial Epidermoid Cysts Reveals Immune-Associated Mechanistic and Potential Actionable Targets

1. Introduction

Intracranial Epidermoid Cysts (IECs) are congenital inclusions arising from the neuroectodermal epithelial rests that result from the defective closure of the neural tube between the third and fifth week of embryonic life[1], [2]. IECs are extremely rare, representing only 0.3-1.8% of all intracranial tumors[2], [3].They typically occur in both genders, between the ages of 20–60 years, with a peak incidence in the fourth decade[4]. IECs are often located in the cerebellopontine angle (37.3%) and parasellar region (30%), spreading in the subarachnoid space of the basal cisterns[5]. They may also be found in the middle cranial fossa (18%), diploe (16%), spinal canal (5%)[5], [6], and rarely in the brain stem or medulla[7], [8], [9]. Magnetic resonance imaging (MRI) is the diagnostic approach to identify IECs. Diffusion-weighted imaging allows greater accuracy in the preoperative differential diagnosis[10].

Despite their benign intercourse, tumor growth can be driven by the division of the stratified squamous epithelium lining its cavity. Cyst contents are largely composed of acellular keratin debris and cholesterol inclusions. Liquefaction of the cyst contents is associated with infection or loss of vascularity[11].

IECs have a strong propensity to adhere to critical neurovascular structures, resulting in significant morbidity and neurologic impairment[12], [13]. In addition, recurrent malignant transformation of IECs has been reported [14], [15], [16], including the development of secondary lesions such as malignant melanoma[17] and Squamous Cell Carcinoma (SCC)[18]

The primary treatment for IECs is surgery. Gross total resection is critical to minimize progression[19]; however, cyst adherence results in difficult removal. A subtotal resection yields recurrence rates of up to 93%[20]. Radiation does not seem to have a role in the management of IECs but it has been used as adjuvant therapy for cysts that exhibit malignant transformation[21]. Due to the rarity of IECs, limited research has been conducted focusing on cyst prognosis, diagnosis, and treatment. There is a lack of timely access to molecular testing to determine eligibility for treatment with targeted therapies, difficulty enrolling sufficient numbers of patients in clinical trials, and limited incentives for drug development. Due to the insufficient understanding of the pathophysiology and the molecular characteristics of IECs, no effective systemic or targeted therapy for these patients is available.

Genomic instability is a characteristic of human cancers[22]. With the emergence of exome sequencing, the genetic landscape of many diseases has been unveiled, enabling the identification of actionable targets for new drug development[23] and for which FDA-approved options already exist, highlighting the potential repositioning of these drugs to a new indication[24], and reshaping cancer treatment.

The genetics of IECs remain underexplored. There have been suggestions of a familial inheritance pattern observed for epidermoid cysts of the spleen[25]. Mice lacking the IL-1 receptor (IL-1R) (IL1r^–/–) or deficient in IL1-β developed immunosuppression and multiple epidermal cysts after chronic UVB[26], suggesting that induced somatic events and an altered innate immune response may be involved in the initiation of epidermoid cysts. However, the role of somatic genetic variants in driving trapped fibroblasts to form skull base epidermoid cysts is completely unknown.

In this study, we applied Whole Exome Sequencing (WES) to establish the somatic signature of IECs. A particular focus of our efforts was to gain insights into the biology of IECs. Improved knowledge about the mechanisms of cyst development and progression could enable the identification of potentially actionable variants of clinical impact, and ultimately, bring treatment advances to these patients who still depend on a very limited standard of care.

2. Materials and Methods

Specimens, Patients, and Clinical Data

Studies were conducted following the U.S. Common Rule ethical guidelines. Tumor tissue was resected from participants with the diagnosis of IEC who underwent surgery between 1995 and 2021 at the University of Washington hospitals (Seattle, WA, USA). The respective clinical data was extracted from the University of Washington, School of Medicine clinical database. Data and specimen collection were reviewed and approved by the University of Washington Institutional Review (STUDY R1813). Written informed consent was obtained from all subjects. Samples were collected and stored in tumor biobank (STUDY 00002162) until further processing. Specimens were reviewed by three neuropathologists and neurosurgeons. Clinical data was gathered regarding date, history, demographics, imaging, neuropathology reports, operative information, pre-and post-operative symptoms, adjuvant treatment, and outcomes.

Whole-Exome Sequencing

Genomic DNA (gDNA) was isolated from six fresh-frozen tissue and matching blood (when available) using the QIAmp DNA Mini Kit (Qiagen, Hilden, Germany) and QIAmp DNA Blood Mini Kit (Qiagen, Hilden, Germany), respectively, following the manufacturer’s recommended protocol. DNA was quantified using the Qubit fluorometer (Invitrogen, Carlsbad, CA). Sequencing of gDNA was carried out at the Northwest Genomics Center, Seattle, WA. All sequencing and library preparation were performed at the University of Washington Northwest Genomics Center. Automated library construction and exome capture were carried out in 96-well plate format, using Perkin-Elmer Janus II (PerkinElmer, Waltham, MA) equipment. Five hundred nanograms of genomic DNA were subjected to a series of shotgun library construction steps, including fragmentation through acoustic sonication (Covaris, Woburn, MA), end-polishing, and A-tailing ligation of sequencing adaptors, followed by PCR amplification with dual 10bp barcodes for multiplexing. Libraries underwent exome capture using the Twist+RefSeq exome (36.5MB target) (Twist Biosciences, San Francisco, CA). Briefly, 187.5ng of shotgun library was hybridized to biotinylated capture probes for 16-18 hours. Enriched fragments were recovered via streptavidin beads and PCR amplified. Since each library was uniquely barcoded, samples were captured in multiplex. Prior to sequencing, the pooled library concentration was determined by fluorometric assay, and molecular weight distributions were verified on the Agilent Bioanalyzer (Agilent, Santa Clara, CA) (consistently 180 ± 15bp). Parallel sequencing-by-synthesis with fluorescently labeled, reversibly terminating nucleotides was carried out on the NovaSeq6000 instrument (RTA 3.1.5) (Illumina, San Diego, CA).

Read Processing, Quality Control, and Somatic Variant Calling

Base calls generated in real-time on the NovaSeq6000 were demultiplexed. Raw reads were assessed for Phred score quality using the FastQC tool kit v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Trimmomatic[27] was used to detect and remove sequencing adapters, primers, and low-quality nucleotides. Trimmed reads were aligned against the human reference genome GRCh38 using the Burrows-Wheeler Aligner (BWA) (v0.7.15)[28]. All aligned read data were sorted and subject to duplicate removal using Picard (v2.6.0) (https://broadinstitute.github.io/picard/). Base qualities were recalibrated with GATK BaseRecalibrator (v4.2.6.1) (https://gatk.broadinstitute.org/). Somatic variant calling in unmatched samples was performed against a panel of normal (PoN), using three variant callers: Mutect2 (v4.2.6.1) (GATK), Varscan (v2.4.2)[29], and Vardict[30]. Tumor/normal variant calling was conducted with Mutect2(v4.2.6.1) (GATK), Varscan (v2.4.2)[29], MuSe (v1.0) (MD Anderson Cancer Center), and Strelka[31]. To identify high-confidence mutations, a joint analysis was applied by combining at least three softwares (three-caller pipeline approach)[32]. Mutation calls were selected through a stringent filtering process and functionally annotated with ANNOVAR[33].

Somatic Copy-Number Inference

Somatic Copy Number Variations (CNVs) were identified using CNVkit (v0.9.9)[34]. Whole-exome sequencing read alignments in BAM format, the capture bait locations, and a pre-built human reference were the inputs to the program. All additional data files used in the workflow, such as GC content and the location of sequence repeats, were extracted from the genome sequence in FASTA format using scripts included with the CNVkit distribution. Log₂ copy ratios across the genome for each sample were calculated based on on-target reads and the nonspecifically captured off-target reads. The baseline of normalized sequencing depth on targeted regions was first constructed based on normal samples. For each tumor sample, log2 copy number ratio of normalized sequencing depth on each targeted region between the tumor sample and the baseline was then calculated. A built-in segmentation algorithm was applied to the log₂ ratio values to infer discrete copy number segments[34]. The log₂ ratios and segments were used as inputs for the genomic Identification of Significant Targets in Cancer (GISTIC) v2.0.23[35], allowing the identification of copy number regions that were significantly amplified or deleted across the set of samples.

Data Analysis and Visualization

Search Tool for the Retrieval of Interacting Genes (STRING, http://www.string-db.org/) was used to construct the Protein-Protein Interaction (PPI) network. Genomic data analysis, visualization, and annotation were conducted in R through the usage of Bioconductor packages. Functional enrichment analysis was performed by enrichR[36] with the inclusion of the General Ontology (GO) Molecular Function, Biological Process, and Molecular Component databases[37]. Pathway enrichment analysis was carried out using the Kyoto Encyclopedia of Genes and Genomes (KEGG) 2021 Human database[38]. We utilized the sigminer[39], clusterProfiler[40] , and Maftools[41] packages to extract, analyze, and visualize copy-number aberrations and mutational signatures. For all analyses, p < 0.05 was considered as statistically significant, unless stated otherwise. Biorender software (Toronto, CA) was used for scientific illustration.

Data will be deposited in a publicly available database prior to publication.

3. Results

3.1. Samples and Clinical Data

We performed WES in six tissue samples resected from patients who underwent surgical resection of an IEC from 1995 to 2021 at the University of Washington hospitals (Table 1). WES was also performed in the gDNA of matching blood samples, when available. The mean age of the patient cohort was 36.8 years, with all but one individual being of female gender. Fifty percent of the patients presented with symptoms of headaches, while 33% experienced changes in vision before surgical intervention. Subsequent to diagnosis, all patients underwent craniotomy for tumor resection. Among them, 50% of the cases underwent Gross Total Resection (GTR), two patients (33%) received Subtotal Resection (STR), and one individual (17%) underwent a near Gross Total Resection (nGTR) (Table 1). In H766, an STR posed recurrence/progression risk, which was diagnosed ~7 months post-STR. For 4316, two resections (2011, 2016) resulted in recurrence/progression. The complete clinical information is presented in the Supplementary Table S1.

3.2. Overview of the Somatic Mutational Signature of IECs

The WES workflow is shown in Figure 1A. We performed the gDNA isolation of six fresh-frozen epidermoid samples resected during surgery and respective matching blood, when available. After library preparation, the WES was carried out. Read alignment, somatic variant calling, and VCF filtration were carried out as described in methods. For translational purposes, only exonic non-synonymous alterations and indels were considered for downstream analysis. A total of 1,221 variants affecting 737 unique genes were identified (Supplementary Table S2). From these, ~ 13.00% (158/1221) have been reported in the Catalogue of Somatic Mutations in Cancer - COSMIC database (https://cancer.sanger.ac.uk/cosmic) (Figure 1A).

The median number of variants per sample was 178. The majority of variants identified were missense mutations (1065/1221, 87.22%). In-frame deletions and insertions corresponded to 72/1221 (6.47%) and 15/1221 (1.23%), respectively. Frame shift mutations constituted 2.7% (33/1221) of all somatic alterations. Nonsense and nonstop mutations corresponded to 22/1221 (1.80%) and 5/1221 (0.41%), respectively. Two (0.16%) translation start site alterations were also identified. According to the variant type, 89.43% (1092/1221) were Single Nucleotide Polymorphisms (SNPs), 2.21% (27/1221) were insertions, and 8.35% (102/1221) were deletions. The most predominant Single Nucleotide Variant (SNV) class alterations were C>T (n=408) and T>C (n=231) (Figure 1B). Hypermutated regions were explored by plotting inter-variant distances on a linear genomic scale. Visualization through a rainfall plot revealed changes in inter-event distances located on chromosomes 7, 12, 16, and 19 (Figure 1C). Potential cancer driver candidates identified by the algorithm oncodriveCLUST[42] included the genes PDPR, DNAH2, UBXN11, KIR2DL1/3, FCGBP, PRSS2, MTCH2, LILRB1, HLA-DQB1, HLA-DRB1, and TPTE (Figure 1D). Pair-wise Fisher’s Exact of the top fifty altered genes of the dataset detected pairs of mutually exclusive or co-occurring sets of mutated genes. Relationships were identified when p < 0.1 (Figure 1E).

Figure 1. Overview of the Somatic Mutational Signature of IECs. (A) WES workflow: whole-exome sequencing was performed in six fresh-frozen epidermoids. A total of 1,221 alterations on 737 unique genes were identified. Nearly 13.0% (n=158/1221) mutations have been reported in COSMIC. (B) Data Summary: The median number of variants per sample was 178. Most variants identified were missense mutations (1065/1221, 87.22%) and SNPs (1,092/1,221, 89.43%). Predominant SNV alterations included C>T (n=408) and T>C (n=231). (C) Rainfall plot of hypermutated regions: each point is color-coded according to the SNV class. Hypermutated genomic segments are highlighted by black arrowheads (chromosomes 7, 12, 16, and 19). (D) Prediction of Disease-Associated Driver Genes: driver candidates were predicted by the algorithm oncodriveCLUST. The number of closely spaced mutational clusters is highlighted within brackets. (E) Mutually Exclusive or Co-occurring genes: the triangular matrix displays relationships among the top mutated genes (p < 0.1). Green indicates a tendency toward co-occurrence, whereas pink indicates a tendency toward exclusiveness.

Figure 1. Overview of the Somatic Mutational Signature of IECs. (A) WES workflow: whole-exome sequencing was performed in six fresh-frozen epidermoids. A total of 1,221 alterations on 737 unique genes were identified. Nearly 13.0% (n=158/1221) mutations have been reported in COSMIC. (B) Data Summary: The median number of variants per sample was 178. Most variants identified were missense mutations (1065/1221, 87.22%) and SNPs (1,092/1,221, 89.43%). Predominant SNV alterations included C>T (n=408) and T>C (n=231). (C) Rainfall plot of hypermutated regions: each point is color-coded according to the SNV class. Hypermutated genomic segments are highlighted by black arrowheads (chromosomes 7, 12, 16, and 19). (D) Prediction of Disease-Associated Driver Genes: driver candidates were predicted by the algorithm oncodriveCLUST. The number of closely spaced mutational clusters is highlighted within brackets. (E) Mutually Exclusive or Co-occurring genes: the triangular matrix displays relationships among the top mutated genes (p < 0.1). Green indicates a tendency toward co-occurrence, whereas pink indicates a tendency toward exclusiveness.

3.3. IECs Are Characterized by an Altered Immune Repertoire

GO annotation and pathway analyses were performed on the 737 altered genes. The top five cellular components involved the MHC protein complex, the luminal side of the endoplasmic reticulum membrane, the collagen-containing extracellular matrix, and the MHC class II protein complex (Figure 2A). The molecular function annotation showed extracellular matrix structural constituent, carbohydrate binding, inhibitory MHC class I receptor activity, peptide binding antigen, and MHC protein complex binding (Figure 2B). The main biology processes were associated with the MHC protein complex assembly, peptide antigen assembly with the MHC protein complex, and MHC class II protein complex assembly (Figure 2C) (Supplementary Table S3). The top five pathways included antigen processing and presentation, graft-versus-host disease, cell adhesion molecules, ECM-receptor interaction, and human papillomavirus infection. The KEGG pathway enrichment and respective associated genes are presented in Figure 2D and Supplementary Table S4. Both GO and pathways enrichments demonstrated that many altered genes in IECs are closely associated with immunity, molecular signaling, and interaction, suggesting the involvement of the immune system and microenvironment as a mechanism of intracranial epidermoid formation.

To better visualize alterations in the immune repertoire, we visually assessed the genetic alterations in immune-related genes through an oncoplot (Figure 2E). Notably, somatic mutations on these genes affected 100% of the cohort. The top altered immune-related genes in IECs were Killer cell immunoglobulin-like receptor 2DL1 (KIR2DL1) and Killer cell immunoglobulin-like receptor 2DL3 (KIR2DL3), which comprised multiple missense and multi-hit mutations, affecting ~70% of IECs patients. We also observed DNA variations in other Killer cell immunoglobulin-like receptors (KIR) (KIR2DS1/2/4, KIR3DL1/2, KIR2DL4), HLA genes (HLA-A/B/C, HLA-DRB1/5, HLA-DQB1, HLA-DPB1, HLA-DQA1/2), and additional immune-associated genes which were affected in lower frequencies within the study cohort (Figure 2E, Supplementary Table S2). Taken together, these data indicate an altered immune repertoire associated with IECs, strongly suggesting tumor mechanistic associated with immune evasion.

Figure 2. GO and pathway enrichment. (A) Cellular Component. (B) Molecular Function. (C) Biological Process. (D) Pathway analysis and respective gene networks. (E) Oncoplot of altered immune-associated genes.

Figure 2. GO and pathway enrichment. (A) Cellular Component. (B) Molecular Function. (C) Biological Process. (D) Pathway analysis and respective gene networks. (E) Oncoplot of altered immune-associated genes.

3.4. Oncogenic Driver Candidates of IECs

To identify highly mutated genes that are likely to be the oncogenic drivers of IECs, we considered genes that were found altered in at least 50% of the cohort (>= 3 patients). A set of thirty-three genes met the criteria (Figure 3A, Supplementary Table S5). Except for the Basic salivary proline-rich protein 1 (PRB1) which, according to the literature, does not present a cancer hallmark, all the other genes have been previously related to oncogenic transformation. In addition, over 50% of them have been previously associated with tumor immune infiltration and microenvironment (Table 2), supporting their role as oncogenic drivers. The malignant transformation of IECs into Squamous Cell Carcinoma (SCC) and melanoma has been previously reported[17], [18]. Notably, ~36% (n=12/33) of the most affected genes were also found altered in SCC and melanoma (Table 2 – genes highlighted in bold letters), strongly suggesting a role in the progression of IECs; therefore, may represent interesting targets to further explore in the context of therapeutics to block and/or delay epidermoid malignant transformation.

Variants identified in the driver candidates were frequently found located in specific DNA sites suggesting disease-related mutation hotspots. The hotspots were often positioned in conserved protein regions, domains, and motifs (Supplementary Figure S1), strengthening a potential involvement in the disease.

The functional protein association network analysis revealed three interaction clusters associated with the MHC protein complex, antigen processing, and presentation, suggesting again an altered mechanism to circumvent immune control (Figure 3B, Supplementary Table S6).

Table 2. Characteristics of potential drivers of IECs.

Hugo Symbol	Status	Cancer Type	Oncogenic Function	References
DSPP	Upregulation	Head Neck SCC	ER stress invasion Metastasis Angiogenesis EMT	[43], [44]
FCGBP	Aberrant regulation	Glioma Colorectal SCC	Immune infiltration	[45], [46], [47]
FMN2	NR	Osteosarcoma	Regulation of tumor suppressor pathway mediated by p14ARF	[48]
FRAS1	Upregulation	SCC	Cell growth Metastasis	[49]
GGT2	Upregulation	Endometrial Gastrointestinal SCC Melanoma	Metabolic signaling Cell growth	[50], [51], [52]
GSTT4	Mutation	Acute leukemia SCC	Cellular detoxification	[53], [54], [55], [56]
HLA-B	Mutation Downregulation	Breast Colon Rectal SCC Melanoma	Immune infiltration	[57], [58], [59], [60]
HLA-DRB1	Upregulation Fucosylation	Melanoma Lung SCC	Immune infiltration	[61], [62], [63], [64], [65]
HLA-DRB5	Upregulation	Melanoma Lung SCC	Immune infiltration	[63], [64], [65]
IGFN1	Aberrant expression Mutation	NSCL Renal	Metastasis	[66], [67]
KCNJ12	Mutation	Glioblastoma Prostate SCC	Ion channel	[68], [69], [70]
KIR2DL1	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR2DL3	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR2DS1	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR2DS2	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR3DL1	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
LILRB1	Upregulation	Melanoma Multiple cancers	Immune infiltration	[73], [74], [75]
MAGEC1	Aberrant expression	Cholangiocarcinoma Hepatocellular	Immune infiltration Apoptosis Cell cycle	[76], [77]
MST1L	Downregulation	Breast	Cell differentiation Adhesion Migration Apoptosis	[78]
NOTCH2	Aberrant expression Mutation	SCC Melanoma Multiple cancers	Immune infiltration Stem-like proliferation	[79], [80], [81], [82], [83]
OGFR	Aberrant expression Mutation	Gynecological Ovarian SCC	Cell proliferation Cell cycle	[84], [85], [86]
OR10G2	Aberrant expression	Multiple cancers	Cell differentiation Invasion Metastasis	[87], [88]
OR9G1	Aberrant expression	Multiple cancers	Cell differentiation Invasion Metastasis	[87], [88]
PABPC3	Aberrant expression	Pancreas Osteosarcoma	Immune infiltration Cell proliferation Metastasis	[89], [90]
PDPR	Downregulation	Thyroid	NR	[91]
PRB1	NR	NR	NR	NR
PRSS2	Upregulation	Gastric	EMT Metastasis Microenvironment	[92], [93]
RBMX	Downregulation	Multiple cancers	Tumor development	[94], [95]
SLC35G4	Upregulation	Glioma Breast	Immune infiltration	[96], [97]
TPTE	Upregulation	Lung Prostate	Immune infiltration	[98], [99]
UBXN11	Upregulation Mutation	Multiple cancers	Tumor progression	[100]
USP8	Upregulation Mutation	SCC Melanoma Multiple cancers	Immune infiltration Tumor progression Drug resistance	[101], [102], [103]
ZNF717	Mutation	Colorectal Hepatocellular	NR	[104], [105]

Gene alterations previously reported for Squamous Cell Carcinoma and or melanoma are highlighted in bold letters. SCC: Squamous Cell Carcinoma; NSCL: Non-small-cell lung cancer; EMT: Epithelial-Mesenchymal Transition; ER: Endoplasmic Reticulum; NR: Not Reported.

Table 2. Characteristics of potential drivers of IECs.

Hugo Symbol	Status	Cancer Type	Oncogenic Function	References
DSPP	Upregulation	Head Neck SCC	ER stress invasion Metastasis Angiogenesis EMT	[43], [44]
FCGBP	Aberrant regulation	Glioma Colorectal SCC	Immune infiltration	[45], [46], [47]
FMN2	NR	Osteosarcoma	Regulation of tumor suppressor pathway mediated by p14ARF	[48]
FRAS1	Upregulation	SCC	Cell growth Metastasis	[49]
GGT2	Upregulation	Endometrial Gastrointestinal SCC Melanoma	Metabolic signaling Cell growth	[50], [51], [52]
GSTT4	Mutation	Acute leukemia SCC	Cellular detoxification	[53], [54], [55], [56]
HLA-B	Mutation Downregulation	Breast Colon Rectal SCC Melanoma	Immune infiltration	[57], [58], [59], [60]
HLA-DRB1	Upregulation Fucosylation	Melanoma Lung SCC	Immune infiltration	[61], [62], [63], [64], [65]
HLA-DRB5	Upregulation	Melanoma Lung SCC	Immune infiltration	[63], [64], [65]
IGFN1	Aberrant expression Mutation	NSCL Renal	Metastasis	[66], [67]
KCNJ12	Mutation	Glioblastoma Prostate SCC	Ion channel	[68], [69], [70]
KIR2DL1	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR2DL3	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR2DS1	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR2DS2	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
KIR3DL1	Mutation	Glioblastoma Melanoma	Immune infiltration	[71], [72]
LILRB1	Upregulation	Melanoma Multiple cancers	Immune infiltration	[73], [74], [75]
MAGEC1	Aberrant expression	Cholangiocarcinoma Hepatocellular	Immune infiltration Apoptosis Cell cycle	[76], [77]
MST1L	Downregulation	Breast	Cell differentiation Adhesion Migration Apoptosis	[78]
NOTCH2	Aberrant expression Mutation	SCC Melanoma Multiple cancers	Immune infiltration Stem-like proliferation	[79], [80], [81], [82], [83]
OGFR	Aberrant expression Mutation	Gynecological Ovarian SCC	Cell proliferation Cell cycle	[84], [85], [86]
OR10G2	Aberrant expression	Multiple cancers	Cell differentiation Invasion Metastasis	[87], [88]
OR9G1	Aberrant expression	Multiple cancers	Cell differentiation Invasion Metastasis	[87], [88]
PABPC3	Aberrant expression	Pancreas Osteosarcoma	Immune infiltration Cell proliferation Metastasis	[89], [90]
PDPR	Downregulation	Thyroid	NR	[91]
PRB1	NR	NR	NR	NR
PRSS2	Upregulation	Gastric	EMT Metastasis Microenvironment	[92], [93]
RBMX	Downregulation	Multiple cancers	Tumor development	[94], [95]
SLC35G4	Upregulation	Glioma Breast	Immune infiltration	[96], [97]
TPTE	Upregulation	Lung Prostate	Immune infiltration	[98], [99]
UBXN11	Upregulation Mutation	Multiple cancers	Tumor progression	[100]
USP8	Upregulation Mutation	SCC Melanoma Multiple cancers	Immune infiltration Tumor progression Drug resistance	[101], [102], [103]
ZNF717	Mutation	Colorectal Hepatocellular	NR	[104], [105]

Gene alterations previously reported for Squamous Cell Carcinoma and or melanoma are highlighted in bold letters. SCC: Squamous Cell Carcinoma; NSCL: Non-small-cell lung cancer; EMT: Epithelial-Mesenchymal Transition; ER: Endoplasmic Reticulum; NR: Not Reported.

3.5. Recurrent Variants in NOTCH2 and USP8

After the investigation of the most altered genes, we conducted the analysis focusing on recurrent mutations observed in IECs. The most frequently observed alterations, affecting at least 50% (n=3/6) of the cohort were identified in GSTT4 (p.S227P), KIR2DL1 (p.K237E), KIR2DL3 (p.E295D), MAGEC1 (p.T221S), OR9G1 (p.T83I), PABPC3 (p.H144R), PDPR (p.I106V), PRSS2 (p.R68C and p.C171Y), SLC35G4 (p.T263N), and USP8 (p.R657W). None of these variants has been previously associated with cancer (Figure 3C).

Next, from all the recurrent variants identified, we investigated alterations that are likely to be deleterious, according to the Combined Annotation-Dependent Depletion (CADD) scores[106]. (Figure 3D, Supplementary Table S5). Recurrent and deleterious alterations included the two missense mutations p.R208W and p.R535W in GGT2, involving the gamma-glutamyl transpeptidase domain, p.S149N and p.V209M in HLA-DRB1, both positioned in the Ig-like C1-type domain and Beta-2 region, the missense p.Q677P in the EGF-like calcium-binding domain of NOTCH2, and p.P381S and p.G384C in UBXN11, both located adjacently to the UBX domain. In PDPR, the variants p.I106V and p.N562K are positioned in the chain and the GcvT domain of the protein, respectively. PDPR p.I106V, PABPC3 p.H144R, and USP8 p.R657W occurred in the highest frequency within the cohort, affecting 50% of patients. PDPR p.I106V is in the chain-pyruvate dehydrogenase phosphatase regulatory subunit, which is involved in the regulation of pyruvate metabolism. PABPC3 p.H144R is positioned in the RNA recognition motif (RRM) domain, which is the most prevalent RNA binding domain in eukaryotes, and it is responsible for the high-affinity binding to homopolymeric adenosines for the regulation of mRNA stability and translation initiation[107]. USP8 p.R657W is adjacent to the peptidase C19 domain of the deubiquitinase, which is important to recognize and hydrolyze the peptide bond at the C-terminal glycine of the ubiquitin (Figure 3E, Supplementary Figure S1).

Among the recurrent and predicted deleterious variants, only the GGT2 p.R535W has been previously reported in cancer (COSM3766685) (Figure 3D). USP8 and NOTCH2 are well-known oncogenes that are frequently overexpressed in human cancers. The high expression of NOTCH2 and USP8 has been correlated to a worse prognosis in many tumor types [102], [103], [108], [109], including SCC[82], [108], [109] and melanoma[83], [103]. Consistently, USP8 and NOTCH2 mutations were characterized as gain-of-function alterations [110], [111]. These observations suggest that the USP8 and NOTCH2 variants identified in IECs may lead to an increased protein function to exercise tumorigenesis; nonetheless, additional experimental investigation is needed to fully characterize their functions in the epidermoid formation and progression.

Figure 3. Potential oncogenic drivers and recurrent somatic variants in IECs. (A) Top altered genes: list of the 33 genes found altered in at least 50% of the cohort (>=3 patients). (B) Interaction network: string functional protein association network among the top altered genes. (C) List of top recurrent variants: most recurrent somatic variants, affecting at least 50% (n=3/6) of the cohort. (D) List of recurrent and predicted deleterious variants: recurrent somatic variants with CADD score > 10. (E) Lollipop plots of recurrent and predicted deleterious variants on protein domains. The circle indicates the most frequently occurring mutations in the protein sequence.

Figure 3. Potential oncogenic drivers and recurrent somatic variants in IECs. (A) Top altered genes: list of the 33 genes found altered in at least 50% of the cohort (>=3 patients). (B) Interaction network: string functional protein association network among the top altered genes. (C) List of top recurrent variants: most recurrent somatic variants, affecting at least 50% (n=3/6) of the cohort. (D) List of recurrent and predicted deleterious variants: recurrent somatic variants with CADD score > 10. (E) Lollipop plots of recurrent and predicted deleterious variants on protein domains. The circle indicates the most frequently occurring mutations in the protein sequence.

3.6. Somatic Copy-Number Inference Strengthens the Involvement of Immune-Associated Genes in the Pathogenesis of IECs

The CNV analysis of each sample was performed by using CNVkit[34] (Supplementary Table S8). CNVs stratified by sample, CNV length, and hit counts by chromosome are shown in Supplementary Figure S2. The chromosome 19 was the most affected (Supplementary Figure S2).

To identify somatic copy number regions that are significantly amplified or deleted across all six samples, the log₂ ratios and segments generated with CNVkit[34] (Supplementary Table S7) were used as inputs for GISTIC[35] (Supplementary Table S8). The cohort analysis of the distributions of amplification and deletion events revealed seven cytobands significantly altered in the cohort, among them, amplifications of 7q35, 4p16.1, 17q12 (AP_5 and AP_6), 8p23.1, and 14q11.2, and the deletion of 2p11.1 (Figure 4A). Of these, the events with the lowest False Discovery Rate qvalues calculated for aberrant regions included the amplification of the cytobands 7q35 (AP_5 and AP_6), 4p16.1, and 17q12 (Figure 4A), affecting 67% (n=4/6), 50% (n=3/6), and 33% (n=2/6) of the cohort, respectively (Figure 4B). The list of genes located in each wide peak region is shown in Figure 4C. Interestingly, the pattern of alterations in deubiquitinases and immune-related genes is also present at copy number levels. Recurrent amplifications were predominantly noted in DUB/ubiquitin-specific protease 17 (USP17) family, TBC1 domain protein family, and beta-defensins genes (DEFB4B, DEFB103A, DEFB104A, DEFB105A, DEFB106A, DEFB107A, DEFB131) within the cohort, which overexpression plays a major role in the immune response to cancer[112], [113], [114], [115], [116], [117].

3.7. Potential Actionable Targets in IECs

To generate hypotheses about how the potential driver candidates of IECs might be targeted therapeutically or prioritized for drug development, we matched the list of the thirty-three gene driver candidates against the Drug-Gene Interaction database (DGIdb), a compendium of drug-gene interactions and potentially druggable genes, to prioritize drug-gene interactions[118], [119]. The categorical analysis revealed NOTCH2, which presented the mutations p.Q677P, p.H1024Q, p.N2265K affecting 50% of the IECs, as a clinically actionable target for IECs (Figure 5A). Drug-gene interaction analysis revealed inhibitors and antibody candidates for potential targeted therapies (Figure 5B). Of these, the NOTCH2 inhibitor Nirogacestat is a Food and Drug Administration (FDA) approved chemotherapeutic under the commercial name of Ogsiveo® (SpringWorks Therapeutics). Nirogacestat (Ogsiveo®, SpringWorks Therapeutics) has been indicated to treat adult patients with progressing desmoid tumors who require systemic treatment. The NOTCH2 inhibitors RO4929097 and MK0752, and the NOTCH2 antibodies Tarextumab, and REGN-421 are currently in trial for solid tumors. Additional FDA-approved options that are being currently used to treat cancer include Pazopanib (Votrient®, Novartis) and Lapatinib (Tykerb®, Novartis) targeting HLA-B and HLA-DRB1, respectively (Figure 5B).

4. Discussion

There is an urgent need to understand the mechanisms as well as identify oncogenic drivers underlying the progression of IECs to improve treatment options. Due to the rarity of these tumors, limited research has been done; consequently, these patients continue to rely on a very limited standard of care, resulting in neurological impairment and poor outcomes.

Next-generation sequencing has dramatically reshaped oncological treatment through the identification of genetic variants that provide prognostic information, and aid in therapeutic selection, resulting in hypothesis-driven clinical trials[120]. Since the genetics of IECs was completely unexplored, we sought to apply WES on resected IECs to improve our knowledge of the mechanisms of IEC's oncogenic transformation and identify possible actionable drivers of potential clinical value.

Tumor cells stimulate significant molecular, cellular, and physical changes in their surrounding tumor microenvironment (TME) to keep it continuously evolving as the tumor develops, acting as an active promoter of cancer progression[121]. The GO analysis of the altered gene set identified in IECs revealed a strong association with the MHC complex. The pathway enrichment indicated a tight involvement in antigen processing and presentation, ECM-receptor, and cellular adhesion. Since the TME is composed of blood vessels, extracellular matrix, stromal and immune cells[121], [122], our results demonstrate that genetic alterations in these immune, ECM-receptor, and adhesion genes may have a role in shaping an epidermoid cyst-permissive TME.

IECs were characterized by an altered immune repertoire. The immune cells of the TME are the ones involved in the acquisition of immune escape mechanisms for cancer progression and development[121], [122]. Tumor cells with stronger immunogenicity can be recognized and eliminated by the immune system, while some others can escape from it by several mechanisms and then develop into cancers[123]. Somatic mutations in immune genes affected 100% of the cohort. Top altered immune-related genes included HLA and KIR genes. HLA molecules are fundamental for triggering anti-tumor immunity[124]. KIR is a superfamily of immunoglobulins located on the surface of natural killer (NK) cells, consisting of Ig-like domains that bind to HLA class-I molecules, which help them distinguish between what is “the self” and what is “the non-self”[125], [126]. The innate immune response is the first mechanism activated for tumor immune response, where NK cells play a significant function. The interaction between HLA class-I molecules on normal tissue cells, and inhibitory KIRs on the surface of NK cells shapes the autoimmune tolerance[127]. Taken together, a mechanism of immune evasion in IECs may be based on the HLA-KIR interaction, where the tumor can insert mutations on its LA molecules to decrease antigen presentation[123], [128]. Interestingly, preliminary data demonstrated that mice lacking the IL-1 receptor (IL-1R) (IL1r^–/–) or deficient in IL1-β developed immunosuppression and multiple epidermal cysts after chronic UVB[26], suggesting that induced somatic events and an altered innate immune response may be involved in the initiation of epidermoid cysts, strongly corroborating our findings. The blockade of, for example, the inhibitory KIRs might be a promising immunotherapy strategy to optimize anti-tumor response.

To identify potential oncogenic drivers of IECs, we considered genes that were found altered in at least 50% of the cohort (>=3 patients), resulting in a set of thirty-three driver candidates. The functional interaction network also strengthened immune evasion as the leading oncogenic mechanism of IECs. Our review of the literature revealed that, except for the PRB1, all the other driver candidates have been previously associated with cancer immune infiltration, TME, and/or endoplasmic reticulum (ER) stress, supporting their role as oncogenic drivers in IECs. Interestingly, many members of the list of IEC driver candidates have been reported in SCC and melanoma, probably due to the epidermoid inclusions consisting of ectodermal epithelial/stratified squamous epithelial tissue applied to connective tissue with accumulation of keratin over the epithelium[12]. The transformation of IECs into malignant SCC[18] and melanoma[17] has been reported. Targeting these overlapping genes may be an interesting approach to potentially block the malignant transformation of IECs.

Recurrent mutations on NOTCH2 and USP8 were identified in IECs, affecting 50% of the cohort. The most frequent alterations included USP8 p.R657W and NOTCH2 p.Q677P, which were observed in 50% and 30% of the cohort, respectively. Both USP8 and NOTCH2 are frequently overexpressed in human cancers[102], [103], [108], [109]. Consistently, mutations on these genes confer upregulated protein functionality [110], [111]. These observations indicate that the USP8 and NOTCH2 variants identified in IECs may display also enhanced protein function; therefore, they may be potential candidates for targeted therapy.

The missense p.Q677P in NOTCH2 is positioned in the EGF-like calcium-binding domain of the protein. USP8 p.R657W sits adjacent to the 14-3-3 binding motif (RSYSS) of the protein. Somatic activating mutations in USP8 located between the amino acids 713 and 720 have been identified in (ACTH)-secreting neuroendocrine tumors. In ACTH pituitary adenomas, USP8 mutations are also positioned within or adjacent to the 14-3-3 binding motif (RSYSS). The characterization of USP8 mutations in neuroendocrine tumors revealed that they promote oncogenic transformation via activation of EGF receptor signaling[129]. Together, these observations indicate that USP8 and NOTCH2 alterations may result in downstream dysregulation of EGF receptor signaling in IECs. As in ACTH-pituitary adenomas, USP8 mutations identified in the present work may represent a potential diagnostic marker for IECs. Finally, the inhibition of USP8, NOTCH2, or EGFR should be further explored to tailor treatment strategies for patients with the diagnosis of IECs.

Many challenges are involved in drug development for IECs and other ultra-rare cancers, including the insufficient understanding of the cancer pathophysiology, molecular characteristics, natural history, and lack of timely access to molecular testing to determine eligibility for treatment with targeted therapies, difficulty enrolling sufficient numbers of patients to clinical trials. The use of FDA-approved drugs for new indications may promote the development of safe and effective therapeutic approaches to treat patients with ultra-rare cancers[24]. We performed drug-gene interaction analysis to generate hypotheses about how the potential driver candidates for IECs might be therapeutically targeted. NOTCH2 was highlighted as a clinically actionable target for IECs. The Nirogacestat (Ogsiveo®, SpringWorks Therapeutics) is an FDA-approved NOTCH2 inhibitor indicated for the treatment of adult patients with progressing desmoid tumors. The patient’s treatment with Nirogacestat has provided progression-free survival benefit[130]. It is worth noting that dermoid cysts and epidermoid cysts are similar in structure and origin; therefore, it may represent a very interesting approach for drug repurposing towards IECs. Finally, considering the high genomic instability in deubiquitinases and the indication of a tumor immune evasion mechanism (both also observed at copy-number level), proteasome inhibitors and immunotherapy may also represent treatment options that deserve further investigation in intracranial epidermoid cysts.

5. Conclusions

We acknowledge the limitations of our study. IECs are considered ultra-rare cancers. An experimental design involving samples from a cohort of six patients may be prone to bias due to the low sample sizes. It also makes it difficult to correlate the molecular findings with clinical data, challenging the identification of biomarker candidates. Experimental validation is challenging due to the lack of in vitro or in vivo models. Nonetheless, the molecular information about these rare intracranial cysts is scarce. The relationship between immune system modulation and tumorigenesis in IECs has not been yet characterized. The genetic information demonstrated herein provides a list of driver candidates that can be potentially targeted with existing drugs. It will also contribute to a better understanding of the molecular mechanisms underlying IEC oncogenesis. Additional investigation is needed to better characterize how these mutations may induce changes in surrounding non-tumor immune cells of the TME to confer the tumor immune escape mechanisms. However, the highly reproducible somatic pattern of an altered immune repertoire suggests that the immune escape in IECs is driven by similar underlying molecular mechanisms across most of the patients for tumor progression. These results highlight that the blockade of these targets might be a promising strategy to reverse tumor immune escape and consequently, its progression. Finally, the identification of NOTCH2 and USP8 variations may provide opportunities for novel investigations and the development of tailored approaches to treat IECs.