Comprehensive Molecular Characterization of the Hippo Signaling Pathway in Cancer.

Hippo signaling has been recognized as a key tumor suppressor pathway. Here, we perform a comprehensive molecular characterization of 19 Hippo core genes in 9,125 tumor samples across 33 cancer types using multidimensional "omic" data from The Cancer Genome Atlas. We identify somatic drivers among Hippo genes and the related microRNA (miRNA) regulators, and using functional genomic approaches, we experimentally characterize YAP and TAZ mutation effects and miR-590 and miR-200a regulation for TAZ. Hippo pathway activity is best characterized by a YAP/TAZ transcriptional target signature of 22 genes, which shows robust prognostic power across cancer types. Our elastic-net integrated modeling further reveals cancer-type-specific pathway regulators and associated cancer drivers. Our results highlight the importance of Hippo signaling in squamous cell cancers, characterized by frequent amplification of YAP/TAZ, high expression heterogeneity, and significant prognostic patterns. This study represents a systems-biology approach to characterizing key cancer signaling pathways in the post-genomic era.

INTRODUCTION

The Hippo signaling pathway is an evolutionarily conserved pathway that controls organ size and cell differentiation across diverse organisms through the regulation of cell proliferation and apoptosis (Di Cara et al., 2015; Pan, 2010; Yu and Guan, 2013; Zhao et al., 2007). The core of the Hippo pathway is a kinase cascade consisting of STE20-like protein kinase 1 (STK3, also known as MST2, and STK4, also known as MST1), the large tumor suppressors (LATS1 and LATS2), and adaptor proteins Salvador homolog 1 (SAV1) and MOB kinase activators (MOB1A/MOB1B) (Tapon et al., 2002). These kinases directly phosphorylate the major Hippo pathway downstream effectors, the yes-associated protein (YAP) and the transcriptional coactivator with PDZ-binding motif (TAZ), and thus inhibit the transcription of their downstream target genes (Cho et al., 2006; Hansen et al., 2015; Johnson and Halder, 2014; Oka et al., 2008). In recent years, the critical role of the Hippo pathway in cancer development has been increasingly recognized (Bonilla et al., 2016; Bueno et al., 2016; Gao et al., 2014; Yu et al., 2015; Zanconato et al., 2016). Dysregulated signaling by the Hippo pathway has been reported in several cancer types such as breast, liver, lung, prostate, gastric, and colorectal tumors (Cordenonsi et al., 2011; Harvey et al., 2013; Jiao et al., 2014; Lau et al., 2014; Nguyen et al., 2015; Yimlamai et al., 2014; Zhang et al., 2017; Zhou et al., 2009). Extensive studies have experimentally established its tumor suppressor function. For example, STK3/STK4 loss in a mouse liver leads to uncontrolled cell proliferation and differentiation (Zhou et al., 2009), and YAP/TAZ overexpression leads to tissue overgrowth and cancer (Camargo et al., 2007; Dong et al., 2007).

Previous studies on the Hippo pathway have focused on the effect of its individual components in a small set of tumor contexts, which provides limited or biased views of this important pathway. So far, a comprehensive molecular portrait of the Hippo pathway in cancer has not been characterized, resulting in important knowledge gaps for utilizing this pathway in cancer medicine. First, given the tremendous molecular diversity between cancer types, the relative importance of the Hippo pathway in different tumor contexts remains poorly understood. Second, because the Hippo pathway receives upstream signaling from multiple sources and involves many components, it remains unclear how to effectively characterize Hippo pathway activity as a molecular signature that can confer some clinical utility. Third, molecular drivers and regulators that affect Hippo pathway signaling in cancer development have not been fully characterized. Combining the bioinformatics analysis of multidimensional molecular profiling data from The Cancer Genome Atlas (TCGA), functional perturbation approaches, and cancer-type-specific predictive modeling, we focused on the YAP/TAZ-centered Hippo pathway (19 core genes) and aimed to address these critical questions in a systematic way.

RESULTS

Somatic Alteration Landscape of the Hippo Pathway

Curated from the literature, we defined 19 genes that function mainly through the Hippo pathway as Hippo core genes (Figure 1A). Focusing on these genes, we calculated the somatic copy number alteration (SCNA) and mutation frequency in the pan-cancer cohort of 9,125 patients (Figure 1B; Tables S1 and S2). The overall DNA aberration level was low, ranging from 1% to 5%. STK3 and TAZ showed the highest amplification frequency, followed by TEAD4, YAP1, and STK4. Deep deletions occurred mostly in LATS1/2, consistent with their tumor suppressor role in cancer development. To identify somatic driver candidates among Hippo genes, we identified background-corrected significantly recurrent SCNAs and significantly mutated genes (SMGs) in each cancer type using GISTIC2 (Mermel et al., 2011) and MutSigCV (Lawrence et al., 2013), respectively (q < 0.25; Figures 1C and 1D). For oncogenes, the most significant amplification peak was an amplicon of 11q22.1 in cervical squamous cell carcinoma (CESC), driven by the amplification of YAP1 (Figure S1A). For tumor suppressors, the most significant deletion peak was 17p in sarcoma (SARC), where TAOK1 resides (Figure S1B). Regarding the mutational profile, NF2 (23.2%) and LATS2 (9.8%) showed the highest mutation frequencies in mesothelioma (MESO). NF2 showed a striking pattern in MESO: all of the mutations were truncating mutations and led to reduced protein expression, indicating its loss-of-function effects (Figure S1C). These results showed that NF2 loss is a major cancer driver in silencing the Hippo pathway in this disease (Bianchi et al., 1995; Li et al., 2014; Sekido et al., 1995).

Figure 1.

Somatic Alteration Landscape of the Hippo Pathway

(A) Diagram of 19 Hippo pathway core genes. Red depicts oncogene, and blue depicts tumor suppressor.

(B) Waterfall plots of gene mutation and copy number alteration of 19 Hippo core genes. Each row represents a gene, and each column represents a sample. Genes are ranked from high to low somatic alteration frequency. Oncogenes are highlighted in red.

