The mutational landscape of adenoid cystic carcinoma.

Adenoid cystic carcinomas (ACCs) are among the most enigmatic of human malignancies. These aggressive salivary gland cancers frequently recur and metastasize despite definitive treatment, with no known effective chemotherapy regimen. Here we determined the ACC mutational landscape and report the exome or whole-genome sequences of 60 ACC tumor-normal pairs. These analyses identified a low exonic somatic mutation rate (0.31 non-silent events per megabase) and wide mutational diversity. Notably, we found mutations in genes encoding chromatin-state regulators, such as SMARCA2, CREBBP and KDM6A, suggesting that there is aberrant epigenetic regulation in ACC oncogenesis. Mutations in genes central to the DNA damage response and protein kinase A signaling also implicate these processes. We observed MYB-NFIB translocations and somatic mutations in MYB-associated genes, solidifying the role of these aberrations as critical events in ACC. Lastly, we identified recurrent mutations in the FGF-IGF-PI3K pathway (30% of tumors) that might represent new avenues for therapy. Collectively, our observations establish a molecular foundation for understanding and exploring new treatments for ACC.

Adenoid cystic carcinomas (ACCs) are malignancies that cause significant morbidity and mortality. Typically arising from salivary glands, they are characterized by unpredictable growth, extensive perineural invasion, and high rates of metastasis, ultimately resulting in low survival rates. Treatment remains limited to surgery and radiation, and no systemic agent has been found to be effective.[1] Although high EGFR and c-kit expression have been identified[2,3], drugs against these targets have not led to substantive clinical responses.[4] Other initial findings include the identification of MYB-NFIB fusion proteins but their significance remains unclear.[5] A deeper grasp of the underlying genomic foundations of ACC is needed to understand its molecular basis and guide the development of effective therapies.

We sequenced the entire exome (n=55) or genome (n=5) of 60 ACC samples with matched normal DNA (Table 1). This approach generated 2,221 gigabases of mapped sequence with 92.4% of the target sequence covered to at least 10x depth. Sequencing generated a mean exome and genome coverage of 106x and 37x, respectively. To ensure the accuracy of our massively parallel sequencing data, we conducted extensive validation of nearly every candidate somatic mutation identified (2,751 variant calls) (Supplementary Fig. 1) using targeted re-sequencing (Supplementary Fig. 2, Supplementary Table 1). In addition, we also performed FISH analysis for the MYB-NFIB translocation (Supplementary Fig. 3, Supplementary Table 2).

Table 1

ACC whole exome and whole genome sequencing statistics.

Statistic	Whole exome	Whole Genome
Tumor/normal pairs sequenced	55	5
Total Gb sequenced	971.5	1249.4
Mean fold tumor target coverage	100 (23–167)	38 (34–45)
Mean fold normal target coverage	112 (25–185)	35 (33–37)
Mean somatic mutation rate per Mb	0.4 (0.06–0.8)	0.3 (0.2–0.5)
Mean number of validated nonsynonymous mutations/patient	11 (1–36)	11 (5–17)
Mean number of transcribed noncoding mutations/patient (UTRs)	6 (1–13)	3 (1–6)
Total number of structural rearrangements	N/A	17
Total number of frame-preserving genic rearrangements	N/A	1

We identified a mean of 22 somatic mutations per sample, corresponding to approximately 0.31 non-silent mutations per MB. This rate is quite low compared to most adult solid tumors such as head and neck squamous cell carcinoma[6,7] and colon cancer[8] yet similar to hematologic malignancies and neuroblastoma.[9-11] The transition/transversion ratio (Ti/Tv) was 1.1, similar to some carcinogen-driven malignancies[6,7,12] but unlike most described cancer types.[13] The somatic mutation frequency correlated with solid histology (Wilcoxon rank-sum test, p = 4.0 × 10−2), and MYB translocations occurred in 57% of samples (34/60).

