Genotype correlates with clinical severity in PIK3CA-associated lymphatic malformations

Lymphatic malformations (LMs) are the most common vascular malformation of the head and neck affecting 1 in 4000 live births (1, 2). Previously known as cystic hygromas and lymphangiomas, isolated LMs are recognized as nonneoplastic malformations of lymphatic channels thought to arise from anomalous lymphovasculogenesis (1). They predominantly affect the head and neck (>75%) and may cause significant disfigurement and functional compromise, ranging from impaired speaking and swallowing to severe airway obstruction (3). The clinical course varies greatly. Large LMs are treated with multiple invasive, high-risk procedures including tracheotomy, while smaller neck LMs may resolve without invasive therapy (4). Lesion location and cystic structure are related to clinical severity (Figure 1). For surgeons, location above the hyoid bone distinguishes difficult-to-treat LMs in the oral-facial-pharyngeal region from those in the neck and body (5). The de Serres classification uses suprahyoid-infrahyoid location and laterality to create a 5-stage system that correlates well with posttreatment disease persistence (6). Radiographically determined cystic structures are divided into macrocystic (large cysts >2 cm diameter, cystic hygroma), microcystic (many small cysts, all <2 cm, lymphangioma), and mixed lesions (7). Macrocystic infrahyoid LMs are more responsive to therapy and more likely to resolve spontaneously, while orofacial microcystic lesions are highly infiltrative and persistent (5). Treatments are usually invasive, including sclerotherapy and surgery (8–10). However, recent discovery of mosaic PIK3CA mutations in LMs opened the possibility for targeted inhibition of the PI3K/AKT/mTOR pathway and prompted the use of sirolimus as a noninvasive therapeutic option (9, 11–15).

Figure 1

Clinical examples of lymphatic malformations. Preoperative photos and MRI imaging of individuals with lymphatic malformations of varying severity. Arrows indicate diseased area on MRI. (A) Clinical photo and coronal T2-weighted MRI image of infant with unilateral, macrocystic LM isolated to the right lateral neck. (B) Clinical photo and axial T1-weighted MRI of teenage male with unilateral microcystic LM of the left upper lip. Notably, this lesion is much more challenging to assess on MRI compared with a macrocystic lesion. (C) Clinical photo and axial T2-weighted MRI of bilateral microcystic LM of the anterior tongue. (D) Clinical photo and coronal T2-weighted MRI of an infant with a large, bilateral mixed cystic LM in a “beard distribution.” This severe level of disease burden markedly impacted this patient’s ability to breathe and he required a tracheostomy to maintain a safe airway. LM, lymphatic malformation; MRI, magnetic resonance imaging.

PIK3CA is the most widely expressed catalytic subunit of phosphatidylinositol 3-kinase (PI3K), an upstream element of the PI3K/AKT/mTOR pathway that regulates multiple cellular functions including proliferation, growth, and angiogenesis (16). As a proto-oncogene, PIK3CA mutations seen in cancer demonstrate significant gain of function (GOF) in vitro (17). The 3 most common mutations have been designated hotspot (HS) mutations (18–21). These 3 PIK3CA mutations are responsible for the majority of LMs and other PIK3CA-related overgrowth disorders (11). Prior studies of LMs have not examined correlations between clinical severity and specific mutation or level of mosaicism. Here we report genetic and clinical analysis of a large series of isolated LMs. We propose a metric that adjusts for differences in mutation “strength” — the genotype-adjusted variant allele fraction (GVAF) — which identified significant associations between mosaic genotype and disease severity.

Experimental design and cohort. We sequenced PIK3CA in 271 affected and unaffected tissue samples from 81 individuals with isolated LMs using single-molecule molecular inversion probes (smMIPs) (Figure 2A and Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.129884DS1). Tissue sources included 180 LM lesion samples, 14 LM cyst fluid samples, 7 blood samples, and 70 skin, mucosa, or nearby unaffected tissue samples. Individuals had an average of 3.3 total and 2.3 lesion samples, obtained either from distinct portions of the lesion or at separate surgeries (range: 1–13 samples per individual). Our cohort included 16 individuals with LMs isolated to the neck, 9 in the body, and 56 in the orofacial region; phenotype data are summarized in Table 1. Due to the relatively small number of isolated neck and body LMs, these individuals were combined into a single group for analysis purposes. Twenty-six had macrocystic LMs, 14 were microcystic, and 41 were mixed cystic LMs. Individuals in our cohort had 4.6 invasive procedures on average (range: 0–22). One patient had aspiration of cyst fluid only and was therefore considered to have zero invasive procedures. As expected from prior reports, macrocystic lesions were more commonly unilateral (P < 0.00001, Fisher) and present in the neck and body (P = 0.0002, Fisher) (Supplemental Figure 1A). Increased numbers of procedures were associated with orofacial location (P = 0.003, Fisher), bilateral lesions (P < 0.00001, Fisher), and microcystic structure (P < 0.0001, Fisher) (Supplemental Figure 1B).

