Machine-Learning and Chemicogenomics Approach Defines and Predicts Cross-Talk of Hippo and MAPK Pathways

One challenge of cancer precision medicine is the heterogeneity of genetic and nongenetic alterations that result in aberrant pathway signaling. Recurrent mutations and genetic alterations have been identified in many oncogenic signaling pathways, including MAPK and PI3K (1, 2), whereas other signaling pathways such as Hippo lack canonical hotspot mutations. Yet dysregulation in Hippo pathway signaling is known to drive oncogenesis across numerous cancer types. The Hippo pathway is emerging as the target of drug discovery efforts, but it lacks hotspot mutations; identifying relevant Hippo pathway–dependent patient population(s) and biomarker(s) of response is a prerequisite for precision medicine in tumors that leverage this pathway.

The Hippo pathway controls multiple cellular functions that drive oncogenesis, including proliferation, cell fate determination, and cell survival. Perturbation of the pathway has been shown to trigger tumorigenesis in mice (3). The pathway is evolutionarily conserved across diverse species and was first identified in Drosophila melanogaster through multiple genetic screens for gene mutations that cause overgrowth phenotype (4–6). These led to the discovery of the conserved Hippo pathway core components consisting of serine/threonine kinases named Mammalian STE20-like 1/2 (MST1/2) with adaptor protein SAV1 that directly phosphorylate the large tumor suppressors (LATS1/2). Together with the kinase activator MOB1, LATS1/2 can phosphorylate the two major downstream coactivators YAP (YAP1) and TAZ (WWTR1; Fig. 1A, inset). When the pathway is deregulated, unphosphorylated YAP and TAZ are translocated to the nucleus and activate downstream target gene expression by binding to TEAD family transcription factors (TF; refs. 7–13; Fig. 1A, inset). Widespread dysregulation of the Hippo pathway components has been observed in multiple human cancer types including glioma and breast, liver, lung, prostate, colorectal, and gastric cancers (14–17). Furthermore, tumors with dysregulated Hippo components are not only insensitive to the intrinsic cellular death barriers (3, 18) but also resistant to chemo and molecular targeted therapies (19–21).

Extensive studies have established the importance of the Hippo pathway in biology and cancers. As drug development interest in targeting the pathway continues to grow (22–25), one key clinical challenge is to identify patient populations that would benefit from such a therapy. Previous studies on the Hippo pathway have either defined broad genetic alterations in pathway component(s) or focused on individual cancer types or cancer cell lines. This experimental strategy has established the role of the Hippo pathway in cancers; however, regulation of Hippo pathway signaling can be highly complex with many linked signaling inputs from the orthogonal pathways (26). Here, we employ an integrated experimental–computational strategy to identify a lineage-independent signature for the Hippo pathway in cancers. By evaluating transcriptional profiles, we observed a relationship between YAP/TAZ dependency and MAPK pathway activity, leading us to nominate potential combination therapies with Hippo pathway inhibition.

To understand the role of the Hippo pathway in human cancers, we first examined the pathway alternations using data from The Cancer Genome Atlas (TCGA; ref. 27). Although genetic alterations in the Hippo pathway are infrequent (1%–15% across individual cancer types), YAP1 amplifications are among the most frequent alterations pan-cancer in the Hippo pathway (Fig. 1A) and most frequently observed in patients with cervical and head and neck squamous cell cancer. As expected, genetic YAP1 amplifications, but no other Hippo pathway alterations, were exclusively associated with YAP1 RNA overexpression in multiple cancer types (Fig. 1B). Furthermore, genetic YAP1 amplifications, along with alterations in other Hippo pathway members, were mutually exclusive across samples of patients with cancer (Fig. 1B), suggesting these low-frequency mutations may function similarly to deregulate the Hippo pathway.

As YAP is a transcriptional coactivator, the most frequently altered regulator of the Hippo pathway, and previously associated with treatment resistance (19, 20), we hypothesized that its oncogenic potential must be mediated by its downstream transcriptional target genes. We aimed to develop a first-principle approach to map a lineage-independent transcriptional signature for Hippo deregulation. We first identified seven cell lines originating from different tissues, but all carry YAP1 amplification (copy number, 6.29 ± 1.50) with marked YAP1 mRNA overexpression (Fig. 1C). We performed knockdown of YAP1 and its paralog WWTR1 and then performed RNA sequencing (RNA-seq) on the parental and YAP1/WWTR1 knockdown lines. YAP1/WWTR1 knockdown resulted in broad transcriptomic deregulation in all cell lines (Fig. 1C). Although CTGF expression (a canonical Hippo pathway target gene) was significantly decreased in all cell lines (Fig. 1C), there was no clear association between CTGF expression and magnitude of global gene expression changes to a cell line's sensitivity to YAP1/WWTR1 knockdown (Supplementary Fig. S1A and S1B). Taken together, this suggested that YAP/TAZ dependency may be more complex, which necessitates expanding beyond a single marker of pathway activity to capture the pathway dependency over different cell lineages.