We validated 710 distinct nonsynonymous mutations across 643 genes (1–36 per tumor) (Fig. 1a, Supplementary Table 3, Supplementary Fig. 4). This represents substantial mutational heterogeneity across tumors (Fig. 1b, Fig. 2). We employed CHASM, a widely used approach for distinguishing driver from passenger mutations,[14] to identify multiple potential driver mutations, including those in PIK3CA, TP53, PTEN, SMARCA2, KDM6A, and CREBBP (Supplementary Tables 4, 5). Analysis of these driver genes demonstrated marked enrichment in pathways involved in chromatin remodeling, DNA damage, MYB, protein kinase A (PKA) signaling, and PI3K signaling (Fig. 1c–e). For example, despite low overall mutation frequency, 35% of ACC tumors were mutated in chromatin regulators. These data suggest that the ACC mutational landscape is characterized by combinations of MYB pathway alterations and mutations in specific biological processes (Fig. 2). Interestingly, a small subset (n=8) were observed with no CHASM-designated driver mutations. It is possible that some mutations in these tumors are drivers not called by CHASM or that other, non-exonic alterations are important in these tumors.

Figure 1

Mutational landscape of adenoid cystic carcinoma

(a). Number of validated nonsynonymous somatic mutations per sample across 60 ACC cases. (b). Representative list of recurrent nonsynonymous somatic mutations. Multiple mutations within a given sample in the same gene were only counted once. (c). Pathways affected by driver mutations identified by CHASM. Bonferroni FDR-corrected p-values for pathway enrichment are shown. (d). Location of ACC mutations in key chromatin remodeling genes. HSA, helicase-SANT-associated; Zf-TAZ, TAZ zinc finger; DUF, domain of unknown function; KAT11, histone acetylation protein; ZZ, ZZ-type zinc finger; SNF2_N, SNF2 family N-terminal domain; Helicase_C, Helicase conserved C-terminal domain; TPR, tetratricopeptide repeat; JmjC, JmjC hydroxylase domain; FtsJ, FtsJ-like methyltransferase. (e). Location of ACC mutations in key established cancer genes. PI3K_p85B, PI3K p85-binding domain; PI3K_rbd, PI3K ras-binding domain; PI3Ka, PI3K accessory domain; EGF_Ca, calcium-binding EGF domain; hEGF, human growth factor-like EGF; Ank repeat, ankyrin repeat; DSPc, dual specificity phosphatase, catalytic domain; P53 TAD, P53 transactivation motif; P53_tetramer, P53 tetramerization motif.

Figure 2

Integrated analysis of adenoid cystic carcinoma genetic alterations

Clinical data, validated somatic mutations, validated structural variants, and copy number alterations for 60 ACC cases. Genes highlighted in bold contain significant driver mutations identified by CHASM.[14] Remaining genes contain selected mutations of interest that did not reach significance by CHASM but which are either altered in other malignancies or are integral components of affected pathways. Orange, nonsynonymous somatic mutation (missense, nonsense, splice site, indel); blue, homozygous copy number deletion; red, high-level copy number amplification.

We used exome and genome sequencing data to characterize the copy number landscape of these tumors. We analyzed somatic copy number variations (CNVs) using ExomeCNV[15] and found high concordance with a subset (n=12) that underwent array-based analysis. GISTIC2.0[16] identified recurrent high-level losses in 6q24, 12q13, and 14q (Fig. 3a, Supplementary Table 6).[5] Samples with 14q loss were more likely to be of solid histology (Fisher’s exact test, p = 2.0 × 10−4), while samples with 6q24 loss were enriched for advanced stage (p = 4.0 × 10−2). Expression array analysis on 23 ACC tumors found no distinct subgroups (Supplementary Fig. 5). Genes harboring driver mutations were confirmed to be generally expressed in ACC tumors (Supplementary Table 7).

Figure 3

Structural variations and copy number landscape of adenoid cystic carcinoma

(a). Clustered copy number profiles of 60 ACC cases. Clustering revealed 4 sub-clusters of tumors with one each defined by loss of 14q, loss of 12q13, a cluster with two samples containing a large number of copy number alterations, and most tumors in a large group with minimal copy number changes. White, normal (diploid) copy number log-ratio; blue, copy number loss; red, copy number amplification. (b). Circos plots[73] of genetic alterations in 5 ACC cases. Plots depict validated structural genetic variants, DNA copy number alterations, intra- and inter-chromosomal translocations, and sequence alterations. Loss-of-heterozygosity, green track; amplification, red track; copy number loss, blue track; missense mutation, black gene name; nonsense mutation, blue gene name; splice site mutation, orange gene name; inter-chromosomal translocation, red line; deletion, gray line; intra-chromosomal translocation, blue line; insertion, purple line.