(C) Significant amplification peaks of oncogenes and deletion peaks of tumor suppressors in each cancer type identified by GISTIC2 (q < 0.25). Dot size shows level of significance; color depicts peak status (red: amplification peak; blue: deletion peak). Oncogenes are highlighted in red.

(D) Mutation frequency heatmap of Hippo pathway core genes in each cancer type. Color depicts mutation frequency, with higher mutation frequency in a cancer type indicated by darker color, and significantly mutated genes identified by MutSigCV (q < 0.25) are highlighted in black boxes.

(E) Hippo pathway DNA alteration load relative to the genomic background of each cancer type. The bar plots show the top percentile of the observed mutation load in the distribution of 1,000 random sample sets with the same number and same gene length of Hippo pathway genes.

See also Figures S1 and S2 and Tables S1 and S2.

Several cancer types showed little somatic alterations on the Hippo pathway compared with other cancers, such as pheochromocytoma and paraganglioma (PCPG), acute myeloid leukemia (LAML), and diffuse large B cell lymphoma (DLBC). To assess the cancer types in which Hippo pathway alterations are more (or less) prevalent than background expectation, we employed a random sampling approach (STAR Methods). Given the background DNA aberration frequencies in a specific cancer type, we found that the Hippo pathway was hyper-altered in five cancer types (top 5% among 1,000 permutations) and hypoaltered in LAML (bottom 5%) (Figure 1E; Figure S2). Among the hyper-altered cancer types, MESO, kidney renal papillary cell carcinoma (KIRP), and CESC were the most significant, which was probably due to the high frequency of NF2, LATS2, and SAV1 mutations and YAP1 amplification. These results reveal a highly heterogeneous somatic alternation landscape across 33 cancer types, suggesting the relative importance of Hippo pathway dysregulation in different tumor contexts.

YAP1/TAZ Somatic Alterations and Their Functional Effects

YAP1 and TAZ are the major downstream effectors of the Hippo pathway. To understand the SCNA patterns of YAP1 and TAZ across 33 cancer types, we first calculated their combined amplification frequency, which ranged from 0% to 19%. Interestingly, the top six cancer types with the highest amplification frequency included all five squamous cell-involved cancers (CESC, lung squamous cell carcinoma [LUSC], esophageal squamous cell carcinoma [ESCA], head and neck squamous cell carcinoma [HNSC], and bladder urothelial carcinoma [BLCA]), suggesting a driving role in squamous cell cancers (Figure 2A). Next, focusing on these squamous cell cancer types, we examined the mutually exclusive pattern of YAP1 and TAZ. Interestingly, the high-level amplifications of these two genes occurred in a mutually exclusive pattern in two squamous cell cancers (HNSC and CESC), collectively contributing to 14% and 19% of the tumor samples (Figure 2B). The patterns suggest functional complementarity of YAP1 and TAZ in activating down-stream transcriptional targets.

Figure 2.

YAP1/TAZ Amplification and Mutation Patterns

(A) The combined amplification frequency of YAP1/TAZ in different cancer types. Squamous cell cancer types are labeled in red.

(B) Significant mutually exclusive pattern of YAP1 and TAZ in squamous cell cancer. Each bar represents a gene, and each column is the copy number status of a patient sample. Red depicts high-level amplification, pink depicts low-level amplification, light blue depicts heterozygous deletion, and dark blue depicts homozygous deletion. p values of each pair are shown.

(C) Functional annotation of YAP1 somatic mutations based on cell viability assays. Red depicts moderately activating, orange depicts weakly activating, yellow depicts no difference from wild-type, and blue depicts inactivating. The relative mutation positions are shown at the top.

(D) Functional annotation of TAZ somatic mutations.

(E) Cell images of three functional YAP1 mutants compared with wild-type in 4 (upper row) and 11 days (lower row).

(F) Cell images of three functional TAZ mutants compared with wild-type in 4 (upper row) and 11 days (lower row).

YAP1/TAZ mutations have been rarely studied because of their low frequencies. To understand the impact of YAP1/TAZ somatic mutations, we experimentally assessed the effect of each of 22 YAP1 and 19 TAZ non-silent mutations observed in the pan-cancer cohort using moderate-throughput, sensitive cell viability assays (Figures 2C and 2D; STAR Methods). We transfected wild-type genes and mutants into MCF10A cells and measured their effects. Upon overexpressing the wild-type genes, we confirmed their stimulating effects on cell viability. Interestingly, compared with the wild-type, we identified three activating and three inactivating mutations in YAP1 and three activating and two inactivating mutations in TAZ. All of the identified inactivating mutations were truncating mutations, consistent with the finding that truncated proteins mainly lose their functions. To confirm our results, we further measured the cell numbers of selected activating and inactivating mutations, which showed the same effects (Figures 2E and 2F). These results provide large-scale experimental evidence of YAP1/TAZ mutation effects, highlighting an underappreciated functional diversity of somatic mutations observed in human cancers.

The Pan-Cancer Expression Patterns of the Hippo Pathway

To unravel the expression pattern of Hippo pathway genes, we first performed an analysis of mRNA expression in 15 cancer types with a sufficient number of matched tumor and normal samples. Among them, kidney renal clear cell carcinoma (KIRC), kidney chromophobe (KICH), LUSC, lung adenocarcinoma (LUAD), and thyroid cancer (THCA) had the largest numbers of differentially expressed Hippo genes (Figure S3A). The dysregulated pattern of individual genes varied greatly between cancer types, and the most consistently inhibited genes were the tumor suppressors LATS1/2 (Figure S3B). These results reveal perturbed Hippo pathway activity in cancer patients.

To gain a global view of Hippo gene expression across cancer types, we next performed an unsupervised consensus clustering of 9,125 samples of 33 cancer types based on the mRNA expression data. In addition to the 19 core genes, we included another 31 genes involved in regulating Hippo pathway signaling in order to obtain robust clustering patterns. Based on 50 Hippo-related genes, we identified five robust clusters (Figure 3A; STAR Methods). As expected, the samples in most cancer types almost exclusively fell into one cluster, suggesting the dominance of tissue effects (Figures 3A and 3B). Strikingly, only five squamous cell cancers (BLCA, CESC, ESCA, HNSC, and LUSC) were split into two clusters (clusters 1 and 3), whereas other tumor types resided mainly in a single cluster. Thus, squamous cell cancer types appear to have very distinct Hippo pathway activities that characterize the molecular diversity of these diseases.