We performed an unbiased weighted correlation network analysis (28, 29) among a consensus set of genes that were broadly expressed across all tissues in addition to being significantly and consistently downregulated (in at least three of seven cell lines) upon YAP1/WWTR1 knockdown (Fig. 1D). This identified four distinct gene clusters of coexpressed genes and one cluster of noncorrelated genes (Fig. 1E). Interestingly, we noted that many of the canonical Hippo pathway–regulated genes (e.g., CTGF, CYR61, etc.) were found within gene Cluster 2, suggesting that Cluster 2 may be most proximal to Hippo pathway signaling. Among the 145 genes in Cluster 2, 86% (n = 124) have not been reported in previous YAP/TAZ gene signatures (Supplementary Fig. S1C) in which we orthogonally validate several genes using RT-PCR (Supplementary Fig. S1D). To further validate Cluster 2, we leveraged recent systematic CRISPR and RNAi dependency screens (30). Although these datasets utilize only single-gene knockout, nevertheless, we performed gene-wide regression analysis with overlapping cell lines to assess whether the new gene set is associated with a given gene knockout/knockdown. Among the most significant gene dependencies, this analysis confirmed many Hippo pathway effectors included WWTR1, YAP1, and TEAD1 (Supplementary Fig. S1E and Supplementary Table S1). We then performed an unbiased analysis of somatic genetic predictors of Cluster 2 single-sample gene set enrichment analysis scores in the TCGA pan-cancer cohort (Supplementary Fig. S1F). The most significant results included NF2 loss-of-function mutations and homozygous deletions (Supplementary Fig. S1F and S1G), consistent with Hippo pathway regulation. Furthermore, we performed RNA-seq on three independent, NF2-null (Hippo pathway altered) cell lines [i.e., GOS-3 (glioma), MDA-MB-231 (triple-negative breast cancer), and MS751 (cervical)] after YAP1/WWTR1 knockdown. Consistent with the original seven YAP1-amplified cell lines, we observed similar numbers of overlapping, significantly downregulated genes in each of the three independent NF2-null cell lines (Supplementary Fig. S1H).

Last, we performed Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) on Detroit 562 and PA-TU-8902, a pancreatic adenocarcinoma cancer cell line with TEAD4 amplification, and confirmed that Cluster 2 genes were most strongly associated with loss of chromatin accessibility upon YAP1/WWTR1 knockdown (Fig. 1F; Supplementary Fig. S1I and S1J). Taken together, Cluster 2 genes included the most well-known canonical Hippo pathway marker genes, mostly correlated with previously reported Hippo pathway activity genes, and were associated with loss of chromatin accessibility upon knockdown.

Aberrant Hippo pathway signaling has been known to drive oncogenesis in several cancer types, many of which lack a known Hippo pathway genetic alteration. To better identify potential Hippo pathway–dependent populations, we sought to predict YAP/TAZ dependency using the cluster genes we have identified here. We identified and performed RNA-seq on a broader set of 42 cancer cell lines exhibiting a spectrum of Hippo pathway activity. Next, we assessed each cell line's sensitivity to YAP1/WWTR1 knockdown to train a machine-learning (ML) computational model to predict YAP/TAZ dependency given a cell line's parental cluster gene expression profile. The ensemble-based algorithm learned a combination of gene expression values to predict the change in viability after YAP1/WWTR1 knockdown (Fig. 2A). Cluster 2 score was the most correlated to predicted dependency, further supporting that Cluster 2 is most proximal to aberrant Hippo pathway signaling, the primary driver of YAP/TAZ dependency (Supplementary Fig. S2A and Supplementary Table S2). Given Cluster 2 included many novel genes not reported in previous gene sets, we next benchmarked our gene set compared with previously published gene sets (21, 31). We observed that our gene set performed better, independent of algorithm or training data (Supplementary Fig. S2B and S2C; Methods). Although we see that the known genes (e.g., CTGF, CYR61, etc.) have high importance/weight in the ML model (Supplementary Fig. S2D), many of the genes not found in previous gene sets (21, 31) were among the greatest importance/weight in the ML model's predictive power (Supplementary Fig. S2E and Supplementary Table S3), including CCDC42EP1, TNFRSF12A (32), and PHLDB2I (33).

Certain tissue lineages and histologic cell types were significantly associated with YAP/TAZ dependency, including hematologic cell lines that are predicted to be not dependent on YAP1/WWTR1 knockdown (P value < 10⁻⁴⁷), whereas mesothelioma histologic subtype was among the most predicted to be YAP/TAZ dependent (Fig. 2A and B). We then validated our ML model by selecting 12 additional cell lines to confirm the predicted YAP/TAZ dependency: 6 cell lines which were predicted to be YAP/TAZ dependent and 6 predicted independent from a variety of different lineages (Fig. 2C–E; Supplementary Fig. S2F–S2J). This provides a landscape of YAP/TAZ dependency across cancer cell line models and enables nomination of cell line models as well as prioritization of cancer indications that would potentially benefit from a Hippo pathway inhibitor.

To functionally annotate gene clusters, we performed a systematic gene signature correlation analysis of the gene clusters with other previously published gene signatures (Supplementary Table S4). As expected, Cluster 2 scores were highly correlated with the previously published YAP gene signature (CORDENONSI_YAP_CONSERVED_SIGNATURE, Pearson ρ: 0.91; Supplementary Fig. S2K). Cluster 1 and 3 scores were associated with several proliferation-associated gene signatures (HALLMARK_MYC_TARGETS_V2, Pearson ρ: 0.95)or epithelial-to-mesenchymal transition (HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, Pearson ρ: 0.87; Supplementary Fig. S2L and S2M), respectively, both of which have been previously implicated in aberrant Hippo pathway signaling (16, 34–36). Interestingly, Cluster 4 was strongly associated with a KRAS dependency gene signature(SINGH_KRAS_DEPENDENCY_SIGNATURE, Pearson ρ: 0.87;Fig. 2F), and although previous reports have suggested YAP1 overexpression as a bypass mechanism to KRAS activation (19), this result suggests that the MAPK pathway may play a role in the context of Hippo signaling. We hypothesized the other gene clusters may also be associated with Hippo signaling although not directly downstream. Beyond Cluster 2 scores as the strongest single predictor of YAP/TAZ dependency, we noted that cell lines with the largest magnitude of decrease in Cluster 4 genes were also those that were most dependent on YAP1/WWTR1 knockdown (Fig. 2G). Taken together, this suggests that additional suppression of MAPK pathway may serve to further enhance therapeutic efficacy of a Hippo pathway inhibitor.

