Drug-Sensitivity Screening and Genomic Characterization of 45 HPV-Negative Head and Neck Carcinoma Cell Lines for Novel Biomarkers of Drug Efficacy

Head and neck squamous cell carcinoma (HNSCC) is a common type of human cancer worldwide (1). Local HNSCCs are often treated with surgery, radiotherapy, and/or chemoradiotherapy, but cancer recurrence is frequent, and emergence of distant metastases invariably leads to death (2, 3). Systemic treatments for advanced HNSCC often consist of chemotherapeutic drugs, such as platinum salts and taxanes, but these drugs have only moderate efficacy (4). The only approved targeted drugs for advanced HNSCC are cetuximab, a monoclonal antibody that binds to the epidermal growth factor receptor (EGFR), and, recently, the anti-PD-1 antibodies nivolumab and pembrolizumab (5). Despite EGFR is frequently amplified and overexpressed in HNSCC, cetuximab has shown only moderate efficacy (6, 7), and there are no clinically approved predictive biomarkers for cetuximab response. There is a clear unmet need for novel, effective targeted treatments and accompanying biomarkers for efficacy.

Recent comprehensive genomic analyses of HNSCCs have greatly advanced our understanding of the disease and the key molecular alterations that drive these cancers (8–13). However, the translation of the molecular aberrations into treatment strategies for patients is challenging as the efficacy of many anticancer drugs is influenced by the genomic and transcriptomic heterogeneity of tumors (14). Numerous novel anticancer drugs are entering the clinical phase of testing, and testing all of them and their potential combinations in clinical trials is challenging. Combining molecular and genetic analyses with comprehensive preclinical compound screening in vitro might reveal new drug targets, and could facilitate rational drug testing in clinical trials.

HNSCCs comprise both human papillomavirus (HPV)-positive and HPV-negative cancers, the latter being more common. Patients with HPV-negative cancer have, in general, less favorable survival as compared with those with HPV-positive HNSCC (15). Recurrent or metastatic HPV-positive cancers resemble genetically HPV-negative tumors (16). In the current study, we performed a large-scale drug screen on 45 HPV-negative HNSCC cell lines, and integrated the compound sensitivity landscapes with detailed genomic and gene expression analyses to identify novel drugs that might be effective for HNSCC along with predictive biomarkers for these drugs.

The series consists of 45 HPV-negative HNSCC cell lines established from HNSCC (Table 1). The cell lines were established at the Department of Otorhinolaryngology-Head and Neck Surgery, Turku University Hospital (Turku, Finland) as described elsewhere (17, 18). Thirty-one (69%) of the 45 patients whose cancer gave rise to a cell line were male, and the mean age was 62 years at the time of tissue sampling. A majority (36 of 45, 80%) of the cell lines originated from a primary tumor, 4 (9%) from a primary tumor that persisted after treatment, and 5 (11%) from recurrent cancer. Most of the cell lines (n = 33, 73%) originated from an oral cavity or oropharyngeal tumor, 11 (24%) from larynx cancers, and 1 (2%) from a maxillary sinus cancer. All cell lines tested negative for mycoplasma before the experiments using a VenorGeM Mycoplasma PCR Detection Kit (Minerva Biolabs GmbH).

The cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 2 mmol/L L-glutamine, 10% fetal bovine serum, nonessential amino acids solution, penicillin (100 U/mL), and streptomycin (100 μg/mL) at 37°C in an atmosphere of 5% CO₂ (all reagents were purchased from Lonza BioWhittaker). Genotyping for HPV 6, 11, 16, 18, 26, 31, 33, 35, 39, 42, 43, 44, 45, 51, 52, 53, 56, 58, 59, 66, 68, 70, 73, and 82 strains was performed using a Multiplex HPV Genotyping Kit (DiaMex) in PCR-amplified samples of genomic DNA isolated from the cell lines. All cell lines were negative for all the HPV strains tested. Normal skin fibroblasts were obtained from 3 patients whose tumors gave rise to the UT-SCC-19A, UT-SCC-29, and UT-SCC-45 cell lines, and they were cultured as the cancer cells except that 20 mmol/L HEPES (Invitrogen) was added to the cell culture media. CCD 1106 KERTr (ATCC CRL-2309) keratinocytes were grown in keratinocyte-serum free medium supplemented with epidermal growth factor and bovine pituitary extract (Invitrogen).

