Drug Sensitivity Prediction Models Reveal a Link between DNA Repair Defects and Poor Prognosis in HNSCC

DNA repair pathways, such as homologous recombination (HR), Fanconi anemia (FA), and associated cellular repair pathways govern genomic integrity and maintain chromosomal stability. Recent genetic studies highlight the identification of an increasing number of DNA damage response (DDR) and DNA repair gene defects in cancer (1), but the role of these defects in cancer outcome or treatment response is not yet fully understood. Early in vitro studies showed how defects in DDR and DNA repair pathways cause increased sensitivities to chemotherapeutic agents or radiation. Thus, given the DNA-damaging nature of many cancer treatment agents, tumor repair defects are expected to be beneficial for overall patient outcome. However, evidence, mostly assessed in breast cancer, is inconclusive (2). Studies in other cancer types are less numerous and are greatly hampered by the lack of biomarkers that identify functional repair defects. Here we set out to develop DNA repair defect markers to determine the clinical impact of DNA repair defects in head and neck squamous cell carcinoma (HNSCC).

DNA repair gene variants may be used to depict DNA repair defects, such as in the case of pathogenic BRCA1/2 variants in breast or ovarian cancer (1, 3). Although mutations in such genes have been shown to predict platinum (4) and olaparib (5) response, determining which mutations affect the function of a given protein or pathway remains difficult (6). In addition, despite the high number of variants and mutations, the frequency of specific variants is very low. Other attempts successfully applied indirect repair defect measures such as “genomic scars” (7). The genomic scar biomarkers were developed in BRCA1/2-mutant breast cancer and show promise for patient stratification for olaparib treatment in clinical trials in breast cancer (8). Such scar-based biomarkers depict past repair defect events. Multiple resistance mechanisms have recently been described and some show DNA repair restoration (9, 10).

The quest to assess the clinical impact of DDR and repair defects in other cancer types faces many challenges. Numerous DDR genes have been shown to be mutated in human cancer (11). However, the functional significance of most of these mutations is largely unknown. A multitude of genes are affected and the mutated genes or pathways vary largely among different cancer types. The oncogenetic context in these cancers also impacts DNA repair pathway engagement or performance. Detection of DNA repair defects therefore requires broadening the scope beyond breast- and ovarian cancer–specific BRCA mutations. A gene expression profile for mutations that phenocopy p53 inactivation was developed as a proxy for a deficiency in DNA damage response, and showed an association with poor prognosis in some cancers (12). However, an exclusive focus on TP53 mutations limits the prediction to altered p53 response rather than DNA repair defect prediction. Functional analyses, such as the inability of cells to recover from DNA crosslinker exposure, reliably reveal DNA crosslink repair defects, independent of the underlying genetic cause. Similar to genomic scars, functional repair endpoint studies aim to improve the association with true yet undetected DNA repair defects (13–15). Previous studies by us and others exposed functional DNA repair defects in the HR/FA pathway in HNSCC (16–18). These defects were identified by determining sensitivity to the two DNA-damaging agents: mitomycin C (MMC) or olaparib in patient-derived HNSCC cell lines (16). MMC is a potent DNA crosslinking agent. Cellular survival after exposure to DNA crosslinking agents requires the concerted action of multiple repair pathways, including nucleotide excision repair, FA and HR (19). MMC hypersensitivity is therefore a hallmark of FA pathway defects and has been used in the identification of many members of the FA and HR pathway, such as the BRCA1/2 genes (20). The PARP inhibitor olaparib interferes in PARP-mediated DNA repair, primarily single-strand break repair and inhibitor-induced trapped PARP poses replication blocks (21). Together with the increased load of unresolved DNA damage, such blocks also contribute to the induction of secondary DNA double-strand breaks (22). Cells with HR/FA defects and with high levels of replication stress struggle to survive PARP inhibition. As with MMC, hypersensitivity to PARP inhibition can therefore reveal functional HR and FA pathway defects, also of unknown genetic causes. Notably, a shift in the engagement from HR toward other DSB repair pathways can circumvent defects that cause PARP inhibitor sensitivity, while maintaining crosslinker sensitivity (9, 20, 23). Functional pathway defects were accompanied by genetic variants in multiple FA and HR genes in some of the HNSCC cell lines, but not all, in our in vitro study, further arguing for functional endpoint-based assays (16).