Figure 3.

Gene Expression Pattern of Hippo Pathway Genes

(A) Unsupervised consensus clustering of Hippo component gene expression showing five distinct clusters, each marked by a different color in the top row bar.

(B) Sample distribution in the five clusters. Each row is a cancer type, and each column represents a cluster. The number in each cell shows the number of samples classified in the corresponding cluster; red intensity shows the percentage of samples falling into one cluster.

(C) Comparison of clustering pattern of Hippo pathway with the other nine major biological pathways. y axis shows the fraction of the second largest cluster in a given cancertype. The middle line in the box is the median, and the bottom and top of the box are the first and third quartiles, and the whiskers extend to 1.5× interquartile range of the lower quartile and the upper quartile, respectively. Only the Hippo pathway shows a statistically significant higher ratio for squamous cell-involved cancers than for non-squamous cell cancers, indicating the uniqueness of the observed clustering pattern of Hippo pathway. Wilcoxon rank-sum test p value and q are shown for significant comparisons.

See also Figure S3.

To determine whether the clustering pattern of squamous cell-involved cancer types is specific to the Hippo pathway, we performed the same consensus clustering analysis for another nine major biological pathways. We used the fraction of the second largest cluster as an index to quantitatively characterize the two-main-cluster pattern of squamous cell cancers (STAR Methods). Among all of the pathways surveyed, this index was significantly higher for squamous cell-involved cancer types than non-squamous cell cancers in the Hippo pathway only (p = 0.005; q = 0.05; Figure 3C). These results indicate that the high expression heterogeneity of Hippo-related genes is unique, rather than a common property of squamous cell cancer types.

miRNA Regulation of Hippo Pathway Genes

To delineate microRNA (miRNA)-mediated regulation, we integrated the target sequence information and co-expression pattern of mRNA-miRNA to identify potential major regulators for the Hippo pathway (Figure 4A). With a stringent selection criterion, we identified a total of 69 high-confidence, recurrent mRNA-miRNA regulatory pairs between 5 Hippo core genes and 18 miRNAs across 9 cancer types. The most frequently identified miRNA regulator (32 out of 69) was the miR-200 family (hsa-miR-200a, hsa-miR-200b, hsa-miR-200c, hsa-miR-141, and hsa-miR-429). They were predicted to regulate FRMD6, LATS2, and TAZ. The miR-200 family is best known for its regulatory role in epithelial-mesenchymal transition (EMT), and it has been reported to inhibit YAP1 in human breast cancer (Yu et al., 2013). In contrast, for the other Hippo major effector TAZ, few miRNAs have been reported with strong experimental evidence. To validate the predicted miRNA regulators of TAZ, we focused on miR-200a-3p and miR-590-3p, each of which contains at least two perfect binding sites in the 3′ UTR of TAZ (Figure 4B). Upon transfection of miR-200a-3p and miR-590-3p mimics in two colon cancer cell lines (HCT116 and KM12), we found that the two miRNAs significantly suppressed TAZ expression at both mRNA and protein levels (Figures 4C–4F). Together, these results not only present a systematic view of miRNA regulators for Hippo pathway genes but also suggest that the related miRNA analogs might be exploited clinically.

Figure 4.

miRNA Regulation of the Hippo Pathway

(A) High-confidence miRNA-Hippo pathway gene interaction network. miRNA-Hippo gene regulation is identified by integrating sequence information and co-expression pattern. Significant pairs recurring in at least three cancer types are shown (q < 10−5 and Rs < −0.5).

(B) TAZ 3′ UTR binding sites with miR-200a-3p and miR-590-3p. Red depicts the binding regions for miR-590-3p; purple depicts the binding regions for miR-200a-3p. Detailed sequences are magnified for each binding motif.

(C) Normalized TAZ gene expression level after miR-590-3p transfection in HCT116 and KM12 cell lines (**p < 0.01; ***p < 0.001).

(D) Normalized TAZ gene expression level after miR-200a-3p transfection in HCT116 and KM12 cell lines. (C and D) Data are represented as mean ± SD.

(E) Relative protein expression level of TAZ upon overexpression of miR-590-3p in HCT116 and KM12 cell lines. Western blot gel (upper panel) and quantification (lower panel) are shown.

(F) Relative protein expression level of TAZ upon overexpression of miR-200a-3p in HCT116 and KM12 cell lines. Western blot gel (upper panel) and quantification (lower panel) are shown.

Robust Prognostic Power of a YAP/TAZ Target Gene Signature

To assess the clinical impact of a dysregulated Hippo pathway in cancer, we sought to examine the correlation of its activity with patients’ overall survival times. We considered the following indexes: (1) TAZ mRNA expression; (2) TAZ protein expression; (3) YAP1 mRNA expression; (4) phosphorylated YAP protein expression; (5) the reverse-phase protein array (RPPA)-based Hippo pathway score (containing three protein markers); and (6) the mRNA-based Hippo pathway score (containing 19 core genes). In addition, because the Hippo pathway receives upstream signaling from multiple sources and its components undergo extensive post-transcriptional and post-translational regulations, the expression of pathway genes may not accurately capture the pathway activity. An alternative approach is to evaluate the consequences of the Hippo pathway: the expression of the YAP/TAZ direct downstream target genes. To achieve this, we manually curated 22 downstream target genes of the Hippo pathway from the literature and experimental data, and included the expression of this gene signature as another index (Figure 5A; Figure S4; STAR Methods). Based on independent cell line data, the selected 22 target genes showed high-level expression changes upon YAP/TAZ perturbations and strong positive correlations with YAP and TAZ protein expression (Figure 5A; Figure S4).

Figure 5.

Prognostic Power of the YAP/TAZ Target Gene Signature