Figure 2

Experimental workflow, smMIP coverage, RFLP results, and mutation results. (A) Experimental workflow for samples in our cohort. All samples initially went through an smMIP panel, after which 41 individuals had a mutation detected. The remaining 40 individuals underwent screening with 5 ddPCR assays as listed. Following ddPCR screening, 23 more individuals had a mutation detected and 17 remained unsolved. (B) Box plot of reads per smMIP for 271 samples and 42 smMIPs covering the entire coding region of PIK3CA. Boxes are defined by the first quartile inferiorly, median, and third quartile superiorly, with whiskers extending to the farthest non-outlier point (defined as within 1.5 times the interquartile range). Green indicates smMIPs in which a mutation was detected. The right y axis indicates the minimum possible VAF detectable according to read depth. A schematic of PIK3CA cDNA (NM_006218) is along the x axis aligned with the corresponding smMIPs. Mutations detected in our cohort are indicated beneath the cDNA schematic with an asterisk (*) to denote hotspot mutations. (C and D) smMIP and RFLP results for non-hotspot mutations not screened by ddPCR, (C) p.Gly106Val and (D) p.Gln546Lys. Variable mutation detection within an individual is demonstrated in D, as 2 lesion samples were tested but only 1 was positive by both smMIPs and RFLP. (E) Pie chart summarizing mutations detected in the cohort. ddPCR, droplet digital polymerase chain reaction; N/A, not applicable; NHS, non-hotspot mutations; RFLP, restriction fragment length polymorphism; smMIPs, single-molecule molecular inversion probes; VAF, variant allele fraction; WT, wild type.

Table 1

Summary of clinical characteristics by mutation

Molecular diagnosis: smMIPs and ddPCR. We identified pathogenic PIK3CA variants in 41 of 81 individuals (51%) via smMIPs (Figure 2A). No pathogenic variants were found in the 43 remaining genes in the smMIP panel (Supplemental Table 2). The average number of collapsed reads per MIP was 113, although the average reads per MIP varied widely across PIK3CA from 0 to 277 reads (Figure 2B). Variants were deemed diagnostic if found in 3 or more independent reads out of 30 or more total reads in a single sample. Samples with fewer than 30 average reads across PIK3CA were considered negative. All detected variants that met diagnostic criteria were previously reported as pathogenic in the literature and will be referred to as mutations for the remainder of this text. The 40 individuals (91 samples) with no variant detected via smMIP sequencing underwent screening for the most common recurrent GOF mutations via droplet digital polymerase chain reaction (ddPCR). Mutations were detected via ddPCR in 24 additional individuals (Figure 2, A and E).

We detected p.Glu545Ala substitutions (NM_006218, c.1634A>C) in 63 individuals (154 of 271 samples). This is a known artifact consistent with amplification of a noncoding pseudogene on chromosome 22q11.2 that is 95% homologous to PIK3CA coding exons 9–13 (22). All samples with detection of p.Glu545Lys and p.Glu542Lys via smMIPs were manually inspected in Integrated Genome Viewer (IGV) and the VAF was recalculated, excluding any p.Glu545Ala mismapped reads from the total number of reads. p.Glu545Ala reads contributed 0%–12% of total reads (0–45 reads) in individuals with p.Glu545Lys and p.Glu542Lys mutations.

In total, we detected HS mutations p.His1047Arg (c.3140A>G), p.Glu545Lys (c.1633A>G), and p.Glu542Lys (c.1624G>A) in 22, 18, and 18 individuals, respectively. Non-hotspot (NHS) mutations included p.His1047Leu (c.3140A>T) in 4 individuals, and p.Gly106Val (c.317G>T) and p.Gln546Lys (c.1636C>A) in 1 individual each (Table 1 and Figure 2E).

p.Gln546Lys (15 of 173 reads, 8.7% VAF) and p.Gly106Val (9 of 212 reads, 4.7% VAF) were confirmed via restriction fragment length polymorphism (RFLP) assays in the same tissue samples that screened positive on smMIPs (Figure 2, C and D). p.Gln546Lys is within the helical domain of p110α and has previously been reported in an isolated LM (12). p.Gly106Val has not previously been reported in an LM; however, it has previously been documented as a pathogenic GOF mutation in the ABD domain of p110α (23).