HNSCC are in particular interesting in that they have been shown to exhibit functional DNA repair defects that are caused by a multitude of genetic aberrations in different DNA repair pathways. FA and HR gene variants were identified in larger HNSCC sequencing efforts (18) and found to be common in young patients with HNSCC (24), but the relation to treatment outcome is less well explored. A small study compared HR pathway engagement in biopsies of 13 patients with HNSCC before and after the start of chemoradiotherapy by analyzing the ability to form Rad51 foci; this points to a poorer patient outcome of carriers of DNA repair–affected tumors (14). In line with such an influence, a recent study highlights the association between chromosomal instability and lymph node status in HNSCC (25). Larger studies investigating the effect of functional DNA repair defects on patient outcome in HNSCC are lacking.

In this study, we aimed to determine the impact of such DNA repair defects on outcome and treatment response in HNSCC. HR and FA repair defects provoke replication stress in proliferating cells and a network of various compensatory DNA repair options are engaged to assure cellular survival. We therefore posit that DNA repair–defective cells will show a distinct transcriptional pattern that indirectly reflects the DNA damage response to unrepaired damage from endogenous sources. To embrace the variety of repair defect types and genetic contexts, we employed a large panel of patient-derived HNSCC cell lines, of which, a subset showed functional repair defects as determined by hypersensitivity to MMC and/or the PARP inhibitor olaparib and used machine learning models to predict a biological characteristic, DNA crosslink repair status. Independent in vitro datasets were used to test the performance in identifying DNA repair defects. The impact of DNA repair defects, as identified by these expression models, in clinical outcome was determined in two independent HNSCC patient cohorts uniformly receiving standard-of-care cisplatin-based chemoradiotherapy.

UT-SCC are primary tumor cell lines from patients with HNSCC and were established by Prof. Grenman (University of Turku, Turku, Finland) from whom we obtained the cell lines. Cell lines were genotyped (Eurofins) and further characterized by photographs, cell doubling times, and assayed at low passage numbers. Cells were cultured in DMEM + GlutaMAX (Gibco, Life Technologies), supplemented with 10% FBS (Gibco, Life Technologies), 1% penicillin–streptomycin (Gibco, Life Technologies) and 1× MEM non-essential amino acids solution (NEAA; Gibco, Life Technologies) at 37°C, 5% CO₂, and 3% O₂. Cell lines were routinely tested for Mycoplasma contamination.

RNA was isolated using the DNA AllPrep DNA/RNA Mini Kit (Qiagen). Strand-specific libraries were generated using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina Inc., RS-122-2101/2) according to the manufacturer's instructions. Reads were mapped against the GRCh38 human genome with TopHat2.1, using the options fr-firststrand, transcriptome-index, and prefilter multi-hits. Gene expression was determined using HTSeq-count with options stranded and mode union. Differential expression was performed using DESeq2. Genes with log₂-fold change higher than one were classified as differentially expressed, with a false discovery rate of 10%.

Data are available at the European Genome-phenome Archive (EGA).

For modelling, RPKM values were calculated from the read counts and were scaled to have mean 0 and SD 1. Only genes whose expression was measured in the CGP (26) and Peng (27) datasets (15,383 genes) were considered to allow model validation in these datasets. In addition, expressions of 13,886 signatures from the Enrichr website (Supplementary Table S1) were calculated using the GSVA BioConductor package. Variables were selected as explained in Supplementary Fig. S1A. Models were fit using the train function from the Caret R package. Options used were trainControl (boot = 10), tuneLength = 10, and preprocessing (centre, scale). pROC package was used to calculate performance. Each of the best performing combinations of feature selection and classifier was trained on the whole dataset 20 times. For predictions, class probabilities were calculated as the average of 20 model predictions. For the model ensemble predictions, class probabilities of the selected models of each type (mmcOnly, olaOnly, or threeClass) were averaged to obtain the final prediction value. The mmcOnly model is available as an R package at bitbucket.org/paulessers/dna-repair-model-ensemble/.

Biopsies and collection of fresh-frozen tumor material were approved by the Institutional Review Board and all patients granted written informed consent. The retrospective “NKI-CRAD” cohort provided pretreatment tumor biopsy material from 98 patients with hypopharyngeal, laryngeal, or HPV-negative oropharyngeal HNSCC who were treated at the Netherlands Cancer Institute between 2001 and 2014. Samples with tumor percentage <40%, as determined by hematoxylin and eosin (H&E) staining, were excluded. HPV status was determined using p16 and p53 IHC and validated using PCR and RNA-sequencing. HPV-positive patients were excluded. All patients were treated with cisplatin-based chemoradiotherapy. Radiotherapy was administered according to a conventional (35 fractions of 2 Gy over 7 weeks) or DAHANCA scheme (35 fractions of 2 Gy over 6 weeks). Chemotherapy consisted of daily (25 times 6 mg/m²), weekly (6 × 40 mg/m²) or 3-weekly (3 × 100 mg/m²) intravenous administration of cisplatin. Eight patients received intra-arterial cisplatin (4 × 150 mg/m²) according to the RADPLAT trial (28). Not all patients completed chemotherapy and the individual total cumulative cisplatin dose was recorded. Overall survival was calculated from the start of the treatment until the event or censored at the time of the last follow-up. Metastasis and locoregional recurrence-free survival were additionally censored at the time of death. The DESIGN cohort consisted of 82 advanced HPV-negative patients with HNSCC who were included in the Dutch multicenter DESIGN study (KWF-A6C7072) and were treated at the VUmc Cancer Center Amsterdam or in the Netherlands Cancer Institute (Amsterdam, the Netherlands). Patient selection criteria were equal to those of the NKI-CRAD cohort; patients who were previously treated with radiotherapy in the head and neck area or chemotherapy were excluded.

