Ovarian cancer-derived copy number alterations signatures are prognostic in chemoradiotherapy-treated head and neck squamous cell carcinoma.

DNA copy number alterations (CNAs) are frequent in cancer, and recently developed CNA signatures revealed their value in molecular tumor stratification for patient prognosis and platinum resistance prediction in ovarian cancer. Head and neck squamous cell carcinoma (HNSCC) is also characterized by high CNAs. In this study, we determined CNA in 173 human papilloma virus-negative HNSCC from a Dutch multicenter cohort by low-coverage whole genome sequencing and tested the prognostic value of seven cancer-derived CNA signatures for these cisplatin- and radiotherapy-treated patients. We find that a high CNA signature 1 (s1) score is associated with low values for all other signatures and better patient outcomes in the Dutch cohorts and The Cancer Genome Atlas HNSCC data set. High s5 and s7 scores are associated with increased distant metastasis rates and high s6 scores with poor overall survival. High cumulative cisplatin doses result in improved outcomes in chemoradiotherapy-treated HNSCC patients. Here we find that tumors high in s1 or low in s6 are most responsive to a change in cisplatin dose. High s5 values, however, significantly increase the risk for metastasis in patients with low cumulative cisplatin doses. Together this suggests that the processes causing these CNA signatures affect cisplatin response in HNSCC. In conclusion, CNA signatures derived from a different cancer type were prognostic and associated with cisplatin response in HNSCC, suggesting they represent underlying molecular processes that define patient outcome.

copy number alterations

Cox proportional hazards

distant metastasis

head and neck squamous cell carcinoma

human papilloma virus

hazard ratio

ovarian cancer

overall survival

progression‐free survival

CNA signature 1–7

squamous cell carcinoma

The Cancer Genome Atlas

laryngeal and pharyngeal tumors of the TCGA HNSCC cohort

oral cavity tumors of the TCGA HNSCC cohort

HPV‐positive HNSCC of the TCGA HNSCC cohort

whole genome sequencing

Introduction

DNA damage and aberrant or failed DNA repair promote the mutational processes that drive DNA copy number alterations (CNAs) and shape CNA patterns in human cancer. Thus, similar to mutational signatures,1 CNAs reflect the involvement of crucial mutational processes and have the potential to assist molecular stratification of tumors for precision medicine. Studying high‐grade serous ovarian cancer (OC), MacIntyre et al.2 derived and validated seven CNA signatures to decode the complexity of the genomic changes embedded in the copy number change patterns. Importantly, these CNA signatures were found to be associated with patient outcome2 and platinum resistance relapse, thereby showing both, the relevance of such mutational processes in treatment outcome and the potential as biomarker. However, the application and value in other cancer types remain elusive.

Among all cancer types, head and neck squamous cell carcinoma (HNSCC) is, similar to OCs and BRCA1/2‐mutated breast cancers, predominantly characterized by high CNAs.3, 4, 5 Although certain frequently found gene amplifications and losses are involved in the progression of healthy tissue to cancerous lesions,6 only few of these particular CNAs have a prognostic value in established tumors.7, 8 Like ovarian carcinoma, HNSCC exhibits frequent alterations in DNA repair pathways that govern genomic stability8, 9, 10 and respond to cisplatin combination treatment in the clinic.11 Using novel gene expression models that identify DNA crosslink repair defects, our previous studies revealed the significance of such repair defects in HNSCC and a link to cisplatin‐based chemoradiotherapy treatment response in advanced stage human papilloma virus (HPV)‐negative HNSCC.12

In this study we set out to investigate whether the OC‐derived CNA signatures are discernible and prognostic in HNSCC.

Material and Methods

Patient data and material

Primary pretreatment tumor material from 173 patients with hypopharyngeal, laryngeal or HPV‐negative oropharyngeal HNSCC treated at the Netherlands Cancer Institute (NKI) and at the Amsterdam Universitair Medische Centra (UMC) location VUmc, between 2001 and 2014 (Supporting Information Fig. S1 and Table S1) was available for this retrospective study. Fresh‐frozen tumor material from pretreatment biopsies was analyzed, and the study was approved by the Institutional Review Boards of the two institutions, NKI and VUmc. All patients signed informed consent for the use of material and data. Samples with tumor percentage <40%, as determined by hematoxylin and eosin staining, were excluded. HPV‐positive patients were excluded using immunohistochemistry for p16 and p53, subsequent PCR and RNA‐sequencing.12 All patients received cisplatin‐based chemoradiotherapy with a total dose of 70 Gy that followed the conventional (35 fractions of 2Gy over 7 weeks) or the Danish Head and Neck Cancer Group scheme (35 fractions of 2 Gy over 6 weeks). Cisplatin treatment in this cohort encompassed daily (25 × 6 mg/m2 body surface area), weekly (6 or 7 × 40 mg/m2) or 3‐weekly (3 × 100 mg/m2) intravenous administration. Eight patients received intraarterial cisplatin (4 × 150 mg/m2) according to the RADPLAT trial.13 The total cumulative cisplatin dose received was available for 170 patients and used for outcome association analyses. Overall survival (OS) 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.