As Hippo pathway inhibitors are under active development, identifying clinically actionable combinations becomes an important next step in augmenting therapeutic response. To determine whether MAPK is a uniquely actionable pathway that cross-talks with Hippo pathway dysregulation, we undertook a chemical genetic screening approach. We screened a drug library of 487 small-molecule compounds in Detroit 562 cells stably transfected with an inducible YAP1 shRNA. Detroit 562 cells are very sensitive to YAP1 knockdown alone, so we decided to knock down only YAP1 (Supplementary Fig. S3A) in our screen to reduce any small hairpin–related RNA toxicity. We assessed whether addition of each compound to the doxycycline-induced knockdown of YAP1 had a greater effect on cell viability than the noninduced shYAP1 arm. We observed that MEK and ERK inhibitors were among the highest-scoring hits showing the largest impact on viability in combination with YAP1 knockdown (adjusted P < 0.1), whereas broad-spectrum cytotoxic chemotherapies did not modulate the effect of YAP1 knockdown, suggesting abrogating MAPK signaling further sensitizes cells to YAP1 knockdown (Fig. 3A; Supplementary Fig. S3B–S3D; and Supplementary Table S5). We sought to expand this observation to a larger panel of YAP1-amplified cell lines treated with several MEK inhibitors including cobimetinib, selumetinib, and PD-901 (Fig. 3B–E; Supplementary Fig. S3E and S3F). Although this sensitization was observed in Hippo pathway–deregulated, YAP1-amplified cell lines, this combination did not show further sensitization in squamous cell cancer lines that lack Hippo pathway alterations (Fig. 3C). The combination of MEK inhibition and YAP1 knockdown promoted caspase-mediated cell death that was measured by increase in caspase 3/7 activity, which was reversed upon treatment with a pan-caspase inhibitor (QVD; Fig. 3F). Clonogenic assays also confirmed the cooperation of YAP1 knockdown and MEK inhibitors in all three YAP1-amplified cell lines (Fig. 3G and H). This combination is unlikely due to general toxicity, as both SK-N-FI (a Hippo-independent model) and MCF10-A (a nonmalignant breast epithelial model) did not show further sensitization (Supplementary Fig. S3G). Importantly, inducible YAP1 depletion in vivo and cobimetinib combination exhibited significant tumor regression in Detroit 562 xenograft model (Fig. 3I; Supplementary Fig. S3H). Taken together, this small-molecule drug library screen provided an orthogonal validation of the role of the MAPK pathway in YAP-dependent cancers.

To further assess the impact of cobimetinib and YAP1 knockdown, in Detroit 562 cells, we performed RNA-seq and ATAC-seq comparing YAP1 knockdown, or cobimetinib treatment, and the combination of YAP1 knockdown and cobimetinib treatment with the control treatment. Together, the combination treatment significantly decreased expression of proliferation genes (Fig. 4A) compared with each individual treatment, consistent with in vitro and in vivo observations (Fig. 3G–I; Supplementary Fig. S3C–S3F). Furthermore, Cluster 2 genes were downregulated in the YAP1 knockdown and the combination treatment but not in cobimetinib treatment alone (Fig. 4A). Conversely, MAPK pathway genes were downregulated upon MEK inhibition and combination treatment but not after YAP1 knockdown alone (Fig. 4A). Although neither Hippo nor MAPK pathway genes were further suppressed upon the combination treatment, we hypothesized that YAP1 knockdown and cobimetinib jointly affect a common node rather than further downregulating their individual pathways. We noted an overlapping 3,479 peaks exhibiting loss of chromatin accessibility in the combination of YAP1 knockdown and cobimetinib treatment were also found in the single treatments alone (Fig. 4B). Motif enrichment analysis revealed decreased chromatin accessibility at TEAD- and AP1-binding sites upon YAP1 depletion and MEK inhibition (Fig. 4C), respectively. Although enrichment was significant upon MEK inhibition, the combination treatment exhibited greater significance of AP1 motif in peaks with loss of chromatin accessibility (Fig. 4C). Previous studies have shown that TEADs and AP1 can coregulate gene transcription through changes in enhancer and promoter regions (33, 37). Taken together, these findings suggest that concomitant YAP1 depletion and MEK inhibition serve to further enhance the loss of AP1-binding sites through Hippo and MAPK pathways, respectively.

As cobimetinib inhibitors affect MEK kinase activity, we performed global phosphoproteomics analysis in Detroit 562 to identify changes in phosphorylation sites across proteins upon YAP1 depletion and/or MEK inhibition. This identified significant changes in 18,800 phosphopeptides across 8,500 proteins. Consistent with the enriched AP1 motif in peaks with loss of chromatin accessibility, we noted 2-fold decrease in FOSL1 phosphorylation (Fig. 4D–F; Supplementary Fig. S4A; Supplementary Tables S6 and S7). In addition, we observed a significant decrease in AP1/FOSL1 target gene expression in the combination treatment (Supplementary Fig. S4B) compared with YAP1 knockdown or cobimetinib treatment alone, suggesting the combination may further contribute to loss of FOSL1 activity. As previous studies have suggested that TEADs (38) may have numerous dimerization partners, we hypothesized that FOSL1 may interact with TEADs (33, 37). Consistent with previous reports, we confirm that TEAD directly interacts with FOSL1 (Supplementary Fig. S4C and S4D) through coimmunoprecipitation. Furthermore, YAP1 is required for FOSL1–TEAD interaction, whereas YAP–TEAD interaction is independent of FOSL1 (Supplementary Fig. S4C). Consistent with these observations, FOSL1–TEAD interaction is abolished upon cobimetinib treatment (Supplementary Fig. S4D), and deletion of FOSL1, together with YAP1, mimics the synergistic effect observed with cobimetinib (Supplementary Fig. S4E and S4F).