For migration assays, cells were seeded in DMEM + GlutaMAX in an IncuCyte ImageLock 96-wells Plate (Essen Bioscience) at 100% confluence. After 24 hours, a homogeneous scratch was made with the use of the IncuCyte WoundMaker (Essen Bioscience) and washed twice with culture medium. Relative wound density was measured every 4 hours using the IncuCyte (Essen Bioscience). For invasion assays, the same protocol was used with modifications: 100 μg/mL growth factor reduced (GFR) Matrigel (Corning) was added to an IncuCyte ImageLock 96-well plate. After 24 hours, nonsolidified Matrigel was aspirated and cells seeded. After 4 hours, the cells were scratched and washed twice and 8 mg/mL GFR Matrigel was added. Thirty minutes later, culture medium was added and plates were placed in the IncuCyte. Rad51 experiments were performed by adding 500 nmol/L Rad51 inhibitor (B02, SML0364, Sigma) after cell seeding and again after the scratch was made. Inhibitors of DNA-PK NU7026, ATM Ku-55933, and WNT XAV-939 (all from SelleckChem) were added after the scratch was made.

For migration assays, multiple UT-SCC cell lines were seeded (1.11 × 10⁴ cells/insert) in culture medium with low serum concentration (0.5% FBS) in the top of an 8-μm pore size Transwell insert (Falcon, Thermo Fisher Scientific). Culture medium with 10% FBS was placed in the bottom of a 24-well plate. After 24 hours, membranes were fixed for 2 minutes at room temperature with 3.7% formaldehyde and washed twice with PBS. Cells were permeabilized with 100% methanol for 20 minutes at room temperature, washed twice with PBS, and stained with Hoechst (16 μg/mL; Thermo Fisher Scientific) for 10 minutes at room temperature. The membrane was washed with PBS and cells were removed from the bottom or top of the membrane with a cotton swab. Images were made (10× objective) with the Axiovert 200M (Zeiss) microscope and analyzed using ImageJ software. For invasion experiments, 125 μg/mL GFR Matrigel was added to the top of the insert and left open overnight, allowing the Matrigel to dry onto the membrane. The Matrigel was reconstituted with culture medium for 1 hour.

Analyses were performed in the R statistical computing environment. Scripts used for visualizing patient survival are available at https://github.com/PaulEssers/SurvivalPlots/. Spindle index, scratch assays, and transwell assays were analyzed by linear-mixed effects models using the nlme R package. Models were fit with “Group” as the fixed effect and the nested term “1|cellLine/Experiment” as the random effect. P values reported are those for the fixed effects.

A dual drug sensitivity analysis in patient-derived HNSCC cell lines revealed three classes of HNSCC: (i) MMC and olaparib sensitive, consistent with pronounced crosslink repair, i.e., HR and FA, pathway defects; (ii) MMC sensitive that lack concomitant olaparib sensitivity, likely due to compensatory repair mechanisms and a shift in replication associated repair pathway preference and (iii) HNSCC that do not show strong repair defects as identified by drug sensitivity (Fig. 1A; ref. 16). We will refer to these groups as SensMO (i), SensM (ii), and Normal (iii), respectively. Upon further analysis, HR and FA gene mutations were found in the MMC-sensitive HNSCC cell lines, providing the potential genetic causes of the functional repair defects (Supplementary Fig. S1; refs. 16, 29).

We reasoned that changes in the transcriptome will reflect cellular DNA repair status via activation of the DNA damage stress response and downstream mediators. We therefore used RNA-sequencing data from HNSCCs assigned to these three different functional endpoint classes. By pairwise comparisons, 82 genes were found to be differentially expressed between the groups. However, expression levels of these genes varied considerably within groups, and there were no genes that were expressed exclusively and ubiquitously in each group (Fig. 1B). Consequently, these traditional gene signatures, derived from these differentially expressed genes, performed poorly as predictors of drug sensitivity as determined by the AUC [area under the receiver operating characteristic (ROC) curves; Fig. 1C)]. Therefore, we used machine learning algorithms to capture more complex relationships between gene expression and the functional DNA repair endpoint.