Correlation of number of samples with variant detection rate. We observed that unsolved individuals appeared to have fewer samples than those individuals with detected mutations. To test this observation, we employed a bootstrap analysis to assess variant detection sensitivity based on the number of lesion samples per individual. As the number of lesion samples per individual increased, the likelihood of detecting a pathogenic variant in at least 1 sample per individual also increased (Figure 3). The majority of unsolved individuals (94%) had only 1 or 2 lesion samples available for testing, and the interval increase in detection rate as the number of sampled lesions increases begins to decrease around 4 samples (Figure 3). Interestingly, orofacial LMs and microcystic LMs were more likely to be solved than neck/body or macrocystic LMs (P < 0.00001 and P = 0.013, Fisher, respectively), which is likely due to a lower number of lesion samples in the latter 2 groups (location: P = 0.0095, Wilcoxon; cystic structure: P = 0.018, Kruskal-Wallis).

Figure 3

Number of lesion samples correlates with likelihood of detecting a pathogenic variant. The results of a bootstrap analysis looking at likelihood of detecting a pathogenic variant given sampling of different number of lesion samples per individual. Mean detection rates for smMIPs only and smMIPs with ddPCR are shown, with error bars indicating the 2.5 and 97.5 percentiles. The number of individuals with adequate sample numbers for each analysis is indicated at the bottom of the chart. For smMIPs alone, the experimental rate of detection corresponded to the cumulative distribution function of the geometric distribution for P = 0.32, y = 1 – (0.68)^x, which demonstrates the likelihood of achieving a diagnostic result with a 32% mutation detection rate for 1 lesion sample. ddPCR, droplet digital polymerase chain reaction; smMIPs, single-molecule molecular inversion probes.

To further understand why increased variant detection was seen with increasing sample numbers, mean VAF and standard deviation were calculated for individuals with 3 or more lesion samples with detected variants. The average VAF across individuals with 3 or more lesion samples was 2.3%, with an average standard deviation of 1.6%. This variance between samples is exemplified by LR16-263, who had multiple lesion samples obtained at different surgeries and underwent RFLP confirmation of a p.Gln546Lys mutation (Figure 2C). One lesion sample from this individual had 432 reads but no detectable mutation, while 2 other lesion samples had 15 of 186 and 5 of 182 variant reads, respectively. This demonstrates the variability seen in a single lesion and supports the bootstrap analysis findings.

Analysis of phenotype by mutation and VAF. We analyzed genetic data by primary location, laterality, cystic structure (macrocystic versus microcystic, excluding mixed), number of invasive procedures, and volume. Given the small number of NHS mutations, p.His1047Leu, p.Gln546Lys, and p.Gly106Val were combined into a single group referred to as NHS mutations. Prior in vitro studies have demonstrated GOF effects for these NHS mutations (21, 23).

p.His1047Arg, p.Glu545Lys, p.Glu542Lys, and NHS mutations were all detected in orofacial and neck/body LMs and in both unilateral and bilateral lesions (Figure 4, A and B). NHS mutations were not detected in microcystic LMs. All 3 HS mutations and at least 1 NHS mutation were present in macrocystic and mixed LMs. The mutations detected in macrocystic and microcystic lesions were significantly different, with more p.His1047Arg mutations in microcystic lesions and more NHS mutations in macrocystic lesions (P = 0.015, Fisher; Figure 4C). A nonsignificant trend was noted when looking at the distribution of mutations by location, with orofacial LMs having a smaller proportion of NHS mutations than neck and body LMs (P = 0.091, Fisher; Figure 4A). No associations were detected between mutation and lesion laterality (Figure 4B), volume (Supplemental Figure 2), or number of invasive procedures (Figure 4D).

Figure 4

Mutation spectrum varies by cystic structure. Histograms depicting interactions between mutation and phenotype for (A) primary location, (B) laterality, (C) cystic structure, and (D) number of invasive procedures. Percentages are rounded to the nearest whole number. P values in the left upper corner indicate results of Fisher’s exact test when 2 categorical variables were present (location and laterality) or Kruskal-Wallis rank-sum test when 3 categorical variables were present (cystic structure and invasive procedures). *P < 0.05. In C, post hoc Fisher’s exact test showed a significant difference (P = 0.015) in mutations between macro- and microcystic lesions: **9%, ***4%, ****6%. NHS, non-hotspot mutations.

VAF ranged from 0.11% to 13.0%, with 98% of individuals having less than 10% VAF and 69% of individuals having less than 5% VAF. To assess for VAF patterns within specific mutations, we separated individuals by their PIK3CA mutation (combining NHS into 1 group) and compared VAF within each group to location, laterality, and cystic structure. Visual inspection of VAF plotted against location suggested a trend toward higher VAF with orofacial lesions and lower VAF with neck/body lesions within each mutation; however, the results of these analyses were not statistically significant, likely due to small sample size (Supplemental Figure 3, A and B).

