Overlapping variants in the blood, tissues and cell lines for patients with intracranial meningiomas are predominant in stem cell-related genes.

OBJECTIVE: Bulk tissue genomic analysis of meningiomas identified common somatic mutations, however, it often excluded blood-related variants. In contrast, genomic characterisation of primary cell lines that can provide critical information regarding growth and proliferation, have been rare. In our work, we identified the variants that are present in the blood, tissues and corresponding cell lines that are likely to be predictive, tumorigenic and progressive.
METHOD: Whole-exome sequencing was used to identify variants and distinguish related pathways that exist in 42 blood, tissues and corresponding cell lines (BTCs) samples for patients with intracranial meningiomas. Conventional sequencing was used for the confirmation of variants. Integrative analysis of the gene expression for the corresponding samples was utilised for further interpretations.
RESULTS: In total, 926 BTC variants were detected, implicating 845 genes. A pathway analysis of all BTC genes with damaging variants indicated the 'cell morphogenesis involved in differentiation' stem cell-related pathway to be the most frequently affected pathway. Concordantly, five stem cell-related genes, GPRIN2, ALDH3B2, ASPN, THSD7A and SIGLEC6, showed BTC variants in at least five of the patients. Variants that were heterozygous in the blood and homozygous in the tissues or the corresponding cell lines were rare (average: 1.3 ± 0.3%), and included variants in the RUNX2 and CCDC114 genes. An analysis comparing the variants detected only in tumours with aggressive features indicated a total of 240 BTC genes, implicating the 'homophilic cell adhesion via plasma membrane adhesion molecules' pathway, and identifying the stem cell-related transcription coactivator NCOA3/AIB1/SRC3 as the most frequent BTC gene. Further analysis of the possible impact of the poly-Q mutation present in the NCOA3 gene indicated associated deregulation of 15 genes, including the up-regulation of the stem cell related SEMA3D gene and the angiogenesis related VEGFA gene.
CONCLUSION: Stem cell-related pathways and genes showed high prevalence in the BTC variants, and novel variants in stem cell-related genes were identified for meningioma. These variants can potentially be used as predictive, tumorigenic and progressive biomarkers for meningioma.

Introduction

Meningiomas, the most common central nervous system tumours (CNSTs), occur within the arachnoid membrane in multiple extra-axial locations [1, 2, 3]. The World Health Organization (WHO) characterizes meningiomas into 15 histopathological variants, graded I to III. Despite recent developments in targeted therapy, surgery and radiation therapy remain the main methods of treatment for meningiomas, even though both methods, depending on the tumour location, pose post-treatment challenges [4].

Early analysis of meningioma cells identified the tumour-promoting role of the Type 2 neurofibromatosis (NF2)/Merlin gene located on chromosome 22 [5]. Bulk tissue genomic analysis has indicated that approximately 50% of patients who experience biallelic inactivation of NF2 develop meningiomas [6]. However, this association is not necessarily as strong in non-Caucasian populations [7]. Other genetic changes that implicate different genes have also been identified, including phosphoinositide 3-kinase (PI3K), G protein-coupled receptor smoothened (SMO) and, more recently, Forkhead Box M1 (FOXM1) [8, 9, 10, 11].

To date, up to 20% of Grade I tumours reoccur; however, unlike Mib-1, molecular markers approved by the WHO to predict recurrence have not been established [3, 12, 13, 14]. Thus, identifying targetable, predictable and progressive mutations that are highly tumorigenic in meningiomas remains a challenge. While bulk tissue genomic analysis of meningiomas has been useful for the identification of common mutations [6], such analyses have often involved the exclusion of blood-related variants, even though blood is likely to contain cancer derived nucleic acids [15, 16]. Deoxyribonucleic acid (DNA) from commercially available cell lines (HBL-52, Ben-Men-1, IOMM-Lee and CH157-MN) has also been sequenced in an effort to identify mutations that contribute to the progression of meningiomas [17]. However, it is commonly accepted that early passaged primary cell lines derived from patients can represent the “central dogma” of those patients’ tumours more accurately [18, 19]. Thus, genomic characterization of primary cell lines can provide critical information regarding the nature of variants important for growth and proliferation.

We previously published gene expression profiles for meningioma patients’ tissues collected for our cohort, analysed the hetero-dynamic characteristics of cancer stem cells in situ, and characterized the drug-resistant nature of the corresponding cell lines [20, 21]. In this work, we aim to determine rare tumorigenic variants that could potentially contribute to high proliferation and cell culture survival. Using whole-exome sequencing, we identify variants that exist in the blood, tissues, and corresponding primary cell lines (BTCs) for 14 meningioma patients, and we also distinguish their related pathways.

Materials and methods

Patient cohort

Meningioma specimens were collected between February 2013 and December 2015 from 14 patients. The work was approved by the Ethical Board of King Abdulaziz University Hospital (KAUH) (board registration number at the National Committee of Bio. and Med. Ethics is HA-02-J-008) (Project Reference No. 976-12). A signed informed consent form designed according to the Declaration of Helsinki was obtained for each donated tumour sample. Neuropathologists diagnosed the surgical specimens according to WHO classification. The clinical profiles for the patients and their tumours’ histopathological features are shown in Table S1. Haematoxylin and eosin representative sections of histological variants of meningiomas had been previously published [22, 23].

Tumour sampling and cell culture methods

The meningioma specimens were obtained within 30 min after the tumour removal, dissected into three portions, and processed for DNA/RNA extraction [20], tissue freezing at -80 °C [23], and cell culture initiation [22], according to previously published studies. Genomic DNA from the tumour and cell lines was extracted from an average of 30 mg of frozen tissue or 5×106 cells per cell line using All PrepDNA/RNA kits (Qiagen, Hilden, Germany) according to the manufacturer's instructions. DNA from the blood samples was extracted from 3 mL of whole blood using QIAamp DNA Blood Maxi kits (Qiagen), as per the manufacturer's instructions. For the cell lines, cells were collected within an average of 5.1 days (±1.2 days) after the tumour retrieval. All cells were harvested in passage zero before flask expansion and before completing the time taken to double the cells' number.