Whole genome paired-end sequencing identified numerous structural variants (SVs) (Fig. 3b, Supplementary Table 8), with the existence of 17 SVs across 5 samples confirmed using PCR (Supplementary Fig. 6, Supplementary Table 9). Combined with FISH data, MYB translocations were the only recurrent SVs detected. However, we cannot exclude the presence of less common recurrent translocations.[17] Intriguingly, one sample harbored a tandem duplication within FGFR2. Similar alterations comprised of dimerized, constitutively active variants have been described in hematologic malignancies.[18] We did not detect a fusion transcript from this tumor using RT-PCR, but this does not rule out an alternative transcript configuration. Larger numbers will need to be analyzed to further characterize its incidence and impact. The remaining variants were not predicted to be in-frame. Collectively, few CNVs or SVs were found across the cohort, demonstrating the relatively “quiet” nature of the ACC genome. Thus, our data indicate that ACC’s major SVs are MYB translocations and recurrent deletions on 6q24, 12q13, and 14q.

A prominent feature of the ACC mutational landscape is the presence of multiple mutations targeting chromatin remodeling genes (35% incidence) (Fig. 1d, Fig. 4a). Such alterations are increasingly recognized as playing key roles in oncogenesis[19,20] and have been reported in various other tumors[10,21] but not ACC. Among ACC alterations, chromatin state modifiers were substantially enriched for somatic mutations (q = 4.5 × 10−3). We identified multiple aberrations in the SWI/SNF-related, matrix associated, actin dependent regulator of chromatin (SMARC) family, including SMARCA2 (5%) and single mutations in SMARCE1 (2%), ARID1A (2%), and ATRX (2%). SMARC mutations contribute to the development of both cancers and genetic diseases.[21-25] SMARCA2 encodes a core catalytic subunit of the SWI/SNF complex involved in regulating gene transcription.[26] All SMARCA2 mutations were clustered within the Helicase C family domain (T1126I, G1132V, G1164W). Notably, mutated helicases have been shown to increase cancer susceptibility, likely by disturbing core repair mechanisms.[27] Similarly, we identified likely loss of function in SMARCE1 (2%), another SWI/SNF complex gene with a HMG box domain mutation (Y73C) that likely affects DNA binding.[28] Additional recurrent mutations include CREBBP (7%), a histone acetyltransferase that enables transcription across multiple critical signaling pathways.[24,25,29] The majority of CREBBP mutations clustered in the critical KAT11 histone acetylation domain (R1446C, I1453N, W1472S).[30] We also observed numerous KDM6A missense mutations (7%), a histone demethylase also termed UTX.[31] The functional impact of these KDM6A mutations was evaluated via a well-characterized H3K27me3 assay. We observed abrogation of demethylase activity in cells overexpressing mutant but not those overexpressing wildtype KDM6A (Fig. 5a–b). Moreover, whereas wildtype KDM6A suppressed growth, mutants either lost the ability to suppress growth or in some cases, augmented it (dominant phenotype).

Figure 4

Diversity and frequency of ACC genetic changes culminating in aberrant signaling pathways

Alterations are defined by somatic mutations, homozygous deletions, high-level focal amplifications, and structural variants validated by FISH or PCR. Frequencies are expressed as a percentage of all cases. Red background denotes activating alteration, blue background denotes inactivating alteration, white background denotes numerous alterations, and purple background denotes alteration of unclear significance. Interactions based upon IPA or Pathway Commons.[38] (a). Epigenetic modification. (b). DNA damage/checkpoint signaling pathway. (c). MYB/MYC signaling pathway, with fusion denoted by recurrent t(6;9) translocation between MYB and NFIB. (d). FGF/IGF/PI3K signaling pathway. (e). Notch signaling pathway.

Figure 5

Functional consequences of genetic alterations in KDM6A, TP53, PI3K, and FAT4 pathway genes

(a). Overexpression of KDM6A mutants identified in ACC compared to overexpressed wildtype demonstrates increased growth. Experiments performed in triplicate. (b). KDM6A tumor-specific mutants exhibit abrogated H3K27me3 demethylase activity relative to wildtype KDM6A. Representative immunostaining of cells overexpressing wildtype or mutant KDM6A is shown with corresponding quantitation. Arrows represent cells expressing FLAG-tagged KDM6A construct. Experiments performed in triplicate. Scale bars, 20µm. (c). Heatmap shows differential expression of key TP53 transcriptional targets in TP53-pathway altered ACC samples (TP53-mutant or MDM2-amplified) compared to TP53-wildtype ACC samples and normal salivary tissue. (d). PI3K-mutant ACCs demonstrate marked downstream p-AKT and p-PRAS40 immunohistochemistry staining relative to PI3K-wildtype ACCs. Scale bars, 120µm. (e). FAT4 knockdown demonstrates increased growth in human foreskin fibroblasts (HFF-1), immortalized human salivary cells (HSG), and human salivary adenocarcinoma cells (HSY). Experiments performed in triplicate. KD, knockdown. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001, ns, not significant (ANOVA). Error bars, ± 1 S.E.M.

