The molecular landscape of glioma in patients with Neurofibromatosis 1

Abstract

Neurofibromatosis type 1 (NF1) is a common tumor predisposition syndrome in which glioma is one of the prevalent tumors. Gliomagenesis in NF1 results in a heterogeneous spectrum of low- to high-grade neoplasms occurring during the entire lifespan of patients. The pattern of genetic and epigenetic alterations of glioma that develops in NF1 patients and the similarities with sporadic glioma remain unknown. Here, we present the molecular landscape of low- and high-grade gliomas in patients affected by NF1 (NF1-glioma). We found that the predisposing germline mutation of the NF1 gene was frequently converted to homozygosity and the somatic mutational load of NF1-glioma was influenced by age and grade. High-grade tumors harbored genetic alterations of TP53 and CDKN2A, frequent mutations of ATRX associated with Alternative Lengthening of Telomere, and were enriched in genetic alterations of transcription/chromatin regulation and PI3 kinase pathways. Low-grade tumors exhibited fewer mutations that were over-represented in genes of the MAP kinase pathway. Approximately 50% of low-grade NF1-gliomas displayed an immune signature, T lymphocyte infiltrates, and increased neo-antigen load. DNA methylation assigned NF1-glioma to LGm6, a poorly defined Isocitrate Dehydrogenase 1 wild-type subgroup enriched with ATRX mutations. Thus, the profiling of NF1-glioma defined a distinct landscape that recapitulates a subset of sporadic tumors.

Extended Data Fig. 1 |. Data analysis workflow.

Fifty nine tumor samples from 56 NF1-glioma patients with 43 matched normal were profiled with WES, DNA Methylation profiles (31 tumors) and RNA sequencing (29 tumors). WES was used to call NF1 germline mutations using HaplotypeCaller and Somatic-germline log odds filter. Somatic SNVs were called from WES data by integrating the results of five algorithms (Freebayes, MuTect, Strelka, VarDict and VarScan). Recurrent CNVs were detected by GATK and GISTIC2. SNVs and CNVs were validated by Sanger sequencing (93% validation rate) and genomic qPCR (96% validation rate), respectively. Neoantigen prediction was obtained using netMHCpan and HLA genotype was determined by Polysolver, Optitype, Phlat and Seq2hla and validated by affinity binding kinetics. COSMIC cancer mutation signatures were identified by deconstructSig and compared to those occurring in sporadic glioma. DNA Methylation arrays were used to classify NF1 glioma in the methylation subtypes of sporadic glioma form the TCGA pan-glioma dataset (KNN). RNAseq was used to define gene expression clusters and immune subtypes of low-grade NF1-glioma and results were confirmed by RT-qPCR and immunohistochemistry. Integrative analysis of gene expression and DNA Methylation identified epigenetic signatures characterizing immune subtypes of low-grade glioma. A pan-glioma gene regulatory network was used to identify MRs of the ATRX-mutant phenotype in LGm6 sporadic and NF1-glioma (RGBM). Finally, the impact of ATRX mutation on survival was assessed using TCGA pan-glioma and NF1-glioma data.

Extended Data Fig. 2 |. Fingerprint analysis of WES NF1 samples.

Dendrogram of hierarchical clustering of 59 tumor and 43 normal samples based on Pearson correlation coefficients of SNPs allele fractions. Case ID and the tissue specimen are indicated (blood DNA, red; tumor with available matched blood DNA, blue; tumor without matched normal DNA, yellow). The analysis confirmed proper matching of samples for each of the 43 tumor-blood DNA pairs. Thirteen tumors without available paired normal DNA (yellow) showed individual branches in the clustering dendrogram.

Extended Data Fig. 3 |. Validation of recurrent CNVs.

Genomic qPCR was performed to assay copy number changes for TERT (n = 10 glioma samples), b, IL-15 (n = 8 glioma samples), c, FGF1 (n = 17 glioma samples) and d, CDKN2A (n = 11 glioma samples). Red and blue bars indicate WES-inferred gene gain and loss, respectively. Analysis of normal DNA (green bars) was included to define diploidy (dotted line). Tumor samples diploid for the tested gene were included as control (white bars). Bar graphs show mean± s.d. of 3 technical replicates. Experiments were repeated three times with similar results.

Extended Data Fig. 4 |. Somatic mutation burden of NF1-glioma and pediatric and adult cancer genomes.