(A) Flow diagram of the curation of the YAP/TAZ target score. Venn diagram of shared genes from published YAP/TAZ chromatin immunoprecipitation (ChIP)-sequencing data (Galli et al., 2015; Stein et al., 2015; Zanconato et al., 2015). Gene list was then filtered for significantly upregulated or downregulated genes (≥2-fold) in YAP/TAZ knockdown or YAP overexpression or activation transcriptomic data (Elbediwy et al., 2016; Galli et al., 2015; Kim et al., 2015; Stein et al., 2015; Zanconato et al., 2015; Zhang et al., 2008), following selection by significant correlations with YAP and TAZ protein expression in 891 cell lines from the Cancer Cell Line Encyclopedia (CCLE; p < 10−4) and manual curation leading to a 22-gene score.

(B) Summary of correlations of (1) mRNA-based 22-gene YAP/TAZ target score; (2) mRNA-based 22-gene YAP/TAZ target score after adjusting for tumor purity; (3) YAP1 mRNA; (4) YAP protein; (5) TAZ mRNA; (6) TAZ protein; (7) protein-based Hippo pathway gene score; and (8) mRNA-based Hippo pathway gene score with patient overall survival in different cancer types. Circle size represents statistical significance; color represents direction. Overall, the YAP/TAZ target score shows the most consistent and extensive significant correlations with patient survival times, thereby representing the best indicator for Hippo pathway activity.

(C) Kaplan-Meier survival plots of patients grouped by the YAP/TAZ target score in individual cancer types. Both p values for log rank and Cox tests are shown.

See also Figure S4.

For the individual gene expressions of TAZ and YAP1, we observed compatible patterns with their assumed functional roles as oncogenes in multiple cancer types, but occasionally observed the opposite patterns. For the mRNA-based Hippo pathway score, we observed five significant but mixed correlations (three versus two); for the RPPA-based Hippo pathway score, we observed more consistent positive correlations with patient survival in six cancer types (adrenocortical carcinoma [ACC], glioblastoma [GBM], LUSC, KIRP, brain lower grade glioma [LGG], and KIRC), consistent with the tumor suppressor role of the pathway. The most consistent and extensive correlations came from the YAP/TAZ target score based on the integrated mRNA expression of their 22 target genes (Figure 5A; STAR Methods). We found a significant negative correlation between target gene expression and overall patient survival times in nine cancer types (Figure 5B). Indeed, the survival patterns were quite striking in multiple cancer types (Figure 5C). Notably, four out of the five squamous cell cancers (BLCA, LUSC, HNSC, and CESC) showed decreased overall patient survival time with increased YAP/TAZ transcriptional activity. In addition, the significant survival patterns remained robust after adjusting for tumor purity in the analysis: two more cancer types, ovarian serous cystadenocarcinoma (OV) and CESC, showed significant correlations with survival time, whereas MESO lost significance mainly because of incomplete tumor purity data. These results indicate that the 22-gene target signature we defined represents a robust index that can effectively capture Hippo pathway activity in cancer development, and suggest its potential utility as a prognostic marker.

Biological Pathways and Infiltrating Immune Cell Compositions Associated with Hippo Signaling

Given the robust prognostic power of the YAP/TAZ target score across cancer types, we next focused on this index as a functional readout of Hippo pathway activity to understand the associated biological pathways. We first performed gene set enrichment analysis (GSEA) (Subramanian et al., 2005) based on the rank of mRNA expression levels in correlation with the target score. Quite a few hallmark pathways showed consistent correlations across all cancer types: the positively correlated pathways included EMT, tumor necrosis factor alpha (TNF-α) signaling, inflammatory response, angiogenesis, and more, whereas the negatively correlated pathways included fatty amino acid metabolism, DNA repair, MYC targets, and oxidative phosphorylation (Figure 6A). In parallel, we examined the correlations of the YAP/TAZ target score with 11 major signaling pathway scores defined by RPPA data in 33 cancer types, and observed highly consistent patterns (Figure S5). Cross-talk among several pathways has been reported and discussed, such as between the EMT pathway and nuclear factor-κB (NF-κB) signaling (Cho et al., 2010; Lei et al., 2008; Overholtzer et al., 2006). Therefore, it is worth looking at other networks that may indirectly contribute to dysregulation of the Hippo pathway in human cancers.

Figure 6.

Biological Pathways and Tumor-Infiltrating Immune Cell Abundance Associated with the YAP/TAZ Target Signature

(A) Bar plots showing the most consistent hallmark pathways correlated with the YAP/TAZ target score (q < 0.05). Color indicates the correlation direction.

(B) Heatmap showing the correlation of significant tumor-infiltrating immune cell abundance with YAP/TAZ target scores (q < 0.05 and |corr| > 0.3). Red depicts positive correlations, and blue depicts negative correlations.

See also Figure S5.

The Hippo pathway is primarily known for its role in cell proliferation and differentiation, but emerging evidence implicates its importance in cancer immunity (Moroishi et al., 2016; Wang et al., 2016). To obtain an overview of the potential interactions between Hippo pathway activity and tumor immunogenicity, we calculated the correlations of YAP/TAZ target scores with tumor-infiltrating cell abundance (inferred by TIMER), including B cells, CD8 T cells, CD4 T cells, neutrophil, macrophage, and dendritic cells (Li et al., 2016). We observed strong positive correlations between all of these types of tumor-infiltrating cells with YAP/TAZ activity in the vast majority of cancer types except for testicular germ cell tumors (TGCTs) and thymomas (THYMs) (Figure 6B). This pattern is concordant with previous findings that activated YAP/TAZ signaling promotes immune cell recruitment (Moroishi et al., 2016; Wang et al., 2016). The strong correlations observed in this comprehensive analysis also highlight the importance of the Hippo pathway in the tumor microenvironment and its potential applications in cancer immunotherapy.

Cancer-Type-Specific Hippo Pathway Regulation Revealed by Integrative Modeling