Genomic DNA was isolated from the cell lines using a DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's protocol. The concentration and purity of the genomic DNA was estimated with a NanoDrop spectrophotometer (Thermo Scientific). Total RNA was isolated using a Macherey-Nagel (Düren) Nucleospin RNA II Kit. The quality of total RNA was assessed with a 2100 Agilent Bioanalyzer (Agilent Technologies), and only samples of high quality (2100 Bioanalyzer RNA integrity value >7) were included in the analyses.

An arrayCGH was carried out using an Agilent Human Genome CGH 244A Oligo Microarray. The samples were preprocessed according to an Agilent protocol (version 6.0, November 2008; Agilent Oligonucleotide Array-Based CGH for Genomic DNA Analysis). The gene copy-number data were processed as described earlier (19). Briefly, the data were LOESS normalized, segmented with circular binary segmentation using an R-package DNAcopy (20), and the regions with segment means more than two standard deviations from the mean (±0.67) of three self-versus-self control cell lines were called altered. The arrayCGH data are available under a Gene Expression Omnibus (GEO) accession number GSE108062. The significance of frequent alterations was assessed with the GISTIC v.2.0.21 software (the parameters used were q < 0.25, tamp = 2^0.67, tdel = 2^−0.67, cap = 2, join segment size = 20). Copy-number polymorphic loci were excluded before data segmentation. The data were analyzed with the Anduril software framework (21). As due to the stringent criteria in arrayCGH analysis some previously reported alterations might be missed, the data were checked manually for consistency for selected regions with reported alterations.

Gene expression analysis was carried out on Illumina HT-12 human microarrays (Illumina Inc.). The RNA samples from the 45 cell lines were processed according to the manufacturer's recommendations. The gene expression data were quantile normalized using the lumi R-package (22). The gene expression data are available under a GEO accession number GSE108062. Downstream analysis of the gene expression data is described in detail below.

A custom head and neck cancer panel was created using a Roche NimbleGen (Roche NimbleGen, Inc) Comprehensive Cancer Design consisting of 578 genes as a baseline. We customized this panel by analyzing these 578 genes in the cBioPortal (23, 24) for published mutations in head and neck cancers. In addition, information about the gene sequence alterations in head and neck cancers was searched from the PubMed. The genes with no known mutations in head and neck cancer were removed from the panel, and genes related to head and neck cancer were added, resulting in a head and neck cancer customized panel consisting of 694 genes (Supplementary Table S1). The gene list of the 694 genes was provided to Roche NimbleGen (Roche NimbleGen, Inc) for a SeqCap EZ Choice custom panel design. The exon locations were collated from ENSEMBL (v75), CCDS (29Nov2013), VEGA (v54), and RefSeq (13Aug2013). Overlaps were collapsed, and a single unified list of target exons was generated. These coordinates were further adjusted by expanding any targets less than 100 bp to a minimum length of 100 bp, centered in the middle of the exon. Any overlaps were then consolidated a second time. This coordinate file was used as the input for capture probe selection for a unique probe set containing unique probes as determined by the SSAHA (Sequence Search and Alignment by Hashing Algorithm). A probe was considered to match to the genome if there were five or fewer single-base insertions, deletions, or substitutions between the probe and the genome.

The genomic DNA was sheared with an S2-focused ultrasonicator (Covaris, Inc.). The library was prepared with a ThruPLEX DNA-seq Kit (Rubicon Genomics), and the targets were enriched using the custom SeqCap EZ Choice Library (Roche NimbleGen, Inc.). The samples were sequenced with an Illumina HiSeq sequencing system (Illumina Inc.) using a 100-nt paired ends read length. We reached an average coverage of 197 at the target regions.