Distribution of somatic non-synonymous coding mutation rate is represented on a logarithmic scale for NF1- and sporadic glioma (bold) and other frequent cancer types, including pediatric tumors. Cancer types and subgroups are ordered by increasing mutation frequency median, with the lowest frequencies (left) found in pediatric tumors and low-grade NF1-glioma. Somatic mutations used to calculate the mutational burden for different cancer types were retrieved from TCGA (adult tumors) and TARGET (pediatric tumors) databases.

Extended Data Fig. 5 |. Mutational clonality.

Analysis of mutational clonality in 55 NF1-glioma samples. a, Number of mutation clones relative to age (Pearson correlation coefficient = −0.126 and p = 0.363), and b, tumor grade (Pearson correlation coefficient = 0.031 and P = 0.820). Blue line: linear regression; shaded area: 95% confidence interval.

Extended Data Fig. 6 |. Analysis of DNA Copy Number Variations.

Schematics of chromosome location peaks (gain, red; loss, blue) identified using GISTIC2. Peaks are designated by candidate targets for each region, selected according to criteria described in Methods. The complete list of chromosome location peaks is included in Supplementary Table 6a, b.

Extended Data Fig. 7 |. Mutual exclusivity and co-occurrence of genetic alterations in NF1-glioma.

a, Mutually exclusive and b, co-occurring genetic alterations in NF1-glioma were evaluated using CoMEt and two-sided Fisher’s exact test, respectively. Significant mutual relationships between two gene alterations are indicated by a line (green, exclusion; red, co-occurrence) whose thickness represents −log₁₀ of p-value (reported in Supplementary Table 7).

Extended Data Fig. 8 |. Distribution of somatic mutation spectrum in NF1-glioma.

Dirichlet multinomial regression test for ATRX status (n = 10 and n = 46 ATRX mutant and ATRX wild-type samples, respectively), age (n = 22 pediatric glioma; n = 33 adult glioma) and glioma grade (n = 24 high-grade glioma; n = 32 low-grade glioma). b, The relative proportions of the six different possible base-pair substitutions are represented by barplots for ATRX mutant (n = 10, solid fill) and ATRX wild-type (n = 46, patterned fill). The relative frequency of C > T transition was significantly higher in ATRX mutant tumors (p = 5.1 × 10–3, two-sided Fisher’s exact test).

Extended Data Fig. 9 |. Somatic alterations in PI3K and Transcription/Chromatin regulation pathways in NF1-glioma.

Integrated matrix of 59 NF1-glioma samples (56 patients) and somatic alterations (SNVs and indels, and significant copy number variations) occurring in genes linked to PI3K and transcription/chromatin regulation pathways (left panel, high-grade glioma; right panels low-grade glioma). Rows and columns represent genes and tumor samples, respectively. NF1-glioma samples are sorted in the same order of Fig. 2. Genes are grouped by PI3K (purple) and transcription/chromatin regulation (blue) pathways. Genomic alterations, age, the histology of glioma and the identification of NF1 germline mutation are shown by the indicated colors. Validation by Sanger sequencing (SNVs) and quantitative-genomic PCR (gains and losses) are indicated by yellow and green triangles, respectively.

Extended Data Fig. 10 |. Somatic alterations in splicing, MAPK and cilium/centrosome pathways in NF1-glioma.

Integrated matrix of 59 NFl-glioma (56 patients) and somatic alterations (SNVs and indels, and significant copy number variations) occurring in genes included in splicing, MAPK and cilium/centrosome pathways (left panel, high-grade glioma; right panels low-grade glioma). Rows and columns represent genes and tumor samples, respectively. NFl-glioma samples are sorted in the same order of Fig. 2. Genes are grouped by splicing (red), MAPK (yellow) and cilium/centrosome (green) pathways. Genomic alterations, age, the histology of glioma and the identification of NF1 germline mutation are shown by color as indicated. Validation by Sanger sequencing (SNVs) and quantitative-genomic PCR (gains and losses) are indicated by yellow and green triangles, respectively.

Fig. 1 |. Analysis of germline and somatic mutations in NF1-glioma patients.