Exome sequencing

Libraries were generated in the Galician Foundation of Genomic Medicine-SERGAS centre, Spain, using Ion Ampli-Seq™ Exome technology (ThermoFisher Scientific, Massachusetts, USA). Emulsion PCR was performed on diluted and pooled libraries for template preparation using the Ion OneTouch (Life Technologies, California, USA). The templated ion sphere particles were enriched and sequenced using Ion PI chip v2 (Life Technologies) on an ion proton sequencer. Base calling and sequence alignment were performed using Ion Torrent Suite software v.4.2.1. For the text manipulation, concatenate datasets tail-to-head (cat) was used. Sequence mapping was performed with BWA for Ion Torrent using customized build HG19 as the reference genome. BAM file generation and manipulation was performed with Picard application-using tools such as SortSam, Markduplication and BuildBamIndex. HaplotypeCaller was employed in the variant calling GATK. ANNOVAR was utilized for the annotation. The samples had an average total read depth of 106.4 X, and variants were accepted with mutant alleles of at least 20% in the germline. The pipeline used for the VCF analysis is shown in Figure 1. The files were annotated in the BaseSpace Variant Interpreter (accessed on the 28/10/2019). Unfortunately, there is no registered Saudi genome variant source; thus, blood samples from 57 local individuals with no clinical diagnosis of cancer at the time of the sample collection were used to filter out possible ethnicity-related common variants (minor allele frequencies (MEF) > 0.01) (curtsey of Roya Specialised Medical Laboratory at King Abdulaziz University). A preliminary check for COSMIC variants of genes previously known to be associated with meningiomas was carried out (Table S2). PolyPhen-2 Wiki [24], was used to predict damaging consequences for the variants that were reported as unknown by BaseSpace.

Figure 1

The analysis workflow for the processing of files generated for the whole exome sequencing of all samples.

Conventional sequencing

The top seven damaging variants detected in the BTCs for samples with aggressive features were checked using conventional sequencing. Mutational analysis was performed on genomic DNA templates essentially as described in our standard protocol [25]. The primers used for the PCR and sequencing are listed in Table S3. The DNA size of the PCR products was verified through gel electrophoresis and then purified and subjected to cycle sequence reactions using a BigDye Terminator V3.1 Cycle Sequencing kit (Thermo Fisher Scientific Inc., Waltham, MA). The sequencing products were subsequently resolved by capillary electrophoresis on a 3130 Genetic Analyzer (Thermo Fisher Scientific Inc., Massachusetts, USA).

Pathway analysis

Biological functions and implicated pathways were interpreted by employing both Panther [26] and Metascape platforms [27].

Gene expression comparisons for NCOA3 target and partner genes

The PolyPhen damaging rare variant (Gln1274_Gln1276del), within the oncogenic transcription coactivator NCOA3 gene, was the top commonly detected BTCs variant present in only aggressive tumours. Since NCOA3 is a transcription coactivator, a selection of previously published gene expression profiles for the corresponding meningioma samples from GEO submission GSE77259, and three meninges samples from GEO submission GSE100534 were analysed [20, 21] to investigate any possible related effects of this variant. Transcriptome Analysis Console 4.0.2 (Thermo Fisher Scientific, Waltham, MA) was employed to calculate algorithms for average (Log2) expression. To compare means, the samples were grouped into two groups. The first was the NCOA3WT group, which included the low-grade tumours with no NCOA3 variants, the high-grade Jed13_MN tumour with no NCOA3 variants, and three meninges from healthy participants, and the second group was the NCOA3poly−Q group, which included aggressive tumours with the NCOA3 poly-Q variant.

A T-test was applied for 22,216 genes, and 810 genes showed significantly different gene expression when comparing the two groups. The values for these genes were then analysed using the NetworkAnalyst 3.0 platform [28]. In this analysis, files were filtered to remove data that were unlikely to be informative or were simply erroneous. The variance filter was set to a percentile rank of 15, and the relative abundance tab was set to 5. No further normalization was applied since the data had been normalized initially. A Limma statistical analysis was then used for a specific comparison of the differential gene expression values between the two groups.

Results

Basic features of the variants present in BTCs

The basic features of the BTC variants for each meningioma patient in this study are shown in Figure 2. In total, 926 BTC variants were detected that implicated 845 genes, each of which was seen in at least one patient (Table S4). The most variants were detected within chromosome 1 (11.2%), chromosome 19 (8.8%) and chromosome 17 (7.4%). The most common transversions seen in the patients’ exomic regions were C > T and G > A transversions, followed by A > G and T > C transversions. The variants were mainly missense mutations (average: 77.9 ± 12.2%) and inframe deletion (average: 5.4 ± 3.1%). The most common type of zygosity observed when comparing the BTC variants was heterozygous for all (btc) (average: 88.4 ± 15.8%), followed by homozygous for all (BTC) (average: 4.5 ± 3.3%). Overall, variants that were heterozygous in the blood and homozygous in the tissues or the corresponding cell lines were rare (average: 1.3 ± 0.3%). These variants were detected in ficolin-3 (FCN3), transient receptor potential cation channel subfamily m member 8 (TRPM8), ankyrin repeat domain-containing protein 65 (ANKRD65), transmembrane protein 232 (TMEM232), angiomotin (AMOT), cysteine-rich with EGF-like domain protein 2 (CRELD2), LINE-1 type transposase domain-containing protein 1 (L1TD1), apolipoprotein L6 (APOL6), anion exchange protein 3(SLC4A3), cadherin-23 (CDH23), probable tubulin polyglutamylase TTLL2 (TTLL2), apolipoprotein(a) (LPA), E3 ubiquitin-protein ligase ZNRF3 (ZNRF3), MIR205 host gene (MIR205HG), runt-related transcription factor 2 (RUNX2), and coiled-coil domain-containing protein 114 (CCDC114). While the variants detected in RUNX2 and CCDC114 were each observed in two patients, the other variants were each detected in only one patient.