The samples were aligned to hg19 using “bwa mem” (v. 0.7.10) with default parameters (25) and sorted with the Picard software (http://broadinstitute.github.io/picard). The targeted sequencing protocol greatly increases the likelihood of pseudoduplicated reads (i.e., reads that appear duplicated) and distinguishes real and false PCR duplicates. PCR duplicates were therefore not removed. Mutations were called using the VarScan v.2.3 software using parameters P = 0.95, somatic P = 0.05, strand-filter = 1, and min-var-frequency = 0.2 (26). Each cancer cell line was compared against the pooled sequences from the three control cell lines (normal skin fibroblasts from the same patient as UT_SCC_19A, 29 and 45). The germline variants included in the dbSNP and the 1000Genomes databases were removed after calling.

We used a compound library that consisted of 220 FDA-approved drugs as well as investigational compounds. A large proportion of the compounds were classified as signal transduction inhibitors (n = 94), mostly targeting either receptor or cytoplasmic kinases (Supplementary Table S2). Other classes that contained over 10 different compounds included alkylating agents, hormonal agents, antimetabolites, angiogenesis inhibitors, and topoisomerase inhibitor. Plerixafor and platinum salts including cisplatin and carboplatin were dissolved in water as dimethyl sulfoxide (DMSO) inactivates them while the rest of the compounds in the drug-sensitivity and resistance testing (DSRT) screening platform were dissolved in DMSO. The DSRT was carried out as described previously (27). Briefly, all compounds were tested over a 10,000-fold concentration range in five different concentrations in 384-well plates to generate quantitative and reliable dose-response data. Drugs were dispensed on plates by acoustic liquid handler (Echo 550, Labcyte) and dissolved in 5 μL of media. The cells were plated at a density of 1,000 cells/well in 20 μL volume by MultiDrop Combi peristaltic dispenser (Thermo Scientific) and incubated at 37°C and 5% CO₂ for 72 hours. After 3 days cell viability was measured by CellTiter-Glo luminescent assay (Promega). DMSO and benzethonium chloride (100 μmol/L) were used as a negative and positive control, respectively.

The dose-response curves were modeled using a four-parameter logistic function, defined by the top and bottom asymptotes, the slope, and the half-maximal inhibition constant (IC₅₀). In the curve fitting, the top asymptote of the curve was fixed to 100% viability, while the bottom asymptote was allowed to float between 0 and 75% (i.e., compounds causing less than 25% inhibition were considered inactive). To quantitatively profile the individual HNSCC cells in terms of their compound responses, we calculated the model-based drug-sensitivity score (DSS), which has been shown to improve the identification of cellular addictions and other druggable vulnerabilities in individual cancer patients over IC₅₀ and other univariate response parameters (27, 28). The differential DSS (dDSS) was calculated by comparing the DSS of particular HNSCC cells to the mean DSS over all the cell lines. Unsupervised clustering of the DSS profiles across the tested compounds and HNSCC cell lines was performed using the Ward hierarchical complete-linkage clustering algorithm using Spearman and Manhattan distance measures for the drug and sample profiles, respectively.

Drug classes were defined as follows: EGFR inhibitors consisted of afatinib, canertinib, gefitinib, and erlotinib; MEK inhibitors of refametinib, binimetinib, selumetinib, trametinib, pimasertib, and TAK-733; and mTOR inhibitors of OSI-027, AZD8055, sirolimus, temsirolimus, ridaforolimus, and everolimus. For each of these drug classes, we grouped the six cell lines with the highest class-specific DSS (median drug sensitivity calculated from the compound-specific DSSs in the drug class) into a high-sensitivity group, and the six cell lines with the lowest class-specific DSS into a low-sensitivity group. We compared the difference in the gene expression between the high and low-sensitivity groups, and the genes with a fold change of > 2 (group medians) were considered significantly differentially expressed (t test, P < 0.01). The P values for the significantly differentially expressed genes were adjusted by the Benjamini–Hochberg method (29) using “stats” R-package (R Foundation for Statistical Computing).

To assess the influence of a gene-specific mutation on drug sensitivity, the cell lines with a mutation were compared with the rest of the cell lines using the Wilcoxon rank-sum test (two-sided). The difference between two groups was assessed with permutation test using “coin” R-package (30). The P values were adjusted with the Benjamini–Hochberg (29) method as described above. Due to the functional relevance of stop-gain mutations to protein product, the identical analysis was performed to assess the influence of gene-specific stop-gain mutations on drug sensitivity. We restricted our analysis only to the most common mutations we observed in a gene panel of 694 genes.

The difference in the drug sensitivity between the cell lines with or without a copy-number alteration (CNA) was tested using the Wilcoxon rank-sum test, and the difference was considered significant at a false discovery rate (FDR) level q < 0.05 (Storey correction for multiple testing; ref. 31). The differences in drug sensitivity were quantified by subtracting the median DSS of the samples with a CNA from the median DSS of the samples without the CNA. The statistical testing was done for each drug separately. We further compared the gene expression data with the CNA data to find the genes whose expression is modified due to a CNA. By using the CNA data, we categorized the genes located in the same chromosome region into two groups, altered and nonaltered, and calculated the signal-to-noise ratio (32). The significance of the difference in gene expression between the two groups was assessed with a permutation test and corrected for multiple hypothesis using the Benjamini–Hochberg correction (29). Genes with an FDR q < 0.05 were considered significantly altered. Only CNAs that occurred in at least 10% of the cell lines were tested. Finally, the genes whose CNA status was significantly associated with both gene expression and drug sensitivity were further tested for an association between gene expression and drug-sensitivity profiles using the bootstrapCor function in the maigesPack R-package (−0.3 ≥ |r| ≥ 0.3; ref. 33) and adjusted with the Benjamini–Hochberg (29) method.

mRNA was reverse transcribed to cDNA using an iScript cDNA synthesis kit (Bio-Rad) following the manufacturer's instructions. RT-qPCR analyses were done using the SYBR Green dye (LightCycler 480 DNA SYBR Green I Master, Cat. No. 04707516001, Roche Diagnostics GmbH). Quantitation was performed using LightCycler 480 384-multiwell plates and the LightCycler 480 instrument (Roche Diagnostics GmbH). The analyses were done in a triplicate for each gene. The primers used are listed in Supplementary Table S3. The data were analyzed with a Basic relative quantification program of the LightCycler 480 Software v.1.5 (Roche Diagnostics GbmH), based on the 2-ΔΔCT-method. A keratinocyte cell line, CCD 1106 KERTr (ATCC CRL-2309), served as the reference sample for gene expression normalization. GAPDH was selected as the reference gene for normalization.

To validate selected results from the drug screen, we silenced FAM83H (pool of siRNAs s50090, s50091, and s50092, Ambion) and RAB25 (pool of siRNAs SI03036544, SI02644089, SI02644082, and SI00126084, Qiagen) in UT-SCC-14 and UT-SCC-8 cell lines, respectively. A total of 12 nmol/L pooled siRNAs were used in three replicate experiments. Qiagen All Stars Negative (#1027281) served as a negative and Qiagen All Stars Hs Death Control (#1027299) as positive control. Cells were reverse transfected using Lipofectamine RNAiMAX (Thermo Fisher Scientific), and 24 hours later supplemented with nine different afatinib or erlotinib concentrations (0.1, 0.3, 1, 3.2, 10, 31.6, 100, 316.2, and 1,000 nmol/L or 1, 3.2, 10, 31.6, 100, 316.2, 1000, 3162.3, and 10,000 nmol/L, respectively). After 72 hours of drug treatment, the number of viable cells was measured using the CellTiter-Glo luminescent cell viability assay (Promega). Percent inhibition was calculated using the negative and positive controls using formula:

To account for the siRNA effect on the cells, the percent inhibition value at the lowest concentration of each dose response was treated as the siRNA background, and this value was subtracted from values of the other concentrations. The percent inhibition data for each drug were then curve fitted with Levenberg–Marquardt method. The replicates were averaged for each concentration during the curve fitting. The IC₅₀, slope, maximal and minimal response coming out of the curve fitting are then used to calculate the DSS (28).

A high-throughput DSRT platform of 220 compounds showed selective responses across the 45 HNSCC cell lines (Supplementary Table S4) as measured with the DSS, a modified form of the area under the curve calculation that integrates multiple dose-response parameters for each drug. DSS therefore captures both the potency and the efficacy of the drug effect (Supplementary Table S4). Unsupervised hierarchical clustering of drug sensitivity over all 45 cell lines yielded eight clusters (Supplementary Fig. S1). Of these, group I was unresponsive to most of the 220 drugs tested. Three groups showed selectivity to EGFR and MEK inhibitors (groups II, VI, and VIII) as single agents, proposing synergistic effects for selected cell lines or in combination with an mTOR inhibitor (group VIII). Group V cell lines were responsive to mTOR inhibitors and mostly unresponsive to EGFR and MEK inhibitors.

The activity of EGFR tyrosine kinase inhibitors afatinib, canertinib, erlotinib, and gefitinib had a median DSS of >10, with canertinib (a pan-ErbB inhibitor) showing the highest activity (Supplementary Fig. S2; Supplementary Table S4). Many of the genomic alterations found in HNSCC may result in a persistent activation of the PI3K–AKT–mTOR pathway (11). In accordance with this, the cell lines often responded well to five of six tested mTOR inhibitors (AZD8055, everolimus, ridaforolimus, sirolimus, and temsirolimus) and to all tested PI3K/mTOR inhibitors (apitolisib, dactolisib, omipalisib, and PF-04691502). The median DSS to these drugs was >10 (Supplementary Figs. S1 and S2; Supplementary Table S4). Interestingly, all the other PI3K inhibitors except pictilisib did not show substantial efficacy (the median DSS <10). MEK inhibitors had variable efficacy, with a DSS >10 obtained with pimasertib, refametinib, TAK-733, and trametinib.

The cell lines were subjected to a genome-wide analysis of gene copy-number gains and losses using an Agilent Human Genome CGH 244A Oligo Microarray and targeted to exome sequencing of 694 genes to correlate the genomic alterations with drug sensitivity. In general, the detected mutations, which occurred frequently in TP53, NOTCH1, FAT1, CSMD3, and CDKN2A, and the commonly observed CNAs (amplifications at 3q26, 7p11.2, 8q24, and 11q13 and deletions at 3p14, 8p23, 9p21, 9p23, and 18q), resemble those observed in clinical HNSCC (refs. 10–14, 34; Fig. 1; Supplementary Tables S5 and S6). In total, we detected 552 variants in 279 genes in the targeted HNSCC gene panel by next-generation sequencing (Supplementary Table S5). The number of variants among the studied 694 genes differed between the cell lines (range, 3–29; mean, 12). The UT-SCC-60A cell line had the largest and the UT-SCC-99 and UT-SCC-36 the lowest number of mutations. The number of mutations did not correlate with the primary tumor tumor–node–metastasis classification size category (P = 0.379; unpaired t test) or with cancer regional lymph node involvement (P = 0.887; unpaired t test).

The genes that were mutated in at least five (≥10%) of the cell lines were compared with drug sensitivity (using the DSS) across the 45 cell lines. The results are shown in Supplementary Table S7. The analysis between CNA alterations and drug responses revealed 773 CNAs that were significantly associated with the drug sensitivity (Wilcoxon rank-sum test, Q < 0.1; Supplementary Table S8).

The influence of a CNA to gene expression was assessed by expression profiling using Illumina HT-12 microarrays. Of the 4,298 genes (1,952 deleted, 2,347 gained) that were located in the regions altered in at least 10% of the cell lines, 472 showed significant downregulation and 445 significant upregulation based on CNA cis effects at their loci (Supplementary Tables S9 and S10). The amplified and upregulated genes included PIK3CA at 3q26.3, EGFR at 7p11.2, CCND1, ORAOV1, PPFIA1, CTTN, and FADD at 11q13, and BIRC2, BIRC3, and YAP1 at 11q22. All of these genes are known to be located in frequently amplified chromosomal regions in HNSCC (10–13). The deleted and downregulated genes included CDKN2A at 9p21.3 and SMAD4 at the 18q21.2 region. A real-time quantitative PCR reaction of eight genes validated the gene expression microarray data (Supplementary Table S3).

Finally, to identify the genes that might be biomarkers for drug response, within the amplified and deleted regions, we integrated drug-sensitivity data with the gene expression data. Only those genes that showed marked CNAs in-cis and substantial changes in mRNA expression were included in this analysis. Significant drug-sensitivity-gene expression-CNA triads are presented in Supplementary Table S11. The most common genetic alterations and their association with drug response to EGFR, MEK and mTOR inhibitors for the 45 HNSCC cell lines are presented in Fig. 2.

To find potential biomarkers for mTOR and MEK inhibitor sensitivity, we searched for somatic mutations and copy-number aberrations that associated with drug sensitivity from the analyses described above. In addition, the gene expression patterns of the six most sensitive and the six most resistant cell lines for these inhibitor classes were compared. These drugs were selected, because they showed high selectivity across all cell lines, and they are of clinical interest in the treatment of HNSCC.

A 24-gene signature discriminated the high and poorly responsive cell lines to six mTOR inhibitors (Fig. 3A; Supplementary Table S12), from which everolimus and temsirolimus are currently being evaluated as therapeutic options in combination with chemotherapy, radiotherapy, or other molecularly targeted drugs in HNSCC (35). One of these genes was oxysterol binding protein like 1A (OSBPL1A), a gene involved in lipid metabolism. The integrative analysis of CNA-gene expression and DSS revealed that the expression of OSBPL1A was affected by deletion in a group of poorly responsive cell lines. In comparison with the gene expression profiles of the six most sensitive and resistant cell lines for six MEK inhibitors, 11 genes, including EFNB1 and STAT6, showed higher levels of expression in the cell lines sensitive to MEK inhibitors as compared with cell lines resistant to these drugs (Fig. 3B; Supplementary Table S12).

Regarding CNAs, deletion of PTPRD at 9p24.1, present in nine cell lines, was statistically significantly associated with poor sensitivity to MEK inhibitor refametinib as assessed by Wilcoxon rank-sum test (Q < 0.1; Supplementary Table S8; Fig. 3C).

We found no significant associations between responses to mTOR or MEK inhibitors with mutations after correcting for multiple hypothesis testing (permutation test Q > 0.1). However, as specific genes have mutations in relatively small number of cell lines in our cohort, a larger series of samples is needed to further assess the statistical significance of the alterations. As an example, UT-SCC cell lines that harbored a stop-gain mutation in FAT1 (present in 9% of the cell lines) showed a tendency for a higher sensitivity to temsirolimus as compared with the other cell lines (Wilcoxon rank-sum test P < 0.05, Fig. 3D; Supplementary Table S7).

Mutations in NOTCH1, which is one of the most commonly mutated genes in HNSCC, were associated with a poor response to refametinib (Wilcoxon rank-sum test P < 0.05; Fig. 3C and D). Interestingly, mutations in NAV3 were associated with poorer efficacy of four of the six MEK inhibitors tested (Wilcoxon rank-sum test P < 0.05; Supplementary Table S7; Fig. 3C). Cell lines with CDKN2A mutation, which occurred in 20% of our cell lines, had tendency for poor sensitivity to the MEK inhibitor trametinib (Wilcoxon rank-sum test P < 0.05, Fig. 3D; Supplementary Table S7). On the other hand, CDKN2A deletions or decreased expression did not show significant associations with drug response.

EGFR status has not been considered as a reliable biomarker for response to EGFR inhibitors in HNSCC (36). In our study, cell lines with EGFR amplification and overexpression were weakly associated with sensitivity to selected inhibitors, such as erlotinib (r = 0.36) and canertinib (r = 0.30; Supplementary Table S11). Most cell lines with increased EGFR copy number showed sensitivity to EGFR inhibitors, but in line with previous studies, the result was not conclusive (Fig. 4A).

Interestingly, when we compared the gene expression between cell lines responsive and resistant to EGFR inhibitors, the cell lines sensitive to EGFR inhibitors showed high mRNA levels for several genes associated with membrane trafficking and transport, such as RAB25, CLIC3, SYK, SLC31A2, LYPD3, TMEM125, and TSPAN1 (Fig. 4B and C). Knockdown of RAB25 slightly reduced the sensitivity of UT-SCC-8 cells to erlotinib and afatinib (Supplementary Fig. S3A–S3D and S3I).

In the integrative analysis of drug sensitivity, gene expression, and gene copy numbers, FAM83H overexpression, which was caused by gene amplification, correlated with enhanced sensitivity to EGFR inhibitors afatinib (r = 0.53, permutation P = 0.0001, Q = 0.2, Fig. 4D), erlotinib (r = 0.53), canertinib (r = 0.51), and gefitinib (r = 0.49; permutation P < 0.001, Q > 0.2; Supplementary Table S11). This gene is amplified in 12% of HNSCCs (13). Moreover, FAM83H silencing made UT-SCC-14 cells more resistant to afatinib and erlotinib (Supplementary Fig. S3E–S3I).

To validate the gene expression results for RAB25, TMEM125, and FAM83H in an independent data set, we determined drug-sensitivity calls (sensitive, intermediate, or resistant) in the Cancer Cell Line Encyclopedia (CCLE; ref. 37), which provides genetic and expression data for 490 cell lines together with their response to EGFR inhibitor erlotinib. We used 3-bin k-means clustering of the −logIC₅₀ values. The sensitivity to erlotinib and expression of RAB25, TMEM125, and FAM83H were significantly associated in the CCLE data set (Fig. 4E).

Similar to MEK inhibitors, mutations in NAV3 and NOTCH1 showed tendency with a poor response to EGFR inhibitors (P < 0.05; Supplementary Table S7; Fig. 4D) discriminating a group of cell lines with poorer responses to both EGFR and MEK inhibitors. In addition, several CNA alterations were associated with response to EGFR inhibitors afatinib (Fig. 4D), canertinib, gefitinib, and erlotinib (Supplementary Table S8). Cell lines with MYC amplification showed tendency for better response to EGFR inhibitors gefitinib, afatinib, and erlotinib as compared with the other cell lines (Supplementary Table S8). FAM83H was coamplified with MYC (Fig. 4D), which might explain the association between MYC amplification and EGFR inhibitor response. In line with these results, there was also a significant correlation between FAM83H and MYC amplification in The Cancer Genome Atlas (TCGA) HNSCC study (13) when we analyzed the cBioPortal database (Fisher exact test, P < 0.001).

Apart from cetuximab and immune checkpoint inhibitors (6, 7, 38), no other targeted drugs to treat advanced HNSCC are currently available. We studied 45 HNSCC cell lines established from HPV-negative HNSCCs on a high-throughput drug-testing platform (DSRT) and combined the drug-sensitivity information with detailed genomic and transcriptomic data, which allowed us to find novel molecular associations for drug responses. We found several genomic alterations such as mutations in NOTCH1 and FAT1 and deletions in PTPRD1, which showed tendency with poor efficacy of EGFR and MEK inhibitors. In addition, overexpression of FAM83H and amplification of MYC showed tendency for improved sensitivity to EGFR inhibitors. To our knowledge, the present study is the largest study on HNSCC that integrates genetic information with drug responses.

A limitation of this study is that only cell lines were tested, and the number of single alterations was low, affecting the power of the statistical analyses. However, testing of a large number of compounds and biomarkers is challenging in clinical trials or in xenograft models due to financial and other constraints. Repeated expansion of the cell lines likely results in adaptive changes in gene expression, but many of the cell lines used in this study had low passage numbers. Prior studies on cancer cell lines generally agree with their value in drug discovery for molecularly targeted therapies (37, 39–43). Genomic characterization of the cell lines corresponded well with the mutation spectrum observed in human HNSCCs (13), suggesting that the cell lines recapitulate oncogenic alterations observed in clinical tumor samples. Recognizing that the translation of potential associations detected in preclinical models to clinically useful biomarkers may be challenging and requires validation (42, 43), the integration of data obtained from DSRT and genomic characterization appears a powerful strategy to identify systemic treatment options and potential biomarkers (27, 28).

EGFR inhibitors are currently the only targeted kinase inhibitors that have been approved for the treatment of advanced HNSCC and have modest efficacy. In HNSCC, EGFR expression or an increased EGFR gene copy number has been extensively studied as biomarkers for response to EGFR inhibitors, but they are not predictive or prognostic (36). In our study, EGFR amplification or overexpression was weakly associated with sensitivity to selected EGFR inhibitors, but in line with the previous studies, the results were not conclusive. On the other hand, cell lines sensitive for EGFR inhibitors had high mRNA levels for genes involved in membrane trafficking and receptor and membrane transport, suggesting a role for the plasma membrane protein organization and EGFR trafficking in drug response. For example, the EGFR inhibitor–sensitive cell lines showed 8 times higher RAB25 mRNA levels as compared with the insensitive cell lines. In agreement with this finding, Rab25 was found to influence EGFR endocytosis and gefitinib efficacy (44). Interestingly, the cell lines with MYC amplification responded significantly better to EGFR inhibitors as compared with other cell lines. In a series with non–small cell lung cancer patients, a high MYC copy number correlated with high sensitivity to EGFR inhibitors (45).

An example of a novel molecular association with EGFR inhibitor response in the present study is amplification and mRNA overexpression of FAM83H. The FAM83 proteins were recently identified as novel transforming oncogenes that promote resistance to EGFR tyrosine kinase inhibitors in breast cancer (46). However, suggesting that the effects of FAM83 proteins on drug sensitivity might vary for different FAM83 family members and might be tissue-type dependent, the detailed molecular mechanisms that explain these associations remain to be elucidated.

Head and neck cancers harbor a wide spectrum of genetic alterations, and observations from several studies point to the importance of mutations in TP53, NOTCH1, CDKN2A, PIK3CA, FBXW7, and HRAS in the molecular pathogenesis of HNSCC (10–13). In the present study, NOTCH1 was the second most commonly mutated gene (present in 29% of the cell lines), and the mutations were mostly found in cell lines with poor sensitivity to EGFR and MEK inhibitors. NOTCH1 is regulated by TP53, and it has both oncogenic and tumorigenic functions (47–49). NOTCH1 is involved in squamous cell differentiation, and in keratinocytes increased NOTCH activity reduces proliferation and promotes differentiation, whereas loss of NOTCH1 promotes carcinogenesis (47–49). NOTCH1 activity and concurrent downstream activation of MYC are associated with resistance of breast cancer cells to PI3K inhibitors (50). Only 3 of the 24 cell lines harboring mutations either in NOTCH1 or in the PI3K/mTOR pathway had mutations in both, which is in line with previous clinical data showing no correlation between the PI3K/mTOR pathway and NOTCH1 mutations (11). Taken together, the results support an important role for NOTCH1 signaling in the genesis of some HNSCC and as a potential drug target and biomarker.

In summary, the current comprehensive genomic analysis integrated with high-throughput drug testing on a large panel of HNSCC cell lines identified several biomarkers for EGFR, MEK, and mTOR inhibitor efficacy in HNSCC. These biomarkers included genetic alterations such as NOTCH1 and NAV3 mutations, MYC and FAM83H amplifications, and PTPRD deletion. These potential biomarkers deserve further validation in clinical series, because robust biomarkers for EGFR, MEK, and mTOR inhibitor efficacy are lacking.

H. Joensuu is Vice President at and is a consultant/advisory board member for Orion Pharma and has ownership interest (including stock, patents, etc.) in Sartar Therapeutics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: T. Lepikhova, P.-R. Karhemo, H. Joensuu, Outi Monni

Development of methodology: T. Lepikhova, R. Grénman, K. Wennerberg

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): T. Lepikhova, P.-R. Karhemo, A. Murumägi, E. Kulesskiy, H. Sihto, R. Grénman, S.M. Syrjänen, O. Kallioniemi, H. Joensuu, Outi Monni

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): T. Lepikhova, P.-R. Karhemo, R. Louhimo, B. Yadav, E. Kulesskiy, M. Kivento, H. Sihto, T. Aittokallio, K. Wennerberg, H. Joensuu, Outi Monni