a, The relative frequency distribution of age at diagnosis is represented by density plot: the overall distribution of NF1-glioma patients (dashed black line, n = 55) by age identifies two peaks, 13.5 and 38.8 years. Low-grade gliomas (blue line, n = 32) occur more frequently in children, while high-grade gliomas (red line, n = 23) are diagnosed more frequently in adults. b, Germline mutations in the NF1 gene identified in NFI-glioma patients by WES. The spectrum of NF1 germline variants (SNVs and indels) is represented with each mutation shown only once per patient. We identified NF1 germline pathogenic mutation in 51 of 56 (91%) patients. Patients no. 47 and no. 52 had one additional pathogenic germline mutation. c, Scatter plot showing the number of somatic mutations (SNVs and indels) occurring in low-grade and high-grade NFI-glioma (low-grade glioma, n = 32; high-grade glioma, n = 24; P = 7.4×10⁻⁵, two-sided Mann-Whitney-Wilcoxon (MWW) test). d, Scatter plot showing the number of somatic mutations occurring in patients younger than 16 years (pediatric, n = 22) or older than 18 years (adult, n = 33; P = 9.8×10⁻⁴, two-sided MWW test). Mutations in the patient of unknown age were not included in the analysis. e, Scatter plot showing the number of mutations according to grade and age (low-grade glioma, pediatric, n = 17; high-grade glioma, pediatric, n = 5; low-grade glioma, adult, n = 15; high-grade glioma, adult, n = 18). Mutations in the patient of unknown age were not included in the analysis. P = 3.8×10⁻³, high-grade versus low-grade adult gliomas; P = 0.025, adult versus pediatric low-grade gliomas; P = 0.06, low-grade versus high-grade pediatric gliomas; P = 0.07, pediatric versus adult high-grade gliomas (two-sided MWW test). Scatter plots show median and interquartile range. Median and range of mutations are reported below each plot.

Fig. 2 |. Landscape of somatic genomic alterations in NF1-glioma.

a, Integrated matrix of 59 glioma samples from 56 patients and gene variants (SNVs, indels, and significant CNVs) observed in NF1-glioma (left panel, high-grade glioma; right panels, low-grade glioma). Rows and columns represent genes and tumor samples, respectively. Genomic alterations, age, the histology of glioma, and NF1 germline mutations are indicated. NFI-glioma samples are sorted by their mutation profiles, except for patient no. 5, hypermutated high-grade glioma, and patient no. 39, including four spatially distinct glioma samples, which are shown at the last columns of left and right panel, respectively. Recurrently mutated genes are selected for their previously established association with glioma (ATRX, CDKN2A, TP53, PIK3CA/B, PTEN, BRAF, FGFR1 and FGF1, PRKCA, TERT), cancer biology (DOCK2/3/6, FDZ3/8, BCL9/9L, TOP2/3B), and immune functions (IL15, DGKQ). Genes are sorted according to higher frequency (percentage of patients) in high-grade (top, red) or low-grade gliomas (bottom, blue), respectively. Validations by Sanger sequencing (SNVs) and quantitative-genomic PCR (gains and losses) are indicated by yellow and green triangles, respectively. LOH, loss of heterozygosity. b, Function/pathway analysis of damaging somatic mutations and CNVs. Genetic alterations in NF1-gliomas grouped into PI3K, transcription/chromatin regulation, splicing, MAPK, and cilium/centrosome functions. A significantly higher frequency of genetic alterations in PI3K, transcription/chromatin regulation, and splicing pathway were observed in high-grade glioma (n = 24; P = 4.7×10⁻⁵, P = 9.1×10⁻⁴, and P = 0.03, respectively; two-sided Fisher’s exact test), while mutations in the MAPK pathway were more frequent in low-grade glioma (n = 32; P = 0.03, two-sided Fisher’s exact test). The integrated matrices of NF1-glioma and gene pathway alterations are reported in Extended Data Figs. 9 and 10.

Fig. 3 |. Analysis of ATRX somatic mutations in NF1-glioma patients.

a, ATRX mutations were identified by WES. The spectrum of ATRX somatic variants (SNVs and indels) is represented with each mutation shown only once per patient. We identified and validated by Sanger sequencing ATRX pathogenic mutations in 10 patients (1 low-grade glioma, 3%; 9 high-grade gliomas, 37.5%). b, The relative frequency of age distribution is represented for all patients (dashed black line, n = 55), ATRX wild-type (green line, n = 46), and ATRX mutant gliomas (red line, n = 9; for one patient carrying ATRX mutation age was unknown and was not included in the analysis). c, Microphotographs of ATRX immunohistochemistry in gliomas from NF1 patients. Representative images are from n = 7 low-grade gliomas (left) and n = 16 high-grade gliomas (right). Results were validated on more than ten independent samples to ensure the staining pattern on human tissue was reproducible. High-grade glioma samples were negative for ATRX expression whereas low-grade gliomas retained ATRX protein expression. d, Contingency table shows loss of ATRX protein expression in 8 of 16 high-grade and in none of 7 low-grade NF1-gliomas (P= 0.05, two-sided Fisher’s exact test). e, C-circle (CC) assay was performed to measure ALT activity in NF1-glioma samples. The scatter plot reports the normalized CC content for each glioma according to ATRX mutational status: ATRX wild-type (blue, n = 11) and ATRX mutant gliomas (red, n = 10). For each group the median with interquartile range is indicated. All ATRX mutant gliomas but only one ATRX wild-type glioma showed increased ALT activity (normalized CC content greater than 1; P = 2.3×10⁻⁵, two-sided MWW test).