Given that the combination of YAP1 depletion and cobimetinib treatment affect TEAD interaction partners, we hypothesized that modulating these two interaction partners may affect TEAD protein stability. Although neither YAP1 depletion nor cobimetinib treatment alone changed TEAD protein levels (Fig. 5A and B), we observed a significant decrease in TEAD protein half-life upon the combination of cobimetinib and YAP1 depletion following cycloheximide (CHX) chase in YAP1-amplified cell lines (Fig. 5A and B; Supplementary Fig. S5A), whereas TEAD transcript levels were unaltered (Supplementary Fig. S5B). Decrease in TEAD protein half-life was reversed upon MG132 treatment 24 hours after CHX treatment (Fig. 5C and D; Supplementary Fig. S5C), suggesting that the decrease in TEAD half-life is mediated by proteasomal degradation. These data imply that the combination of YAP1 depletion and cobimetinib treatment results in a decrease in TEAD stability. Furthermore, we noted that several proliferation genes such as MYC and FOSL2 (Fig. 5E and F) have nearby TEAD- and AP1-binding sites with the potential to regulate their expression. Together, our findings suggest that the convergence of YAP1 depletion and treatment with cobimetinib is mediated through the cooperative interaction between AP1/FOSL1 and TEAD.

In this study, we developed an ML approach to understand Hippo pathway activity (Fig. 1). This has identified a robust, lineage-independent predictive Hippo pathway signature (Fig. 2) and nominated the MAPK pathway as a potential focus for drug combinations that was orthogonally identified and confirmed through a small-molecule drug screen (Fig. 3). Further investigation revealed a novel mechanism in which both Hippo and MAPK pathways regulated TEAD function through decreasing its stability with observed loss of chromatin accessibility at TEAD-binding motifs (Figs. 4 and 5).

The Hippo pathway is emerging as an important area for targeted drug discovery efforts, but greater understanding of this pathway is warranted. As opposed to previous efforts that have derived a curated list of Hippo pathway target genes, here we utilized an ML approach to systematically define four core target gene clusters that are altered as a direct result of loss of Hippo signaling. This lineage-independent approach identified many novel genes not reported in previous YAP/TAZ gene sets. This has enabled accurate prediction of YAP/TAZ dependency in vitro and yielded a signature that can be used to define and prioritize cell line models and patient populations. We validated our gene set to be robust across different genotypes and cancer types. To our knowledge, this is the first lineage-independent, unbiased method that is predictive of Hippo pathway dependency and thus can serve as a valuable tool to identify biomarkers of interest in a tumor-agnostic manner.

We then uncovered the molecular mechanism underpinning the effects of combined inhibition of MAPK and Hippo pathways through several orthogonal technologies. We performed ATAC-seq after YAP1 depletion and/or cobimetinib treatment; results here suggested that the combination converges on the loss of chromatin accessibility at AP1- and TEAD-binding sites. The combination of modulating both TEAD interaction partners, YAP and FOSL1 (via MEK inhibition), resulted in decreased TEAD protein stability (Fig. 5G). Given previous studies have suggested differential regulation of YAP and TAZ (39, 40), future studies will be necessary to elucidate these mechanistic differences, if any, in the context of MAPK pathway inhibition.

Last, this approach serves as a framework for systematic characterization of signaling pathways. Here, we focused on the Hippo pathway in the context of combinations with MAPK inhibitors. Earlier work has shown the cooperative activity of YAP/TAZ/TEAD and AP1 at many enhancers and promoters (33, 37) as well as a role for Hippo pathway inhibitors to combat resistance to BRAF^V600E inhibition (20). As MEK inhibitors are currently in the clinic and Hippo pathway inhibitors are in development, our studies suggest that cotargeting these pathways may achieve a deeper therapeutic response compared with single-agent treatment alone in Hippo pathway–dependent cancers. Our unbiased characterization through a lineage-independent approach to study Hippo pathway activity and dependency underscores the importance of understanding pathway cross-talk as a strategy to nominate potential treatment combinations (Fig. 5G). Our findings here have translational impact not only on Hippo-dependent cancers, but also on tumors with MAPK pathway alterations in primary as well as resistance settings (20, 41).

In summary, our study has made several significant contributions to understanding the Hippo pathway, and in addition, we developed an approach to identify possible pathway targets. This computational–experimental study defines a framework to establish a new paradigm to apply the ML toolbox to accelerate biology and drug development efforts.

Cell lines used in this study (Detroit 562, COLO-680N, HEp-2, and OVCAR-8) were obtained from the ATCC. Detroit 562 with shYAP1 was generated by transfecting in the shYAP1-pLKO lentiviral vector and selecting for Puromycin-positive cells. Cell line authentication was conducted for short tandem repeat (STR) profiling using the Promega PowerPlex 16 System. This is performed once when receiving a new cell line and compared with external STR profiles of cell lines (when available) to determine cell line ancestry. Cell line authentication was routinely conducted by SNP-based genotyping using Fluidigm-multiplexed assays at the Genentech cell line core facility.

Antibodies used in this study are Pan-TEAD [13295, Cell Signaling Technology (CST)]; YAP (14074, CST); p44/42 MAPK (ERK1/2, 4696, CST); Phospho-p44/42 MAPK (p-Erk1/2, 4370, CST); TAZ (70148, CST); MAX (sc-765 and sc-8011, Santa Cruz Biotechnology); α-tubulin (3873, CST); β-actin (3700, CST); cleaved PARP (9541, CST); FRA1(5281, CST); Myc-tag (2278, CST); V5 (680602, Biolegend); anti-rabbit and anti-mouse horseradish peroxidase linked (7074 and 7076, CST), and IRDye anti-rabbit and anti-mouse (68070 and 32211, LI-COR). Cobimetinib was synthesized at Genentech. Selumetinib (S1008), PD0325901 (S1036), and gemcitabine (S1714) were purchased from Selleckchem. CHX solution (C4859) and MG132 solution (M7449) were purchased from Sigma-Aldrich.