Low‐coverage whole genome sequencing and copy numbers calculations

DNA was extracted using the AllPrep DNA/RNA kit (QiaGen, Hilden, Germany) and libraries were constructed using the KAPA HTP library preparation kit (Ref no. 7961901001, KAPA/Roche, Wilmington, DL) according to manufacturers' instructions and sequenced on an Illumina HiSeq 2500 to a depth of approximately 0.5×. Reads were mapped to the hg19 genome using bwa‐0.7.17.14 Relative copy numbers were obtained for predefined 10‐kb bins using QDNASeq.15 Absolute copy numbers were calculated using the squaremodel function in the ACE BioConductor package,16 with the settings “penalty = 0.5” and “penploidy = 1.5.” Of the proposed models, the best fit that did not imply a 100% pure tumor was selected. The MacIntyre data set was used to benchmark absolute copy number estimation (ACE) in Poell et al.,16 showing a high concordance of ACE‐derived absolute copy numbers with those derived from deep sequencing data.16 ACE performance in our data set was further evaluated using exon capture sequencing data for 556 genes available for a subset of 54 HNSCC patients, as described in Verhagen et al.8 Overall tumor purity and gene ploidy estimates for those samples were calculated using PureCN.17 ACE and PureCN tumor purity estimates were highly correlated (R 2 = 0.87; Supporting Information Figs. S2). Inspection of ploidy estimates for individual genes confirmed that ACE estimates are largely in agreement with PureCN estimates, indicating accurate estimation of absolute copy numbers (Supporting Information Fig. S2).

CNA signatures were calculated using the R scripts available from MacIntyre et al.2 Low numbers of breakpoints and large segment sizes characterize signature 1 (s1) that oppose the high numbers with small segment sizes of s2. Copy number changes to 3 or single copies and high segment copy numbers distinguish s2, s3 and s4 respectively. However, s5 captures the presence of subclonal copy numbers reductions and s6 is characterized by large and pronounced copy number changes and gains. Common losses in tetraploid genomes are the most characteristic feature for s7.

The Cancer Genome Atlas data sets

Copy number data for The Cancer Genome Atlas (TCGA) HNSC cohort (n = 522) were downloaded from https://portal.gdc.cancer.gov/ (on June 19, 2019). Surgical resection was the primary treatment for all patients in this cohort. HPV‐negative laryngeal and hypopharyngeal or oropharyngeal cases were selected for the laryngeal and pharyngeal tumors of the TCGA‐HNSCC (TCGA‐HNSC‐LP) subcohort (n = 128) analyses and compared to the oral cavity cases in the oral cavity tumors of the TCGA‐HNSCC (TCGA‐HNSC‐OrC) cohort of this study (n = 293), to HPV‐positive HNSCC of the TCGA HNSCC (TCGA‐HNSC‐HPV+) cohort (n = 99) or the oropharyngeal TCGA‐HNSCC (TCGA‐OPSCC‐HPV+) (n = 40) with all HPV‐positive HNSCC or solely oropharyngeal, respectively. Clinical data containing OS and progression‐free interval was obtained from Liu et al.18 Absolute copy numbers and CNA signatures were calculated as described above.

Statistics

All statistics were carried out in the R environment for statistical computing19 as specified.

Data availability

The Dutch multicenter cohort CNA and RNA‐Seq data that support the findings of this study are available in European genome‐phenome archive (EGA) at https://ega-archive.com under the EGA study numbers EGAS00001004090 (RNA‐Seq data), EGAS00001004091 (low‐coverage whole genome sequencing [WGS]) and data set numbers EGAD00001005716, EGAD00001005715 and EGAD00001005719, EGAD00001005718 for the RNA‐Seq and low‐coverage WGS, respectively.

Results and Discussion

CNV signatures in HNSCC

We performed low‐coverage WGS on pretreatment tumor biopsy material from 173 patients from two previously described Dutch multicenter cohorts (NKI‐CRAD/DESIGN).12 All patients had advanced pharyngeal and laryngeal squamous cell carcinoma (SCC) and were treated with cisplatin‐based chemoradiotherapy (Supporting Information Fig. S1 and Table S1). Given the distinct biology of HPV‐positive oropharyngeal cancers that are associated with a good prognosis, these were excluded.20