Writing, review, and/or revision of the manuscript: T. Lepikhova, P.-R. Karhemo, R. Louhimo, B. Yadav, A. Murumägi, E. Kulesskiy, M. Kivento, H. Sihto, R. Grénman, S.M. Syrjänen, O. Kallioniemi, T. Aittokallio, K. Wennerberg, H. Joensuu, Outi Monni

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): P.-R. Karhemo, M. Kivento, O. Kallioniemi, H. Joensuu, Outi Monni

Study supervision: H. Joensuu, Outi Monni

Other (coordination and detailed planning of the project): T. Lepikhova

Other (data integration): B. Yadav

Other (HPV genotyping): S.M. Syrjänen

We thank the Biomedicum Functional Genomics Unit (FuGU), FIMM High-Throughput Biomedicine Unit, and FIMM Sequencing Unit at the University of Helsinki, Finland, for excellent technical assistance.

This work was financially supported by Academy of Finland (O. Monni, H. Joensuu, T. Aittokallio, and K. Wennerberg; Center of Excellence for Translational Cancer Biology, O. Kallioniemi, H. Joensuu), Jane & Aatos Erkko Foundation (O. Monni, H. Joensuu, and K. Wennerberg), Finnish Cancer Society (O. Monni, H. Joensuu, T. Aittokallio, K. Wennerberg, and O. Kallioniemi), Sigrid Juselius Foundation (H. Joensuu, K. Wennerberg, and O. Kallioniemi), and Research Funds of Helsinki University Hospital (H. Joensuu).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.