We found no relationship between VAF and primary location, laterality, cystic structure, number of procedures (Figure 5B), or volume (Figure 5D) when VAFs for all mutations were combined. We noted a trend toward higher VAF for NHS mutations compared with HS mutations (P = 0.062, Wilcoxon’s rank-sum test; Figure 5A). While not statistically significant, this suggests a relationship between genotype and VAF, which would interfere with statistical analyses combining VAF data from multiple mutations.

Figure 5

Variant allele fraction (VAF) and genotype-adjusted VAF (GVAF) comparisons by phenotype. (A) Dot-box plot of VAF by mutation. Wilcoxon’s rank-sum test of hotspots versus NHS showed a trend toward increased VAF with NHS (P = 0.062). (B and C) Dot-box plots with (B) VAF (blue) on the left and (C) GVAF (green) on the right plotted by the 4 phenotype outputs: primary location, laterality, cystic structure, and number of procedures. Boxes are defined by the first quartile inferiorly, median, and third quartile superiorly, with whiskers extending to the farthest non-outlier point (defined as within 1.5 times the interquartile range) (A–C). Statistically significant results were seen for GVAF when analyzed by location and macro- versus microcystic structure. *P < 0.05 via Wilcoxon’s rank-sum test. (D and E) Dot plots with dots colored by cystic structure for volume by VAF (D) and GVAF (E) showed no correlation by Pearson’s correlation (n = 31). NHS, non-hotspot mutations.

GVAF. HS and NHS mutations upregulate p110α activity to different degrees in vitro (21, 23–29), suggesting that VAF may need adjustment prior to analysis to account for differences between mutations. No published in vitro experiments provide sufficient quantitative comparisons between HS and our detected NHS mutations. However, substantial data on PIK3CA mutation frequency are available in cancer databases. The COSMIC database (Catalogue of Somatic Mutations in Cancer, cancer.sanger.ac.uk/cosmic) contains data on 4 million mutations in cancer, including 13,969 mutations in PIK3CA (30). In the COSMIC database, the 3 HS mutations account for 67% of all PIK3CA variants (Table 2), while the 10 most common variants (including 5 mutations found in our cohort) account for 81% of all PIK3CA variants.

Table 2

Summary of the most common PIK3CA mutations in the COSMIC database and mutation strength adjustment factor

The relative frequency of individual PIK3CA mutations within the COSMIC database generally correlates with the level of pathway activation by in vitro studies (21, 23–27), but few studies have systematically analyzed across a large number of mutations. To address this, we analyzed previously reported quantitative data from growth assays for 27 PIK3CA pathogenic variants against the relative frequency of variants in the COSMIC database across all cancers (31). Pearson’s correlation revealed very strong correlations between COSMIC variant frequency and fold-change increase in growth in both an IL-3–restricted assay of Ba/F3 cells (a murine cell line whose growth is IL-3 dependent) and an insulin- and epidermal growth factor–restricted assay of MCF10A cells (a normal human breast epithelial cell line whose growth is insulin and epidermal growth factor dependent) (r = 0.72, P = 2.45 × 10^–5, and r = 0.82, P = 1.33 × 10^–7, respectively) (Supplemental Table 3). The observation of a correlation between COSMIC variant frequency and GOF strength in 2 distinct functional assays across 27 PIK3CA pathogenic variants supports the hypothesis that COSMIC mutation frequency could represent a surrogate for mutation GOF strength and be used to adjust VAF for the different GOF mutations observed in our study. We tested this hypothesis further by constructing genotype-specific adjustment factors based on COSMIC database mutation frequency (Table 2), to determine if the product of mutation adjustment factor and maximum VAF (a parameter we refer to as GVAF) correlated with measures of clinical severity.

When GVAF was compared to measures of clinical severity we found significant associations with location and cystic structure (Figure 5C). Orofacial location and microcystic structure were associated with higher GVAF, while neck/body and macrocystic lesions had significantly lower GVAF (P = 0.016 and P = 0.045, respectively, Wilcoxon’s rank-sum test). We also noted a trend toward increased GVAF with increased procedure numbers (P = 0.063, Kruskal-Wallis). An association was present between GVAF and de Serres stage (P = 0.029, Kruskal-Wallis) which was not seen when genotype or VAF alone was analyzed by stage (Supplemental Figure 4). No correlation was seen between GVAF and LM laterality or volume (Figure 5, C and E).

