Somatic mosaicism in the MAPK pathway in sporadic brain arteriovenous malformation and association with phenotype.

All patients were prospectively enrolled in the University of California, San Francisco (UCSF) BAVM Project between 2002 and 2017 and had no family history of BAVM, hereditary hemorrhagic telangiectasia, or other known genetic syndromes. Diagnosis of BAVM was confirmed with catheter angiography and/or pathology evaluation. Informed consent for genetic studies was obtained, and the study was approved by the Committee on Human Research at UCSF. For the discovery cohort, we analyzed data from 14 White patients with a diagnosis of sporadic BAVM and both BAVM tissue that had been frozen in optimal cutting temperature (OCT) compound and blood available for analysis. All patients who underwent microsurgical resection of BAVM had not received prior interventional treatment, except 1 patient whose first treatment was embolization. The patients were evenly split in terms of sex and between those with ruptured and unruptured BAVMs. All patients were White, the most common race/ethnicity in the UCSF BAVM cohort, to minimize confounding due to genetic differences. For the replication cohort, we analyzed an additional 56 patients with sporadic BAVM, of whom 8 had undergone embolization and 1 had undergone radiosurgery prior to microsurgical resection. The replication cohort included patients with BAVM with full clinical data, and sufficient OCT-frozen tissue available for analysis. Thus, the replication cohort included patients with multiple ethnic backgrounds and is more representative of the patients included in the overall UCSF BAVM Project cohort. Of 70 patients included in the entire cohort, 43 had ICH. Demographic and phenotypic data for all patients are shown in Supplemental Table 1 . Blood genomic DNA was previously isolated and banked with standard DNA extraction techniques. For BAVM tissue, samples were snap-frozen in liquid nitrogen after resection and banked at −80°C in OCT compound. Frozen tissue was cryosectioned in bulk, and genomic DNA was isolated with the DNeasy Blood & Tissue Kit (QIAGEN). We used Nextera Rapid Capture Exome kits (Illumina) for exome enrichment and library preparation, with 250 ng DNA per sample. Paired-end sequencing (2 × 76 bp) was performed on the Illumina HiSeq 2500 platform at the Broad Institute, Cambridge, Massachusetts. Sequencing reads were mapped to the human reference genome assembly (hg19) with the Burrows-Wheeler Alignment tool, 25 then preprocessed according to the Genome Analysis Toolkit (GATK) best practice workflows for somatic short-variant discovery. 26 , 27 Mutect2 in GATK version 3.6 was used for calling somatic single-nucleotide variants (SNVs) with default parameters. We devised a customized workflow to prioritize candidate somatic variants associated with BAVM, with filters applied sequentially ( Supplemental Fig. 1 and Supplemental Table 2 ). Because the TLOD (theta logarithm of the odds) values (i.e., initial LOD threshold for calling tumor variant) of both KRAS mutations (G12V and G12D) were greater than 4.2, we used TLOD ≥ 4.2 as the first filter. Finally, SnpEff 28 and PolyPhen-2 29 (HumDiv-trained) were used for variant annotation and functional effect prediction, respectively. Variants predicted to be “probably damaging” were retained. We applied MutSigCV 30 and MUFFINN 31 to identify significantly mutated genes harboring mutations associated with BAVM. Because MutSigCV requires both silent and nonsilent variants to evaluate the background mutation rate, we used the 12,272 SNVs that passed the first four filters of the analysis ( Supplemental Fig. 1 ). Oncotator 32 was used to generate the input mutation file required by MutSigCV. For MUFFINN, we used the 3624 SNVs that passed all filters of the analysis ( Supplemental Fig. 1 ), network algorithm NDsum, and functional network HumanNet. We used ExomeCQA 33 to assess the coverage distributions of the sequencing reads in the study cohort, and we used bedtools to calculate the coverage files required. Cohort Coverage Sparseness was used to assess the coverage of all exons across the genome, and Unevenness was used to evaluate the coverage of a given exon across all samples. For validation and replication, we assayed KRAS G12V and G12D by using commercially available digital droplet polymerase chain reaction (ddPCR) assays according to the manufacturer’s instructions (Bio-Rad Unique Assay ID dHsaCP2000005 for mutant G12V, dHsaCP2000001 for mutant G12D, dHsaCP2000006 for wild-type G12V, and dHsaCP2000002 for wild-type G12D). The QX100 ddPCR system (Bio-Rad) was used for droplet generation and data acquisition. A total of 40 ng DNA per reaction was used for a maximum sensitivity of 0.025%. We included a positive control (somatic mutation positive), negative control (wild type only), and nontemplate control on all plates. Data analysis included determination of allelic burden with Bio-Rad QuantaSoft Analysis Pro software. All assays were performed in duplicate or triplicate. Furthermore, we required a minimum of 250 signal-positive droplets per patient and 100 signal-positive droplets per well. We included samples in the statistical analysis if they passed quality control for both ddPCR assays for G12D and G12V. We calculated summary statistics for patient and BAVM characteristics as mean ± standard deviation (SD) for numeric variables or as count and percentage for categorical variables. For our primary analysis, we tested whether the percentages of G12D or G12V alleles were associated with three outcomes: 1) age at tissue collection; 2) BAVM lesion size (maximum lesion diameter in millimeters); and 3) age at first ICH. For outcomes 1 and 2, we used linear regression models, adjusted for sex and prior ICH, and report exponentiated regression coefficients, or proportional increase (PI), with 95% confidence intervals (CIs). For outcome 3, we used a Cox proportional hazards model, adjusted for sex and BAVM size, and report hazard ratio (HR) and 95% CI; patients without ICH at age of tissue collection were censored. Secondary analyses evaluated individually G12D and G12V allelic burden as a percentage, and we also dichotomized patients according to the presence or absence of either somatic mutation (allelic burden threshold of 0.5%). Data analysis was conducted with Stata version 15.1 (StataCorp LLC).