To elucidate the key Hippo pathway components and cancer drivers affecting or associated with the pathway activity in various cancer types, we developed an elastic-net-based machine learning approach to build predictive models for the YAP/TAZ target score (Figure 7; STAR Methods). For each cancer type, we incorporated somatic mutation, SCNA, DNA methylation, mRNA expression and RPPA-based protein expression, and miRNA regulators, combined with cancer-type-specific somatic drivers (e.g., SMGs and recurrent SCNAs), and selected the top weighted features to construct a cancer-type-specific Hippo pathway regulatory network (Figure 7A). These networks highlight the diverse regulatory mechanisms in different cancer types (Figure 7B; Figures S6 and S7). Among Hippo genes, the most important regulators were mRNA expression levels of YAP1, TAZ, and TEADs. Specifically, elevated YAP/TAZ expression levels caused by SCNAs were selected with strong preference in our models of squamous cell-involved cancers (Figure 7B), supporting YAP/TAZ amplification as a driving force in these cancer types. Decreased NF2 expression level caused by NF2 truncating mutation was a dominant regulator of YAP/TAZ activity in MESO. Hyper-methylation of STK3, WWC1, and TAOK2 was predicted to strongly promote YAP/TAZ amplification in colorectal cancer, MESO, and BLCA, respectively. Furthermore, copy number loss for the TAOK family was a major cause of elevated YAP/TAZ activity in KIRP and LUSC. As for miRNA regulators, the miR-200 family was the most frequently identified negative regulator in seven out of the nine cancer types with a prognostic YAP/TAZ target score (Figure 5C).

Figure 7.

Integrative Modeling of Cancer-Type-Specific Hippo Regulatory Networks

(A) Workflow of integration modeling based on elastic net regression. For each cancertype, somatic mutation, SCNA, mRNA expression, DNA methylation, and RPPA protein expression data of 19 Hippo core genes, and cancer-specific predicted miRNA regulators were used as features to build an elastic net regression model to predict YAP/TAZ target score. In parallel, the modeling for cancer-type-specific significantly mutated genes and SCNAs were carried out. After 100 rounds of 10-fold cross-validation, the best model yielding the lowest mean squared error was selected, and the top five weighted predictors were used to construct the Hippo regulatory network.

(B) Hippo regulatory networks of the nine cancer types with YAP/TAZ target scores that showed significant survival correlations. Each color depicts a molecular data type, and each shape depicts a molecular type. Blue line depicts negative regulations, and red line depicts positive regulations.

See also Figures S6 and S7.

For somatic drivers associated with Hippo pathway activity, the PTEN copy number was negatively associated with YAP/TAZ activity in BLCA, stomach adenocarcinoma, LUSC, and HNSC. Other recurrently identified somatic driver events included ERBB2, RB1, and SOX2 copy number changes. Besides SCNAs, IDH1 and FGFR3 mutations were predicted to be associated with YAP/TAZ activity in LGG and BLCA, respectively. Other identified associated somatic driver events included EGFR, KRAS, and PDGFR copy number changes, consistent with the reported functional coupling of these pathways with Hippo pathway signaling (Chaib et al., 2017; Kapoor et al., 2014; Lin et al., 2015; Shao et al., 2014). We note that the cancer driver events identified may not directly affect Hippo pathway signaling, but rather simply reflect the tumor contexts in which the Hippo pathway tends to be perturbed; further efforts are required to elucidate the mechanisms underlying the observed associations. Importantly, we found that Hippo pathway activity was associated with aberrations in a few clinically actionable genes with Food and Drug Administration (FDA)-approved drugs or compelling clinical evidence (Chakravarty et al., 2017), which provides insights into clinical applications targeting the Hippo pathway in cancer therapy.

DISCUSSION

Using the latest TCGA multi-dimensional molecular profiling data, we performed a comprehensive molecular characterization of the Hippo pathway across >9,000 samples of 33 cancer types. In addition to computational analyses, we have experimentally assessed the effects of YAP1/TAZ somatic mutations using cellular assays. In total, we identified 11 non-silent mutations that affect cell viability and proliferation, filling a critical knowledge gap about YAP1/TAZ somatic mutations. We inferred miRNA regulation by considering both sequence information and expression pattern, and validated two TAZ-negative regulators, miR-200a and miR-590. This combined computational and experimental approach represents the direction for characterizing key cancer pathways in the post-genome era.

Our analysis reveals that diverse mechanisms drive pathway dysregulation in tumor contexts. We not only confirm the mutation-driven mechanisms in MESO as previously reported, but also highlight the critical role of Hippo pathway signaling in squamous cell cancers. Through our integrated analyses of BLCA, CESC, ESCA, LUSC, and HNSC, we reveal frequent YAP1/TAZ amplifications, high expression heterogeneity, and significant prognostic correlations. Given the strong selective pressure to alter the Hippo pathway and activate YAP/TAZ in these cancers, Hippo pathway signaling appears to play a major role in shaping tumor progression and may represent a promising target for future development of novel therapies. Another cancer of particular interest is LGG, for which Hippo pathway activity shows the strongest correlation with patient survival times and the most extensive associations with known drivers such as IDH1 mutations. However, this disease shows low-frequency aberrations among Hippo pathway genes in terms of both mutations and SCNAs, suggesting a role of unknown epigenetic mechanisms controlling Hippo pathway signaling.

Another significant contribution of this study is that we developed a YAP/TAZ transcriptional target signature of 22 genes and demonstrated the superior performance of the expression score based on this signature over that of individual pathway genes in predicting patient prognosis. This represents an alternative strategy for robustly inferring pathway activity: focus on the genes affected downstream of a pathway rather than the molecular aberrations within the pathway itself. We suspect that similar target gene signatures can be developed for other major cancer signaling pathways and used for pathway activity inference. The robust prognostic power and the relatively small number of genes involved also make this gene signature highly valuable for potential clinical applications.

We systematically constructed cancer-type-specific Hippo regulatory networks by integrating various types of molecular data. The associations between somatic driver events with Hippo pathway signaling may coordinate the cross-talk between the Hippo pathway and other signaling pathways. In the modeling results, we identified the miR-200 family as master regulators for YAP/TAZ activity. The miR-200 family is best known for its regulatory role in the EMT pathway. Although it has been demonstrated that miR-200 targets YAP1 in breast cancer (Yu et al., 2013), its functional impact on the Hippo pathway remains largely unexplored and is worth further exploration in the future. Taken together, our focused and systematic analyses of this important molecular pathway will serve as a valuable resource for understanding its dysregulation in cancers and how to maximize its clinical utility.