We report the largest cohort of individuals with isolated LMs who underwent PIK3CA sequencing to date. Taking advantage of 2 independent DNA sequencing methods (smMIPs and ddPCR) and a retrospective cohort of detailed clinical data, this study identifies factors correlating with the variable pathogenesis of LMs through a genotype/phenotype analysis for mosaic disorders.

Prior reports detected PIK3CA mutations in 16 of 17 (94%) individuals from a Boston cohort and 23 of 31 (74%) from a portion of this study’s cohort (11). Using a larger cohort and 2 deep targeted sequencing methods, we detected PIK3CA mutations in LM tissue in 64 of 81 (79%) individuals with isolated LMs, including 58 HS and 6 NHS mutations, with VAFs ranging from 0.1% to 13%. Furthermore, 44 of 64 (69%) individuals had a maximum VAF below 5%, which demonstrates the need for DNA sequencing methods optimized to detect very low VAFs within LM tissue–derived DNA.

Our data show a robust association between the number of samples tested and the mutation detection rate (Figure 3). Given the importance of accurate molecular diagnosis in the age of personalized therapeutics, these results strongly suggest that genetic testing should be performed on multiple samples. The samples comprising our data set either come from different locations within an individual LM or from different surgical procedures. With this in mind, we recommend that 3 to 4 separate lesion samples be obtained for DNA isolation and testing for LMs. Because approximately 10% of our PIK3CA-positive patients possessed NHS mutations, we recommend performing whole-gene sequencing to a depth of at least 300 reads. Alternatively, a tiered approach using HS allele–specific assays (such as ddPCR) first, reflexing to whole-gene sequencing if negative, may be more cost effective. This recommendation may also apply to other vascular anomalies and overgrowth disorders that occur only with mosaic mutations (32). Recently, 19 individuals with syndromic, nonneoplastic PIK3CA-related overgrowth disorders were reported to respond positively to treatment with the PI3K-specific inhibitor BYL719 (alpelisib), which has received FDA approval for treatment of breast cancer (33, 34). Identifying sensitive and efficient diagnostic approaches for individuals with LMs is critical given the interest in treating individuals with LMs with this class of drugs.

Correlations between PIK3CA mutations and disease severity have been previously noted. In a cohort of individuals with PIK3CA-associated neurodevelopmental disorder, HS mutations were found in subjects with hemimegalencephaly, while only NHS mutations were found in individuals with megalencephaly-capillary malformation syndrome (MCAP), a less severe phenotype (35). Another study of 162 individuals with heterogeneous brain and/or body overgrowth phenotypes found more strong oncogenic mutations (HS mutations and other strongly activating mutations) in individuals without brain overgrowth, while those with brain overgrowth were more likely to have intermediate or weak activating mutations (36). Consistent with prior studies, we found a correlation between GOF strength of PIK3CA mutation and the cystic structure of the LM (Figure 4). Our report is the first to our knowledge to demonstrate a genotype-phenotype correlation in a cohort that consists entirely of individuals with LMs, the most common subtype of PIK3CA-related overgrowth disorders (37).

Detection of genotype-phenotype associations in diseases associated with postzygotic mutations is challenging, due in part to the additional variables (timing of mutation in development, affected cell type, and VAF) that are not present in nonmosaic disorders. The studies described above (35, 36) did not find a difference in affected tissue VAF between disease groups, but Kuentz et al. found increased VAF in affected tissue compared with blood and increased VAF in blood in patients with brain overgrowth compared with those with only body overgrowth. A recent study on arteriovenous malformations (AVMs) with mosaic KRAS and BRAF mutations found a negative correlation between VAF and AVM nidus size, with larger AVMs associated with lower VAFs (38). These studies did not, however, provide analysis of affected gene or individual mutations and their relationship to VAF.

Our large cohort of individuals with a single disease allowed for genotype-phenotype analysis not previously possible for disorders associated with mosaic mutations. We detected a statistically significant association between HS/NHS mutations and cystic structure, with more frequent NHS mutations in macrocystic lesions compared with microcystic lesions (Figure 4C). We also detected a trend toward lower VAF with HS mutations compared with NHS mutations (even with only 6 NHS mutations) and noted a trend toward increased proportion of NHS mutations in neck/body lesions. Our results suggest that combining genotype and VAF together may capture associations that would not be detected by analyzing genotype and VAF separately. This led us to develop the GVAF metric.