Means of 123 million and 119 million reads were generated for the BAVM lesion and blood DNA samples, respectively. Averages of 99.04% and 98.80% of reads were mapped for the lesion and blood DNA samples, respectively. The average (range) sequencing coverage was 108× (88× to 142×) for lesion DNA and 97× (70× to 125×) for blood DNA ( Supplemental Table 3 ). We observed variations in coverage across all samples ( Supplemental Table 3 ) and over the targeted regions within each sample (data not shown). A total of 37,591 somatic variants were identified exome-wide (by using the default parameters of Mutect2). After application of sequential filters ( Supplemental Fig. 1 ), the final number of candidate somatic variants decreased to 3624. We cross-referenced the list of candidate somatic mutations with 295 MAPK pathway genes (included in the Kyoto Encyclopedia of Genes and Genomes database, https://www.genome.jp/kegg/ ). KRAS was the only MAPK pathway gene with a recurrent mutation (G12V) that was observed in our study ( Fig. 1 ). There were 10 other MAPK pathway genes with candidate somatic variants, each variant unique to 1 patient ( Fig. 1 ). A total of 22 “probably damaging” (PolyPhen2) variants were identified in these 10 MAPK genes (2 or 3 variants per gene) ( Fig. 2 ). Only 1 variant had allelic frequency greater than 6%. Each of the two previously reported recurrent KRAS mutations 16 – 21 was identified in two different BAVM samples by utilizing Mutect2. The KRAS G12V mutation was among the final group of variants that passed all filters and was identified as a recurrent mutation in our analysis. The KRAS G12D mutation was predicted as “possibly damaging” and was filtered out with functional effect prediction analysis. When we counted the number of sequencing reads that contained these mutations without using any filters, G12V was found in 3 of 14 patients and G12D in 5 patients ( Fig. 3A and B ); all positive results were subsequently validated with ddPCR ( Fig. 3E and F ). No samples had positive results for both KRAS mutations. The two KRAS mutations were not detected in any blood samples. No other candidate somatic mutations were identified in KRAS . We applied two different approaches, MutSigCV and MUFFINN, to evaluate whether the rate of identified somatic mutations in candidate genes was higher than the expected background mutation rate. Exome-wide, a total of 218 genes were identified as significantly mutated (p < 0.1) with MutSigCV. Six of these were MAPK signaling pathway genes ( ATF2 , CRKL , DUSP6 , KRAS , PDGFB , and PDGFRB) . KRAS was the only MAPK pathway gene with recurrent somatic mutation, whereas CRKL and PDGFRB had different mutations in at least 2 patients. By using MUFFINN, we allocated a score based on the number of candidate mutations and known functional network information to 458 genes. Eleven of these were MAPK pathway genes, including KRAS and all genes listed in Fig. 1 . Overall, KRAS , CRKL , and PDGFRB were the only MAPK pathway genes identified as bearing candidate somatic variants associated with BAVM according to both MutSigCV and MUFFINN. Among all MAPK pathway genes, KRAS and PDGFRB were the top two candidate genes containing BAVM-associated variants. MutSigCV analysis identified KRAS as the most significantly mutated gene, followed by PDGFRB ; MUFFINN analysis assigned the highest score to PDGFRB , followed by KRAS . By combining the discovery and replication cohorts, we used ddPCR to assay KRAS G12D and G12V in 69 BAVM lesion samples and 32 paired blood samples ( Supplemental Table 1 ). Fifty-six of 69 (81%) samples passed quality control and were included in the statistical analysis. The mean age ± SD at tissue collection was 37.1 ± 17.7 years, and 28 (50%) patients were female. The mean BAVM diameter was 19.2 ± 10.7 mm. Thirty-five (63%) patients had ICH prior to resection at a mean age of 35.7 ± 19.7 years. G12V or G12D somatic mutation was detected in 31 (55%) BAVM DNA samples and in zero blood DNA samples. The somatic mutation allelic burden of KRAS G12D or G12V exceeded 0.5% in 18 (32%) BAVM samples; the highest level detected was 5.0%. We evaluated the associations between allelic burden of somatic KRAS mutation and three phenotypes of BAVM severity: age at resection, lesion size (maximum diameter), and age at first ICH ( Table 1 ). We found no evidence of an association between somatic KRAS mutation burden and age at BAVM resection (PI 0.7 [95% CI –2.9 to 4.3], p = 0.69), BAVM size (PI 1.4 [95% CI 0.7–3.6], p = 0.18), or age at first ICH (HR 0.87 [95% CI 0.64–1.20], p = 0.40). Individual G12V or G12D mutation burden, or presence of each mutation, was also not significantly associated with any BAVM phenotype ( Table 1 ).

We used WES to evaluate paired lesion and blood DNA samples from 14 patients with BAVM for somatic mutations in 295 MAPK pathway genes. Two KRAS mutations (G12V and G12D), which were associated with BAVM in recent studies, were also identified in our study cohort and further validated in a larger cohort with ddPCR. By using these two mutations as a standard and applying a pipeline of stringent filters, we identified candidate BAVM-associated somatic mutations in 10 other MAPK pathway genes. Overall, KRAS mutations (G12V and G12D) have been identified in at least seven BAVM cohorts, 16 – 21 including ours, which is one of the largest cohorts published to date. Nikolaev et al. 16 first reported that 60% of patients with sporadic BAVM carried these KRAS mutations, whereas Hong et al. 17 reported that 90% of patients had positive results for KRAS mutation. In our study, we found somatic KRAS mutations in 43% of patients with BAVMs. The lower percentage in our study may be due to, in part, the lower sequencing coverage (70×) of our WES data compared with those of other studies (> 1000× 17 or 300× 16 ). However, the reported allelic burden of somatic KRAS mutations was similar in all studies (as high as approximately 5%), including ours. An interesting observation across multiple studies is that G12D was found in more patients with BAVM than G12V, but G12V had a greater average allelic burden (2.01%) than G12D (1.14%). This suggests that a lower allelic burden of the KRAS G12D mutation may be sufficient to produce the BAVM phenotype in comparison with that of G12V, and that G12D may have a stronger functional impact on the KRAS protein in the cerebrovascular context. Functional studies of these KRAS mutations, mostly in epithelial cancer cell lines, suggest that KRAS G12V and G12D mutations may have differential effects on proliferation, selection, and downstream signaling. 34 , 35 Across the seven studies on sporadic BAVM that have been conducted to date, KRAS is the only gene that has been consistently found with recurrent somatic mutations. 17 – 21 This implies that KRAS may play an important role in the etiology of BAVM and may be a potential therapeutic target in the majority of patients with sporadic BAVM. KRAS mutations are commonly found in patients with many types of cancer (including colorectal, pancreatic, and lung cancer 36 ) and other angiogenesis-related disorders such as endometriosis. 37 We also evaluated the associations between somatic KRAS mutations (allelic burden or presence) and three phenotypes of BAVM severity (age at resection, lesion size, or age at first hemorrhage) that we hypothesized could be affected by somatic mutation burden; however, we did not detect any statistically significant associations. A potential explanation for the lack of association may be confounding by the variable presence of nonlesion cells in the grossly cryosectioned samples. A more precise evaluation of allelic burden in BAVM lesion cells may be necessary to detect an association: for example, the endothelial cells of BAVM lesions may need to be isolated first with laser capture microscopy. The lack of a statistically significant association between somatic KRAS mutation allelic burden and phenotype could be due to low precision because the 95% CIs are wide ( Table 1 ). For example, the 95% CIs for age at tissue collection and presence of a mutation all have spans greater than 20 years. We identified 10 other genes in the MAPK pathway with candidate somatic mutations that were present in multiple patients included in the discovery cohort. Among these, PDGFRA and PDGFRB have differential roles in angiogenesis and vessel stability. 38 PDGFRB was one of the three genes (alongside KRAS and CRKL ) with candidate somatic mutations that was detected more often than expected. Experimental animal models of BAVM show reduced Pdgfr-β expression, 39 and a significant reduction in PDGFRB-positive pericyte coverage has been observed in BAVM vessels and adjacent perivascular space, suggesting a potential mechanism by which damaging somatic mutations in PDGFRB may contribute to BAVM development. 40 Recently, activating PDGFRB somatic mutations have also been identified in fusiform (but not saccular) cerebral aneurysms, further supporting a role in cerebrovascular anomaly development. 41 CRKL , the third gene with an excess of candidate somatic mutations in BAVM lesions, is a putative proto-oncogene that is upregulated in multiple types of cancers. 42 In vitro, amplification of CRKL in EGFR -mutant cells induces drug resistance via activation of ERK and AKT signaling, 43 whereas CRKL overexpression promotes cell invasion via upregulation of MMP-9 expression. 44 This study has several limitations. Subfractionation was not used to isolate BAVM lesion cells or endothelial cells, and this may have reduced the amount of mutant DNA molecules in the sample. As discussed above, the relatively low sequencing coverage of our original WES sequencing data set (70× compared with 300× to > 1000× in other studies 16 , 17 ) reduced sensitivity to detect candidate mutations with allelic representation < 2%. We also observed variation in coverage across samples ( Supplemental Table 2 and Supplemental Fig. 2 ) and uneven coverage of reads over the targeted regions within samples (data not shown), and these further reduced the detection and the accuracy of variant calling. BAVM is a complex lesion with likely genetic heterogeneity. In all previous studies, as well as in ours, somatic KRAS mutations were not found in every patient. We identified several other candidate mutations in MAPK pathway genes; in some cases, somatic mutations in multiple genes were identified in a single patient. Although these mutations remain to be validated, they suggest that a single mutation may be insufficient to induce BAVM formation, and instead multiple hits to genes in the same pathway may be required. These gene hits may come from other mechanisms, e.g., epigenetic inactivation, or BAVM may be related to aberrations in other genes and pathways. This observed genetic heterogeneity ultimately suggests that different therapeutic approaches along with reliable biomarkers may be needed for different patients. In the future, diagnostic sequencing of patient samples (e.g., via targeted sequencing of cell-free DNA or endovascular sampling of lesion cells 45 ) could be used to match patient subgroups with appropriate targeted therapeutic strategies.

Our study suggests that BAVM exhibits somatic mutational heterogeneity, because different combinations of rare somatic mutations in MAPK pathway genes were observed in BAVM lesions (but not in blood DNA). In our study, we replicated the high prevalence of two previously identified KRAS mutations (G12V and G12D), confirming that these two mutations are the most prevalent somatic mutations present in BAVM tissue. We did not find a significant association between somatic KRAS mutation burden and three phenotypes of BAVM severity, but more sensitive detection of mutation burden may be needed. We also identified several new candidate somatic mutations associated with sporadic BAVM, including mutations in PDGFRB and CRKL . Validation and functional studies of these mutations are needed to elucidate the biological processes leading to BAVM.