Figure 2

The basic characteristics of the variations present in the BTC for each meningioma patient. A) The number of variants detected in each chromosome; B) the types of detected transversions in the exomic regions; C) the number of all detected variant types based on their consequence; and D) the number of variants for each type of zygosity observed when comparing the variants in the BTCs of each patient. All selected variants had PolyPhen-damaging consequences, according to BaseSpace or PolyPhen-2 Wiki.

Pathways associated with BTC genes that have damaging variants

To determine the possible functional characteristics of BTC genes that have damaging variants, the associated pathways for all 845 genes were identified using the Metascape platform. The top 20 significantly implicated pathways are shown in Figure 3 and Table S5. The top pathway with the highest number of grouped genes (49 genes) was the ‘cell morphogenesis involved in differentiation’ pathway (GO:0000904). For this pathway, on average, 5 variants (±3) per patient were included. The most frequent BTC genes with variants, which occurred in at least three patients, were cholinergic receptor nicotinic alpha 3 (CHRNA3), C-type lectin domain family 1 member B (CLEC1B), EPH receptor A1 (EPHA1), EPH receptor A5 (EPHA5), GLI family zinc finger 3 (GLI3), and protein tyrosine phosphatase receptor type D (PTPRD).

Figure 3

Enriched pathways for all 845 affected genes detected in the BTCs. A) A bar graph for the top 20 significantly enriched pathways based on a combined analysis for all affected genes. Pathways are ranked in a descending order based on the number of genes implicated per pathway. B) Fifteen densely connected network components were calculated using Molecular Complex Detection (MCODE) algorithm9 and were based on three databases: BioGrid6, InWeb_IM7, OmniPath8. The significantly enriched functional descriptions for each complex are as follow: MCODE1 include regulation of transcription by RNA polymerase II (GO:0006357), transcription by RNA polymerase II (GO:0006366), and regulation of transcription, DNA-templated (GO:0006355); MCODE2 include Tight junction (hsa04530), ribonucleoprotein complex biogenesis (GO-0022613), and ribosome biogenesis (GO-0042254); MCODE3 include epoxygenase P450 pathway (GO-0019373), Phase I - Functionalization of compounds (R-HSA-211945), and arachidonic acid metabolic process (GO-0019369); MCODE4 include Fanconi Anemia Pathway (R-HSA-6783310), FA core complex (Fanconi anemia core complex) (CORUM-1625), and FA complex (Fanconi anemia complex) (CORUM-1152); MCODE5 include 55S ribosome mitochondrial (CORUM-320), Mitochondrial translation elongation (R-HSA-5389840), and Mitochondrial translation termination (R-HSA-5419276); MCODE6 include MET activates PTK2 signaling (R-HSA-8874081), MET promotes cell motility (R-HSA-8875878), and Collagen chain trimerization (R-HSA-8948216); MCODE7 include G alpha (q) signalling events R-(HSA-416476), Class A/1 (Rhodopsin-like receptors) (R-HSA-373076), and GPCR ligand binding (R-HSA-500792); MCODE8 include transforming growth factor beta receptor signaling pathway (GO-0007179), metabolic process (GO-0046710 GDP), glomerular visceral epithelial cell development (GO-0072015), and positive regulation of mRNA catabolic process (GO-0061014); MCODE9 include Resolution of D-loop Structures through Synthesis-Dependent Strand Annealing (SDSA) (R-HSA-5693554), Resolution of D-loop Structures through Holliday Junction Intermediates (R-HSA-5693568), and Resolution of D-Loop Structures (R-HSA-5693537); MCODE10 include: rRNA modification in the nucleus and cytosol (R-HSA-6790901), Major pathway of rRNA processing in the nucleolus and cytosol (R-HSA-6791226), and rRNA processing in the nucleus and cytosol (R-HSA-8868773); MCODE11 include Major pathway of rRNA processing in the nucleolus and cytosol (R-HSA-6791226), rRNA processing in the nucleus and cytosol (R-HSA-8868773), and rRNA processing (R-HSA-72312); MCODE12 include Striated Muscle Contraction (R-HSA-390522), muscle filament sliding (GO-0030049), and actin-myosin filament sliding (GO-0033275); MCODE13 include positive regulation of metanephric cap mesenchymal cell proliferation (GO-0090096), cotranslational protein targeting to membrane (GO-0006613), and regulation of post-translational protein modification (GO-1901873); MCODE14 include: positive regulation of protein insertion into mitochondrial membrane involved in apoptotic signaling pathway (GO-1900740), mitotic cytokinetic process (GO-1902410), and negative regulation of angiogenesis (GO-0016525); MCODE15 include G alpha (i) signalling events (R-HSA-418594), and GPCR ligand binding (R-HSA-500792).

The second most common pathway (48 genes) was the microtubule-associated pathway GO:0007017 ‘microtubule-based process’. This pathway included an average of 8 variants (±3) per patient, and the most frequent genes with variants that were detected in at least three patients were Dynein Axonemal Heavy Chain 3 (DNAH3), RPGRIP1 Like (RPGRIP1L), centrosomal protein 350kDa (CEP350), FERM domain containing 7 (FRMD7), coiled-coil domain containing 114 (CCDC114), cilia- and flagella-associated protein 100 (CFAP100), kinesin family member 1A (KIF1A), kinesin family member 1B (KIF1B), pericentrin (PCNT) and tubulin beta 8 class VIII (TUBB8). The third pathway, which included 33 genes, was the ‘NABA CORE MATRISOME’ pathway (Canonical Pathways M5884). On average, four variants (±1) per patient were included.