Cells were seeded at 500 to 1,000/well on a 96-well plate for 24 hours. They were then treated with indicated siRNAs (final concentration of 25 nmol/L) or indicated inhibitors (with indicated concentration). Cell growth was assessed using Cell Titer-Glo Luminescent Cell Viability Assay (Promega). All cell viability data were collected and calculated for at least six replicates per time point, per condition. IC₅₀ for the inhibitors was determined by fitting the nonlinear regression curve generated by GraphPad Prism.

Cells were seeded at 50,000 to 70,000/well on a 6-well plate for 24 hours. They were then treated with indicated siRNAs (final concentration of 25 nmol/L) or indicated inhibitors (with indicated concentration) for 6 to 10 days. Colony formation was then accessed with crystal violet stain.

In a 96-well plate, cells were treated with indicated siRNAs, cobimetinib, or Q-VD-OPH, a pan-caspase inhibitor (treatment can be alone or in the indicated various combinations). Caspase 3/7 activation was measured using Caspase-Glo 3/7 assay reagent containing proluminescent caspase 3/7 substrate, tetrapeptide sequence DEVD (Promega). The amount of luminescence is displayed as fold changes of treatments to siNTC/DMSO-treated control.

A library comprising 485 compounds including targeted agents, chemotherapeutics, and tool compounds was used to screen for inhibitors that exhibit enhanced efficacy in the context of YAP1 knockdown. Compounds were obtained from in-house synthesis or purchased from commercial vendors and managed as previously described (1). Detroit 562 cells harboring a doxycycline-inducible shRNA targeting YAP1 were treated with and without doxycycline for 72 hours prior to seeding into 384-well plates (BD Falcon) at a density of 1,000 cells per well. Cells were maintained at 37°C with 5% CO₂ and ± doxycycline throughout the course of the experiment. Twenty-four hours after seeding, cells were treated with a nine-point dose titration of each compound or DMSO control (0.1%) for 96 hours. Cell viability was then assessed using CellTiter-Glo (Promega). IC₅₀ (concentration yielding 50% reduction in viability) and mean viability (roughly equivalent to the area under the dose–response curve; ref. 42) values were determined by fitting curves using Genedata Screener software (Genedata). Compounds exhibiting more activity in the context of YAP1 knockdown were determined by calculating both the difference in mean viability and fold change in IC₅₀ between doxycycline-treated and nontreated arms. Screening results can be found in Supplementary Table S3.

For tumor xenograft models, Detroit shYAP1 cells (4.25 × 10⁶) were injected s.c. into right thoracic regions of 6- to 10-week-old Nu/Nu (nude-CRL) mice. When tumors reach a mean volume between 150 and 300 mm³, mice were then grouped into 10 per group for further intervention treatments including doxycycline (1 μg/μL), cobimetinib (7.5 mg/kg), and combination of doxycycline and cobimetinib. Tumor volume was collected and calculated after 3 days and up to 28 days. A small group of mice was euthanized at 6 weeks after treatment for immunoblotting. All animal experiments were approved by Genentech Animal Care and Use Committee.

Cells were lysed in RIPA lysis buffer (89900, Thermo Fisher Scientific) containing protease inhibitor (Roche) and phosphatase inhibitor (Roche). Lysates were prepared by taking supernatants from centrifugation at 12,000 × g for 15 minutes at 4°C. Equivalent amounts of proteins were loaded and separated by SDS-PAGE followed by transferring to membranes.

For endogenous coimmunoprecipitation experiments, 1 × 10⁷ cells were lysed using RIPA buffer (Thermo) and immunoprecipitation with indicated antibody overnight at 4°C. After washing with RIPA buffer (Thermo), coimmunoprecipitated endogenous proteins were then detected by immunoblotting.

Cell pellets were lysed in 8 mol/L urea, 150 mmol/L NaCl, 50 mmol/L HEPES (pH 7.2), complete-mini (EDTA free) protease inhibitor (Roche), PhosSTOP phosphatase inhibitor (Roche), 1 mmol/L Na-orthovanadate, 2.5 mmol/L Na-pyrophosphate, and 1 mmol/L beta-glycerol-phosphate by 15× passages through a 21 g needle. Protein concentrations were then estimated by BCA assay (Thermo Fisher Pierce), and 2.5 mg of protein material was used for further sample preparation. Disulfide bonds were reduced by incubation with 5 mmol/L DTT (45 minutes at 37°C), followed by alkylation of cysteine residues by 15 mmol/L IAA (30 minutes, RT Dark), and finally capped by the addition of 5 mmol/L DTT (15 minutes, RT Dark). Proteins were then precipitated by chloroform/methanol precipitation and resuspended in digestion buffer (8 mol/L urea, 150 mmol/L NaCl, and 50 mmol/L HEPES, pH 7.2). Initial protein digestion was performed by the addition of 1:100 LysC followed by incubation at 37°C for 3 hours. Samples were then diluted to 1.5 mol/L urea with 50 mmol/L HEPES (pH 7.2) before the addition of 1:50 Trypsin and incubation overnight at 37°C. Peptide mixtures were acidified and desalted via solid-phase extraction (SPE; SepPak–Waters).

Before phosphopeptide enrichment, a 100 μg aliquot of peptides was saved for global proteome analysis. Phosphopeptides were enriched from the remaining material utilizing iron-IMAC–based magnetic beads as previously described (43). To enable multiplexed quantitation, both the phosphopeptide-enriched and unenriched peptide pools were resuspended in 200 mmol/L HEPES (pH 8.5) and mixed with tandem mass tags (TMT; Thermo Fisher Pierce) at a label to protein ratio of 2:1. After 1 hour of labeling, the reaction was quenched by the addition of 5% hydroxylamine and incubated at room temperature for 15 minutes. Labeled peptides were then mixed, acidified, and purified by SPE.