Fig. 4 |. Transcriptomic analysis of NF1-glioma.

a, Consensus clustering on the Euclidean distance matrix based on the most variable genes among 29 NF1-glioma samples (1,330 genes). The consensus matrix is obtained from 10,000 random samplings using 70% of the 29 samples. The 10 high-grade samples fall in 1 cluster (red) and all low-grade samples (n = 19) fall in a different cluster (blue). b, Hierarchical clustering of 29 NF1-gliomas by Euclidean distance with the Ward linkage method was based on the 100 most differentially expressed genes (two-sided MWW test, top and bottom 50 genes). c,d, Enrichment map network of statistically significant gene ontology categories in (c) ten high-grade and (d) 19 low-grade NF1-gliomas (two-sided MWW-GST; q < 0.001, absolute NES >0.6). Nodes represent gene ontology terms and lines their connectivity. Node size is proportional to number of genes in the gene ontology category and line thickness indicates the fraction of genes shared between groups. Gene network categories in NF1 high-grade gliomas are linked to mitotic progression, chromosome organization, and RNA biogenesis/regulation. Gene network categories in NF1 low-grade gliomas converge on proinflammatory immune response enriched for T lymphocyte effector functions. e, Unsupervised clustering of single-sample MWW-GST enrichments of the categories in d. Low-grade NF1-gliomas are divided into two clusters (red and green), characterized by high- and low-immune gene set enrichments, respectively (two-sided MWW test; only statistically significant categories are shown; the complete list is presented in Supplementary Table 8). f, Tumor purity scores of low-grade/high-immune, low-grade/low-immune, and high-grade tumors computed by ESTIMATE. The low-grade/high-immune group has significantly lower tumor purity when compared with either the low-grade/low-immune or the high-grade glioma groups (P = 0.001, high-immune versus low-immune low-grade gliomas; P = 0.01 high-immune low-grade gliomas versus high-grade gliomas; P = 0.762 for low-immune low-grade gliomas versus high-grade gliomas; two-sided t-test). g, Immune scores of low-grade/high-immune, low-grade/low-immune, and high-grade tumors computed by ESTIMATE (P = 9.6×10⁻⁵, high-immune versus low-immune low-grade gliomas; P = 0.005 high-immune low-grade gliomas versus high-grade gliomas; P = 0.852 for low-immune low-grade gliomas versus high-grade gliomas; two-sided t-test). h-l, Enrichments of CD8⁺ T cell functions in low-grade/high-immune compared with low-grade/low-immune and high-grade gliomas. Boxplots report the z-scores and P values (two-sided MWW test) for published CD8⁺ T cell signatures^–. Scatter plots show median, interquartile, and minimum to maximum range.

Fig. 5 |. T cell infiltration and neoantigen analysis in low-grade NF1-glioma subclusters.