A large number of functional studies have already demonstrated that different PIK3CA mutations lead to activation of downstream pathways to different degrees (21, 23–29), but most studies have looked at only a handful of PIK3CA mutations, or have used qualitative assays to rank mutation strength. In the most systematic and comprehensive study of the effect of oncogenic PIK3CA mutations to date, Dogruluk et al. modeled 27 PIK3CA mutations in growth factor–free cell survival assays in 2 distinct cell types, and found strong correlations between the magnitude of cell survival and PIK3CA mutation frequency in breast cancer (31). We also detected strong correlations between magnitude of cell survival and COSMIC frequency of PIK3CA mutations across all cancers (Supplemental Table 3). These findings support the use of COSMIC frequency of PIK3CA mutation as a surrogate for PIK3CA GOF strength. A limitation of this approach is the fact that it assumes the effects of activating PIK3CA mutations in cell lines and cancers will be similar to that in lymphatic endothelial cells — the cell type most likely to be driving LMs. However, the presence of correlations between clinically important variables in lymphatic malformations (cystic structure and anatomic location, Figure 5) and our GVAF metric support the utility of this approach. Future studies modeling these specific PIK3CA mutations in lymphatic endothelial cell culture, or in animal models of LM, will be needed to provide further support for, and perhaps modification of, the GVAF metric.

We identified statistically significant associations between GVAF and 2 important indicators of clinical severity: location and cystic structure. Higher GVAFs were seen in orofacial LMs compared with neck/body, and higher GVAFs were associated with microcystic lesions when compared with macrocystic LMs (Figure 5C). We also found statistically significant differences in GVAF by de Serres stage (Supplemental Figure 4) and a trend toward higher GVAF with increasing number of procedures (Figure 5C). These measures of severity are not independent, as demonstrated by our initial analysis of the cohort; however, this remains an important development in the understanding of LM severity in relation to genotype and provides the first demonstration to our knowledge of combined mutation and VAF data analysis compared to clinical severity.

The first LM clinical classification divided lesions into 2 categories: macrocystic infrahyoid lesions (type I LMs) and microcystic suprahyoid lesions (type II LMs) (5). Decades of clinical experience have confirmed the initial observation that type II LMs persist and are more difficult to treat than type I LMs. Type I LMs respond well to sclerotherapy and are more likely to resolve spontaneously, while type II LMs have a poor response to sclerotherapy and greater impact on breathing and swallowing (2, 4, 39–41). Consistent with these observations, this study identified higher GVAFs in the more serious orofacial and microcystic lesions, providing a genetic explanation for disease severity.

The 2-tiered description does not encompass the full range of LM phenotypes defined in more recent classification systems (6), but provides insight into what may be 2 distinct biological processes. In a recent study, macrocystic and microcystic LMs were found to have different expression profiles of genes in the PI3K/AKT/mTOR pathway (42). Oncogenic transcripts including RAN and CDC20 were overexpressed in microcystic lesions compared with controls, while macrocystic lesions had overexpression of genes related to cellular respiration, epithelial-mesenchymal transition, and adaptation to hypoxia. These observations further support the hypothesis that cystic structure is a by-product of upstream developmental differences between lesions, as opposed to a purely anatomic difference, which is also seen in our data.

Our study has several limitations that are common when studying mosaic disorders. Since genetic analysis in mosaic disorders required affected tissue, we only included subjects that underwent surgical intervention, biasing our cohort towards more severe phenotypes. We also cannot determine how spatial or temporal differences might contribute to variation in VAF because samples were obtained from multiple locations within a lesion or from different surgeries. While our cohort is the largest published for a single vascular malformation, it remains relatively small, with a skewed distribution of phenotypes. Although LMs have a propensity to occur in the head and neck, the majority of our subjects were ascertained through the head and neck surgery team, resulting in more orofacial than neck or body LMs. Fewer samples were available per individual for neck and body LMs compared with orofacial lesions, which led to a lower detection rate in those groups. This is likely a consequence of neck and body lesions having fewer surgical interventions due to decreased clinical severity. Nevertheless, we were able to detect a strong association between lower GVAF and neck/body and macrocystic lesions despite the unequal distribution of phenotypes. Although we provide evidence that GOF strength is correlated with rank in COSMIC, we recognize that the genotype adjustment factor is not a perfect surrogate for GOF strength. Nonetheless, the validity of this approach is supported by the GVAF correlation with clinical variables. Testing these variants in lymphatic endothelial cell culture or some other model of lymphatic malformations will be needed to further address the GVAF hypothesis.

As sequencing technologies have improved, mosaic diseases have become a more prevalent topic. However, the study of these diseases remains challenging due to the need for affected tissue samples to uncover a molecular diagnosis. Within the spectrum of congenital PIK3CA-related overgrowth disorders, LMs are the most common phenotype. They are associated with a narrow spectrum of PIK3CA mutations and recent studies suggest a single cell type, endothelial cells, may be mutated in LMs (12, 14). Thus, LMs may serve as a model phenotype in which we can develop tools to better understand mosaic disorders in general. Traditional genotype-phenotype analyses have limited utility in understanding mosaic disorders due to the complexity of multidimensional genotype factors. We propose that the “mosaic genotype” encompasses more than just the mutation or the VAF, but that these 2 previously distinct data points should be combined with adjustment for the mutation’s effect on cell function. In this study, using the GVAF metric that integrates independent data on mutation severity and VAF, we demonstrate a consistent correlation between mosaic genotype and indicators of LM clinical severity, including location and cystic structure.