Other related mutations include those involved in histone acetyltransferase/deacetylase activity (EP300, ARID4B, ARID5B, BRD1), histone methyltransferase/demethylase function (FTSDJ1, MLL3), as well as histones themselves (HIST1H2AL, HIST1H1E). EP300 binds CREBBP as a transcriptional co-activator that regulates cell proliferation and differentiation[32], and harbored a splice site mutation within its KAT11 histone acetylation domain (G1429_splice). HIST1H2AL encodes a core component of the nucleosome, while HIST1H1E is known to link histones for compaction into higher-order structures. Importantly, histone mutations have been implicated in cancers such as pediatric gliomas.[8,33,34] While mutations in our study were observed within known domains, (R18L in core histone H2A/H2B/H3/H4 and K75N in linker histone H1/H5 domains, respectively), further investigation will be necessary to determine their functional impact. Overall, a diversity of chromatin genes was involved, suggesting that disparate means of chromatin dysregulation may promote ACC development through convergent pathways. Given the emerging importance of epigenetic agents on other disorders[35,36], similar treatments for ACC cases with chromatin derangements may hold promise.

We observed a second group of mutations in genes involved in DNA damage response (q = 5.6 × 10−3) (Fig. 4b). In TP53 (5% incidence), one missense (P151S) and 2 nonsense mutations (R213*, R342*) in the central binding and tetramer motifs were found (Fig. 1e). Similarly UHRF1, a ubiquitin ligase involved in the p53-dependent DNA damage checkpoint[37], was altered in 8% of patients (2 mutations, 3 homozygous deletions). TXNIP, a protein frequently repressed in cancers[38,39], was noted to have a frameshift insertion in its arrestin domain (L129fs). We also discovered 2 missense mutations each within individual samples in ATM, BRCA1, and DCLRE1A, a DNA cross-link repair gene. Notably, markedly decreased expression in key TP53 transcriptional targets was specifically identified in tumors with TP53 pathway alterations (binomial test, p = 1.0 × 10−4) (Fig. 5c).

Previous work has implicated MYB as an ACC fusion oncogene with NFIB via a t(6;9) translocation[40] (Fig. 4c). Indeed, we verified 57% of our cohort as MYB-NFIB(+) by FISH analysis. We extend this initial finding and observe further MYB pathway dysregulation with alterations found in an additional 8% of cases. Specifically, splice site and coding mutations involving exon 10 of MYB were seen, potentially disrupting its leucine-rich negative regulatory domain[40] and triggering constitutive activation. Five additional homozygous deletions were identified on CNV analysis (4/5 were MYB-NFIB(+)), though these may be secondary alterations around the translocation. Other mutations identified in the MYB pathway included MYCBP2, MGA, and MCM4. In NFIB, we observed two truncating mutations in the CTF/NFI transcription modulation domains (Y249*, P390fs), in addition to 4 homozygous deletions (all were MYB-NFIB(+)). Our data highlight MYB in ACC as an active oncogenic partner in fusion transcripts but also suggest a separate role for NFIB given the presence of specific mutations in this gene.

Multiple genes involved in the PKA pathway were found to be mutated (27% incidence, q = 4.2 × 10−3). We identified recurrent mutations in RYR3 (7%) and RYR2 (2%), both tetrameric, intracellular calcium channels implicated in breast cancer progression.[41,42] Mutations were also observed in PKA pathway-associated genes that are known or potential tumor suppressors, including tyrosine phosphatases (PTPRG, PTPRH, PTPRJ, PTPRK). Nonsense (PTPRG, PTPRH) and frameshift mutations (PTPRK) abrogate their phosphatase domains (E736*, W602*, L457fs, respectively), while an additional 4 samples had PTPRK homozygous deletions. PKA-related genes are known to facilitate CREBBP recruitment via calcium influx[43] and CREB phosphorylation, regulating multiple critical developmental processes involved in tumorigenesis.[24,44]