Genes with BTC variants that occurred in at least four patients

Thirty-one BTC variants were present in at least four patients (4/14 patients, 29%; Figure 4). The most frequently detected BTC variant, which was the inframe insertion c.727_728insAGGTGGGGG,p.(Arg242_Ala243insGluValGly) present in the G protein regulated inducer of neurite outgrowth 2 (GPRIN2) gene, was seen in 10 patients. Three variants were detected in the dynein axonemal heavy chain 3 (DNAH3) gene. The variants c.10933C > T,p.(Arg3645Cys) and c.8846A > C,p.(Lys2949Thr) were detected in six patients, and the variant c.11230C > T,p.(Arg3744Trp) was seen in five patients. The variant c.2231G > A,p.(Arg744Gln) in RPGRIP1 like (RPGRIP1L) was detected in samples from six patients.

Figure 4

Critical variants detected for all meningioma patients' BTCs and frequent in at least 4 patients. Only variants with a possible/probable damaging PolyPhen effect were included, as per data annotated by BaseSpace or as detected manually using PolyPhen-2 Wiki. Group 1 tumours are all WHO classified as Grade I.

Other variants that were observed in five patients were detected in genes for heparan-sulfate 6-O-sulfotransferase 3 (HS6ST3), aldehyde dehydrogenase family 3 member B2 (ALDH3B2), Asporin (ASPN), LIM domain only protein 7 (LMO7), uromodulin-like 1 (UMODL1), thrombospondin type-1 domain-containing protein 7A (THSD7A), leucine-rich repeat, immunoglobulin-like domain and transmembrane domain-containing protein 1 (LRIT1) and sialic acid-binding Ig-like lectin 6 (SIGLEC6). Two variants were detected in the ATP-binding cassette sub-family A member 10 (ABCA10) gene; c.3964C > T,p.(Arg1322Ter) was detected in five patients, and c.1331_1334delCTGT,p.(Ser444PhefsTer17) was detected in four patients.

Variants detected only in tumours with aggressive features

Five tumours were characterised as having aggressive features. These were the Grade II Jed13_MN and Jed18_MN tumours, the recurrent Jed49_MN tumour, the metastatic Jed45_MN tumour and the drug-resistant Jed38_MN tumour; the latter was previously shown to be highly proliferative in culture [22]. Damaging variants present in the BTCs in only this group of samples, and never in the Grade I tumours, were selected in order to identify variants that are likely associated with aggressive meningioma phenotypes. In total, 240 BTC genes with damaging variants were detected (Table S6). The most common variants present in at least two patients and that were conventionally sequenced to confirm their presence are shown in Table 1 and Figure S1. These included variants in the genes for nuclear receptor coactivator 3 (NCOA3), collagen alpha-1 (XXVII) chain (COL27A1), collagen alpha-1(V) chain (COL5A1), kinesin-like protein KIF13B (KIF13B), kinetochore-associated protein 1 (KNTC1), protein spinster homolog 1 (SPNS1), and somatostatin receptor type 5 (SSTR5).

Table 1

Rare damaging variants detected in the BTCs only in patients with aggressive tumours, and in at least two patients within that group.

Gene Symbol	Variant Details	Consequence	Jed13 _MN∗	Jed18_MN∗	Jed38_MNˆ	Jed45_MN∗∗∗	Jed49_MN∗∗
NCOA3	NM_181659.2, c.3753_3761delGCAGCAGCA, p.(Gln1274_Gln1276del)	ID		btc√	btc		Btc
COL27A1	NM_032888.2, c.3878C > T, p.(Thr1293Met)	MV	btc	btc√			Btc
COL5A1	NM_000093.4, c.1588G > A, p.(Gly530Ser)	MV	btc	btc√
KIF13B	NM_015254.3, c.1105C > T, p.(Arg369Trp)	MV	btc				btc√
KNTC1	NM_014708.4, c.6062T > G, p.(Val2021Gly)	MV		btc√		btc
SPNS1	NM_032038.2, c.1550G > A, p.(Arg517His)	MV	btc	btc√
SSTR5	NM_001172560.1, c.142C > A, p.(Leu48Met)	MV				btc	btc√

B:Blood, T:Tissue, C:Cell line. Small cap: Heterozygous, Large cap: Homozygous. MV: missense variant, ID: inframe deletion.√: indicates the variant was checked using conventional sequencing.

Grade II.

Recurrent.

Metastatic.

High mitotic index in culture.

Differential pathways in grade I tumours and tumours with aggressive features

To identify the pathways that differentiate between low-grade tumours and tumours with aggressive features, all the genes in the BTC variants for each group were assembled into a list, and the two corresponding lists were compared using a Metascape multiple list analysis (Figure 5). The top 20 clusters, with their representative enriched terms for low-grade tumours and tumours with aggressive features, are shown in Table S7. Common pathways included pathways for ‘microtubule-based process’ (GO:0007017), ‘cell morphogenesis involved in differentiation’ (GO:0000904), ‘NABA CORE MATRISOME’ (M5884), and ‘diseases associated with O-glycosylation of proteins’ (R-HSA-3906995). Those that were particularly associated with aggressive tumours included pathways for ‘homophilic cell adhesion via plasma membrane adhesion molecules’ (GO:0007156) and ‘rRNA catabolic process’ (GO:0016075).

Figure 5