We generated a predictive model, exploiting the combined drug sensitivity data, so as to be able to distinguish between the three groups (threeClass). Alternatively, we generated models that were based on the sensitivity to each drug separately with the goal to predict olaparib or MMC sensitivity only (mmcOnly and olaOnly). Next to gene expression values, we also calculated expression values for established gene signatures (Supplementary Table S1; ref. 30). Gene expression values, signature expression values, or both were then used as input to train the models. As our ultimate aim was to identify and investigate the effect of DNA crosslink repair defects on patient outcome, we considered a variety of methods to minimize potential method–specific bias and repair-unrelated associations. Relevant genes and gene signatures were preselected using several techniques and were used to train multiple predictive models (see Supplementary Fig. S2 for schematic representation and description). To assess the performance of all the different models derived from the individual combinations of feature selection method and classifier, we computed the mean AUC values from 10 rounds of a 5-fold cross-validation, or 16 rounds of a 3-fold cross-validation of each (Fig. 1D and E; Supplementary Table S2). Nearly all combinations predicting olaparib or MMC sensitivity (i.e., olaOnly and mmcOnly) performed well. We assessed whether the IC₅₀ cutoffs that were used to classify the cell lines were crucial for model performance and found that minor variations had little effect on the AUC (Supplementary Fig. S3). The threeClass approach performed less well in this analysis (Fig. 1E); however, this could be a result of low sample size when dividing the data into three groups. This therefore does not necessarily preclude an ability to identify repair defects. Consequently, the threeClass combinations were included in the subsequent validation studies.

For further external validation, we selected the best performing combinations of feature selection and classifier, defined as those within the top 25 of both the 3- and 5-fold cross-validations, resulting in 9 mmcOnly, 13 olaOnly, and 13 threeClass combinations (Fig. 1F; Supplementary Table S2). We used these to train models on the full HNSCC cell line expression dataset and combined the resulting models into ensembles (Fig. 1F). Scores of these ensembles that represent the probability of each cell line to fall into any of the sensitivity groups (normal, SensM, SensMO) are shown in Fig. 1G.

To verify that our models can be used to predict DNA repair defects, we next determined whether our models are predictive for sensitivity of DNA-damaging agents in other datasets. HNSCC cell line responses to a large number of compounds were determined as part of the Cancer Genome Project (CGP; ref. 26). The predictions from both MMC sensitivity relevant ensembles, mmcOnly, and threeClass (SensM and SensMO combined), correlate with MMC IC₅₀ in this dataset (Fig. 2A). This also translated into an ability to identify DNA repair–defective cell lines, as classified by their MMC sensitivity (Fig. 2B). We did not observe a correlation between the olaOnly ensemble scores and MMC sensitivity. The olaOnly model ensemble was consequently unable to distinguish MMC-sensitive cells. For olaparib sensitivity, we did not observe a correlation with model prediction scores for any ensemble, likely due to uncertainty in the computationally inferred IC₅₀ values in the CGP HNSCC cell lines (Supplementary Fig. S4A–S4C). Importantly, all three model ensembles also predicted the presence of mutations in DNA crosslink repair genes (Supplementary Fig. S4B; Supplementary Table S3), with the mmcOnly and threeClass model ensembles showing the best performance (Fig. 2C).

As a complementary validation strategy that these models predict DNA repair defects beyond drug sensitivity, we determined the capability of our models to detect alterations in DDR and repair. In a study by Peng and colleagues (27), components of the DDR and repair pathways were knocked down using shRNA in breast cancer cells and gene expression was measured using microarrays. We tested whether the model ensembles can distinguish cells treated with no or control shRNA from those treated with the tested shRNAs against DDR genes. We find that the ensembles accurately predicted DDR defects caused by the loss of ATM, ATR, CHK1 CHK2, and 53BP1 (Fig. 2D). The performance for BRCA1 and RAD51 knockdown was more heterogeneous (Fig. 2D), a characteristic that could be explained by the differences in downstream signaling and transcriptomics that are caused by the complete loss of the protein compared with deleterious gene mutations. Interestingly, loss of the DNA damage response mediator MCPH1 (BRIT1) results in a consistently and greatly reduced signal for DNA repair defects, suggesting an essential role of the chromatin remodeling cyclin response in establishing the distinctive gene expression patterns that are detected by these models.

On the whole, we conclude that both the mmcOnly and the threeClass model ensembles capture features of DDR and DNA crosslink repair defects, i.e., drug sensitivity and DDR gene mutations, and repair defects resulting from DDR gene knockdowns.