STAR★METHODS

CONTACT FOR REAGENT AND RESOURCE SHARING

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Han Liang (HLiang1@mdanderson.org).

EXPERIMENTAL MODEL AND SUBJECT DETAILS

LentiX-293T cells (Clontech) were cultured in DMEM (with high glucose, glutamine and sodium pyruvate) with 5% fetal bovine serum and 1 × non-essential amino acid. MCF10A cells (ATCC® CRL-10317) were cultured in specific growth medium (DMEM/F12 medium with 5% HS, 20 ng/ml epidermal growth factor, 0.5 mg/ml hydrocortisone, 100 ng/ml cholera toxin, 10 μg/ml insulin, 1 × Pen/strep). Assay medium for MCF10A cells was MEBM basal medium (Lonza #CC-3151) with 100 ng/ml cholera toxin and 52 ng/ml bovine pituitary extract (Lonza #CC-4009).

HCT116 and KM12 colon carcinoma cell lines were purchased from MD Anderson Characterized Cell Line Core facility. Cell lines were confirmed using Short Tandem Repeat (STR) analysis and were tested and negative for mycoplasma contamination. Cells were cultured in RPMI 1640 medium supplemented with 10% fetal bovine serum. MiR-590-3p and miR-200a-3p mimics, and MISSION miRNA Negative Control 1 (HMC0002) were purchased from Sigma Aldrich. Cells were transfected with 50nM of miRNA mimics or miRNA Negative Control using Lipofectamine RNAiMAX reagent (Thermo Fisher).

METHOD DETAILS

Somatic Copy-Number Alteration (SCNA) Analysis

We obtained TCGA thresholded SCNA scores for 9,125 patient samples from Genome Data Commons. Values equal to 2 and −2 were considered to indicate high-level amplification and homozygous deletion, respectively. Based on the copy number segmentation file, we identified significantly amplified and deleted genes (q < 0.25) in each cancer type using GISTIC2.0 with the default parameters. We performed mutual exclusivity analysis using the CoMET module (Leiserson et al., 2015) in R for YAP1 and TAZ in squamous-cell-involved cancers, and reported significant pairs with p < 0.05.

Somatic Mutation Analysis

We obtained somatic mutation data from Genome Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). Further filtering steps were used to eliminate artifacts and prevent false-positive calls. i) Only mutations with “PASS” in the “FILTER” column were retained for all cancer types except OV and LAML, and we allowed “wga” for these two cancer types. ii) Hypermutated samples with > 1,000 somatic mutations were removed, resulting in a cohort of 8,811 patient samples for downstream analyses. We identified significantly mutated genes using MutSigCV 1.4 with a q-value cutoff of 0.25. To assess the Hippo pathway genetic alteration load, we considered the tumor context and gene length by randomly selecting gene sets with similar length distributions in each cancer type. First, we divided all protein coding genes into bins based on their gene lengths. Then, for each sampling iteration, we randomly selected 19 genes from their corresponding bins to generate a random gene set with a similar length distribution as the Hippo components. We calculated the frequency of somatic mutation and copy number changes of this randomly sampled gene set. Then, by repeating the sampling 1,000 times, we calculated the rank of the observed Hippo gene DNA aberration frequency given the random sampling background.

Functional Assays of YAP1/TAZ Mutations

Functions of YAP1 and TAZ mutations were assayed in MCF10A cells, human non-tumorigenic mammary epithelial cells. These cells depend on exogenous epidermal growth factor and insulin for proliferation. We hypothesized that a “driver” mutation confers proliferation advantages to the cells in the absence of required growth factor(s), but a “passenger” mutation does not. The mutations, wild-types and silent controls of YAP1 (YAP) and TAZ (WWTR1) genes were cloned into pHAGE-PURO lentivirus vector by HiTMMoB technique as described previously (Dogruluk et al., 2015; Tsang et al., 2016). Lentivirus was generated in the LentiX-293T cells by transfecting the pHAGE and two packaging plasmids (psPAX2 and pMD2.G). The medium of the transfected cells was refreshed at 16 hours post-transfection. The virus was harvested at 3 days post-transfection by filtering with 0.45 μM filter and used for transduction of MCF10A cells. For cell viability assay, 5,000 cells were seeded in 96-well plates in triplicate a day before transduction and transduced by spinoculation at 906 × g for 2 hours in the presence of polybrene (final concentration: 2.7 μg/ml). The medium was refreshed after spinoculation with the MCF10A assay medium. Transduced cells were incubated at 37°C for 2 weeks. Cell viability of MCF10A cells was performed using CellTiter-Glo (Promega) on 4, 8, 11 and 15 days post-transduction. The functional annotations of mutations were made based on a comparison with the corresponding wild-type clones. For cell count assays, 1,000 cells were seeded into 24-well plates in triplicate one day before transduction. The cells were transduced and the medium was refreshed as described above. An entire-well image of each well was taken by IncuCyte ZOOM® system on 4, 8, 11 and 15 days post-transduction.

Gene Expression Analysis

We obtained normalized gene expression data from Genome Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). For comparison of tumor tissue to normal tissue, we used a paired t test to detect genes that were differentially expressed in tumor versus normal samples, and considered a false discovery rate (q) < 0.05 as statistically significant. All negative values were considered as missing values (NA), and we took the log-transformed values in all downstream analyses. For cancer types with a sufficient number of tumor-normal matched pairs (n ≥ 10), we used the paired t test to identify expression pattern differences. Unsupervised consensus clustering was performed on 9,125 patient samples across 33 cancer types based on the expression of 50 Hippo core or regulatory genes with less than 0.5% invalid values (NA or 0) using the Pearson distance and Ward’s method. For each cancer type, we used the Wilcoxon rank-sum test or Kruskal-Wallis nonparametric analysis of variance test to detect genes that were differentially expressed among the cancer subtypes, and considered a false discovery rate (q) < 0.05 as statistically significant. To compare the clustering pattern of the Hippo pathway and the other 9 major biological pathways, we used the gene list summarized by TCGA Pancan pathway group, and performed the same clustering method. We quantified the squamous-cell-involved separation pattern using the ratio of the second largest cluster size divided by the entire cohort sample number, and evaluated whether there is a statistical significant difference between squamous-cell-involved and non-squamous-cell cancers. We used the Wilcoxon rank-sum test, considering q < 0.05 as statistically significant.