Differential pathways between low-grade and aggressive tumours. A) A Metascape heatmap of enriched terms across lists for the total genes with BTC variants for low grade and aggressive tumours. Pathways are coloured by p-values. B) MCODE components identified from combined networks. Each network node is displayed as a pie and colour coded red for genes identified in aggressive tumours and blue for genes identified in low grade tumours. The pathways for each component are as follows: MCODE1 include the PID FANCONI PATHWAY (M1), the epoxygenase P450 pathway (GO:0019373) and the DNA Repair pathways (R-HSA-73894); MCODE2 include rRNA processing (GO:0006364), rRNA metabolic process (GO:0016072), and ribosome biogenesis (GO:0042254): MCODE3 include Collagen biosynthesis and modifying enzymes (R-HSA-1650814), Collagen formation (R-HSA-1474290) and MET activates PTK2 signaling (R-HSA-8874081); MCODE4 include axonogenesis (GO:0007409), axon development (GO:0061564) and cell morphogenesis involved in neuron differentiation (GO:0048667). C) MCODE components identified from the top enriched pathways for aggressive tumours. The pathways for each component are as follows: MCODE1 include Major pathway of rRNA processing in the nucleolus and cytosol (R-HSA-6791226), rRNA processing in the nucleus and cytosol (R-HSA-8868773), and rRNA processing (R-HSA-72312); MCODE2 include regulation of transcription by RNA polymerase II (GO:0006357), transcription by RNA polymerase II (GO:0006366), and regulation of transcription, DNA-templated (GO:0006355); MCODE3 include 55S ribosome, mitochondrial (CORUM:320), Mitochondrial translation initiation (R-HSA-5368286), and Mitochondrial translation elongation (R-HSA-5389840); MCODE4 include axonogenesis (GO:0007409), axon development (GO:0061564), and cell morphogenesis involved in neuron differentiation (GO:0048667); MCODE5 include G alpha (q) signalling events (R-HSA-416476), Class A/1 (Rhodopsin-like receptors) (R-HSA-373076), and GPCR ligand binding (R-HSA-500792); MCODE6 include Antigen processing: Ubiquitination & Proteasome degradation (R-HSA-983168), protein polyubiquitination (GO:0000209), and Class I MHC mediated antigen processing & presentation (R-HSA-983169).

Pathways that were particularly common in the low-grade tumours included ‘positive regulation of cellular component biogenesis’ (GO:0044089), and ‘cell-substrate adhesion’ (GO:0031589). The protein–protein interaction enrichment analysis generated a Molecular Complex Detection (MCODE) algorithm for four nodes of proteins implicated in pathways that overlapped between the two groups (Figure 5B). In the aggressive tumours, six nodes of MCODE were identified (Figure 5C).

Expressed genes potentially affected by the poly-Q variant detected in the NCOA3 gene

The PolyPhen damaging rare variant (Gln1274_Gln1276del) within the oncogenic transcription coactivator NCOA3 gene was detected in the BTCs of three patients with aggressive tumours and zero patients with low-grade tumours. This variant shortens the poly-Q region of the gene, which is important for the acetylation activity of proteins and protein–protein interactions [29]. To investigate whether this variant might have a functional effect in meningiomas, we analysed the gene expression profiles for NCOA3, its target genes, and its binding partners that were previously published for the corresponding meningioma samples [20]. The expression levels in the low-grade tumours (Jed04_MN, Jed36_MN, Jed40_MN, Jed43_MN and Jed64_MN), Jed13_MN and three meninges samples referred to as NCOA3wt, which did not have any variants present in the NCOA3 gene, were compared to the expression levels in the aggressive tumours (Jed18_MN, Jed38_MN, Jed49_MN), which had the variant referred to as NCOA3poly−Q. No expression data for Jed45_MN was available. The mean level of expression of the NCOA3 gene was significantly higher in the low-grade and aggressive meningioma samples than in the meninges (meninges: 8.68 ± 0.40; low-grade tumours: 10.43 ± 0.18; all aggressive tumours: 9.96 ± 0.20; significance: P = 0.01 and P = 0.02, respectively). However, this effect was not observed when comparing the NCOA3 gene expression levels for the NCOA3wt and NCOA3poly−Q groups (averages: 10.09 ± 0.77 and 10.02 ± 0.31, respectively; P = 0.841). To identify potential target genes that significantly differ in their gene expression between NCOA3wt and NCOA3poly−Q, a T-test was applied to 22,216 genes. A total of 815 significantly differential genes that showed significance below 0.05 (Table S8) were selected to be analysed further on the NetworkAnalyst platform using the Limma statistical method in order to identify the top significantly differential genes. Fifteen genes were calculated to show significant differential expression (Table 2). The genes that were upregulated by at least one fold were eyes shut homolog (drosophila) (EYS); sema domain; immunoglobulin domain (Ig); short basic domain; secreted (semaphorin) 3D (SEMA3D); family with sequence similarity 20, member A (FAM20A); gamma-aminobutyric acid (GABA) A receptor, beta 2 (GABRB2); C3 and PZP-like, alpha-2-macroglobulin domain containing 8 (CPAMD8), stimulated by retinoic acid 6 (STRA6); vascular endothelial growth factor A (VEGFA); ATP binding cassette subfamily A member 6 (ABCA6); guanine nucleotide binding protein (G protein), gamma 11 (GNG11); leucine rich repeat containing 31 (LRRC31); and early B-cell factor 4 (EBF4). The genes that were downregulated by at least one fold were NCK associated protein 5 (NCKAP5), ornithine decarboxylase 1 (ODC1), BCL2 apoptosis regulator (BCL2), and sortilin related receptor 1 (SORL1).

Table 2

NCOA3wt and NCOA3polyQ significantly differential genes identified by the NetworkAnalyst platform.