The distributions of the six features that, in different linear combinations, compose the seven MacIntyre copy number signatures2 were calculated: segment length, segment copy number, number of breakpoints per chromosome arm, number of breakpoints per 10 kb, length of oscillating copy number segment chains and the difference in copy number between adjacent segments in the two available HNSCC data sets (NKI‐CRAD/DESIGN and TCGA‐HNSC). The distributions of these features were similar to those in the BriTROC ovarian cancer data set used to construct the signatures (Supporting Information Fig. S3). Some notable differences were lower numbers of breakpoints per chromosome arm and larger numbers of deletions in the HNSCC cohorts as compared to the ovarian cancer cohorts (Supporting Information Fig. S3).

The seven signatures were calculated from the distributions of these features for all patient samples in our cohorts. No clear differences in signature scores between anatomical subsites (Fig. 1 a) were apparent. HNSCC and OC are characterized by high average s1 scores, which are not observed in carcinoma from other sites (Supporting Information Fig. S4), a result of distinctly different CNA feature distributions (Supporting Information Fig. S5). On the other hand, s1 high tumors are more frequent in the HPV‐positive oropharyngeal SCC (Supporting Information Fig. S4). Apart from gradually decreasing s7 scores with increasing T‐stage, no strong relations were found between signature scores and other clinical variables (Supporting Information Fig. S6). Overall copy number signature score distributions were similar between the cohorts, with the TCGA‐HNSC dataset displaying more samples with a high s5 score and the NKI‐CRAD/DESIGN cohort more s1 and s6 high samples (Fig. 1 b). Tumors scoring high for s1 typically scored low for all other signatures (Fig. 1 c). Hierarchical clustering confirmed this pattern in both HNSCC cohorts as it resulted in a primary Cluster 1, with a strong “s1‐high/others‐low” CNA signature characteristic (Fig. 1 c). Principal component analysis confirms this pattern (Supporting Information Figs. S7). Analyzing patient survival association for this CNA signature pattern‐characterized HNSCC subgroup, we find that the NKI‐CRAD/DESIGN patients within the “s1‐high” cluster had a significantly reduced risk of distant metastases (Fig. 1 d), although this did not lead to a significant increase in OS (hazard ratio [HR] = 0.76; p = 0.23). A similar distribution was observed in the TCGA‐HNSC‐LP data set patients (Figs. 1 c and 1 e and Supporting Information Figs. S8 for all TCGA‐HNSC) confirming a trend toward better prognosis for the “s1‐high” cluster patients.

Figure 1

Distribution of CNA signatures in HNSCC. (a) Mean signature score proportions across the NKI‐CRAD/DESIGN cohort, in the full cohort and split by anatomical subsite. (b) Comparison of signature scores across HNSCC cohorts used in this study and the ovarian cancer (OC) cohorts used for signature generation and validation in MacIntyre et al.2 Note that signature scores add up to 1 in each patient. (c) Hierarchical clustering of the NKI‐CRAD/DESIGN cohort patients (upper panel) and the TCGA‐HNSC‐LP, that is, restricted to pharyngeal and laryngeal, by the seven CNA signature scores. (d) DM‐free survival for the “s1‐high” cluster NKI‐CRAD/DESIGN patients (labeled in red) from c compared to the rest with multivariable CoxPH ratio values including tumor site and stage and associated p‐value. (e) OS for the “s1‐high” cluster TCGA‐HNSC‐LP patients from c compared to the other two primary clusters patients.

Overall these data show that HNSCC CNA signatures patterns could be applied in molecular stratification of tumors and suggest that these signatures have potential as prognostic markers in HNSCC.

Prognostic value of signature scores

To further investigate the prognostic potential, we tested the individual CNA signatures in the NKI‐CRAD/DESIGN cohort. Cox proportional hazards (CoxPH) models for all available survival endpoints were constructed for each signature separately, which included the signature score, and the known prognostic markers, tumor subsite and tumor stage, as variables (Fig. 2 a). Confirming the cluster‐derived involvement of s1, we find a trend toward an association with decreased risk for all outcomes, most notably distant metastasis (DM) for s1. In contrast, s6 showed a trend toward an association with increased risk of disease progression. Importantly, s5 and s7 were significantly associated with increased DM risk (Fig. 2 a). As HRs for continuous scores can be difficult to interpret, Kaplan–Meier plots are shown in Figure 2 b illustrating these associations.