miRNA-Hippo Pathway Interaction Prediction

We obtained miRNA sequencing data from Firehose (version 2016/1/28). We used two criteria to identify the miRNA regulators for Hippo pathway genes: i) the 3′UTR of a Hippo gene should contain at least one target site matched to the seed region (position 2-8) of the miRNA; and ii) the Spearman correlation between the miRNA and the gene expression levels should be significantly negative (q < 10−5 and Rs < −0.5). Only interaction pairs that recurred in ≥ 3 cancer types were considered as high-confidence miRNA-Hippo regulators.

miRNA Target Experimental Validation

RNA Isolation and Quantitative Real-Time PCR

miRNA expression was assessed 48 h after transfection. Total RNA was extracted using the mirVana miRNA Isolation Kit, with phenol (Ambion) and reverse transcribed using the TaqMan Advanced miRNA cDNA Synthesis Kit (Applied Biosystems). Mature miRNAs were detected and quantified using the TaqMan Fast Advanced Master Mix according to manufacturer’s instructions (Applied Biosystems). TaqMan miRNA assays for hsa-miR-590-3p (Assay ID 478168_mir), hsa-miR-200a-3p (Assay ID 478490_mir), and hsa-miR-26a-5p internal control (Assay ID 477995_mir) were purchased from Thermo Fisher. To quantify the expression of mRNA targets, RNA was extracted using the RNeasy kit (QIAGEN) and cDNA synthesis was carried out using High-Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (Applied Biosystems). Real-time PCR was performed using SYBR Select Master Mix (Applied Biosystems). Primer pairs (TAZ F-5′-AATAGCTCAGATCCTTTCCTC-3′; R-5′-TCTCCTGTATCCATCTCATC-3′; β-actin F-5′-ATTGG CAATGAGCGGTTC-3′; R-5′-CGTGGATGCCACAGGACT-3′) were purchased from Sigma Aldrich. TAZ mRNA level was normalized to beta-actin and reported as fold-change compared with the negative control.

Western Blot Experiments

Total protein was isolated from cells 72 h after transfection. Cells were lysed with RIPA buffer (25 mM Tris-HCl pH 7.6,150 mM NaCl, 1% NP-40, 1% sodium deoxycholate, 0.1% SDS) containing protease and phosphatase inhibitor cocktail (Thermo Fisher). Protein concentrations were determined using bicinchoninic acid assay (Pierce). Cell lysates (20ug) were resolved on 4%–20% SDS-PAGE gels (Bio-Rad) and transferred to polyvinylidene fluoride membrane. Membranes were probed with TAZ (Cell Signaling Technology, #4883) and GAPDH (Santa Cruz Biotechnology, sc-25778) antibodies, and visualized using a chemiluminescence western blot detection kit (GE Healthcare). Protein levels for TAZ were quantified using ImageJ and normalized to GAPDH.

RPPA Analysis

We obtained TCGA RPPA data from The Cancer Proteome Atlas (version 4.2). Among 19 Hippo core genes, three of them had corresponding RPPA data. Based on the predefined role of each gene, we calculated the Hippo pathway RPPA score for each patient using an equal-weighted summation of these three proteins: Hippo RPPA score = RPPANF2- RPPAYAP - RPPATAZ.

YAP/TAZ Target Score

We curated the 22-gene YAP/TAZ target signature (MYOF, AMOTL2, LATS2, CTGF, CYR61, ANKRD1, ASAP1, AXL, F3, IGFBP3, CRIM1, FJX1, FOXF2, GADD45A, CCDC80, NT5E, DOCK5, PTPN14, ARHGEF17, NUAK2, TGFB2, RBMS3) using published RNA-sequencing and ChIP-sequencing data from human cholangiocarcinoma, breast carcinoma, glioma, mesothelioma, and normal diploid fibroblast cell lines (Galli et al., 2015; Stein et al., 2015; Zanconato et al., 2015). We considered genes to be direct YAP/TAZ targets if their expression was robustly modulated upon YAP/TAZ knockdown (RNA-sequencing experiments) and their transcription start site or enhancers were simultaneously YAP- and TEAD-bound (ChIP-sequencing experiments). We excluded genes that were well-established targets of other signaling pathways (e.g., WNT, cell cycle, MAPK). We calculated the YAP/TAZ target score within each cancer type using TCGA RNA-sequencing data (summarizing the Z-normalized log2RSEM of the expression data for the 22 curated YAP/TAZ downstream transcription target genes). The workflow and number of genes filtered in each step is summarized in Figure 5A.

Patient Survival Analysis

We obtained the patients’ clinical data from Genome Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). We assessed the correlation between the YAP/TAZ target score and patient survival times. We used: i) the univariate Cox proportional hazards model to assess the correlation with the patient’s overall survival time; and ii) the log-rank test to compare patient survival curves between the groups with high- and low target score (separated by the median value), considering p < 0.05 as statistically significant. We obtained tumor purity data from Genome Data Commons (Carter et al., 2012), and repeated the above analysis by taking tumor purity as a Cox model covariate.

Biological Pathway and Cancer Immunity Analysis

We identified biological pathways associated with the YAP/TAZ score using the GSEA pre-ranked tool (Subramanian et al., 2005). Briefly, each gene was ranked based on its expression correlation coefficient with the YAP/TAZ target score within each cancer type, and the pre-ranked gene lists were then run against using GSEA Java (version 2.2.3). We reported the pathways that showed consistent significant correlations (q < 0.05) in ≥ 2/3 of all cancer types. We integrated the normalized protein expression data (Z-score) from TCGA RPPA platform into 11 core cellular pathways (Akbani et al., 2014). We calculated the Spearman rank correlations of the YAP/TAZ target score with these pathway scores, and considered q < 0.05 as statistically significant. We obtained tumor-infiltrating immune cell abundance from TIMER (Li et al., 2016), and for each cancer type calculated Spearman rank correlation between YAP/TAZ target score with each immune cell component score with q < 0.05 and |Rs| > 0.3.