Study design. This study was a retrospective analysis of clinical data from individuals with LMs combined with deep targeted sequencing of banked tissue. The goal of the project was to characterize the PIK3CA mutations contributing to LMs and correlate the genetic data to clinical and radiographic attributes from the same LM. We included all individuals with isolated LMs treated at Seattle Children’s Hospital between 2000 and 2015 who had LM tissue available for analysis. Individuals with accompanying overgrowth syndromes such as fibro-adipose vascular anomaly (FAVA), Klippel-Trenaunay syndrome (KTS), and congenital lipomatosis, overgrowth, vascular malformations, epidermal nevi, and skeletal/spinal anomalies (CLOVES) were excluded.

Clinical features. LM features were obtained from imaging review or clinical reports of imaging, physical exam, and operative reports if images were not available. LM features included primary location, laterality, and cystic structure. Primary location groups were based on prior clinical studies that separated head and neck LMs by their relation to the hyoid bone (6). Lesions isolated to the neck (no suprahyoid involvement) were considered neck LMs and any involvement of suprahyoid structures was grouped into a single category of orofacial LMs based on established greater severity of suprahyoid lesions (5, 6). Any lesions outside of the head and neck were categorized as body LMs. Individuals were divided into 3 cystic structure groups: macrocystic, microcystic, and mixed. Macrocystic was defined as lesions with more than 90% large cysts (≥2 cm cyst diameter), microcystic as 100% small cysts (<2 cm), and mixed as having both large and small cysts. Each head and neck LM was staged according to the de Serres clinical staging system (6). Other clinical data were collected retrospectively, including age at diagnosis and intervention, family history of vascular anomaly, and number of invasive procedures. Procedures related to the LM and requiring an operating room were counted as invasive procedures, including sclerotherapy, laser treatments, resections, direct laryngoscopy for airway evaluation, tracheostomy, etc. Numbers of procedures were grouped into 0–1, 2–5, and ≥6 procedures.

Volumetric analysis. Subjects with preoperative MRIs performed at Seattle Children’s Hospital underwent volumetric analysis of LMs to determine total lesion volume. Volume calculations were performed by a board-certified radiologist using hand-drawn regions of interest and OsiriX (Pixmeo SARL) to calculate volume by multiplying slice thickness by the outlined area and adding up individual slice volumes (43, 44).

Sample acquisition and DNA extraction. After obtaining IRB approval and consent, LM tissue was prospectively collected at clinically indicated surgical procedures, flash frozen, and stored in a biorepository at –80°C. DNA isolation from blood and serum was performed with a Puregene Blood Core Kit (QIAGEN). DNA from lesions, mucosa, and skin was isolated with a PureLink Genomic DNA Mini Kit (Invitrogen) or an Oragene Saliva Kit (DNA Genotek).