Figure 2

Prognostic CNA signatures in HNSCC. (a) Forest plot showing the results from CoxPH model fits for OS, PFS, locoregional control and DM‐free survival. Signature scores were scaled to an average of 0 and standard deviation of 1, so that HRs shown represent a 1 standard deviation increase of each score. Each signature was then individually assessed as a variable in a CoxPH model containing anatomical subsite and disease stage as additional variables. Only the estimates for the signatures are shown and compared to each other in this forest plot. (b) Kaplan–Meier plots illustrating the main findings in a and based on a median cohort split. Hence, s5 or s7 scores were 0 in many patients and these are therefore shown as a separate group. The remaining patients were divided into equal sized groups. (c) Illustration of cohort division cutoff dependence and prognostic strength. Patients were split into a low and high group at each possible cutoff for each signature as indicated by the color code in the legend and a CoxPH model was fit with the same variables as in a. At a cutoff of 0.2, 20% of patients are considered “low” and 80% are considered “high.” At each cutoff, the resulting HR for each signature is depicted. HRs associated with a p < 0.05 are indicated with dots. Scripts for these HR plots are available from https://github.com/PaulEssers/SurvivalPlots.

The use of continuous variables in a CoxPH model assumes a linear relationship between signature scores and survival. To circumvent this assumption and avoid choosing an arbitrary cutoff, patients were repeatedly split into a high and low score group at each possible cutoff and for each signature. The resulting multivariable CoxPH‐derived HRs are depicted in Figure 2 c. This confirms the findings shown in Figure 2 a: high s5 and s7 are associated with increased DM risk, regardless of the exact cutoff. In contrast, s6 is associated with OS and progression‐free survival (PFS) at a few cutoff points only. However, s1 outcome associations behave opposite to those with s5, s6 and s7. As noted above, in general, high scores of s1 correlate with low scores of the other signatures. In contrast to HNSCC, s1 and s2 were associated with poor prognosis in the OC study, and s3 and s7 with good prognosis.2 Tissue origin, oncogenomic context, patient treatment or the authors' methodological approach to calculate each signature relative to s5 may have influenced the prognostic role of the CNA signatures.

To confirm associations with prognosis, we used the independent TCGA‐HNSC cohort, which, however, consists of surgically treated cases only, with a plethora of secondary treatments. Nonetheless, the progression prognosis association for s1 and s6 could be confirmed by the TCGA‐HNSC‐LP data set (Supporting Information Figs. S8). Due to the lack of DM data, we were not able to test a metastasis association for s5 and s7 in resected HNSCC. However, we do observe a trend toward poor prognosis for s7‐high TCGA‐HNSC‐LP patients in a different clinical endpoint. In addition, we find a role for the CNA s4 in TCGA‐HNSC patients (Supporting Information Figs. S8), which is not present in the NKI‐CRAD/DESIGN cohort but could result from the different treatment. Interestingly, oral cavity HNSCC does not show an outcome association for s6 (Supporting Information Figs. S8) further suggesting tissue and oncogenetic context specificity.21, 22 Together, we observe a consistent prognosis association for s1, s6 and s7 in both cohorts. The trend toward better prognosis for s1‐high or poorer prognosis after resection when high in s4 was also visible in the HPV‐positive TCGA‐HNSC‐HPV+ cohort (Supporting Information Fig. S8). However, s1 and s7 high tumors are rare in esophageal carcinoma (ESCA) and lung SCC (Supporting Information Fig. S4), yet a poor outcome association is apparent for s6‐high in the TCGA‐ESCA cohort (Supporting Information Fig. S9).

Association with cisplatin treatment response

CNA signatures predicted platinum‐resistant relapse in the OC study of MacIntyre et al.2 High cumulative cisplatin doses of at least 200 mg/m2 were previously reported to improve treatment outcome in chemoradiotherapy‐treated HNSCC23 independent of administration mode (Supporting Information Fig. S10) and result in a HR of 0.57 (p = 0.0053) in the NKI‐CRAD/DESIGN cohort.12 We therefore investigated the CNA signature outcome associations within patient cohorts with different cumulative cisplatin doses. Performing multivariable CoxPH as in Figure 2, we find that the outcome association of s1 and s6 are exclusive to the high‐dose cisplatin group (Fig. 3 a) showing a greater benefit from high cumulative cisplatin doses in the s1‐high and s6‐low patient groups (Fig. 3 b). Similarly, high s5 scores were only prognostic in the low‐dose cisplatin group (Figs. 3) further revealing that the mutational processes linking s5 with a risk of DM in HNSCC (Fig. 2) may influence responsiveness to cisplatin.