Labeled samples were separated by offline high pH reverse-phase fractionation using an ammonium formate–based buffer system delivered by a 1100 HPLC system (Agilent). Peptides were separated over a 2.1 × 150 mm, 3.5 μm 300Extend-C18 Zorbax column (Agilent) and separated over a 75-minute gradient from 5% ACN to 85% ACN into 96 fractions. The fractions were concatenated into either 24 or 12 samples for proteome and phosphoproteome samples, respectively. Fractions were concatenated by mixing different parts of the gradient to produce samples that would be orthogonal to downstream low pH reverse-phase LC-MS/MS. Samples were dried, desalted by SPE, and dried again.

NanoLC-MS/MS analysis was performed on an Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific) coupled to a Dionex Ultimate 3000 RSLC-nano (Thermo Fisher Scientific). Peptides were separated over a 100 μm × 250 mm PicoFrit column (New Objective) packed with 1.7 μmol/L BEH-130 C18 (Waters) at a flow rate of 450 nL/min for a total run time of 180 minutes. The gradient spanned from 2% Buffer B (0.1% FA/98% ACN/2% water) to 30% B over 155 minutes and then to 50% B at 160 minutes. For mass spectrometry analysis, peptides were surveyed within FTMS¹ analyses [120,000 resolution, AGC = 1 × 10⁶, and maximum injection time (max IT) = 50 ms], and the top 10 peaks were selected for MS/MS, ensuring that any given peak was selected only once in a 45-second window (ppm tolerance = 5 ppm). For dynamic exclusion, the “one precursor per charge state” was ON for proteome analysis and OFF for phosphoproteome analysis. For MS² analysis, precursors were isolated using the quadrupole (0.5 Th window), fragmented using CAD (NCE = 35, AGC = 2 × 10⁴, and max IT = 100 ms), and analyzed in the ion trap (scan speed = Turbo). Following MS² analysis, the top 8 (proteome) or 6 (phosphoproteome) ions were simultaneously selected [synchronous precursor selection – SPS, AGC = 2.5 × 10⁵ (proteome) or 3.0 × 10⁵ (phosphoproteome), and max IT = 150 (proteome) or 200 (phosphoproteome) ms] and fragmented by HCD (NCE = 55) before analysis in the Orbitrap (resolution = 50,000). Raw data files are available via the MASSIVE data repository using the identifier MSV000085111.

All mass spectrometry data were searched using Mascot against a concatenated target-decoy human database (downloaded June 2016) containing common contaminant sequences. For the database search, a precursor mass tolerance of 50 ppm, fragment ion tolerance of 0.8 Da, and up to 2 missed cleavages were used. For global proteome and phosphoproteome analysis, carbamidomethyl cysteine (+57.0214), and TMT-labeled n-terminus and lysine (+229.1629) were applied as static modifications, whereas methionine oxidation (+15.9949) was set as a dynamic modification. For phosphoproteome analysis, TMT-labeled tyrosine (+229.1629) and phosphorylation of serine, threonine, and tyrosine (+79.9663) were also set as dynamic modifications. Peptide spectral matches for each run were filtered using line discriminant analysis to an FDR of 2% and subsequently as an aggregate to a protein level FDR of 2%. Localization of phosphorylation sites was performed using a modified version of the Ascore algorithm. TMT-MS3 quantification was performed using Mojave, with only those peptide-spectrum matches (PSM) possessing isolation specificities greater than or equal to 0.5 considered for the final dataset. Intensities of each PSM were added to the peptide and then protein (proteome) or phosphoisoform (phosphoproteome) level. Expression is reported as relative abundance, which is the measured intensity of any given channel divided by the total intensity for that protein or phosphoisoform.

Tumors and cell lines were dissociated and lysed for RNA isolation using RNAeasy Mini Kit (QIAGEN) with the on-column DNA elimination. cDNA was prepared by reverse transcription using the iScript cDNA Synthesis Kit (Bio-Rad) as the manufacturer's protocol. The quantitative RT-PCR was performed using QuantStudio 7 Flex machine with TaqMan probes for YAP1 (Hs00902712_g1), WWTR1 (Hs00210007_m1), TEAD1 (Hs00173359_m1), TEAD2 (Hs01055894_m1), TEAD3 (Hs00243231_m1), TEAD4 (Hs01125032_m1), GAPDH (Hs02786624_g1), CTCF (Hs00902016_m1), CYR61 (Hs00155479_m1), CPA4 (Hs00275311_m1), THBS1 (Hs00962908_m1), PTRF (Hs00396859_m1), EPHA2 (Hs01072272_m1), EDN1(Hs00174961_m1), TSPAN3 (Hs00170681_m1), NDRG1 (Hs00608387_m1), and EMP2 (Hs00171315_m1; Applied Biosystems). Relative expression of each gene to GAPDH of target genes was assessed for at least two to three biological replicates.

Total RNA was extracted using the QIAGEN RNAeasy Mini Kit (QIAGEN) with the on-column DNA elimination. The concentration was identified using NanoDrop 8000 (Thermo Fisher Scientific). Quality control was done by determining RNA integrity with Bioanalyzer 2100 (Agilent Technologies). About 500 ng of RNA was used for library synthesis using the TrueSeq RNA Sample Preparation Kit v2 (Illumina). Size of the libraries was confirmed using 2200 TapeStation and High Sensitivity D1K screen tape (Agilent Technologies), and the concentration was determined by a qPCR-based method using the Library Quantification Kit (KAPA). The libraries were multiplexed and then sequenced on Illumina HiSeq2500 (Illumina) to generate 30M of single-end 50-base pair reads.