PIK3CA targeted sequencing using smMIPs. A pool of 1707 smMIPs was designed targeting the coding sequences of 44 genes using MIPgen (http://shendurelab.github.io/MIPGEN/), as previously reported (45, 46). Forty-two smMIPs in this pool targeted 100% of the coding sequences in PIK3CA. We used 200 ng of genomic DNA in each capture reaction, with a 5000:1 ratio of phosphorylated smMIPs to DNA. Massively parallel sequencing was performed with the Illumina HiSeq to obtain 151-bp paired-end reads (Illumina).

Variant analysis pipeline. Reads were aligned to human assembly hg19 using BWA (47). Reads with the same molecular barcode were collapsed to form independent reads. Variant analysis was performed by 2 pipelines: a standard pipeline comprising consensus calls between GATK and FreeBayes (45) and a mosaic variant caller, PISCES (48). Samples negative by either pipeline were visually inspected with IGV (Broad Institute) at the 3 HS residues in PIK3CA (49). All variants were reported on NM_006218. We required variants to be present in 3 or more independent reads with coverage of at least 30 independent reads. Variants were only considered if they met standard guidelines (50).

RFLP. NHS variants other than p.His1047Leu were confirmed with RFLP assays. Primers for the 521-bp product containing the p.Gln546Lys (c.1636C>A) mutation were forward, 5′-TGGTTCTTTCCTGTCTCTGAAAA-3′, and reverse, 5′-GCATTTAATGTGCCAACTACCAATG- 3′. The c.1636C>A mutation created an HpyAV digestion site resulting in 2 bands at 267 and 254 bp. Primers for the 526-bp product containing the p.Gly106Val (c.317G>T) mutation were forward, 5′-ATGCCTCCACGACCATCATC-3′, and reverse, 5′-GGTGTTAAAAATAGTTCCATAGTTCGA-3′. The c.317G>T mutation created a HincII digestion site resulting in 2 bands at 208 and 318 bp. All reactions were performed using HotStarTaq (QIAGEN) and 12 ng of DNA per reaction. The PCR reactions were performed using the following cycling conditions: 95°C for 5 minutes, 94°C for 30 seconds, 52°C for 30 seconds, and 72°C for 30 seconds, 35 cycles, followed by a final extension at 72°C for 5 minutes. Following the PCR, 10 μL of PCR product was digested with 10 units of restriction enzyme in 30 μL at 37°C overnight. DNA fragments were analyzed by 3% agarose gel electrophoresis.

Targeted mutation detection with ddPCR. Bio-Rad–designed ddPCR assays were used to target 5 PIK3CA mutations (p.Cys420Arg, p.Glu542Lys, p.Glu545Lys, p.His1047Arg, and p.His1047Lys) in samples undiagnosed by smMIPs. Each 24-μL PCR mixture contained 10 μL of ddPCR SuperMix, 1 μL of assay, 8 ng of DNA, and 12 μL of water. Twenty microliters of the PCR mixture was used to generate droplets using a Bio-Rad QX200 Droplet Generator. PCR was performed with the following parameters: 10 minutes at 95°C, followed by 40 cycles of 30 seconds at 94°C, 60 seconds at 55°C, and 10 minutes at 98°C, and held at 4°C. We used the Bio-Rad QX200 Droplet Reader within 24 hours and analyzed the data with Bio-Rad QuantaSoft software. Each run included a no-template control, WT control, and a positive mutant control. Samples were positive for the variant if (i) there were 4000 or more WT-positive droplets, and (ii) if the mutant fluorescence was significantly different from the WT control, given the calculated 95% confidence intervals. VAFs were calculated as the concentration of mutant droplets out of the total concentration of droplets containing at least 1 copy of DNA.

GVAF. To adjust for the strength of activation associated with each mutation, an adjustment factor was calculated using the relative presence of the mutation in the COSMIC database, accessed on December 11, 2018 (cancer.sanger.ac.uk) (30). We compiled the number of mutations recorded in COSMIC for each mutation seen in our cohort. Data on cell growth from 2 growth factor–dependent cell line assays, Ba/F3 cells (a murine cell line that dies in the absence of IL-3) and MCF10A cells (a human breast epithelial cell line that dies in the absence of insulin and epidermal growth factor), and correlation analysis (Pearson’s) between measures of cell growth (fold-change increase relative to vector-only cells) and relative COSMIC frequency was performed (Supplemental Table 3) (31). The frequency of each mutation in comparison to the total number of PIK3CA mutations was calculated and rounded to the nearest whole number to become our adjustment factor. Any mutation seen at less than 1% frequency was assigned an adjustment factor of 1. This frequency, or the genotype adjustment factor, was multiplied by the individual’s maximum VAF to calculate the GVAF (Table 2).

Statistics. Data were analyzed to detect differences between groups using Fisher’s exact test for discrete and ordinal data, Wilcoxon’s rank-sum test for 2 group comparisons of continuous variables, or the Kruskal-Wallis test for >2 group comparisons of continuous variables. Nonparametric tests were used due to small sample size and non-normal distributions. A P value less than 0.05 was considered significant. All analysis was performed using R software and graphs generated with ggplot2.

Bootstrap analysis of mutation detection likelihood. We compiled the number of lesion samples and corresponding mutation results for each individual and performed a bootstrap analysis using R. The first iteration calculated the likelihood of mutation detection for our cohort given 1 lesion sample per individual. The program randomly selected a lesion sample from each subject and calculated the diagnostic rate for the cohort. This process was repeated 10,000 times and the average diagnostic rate, 2.5 percentile, and 97.5 percentiles were recorded. The bootstrap analysis was then performed for all subjects with 2 or more lesion samples and 2 lesion samples were randomly selected. The individual was considered to have a mutation detected if one or more samples had a mutation and the diagnostic rate for the population was again calculated using 10,000 iterations. This process was repeated through 6 samples per subject. The program was then run using only the data obtained from smMIPs and also only data obtained from ddPCR.

Study approval. This study was approved by the Institutional Review Board at Seattle Children’s Hospital. Written, informed consent, including use of patient photos, was obtained for each individual involved in this study prior to tissue collection.