Symbol	Average meninges	Average NCOA3wt	Average NCOA3polyQ	logFC	adj.P.Val	Main functions	Possible link with NCOA3#
EYS	6.37	6.44 ± 1.82	10.71 ± 1.09	3.446	0.033	Calcium ion binding, skeletal muscle tissue regeneration	Regulated by NCOA1ˆ
SEMA3D	5.8	5.55 ± 1.37	8.68 ± 1.29	3.133	0.047	Regulation of neuron differentiation, cell migration	Repressed by AHRˆ
FAM20A	7.91	7.6 ± 0.81	10.38 ± 0.52	2.779	0.02	Phosphotransferase activity, biomineral tissue development	Repressed by NCOA2ˆ
GABRB2	7.94	7.37 ± 1.17	9.95 ± 0.49	2.577	0.04	Chloride channel activity, nervous system process	Regulated by AHRˆ
CPAMD8	8.14	7.64 ± 0.67	9.59 ± 0.64	1.947	0.033	Serine-type endopeptidase inhibitor activity, eye development	Repressed by Junˆ
STRA6	6.92	6.91 ± 0.73	8.66 ± 0.66	1.744	0.043	Transmembrane transport	Regulated by AHRˆ and Jun
VEGFA	8.72	8.04 ± 0.68	9.71 ± 0.45	1.667	0.037	Growth factor receptor binding, positive regulation of angiogenesis, response to hypoxia	Repressed by FOXA1ˆ
ABCA6	6.03	5.15 ± 0.59	6.74 ± 0.12	1.589	0.033	ATPase-coupled transmembrane transporter activity	Regulated by CCND1ˆ, AHRˆ and FOXA1ˆ
GNG11	8	7.85 ± 0.51	9.24 ± 0.53	1.387	0.037	G protein-coupled receptor signaling pathway for Dopamine	Regulated by CCND1ˆ, NOCA2ˆ and FOXA1
LRRC31	4.28	4.65 ± 0.44	5.79 ± 0.53	1.14	0.046	Influencing telomere length, associated with high-grade glioma risk	Regulated by B-Catenin and FOXA1
EBF4	7.7	7.37 ± 0.18	8.41 ± 0.26	1.046	0.017	Positive regulation of transcription by RNA polymerase II	Repressed by AHRˆ
NCKAP5	6.43	6.85 ± 0.29	5.82 ± 0.23	-1.027	0.033	Microtubule bundle formation, microtubule depolymerization	Regulated by CCND1ˆ and NOCA2ˆ
ODC1	9.91	10.5 ± 0.4	9.34 ± 0.27	-1.165	0.033	Ornithine decarboxylase activity, cellular nitrogen compound biosynthetic process	Regulated by AHRˆ, FOXA1ˆ, NCOA1 and NCOA2
BCL2	7.85	8.5 ± 0.48	7.23 ± 0.27	-1.267	0.036	Regulation of intrinsic apoptotic signaling pathway, negative regulation of apoptotic process	Regulated by AHRˆ, and NCOA2
SORL1	10.05	10.05 ± 0.36	8.39 ± 0.49	-1.666	0.017	Post-Golgi vesicle-mediated transport	Regulated by CEBPAˆ and NCOA2ˆ

The logFC and adj.P.Val are significant at the .05 level.

Possible link with NCOA3 according to [95].

interacts with NCOA3 according to [78].

Discussion

Our understanding of the genomic nature of meningiomas comes from studies that analysed bulk tumour masses collected from Western Europe and American patients and focused on somatic mutations by excluding variants detected in the blood. However, tumour cells and tumorigenic DNA can escape into the blood; thus, DNA collected from blood samples may not necessarily represent a pure, non-tumorigenic profile [15, 16]. In addition, the microenvironment of meningioma tissue is heterogeneous and composed of both tumour and non-tumour cells [30]. Selecting for tumorigenic cells via primary cell lines initiation could be useful, as the conditions used in culture initiation do not permit non-tumorigenic cells to grow. Meningioma cells grown in culture must be able to adapt, attach, and divide; thus, they are often aggressive and progressive [22]. However, since selected cultured cells are genomically unstable, new populations are likely to emerge that harbour new mutations, which might not essentially contribute to the progression of meningiomas. Therefore, in order to overcome the limitations associated with each sample type, whole exome sequencing was used to identify variants that are present in the BTC for each of the 14 patients. Importantly, using this approach enabled the identification of variants in stem cell related genes that have not been previously particularly associated with meningioma. In addition, transforming variants seen in patients who had aggressive tumours were selected, and an analysis was carried out to investigate possible effects of the most frequent variant in the stem cell related NCOA3 gene, present in the aggressive meningiomas.

Most variants were detected in chromosomes 1, 17 and 19, which are consistent with regional changes detected previously in meningiomas [11, 31]. Similar to other tumours, the most common variants here were missense mutations, with a C > T transversion in the exomic regions [32, 33, 34]. Tri-homozygosity among variants detected in the BTCs were very rare, and only 16 variants showed a heterozygous pattern in the blood, along with a homozygous change in the tumour or cell line. Of interest were variants detected in the RUNX2 and CCDC114 genes, as they occurred in samples from two of the patients. RUNX2 has been reported in several cancers as a tumour suppressor [35], and the gene has been associated with the modulation of different oncogenic processes and cancer signalling pathways; thus, it might also play a crucial role in meningiomas. Mutations in the CCDC114 gene have been associated with the development of primary ciliary dyskinesia, deafness, and renal disease [36], and a familial variant has been detected in glioblastoma patients and their parents [37].