Figure 3

Cumulative cisplatin dose dependence of CNA signature association with patient outcomes. (a) Forest plot showing the results from multivariable CoxPH model fits that include tumor site and stage, as in Figure 2a. Separate models were fit in the patient subpopulations that did or did not receive at least 200 mg/m2 cisplatin. (b) Kaplan–Meijer plots illustrating the main findings in a. Patients were divided as in Figure 2b and additionally divided by the cumulative cisplatin dose category received.

Links to candidate pathways

To elucidate links to mutational processes, we next investigated whether the individual signatures corresponded to alterations in the candidate pathways proposed in MacIntyre et al.2 In line with the OC report, we find that high s1 scores are strongly linked to the presence of HRAS hotspot mutations in the TCGA HNSCC data set (Supporting Information Fig. S11). However, as previously reported21 and also evident in Dogan et al.,24 HRAS hotspot mutations were exclusive to oral cavity tumors, explaining their absence from the NKI‐CRAD/DESIGN cohort (Supporting Information Fig. S11) but suggesting other origins for the s1‐type CNA in HNSCC. Further analysis revealed that s1 is strongly anticorrelated to the total number of breaks (Supporting Information Figs. S12) that associate with poor outcomes (Supporting Information Fig. S12), suggesting that high s1 scoring patients represent a group with low numbers of copy number alterations. Interestingly, previous reports linked HRAS mutations to a CNA silent group in HNSCC.3, 20 For the prognostic s6, we did observe a correlation with average segment copy numbers (Supporting Information Fig. S12), which suggests that high copy numbers may drive poor prognosis.

We next turned our attention to pathways with a known relevance for patient outcomes in HNSCC patients treated with chemoradiotherapy.25 Prognostic expression markers and models reflecting these biologic factors were analyzed for correlations to the signatures. DNA crosslink repair defects and cellular proliferation were found to be linked to several CNA signatures (Supporting Information Fig. S11). OCs with mutations in BRCA1/2 and homologous recombination genes have higher s3 scores.2 In line with this observation, s3 scores correlate with an expression based marker for DNA crosslink repair defects12 in HNSCC (Supporting Information Figs. S11). Mutation enrichment (at the gene or hallmark level) or mutational signature analyses did not reveal any robust associations with the prognostic s1 cluster. Notably, increased s1 scores, that is, decreased number of breaks, are linked to reduced proliferation gene expression signature values in HNSCC consistent with a role of proliferation in break induction and propagation (Supporting Information Figs. S11). Also consistent with its good prognosis association, we find that s1 cluster tumors have a trend towards lower stem cell marker SLC3A and acute hypoxia scores in both cohorts (Supporting Information Figs. S11).

The observed outcome associations further highlight the role of the processes underlining the formation of these CNA patterns. Together these data point to molecular tumor stratification opportunities based on CNA signatures that may support multiparametric prediction models for HNSCC outcomes to further improve treatment evaluation and choice.

Conclusions

In this study we apply ovarian cancer‐derived CNA signatures in HPV‐negative HNSCC and find a prognostic role for these signatures. These findings suggest that CNA patterns, reflecting molecular processes leading to genomic rearrangements, can characterize HNSCC with different outcomes. We find a similar overall prognostic role in patient outcome for s1, s6 and s7 in both HNSCC cohorts. The lack of a role for the cisplatin responsive s5 characteristics in the surgically treated TCGA‐HNSC patient cohort and a different role for s4 points to the potential of these signatures to predict treatment response that warrants future studies in additional cohorts.

Conflict of interest

The authors have no conflict of interest to declare.

Table S1 Patient cohort characteristics.

Figure S1. Patient selection flowchart.

Figure S2. Validation of ACE estimation of absolute copy numbers.

Figure S3. Copy number features in HNSCC compared to the ovarian cancer cohorts.

Figure S4. Comparison of CNA signature distribution in different cancers.

Figure S5. Copy number features in HNSCC compared to lung and esophageal cancer.

Figure S6. Relation between signature scores and clinical factors.

Figure S7. Principal Component Analysis (PCA) in both HNSCC cohorts.

Figure S8. Prognostic value of signatures in the TCGA‐HNSC cohort.

Figure S9. Prognostic value of signatures in other TCGA cohorts.

Figure S10. Influence of cumulative cisplatin dose in chemoradiotherapy treated advanced HNSCC.

Figure S11. CNA signatures association with candidate genetic alterations or prognostic factors.

Figure S12. Relation of signatures to average copy number and number of breaks.

Click here for additional data file.