Mutations in genes important in the FGF/IGF/PI3K pathway were also identified in 30% of our cohort (q = 2.4 × 10−2) (Fig. 4d). Three tumors harbored separate missense mutations in PIK3CA (Fig. 1e). Another tumor had two mutations in the catalytic domain of PTEN (R130fs, K144Q). All five lesions occurred in hotspots (COSMIC v61).[45] FOXO3[46] was also altered in 7% of patients. We further identified mutations in receptor kinase/PI3K-associated genes, including FGF16, FGFR4, IGFBP2, and ILR17RD. Fibroblast growth factors and insulin-like growth factors are potent oncogenic PI3K activators[47,48], while IL17RD appears to be a PI3K inhibitor through its interaction with FGFR.[49] The samples with PI3K-associated insults all demonstrated solid histology, the most aggressive ACC variant (Fisher’s exact test, p < 1.6 × 10−3). We further evaluated PI3K-altered samples using immunohistochemistry (p-AKT and p-PRAS40). Functional activation of the PI3K pathway was observed in all ACC tumors harboring PIK3CA or PTEN mutations, but not in wild-type tumors (Fig. 5d, Supplementary Fig. 7). This is supported by the presence of significantly enriched PI3K signatures in PI3K-mutant ACC tumors via GSEA (Gene Set Enrichment Analysis) (p < 1.0 × 10−3) (Supplementary Fig. 8a). Together, these findings delineate a previously undescribed subset of ACC patients which we hypothesize may benefit from agents targeting this pathway.[50]

Despite not meeting significance by CHASM analysis, we observed alterations in the Notch signaling pathway in 13% of samples (Fig 4e). In addition to recurrent NOTCH1 (5% incidence, including 3 missense and one nonsense) and FOXP2 (3%) mutations, predicted functional mutations were found in DTX4, FBXW7, and CNTN6. FBXW7 (R465H) is a tumor suppressor, targeting c-myc and NOTCH1 for degradation.[51] Depending on cellular context, NOTCH1 has a dual role in cancer as an oncogene or tumor suppressor.[6,7,52-55] Analysis of NOTCH1-altered ACC samples using GSEA[56] showed a trend toward Notch signaling enrichment (Supplementary Fig. 8b). These alterations require further characterization, but may have therapeutic implications.[57]

Intriguingly, we observed mutations in various genes involved in cell metabolism processes. Of particular interest are recurrent alterations in HSPG2 or perlecan (7%), a basement membrane-based proteoglycan known to modulate tumorigenic growth factors, including FGF. FGF signaling promotes mitogenesis and angiogenesis[58], while deletion of perlecan heparan sulfate domains has been shown to impair tumorigenesis.[59-61] Importantly, ACC produces high perlecan levels within its pseudocystic cribiform structures[62,63], and maximum perlecan synthesis occurs during the tumor’s proliferative phase.[64] Additionally, we identified a hotspot IDH1 mutation in its catalytic domain (R132H).[33,65,66,67]

ACC is unique in its propensity for perineural invasion and distant metastasis. We noted numerous mutations in cell adhesion proteins (29% incidence). We found mutations in the active domains of genes responsible for neuronal axon guidance cues, including NTNG1, SEMA3G, and SEMA5A. Netrins and semaphorins are implicated in neuronal invasion[68] and metastasis.[69] Additionally, mutations in protocadherins were observed, including truncating mutations in FAT3 and FAT4. FAT4 has been hypothesized as a tumor suppressor in gastric adenocarcinoma.[70] We examined FAT4 through RNAi knockdown in human salivary cells (HSG, HSY, HTB-41) as well as fibroblasts (HFF-1) (Fig. 5e, Supplementary Fig. 9). Significant increases in growth occurred in all cells tested following FAT4 knockdown but not in controls.

In summary, we report a large-scale ACC analysis and highlight potential driver mutations in novel genes and pathways. Common alterations across ACC include pathways involving chromatin remodeling, MYB/MYC, DNA damage, and tyrosine kinase signaling (Fig. 4). Our discovery of genomic alterations in targetable pathways suggests potential avenues for novel therapies for a typically chemoresistant malignancy. Verified ACC cell lines[71,72] are needed to further substantiate the clinical utility of mutations identified here. In total, our data provide insights into the genetic framework underlying ACC oncogenesis and establish a foundation for new therapeutic strategies.