The top identified pathway associated with the BTC genes that have damaging variants was the stem cell-related “cell morphogenesis involved in differentiation” pathway. This pathway includes genes that are responsible for the form change (in cells’ shape and size) that occurs upon the differentiation of embryonic or regenerative stem cells. This pathway has been associated with several tumour types, including those in prostate and colorectal cancers [38, 39], metastatic gastric cancers [40] and high grade gliomas [41]. The association of this pathway with BTC mutated genes is consistent with the notion that stem cell-related germline mutations are core to tumorigenesis in meningiomas [22, 23, 42, 43]. Of interest in this pathway is the variant (c.4609C > T,p.(Arg1537Cys), which was detected in the GLI3 gene. The protein encoded by this gene is a transcription factor, which is active as part of the Hedgehog (Hh/HH) signalling pathway known to be important for tissue development, and its aberration has been associated with several cancers, including GBM and medulloblastoma [44]. The variant in the Gli3 gene is localized within the transcription domain-1, and thus could interfere with the transcription of its target genes, tipping the balance in favour of tumorigenesis.

Five stem cell-related genes GPRIN2, ALDH3B2, ASPN, THSD7A and SIGLEC6 showed BTC variants in at least four patients. The most frequent variant was detected in the GPRIN2 gene. The mouse homologue to GPRIN2 was found to be enriched in neural growth cone membranes of embryonic mouse brain cells collected at day 17 of embryogenesis, and has shown to be important for the control of neurons growth [45, 46]. G protein regulated inducers of neurite outgrowth proteins function as intermediates for the communication between G protein-coupled receptors and sequential intracellular targets [47]. GPRIN2 has been observed to bind to Gαo-protein-activating Cdc42, resulting in changes to its cellular morphology [48]. Gene overexpression has been detected in endometrial and head and neck cancers [49]. Studies have reported GPRIN2 mutations in different types of cancers, including melanomas [50] and familial oesophageal squamous cell carcinomas [51]; however, the protein mechanism of action in tumorigenesis has not been clarified. The ALDH3B2 gene encodes a member of the aldehyde dehydrogenase family, a group of isozymes that play a role in the detoxification of aldehydes generated by alcohol metabolism and lipid peroxidation. Gene overexpression has been associated with breast cancer stem cells [52]. The ASPN gene encodes the Asporin protein, which is thought to regulate chondrogenesis by inhibiting transforming growth factor-beta 1-induced gene expression. In the peritumoral stroma of pancreatic cancer tissues, Asporin overexpression has been associated with poor clinical outcomes and has been shown to promote mesenchymal transition, invasion and migration of pancreatic cancer cells through the activation of the stem cell-related CD44-AKT/ERK-NF-κB pathway [53]. The THSD7A gene encodes a membrane-associated N-glycoprotein with a soluble form found almost exclusively in placenta and umbilical cord endothelial cells [54]. The gain and loss of THSD7A expression in a large cohort of samples for renal cell carcinomas and colorectal, breast and prostate cancers, compared to the expression status in non-tumour tissues, have been linked with tumorigenicity, and the protein expression levels are clinically relevant in that study [55]. The SIGLEC6 gene encodes a transmembrane receptor that is expressed in the placenta [56]. This receptor has also been shown to have an immunomodulatory role by binding to sialyl-TN glycans and leptin [57]. Protein expression has been detected in mast cells present in colorectal cancer and is thought to act as a functionally inhibitory receptor [58]. The detected BTC variant is located in the region that codes for the Immunoglobulin I-set domain and thus could compromise its function in, perhaps, the mast cells present in meningiomas [59]. Unfortunately, the ways in which the aforementioned variants are functionally connected with meningioma cancer stem cells is not yet clear. Of note, all five stem cell-related genes were not automatically identified by any of the Panther, Metascape or NetworkAnalyst platforms as genes associated with stem cell pathways. This indicates the limitation of automated pathway analysis when new functions are continuously identified for genes, and it also indicates a stronger stem cell-related characteristic of the BTC variants.

The second most common variant, along with two other variants, were identified in the DNAH3 gene. This gene, a member of the dynein family, encodes a large protein that forms part of a microtubule-associated motor protein complex [49]. The c.8846A > C,p.(Lys2949Thr) variant is located in the microtubule-binding stalk of the dynein motor region, which is required for protein interaction with microtubule strands, while both the c.10933C > T,p.(Arg3645Cys) and c.11230C > T,p.(Arg3744Trp) variants lay in the dynein heavy chain AAA lid domain. The second variant in this gene has been previously detected in head and neck squamous cell carcinoma [60]. Recent studies have revealed other mutations in the DNAH3 gene in lung adenocarcinoma and breast cancers [61, 62]. Kinesin- and dynein-related proteins are known to function generally as cargo transporters and microtubule organizers [63]. However, precisely which proteins’ transportation is affected by the dysfunction of DNAH3 is not currently known.

Other genes that had BTC variants detected in several patients included ABCA10, RPGRIP1L, HS6ST3, LMO7, UMODL1 and LRIT1. Two variants in the ABCA10 gene were detected: c.3964C > T,p.(Arg1322Ter), which is located in the second ABC_tran region, and c.1331_1334delCTGT,p.(Ser444PhefsTer17), which is located in the first ABC_tran region. The second variant has been detected in breast, lung, thyroid, haematopoietic and glioma cancer cell lines [64]. ABCA10 belongs to the ATP-binding cassettes transporter family, which is involved in several physiological functions, including the movement of various xenobiotics and drugs across the cell membrane [65]. These transporters show higher expression in chemotherapeutics-resistant cancer cell lines [66, 67, 68, 69]. Variations in the ABCA10 gene have also been detected in blood samples of GBM patients and their families through familial exome sequencing [37]. RPGRIP1L is a body basal tumour suppressor that has been shown to suppress anchorage-independent growth in hepatocellular carcinoma, partly through the mitotic checkpoint protein Mad2 [70]. Aberrations in the HS6ST3 gene have been seen in breast cancers, and silencing this gene significantly changes the expression of IGF1R and XAF1 in breast cancer cells [71]. The LMO7 gene, which encodes a protein that contains a calponin homology (CH) domain, a PDZ domain, and a LIM domain, was recently discovered as a BRAF-linked fusion in papillary thyroid carcinoma [72]. In contrast, although the variant detected in the UMODL1 gene is probably damaging, associated research indicates a mild influence of this gene on cancer [73]. Similarly, although the LRIT1 gene is expressed in neuronal cells, and gene knockout in mice impairs adaption to background light, the gene does not particularly influence the morphology or molecular composition of photoreceptor synapses [74]; thus, it is not yet clear how variations in the latter two genes might be connected to tumorigenesis.

WHO classifies meningioma aggressiveness according to each meningioma's histological features, mitotic index and capacity to invade the brain, with no particular reference to molecular features [75]. In this work, we identified potential biomarkers for aggressive characteristics, including Grade I, recurrence, metastasis and fast growth in culture. Although Jed38_MN was initially classified as a Grade I tumour, subsequent work indicated a high Ki67 count and cancer stem cell features [23]. In addition, Jed18_MN was observed to reoccur following the sample collection; thus, it is potentially more aggressive than Jed13_MN. The differential pathway analysis indicated a particular focus on cell adhesion and rRNA catabolic processes. This could indicate that aggressive cells have alterations in their plasma membrane adhesion molecules, which then generate deformed adhesion proteins that enable easier cell-cell detachment and allow more movement and invasion [76]. The increased rRNA breakdown rate could be a response to the high rate of faulty mRNA production, perhaps via the non-functional rRNA decay (NRD) pathway [77]. Both pathways can potentially be utilized for diagnostic, prognostic and therapeutic applications, differentiating between slow- and fast-growing meningioma tumours.

The most frequently identified variant in the aggressive samples was detected in the NCOA3 gene. NCOA3, which is also referred to as SRC3 or AIB1, is a master nuclear coactivator with the potential to interact with 129 proteins, enhance the transcriptional activation of 19 pathways, and influence at least 13 biological processes [78]. Of particular interest is the implication of NCOA3 in pluripotency maintenance, as gene knockdown compromises the expression of several pluripotency markers, including Nanog, Oct4 and Sox2, and impairs the differentiation potential of mouse ESCs [79]. This gene has also been implicated in regulating the synapse formation and plasticity of neuronal cells [80]. Gene deregulation has been reported in several cancers, including astrocytoma, endometrial carcinoma, and breast, pancreatic and prostate cancers, and high expression has been associated with resistance to tamoxifen and poor prognosis [81, 82, 83, 84, 85, 86]. Concordantly, NCOA3 was also overexpressed in our tissue samples compared to the protein expression levels in the meninges. Although protein expression is not differential between Grade I and aggressive tumours, it appears that the BTC variant in the poly-Q region is differential. This domain is important for the protein's acetylation activity and protein–protein interactions [29]. The variant, COSM1483713, has been detected in several tumours, including ER-positive breast carcinoma [87], stomach and large intestine adenocarcinoma [88], small cell lung carcinoma, squamous cell carcinoma [89] and skin cancer [90]. However, none of the aforementioned studies looked into the functional effects of this variant. The increased lengths of the poly-Q region were found to be associated with a higher risk of familial breast cancer in carriers of mutations in exon 11 of the BRCA2 gene [91], while shorter lengths of the region were associated with increased risk of sporadic cancers [29], suggesting that subtle changes in this region impact tumours variably. Our investigation into the possible functional effects of the variant in meningioma led to the identification of 15 potentially affected genes. The 11 genes that were upregulated have shown to be overexpressed in several cancer tissues [92]. Of particular interest is the SEMA3D gene, a stem cell-related gene known to promote cell proliferation and neural crest cell development [93] that was found to contribute to perineural invasion and the metastasis of orthotopic pancreatic tumours in mice [94]. This gene is predicted to be down regulated by nuclear receptor coactivator 2 (NCOA2) [95], a protein that is thought to interact with NCOA3. Another intriguing deregulation is the overexpression of VEGFA, an important protein in angiogenesis. The knockdown of NCOA3 has been shown to attenuate the induction of VEGFA in MCF7 cells [96]. The precise mechanisms by which the poly-Q mutation exerts its effects to deregulate the identified genes remain to be clarified.

A limitation of this study has been the patient low sample size. However, to our knowledge, no published work previously incorporated data from exome sequencing of 14 primary meningioma cell lines. Primary cell lines especially for low grade meningiomas are difficult to induce in culture due to their slow growth and inefficient adaptation. The most commonly studied available meningioma cell line IOMM-Lee, was culturally transformed in order to encourage growth [17]. In addition, our samples included a recurrent and a metastatic sample, which are rare and constitute less than 20% of all meningiomas [3]. Future work will need to incorporate data from a larger cohort of meningioma primary cell lines.

Conclusion

This study presented novel predictive, tumorigenic and progressive variants identified using whole-exome sequencing of BTCs of patients with intracranial meningiomas. The identified BTC variants and their associated pathways were particularly inclined to be stem cell-associated, with the most frequent variant occurring in the GPRIN2 gene and the top pathway being ‘cell morphogenesis involved in differentiation’. Transforming variants seen in patients who had aggressive tumours, such as those in the stem cell-related gene NCOA3, were identified. Analysis in relation to the expression and potential target genes of the NCOA3 gene confirmed a link to stem cell and angiogenesis markers and presented NCOA3 as a strong potential biomarker for aggressive meningiomas. Together, the reported variants in the meningioma samples emphasized a link between stem cells' genetic predisposition and the development of meningiomas, and revealed an intricate and patient-related heterogeneity in meningiomas. More studies are needed to provide functional models that enable a full understanding of the impacts of the stem cell-detected variants on tumorigenesis.

Declarations

Author contribution statement

Deema Hussein, Hans-Juergen Schulten: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Ashraf Dallol, Rita Quintas: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Mona Alomari, Ishaq Khan: Performed the experiments; Wrote the paper.

Saleh Baeesa, Mohammed Bangash, Fahad Alghamdi: Conceived and designed the experiments; Performed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

M-Zaki Mustafa ElAssouli: Performed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Mohamad Saka: Conceived and designed the experiments; Performed the experiments; Wrote the paper.

Angel Carracedo, Adeel Chaudhary, Adel Abuzenadah: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This work was supported by the Dean of Scientific Research, , KSA (HiCi 1434-117-11).

Data availability statement

Exome sequencing data associated with this study has been deposited at NCBI SRA under the accession number PRJNA630560. Gene expression data associated with this study has been deposited at NCBI's Gene Expression Omnibus under the accession number GSE77259.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.