Abstract
Young-Gyu Eun1,2, Dongjin Lee1,3, Young Chan Lee2, Bo Hwa Sohn1, Eui Hyun Kim1,4, Sun Young Yim1,5, Kee Hwan Kwon3 and Ju-Seog Lee1
1Department of Systems Biology, The University of Texas MD Anderson Cancer Center, Houston, Texas, USA
2Department of Otolaryngology-Head and Neck Surgery, School of Medicine, Kyung Hee University, Seoul, Republic of Korea
3Department of Otolaryngology-Head and Neck Surgery, School of Medicine, Hallym University, Seoul, Republic of Korea
4Department of Neurosurgery, Severance Hospital, Brain Tumor Center, Yonsei University College of Medicine, Seoul, Republic of Korea
5Division of Gastroenterology and Hepatology, Department of Internal Medicine, Korea University College of Medicine, Seoul, Republic of Korea
Correspondence to:
Ju-Seog Lee, email: [email protected]
Keywords: YAP1; gene signature; head and neck cancer; prognosis
Received: February 27, 2017 Accepted: July 18, 2017 Published: November 27, 2017
ABSTRACT
By analyzing the genomic data of head and neck squamous cell cancer (HNSCC), we investigated clinical significance of YAP1 activation. Copy number and mRNA expression of YAP1 were analyzed together to assess clinical relevance of YAP1 activation in HNSCC. The clinical significance of YAP1 activation was further validated in four independent test cohorts. We also assessed the correlation of YAP1 activation with genomic alterations such as copy number alteration, somatic mutation, and miRNA expression. The YAP1-activated (YA) subgroup showed worse prognosis for HNSCC as tested and validated in five cohorts. In a multivariate risk analysis, the YAP1 signature was the most significant predictor of overall survival. The YAP1-inactivated (YI) subgroup was associated with HPV-positive status. In multiplatform analysis, YA tumors had gain of EGFR and SNAI2; loss of tumor-suppressor genes such as CSMD1, CDKN2A, NOTCH1, and SMAD4; and high mutation rates of TP53 and CDKN2A. YI tumors were characterized by gain of PIK3CA, SOX2, and TP63; deletion of 11q23.1; and high mutation rates of NFE2L2, PTEN, SYNE1, and NSD1. YA tumors also showed weaker immune activity as reflected in low IFNG composite scores and YAP1 activity is negatively associated with potential response to treatment of pembrolizumab. In conclusion, activation of YAP1 is associated with worse prognosis of patients with HNSCC and potential resistance to immunotherapy.
INTRODUCTION
Head and neck squamous cell carcinoma (HNSCC) is the sixth leading type of cancer worldwide, with an annual incidence of approximately 600,000 cases and a mortality rate of 40% to 50% [1, 2]. The major known risk factors are environmental exposures to tobacco products, alcohol, and infection with high-risk human papillomaviruses (HPV) [3]. Despite advances in our knowledge of the epidemiology and pathogenesis of HNSCC and the use of radiation treatment and functional surgery for this disease, survival rates have not improved over the past 40 years [4].
The Hippo pathway is a major tumor-suppressor pathway in many cancers and, in general, has important regulatory functions in cell proliferation, cell survival, cell competition, and maintenance of a stem cell phenotype [5]. The Hippo pathway consists of a regulatory serine-threonine kinase module and a transcriptional module. The kinase module includes mammalian STE20-like protein kinase 1 (MST1) and MST2, large tumor suppressor 1 (LATS1) and LATS2, together with the adaptor proteins Salvador homologue 1 (SAV1), MOB kinase activator 1A (MOB1A), and MOB1B [6-11]. These inhibitory kinase modules regulate tissue growth by suppressing the transcription module such as the Yes-associated protein-1 (YAP1) and transcriptional co-activator with PDZ-binding motif (TAZ) [12]. YAP1 is a key conduit for Hippo pathway regulation and output [5]. Hyperactivation of YAP1 is widespread in cancers. Previous studies showed that cancer features such as cancer stem cell properties, epithelial-mesenchymal transition (EMT), increased migration, and metastasis are regulated by YAP1 [11]. In a study of the mutational landscape across 12 major cancer types, significantly mutated genes of hippo signaling were found in several cancers included HNSCC [13].
In oral squamous cell carcinoma, nuclear YAP1 accumulation marked premalignant dysplastic regions of the oral epithelium and YAP1 promoted tumorigenic phenotypes and a transcriptional program associated with tumor progression [14]. YAP1 expression in combination with p63 can facilitate identification of HNSCC tumors from hyperplastic and benign tissues [15]. Furthermore, YAP1 upregulation is known to be related with resistance to anticancer drugs such as docetaxel, cisplatin, and cetuximab in several cancers [16-18]. In oral squamous cell carcinoma cell lines, nuclear translocation of YAP1 correlated with the acquisition of cisplatin resistance [19]. In HNSCC, amplification of the YAP1 gene was linked to cetuximab resistance in cell lines [20]. However, the clinical relevance of YAP1 activation has still not been examined in HNSCC.
In the current study, we systematically characterized the genomic data from multiple cohorts of patients with HNSCC and found a molecular subtype characterized by YAP1 activation and poor prognosis. We analyzed the gene alterations associated with YAP1 using multiplatforms. We further showed that YAP1 may play roles in resistant to immunotherapy.
RESULTS
Activation of YAP1 in HNSCC
Although activation of YAP1 has been reported in many cancers [5, 16, 18, 21-24], its relative activity across all cancer types has not been systematically examined. Because amplification is one common activation mechanism of YAP1 [25, 26], we first examined copy number alteration of YAP1 in 26 cancers by using genome copy number data from TCGA. HNSCC is the second most common YAP1-amplified cancer type (8.6%) after cervical cancer (12.6%) (Supplementary Figure 1A), suggesting that YAP1 activation might be a critical genetic event for development of HNSCC. Expression of YAP1 was significantly correlated with copy number alteration (Figure 1A, r = 0.781, P < 0.001), further supporting functional roles of YAP1 in the development of HNSCC. Interestingly, expression of YAP1 was substantially higher in many tumors without YAP1 amplification, suggesting that copy number alteration is not the only mechanism of YAP1 activation in HNSCC (Figure 1A).
Figure 1: Activation of YAP1 in HNSCC. (A) Scatter plots between mRNA expression and copy number alteration of YAP1 in TCGA cohort. The expression of YAP1 was substantially higher in many tumors without YAP1 amplification (in circle). (B) Expression patterns of 292 genes in the YAP1 signature. The data are presented in matrix format, in which each row represents an individual gene and each column represents a tissue specimen. Each cell in the matrix represents the level of expression of a gene feature in an individual sample. The red and green cells reflect relatively high and low expression levels, respectively. YA, YAP1 activated; YI, YAP1 inactivated. (C) Kaplan-Meier plots of the YA and YI patients in TCGA cohort. Overall survival time of patients with the YA subtype was significantly worse than that of patients with the YI subtype. P values were calculated using log-rank tests. +, censored data.
Because the best-known molecular activity of YAP1 is transcription activation, we sought to identify potential downstream targets of YAP1 by finding genes with expression significantly correlated with copy number alterations (P < 0.001 and Pearson correlation coefficient > 0.2 or < –0.2), yielding 652 genes (copy number-associated genes). We next constructed a prediction model with these potential downstream target genes to estimate YAP1 activity in HNSCC by using the BCCP model. Briefly, tumors were first divided into two groups according to copy number alterations: YAP1-high (GISTIC score >2) and YAP1-low (GISTIC score <2). As expectedly, the vast majority of YAP1-high tumors showed a high probability (>0.5) of YAP1 activity when BCCP was applied to gene expression data. Interestingly, many tumors without YAP1 amplification also showed a high probability of YAP1 activity (Figure 1A), further supporting our notion that YAP1 amplification may not be the only mechanism for YAP1 activation in HNSCC. Because many tumors showed high YAP1 expression without YAP1 amplification, we also identified genes with expression significantly correlated with mRNA expression of YAP1 (P < 0.001 and Pearson correlation coefficient > 0.2 or < –0.2). This search yielded 4552 genes (mRNA-associated genes). To identify key downstream targets of YAP1 in HNSCC, we further selected 292 genes that were shared in both gene lists (Supplementary Figure 1B, Supplementary Table 1). To assess the clinical relevance of YAP1 activation in HNSCC, patients were re-stratified according to BCCP probability: YAP1-active (YA >0.5) and YAP1-inactive (YI <0.5) (Figure 1B). Overall survival (OS) time of patients with the YA subtype was significantly worse than that of patients with the YI subtype, suggesting that YAP1 activation is significantly associated with poor prognosis in HNSCC (Figure 1C).
YAP1 signature was associated with the prognosis of HNSCC
Having shown that gene expression signature accurately reflecting activation of YAP1 and had a significant association with prognosis, we next sought to validate its association in the independent cohorts (Leipzig, Greek, MD Anderson, and Seattle cohorts) (Figure 2A). Patients in these cohorts were stratified to the YA or YI group by BCCP classifier (Supplementary Figure 2, Table 1). When this stratification was applied to the Leipzig cohort, 115 (42.6%) of 270 patients were classified as YA (Table 1), and their prognosis was significantly worse than for patients in the YI group (OS 27.7% vs. 38.3% at 5 years, P = 0.021) (Figure 2B). Of 109 patients in the Greek cohort, 30 (27.5%) were classified as YA, and the patients with the YA subtype had a significantly worse disease-free survival than did patients with the YI subtype (37.8% vs. 65.1% at 5 years, P = 0.008) (Figure 2B). In the MD Anderson cohort, 19 (25.7%) of 74 patients were classified as YA and had significantly worse disease-free survival than those with the YI subtype (43.3% vs. 75.8%, P = 0.008) (Figure 2B). Consistent with previous observation, patients classified as YA subtype in the Seattle cohort (35 patients or 36.1%) also showed significantly poorer prognosis (P = 0.019) (Figure 2B).
Figure 2: Construction of a prediction model using YAP1 gene signature and clinical significance of YAP1 activation in HNSCC. (A) Schematic diagram of the strategy used to construct the prediction model and evaluate predicted HNSCC outcomes according to the gene expression signature. (B) Kaplan-Meier plots of the YA and YI patients in four independent test cohorts. P values were calculated using log-rank tests. +, censored data.
Table 1: Patient’s characteristics in 5 cohort
TCGA cohort | Leipzig cohort | Greece cohort | MDACC cohort | Seattle cohort | |
|---|---|---|---|---|---|
Number of patients | 513 | 270 | 109 | 74 | 97 |
Gender | |||||
Male | 370 (73.7%) | 223 (82.6%) | 104 (95.4%) | 58 (78.4%) | 66 (68.0%) |
Female | 132 (26.3%) | 47 (17.4%) | 5 (4.6%) | 16 (21.6%) | 31 (32.0%) |
Age (mean ± SD) | 60.9 ± 11.9 | 60.1 ± 10.0 | 63.3 ± 10.1 | 58.1 ± 13.6 | NA |
Anatomic site | |||||
Oral cavity | 301 (60.0%) | 83 (30.7%) | 0 | 71 (95.9%) | 86 (88.7%) |
Oropharynx | 79 (15.7%) | 102 (37.8%) | 0 | 3 (4.1%) | 11 (11.3%) |
Larynx | 113 (22.5%) | 48 (17.8%) | 109 (100%) | 0 | 0 |
Hypopharynx | 9 (1.8%) | 33 (12.2%) | 0 | 0 | 0 |
others | 0 | 4 (1.5%) | 0 | 0 | 0 |
Primary tumor | |||||
T1 | 33 (6.8%) | 35 (13.0%) | NA | 3 (4.1%) | NA |
T2 | 147 (30.2%) | 80 (29.6%) | NA | 27 (36.5%) | NA |
T3 | 129 (26.5%) | 58 (21.5%) | NA | 28 (37.8%) | NA |
T4 | 178 (36.6%) | 97 (35.9%) | NA | 16 (21.6%) | NA |
Regional lymph node | |||||
N0 | 238 (49.5%) | 94 (34.8%) | 42 (56.8%) | NA | NA |
N1 | 79 (16.4%) | 32 (11.9%) | 13 (17.6%) | NA | NA |
N2 | 155 (32.2%) | 132 (48.9%) | 19 (25.7%) | NA | NA |
N3 | 9 (1.9%) | 12 (4.4%) | 0 | NA | NA |
Stage | |||||
I | 20 (4.1%) | 18 (6.7%) | 12 (11.0%) | 3 (4.1%) | 30 (30.9%) |
II | 96 (19.6%) | 37 (13.7%) | 18 (16.5%) | 16 (21.6%) | 11 (11.3%) |
III | 101 (20.7%) | 37 (13.7%) | 36 (33.0%) | 15 (20.3%) | 15 (15.5%) |
IV | 272 (55.6%) | 178 (65.9%) | 43 (39.4%) | 40 (54.1%) | 41 (42.3%) |
HPV status | |||||
Positive | 68 (19.9%) | 60 (23.4%) | NA | NA | 0 |
Negative | 274 (80.1%) | 196 (76.6%) | NA | NA | 97 (100%) |
Tobacco use | |||||
Never | 114 (23.3%) | 48 (17.8%) | 1 (0.9%) | 15 (20.3%) | NA |
Yes | 376 (76.7%) | 222 (82.2%) | 108 (99.1%) | 59 (79.7%) | NA |
Alcohol use | |||||
Never | 154 (42.1%) | 31 (11.5%) | 51 (46.8%) | NA | NA |
Yes | 212 (57.9%) | 239 (88.5%) | 58 (53.2%) | NA | NA |
YAP1 signature | |||||
YAHSC | 187 (36.5%) | 115 (42.6%) | 30 (27.5%) | 19 (25.7%) | 35 (36.1%) |
YIHSC | 326 (63.5%) | 155 (57.4%) | 79 (72.5%) | 55 (74.3%) | 62 (63.9%) |
We next performed analyses to determine whether the prognostic effect of the YAP1 signature was independent of other clinical variables. Three cohorts (TCGA, Leipzig, and MD Anderson) with available OS data were pooled for a Cox proportional hazards model (n = 808). In univariate analysis, YAP1 signature (YA group vs. YI group), sex, age (<60 years old vs. ≥60 years old), anatomic site (oropharynx vs. non-oropharynx), and primary tumor (T4 vs. non-T4) were significant prognostic factors (Table 2). To address a potential confounding effect, we performed a multivariate risk analysis using the Cox proportional hazards model and found that the YAP1 signature was the most significant predictor of OS (hazard ratio 1.31, 95% confidence interval 1.03-1.67, P = 0.028).
Table 2: Univariate and multivariate Cox Proportinal Hazard Regression Analyses of variables affecting 5-year overall survival rate (patients data: TCGA cohort, Leipzig cohort, MDACC cohort; n=808)
Univariate | Multivariate | |||
|---|---|---|---|---|
HR (95% CI) | P-value | HR (95% CI) | P-value | |
YAP1 signature (YAP1 activated) | 1.41 (1.12-1.77) | 0.004 | 1.31 (1.03-1.67) | 0.028 |
Gender (male) | 0.77 (0.60-1.0) | 0.05 | 0.85 (0.65-1.12) | 0.247 |
Age (≥ 60 y) | 1.30 (1.04-1.64) | 0.024 | 1.26 (0.99-1.61) | 0.066 |
Anatomic site (oropharynx) | 0.54 (0.32-0.92) | 0.025 | 0.63 (0.37-1.10) | 0.104 |
Primary tumor (T4) | 2.03 (1.00-4.12) | 0.049 | 1.81 (0.88-3.73) | 0.108 |
Regional lymph node (N+) | 1.02 (0.81-1.28) | 0.9 | 1.01 (0.76-1.34) | 0.971 |
Stage (stage III & IV) | 1.21 (0.91-1.62) | 0.19 | 1.08 (0.74-1.56) | 0.695 |
Next, we assessed whether YAP1 activation may have a potential association with the patient’s response to radiation therapy. We pooled the patient data from TCGA and MD Anderson cohorts with available radiation data (n = 458). For the patients who underwent radiation therapy, the survival of patients with the YA subtype was significantly worse than for patients with the YI subtype (P = 0.03) (Supplementary Figure 3A). For the YA group, radiation therapy did not improve the survival rate (P = 0.15) (Supplementary Figure 3B). However, radiation therapy improved the survival rate of the YI group (P = 0.04) (Supplementary Figure 3C), suggesting that YAP1 activation might be associated with radioresistance in HNSCC.
Copy number alteration associated with YAP1 activation
To further examine the genomic properties associated with YAP1 activation, we compared the copy number alteration between the YA and YI groups (Figure 3). Somatic copy numbers were obtained for 513 samples in TCGA cohort. When frequency of significantly reoccurring alterations was compared using Fisher’s exact test, the YA group had amplification of 7p11.2 (EGFR) and 8q11.21 (SNAI2). The YA group had deletion of 8p23 (CSMD1), 9p21.3 (CDKN2A), 9q34.3 (NOTCH1), and 18q21.2 (SMAD4). CSMD1, CDKN2A, NOTCH1, and SMAD4 are known to be tumor-suppressor genes. The YI group was distinguished by amplification of 3q26/28 (PIK3CA, SOX2, and TP63) and deletion of 11q23.1 (Supplementary Table 2). Gains in the 3q, 5p, and 8q chromosomal regions were features of HNSCC and lung squamous cell carcinoma [27]. Interestingly, while gain in the 3q region was higher in the YI group, gain in the 8q region was higher for the YA group. There was no difference in 5p between the two groups.
Figure 3: DNA copy number alterations in TCGA HNSCC cohort. (A) Heat map of copy number alterations of tumors stratified by YAP1 subtype. (B) Copy number gains and losses in YA and YI subtypes. For each gene, the mean of the segmented GISTIC copy number values, in YA and YI subtypes were computed and plotted in genomic order.
Somatic mutation associated with YAP1 activation
To assess co-occurring somatic mutations with YAP1 activation in HNSCC, we also analyzed somatic mutation data in TCGA cohort (n = 493). When 30 significantly mutated genes in HNSCC were analyzed for association with YAP1 activation [13, 27], TP53 (82.6% vs. 64.4%, P = 1.4 × 10-5) and CDKN2A (28.3% vs. 19.1%, P = 0.02) were identified as co-occurring mutated genes (Figure 4, Supplementary Table 3). The mutations of NFE2L2 (1.1% vs. 7.8%, P = 7.0 × 10-4), PTEN (0.5% vs. 4.2%, P = 0.022), SYNE1 (14.1% vs. 22.3%, P = 0.025), and NSD1 (8.2% vs. 14.9%, P = 0.033) were mutually exclusive for YAP1 activation. Interestingly, AJUBA has been implicated in the Hippo pathway and was segregated predominantly in HPV-negative tumors [27]. However, AJUBA did not significantly differ between groups (7.1% vs. 5.8%, P = 0.57).
Figure 4: Somatic mutation in HNSCC according to two subgroups. The mutations of TP53 and CDKN2A were identified as co-occurring mutated genes in YA subgroup. The mutations of NFE2L2, PTEN, SYNE1, and NSD1 were mutually exclusive for YAP1 activation.
miRNAs associated with YAP1 activation
To explore potential interactions between YAP1 activation and miRNAs in HNSCC, we analyzed miRNA expression data from TCGA cohort (n = 473 tumors). Expression of 12 miRNAs was significantly associated with YAP1 activation (fold change >2, P < 0.001) (Figure 5A, Supplementary Table 4). miRNA-1-2, miRNA-133a, miRNA-133b, and miRNA-206 are known to be frequently downregulated in cancer [22]. They are also known to be involved in the regulation of molecular networks such as Wnt signaling pathway, tight junction, and MAPK signaling pathway in cancer [22]. Because these miRNAs were significantly upregulated in the YA subtype, they might be regulated through YAP1-associated signaling pathways.
Figure 5: Molecular characteristics of the YA and YI subtypes. (A) miRNAs that are differentially expressed in the YA and YI subtypes. (B) Association between YAP1 signature and the other subtypes of HNSCC. HPV infection status was most significantly associated with the two subtypes. Of 105 YA patients, 104 were HPV-negative (P = 4.29 × 10-7).
Because many tumors in the YA subtype had high expression of YAP1 without YAP1 amplification (Figure 1A), we hypothesized that certain miRNAs might play roles in the regulation of YAP1 in these tumors [copy number-independent (CNI)-YA]. Thus, we sought to find miRNAs whose expression is specifically higher in CNI-YA than in copy number-dependent (CND) YA, YI, or normal head and neck tissues. We first selected miRNAs whose expression significantly differed between CNI-YA and CND-YA (fold ratio >1.5, P < 0.01) and further selected miRNAs whose expression was significantly different between CNI-YA and normal head and neck tissues (P < 0.05), yielding 21 miRNAs (Supplementary Table 5). Of these miRNAs, only expressions of miRNA-187 and miRNA-675 were significantly higher in CNI-YA than in YI or in normal tissues (Supplementary Figure 4). In good agreement with our hypothesis, many of the predicted target genes of both miRNAs are negative upstream regulators of YAP1 in the Hippo pathway (Supplementary Table 6, Supplementary Figure 5), suggesting that YAP1 in CNI-YA tumors might be activated through miRNA-mediated mechanisms.
Clinical characteristics associated with YAP1 activation
We next assessed the association of the two subtypes with clinical characteristics of HNSCC in TCGA cohort. Among clinically recognized features such as lymph-node metastasis and T stages, YAP1 activation was significantly associated with HPV status (P = 4.29 × 10-7, Figure 5B). In TCGA cohort, HPV status data were available from only 279 of the 513 patients. Of 36 HPV-positive patients, most (n = 35) were in the YI subgroup. We also assessed the association of YAP1 activation with four previously recognized molecular subtypes of HNSCC (Supplementary Table 7) [28]. YAP1 activation was most associated with the basal subtype, as reflected in the fact that the highest proportion of patients (58.7%) in the YA subtype was from the basal subtype (P = 2.2 × 10-6, Figure 5B). TheYI subgroup lacked any significant association with the four subtypes as it was found in similar proportions in the four subtypes: atypical subtype (33.3%), classical subtype (24.9%), mesenchymal subtype (22.0%), and basal subtype (19.7%).
Association of YAP1 activity with host immune activity
Recent clinical trial demonstrated clinically meaningful anti-tumor activity of pembrolizumab, humanized antibody against immune checkpoint inhibitor PD-1, in HNSCC [29]. Since study also showed that composite scores from expression of 6 interferon γ (INFG)-related genes (CXCL9, CXCL10, IDO1, IFNG, HLA-DRA, and STAT1) were good predictor for identifying responders of pembrolizumab treatment, we next assessed association of two subtypes with potential response to permbrolizumab treatment by generating INFG composite scores in TCGA cohort. Interestingly, YA subtype has significantly low INFG composite scores (P =8.1 x 10-5) (Figure 6A), suggesting that YAP1 may suppress immune activity related immune checkpoint regulation. Further support of the idea is supplied by significant negative correlation between IFNG composite scores and Bayesian probability of active YAP1 in TCGA cohort (Figure 6B). Consistent with results from TCGA cohort, IFNG scores were significantly lower in YA subtype (P = 0.001) and negatively correlated with Bayesian probability of active YAP1 in Leipzig cohort (Figure 6C and 6D).
Figure 6: Negative association of YAP1 activity with IFNG scores in HNSCC (A, C). Comparison of the IFNG scores between YA subgroup and YI subgroup in TCGA and Leipzig cohort. (B, D) Scatter plots between IFNG scores and Bayesian probability of YAP1 activity is TCGA and Leipzig cohort.
DISCUSSION
By systematically analyzing genomic copy number data and mRNA expression data of YAP1 in HNSCC, we identified a gene expression signature reflecting YAP1 activation (YAP1 signature) that is significantly associated with the prognosis of patients with HNSCC. Robustness of the signature was validated in five independent cohorts of patients with HNSCC, as patients with the YA subtype had worse survival rates than patients with the YI subtype in all examined cohorts. The YAP1 signature was an independent prognostic factor. In multiplatform analysis, YA patients had gain of EGFR and SNAI2; loss of tumor-suppressor genes such as CSMD1, CDKN2A, NOTCH1, and SMAD4; and high mutation rates of TP53 and CDKN2A. YI patients were characterized by gain of PIK3CA, SOX2, and TP63; deletion of 11q23.1; and high mutation rates of NFE2L2, PTEN, SYNE1, and NSD1. To our knowledge, we are the first to report the clinical significance of YAP1 activation in HNSCC and demonstrate that the YAP1 signature can be used as a prognostic biomarker for HNSCC.
Our study results are in good agreement with previous observations in other cancers [18, 30-34]. Activation of YAP1 has been correlated with poor prognosis for colorectal cancer and poor response to cetuximab [18]. Likewise, activation of YAP1 was significantly associated with poor prognosis in ovarian cancer, liver cancer, and gastric cancers [18, 32-34]. Previous studies showed that YAP1 had a role in the onset, progression, and drug resistance in HNSCC [14, 15, 19, 20]. Hiemer et al. reported that YAP1 and TAZ activity drives oral squamous cell carcinoma proliferation, survival, and migration in vitro and vivo [14]. YAP expression was elevated in tumor compared with benign tissues and was associated with nodal metastasis [15]. A potential association of YAP1 with resistance to radiation therapy was also supported by previous studies, as amplification of YAP1 was correlated with cetuximab sensitivity in HNSCC [20] and knockdown of YAP1 increased the sensitivity to cisplatin in vitro [19].
Comparisons of copy number alteration and somatic mutation showed significant differences between the two subtypes. The YA group had more loss of copy number and a higher somatic mutation rate of the cell cycle–related gene CDKN2A. In good agreement with a previous study demonstrating that deletion of CDKN2A and inactivating mutation were associated with HPV-negative tumors [27], the vast majority of tumors with the YA subtype were HPV-negative tumors, suggesting that YAP1 activation and inactivation of CDKN2A might be associated genetic events. Thus, it will be interesting to determine in future investigations whether YAP1 can downregulate CDKN2A. Tumors of the YA subtype were also characterized by gain of EGFR and SNAI2. SNAI2 has five zinc finger domains that play a pivotal role during embryo development and mesenchymal tumorigenesis and has been found to be overexpressed in several cancers, and it also promotes invasion in lung adenocarcinoma, glioma, and ovarian, cervical, and pancreatic cancers and is a prognostic marker in some cancers [24, 35]. The YA subtype was also characterized by a higher mutation rate of the tumor-suppressor gene associated with cell survival (TP53) and loss of copy number of tumor-suppressor genes (CSMD1, NOTCH1, and SMAD4), while the YI subtype was characterized by high mutation rates of NFE2L2, PTEN, SYNE1, and NSD1. Spectrin repeat containing nuclear envelope protein 1 (SYNE1) gene has been implicated in the regulation of nuclear polarity, a process that operates upstream of NOTCH1 in squamous epithelia [35]. SYNE1 mutation was associated with autosomal recessive cerebellar ataxia [24] and is known to be associated with glioblastoma and lung, ovarian, colorectal, and head and neck cancers [35, 36]. The nuclear receptor binding SET domain protein 1 (NSD1), a histone methyltransferase, was found to be frequently mutated in the clear cell variant of renal cell carcinoma and associated with DNA hypomethylation [27].
In assessing the association of the YAP1 subgroup with the four previously discovered molecular subtypes [28], we found that the YA group had the most similarity with the basal subtype. The basal subtype was characterized by inactivation of NOTCH1 and co-amplified 11q13/q22 (21) in which YAP1 resided. Interestingly, 94.5% of atypical subtype was YI subtype (Supplementary Table 7). The atypical subtype was characterized by enrichment of HPV-positive tumor and activating mutations in exon 9 that contain PIK3CA (21). In good agreement with this, the YI subtype included most of the HPV-positive tumors and featured the amplification of copy number of PIK3CA. While the YAP1-associated subtype shares some molecular features with the four previously recognized molecular subtypes, it is important to point out that this subtype is clinically relevant as reflected in prognostic difference and response to radiation therapy, while clinical association of the four molecular subtypes has not been clearly demonstrated yet [28].
Recent advances in our understanding of onco-immunology has led to the development of immunotherapies. Particularly, a blockade of checkpoint molecules such as CTLA-4 and PD-1 has emerged as a novel therapeutic approach in oncology [37-39]. PD-1 is a negative co-stimulatory receptor and a strong inhibitor of T cell response [40]. Pembrolizumab targeting PD-1 has been approved for the treatment of multiple cancers including metastatic melanoma [41]. Six-gene based IFNG composite scores were developed and tested as predictor for response to pembrolizumab treatment in HNSCC [29]. When immune activity reflecting potential response to pembrolizumab in HNSCC tumors was assessed by using IFNG scores, tumors in YA subtype have lower scores and IFNG scores were negatively correlated with YAP1 activity, suggesting that YAP1 may be involved in regulation of host immunity against cancer cells. Thus, YA subtype might be more resistant to immunotherapy. Since YA subtype is poor prognostic, it would be also interested to see if YAP1-mediated low immune activity contributes to aggressiveness of cancer cells in YA subtype in future study.
This study has some limitations. First, all cohorts examined in our study were retrospective cohorts. Second, while patients in the YI subgroup showed better survival than those in the YA subgroup when treated with radiation therapy, an interaction test failed to show significant interaction between subgroups and radiation therapy. Thus, these observations should be validated in future prospective study.
In conclusion, our newly discovered clinically relevant subgroup may help facilitate rational design of clinical studies by stratifying patients according to their prognostic risk and response to therapy. Prospective studies are needed to validate the clinical usefulness of the two subgroups.
MATERIALS AND METHODS
Patients and gene expression data
This study used patient data from five independent sources. Gene expression, mutation, miRNA, and copy number data of The Cancer Genome Atlas (TCGA) cohort were downloaded from the UCSC Cancer Genomics Browser (https://genome-cancer.ucsc.edu/) [27]. mRNA expression, miRNA expression, copy number alteration, and mutation data were available from 513, 473, 513, and 493 patients respectively. We also used the gene expression and clinical data for four independent cohorts available from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). These four cohorts consisted of data from the Institute for Medical Informatics, Statistics and Epidemiology (Leipzig cohort, GSE65858, n = 270) [42], Aristotle University of Thessaloniki (Greek cohort, GSE27020, n = 109) [43], MD Anderson Cancer Center (MD Anderson cohort, GSE42743, n = 74) [44], and Fred Hutchinson Cancer Research Center (Seattle cohort, GSE41613, n = 97) [44]. Table 1 shows the pathologic and clinical characteristics of the patients in all five cohorts.
Statistical analysis
The BRB-ArrayTools software program (http://brb.nci.nih.gov/BRB-ArrayTools/) was used for analysis of gene expression data [45]. The R language environment (http://www.r-project.org) was used for other statistical analyses. Raw data on the patient cohorts were normalized using a robust multiarray averaging method [46]. A stringent threshold was used to minimize the number of false-positive findings. Pearson correlation was used for correlation analysis. We estimated prognoses using Kaplan-Meier plots and the log-rank test. We used univariate and multivariate Cox proportional hazards regression analyses to evaluate independent prognostic factors associated with survival. Independent t-test was used to compare values between two groups. Fisher’s exact test was used to evaluate the frequency difference of copy number alteration and somatic mutation. P values less than 0.05 indicated statistical significance, and all statistical tests were two-tailed. A heatmap was generated using the Cluster and TreeView software programs [47]. Target genes of miRNAs were predicted by using miRWalk 2.0 (http://zmf.umm.uni-heidelberg.de/) [48].
Gene expression data and construction of the prediction model
Gene expression data of TCGA cohort were sequenced by Illumina HiSeq2000, the Leipzig cohort by Illumina HumanHT-12 V4.0 expression beadchip, the Greek cohort by Affymetrix Human Genome U133A Array, and the MD Anderson and Seattle cohorts by Affymetrix Human Genome U133 Plus 2.0 Array. All gene expression data were standardized independently across all samples before they were integrated and analyzed together. The strategy used to develop and validate the prediction model on the basis of the gene expression signature and to estimate predictive accuracy was adopted from previous studies [49-51]. To find YAP1-specific genes in HNSCC, we applied a double correlation approach to gene expression data from TCGA cohort. We first identified YAP1-associated genes by using the correlation between copy number of YAP1 and mRNA expression of each gene (copy number-associated genes). Likewise, we also identified second YAP1-associated genes by using the correlation between mRNA expression of YAP1 and mRNA expression of each gene (mRNA-associated genes). Genes were selected if the P value was less than 0.001 and the correlation coefficient was more than 0.2 or less than –0.2. Expression patterns of 292 shared genes in two YAP1-associated gene lists were considered as the YAP1 signature in HNSCC and used for construction of the prediction model. We used data from TCGA cohort as the training set and data from the Leipzig, Greek, MD Anderson, and Seattle cohorts as test sets. The expression patterns of the 292 genes from TCGA cohort were combined to form a classifier according to a Bayesian compound covariate (BCCP) predictor [52]. The robustness of the classifier was estimated using a misclassification rate determined during leave-one-out cross-validation in the training set. The BCCP classifier estimated the likelihood that an individual patient had either a YAP1-active subtype or YAP1-inactive subtype.
Copy number analysis
We used the HNSCC copy number data (gistic2) in Cancer Browser (https://genome-cancer.ucsc.edu). Threshold copy number at reoccurring alteration peaks from GISTIC analysis of BRB CGH-Tools was used for comparison of focal peak frequency across subgroups by YAP1 signature. Fisher’s exact test was used for frequency comparisons of significantly reoccurring alterations by YAP1 signature. For comparison of amplification, only high-level events were considered (gistic ≥2), and for deletion all events were used as described in a previous study (21).
Association with YAP1 signature and other subtypes of HNSCC
For assessing the association with previously published molecular classifications of HNSCC [28], we evaluated four previously established molecular classifications: atypical, classical, basal, and mesenchymal. The predictor was adjusted from 790 genes. Each sample was assigned an expression subtype using the BCCP predictor.
Association with YAP1 signature and interferon interferon γ signature
A six-gene signature of interferon γ-related genes (CXCL9, CXCL10, IDO1, IFNG, HLA-DRA, and STAT1) was used in TCGA cohort. The interferon γ score was calculated as the average of the value of the six genes. The interferon γ score was compared with t-test between YA group and YI group.
Author contributions
YGE, DL, YCL, BHS, EHK, SYY, KHK and JL conceived and designed the research. YGE, DL, BHS, EHK and SYY collected and analyzed the data. YCL and KHK performed the experiments. YGE wrote the manuscript. JL supervised the research. All authors read and approved the final manuscript.
CONFLICTS OF INTEREST
The authors declare that they have no competing interests.
FUNDING
This research was supported by the Bio & Medical Technology Development Program of the National Research Foundation (NRF) & funded by the Korean government (NRF-2015M3A9E2029189) National Institutes of Health grants CA150229 and CA016672, 2016 Institutional Research Grant (IRG) and 2016 Sister Institute Network Fund (SINF) grant from The University of Texas MD Anderson Cancer Center, and a Research Initiative grant from the Korean Research Institute of Bioscience and Biotechnology (KRIBB). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
REFERENCES
1. Jemal A, Bray F, Center MM, Ferlay J, Ward E, Forman D. Global cancer statistics. CA Cancer J Clin. 2011; 61:69-90.
2. Vokes EE, Weichselbaum RR, Lippman SM, Hong WK. Head and neck cancer. N Engl J Med. 1993; 328:184-194.
3. Argiris A, Karamouzis MV, Raben D, Ferris RL. Head and neck cancer. Lancet. 2008; 371:1695-1709.
4. Gupta S, Kong W, Peng Y, Miao Q, Mackillop WJ. Temporal trends in the incidence and survival of cancers of the upper aerodigestive tract in Ontario and the United States. Int J Cancer. 2009; 125:2159-2165.
5. Harvey KF, Zhang X, Thomas DM. The Hippo pathway and human cancer. Nat Rev Cancer. 2013; 13:246-257.
6. Harvey KF, Pfleger CM, Hariharan IK. The Drosophila Mst ortholog, hippo, restricts growth and cell proliferation and promotes apoptosis. Cell. 2003; 114:457-467.
7. Udan RS, Kango-Singh M, Nolo R, Tao C, Halder G. Hippo promotes proliferation arrest and apoptosis in the Salvador/Warts pathway. Nat Cell Biol. 2003; 5:914-920.
8. Pantalacci S, Tapon N, Léopold P. The Salvador partner Hippo promotes apoptosis and cell-cycle exit in Drosophila. Nat Cell Biol. 2003; 5:921-927.
9. Wu S, Huang J, Dong J, Pan D. Hippo encodes a Ste-20 family protein kinase that restricts cell proliferation and promotes apoptosis in conjunction with salvador and warts. Cell. 2003; 114:445-456.
10. Kango-Singh M, Nolo R, Tao C, Verstreken P, Hiesinger PR, Bellen HJ, Halder G. Shar-pei mediates cell proliferation arrest during imaginal disc growth in Drosophila. Development. 2002; 129:5719-5730.
11. Moroishi T, Hansen CG, Guan KL. The emerging roles of YAP and TAZ in cancer. Nat Rev Cancer. 2015; 15:73-79.
12. Huang J, Wu S, Barrera J, Matthews K, Pan D. The Hippo signaling pathway coordinately regulates cell proliferation and apoptosis by inactivating Yorkie, the Drosophila Homolog of YAP. Cell. 2005; 122:421-434.
13. Kandoth C, McLellan MD, Vandin F, Ye K, Niu B, Lu C, Xie M, Zhang Q, McMichael JF, Wyczalkowski MA. Mutational landscape and significance across 12 major cancer types. Nature. 2013; 502:333-339.
14. Hiemer SE, Zhang L, Kartha VK, Packer TS, Almershed M, Noonan V, Kukuruzinska M, Bais MV, Monti S, Varelas X. A YAP/TAZ-regulated molecular signature is associated with oral squamous cell carcinoma. Mol Cancer Res. 2015; 13:957-968.
15. Ge L, Smail M, Meng W, Shyr Y, Ye F, Fan KH, Li X, Zhou HM, Bhowmick NA. Yes-associated protein expression in head and neck squamous cell carcinoma nodal metastasis. PLoS One. 2011; 6:e27529.
16. Hall CA, Wang R, Miao J, Oliva E, Shen X, Wheeler T, Hilsenbeck SG, Orsulic S, Goode S. Hippo pathway effector Yap is an ovarian cancer oncogene. Cancer Res. 2010; 70:8517-8525.
17. Song S, Honjo S, Jin J, Chang SS, Scott AW, Chen Q, Kalhor N, Correa AM, Hofstetter WL, Albarracin CT. The Hippo coactivator YAP1 mediates EGFR overexpression and confers chemoresistance in esophageal cancer. Clin Cancer Res. 2015; 21:2580-2590.
18. Lee KW, Lee SS, Kim SB, Sohn BH, Lee HS, Jang HJ, Park YY, Kopetz S, Kim SS, Oh SC. Significant association of oncogene YAP1 with poor prognosis and cetuximab resistance in colorectal cancer patients. Clin Cancer Res. 2015; 21:357-364.
19. Yoshikawa K, Noguchi K, Nakano Y, Yamamura M, Takaoka K, Hashimoto-Tamaoki T, Kishimoto H. The Hippo pathway transcriptional co-activator, YAP, confers resistance to cisplatin in human oral squamous cell carcinoma. Int J Oncol. 2015; 46:2364-2370.
20. Jerhammar F, Johansson AC, Ceder R, Welander J, Jansson A, Grafström RC, Söderkvist P, Roberg K. YAP1 is a potential biomarker for cetuximab resistance in head and neck cancer. Oral Oncol. 2014; 50:832-839.
21. Mo JS, Park HW, Guan KL. The Hippo signaling pathway in stem cell biology and cancer. EMBO Rep. 2014; 15:642-656.
22. Shao DD, Xue W, Krall EB, Bhutkar A, Piccioni F, Wang X, Schinzel AC, Sood S, Rosenbluh J, Kim JW. KRAS, YAP1 converge to regulate EMT and tumor survival. Cell. 2014; 158:171-184.
23. Kang W, Tong JH, Chan AW, Lee TL, Lung RW, Leung PP, So KK, Wu K, Fan D, Yu J. Yes-associated protein 1 exhibits oncogenic property in gastric cancer and its nuclear accumulation associates with poor prognosis. Clin Cancer Res. 2011; 17:2130-2139.
24. Liu JY, Li YH, Lin HX, Liao YJ, Mai SJ, Liu ZW, Zhang ZL, Jiang LJ, Zhang JX, Kung HF. Overexpression of YAP 1 contributes to progressive features and poor prognosis of human urothelial carcinoma of the bladder. BMC Cancer. 2013; 13:1.
25. Zender L, Spector MS, Xue W, Flemming P, Cordon-Cardo C, Silke J, Fan ST, Luk JM, Wigler M, Hannon GJ. Identification and validation of oncogenes in liver cancer using an integrative oncogenomic approach. Cell. 2006; 125:1253-1267.
26. Modena P, Lualdi E, Facchinetti F, Veltman J, Reid JF, Minardi S, Janssen I, Giangaspero F, Forni M, Finocchiaro G. Identification of tumor-specific molecular signatures in intracranial ependymoma and association with clinical characteristics. J Clin Oncol. 2006; 24:5223-5233.
27. Cancer Genome Atlas Network. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature. 2015; 517:576-582.
28. Walter V, Yin X, Wilkerson MD, Cabanski CR, Zhao N, Du Y, Ang MK, Hayward MC, Salazar AH, Hoadley KA. Molecular subtypes in head and neck cancer exhibit distinct patterns of chromosomal gain and loss of canonical cancer genes. PLoS One. 2013; 8:e56823.
29. Seiwert TY, Burtness B, Mehra R, Weiss J, Berger R, Eder JP, Heath K, McClanahan T, Lunceford J, Gause C. Safety and clinical activity of pembrolizumab for treatment of recurrent or metastatic squamous cell carcinoma of the head and neck (KEYNOTE-012): an open-label, multicentre, phase 1b trial. Lancet Oncol. 2016; 17:956-965.
30. Sun Z, Xu R, Li X, Ren W, Ou C, Wang Q, Zhang H, Zhang X, Ma J, Wang H. Prognostic value of Yes-associated protein 1 (YAP1) in various cancers: a meta-analysis. PLoS One. 2015; 10:e0135119.
31. Sun PL, Kim JE, Yoo SB, Kim H, Jin Y, Jheon S, Kim K, Lee CT, Chung JH. Cytoplasmic YAP expression is associated with prolonged survival in patients with lung adenocarcinomas and epidermal growth factor receptor tyrosine kinase inhibitor treatment. Ann Surg Oncol. 2014; 21:610-618.
32. Sohn BH, Shim JJ, Kim SB, Jang KY, Kim SM, Kim JH, Hwang JE, Jang HJ, Lee HS, Kim SC. Inactivation of Hippo pathway is significantly associated with poor prognosis in hepatocellular carcinoma. Clin Cancer Res. 2016; 22:1256-1264.
33. Park YY, Sohn BH, Johnson RL, Kang MH, Kim SB, Shim JJ, Mangala LS, Kim JH, Yoo JE, Rodriguez-Aguayo C. Yes-associated protein 1 and transcriptional coactivator with PDZ-binding motif activate the mammalian target of rapamycin complex 1 pathway by regulating amino acid transporters in hepatocellular carcinoma. Hepatology. 2016; 63:159-172.
34. Jeong W, Kim SB, Sohn BH, Park YY, Park ES, Kim SC, Kim SS, Johnson RL, Birrer M, Bowtell DS. Activation of YAP1 is associated with poor prognosis and response to taxanes in ovarian cancer. Anticancer Res. 2014; 34:811-817.
35. Stransky N, Egloff AM, Tward AD, Kostic AD, Cibulskis K, Sivachenko A, Kryukov GV, Lawrence MS, Sougnez C, McKenna A. The mutational landscape of head and neck squamous cell carcinoma. Science. 2011; 333:1157-1160.
36. Masica DL, Karchin R. Correlation of somatic mutation and expression identifies genes important in human glioblastoma progression and survival. Cancer Res. 2011; 71:4550-4561.
37. Garon EB, Rizvi NA, Hui R, Leighl N, Balmanoukian AS, Eder JP, Patnaik A, Aggarwal C, Gubens M, Horn L. Pembrolizumab for the treatment of non–small-cell lung cancer. N Engl J Med. 2015; 372:2018-2028.
38. Hodi FS, O’Day SJ, McDermott DF, Weber RW, Sosman JA, Haanen JB, Gonzalez R, Robert C, Schadendorf D, Hassel JC. Improved survival with ipilimumab in patients with metastatic melanoma. N Engl J Med. 2010; 2010:711-723.
39. Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, Loriot Y, Necchi A, Hoffman-Censits J, Perez-Gracia JL. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet. 2017; 389:67-76.
40. Boussiotis VA. Molecular and biochemical aspects of the PD-1 checkpoint pathway. N Engl J Med. 2016; 375:1767-1778.
41. Lee J, Kefford R, Carlino M. PD-1 and PD-L1 inhibitors in melanoma treatment: past success, present application and future challenges. Immunotherapy. 2016; 8:733-746.
42. Wichmann G, Rosolowski M, Krohn K, Kreuz M, Boehm A, Reiche A, Scharrer U, Halama D, Bertolini J, Bauer U. The role of HPV RNA transcription, immune response-related gene expression and disruptive TP53 mutations in diagnostic and prognostic profiling of head and neck cancer. Int J Cancer. 2015; 137:2846-2857.
43. Fountzilas E, Kotoula V, Angouridakis N, Karasmanis I, Wirtz RM, Eleftheraki AG, Veltrup E, Markou K, Nikolaou A, Pectasides D. Identification and validation of a multigene predictor of recurrence in primary laryngeal cancer. PLoS One. 2013; 8:e70429.
44. Lohavanichbutr P, Méndez E, Holsinger FC, Rue TC, Zhang Y, Houck J, Upton MP, Futran N, Schwartz SM, Wang P. A 13-gene signature prognostic of HPV-negative OSCC: discovery and external validation. Clin Cancer Res. 2013; 19:1197-1203.
45. Simon R, Lam A, Li MC, Ngan M, Menenzes S, Zhao Y. Analysis of gene expression data using BRB-array tools. Cancer Inform. 2007; 3:11-17.
46. Bolstad BM, Irizarry RA, Åstrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003; 19:185-193.
47. Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci U S A. 1998; 95:14863-14868.
48. Dweep H, Gretz N. miRWalk2.0: a comprehensive atlas of microRNA-target interactions. Nat Methods. 2015; 12:697-697.
49. Lee JS, Chu IS, Heo J, Calvisi DF, Sun Z, Roskams T, Durnez A, Demetris AJ, Thorgeirsson SS. Classification and prediction of survival in hepatocellular carcinoma by gene expression profiling. Hepatology. 2004; 40:667-676.
50. Lee JS, Chu IS, Mikaelyan A, Calvisi DF, Heo J, Reddy JK, Thorgeirsson SS. Application of comparative functional genomics to identify best-fit mouse models to study human cancer. Nat Genet. 2004; 36:1306-1311.
51. Lee JS, Heo J, Libbrecht L, Chu IS, Kaposi-Novak P, Calvisi DF, Mikaelyan A, Roberts LR, Demetris AJ, Sun Z. A novel prognostic subtype of human hepatocellular carcinoma derived from hepatic progenitor cells. Nat Med. 2006; 12:410-416.
52. Ramaswamy S, Tamayo P, Rifkin R, Mukherjee S, Yeang CH, Angelo M, Ladd C, Reich M, Latulippe E, Mesirov JP. Multiclass cancer diagnosis using tumor gene expression signatures. Proc Natl Acad Sci U S A. 2001; 98:15149-15154.