URLs. MSKCC cBioPortal for Cancer Genomics, http://www.cbioportal.org/; dbSNP, http://www.ncbi.nlm.nih.gov/projects/SNP/; ESP5400, http://evs.gs.washington.edu/EVS/; 1000 Genomes Project, http://www.1000genomes.org/; IPA, http://www.ingenuity.com/; Pathway Commons, http://www.pathwaycommons.org/.

ONLINE METHODS

Tumor samples

Primary tumor samples and matched normal specimens (peripheral blood or tissue) were obtained with written informed consent per approved institutional review board protocols (Supplementary Note). Specimens were snap-frozen in liquid nitrogen immediately upon surgical resection and stored at −80°C. Hematoxylin eosin-stained sections were prepared and diagnosis confirmed by a dedicated head and neck pathologist. Microdissection was performed to ensure >70% tumor purity. DNA was extracted using the DNeasy Blood & Tissue kit (Qiagen) and quantified via Quant-iTPicogreen dsDNA assay (Invitrogen). RNA was extracted using the RNeasy Mini kit (Qiagen) and quantified using a Nanodrop spectrophotometer (Thermo Fisher). The identities of tumor and normal DNA samples were reconfirmed by mass spectrometric fingerprint genotyping (Sequenom) as previously described.[74]

Exome capture and sequencing

Whole exome capture libraries were constructed via Agilent Sureselect 51MB target enrichment (Agilent). Approximately 2–3µg of genomic DNA from each tumor and normal specimen were sheared and ligated to barcoded sequencing adaptors. Enriched exome libraries were sequenced on a HiSeq 2000 platform (Illumina) for generation of paired-end reads (2×76bp) (Supplementary Table 1).

Analysis pipeline and mutation annotation

Reads were aligned to the hg19 (GRCh37) build using the Burrows-Wheeler Aligner (BWA)[75] (Supplementary Fig. 1). Further indel realignment, base quality score recalibration, and duplicate read removal were performed using Genome Analysis Toolkit (GATK) v2.2.[76] Somatic Sniper v1.0.0[77] and GATK Somatic Indel Detector[76] were used to generate single nucleotide variations (SNVs) and indel calls, respectively, via standard default parameters. MuTect[78] was used to confirm SNV calls with conditions described previously.[79] Baseline filters (depth 3x coverage in both tumor and normal specimens, >97% normal allelic fraction, >10% tumor allelic fraction) were chosen. We excluded germline variants found in the 1000 Genomes Project[80], ESP5400 (NHLBI GO Exome Sequencing Project), and dbSNP132.[81] Further filtering was based upon mapping quality score, average read length, and strand bias as previously described.[77] Resulting mutations were annotated based on RefSeq (Release 55) using Annovar.[82] All candidates were manually inspected via Integrative Genomics Viewer (IGV) v2.1[83] (Supplementary Fig. 4). Mean somatic mutations per sample (16 nonsynonymous, 6 synonymous) are comprised of exonic SNVs and indels. The non-silent mutations per megabase (MB) rate (0.31) refers to nonsynonymous mutations found within the exon capture target regions (51MB).

Validation of mutations

Nearly all SNVs or indels were resequenced via SOLiD (Life Technologies) or MiSeq (Illumina) platforms to validate putative mutations. Regions encompassing putative candidates were amplified with PCR using Kapa Fast HotStartTaq Polymerase (KapaBiosystems) with WGA genomic DNA generated via Repli-G kits (Qiagen). For SOLiD preparation, 3µg of DNA was sheared, with 200–220 bp fragments selected by agarose gel electrophoresis. DNA was purified and ligated to P1 and P2 adaptors per manufacturer’s instructions. Following PCR amplification, DNA was quantified and size confirmed by Bioanalyzer DNA 1000 assay (Agilent Systems). Approximately 300–500 ng per sample was used for 24hr hybridization with custom SureSelect baits (Agilent Technologies) per manufacturer instructions. Captured DNA was separated from noncaptured DNA per SureSelect protocol. Following emulsion PCR and bead preparation, each sample was loaded into one octet of a SOLiD plate for paired-end sequencing. This was followed by BWA[75] alignment, GATK v2.2[76] variant detection, and manual inspection on IGV 2.1.