Using these two validated model ensembles, we determined whether survival is different in patients with tumors that show DNA repair defect–associated gene expression patterns. For this purpose, we used pretreatment tumor biopsy gene expression data from our in-house HNSCC cohort (NKI-CRAD) of 98 patients with advanced stage HPV-negative tumors (characteristics in Supplementary Table S4), who underwent radiotherapy with cisplatin-based treatment regimens. Consistent with previous reports, tumor site (31) and cumulative cisplatin dose (32) are prognostic factors in our NKI-CRAD cohort. DNA repair defect scores were established using the mmcOnly model ensemble for all tumor samples. These scores did not correlate with any of the potentially prognostic clinical factors (NKI-CRAD data in Supplementary Fig. S5 and Supplementary Table S5). We then divided the patient population into two equal-sized groups, based on the ranking of the mmcOnly ensemble scores, and found a trend toward poorer prognosis for the patients with the higher scores (Fig. 3A).

Hypoxia and cumulative cisplatin dose are independent prognostic factors in HNSCC (Supplementary Fig. S6A and S6B; ref. 33). To consider the hypoxia status, the average expression of five published hypoxia signatures (34) was calculated and used to divide the patient population into equally sized hypoxic and normoxic subgroups. Subsequent DNA repair defect classification by the mmcOnly models in these subgroups shows that the poor prognosis association is particularly strong in normoxic tumor samples (Fig. 3B). This interaction is statistically significant (P = 0.011) and was therefore included in subsequent multivariate analyses. Upon closer inspection, we find that overall gene expression changes from repair defects (high in DNA repair defect score) under normoxia correlate highly with those from hypoxia exposure (as determined by the hypoxia score) in repair-proficient tumors in both cohorts, even though the actual patient classification does not (Supplementary Fig. S7A–S7D). Pathway enrichment analyses revealed that the downstream effects of hypoxia and DNA repair defects are partially overlapping (Supplementary Fig. S7E–S7G). This may explain the lack of an endpoint association within the hypoxic tumor group and points to common repair defect sources or downstream signaling effects.

Tumors harboring DNA crosslink repair defects are expected to be sensitive to cisplatin treatment. Variations in the cumulative cisplatin dose in our patients, as a result of discontinuation due to systemic toxicities, gave us the opportunity to investigate this hypothesis. The patients were assigned to a high (>200 mg/m²) or low (≤200 mg/m²) cumulative cisplatin dose group according to previous reports, resulting in a univariate HR of 2.3 for the low-dose patients (P = 0.0014). Overall, patients with predicted DNA repair–defective tumors benefitted most from the high cumulative cisplatin dose (Fig. 3C; Supplementary Fig. S8A and S8B).

To consider both the prognostic clinical factors and the herein uncovered interaction with hypoxia, we combined predicted DNA repair status, hypoxia status, cisplatin dose, and tumor site in a multivariate model (Fig. 4A). Our analysis shows that expression profiles associated with DNA repair defects were strongly indicative of poor prognosis (HR = 4.01, P = 0.00091). The interaction between DNA repair and hypoxia status was also significant in this model (HR = 0.19, P = 0.0066), confirming that hypoxia and DNA repair defects are not additive (Fig. 4A, left). To independently test this poor patient prognosis association, we next computed DNA repair defect model scores for tumor samples from a second independent cohort with 82 patients with advanced and HPV-negative HNSCC who were chemoradiotherapy treated and from which pretreatment material was collected within the multicenter DESIGN project (characteristics in Supplementary Table S4). This independent validation, shown in Fig. 4A (right), confirms that predicted DNA repair defects, as determined by our model ensemble, are significantly associated with poor prognosis (HR = 3.2, P = 0.023).

By splitting the cohorts at their median scores into two equal-sized groups, as is common practice in biomarker studies, the above analyses stipulate that DNA repair defects are present in half of the HNSCC. Yet there is little support for this suggestion (50% DNA repair–defective tumors). Furthermore, prevalence can also differ between cohorts. To verify that our findings are not biased by certain cut-off points, we evaluated the multivariate HR at a range of cutoffs (Fig. 4B). The varied cutoffs for repair defect status assignment are based on the proportion of patients in one or the other category as defined by their model score ranking, resulting in different score cutoffs for both cohorts in this analysis. Doing so, we find that DNA repair defects called by the mmcOnly model ensemble were significantly associated with poor prognosis at a broad range of cutoffs in both the NKI-CRAD cohort and the second independent DESIGN cohort.