a, Representative microphotographs of CD3 immunohistochemistry in low-grade/high-immune (left panels) and low-grade/low-immune (right panels) NF1-gliomas. Results were validated on more than ten independent samples to ensure that the staining pattern on human tissue was reproducible. b, The number of CD3-positive cells was scored in at least 5 pictures from low-grade/high-immune (n = 6, red dots) and low-grade/low-immune (n = 6, green dots) (*P = 0.003, two-sided t-test with Welch correction; scatter plots show mean with minimum to maximum range). c, Representative microphotographs of CD8 immunohistochemistry in low-grade/high-immune (left panels) and low-grade/low-immune (right panels) NFI-gliomas. Results were validated on more than ten independent samples to ensure that the staining pattern on human tissue was reproducible. d, The number of CD8-positive cells was scored in at least 5 pictures from low-grade/high-immune (n = 6, red dots) and low-grade/low-immune (n = 6, green dots) (*P = 0.017, two-sided t-test with Welch correction; scatter plots show mean with minimum to maximum range). e, Representative microphotographs of GZMB immunohistochemistry in low-grade/high-immune (left panels) and low-grade/low-immune (right panels) NFI-gliomas. Results were validated on more than ten independent samples to ensure that the staining pattern on human tissue was reproducible. f, The number of GZMB-positive cells was scored in at least 5 pictures from low-grade/high-immune (n = 6, red dots) and low-grade/low-immune (n = 6, green dots) (*P = 0.004, two-sided t-test with Welch correction; scatter plots show mean with minimum to maximum range). g, Quantification of neoantigens in low-grade NFI-glioma subclusters. The number of neoantigens per somatic mutation was significantly higher in the set of low-grade/high-immune NFI-gliomas (P = 0.034, two-sided MWW test). h, In vitro binding affinity kinetics of neoantigens and corresponding wild-type peptides for their restricted HLA class I allele. Data are shown as counts per second with increasing peptide concentration (log₁₀ M). Data are mean of two independent experiments. MT, mutant peptide; WT, wild-type peptide.

Fig. 6 |. NF1-gliomas resemble LGm6 subgroup of sporadic gliomas.

a, Heat map of DNA methylation data for the TCGA pan-glioma cohort (n = 819) and 31 NFI-gliomas according to the methylation clusters of sporadic gliomas. The methylation profiles of NFI-glioma samples were classified using a nearest neighbor classifier based on 1,233 cancer-specific DNA methylation probes. Thirty-one of 31 NFI-glioma samples were assigned to the LGm6 methylation cluster, one of the methylation clusters that includes both low-grade and high-grade gliomas. b, Oncoprint of selected somatic genomic alterations in the LGm6 group of gliomas from the TCGA data set (ATRX, TP53, CDKN2A, PTEN, PIK3CA, NF1, BRAF). Rows and columns represent genes and samples, respectively. Glioma grade was significantly associated with alterations of ATRX, TP53, CDKN2A, PTEN. Glioma grade IV, n = 40; glioma grade III, n = 12; glioma grade II, n = 13; P = 0.01, P = 0.02, P = 0.04, P = 0.002, respectively; two-sided Fisher’s exact test. c, Barplot of ATRX non-synonymous somatic mutations occurring in phenotypic subtypes of IDH wild-type gliomas (classic-like, mesenchymal-like, and LGm6) and LGm6 gliomas grouped by tumor grade. ATRX mutations were significantly enriched in grade III LGm6 (P = 0.01, two-sided Fisher’s exact test). d, Kaplan-Meier survival analysis of LGm6 gliomas stratified according to histological grade and ATRX status for grade III gliomas: grade II (green curve, n = 13), grade III ATRX mutant (blue curve, n = 5), grade III ATRX wild-type (cyan, n = 7), grade IV (red curve, n = 23). The ATRX mutant grade III subgroup showed a significantly worse survival when compared with ATRX wild-type grade III patients (P = 0.03, two-sided log rank test). No difference in clinical outcome was observed when comparing ATRX mutant grade III with grade IV. e, Master regulators (MRs) in ATRX mutant glioma. Gray curves represent the activity of each of the 10 MRs with the highest (red) or lowest (blue) activity. Red or blue lines indicate individual ATRX mutant samples displaying high or low activity, respectively, of the MRs in ATRX mutant compared with ATRX wild-type (n = 8 and n = 40 ATRX mutant and ATRX wild-type samples, respectively; P value, two-sided MWW test for differential activity (left) and mean of the activity (right)). f, Hierarchical clustering of MR activity in 48 high-grade LGm6 IDH wild-type gliomas (36 grade IV and 12 grade III). Data were obtained using the Euclidean distance and Ward linkage method built on differential activity of MRs in ATRX mutant (8 samples, red) versus ATRX wild-type (40 samples, black) tumors (two-sided MWW-GST q <0.01, absolute NES >0.6, and two-sided MWW test for differential activity P <0.01). The activity of 41 of 89 MRs was increased in ATRX mutant samples. g, Enrichment map network of statistically significant gene ontology categories (two-sided Fisher’s exact test q <0.01) for genes included in the regulons of the 10 MRs with the highest activity in ATRX mutant gliomas. Nodes represent gene ontology terms and lines their connectivity. Node size is proportional to number of genes in the gene ontology category and line thickness indicates the fraction of genes shared between groups.