For MiSeq sequencing, amplified templates were purified using AMPure (Beckman Coulter). The purified PCR reactions were sequenced bidirectionally with M13 forward and reverse primer and Big Dye Terminator Kit v.3.1 (Applied Biosystems). Dye terminators were removed using the CleanSEQ kit (Beckman Coulter), and the resulting amplicons were barcoded and prepared for sequencing using the Nextera DNASample Prep Kit (Illumina). Libraries were sequenced using the paired-end 150-cycle protocol, followed by BWA[75] alignment, GATK v2.2[76] variant detection, and manual inspection on IGV 2.1 (Supplementary Table 3).

Whole genome capture and sequencing

Library construction and whole genome sequencing were performed as previously described.[84] Briefly, 5µg of native DNA from tumor and matched normal sample were separately sheared, standard paired-end adaptors were ligated, and products run on a gel and excised. Bands were purified via mini-elute columns (Qiagen) and enriched by PCR amplification (10 cycles). Paired-end sequencing (2×101bp) was performed on the HiSeq 2000 platform (Illumina) (Supplementary Table 8). The HumanOmni 2.5–8v1 BeadChip Kit (Illumina) was used for genotyping.

Structural variant and copy number analysis

Structural variants (SVs) were identified using CREST (Clipping REveals STructure)[85], which employs soft-clipped reads to directly map the breakpoints of SVs. Samples were analyzed using the paired analysis module. Primers were designed for the 1000bp flanking regions surrounding the predicted SV breakpoint using Primer3[86] (Supplementary Table 9). Using a standard touchdown PCR protocol, products were run on a gel, excised, purified using a NucleoSpin Gel and PCR Clean-up kit (Clontech), and sequenced. Candidates were considered germline if found in both the tumor and normal sample (Supplementary Fig. 6).

Copy number variations (CNVs) were calculated with ExomeCNV[15], a statistical method to detect somatic copy number variation using depth-of-coverage from mapped short sequence reads. A subset was also run on AffymetrixGenomeWide SNP 6.0 arrays (Affymetrix) as confirmation. Analysis was performed with the aroma.affymetrix package[87] and segmentation was performed using the CBS algorithm in the DNAcopy R package.[88] For quality control, the subset of samples genotyped on HumanOmni 2.5–8v1 BeadChips (as part of whole genome sequencing) were also analyzed for CNVs by using tQN[89] for normalization, followed by OncoSNP v1.2 for characterization of CNVs and loss-of-heterozygosity (LOH) events.[90] Copy number clustering was performed on segmented copy number data. A unified breakpoint profile (region by sample matrix) was derived by combining breakpoints across all samples and determining the minimal common regions of change. Unified breakpoint profiles were computed using the Bioconductor package CNTools.[91] Hierarchical clustering was done using the R function hclust, with manhattan distance and Ward’s agglomeration method. Gene level copy number changes across samples were identified using GISTIC2.0 (Genomic Identification of Significant Targets in Cancer)[16] (Supplementary Table 6).

Gene expression profiling

Approximately 300ng of RNA was extracted from fresh frozen tissue, and RNA integrity was assessed using the Eukaryote Total RNA Nano Assay of the Agilent 2100 Bioanalyzer (Agilent Systems). Expression profiling of viable RNA was performed using the Human HT-12 Expression BeadChip array (Illumina). Analysis was performed with Partek Genomics Suite v6.5 (Partek).

Fluorescence In Situ Hybridization (FISH)

Tissue microarrays were constructed for tumors with available formalin-fixed, paraffin-embedded (FFPE) tissue. Unstained frozen slides were generated for all other tumors. FISH was performed using a 3-color probe mix consisting of BAC clones for 5’ MYB (RP11–614H6, RP11–104D9) (Green), 3’ MYB (RP11–323N12, RP11–1060C14) (Orange), and 3’ NFIB (RP11–413D24, RP11–589C16) (Red) (BACPAC Resources and Wellcome Trust Sanger Institute[92]) (Supplementary Fig. 3). DNA was prepared by standard alkaline lysis, and labeled by nick translation with fluorochrome-conjugated dUTPs (Enzo Life Sciences). Paraffin and frozen slide sections were processed as previously described.[93]

Driver/Passenger and pathway analysis

Driver/Passenger analysis was performed with CHASM v1.0.5 (Cancer-specific High-throughput Annotation of Somatic Mutations), a random-forest-based classification algorithm that predicts whether tumor-derived somatic missense mutations are important contributors to cancer cell fitness.[14] The default set of 52 features was used with an FDR cutoff of 0.35 (Supplementary Table 4).