We next assessed the DNA repair defect/poor prognosis association by using the threeClass model ensemble. We find that this threeClass model ensemble showed highly similar results, confirming a DNA repair defect/poor prognosis association and that it also predicts poor survival in both advanced HNSCC patient cohorts largely independent of the chosen cutoffs (Fig. 4B; Supplementary Fig. S9A–S9D). As a comparison, clinical stage was not found to be significantly associated with poor survival in either cohort (Supplementary Fig. S10).

Finally, to further reveal the origin of the poor prognosis association and cisplatin responsiveness, we analyzed the association with different outcome endpoints. We find little evidence for a role in locoregional control (Fig. 5A). However, we observed significantly increased HRs for the occurrence of distant metastasis, in particular in the NKI-CRAD cohort (Fig. 5A), which suggests that DNA repair–defective tumors may be more prone to metastasize. As for overall survival (Fig. 3C), this effect was most pronounced in normoxic tumors that also metastasize less, unless repair defective (Fig. 5B). Given the lack of endpoint association in hypoxic tumors, we further investigated the effect of cumulative cisplatin dose in normoxic tumors only. HRs for DNA repair defects were higher in the patient group treated with low-dose cisplatin at all cutoffs (Fig. 5C and D). Although only statistically significant at one cutoff, likely due to low patient numbers, it points to the benefit of high-dose cisplatin to reduce metastasis occurrence in patients with DNA repair–defective tumors.

These findings prompted us to explore the link between DNA repair defects and metastasis in in vitro studies with the confirmed crosslink repair–impaired HNSCCs. We selected several representative cell lines for each sensitivity class (Normal, SensM, and SensMO; described in Fig. 1A). To investigate the link between the observed metastatic behavior and DNA repair defects, we used two assays to evaluate migratory behavior in these cell lines. In scratch assays, a confluent layer of cells is scratched and the rate at which the remaining cells close the resulting gap is determined (35). In transwell assays, the ability of cells to migrate through a porous membrane is measured (35). We found that cell lines that show MMC sensitivity (SensM and SensMO) migrated significantly more often and faster than Normal classified cell lines in both, the scratch assays (Fig. 6A and B) and transwell migration assays (Fig. 6C). To evaluate invasive capacity in vitro, both assays were modified by depositing extracellular matrix proteins (Matrigel) in the scratch or on the membrane. We observed that in the scratch and transwell invasion assays, the SensMO cell lines were able to invade through Matrigel, while this was not evident in the SensM or Normal cell lines (Fig. 6D–F). To further test whether this is DNA repair defect related, we inhibited HR/FA DNA repair in the repair-proficient Normal cell lines using the RAD51 inhibitor B02 for 24 hours. This resulted in strong induction of migratory (P = 1.49 × 10⁻⁵) and invasive behavior (P = 1.8 × 10⁻⁶ in ANOVA; Fig. 6G and H), further supporting a possible direct link between DNA repair defects and invasiveness in HNSCC. Interestingly, the quick induction also indicates a lack of a requirement for genetic adaptation. Treatment with equal concentrations of RAD51 inhibitor was lethal for both SensM and SensMO cell lines (Supplementary Fig. S11), precluding further epistasis experiments. These in vitro data support the notion derived from the clinical data that functional FA/HR repair pathway defects are associated with a more aggressive, metastasis-prone nature.

In search for potential mechanistic links and prompted by the high specificity in identifying DNA damage response deregulation by ATM knockdown (Fig. 2D; Supplementary Fig. S7), we hypothesized that the activation of DNA damage response pathways induces migration. This is in line with our initial position that unrepaired damage from oxidative and endogenous sources would provoke a DNA damage response and accompanying transcriptional changes that allows us to identify repair defects. Indeed, ATM depletion was previously reported to reduce migration and invasion. This was shown to occur through IL8 signaling regulated by ATM activation and induced by oxidative stress (36). We find that CXCL8 (IL8) is consistently upregulated in the migrating but not in the invading cell lines (Fig. 7A). Inhibition of the DNA damage response protein ATM reduced migration mildly in some cell lines, suggesting that this is not specific to a DNA repair–defective status (Fig. 7B; Supplementary Fig. S12). A second key player in the DNA damage response network is the DNA-dependent protein kinase, DNA-PK. DNA-PK has been shown to promote migration and invasion in prostate cancer, possibly through DNA-PK–mediated gene expression changes (37). Others postulated DNA-PK influences invasion and migration through the control of the secretion of multiple metastasis-associated proteins (38). When tested in our cell lines, DNA-PK inhibition reduced migration in the BRCA1-mutated UT-SCC-60B (Fig. 7C; Supplementary Fig. S12A–S12D). Together, our data suggest that DNA-PK activation and IL8 expression–mediated processes may be responsible for the migration of cells within the migratory SensM group.