The RNA-seq data have been deposited in the Gene Expression Omnibus (GEO) with the accession code GSE161019. Data or other materials are available from the corresponding authors upon request.

For RNA-seq data analysis, RNA-seq reads were first aligned to ribosomal RNA sequences to remove ribosomal reads. The remaining reads were aligned to the mouse reference genome (GRCm38) using GSNAP (44, 45) version “2013-10-10,” allowing a maximum of two mismatches per 75 base sequence (parameters: -M 2 -n 10 -B 2 -i 1 -N 1 -w 200000 -E 1 –pairmax-rna = 200000 –clip-overlap; refs. 44, 46). Transcript annotation was based on the Ensembl genes database (release 77). To quantify gene expression levels, the number of reads mapped to the exons of each RefSeq gene was calculated.

Differential gene expression was performed with DESeq2 (47). A prefilter was applied, so that only genes with at least a median per kilobase per million mapped reads value of 10 in one condition were analyzed. P values for other genes were simply set to 1, and log fold changes were set to 0 for visualization purposes, but such genes were not included in the multiple testing correction. Q values were obtained by correcting P values for multiple hypotheses using the Benjamini–Hochberg procedure. Genes were considered if they had a Q value of less than 0.05 and were protein coding. Counts were transformed to log₂ counts per million, quantile normalized, and precision weighted with the “voom” function of the limma package (48).