Pathway analyses for enrichment were performed using IPA (Ingenuity Pathway Analysis) (Winter Release 2012, Ingenuity Systems), DAVID v6.7[94] (Database for Annotation, Visualization and Integrated Discovery), and hypergeometric distribution, with potential CHASM drivers used as input. Enriched genesets reflecting Gene Ontology[95] biological processes, Kyoto Encyclopedia of Genes and Genomes (Release 64.0), and Biocarta pathways were corrected for multiple hypothesis testing using the Bonferroni method. Aberrant signaling pathways (Fig. 4) were constructed referencing IPA and Pathway Commons.[38] Expression data was analyzed using Gene Set Enrichment Analysis (GSEA) v2.0.10.[56] Functional consequences were also evaluated using the MSKCC cBioPortal for Cancer Genomics.[96,97]

Cell culture

Cell lines were cultured with the following: 293T and COS7 cell lines, DMEM 10% FBS; HSG, DMEM 5% FBS; HSY, DMEM:F12 10% FBS; HTB-41, McCoy 10% FBS; HFF-1, DMEM 15% FBS. All cell lines were obtained from ATCC except for HSG (O. Baker) and HSY cells (R. Wong).

FAT4 knockdown, transfection, and viral transduction

FAT4 knockdown was performed using siRNA or shRNA (Dharmacon). siRNAs were transfected into cells in antibiotic-free medium using Lipofectamine RNAiMAX (Invitrogen), medium was changed after 6h and cells were harvested after 48h. shRNAs were transfected into 293T cells, along with packaging and envelope plasmids. Virus was harvested from cell culture supernatant at 48, 72, and 96h. Cells were transduced with 8 µg/mL polybrene. Two independent RNAi sequences were used for knockdown. Relative knockdown fold change was normalized to respective controls.

KDM6A mutagenesis and overexpression

Flag-tagged KDM6A plasmid was obtained (AddGene, #17438) and site-directed mutagenesis performed with QuikChange II XL (Stratagene). 293T and COS7 cells were transfected with FuGene HD (Promega).

Immunofluorescence

COS7 cells were plated on poly-d-lysine cover slips (BD Biosciences) and transfected with wildtype or mutant FLAG-tagged KDM6A. After a 36h incubation, cells were fixed in 4% paraformaldehyde in PBS, permeabilized in 0.5% Triton X-100 in PBS, and blocked in 10% FBS. Primary antibodies were anti-FLAG (Sigma, F3165; 1:50 dilution) and anti-H3K27me3 (Millipore, 07–449; 1:50 dilution). Secondary antibodies were Alexa Fluor 488 anti-mouse IgG and Alex Fluor 568 anti-rabbit IgG (Molecular Probes; 1:200 dilution). Images were acquired using a SP5 confocal microscope (Leica).

Growth curves

Growth curve assays were performed in triplicate and quantified using the Vi-Cell XR Cell Viability Analyzer (Beckman Coulter) or in real time in quadruplicate with the xCELLigence System (Roche). xCELLigence plates were seeded with 5,000 or 20,000 cells per well, and growth is reported as the cell index, a measure of impedance reflecting viable, adherent cells.

TP53 Signaling Pathway PCR Array

Approximately 200 ng of RNA extracted from fresh frozen tissue was used to make cDNA using SuperScript III (Invitrogen). Gene expression was quantified using the RT2 Profiler PCR Array Human p53 Signaling Pathway (SABioscience). The reactions were performed using a Mastercycler™ ep Realplex 4s and Realplex software (Eppendorf). The 2-(ΔΔCt) method was used to calculate the ΔΔCt values. In all tissue-specific comparisons, individual gene expression levels in tumor specimens were compared with that in normal salivary tissue. Comparisons were made within tissue between TP53-pathway wildtype and TP53-pathway altered tissue. All 84 genes were normalized against the average of two most stably expressed reference genes (GAPDH, HPRT1). Gene expression was considered below the limit of detection if the Ct value was 35.0 or greater.

PI3K immunohistochemistry staining

5µm FFPE tissue sections were stained for pAKT (Cell Signaling Technology; 1:100 dilution) on Benchmark Ultra (Ventana) using the OmniMap DAB anti-Rabbit Detection Kit (Ventana, 760-149) and multimer detection system. 5µm FFPE tissue sections were similarly stained for pPRAS40 (1:40 dilution) using biotinylated secondary antibody (Vector Labs, no. BA-1000; 1:300 dilution) and straptividin detection system (Ventana).