As the invasive SensMO group did not show such mechanistic links, we next evaluated the presence of known metastasis-promoting features. Metastasis is partly associated with epithelial-to-mesenchymal transition (EMT; ref. 39). DNA damage, such as induced by radiation, has been reported to promote EMT and migration in cell lines. We therefore tested whether DNA repair–defective cell lines show mesenchymal characteristics, elongated cell body, and expression of key transcription factors (40). Inspection by light microscopy revealed that both SensM and SensMO cell lines are more spindle-shaped than Normal cell lines, with SensMO cells showing the most elongated phenotype (P < 10⁻¹⁶ for both; Fig. 7D; Supplementary Fig. S13A). However, while some migratory and invasive cell lines exhibited mesenchymal characteristics, we did not find consistent differences between the cell line classes for expression of the classical EMT markers or EMT-related signatures (41) when using all cell lines. Assuming repair-proficient cell lines may also exhibit EMT (unrelated or not caused by DNA repair defects), we performed gene set enrichment analysis using only the cell lines for which migration and invasion data were available. EMT-related genes were significantly enriched (Fig. 7E), pointing to a possible mechanistic role for EMT. When testing EMT-related processes (41, 42), we found cell adhesion and migration genes to be enriched (Supplementary Fig. S13B) and most migrating cells to express mesenchymal markers (Supplementary Fig. S13C). EMT and cellular invasive behavior as mediated by TGFβ requires SMAD4 (43, 44). SMAD4 loss, that can cause spontaneous head and neck cancer in mice, however, promotes squamous cell carcinoma metastasis and was shown to influence metastasis formation in cooperation with E-cadherin in other models (45, 46). Using DNA capture sequencing analysis (16), we find that UT-SCC-24B lacks SMAD4 expression due to a genomic deletion of the first three exons (Fig. 7F and G). This finding links HR/FA repair defects with SMAD4 and invasive behavior in one of the cell lines.

In conclusion, while we find evidence that repair defects can be linked to migratory behavior through the engagement of key DNA damage response players and their links to migration as reported by others, mechanisms in the invasive repair-defective class may depend on genetic context and may be more multifaceted and complex.

In this study, we set out to investigate the effect of DNA repair defects on clinical outcome in HNSCC. We defined DNA repair defects by functional outcomes, hypersensitivity to DNA-damaging agents, which integrates the proficiency of all DNA damage response and DNA repair pathways, including adaptations that have arisen to any defects in these pathways. We reasoned that the transcriptome reflects the sum of all failed or successful repair activities on endogenously generated damage and therefore encompasses defects and the activation of backup repair pathways and other coping mechanisms. Reestablishment of DNA repair, through compensatory repair activities, is a mechanism for acquiring treatment resistance (9, 47), which does not revert genomic scars or mutations, but is presumably reflected in the transcriptome. Together, we consider this an important addition to the current arsenal of DNA repair biomarkers.

Our strategy allowed us to investigate the effect of DNA repair defects in patients with HNSCC. This revealed that patients with tumors predicted to harbor DNA repair defects showed worse overall survival in two independent patient cohorts. We conclude that they constitute a HNSCC group with a more aggressive behavior. Hypoxia is known to alter cellular repair capacity and pathway choice and is a poor prognosis biomarker (48). We find a noticeable interaction with hypoxia that suggests that, in patients with hypoxic tumors, prognosis is not further decreased by repair defects. Hypoxia has been shown to promote metastasis and to downregulate HR/FA repair processes (49–51). A considerable overlap in gene expression changes to the tumors with DNA repair defects indicated common repair defect sources or downstream signaling effects. Alternatively, the models may be unable to predict repair defects in hypoxic tumors, because the model ensembles were generated under nonhypoxic conditions. The RNA expression patterns in the DNA repair defective cells were likely caused (in part) by oxidative stress. In contrast to our findings, in breast, ovarian, and lung cancer, DNA repair defects in tumors were found to be associated with a favorable prognosis, possibly due to the treatment with DNA-damaging agents (52–55). Current treatments for HNSCC may hence not fully exploit the cellular sensitivities of tumors with DNA crosslink repair defects. The apparent increased benefit of high cumulative cisplatin dose in our clinical data supports this interpretation. Given the risk, such patients could be candidates for additional systemic treatment with crosslinkers, PARP inhibitors, or other novel therapies targeting such defects.