Integrative Modeling of Hippo Pathway Regulatory Networks

We obtained DNA methylation data (450k, 27k) from Genome Data Commons. For each Hippo pathway gene, we first selected the probe with the most significant negative correlation with its mRNA expression in each cancer type. Genes with no significant correlated probes were removed from this analysis. We separated the modeling of Hippo core genes with somatic driver events to prevent potential mocking effects. First, we integrated the mutation, somatic copy number variation, gene expression, methylation, and RPPA protein expression data of 19 Hippo core genes, and cancer-specific predicted miRNA regulators to build the elastic net regression model. The most informative dimension of molecular data was preselected to represent the overall predictive power of each gene for YAP/TAZ activity. This was evaluated by the Spearman correlation of the molecular data with the YAP/TAZ target score, and the one with the strongest correlation according to its biological definition in the Hippo pathway was selected. Then, all of the selected 19 features and miRNA regulators were pulled together and trained in elastic net regression. Second, cancer-specific significant mutated genes and copy number alterations were used separately in the modeling process. Each regression analysis was repeated 100 times using 10-fold cross-validation. The best model yielding the lowest mean squared error was used, and the top 5 weighted predictors were combined to construct a Hippo regulatory network in each cancer type.

QUANTIFICATION AND STATISTICAL ANALYSIS

All the analyses were based on 9,125 tumor samples except for miRNA (7,939), and RPPA (6,441) due to limited data availability. Definition of significance of various statistical tests were described and referenced in their respective Method Details sections.

DATA AND SOFTWARE AVAILABILITY

The raw data, processed data and clinical data can be found at the legacy archive of the GDC (https://portal.gdc.cancer.gov/legacy-archive/search/f) and the PancanAtlas publication page (https://gdc.cancer.gov/about-data/publications/pancanatlas). The mutation data can be found here (https://gdc.cancer.gov/about-data/publications/mc3-2017). TCGA data can also be explored through the Broad Institute FireBrowse portal (http://gdac.broadinstitute.org) and the Memorial Sloan Kettering Cancer Center cBioPortal (http://www.cbioportal.org). Details for software availability are in the Key Resources Table.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited Data
TCGA somatic copy number alteration thresholded data	Genomic Data Commons	https://gdc.cancer.gov/about-data/publications/pancanatlas
TCGA somatic copy number segmentation data by Affymetrix SNP 6 array	Genomic Data Commons and Firehose	https://gdc.cancer.gov/about-data/publications/pancanatlas
TCGA somatic mutation data	Genomic Data Commons	https://gdc.cancer.gov/about-data/publications/pancanatlas
TCGA gene expression data	Genomic Data Commons	https://gdc.cancer.gov/about-data/publications/pancanatlas
TCGA reverse-phase protein array (RPPA) data	Genomic Data Commons	https://gdc.cancer.gov/about-data/publications/pancanatlas
TCGA DNA methylation data	Genomic Data Commons	https://gdc.cancer.gov/about-data/publications/pancanatlas
TCGA miRNA-seq data	TCGA TCGA Firehose	https://gdac.broadinstitute.org
TCGA patient clinic data	Genomic Data Commons	https://gdc.cancer.gov/about-data/publications/pancanatlas
Software and Algorithms
CoMET	(Leiserson et al., 2015)	https://bioconductor.org/packages/release/bioc/html/coMET.html
MutSigCV	(Lawrence et al., 2013)	http://software.broadinstitute.org/cancer/software/genepattern/modules/docs/MutSigCV
GISTIC2.0	(Mermel et al., 2011)	http://portals.broadinstitute.org/cgi-bin/cancer/publications/view/216
Gene Set Enrichment Analysis (GSEA)	(Mootha et al., 2003; Subramanian et al., 2005)	http://software.broadinstitute.org/gsea/index.jsp
Experimental Models: Cell Lines
LentiX-293T cells	Clontech	Cat#: 632180
MCF10A	ATCC	ATCC® CRL-10317; RRID: CVCL_0598
HCT116	MD Anderson Characterized Cell Line Core facility	HCT116
KM12	MD Anderson Characterized Cell Line Core facility	KM12
Chemicals, Peptides, and Recombinant Proteins
Recombinant human EGF	R&D Systems	Cat#: 236-EG-200
Hydrocortisone	Sigma-Aldrich	Cat#: H-0888
Cholera toxin	Sigma-Aldrich	Cat#: C-8052
Insulin	Sigma-Aldrich	Cat#: I-1882
MEGM Bullet Kit	Lonza	Cat#: CC-3150
polybrene	Sigma-Aldrich	Cat#: TR-1003-G
Lipofectamine RNAiMAX Transfection Reagent	Thermo Fisher	Cat#: 13778150
RPMI 1640	Corning, NY, USA	Cat#:10-040-CV
Fetal bovine serum	GIBCO	Cat#:16140-071
DMEM with 4.5 g/L glucose, L-glutamine, & sodium pyruvate	Corning, NY, USA	Cat#:10-013-CV
TAZ (V386) Antibody	Cell Signaling Technology	Cat#: 4883; RRID: AB_1904158
GAPDH Antibody	Santa Cruz Biotechnology	Cat#: sc-25778; RRID: AB_10167668
Oligonucleotides
miR-590-3p	Sigma-Aldrich	HMI0815
miR-200a-3p	Sigma-Aldrich	MIRAP00247
miRNA Negative Control I	Sigma-Aldrich	HMC0002
Critical Commercial Assays
CellTiter-Glo 2.0 Assay	Promega	Cat#: G9243
mirVana™ miRNA Isolation Kit, with phenol	Ambion	Cat#: AM1561
TaqMan™ Advanced miRNA cDNA Synthesis Kit	Applied Biosystems	Cat#: A28007
TaqMan™ Fast Advanced Master Mix	Applied Biosystems	Cat#: 4444557
High-Capacity cDNA Reverse Transcription Kit	Applied Biosystems	Cat#: 4374967
RNeasy Plus Mini Kit	QIAGEN	Cat#: 74136
SYBR™ Select Master Mix	Applied Biosystems	Cat#: 4472908