Reads were aligned to the human reference genome (NCBI Build 38) using GSNAP31 version “2013-10-10,” allowing a maximum of two mismatches per read sequence (parameters: -M 2 -n 10 -B 2 -i 1–pairmax-dna = 1000–terminal-threshold = 1000–gmap-mode = none–clip-overlap). Reads aligning to locations in the human genome that contain substantial sequence homology to the MT chromosome or to blacklisted regions identified by the ENCODE consortium were omitted from downstream analyses. Properly paired reads derived from nonduplicate sequencing fragments were used to quantify chromatin accessibility according to the ENCODE pipeline standards with minor modifications as follows. Accessible genomic locations were identified by calling peaks with Macs2 (49) using insertion-centered pseudo-fragments (73 base pairs; community standard) generated on the basis of the start positions of the mapped reads. Accessible peak locations were identified as described: in brief, we called peaks on a group-level pooled sample containing all pseudo-fragments observed in all samples within each group. Peaks in the pooled sample that were independently identified in two or more of the constituent biological replicates were retained for downstream analysis, using the union of all group-level reproducible peaks (https://www.encodeproject.org/atac-seq/#standards). We quantified the level of chromatin accessibility within each peak for each replicate as the number of pseudo-fragments that overlapped the peak in question and normalized these estimates using the Trimmed Mean of M-values (TMM) method (50). We identified differentially accessible peaks between groups in the framework of a linear model implemented with the limma R package (51) and incorporating precision weights calculated with the voom function in the limma R package (49). We identified enriched TF motifs using HOMER v4.7 (52). To evaluate the significance of the TF enrichment, we defined peaks as significantly differentially accessible based on a range of FDR-adjusted P value thresholds between 1 and 0.01 and a |log₂[fold change]| in accessibility ≥ 1 (estimated from the model coefficients). Given the strong enrichment of the top motifs across a wide range of P-value cutoffs, we decided to consider peaks as different across groups for a |log₂[fold change]| ≥ 1 and FDR-adjusted P value ≤ 0.05 in subsequent analyses.

The ATAC-seq data have been deposited in the GEO with the accession code GSE161019. Data or other materials are available from the corresponding authors upon request.

Statistically significant downregulated genes in three or more of the seven cell lines were considered for common downregulated genes. Genes were subsequently filtered based on their expression, and genes that were upregulated in any cell line upon YAP1/WWTR1 knockdown were removed. Signed coexpression networks were built using the WGCNA package in R (minModuleSize = 10, reassignThreshold = 1e–6, deepSplit = 2, mergeCutHeight = 0.15) using RNA expression of pan-cancer cell lines (53).

For each differential accessibility comparison (e.g., siYAP1/WWTR1 vs. siNTC), ATAC-seq peaks were found within a window of 1,000 bp before the start and 1,000 bp after the end of every gene in human GENCODE geneset, release 27 (genome assembly version GRCh38). The peaks were then defined as opening (log fold change > 1, adjusted P value < 0.05), closing (log fold change < -1, adjusted P value < 0.05), or background (everything else), and only peaks that were annotated as “protein-coding” were retained. The total number of closing peaks within the 1,000 bp window of n genes in each Hippo pathway cluster was compared with the total number of closing peaks within 1,000 bp of n randomly sampled genes, repeated for 1 million permutations. The P value was taken as the number of times the number of randomly sampled closing peaks was greater than that of the Hippo cluster's, divided by 1 million.

YAP1/WWTR1 knockdown was performed on 42 cell lines, and change in viability was assessed using Cell Titer Glo after 7 days. Effect of viability was quantified as a ratio of fraction of remaining luminescence compared with nontargeting control (i.e., 1 would mean complete loss of luminescence and 0 as no change). To predict this effect of viability, we constructed a random forest model using the expression of Cluster 1 to 4 genes of 31 (70%) cell lines to predict the effect of viability. Predication accuracy was determined by minimizing root mean square error from 5-fold cross-validation of the training data.

We identified two additional YAP/TAZ gene set: Cordenonsi and colleagues (21) or Wang and colleagues (TCGA; ref. 31). In addition, we included our Cluster 0 gene set (i.e., genes that were significantly downregulated by YAP1/WWTR1 knockdown but not assigned into a gene cluster, likely not proximal to Hippo pathway activity) as a negative control. First, we trained an ML model using either our gene set (Cluster 2), Cordenonsi and colleagues' gene set (21), or Wang and colleagues' gene sets (TCGA; ref. 31) using the same training data and ML algorithm (parallel random forest, caret v6.0-85). The ML model trained using Cluster 2 gene set outperformed the ML models trained from Cordenonsi and colleagues (21) and Wang and colleagues (31) in both root mean square error and R-squared metrics (Supplementary Fig. S2B). This was not dependent on the ML algorithm used to train the ML model—we repeated this using two other common ML algorithms: support vector machines and boosted generalized linear model. For both ML algorithms, our Cluster 2 gene set outperformed the other gene sets. Last, this difference in performance is not due to the initial training data used in model training. We created 500 random iterations of training data (createDataPartition, caret v6.0-85), and our gene set significantly outperformed the other two gene sets (Wang and colleagues, P value < 10⁻⁷¹; Cordenonsi and colleagues, P value < 10⁻⁶⁵; paired t-test). Importance of individual genes in a given ML model was calculated using varImp (caret v6.0-85). Genes within our gene were defined as novel if not previously reported in either Cordenonsi and colleagues (21) or Wang and colleagues (31) and known if present in one or both gene sets.

Besides the RNA-seq and ATAC-seq studies that have their own statistical analysis, all of the statistical analysis for in vivo and in vitro studies was done using the Student t test (two-tailed, unpaired) to compare treatment groups with control group. All of the in vitro experiments were performed at least 3 times. For the in vivo work, the gender, age, and weight of animals were matched, and the sample size is 3 to 10 mice per group. A P value of less than 0.05 was considered as statistically significant (*, P < 0.05). Significant values as well as number of replicates are noted for each experiment in the respective figure legends.

siRNA transfection was performed in Detriot-562 and PATU-8902 cell lines (obtained from the ATCC). Cells were seeded at 15,000/well on a 6-well plate or 1,000/well on a 96-well plate for 24 hours. They were then treated with the siRNAs (final concentration of 20 nmol/L) using Lipofectamine RNAi Max in serum-free RPMI media for 72 hours before collection. siRNAs were purchased from Dharmacon, including siNTC (D-001810-10), siYAP1 (L-012200-00), siWWTR1 (L-016083-00), siTEAD1 (J-012603-05), siTEAD2 (J-012611-09), siTEAD3 (L-012604-00), siTEAD4 (J-019570-09), and siFOSL1 (LQ-004341-00). Knockdown was confirmed by Western blot with the following antibodies from CST: Pan-TEAD (13295), YAP (14074), TAZ (70148), FOSL1 (5281), and β-actin (3700).

ATAC-seq was performed as previously described (54, 55). A total of 1 × 10⁵ cell pellets were washed with PBS, and cells were pelleted by centrifugation. Cell pellets were resuspended in 100 μL of ATAC-resuspension buffer (10 mmol/L Tris HCl, pH 7.4, 10 mmol/L NaCl, and 3 mmol/L MgCl₂), and nuclei were pelleted by centrifugation. The nuclei were resuspended in 50 μL reaction buffer containing Tn5 transposase (2.5 μL Tn5 transposase, 25 μL 2 × TD buffer, 16.5 μL PBS, 0.5 μL 1% digitonin, 0.5 μL 10% Tween-20, and 5 μL H2O; Illumina). The reaction was carried out at 37°C for 30 minutes. Tagmented DNA was isolated by MinElute PCR Purification Kit (QIAGEN). Libraries were later generated and sequenced on NextSeq500 (Illumina).

CRISPR transfection was performed in Detroit 562 and OVCAR-8 cell lines (obtained from the ATCC). Cells were seeded at 10,000/well on a 6-well plate for 24 hours. They were then transfected with the RNP complex that contains designed guide RNA siRNAs (final concentration of 1 μmol/L) using Lipofectamine CRISPR Max Cas9 transfection reagent for 72 hours before collection. Cas9 protein was purchased from Invitrogen (Truecut Cas9 protein v2, A364498). crRNA and transcrRNA were purchased from IDT. Multiple guide RNAs for YAP1 were tested to get to the two optimized sequences of gYAP1: ACATCGATCAGACAACAACA and gYAP2: CCACAGGGAGGCGTCATGGG. Transient knockout/knockdown was confirmed by Western blot with the following antibodies from CST: YAP (14074) and β-actin (3700).

C.M. Rose, R. Piskol, V.C. Pham, M.T. Chang, and A. Dey report they are employees of Genentech and shareholders at Roche. S.E. Martin reports personal fees from Genentech outside the submitted work. Z. Modrusan reports personal fees from Roche stock outside the submitted work. C. Klijn reports he was an employee of Roche and shareholder at the time when the work was performed (currently no longer applicable). S. Malek reports other from Genentech/Roche during the conduct of the study and outside the submitted work, and is an employee and shareholder of Genentech/Roche. No disclosures were reported by the other authors.

T.H. Pham: Formal analysis, validation, investigation, writing-original draft, writing-review and editing. T.J. Hagenbeek: Validation, investigation, writing-original draft. H.-J. Lee: Conceptualization, validation, investigation. J. Li: Software, formal analysis. C.M. Rose: Software, formal analysis, investigation. E. Lin: Formal analysis, validation, investigation. M. Yu: Formal analysis, validation, investigation. S.E. Martin: Formal analysis, supervision, writing-review and editing. R. Piskol: Formal analysis, supervision, methodology. J.A. Lacap: Validation, investigation. D. Sampath: Supervision, investigation. V.C. Pham: Validation, investigation. Z. Modrusan: Data curation, formal analysis, writing-review and editing. J.R. Lill: Supervision, visualization, writing-review and editing. C. Klijn: Formal analysis, supervision, methodology. S. Malek: Supervision, writing-review and editing. M.T. Chang: Conceptualization, software, formal analysis, methodology, writing-original draft, writing-review and editing. A. Dey: Conceptualization, supervision, writing-original draft, writing-review and editing.