The in vitro migratory and invasive capacity that was exclusive to the SensMO cell lines confirmed a proposed propensity of tumors with DNA repair defects to metastasize. In addition, we found that the induction of DNA repair defects, in cell lines that were otherwise repair proficient, results in rapid acquisition of this invasive phenotype. This suggests that a direct link and signaling mechanism underlies this effect, rather than genetic contextual adaptations to such repair defects. In line with these findings, a recent study showed that chromosomal instability induces invasive behavior through activation of NFκB signaling (25). Others have reported that knockdown of FA genes in HNSCC cell lines results in increased invasion through activation of the NHEJ repair pathway and via DNA-PK and Rac1 signaling (56). Our DNA-PK inhibition data and differential pathway expression analyses support the engagement of these migration-promoting mechanisms in DNA repair–deficient HNSCC cell lines and provides a mechanistic link to a higher propensity to metastasize. DNA-PK activation–mediated migration is also in line with an increased engagement of NHEJ as indicated by the PARP inhibitor resistance and connects repair defects to the reports from others that show a role for these proteins in migration, invasion, and metastasis. Consistent with the notion that loss of FA/HR DNA repair may induce metastatic behavior, knockdown of BRCA1 and FANCD2 was reported to cause dedifferentiation and induce EMT in human mammary epithelial cells in vitro (56). Cells that have undergone EMT downregulate apoptosis (57) and show increased resistance to DNA damage (58). At first glance, our data seem inconsistent with such observations. Upon closer inspection, we find hints of an increased mesenchymal phenotype in the gene expression data and cellular morphology of the invasive cell lines. In recent years, however, it has become clear that cancer cells tend not to be in a solely epithelial or mesenchymal state, but that many intermediate states exists (59). Also, DNA repair–defective cell lines may show a partial EMT phenotype, this as part of a defense mechanism against further DNA damage and the effects thereof.

In the cell line panel, we observed that invasive behavior and DNA repair defects overlap to a large degree. Considering the recent data from Bakhoum and colleagues (25), our data further support a causal relationship between defective DNA repair that promotes chromosomal instability and migratory behavior. However, given this overlap in features, in patient material and based on our transcriptome analyses alone, a clear distinction between DNA repair deficiency and the sole and independent presentation of invasive features is therefore formally not possible. The association with cisplatin response and the inability to predict survival changes in hypoxic samples or to identify all metastasis-prone cases points, however, to repair-defect features in contrast to pure metastasis features. For more accurate identification of tumor DNA repair defects in individual patients, integration of transcriptomics, genomic rearrangements, pathogenic mutations, and phospho-proteomics is desirable. More cohorts of homogenously treated HSNCC patients, for whom full transcriptomics, genomic data, and carefully recorded clinical data are available, are required to further optimize the identification of patients with tumors harboring DNA repair defects with the purpose to provide alternative treatment options to this poor prognosis group.

In summary, we find that DNA repair–defective cell lines are more motile and invasive in vitro and we have developed gene expression–based models to detect DNA repair defects in HNSCC. Importantly, using these models, we find that these defects are associated with poor overall survival in two independent patient cohorts thereby confirming their prognostic value and the detrimental role of tumor repair defects.

No potential conflicts of interest were disclosed.

Conception and design: P.B.M. Essers, M.W.M. van den Brekel, H. Bartelink, C. Vens

Development of methodology: P.B.M. Essers, C.V.M. Verhagen, M.W.M. van den Brekel, C. Vens

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P.B.M. Essers, M. van der Heijden, C.V.M. Verhagen, E.M. Ploeg, R.H. de Roest, C.R. Leemans, R.H. Brakenhoff, H. Bartelink, M. Verheij, C. Vens

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P.B.M. Essers, M. van der Heijden, C.V.M. Verhagen, C.R. Leemans, R.H. Brakenhoff, M.W.M. van den Brekel, C. Vens

Writing, review, and/or revision of the manuscript: P.B.M. Essers, M. van der Heijden, C.V.M. Verhagen, R.H. de Roest, C.R. Leemans, R.H. Brakenhoff, M.W.M. van den Brekel, H. Bartelink, M. Verheij, C. Vens

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C.V.M. Verhagen, E.M. Ploeg, C. Vens

Study supervision: P.B.M. Essers, M.W.M. van den Brekel, M. Verheij, C. Vens

This research was funded by the EU 7th framework program (257144 ARTFORCE), the Dutch Cancer Society (KWF-A6C7072, DESIGN consortium, to R.H. Brakenhoff, C. Vens, and M. van den Brekel) and the Brunel and Verwelius funds (to C. Vens and M. van den Brekel). We would like to thank the NKI Genomics Core Facility for performing RNA sequencing, the NKI RHPC facility for providing computational resources, and the Core Facility-Molecular Pathology and Biobank (CFMPB) for collecting and preparing tissue samples. We are thankful for B. Floot's help and the kind gift of UT-SCC lines from Professor Reidar Grénman (Turku University Hospital, Turku, Finland). We thank Robert W. Sobol (University South Alabama, Mobile, AL) for reviewing and providing editorial comments.

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.