Introduction

The liver employs two main types of epithelial cells---hepatocytes, the predominant parenchymal cell type, and cholangiocytes, which line the bile ducts---to provide key metabolic functions, including the synthesis and secretion of bile acids and detoxification (Stanger, 2015). Early embryonic development (E8-E9.5) in mammals integrates signals from multiple signaling pathways, notably Wnt, FGF and BMP, to initiate and drive hepatic development, including the generation of embryonic bi-potent progenitor cells called hepatoblasts (Si-Tayeb et al., 2010). Active Notch signals direct hepatoblasts to differentiate to the cholangiocyte fate whereas an absence of such signals enables hepatocyte differentiation (Adams & Jafar-Nejad, 2019). Downstream of Notch and these other signaling pathways, transcription factor networks control cell fate, proliferation and differentiation; in particular, HNF4A functions centrally in hepatoblast specification, with loss or gain of CEBPA determining cholangiocyte versus hepatocyte fate, respectively (Tachmatzidi et al., 2021).

Liver cancers rank among the world’s most commonly diagnosed cancers and are the third leading cause of cancer mortalities (Sung et al., 2021). Hepatocellular carcinomas (HCCs), which constitute the majority of liver cancers, primarily arise in the context of cirrhosis, which is commonly caused by Hepatitis B or C virus (HBV, HCV) infections, toxin exposure, alcoholism, or non-alcoholic steatohepatitis (NASH) (Llovet et al., 2022). HCC is commonly diagnosed at advanced stages. Curative resection or ablation is possible only in early stage tumors and recurrence rates are high (Llovet et al., 2022). Systemic therapies now include tyrosine kinase inhibitors (da Fonseca et al., 2020) as well as the combination of the therapeutic antibodies atezolizumab and bevacizumab to block immune checkpoints and angiogenesis, respectively (Finn et al., 2020). While HCC management has improved patient outcomes over the past decade, there remains a strong medical need to develop improved adjuvant therapies and increase the efficacy and duration of systemic therapies (Llovet et al., 2022).

The term “HCC” serves as an umbrella classifier for a heterogenous cancer type composed of multiple subtypes that reflect differences in etiology, cell of origin, and driver signaling pathways and mutations. Histological classifications broadly categorize tumors as poorly- or well-differentiated and further sub-group tumors based on content of cholangiocyte or hepatocyte features, including mixed-lineage and progenitor tumors (Calderaro et al., 2019). These features have also been correlated with different types of molecular classifiers. The most frequent genetic mutations highlight dysregulated cell cycle control, telomerase activation, Wnt signaling, epigenetic regulation and oncogenic tyrosine kinase (as evidenced in part by FGF19, VEGFA, and MET chromosomal amplifications) and MAPK signaling (Llovet et al., 2022). Transcriptome meta-analyses of global patient populations have helped classify HCCs into proliferative and non-proliferative groups, which largely correspond to the poorly- and well-differentiated classes, respectively. Within the proliferative population, several groups have molecularly identified a progenitor subclass (Ally et al., 2017; Boyault et al., 2007; Hoshida et al., 2009), characterized by the expression of stem/progenitor cell markers (Hoshida et al., 2009) and an aggressive, poor-prognosis phenotype (Llovet et al., 2022).

Notch signals, conserved throughout metazoans, vitally contribute to specifying the fate and function of neighboring cells, including progenitors during liver development (Siebel & Lendahl, 2017). In mammals, binding of cell surface ligands of the Jagged (Jag) or Delta-like (Dll) families conformationally alters the negative-regulatory region (NRR) of one of four Notch receptors (NOTCH1-4) (Gordon et al., 2015) to ultimately trigger activation through a gamma-secretase-catalyzed proteolytic cleavage (De Strooper et al., 1999). The released NOTCH intracellular domain (NICD) translocates into the nucleus, assembles into a transcription factor complex that includes the DNA-binding protein RBPJ, and activates the downstream transcriptional program. Given that Notch signaling has been linked to numerous human diseases, including cancer, inflammation and metabolic disorders (Allen & Maillard, 2021; Aster et al., 2017; Xu & Wang, 2021), researchers have developed numerous Notch-blocking therapeutics, including gamma-secretase inhibitors and monoclonal antibodies that selectively target individual Notch ligands and receptors (Allen & Maillard, 2021).

Mutations in Notch receptors that activate or prolong Notch signaling point to an oncogenic role in multiple cancers, including T cell acute lymphoblastic leukemia (T-ALL), triple-negative breast cancer and adenoid cystic carcinoma (Aster et al., 2017). However, such mutations have not been broadly described in liver cancers, and conflicting reports obscure a clear understanding of Notch signaling in HCC. Whereas some reports suggest that Notch signals suppress HCC development (Viatour et al., 2011), other studies conclude that hepatocyte Notch signaling can drive fibrosis in NASH (Zhu et al., 2018), oncogenesis in HCC (Dill et al., 2013) and expression of a conserved gene signature in HCC patient populations (Villanueva et al., 2012). Furthermore, achieving the therapeutically vital goal of defining which Notch ligands and receptors function in HCC requires moving past gain-of-function methods to approaches that precisely inhibit individual ligands and receptors in physiologically relevant preclinical models.

By conducting a preclinical trial of a JAG1-selective blocking antibody in a large panel of patient-derived xenograft models, we now report the discovery of HCCs that exquisitely depend on a JAG1-NOTCH2 signaling axis within tumor cells. Molecular characterization shows that these tumors display hallmarks of embryonic liver progenitors in the aggressive, proliferative HCC subclass. Strikingly, selective inhibition of JAG1 or NOTCH2 leads to profound tumor regressions. Single-cell sequencing analyses of these tumors and their responses following treatment reveal a mechanism in which Notch inhibition induces cycle arrest and promotes hepatocyte differentiation, through upregulation of CEBPA expression and activation of existing HNF4A, mimicking normal developmental programs. Our results thus identify a Notch-dependent subset of HCCs within the progenitor subclass, outline the relevant transcriptional and differentiation mechanisms, and highlight JAG1 and NOTCH2 as actionable targets in this poor-prognosis group.

Results

Discovery of a Notch-dependent HCC subclass

Previous work from our group implicated JAG1 and NOTCH2 as a functionally relevant Notch ligand-receptor pair that supported growth of HCC and cholangiocarcinoma (CC) tumor types (Huntzicker et al., 2015a). However, this study largely relied on the in vivo transfection of mouse liver cells with constitutively active N-Ras and Akt. To test the hypothesis that JAG1 and NOTCH2 were drivers of HCC in a pre-clinical system using patient samples, we moved to patient-derived xenografts (PDXs) involving the direct transfer of tumor samples from patient to an immunocompromised mouse (Siolas & Hannon, 2013). Relying on transcriptome data that showed expression of JAG1 and NOTCH2 as well as evidence of pathway activity, such as expression of HES1, we chose two models for initial tests of sensitivity to aJ1.b70, a potent and selective therapeutic antibody that inhibits JAG1 (Lafkas et al., 2015). While JAG1 blockade had no effect on the growth of one PDX model, aJ1.b70 treatment potently inhibited growth of the second model, LIV78, in a dose-dependent manner (Fig. 1A). This striking response---regression of an established solid tumor, induced by a single-agent targeting Notch signaling---was unprecedented in our experience of manipulating Notch signaling in solid cancers.

Figure 1.

Given that we discovered this Notch-dependent tumor in one of two PDX models using simple screening criteria and that the majority of PDX models and HCC tumors express JAG1 and NOTCH2 (Huntzicker et al., 2015a), TCGA Research Network: https://www.cancer.gov/tcga), we initiated a pre-clinical trial to assess the frequency of response to Notch inhibition. We tested aJ1.b70 in a panel of 42 additional PDX liver cancer models (comprising 35 HCC models plus seven CCs and other liver tumor subtypes), chosen without bias or prescreening, including for Notch pathway status. We used five mice per control and anti-JAG1 groups and dosed for approximately three weeks, given the ambitious scope of the study. Examining the percentage of aJ1.b70 tumor growth inhibition (%TGI) revealed that Jag1 blockade had little to no effect on growth of 30 tumors, although %TGI varied significantly between individual tumors, likely reflecting the small group sizes and short growth period (Fig. 1B, white panel). Three models, including a CC and an adenosquamous carcinoma, trended towards enhanced growth following treatment (Fig. 1B, red area). We classified 9 tumors as partial responders, with %TGI ranging from approximately 26 to 57 (Fig. 1B, light green area). Thus, this unbiased in vivo screen revealed only a modest JAG1-dependence (21%) within this cohort of PDX. LIV78 (dark green area), however, strikingly stood out as showing treatment-induced regressions.

We expanded our screen, searching for additional “super-responder” models that shared features of LIV78. In particular, LIV78 not only expressed high levels of JAG1 and NOTCH2, but also expressed multiple markers of embryonic liver progenitor cells. Using this approach, we identified four additional models and tested them for sensitivity to aJ1.b70 in vivo. Three models responded moderately or not at all (Fig. 1C). However, one model, Li1035, responded dramatically (Fig. 1C), showing tumor regressions in a dose-dependent manner (Fig. 1D) and revealing a JAG1 dependence similarly strong as that observed in LIV78.

To identify the Notch signaling components expressed in both the human cancer cells and mouse stromal cells that typically comprise PDX tumors, we performed bulk RNA sequencing. Tumors from LIV78 and Li1035 were composed of approximately 95% human cancer cells and 5% mouse stromal cells, based on the percentages of human versus mouse sequencing reads (Supplemental Fig. 1A). Among the four canonical ligands, JAG1 was the only one expressed by the human cancer cells (Fig. 2A), consistent with growth inhibition following JAG1 blockade. In contrast, Dll4 was the only ligand significantly expressed in mouse cells, likely reflecting the well-characterized DLL4 expression in endothelial cells (Mack & Luisa Iruela-Arispe, 2018) (Fig. 2A). Both LIV78 and Li1035 cancer cells expressed human NOTCH1 and NOTCH2 whereas the mouse stromal cells expressed Notch1, Notch3 and Notch4 (Fig. 2A). To determine whether NOTCH1, NOTCH2 or both were the relevant receptors mediating Notch dependency, we treated LIV78 or Li1035 PDX models with blocking antibodies that selectively and potently inhibit each receptor (Wu et al., 2010). Similar to treatment with aJ1.b70, treatment with a NOTCH2 blocking antibody (aNRR2) induced regressions of LIV78 and Li1035 tumors (Figs. 2B and C) whereas a NOTCH1 blocking antibody (aNRR1) showed only a modest and delayed effect (Fig. 2C). Immunohistochemical staining with an antibody targeting an epitope in the NOTCH2 intracellular domain (NICD2, cleaved and uncleaved forms) revealed clear NOTCH2 staining in the majority of tumor cells in both models, including membranous/diffuse (inactive) and nuclear (active) staining. Active (nuclear) NOTCH2 signaling was observed in approximately 30-50% of tumor cells. (Fig. 2D.i and Supplemental Fig. 1B.i), while nuclear NICD2 staining was not observed in non-regressing models (Fig. 1C, Supplemental Fig. 1C). Most importantly, tumors harvested after treatment with aJ1.b70 or aNRR2 no longer showed nuclear NICD2 staining, consistent with effective blockade of the Notch2 signal (Fig. 2D.ii-iii and Supplemental Fig. 1B.ii-iii). In contrast, Notch1 signaling was detected only in what morphologically appeared to be endothelial cells (Fig. 2D.v-viii and Supplemental Fig. 1B.v-viii); and treatment with aNRR1 inhibited the endothelial Notch1 signal (Fig. 2D.viii and Supplemental Fig. 1B.viii) without affecting the NICD2 nuclear signal in tumor cells (Fig. 2D.iv and Supplemental Fig. 1B.iv).

Figure 2.

To determine the functionally relevant source of JAG1, which could be human JAG1 expressed on the cancer cells or mouse JAG1 expressed on stromal cells, we developed a novel anti-JAG1 blocking antibody. This antibody, aJ1.b6A9, is a mouse anti-human IgG that selectively inhibits human JAG1 but does not inhibit mouse JAG1 (Supplemental Fig. 2D). Treatment with this human-specific antibody induced regressions in both tumor models (Figs. 2E-F, purple line) in the same manner as did aJ1.b70, which cross-reacts with human and mouse JAG1 (Lafkas et al., 2015). Taken together, our results demonstrate that a JAG1-NOTCH2 signaling axis within the human cancer cell population drives the Notch dependency in the super-responder models, with NOTCH2 serving as the functionally relevant receptor.

Co-expression of cholangiocyte and hepatocyte features distinguish the Notch-dependent HCC tumor models

To better understand the hallmarks of the Notch-dependent tumors, we performed RNA sequencing on all HCC models within the PDX panel and focused on the tumor (human) transcriptomes. K-means analysis of the most variably expressed genes yielded 5 clusters (Fig. 3A), with the sensitive models LIV78 and Li1035 forming a sub-cluster within Cluster 3 (Fig. 3A, B). Comparing the expression signatures of these clusters to cell-type specific signatures generated for normal human fetal and adult liver (Segal et al., 2019) revealed that Cluster 3 is enriched for features of fetal hepatocytes and scores positive for signatures of adult hepatocytes and mixed-lineage hepatobiliary hybrid progenitors (HHyPs) (Segal et al., 2019) (Supp. Fig. 2A).

Figure 3.

We then focused the gene signature comparison to the S1, S2 and S3 HCC subtypes (Hoshida et al., 2009). S1 and S2 subtypes represent aggressive and poorly/moderately-differentiated tumors, enriched for TP53 mutations and WNT signaling (but not CTNNB1 mutations); S2 tumors typically express AFP, EPCAM, and GPX3, and are marked by IGF2 signaling, a hepatoblastoma-like gene signature, and cancer stem cell features, while S3 tumors are the least aggressive and hepatocyte-like (Goossens et al., 2015a). We found that Cluster 3 was the most S2-like and least S1-like when compared to the other four clusters (Fig. 3C). Although the S2 classification enriched for Notch-dependent models, it was insufficient to fully predict this dependency, as non-responding models also fell within Cluster 3 and the Hoshida S2 type (Fig. 3C). Histological examination of models in Cluster 3 did not reveal any obvious distinguishing features between sensitive and insensitive models also (Supplemental Fig. 2B). However, we noticed that LIV78 and Li1035 express AFP at a low level, atypical for S2 tumors (Fig. 3D). Furthermore, differential gene expression analysis revealed that the Notch-dependent PDX models show a distinguishable upregulated expression of FGF9 and GAS7 (Figs. 3E, F, Supplemental.Table 1), compared to the insensitive models within Cluster 3. In contrast to the PDX dataset, however, FGF9 and GAS7 expression was readily detectable in the majority of patient samples and did not enable subset classifications (Supplemental Fig. 2C). We speculate that FGF9/GAS7 expression in non-tumor cells, absent in our PDX analysis but present in patient samples, may limit the usefulness of these markers as clinical diagnostics. Nevertheless, our results highlight the discovery of an HCC subtype within the aggressive and poorly differentiated S2 class, that exquisitely depends on Notch signaling and resembles a hybrid progenitor cell expressing FGF9 and GAS7.

Notch inhibition promotes tumor cell differentiation towards the hepatocyte lineage

To elucidate the molecular consequences of Notch inhibition, we performed RNA sequencing of LIV78, Li1035 and three non-responsive models within Cluster 3 (Li0050, Li0240, and Li0271) after treatment with various Notch-blocking and control antibodies. At 72 hours after treatment, JAG1 blockade induced expression changes in hundreds of genes in both the LIV78 and Li1035 tumors; in clear contrast, few or no changes were observed in the three insensitive models (Fig. 4A, Supplemental Fig. 3A, Supplemental Table2). These results likely reflect the presence of an active JAG1-dependent transcriptional program in the regressing models and its absence in the insensitive models. Gene expression changes started as early as 8 hours after treatment, increased in number at 24 and 72 hours, and persisted at least until seven days (Supplemental Fig. 3B, C, Supplemental Table2). NOTCH1 inhibition failed to induce changes in gene expression (Fig. 4A, Supplemental Fig. 3B, Supplemental Table2), consistent with the lack of tumor growth inhibition (Fig. 2C). In contrast, gene expression changes induced by NOTCH2 inhibition were strikingly similar to those caused by JAG1 blockade in both sensitive models (Fig. 4A, Supplemental Fig. 3B, Supplemental Table2), supporting the conclusion that JAG1 and NOTCH2 co-drive the Notch transcriptional program.

Figure 4.

To ascertain the consequences of these changes, we mapped the treatment-induced expression signatures onto published signatures of normal liver development and the HCC signatures defined by Hoshida et al. Blockade of JAG1 or NOTCH2 (but not NOTCH1) drove the signatures away from those that describe hybrid/progenitor states towards those of differentiated hepatocytes, with the magnitude and significance of the shifts increasing over the first three days (Fig. 4B). The gene-specific changes underlying this progenitor-to-hepatocyte shift were nearly identical in LIV78 and Li1035 (Fig. 4C) and included decreases in common Notch transcriptional targets (Fig. 4D-E, Supplemental Fig. 3D-F) and increases in hepatocyte metabolic enzymes (Fig. 4C, E). Together, these results indicate that a JAG1-NOTCH2 signaling axis, not involving NOTCH1, drives a shared Notch transcriptional program in the Notch-dependent super responder HCC models.

Gene set enrichment analysis (Liberzon et al., 2015; Subramanian et al., 2005) of the down-regulated transcriptional program not only reinforced the connection to Notch signaling but also highlighted other developmental signaling pathways (‘Wnt beta-catenin signaling’ and ‘TNFa signaling via NFkB’) (Fig. 4E, Supplemental Fig. 3G, Supplemental Table3). JAG1 inhibition also reduced the expression of SOX4, SOX9, SPP1, and HNF1B (Supplemental Fig. 3I-L)---genes that mark cholangiocytes and bi-potent liver progenitors (Coffinier et al., 2002; Lesaffer et al., 2019; Poncy et al., 2015). Similarly down-regulated were gene sets associated with cell proliferation (‘E2F targets’, G2M checkpoint’, ’Myc targets 1 and 2’) (Supplemental Fig. 3G,H, Supplemental Table3). DNA content analysis confirmed that JAG1 inhibition reduced the number of cells in the S and G₂-M phases (Supplemental Fig. 3M), indicating a treatment-induced cell cycle block in G₀-G₁, further supported by reductions in Ki67 staining and BrdU incorporation (Supplemental Fig. 3N,O).

JAG1 inhibition in LIV78 and Li1035 led to increased expression of ‘Hallmark’ gene sets related to ‘xenobiotic metabolism’ or lipid and bile-acid biosynthesis processes (‘Fatty Acid Metabolism’, ‘Peroxisome’, ‘Cholesterol Homeostasis’, ‘Bile acid Metabolism’) (Fig. 4E, Supplemental Fig 3G,H, Supplemental Table3), all specialized functions of mature liver hepatocytes (Stanger, 2015; Trefts et al., 2017). The top upregulated transcripts included genes encoding glutathione transferases (GSTA1-3), essential for enzymatic detoxification of electrophilic compounds such as xenobiotics and oxidative stress products, as well as several members of the cytochrome P450 family, which catalyze endobiotic oxidation and xenobiotic detoxification (Fig. 4C) (Bachmann, 1996; Ramsay & Dilda, 2014). Hepatocytes are also key to regulating lipid homeostasis through de novo synthesis and fatty acid uptake (Trefts et al., 2017). Notch inhibition upregulated genes functioning in fatty acid metabolism, such as the apolipoprotein encoding genes ApoA4, ApoC3, and ApoL1 (Supplemental Fig. 3P-R) and induced an increase in tumor fat content (Fig. 4F). Furthermore, JAG1 blockade induced expression of albumin (ALB), a classic hepatocyte marker (Supplemental Fig. 3S). Taken together, our results support a model in which JAG1-NOTCH2 inhibition targets sensitive HCCs by promoting differentiation of a progenitor/hybrid lineage to a mature hepatocyte phenotype, which is accompanied by cell cycle exit.

scRNA sequencing analysis reveals heterogeneity of LIV78 tumors and co-expression of JAG1 and NOTCH2

To clarify the cellular composition and the mechanism of efficacy following JAG1/NOTCH2 inhibition, we performed single-cell sequencing of LIV78 tumors isolated from control animals or following inhibition of JAG1 or NOTCH2. Analysis of human cells from all treatment groups yielded six distinct clusters of cells characterized by the expression of markers for: hepatoblasts with high extracellular matrix (ECM) and endoplasmic reticulum (ER) components (cluster 0), oxidative phosphorylation (cluster 1), cell cycle/proliferation (cluster 2), hepatocytes (cluster 3), hypoxia (cluster 4), and ribosomal components (cluster 5) (Fig. 5A). Given the heterogeneity of NOTCH2 activation in the sensitive HCCs (Fig. 2D and Supplemental Fig. 1B) and the separation of ligand- and receptor-expressing populations in other tumor types (Lim et al., 2017), we hypothesized that ligand- and receptor-expressing cells would fall into different clusters. Instead, we found many cells that co-express JAG1 and NOTCH2, with 55% of NOTCH2-positive cells also expressing JAG1 (Supplemental Fig. 4A, B). Ligand and receptor expressing cells co-contributed to multiple clusters, including the largest cluster (0, hepatoblasts), with common Notch targets expressed in a similar distribution (Supplemental Fig. 4C), highlighting cells active for Notch signaling. In contrast, clusters 1 and 5 showed few cells that expressed JAG1 or NOTCH2 (Supplemental Fig. 4A, B). These results suggest that sensitive tumors consist of heterogeneous populations of cells in distinct transcriptional states and that Notch signals emerge from ligand-receptor co-expressing cells.

Figure 5.

Notch inhibition triggers a cell state change regulated by a CEBPA-HNF4A transcriptional program

Given that Notch signaling commonly regulates binary cell fate decisions (Artavanis-Tsakonas et al., 1999; Koch et al., 2013), including self-renewal versus differentiation of progenitors, we investigated whether treatment affected the distribution of cells between the clusters. Notch blockade dramatically increased the fraction of cells most resembling mature hepatocytes (cluster 3, blue), with 18% of cells from treated tumors contributing to this cluster compared to 0.5% of cells from controls. Concomitantly, JAG1/NOTCH2 blockade reduced the cellular contribution to the hepatoblast-like progenitor population (cluster 0, red) from 70% in controls to 53% after treatment (Fig. 5B, C and Supplemental Fig. 4D-H). This shift from the progenitor to hepatocyte state confirms our bulk RNA- sequencing results (Fig. 4B) and further corroborates our hypothesis that Notch inhibition induces hepatocyte differentiation in the sensitive progenitor-like tumors.

To illuminate the transcriptional programs regulating this progenitor-to-hepatocyte transition, we used Latent Dirichlet Allocation (LDA) to infer continuous transcriptional programs (topics) within cells and correlated these topics with transcription factor (TF) activities for 289 TFs across all cells. Among the nine topics that describe the single-cell data (Supplemental Fig. 4I, J), Topic 3 most strongly associated with cells after treatment (Fig. 5D, Supplemental Fig. 4J). Two of the top three TF activities that best correlated with Topic 3 were HNF4A and CEBPA (Fig. 5E, F), both key regulators of hepatocyte differentiation (Flodby et al., 1996; Watt et al., 2003) that codependently localize on chromatin as cooperatively functioning TF partners (Stefflova et al., 2013). The Topic 3 enrichment for HNF4A activity did not stem from increases in HNF4A protein Supplemental Fig. 4K) or mRNA, the expression of which was only modestly modulated by Notch inhibition and indistinguishable across baseline clusters (Fig. 5G andSupplemental Fig. 4L). In contrast, enrichment for CEBPA activity paralleled increased CEBPA mRNA expression after Notch inhibition (Fig. 5G and Supplemental Fig. 4M), consistent with Notch suppressing CEBPA transcription in other cell types (De Obaldia et al., 2013).

To assess how Notch inhibition impacts chromatin accessibility and how chromatin changes relate to HNF4A and CEBPA activities, we performed Assay-for-Transposase-Accessible-Chromatin (ATAC) sequencing of LIV78 tumors, 72 hours after treatment with NOTCH2 blocking or control antibodies. NOTCH2 blockade profoundly changed the chromatin landscape, increasing or decreasing accessibility at 3934 and 3819 peaks, respectively (Fig. 5H, Supplemental Table 4), with HNF4A and CEBPA binding motifs enriched in chromatin regions with increased accessibility and the SOX family TF motifs in regions exhibiting decreased accessibility (Fig. 5I). HNF4A and CEBPA motifs localized in a characteristic pattern, with HNF4A binding close to the centers of opening regions and CEBPA depleted at these center positions but enriched up- and down-stream (Supplemental Fig. 4N).

We integrated chromatin accessibility data with gene expression measurements, using the BETA software to illuminate how these chromatin changes might relate to differences in gene expression. The increased accessibility of chromatin following NOTCH2 inhibition significantly correlated with increases in gene expression, with the greatest increases in genes related to hepatocyte maturation, which harbor HNF4A and CEPBA motifs in their promoter region. (Supplemental Fig. 4O, Supplemental Table 4). Conversely, decreased accessibility of chromatin correlated with down-regulation of genes associated with progenitor or bile duct cells, which had SOX transcription factor motifs present in the promoter region (Supplemental Fig. 4P, Supplemental Table4). Our TF activity and chromatin analyses thus lead to a model in which Notch inhibition leads to increased expression of CEBPA, thus enabling CEBPA to partner with HNF4A to jointly drive transcriptional programs and the underlying chromatin rearrangements that promote differentiation of progenitor-like tumor cells to a mature hepatocyte fate incompatible with tumor growth and maintenance (Fig. 6).

Figure 6.

Discussion

Discovery of a Notch-dependent, super responder HCC subtype

Our previous work in a liver cancer mouse model (Huntzicker et al., 2015b) prompted us to test the efficacy of selective Notch-targeting therapeutic antibodies in 45 liver tumor PDX models—to our knowledge, the largest in vivo liver cancer screen that has been reported. Our studies uncovered an HCC subtype that depends on Notch signaling, and the molecular and cellular mechanisms we describe provide a diagnostic foundation to help identify patient populations predicted to benefit from a Notch inhibitor.

Studies of Notch signaling highlight a complex, multifaceted role in cancer in general and a confusing, contradictory literature in liver cancer in particular. Notch has been described as both oncogenic and tumor suppressive, including in HCC (Dill et al., 2013; Viatour et al., 2011; Villanueva et al., 2012; Zhu et al., 2021). Moreover, studies disagree on the functional relevance of particular ligands and receptors, including NOTCH1 versus NOTCH2 (Dill et al., 2013; Kawaguchi et al., 2016; Villanueva et al., 2012). Some of the conflicting conclusions may reflect choices in model systems, reliance on receptor gain-of-function methods, tumor versus stromal signals, or differences amongst distinct HCC subtypes. The extent of tumor growth inhibition in vivo must also be noted. Whereas some dramatic responses to single-agent Notch inhibition have been observed, particularly in tumor models carrying Notch-activating mutations (Agnusdei et al., 2014; Ferrarotto et al., 2017; Stoeck et al., 2014; Y. Wu et al., 2010), regressions or stasis have been rarely achieved, including in our own work testing Notch inhibitors in multiple oncology indications (e.g., Choy et al., 2017). Our current study relies on a solid foundation, built by combining patient-derived xenografts, considered as “high-bar” preclinical models, with therapeutic antibodies highly selective for individual Notch ligands and receptors. While this experimental approach identified numerous models that displayed modest to moderate sensitivity (26-57% TGI), we also identified two exceptional models that were exquisitely sensitive to Notch inhibition. PDX tumors from LIV78 and Li1035 regressed after Notch blockade in a manner that was consistent, reproducible, rapid, and durable, thus displaying a “super responder” phenotype (Marx, 2015; Wheeler et al., 2021) unique to our experience. By using antibodies that selectively block individual ligands and receptors, we identified JAG1 and NOTCH2 as the functionally relevant ligand-receptor pair. While other Notch receptors, particularly NOTCH1 (Villanueva et al., 2012), have been proposed as HCC drivers, we found that NOTCH1 in the super responders was primarily expressed in tumor stromal cells (likely endothelial) and was not required for tumor maintenance. Our results thus highlight the same ligand-receptor pair that determines cholangiocyte versus hepatocyte fate during liver development and reveal the potential of JAG1 and NOTCH2 for targeted therapies in this HCC subclass.

Intriguing questions remain to explain what drives Notch signaling in the super responder subclass. Receptor mutations in the juxtamembrane negative regulatory region (NRR) that enable ligand-independent signaling or in the PEST domain degron that prolong signaling, as well as other genomic alterations that increase Notch signal strength, have been described in multiple cancer types; however, such changes have not been broadly observed in HCC, and we did not find any in our super responders. Similarly, JAG1 and NOTCH2 were expressed in the majority of the PDX models and thus were insufficient to mark a Notch-dependent phenotype. However, NOTCH2 activity, defined by nuclear localization of NICD2 and expression of Notch transcriptional targets, was exclusive in the LIV78 and Li1035 tumor cells. JAG1 blockade inhibited this NOTCH2 activity, demonstrating its ligand dependence. In our PDX models, the relevant JAG1 source could be mouse endothelial cells, which express JAG1 (Loomes et al., 1999), or human tumor cells, in which we detected JAG1 mRNA and protein. The newly developed aJ1.b6A9 antibody, which inhibits human but not mouse JAG1, blocked tumor cell NOTCH2 activity and induced regressions. Thus, we conclude that JAG1 induction of NOTCH2 signaling within or between human tumor cells is required for tumor maintenance whereas endothelial (mouse) JAG1 is likely not a major contributor. This tumor cell intrinsic ligand-dependent Notch activation mechanism clearly differs from the mutational activation mechanisms previously described.

Molecular analyses showed that the super responder models display fetal progenitor-like properties and associate with the S2/proliferative HCC subclass. Genetically engineered Notch activation in mouse liver cells positive for AFP, a fetal progenitor marker, was sufficient to drive formation of HCCs (Dill et al., 2013; Villanueva et al., 2012). These Notch-driven mouse tumors display transcriptional profiles that include developmentally-restricted imprinted genes like Igf2 and Meg3, and they classify as undifferentiated and proliferative, thus bearing striking resemblance to LIV78 and Li1035. The proliferative classification, particularly within the Hoshida-S2 subclass, has been previously linked to Notch activation and hepatitis infection (Goossens et al., 2015b and references therein), both of which we note in LIV78 and Li1035. Specifically, signatures derived from human HCCs with Hepatitis C etiology and mouse HCCs driven by Notch clustered together as S2/proliferative (Villanueva et al., 2012). While Notch has also been linked to the S1 subtype (Zhu et al., 2021), this link relied on a GSEA signature predominantly consisting of pathway components with undemonstrated relevance to pathway activity.

In anticipation of possible clinical tests, it is important to consider how to best identify the subset of HCC patients most likely to benefit from JAG1 or NOTCH2 inhibition. Super responder status was rare (∼5%) in our PDX collection. However, PDX models have undergone a powerful biological selection during the xenografting of human tumor samples in mice, and thus the frequency of Notch-dependent super responders in the human HCC population remains to be determined. Nevertheless, our results suggest that expression of fetal progenitor markers and the previously described S2 and proliferative subclassifications provide a valuable initial filter. We also note that some of our PDX models that displayed these traits responded weakly to JAG1 inhibition, revealing limitations to this diagnostic approach. Among our entire PDX panel, including the non-responsive models that transcriptionally co-clustered with the super responders, we also identified GAS7 and FGF9 as individual genes that were highly expressed only in the super responders. Furthermore, expression levels of these genes decreased following Notch inhibition, indicating that they are direct or indirect Notch targets. However, GAS7 and FGF9 transcript levels were insufficient to clearly stratify HCC patients using large HCC patient datasets. Understanding this discrepancy requires further study, perhaps including spatial methods (e.g., IHC, spatial transcriptomics) to discern GAS7 and FGF9 expression in cells from the tumor, normal liver, stroma, or immune compartment.

Mechanism of efficacy of NOTCH inhibition in tumors: a replay of cell fate regulation during development?

We found that the efficacy of JAG1 or NOTCH2 inhibition in sensitive hepatocellular carcinomas is based on the transition of a subset of tumor cells from a progenitor-like state to a mature hepatocyte-like phenotype. This trajectory and its regulation closely recapitulate key events during specification of the hepatocyte lineage during embryonic development; in absence of a Notch signal, bi-potent hepatoblasts give rise to hepatocytes, while cholangiocytes are formed only when JAG1-NOTCH2 signaling is active (McCright et al., 2002; Zong et al., 2009 and references therein). Moreover, several lines of evidence highlight the importance of Notch signaling activation in adult hepatocytes for promoting cholangiocyte programs and driving plasticity (Fan et al., 2012; Yimlamai et al., 2014; Zong et al., 2009).

Importantly, our single cell studies enabled in-depth analyses of the NOTCH-regulated molecular interactions that trigger differentiation along the hepatocyte lineage and tumor regressions. These studies identified two closely related tumor cell clusters, one containing progenitor-like cells and one with mature hepatocyte characteristics, which was primarily composed of cells that experienced Notch inhibition. By inferring transcription factor activities, we identified HNF4A and CEBPA as specifically associated with upregulated programs in the treatment-induced hepatocyte-like cluster. Interestingly, both transcription factors are enriched in the liver and essential for hepatocyte formation during embryonic development (Flodby et al., 1996; Parviz et al., 2003; Wang et al., 1995). Moreover, a network comprising both HNF4A, CEBPA, as well as HNF1A, another hepatocyte-lineage expressed transcription factor, and the cholangiocyte-specific factors HNF6 and HNF1B was found to be regulated by NOTCH signaling during lineage commitment of hepatoblasts (Tanimizu & Miyajima, 2004).

Hnf4a deficiency leads to aberrant hepatocyte morphology and function, including changes in cell adhesion that disrupt lobule architecture. Deletion of Cebpa in mice severely disturbs hepatocyte maturation, with parenchymal cells displaying mixed hepatocyte and biliary epithelial cell characteristics and forming pseudoglandular structures (Flodby et al., 1996). These structures, which notably resemble pseudoglandular hepatocellular carcinoma or tissue undergoing damage-repair, not only co-express markers for both lineages (including HNF4A and SOX9), but also show upregulated expression of Jag1 and Notch2 (Akai et al., 2014; Flodby et al., 1996; Yamasaki et al., 2006). Thus, the progenitor-like state we observed in Notch dependent HCC tumor models at baseline is highly similar to the one Cebpa loss creates during hepatogenesis.

Our single-cell studies of the LIV78 model showed that HNF4A transcript was expressed across all clusters, with transcript levels showing little to no treatment-dependent variation. In contrast, activity and expression of CEBPA strongly increased after Notch inhibition, with particularly high levels in cells contributing to a hepatocyte-like cluster that emerged following Notch-inhibition. Interestingly, CEBPA was shown to be negatively regulated by HES1(a bona fide Notch target gene) (De Obaldia et al., 2013) in the context of myeloid fate regulation. In chronic myelomonocytic leukemia, which results from an inactivating mutation in the y-secretase complex component nicastrin, Notch signaling has a tumor suppressive function, that is mediated through direct repression of CEBPA and PU.1 by HES1 (Klinakis et al., 2011). In the liver, HNF4A and CEBPA were shown to extensively co-localize on chromatin and loss of HNF4A reduces the ability of CEBPA to bind DNA and vice versa (Stefflova et al., 2013), rendering HNF4A necessary, but not sufficient for the development of mature hepatocyte features. The hypothesis that Notch inhibition through de-repressing CEBPA transcription allows chromatin engagement of both CEBPA and HNF4A and thus initiation of hepatocyte differentiation programs, is strongly supported by our ATAC sequencing data; chromatin regions that open following NOTCH2 inhibition are strongly enriched for CEBPA and HNF4A motifs which were found in close proximity of each other. Closing chromatin regions on the other hand are associated with motif enrichment for SOX family transcription factors, and HNF1B, re-emphasizing the role of active Notch signaling in maintaining cholangiocyte specific programs active and the progenitor-like tumor cells in a multipotent state.

Our study provides proof-of-concept for a clinical paradigm using therapeutic antibodies targeting JAG1 or NOTCH2 to promote a hepatocyte-like state incompatible with tumor maintenance, and identifies the key regulatory network underlying Notch inhibition-induced tumor cell differentiation. This, together with the finding that therapeutic antibodies against NOTCH2 are well tolerated in patients (Smith et al., 2019), holds great promise for the treatment of Notch-active HCCs and possibly other progenitor-like cancers.

Material and methods

Treatment and tumor growth assessment of liver cancer PDX models

Animals were maintained in accordance with the Guide for the Care and Use of Laboratory Animals (National Research Council 2011). Genentech is an AAALAC-accredited facility and all animal activities in this research study were conducted under protocols approved by the Genentech Institutional Animal Care and Use Committee (IACUC). Mice were housed in individually ventilated cages within animal rooms maintained on a 14:10-hour, light: dark cycle. Animal rooms were temperature and humidity-controlled, between 68 to 79°F (20.0 to 26.1°C) and 30 to 70% respectively, with 10 to 15 room air exchanges per hour.

To generate tumor-bearing animals for the mouse clinical trial and all subsequent studies, tumor fragments from 42 established PDX liver cancer models (comprised of 35 HCC models plus CCs and other liver tumor subtypes) were passaged from tumor-bearing animals and implanted sub-cutaneously (s.c.) into the right flank of either Balb/c nude (NU/J Foxn1^nu (JAX:002019), Jackson Laboratories) or NCr nude (CrTac:NCr-Foxn1^nu;Taconic Biosciences) immunocompromised mice to generate patient-derived xenografts (PDXs). Tumors were allowed to establish until they had engrafted and demonstrated consistent growth (typically between 80-280 mm³). Animals were then grouped into treatment cohorts (between n=5-10 across studies) to normalize across groups based upon tumor volume prior to initiation of treatment. During the mouse clinical trial, mice were treated with either IgG1 isotype anti-Ragweed control antibodies or anti-Jag1 (aJ1.b70) antibodies intraperitoneally (i.p.) at 10 mg/kg, weekly (QW) for 3 weeks.

For subsequent studies, either IgG1 isotype anti-gD, anti-Jag1 (aJ1.b70), anti-huma-Jag1 (ahJ1), anti-Notch1 (anti-NRR1), or anti-Notch2 (anti-NRR2) antibodies were used at the indicated doses (between 1.25 to 30 mg/kg), QW, i.p. either for efficacy assessment or generation of samples for analysis of pharmacodynamics. Animals were monitored routinely to assess health and body condition with tumor volumes and body weights being measured twice weekly. Animals demonstrating excessive body weight loss (>15%) or tumors in excess of 1500 mm3 or bearing ulcerations were routinely euthanized. In general all treatments were well-tolerated with the exception of anti-Notch1 (aNRR1), which has known and expected side-effects resulting in body weight loss.

Analyses and comparisons of tumor growth were performed using a package of customized functions in R (Version 3.4.2; R Foundation for Statistical Computing; Vienna, Austria), which integrates software from open source packages (e.g., lme4, mgcv, gamm4, multcomp, settings, and plyr) and several packages from tidyverse (e.g., magrittr, dplyr, tidyr, and ggplot2) (Forrest et al., 2020). Briefly, as tumors generally exhibit exponential growth, tumor volumes were subjected to natural log transformation before analysis. All raw tumor volume measurements less than 8 mm³ were judged to reflect complete tumor absence and were converted to 8 mm³ prior to natural log transformation. Additionally, all raw tumor volume measurements less than 16 mm³ were considered miniscule tumors too small to be measured accurately and were converted to 16 mm³ prior to natural log transformation. The same generalized additive mixed model (GAMM) was then applied to fit the temporal profile of the log-transformed tumor volumes in all study groups with regression splines and automatically generated spline bases. This approach addresses both repeated measurements from the same study subjects and moderate dropouts before the end of the study.

Estimates of group-level efficacy were obtained using the calculation of percent tumor growth inhibition (TGI). This value represents the percent difference in tumor growth between the treatment and reference groups according to each group fit baseline-corrected area under the curve (AUC). This AUC is obtained after back-transforming tumor volumes to the original scale and correcting for starting tumor burden. Positive values indicate antitumor effects, with 100% denoting stasis and values >100% denoting regression (negative values indicate a pro-tumor effect). The 95% confidence intervals are based on the fitted model and variability measures of the data obtained using parametric bootstrapping.

Generation and characterization of human-specific anti-Jagged1 antibody

For development of anti-huJagged1 hybridoma, five Balb/c mice (Charles River Laboratories International, Inc., Wilmington, MA, USA) were hyperimmunized, in each hind footpad at 3-4 day intervals, with recombinant soluble human Jagged1 ECD in RIBI adjuvant (Ribi Immunochem Research Inc. Hamilton, MT, USA). After 5 injections, serum titers were evaluated by standard enzyme-linked immunosorbant assay (ELISA) to identify mice with positive serum titers to huJagged1 ECD. B cells from popliteal lymph nodes from five mice were pooled and were fused with mouse myeloma cells (X63.Ag8.653; American Type Culture Collection, Manassas, VA, USA) by electrofusion (ECM 2001; Harvard Apparatus, Inc., Holliston, MA, USA). After 10-14 days, hybridoma supernatants were harvested and screened for antibody production with a huJagged1 binding ELISA. All specific clones were then re-screened for binding by FACS (293/3T3-huJagged1 cells) and for blocking Notch1-Jagged1 signaling by co-culture assay. Hybridoma clone, 6A9.7.3, showed high immunobinding and signal blocking after the second round of single cell per well subcloning (FACSAria cell sorter; BD Biosciences, San Jose, CA, USA) and was scaled up for purification in INTEGRA CELLine 1000 bioreactor (INTEGRA Biosciences AG, Zizers, Switzerland). Supernatant was purified by protein A affinity chromatography, sterile-filtered, and stored at 4°C in PBS. Isotype of the mAb was determined using the Isostrip Mouse mAb Isotyping Kit (Roche Applied Biosciences, Indianapolis, IN, USA). The isotype was determined to be IgG1, kappa. Clone 6A9.7.3 was sequenced and chimerized into huIgG1, kappa backbone.

Notch reporter assays were performed as previously described (Lafkas et al., 2015) using U87glioblastoma cells co-transfected with a Notch-response TP-1 (12XCSL) firefly luciferase reporter and a constitutively expressed Renilla luciferase reporter (pRL-CMV, Promega E2261). Antibodies were added with the ligand-expressing cells (NIH-3T3 cells stably transfected with human JAG1 or OP9 cells stably transfected with mouse Jag1).

Generation and analysis of bulk RNAseq data

Tumor tissue was harvested, snap frozen in liquid nitrogen and stored at -80C. Total RNA was extracted from 15-25mg intact or powderized tumor tissue using TRIzol^TM Reagent (Thermo Fisher Scientific Cat# 15596026) according to the manufacturer’s instructions. Samples were homogenized using a TissueLyser II (Qiagen) at 25k rpm for 3min. The optional 5min incubation to allow for complete dissociation of the nucleoprotein complex was performed. Following isopropanol precipitation, samples were loaded onto RNEasy Mini prep columns (Qiagen) and cleanup steps, including on-column DNAse digestion, were performed according to the RNEasy Kit’s instructions. Samples were eluted in 30ul RNAse-free water.

2-6 replicate samples were collected for each treatment condition. The concentration of RNA was determined using NanoDrop 8000 (Thermo Scientific) and the integrity of RNA was determined by Fragment Analyzer (Advanced Analytical Technologies). Approximately 500 ng of total RNA was used as an input for library preparation using TruSeq RNA Sample Preparation Kit v2 (Illumina). Size of the libraries was confirmed using High Sensitivity D1K screen tape (Agilent Technologies) and their concentration was determined by qPCR-based method using Library quantification kit (KAPA). The libraries were multiplexed and sequenced on Illumina HiSeq 2500 (Illumina). We obtained on average 30 million single-end RNA-seq reads (50bp) per sample.

Sequencing reads were processed using a XenoFilter pipeline to obtain human- and mouse-specific gene expression quantifications. Reads were first aligned to ribosomal RNA sequences to remove ribosomal reads. The remaining reads were aligned to the human reference genome (NCBI Build 38) using GSNAP (T. D. Wu & Nacu, 2010) version ‘2013-10-10’, allowing a maximum of two mismatches per 50 base pair sequence (parameters: ‘-M 2 -n 10 -B 2 -i 1 -N 1 -w 200000 -E 1 --pairmax-rna=200000 --clip-overlap’). The same procedure was used to map all reads to the mouse reference genome (NCBI Build 38). Reads that were multi-mappers in either of the two organisms were removed. Subsequently, concordant and unique reads with <= 3 mismatches to the human genome and >3 mismatches to the mouse genome were called human-specific. The reverse mismatch criteria were used to determine mouse-specific reads. Lastly, gene expression levels were quantified by calculating the number of reads mapped to the exons of each RefSeq gene using the HTSeqGenie R package. Transcript annotation was based on the RefSeq database NCBI Annotation Releases 104 and 106 for mouse and human, respectively. Read counts were scaled by library size, quantile normalized and precision weights calculated using the “voom” R package (Law et al., 2014). Subsequently, differential expression analysis on the normalized count data was performed using the “limma” R package (Ritchie et al., 2015) by contrasting anti-JAG1 or anti-Notch antibody treated samples with control samples. Gene expression was considered significantly different across groups if we observed an |log2-fold change| ≥ 1 (estimated from the model coefficients) associated with an FDR adjusted P-value ≤ 0.05. In addition, gene expression was obtained in the form of normalized Reads Per Kilobase gene model per Million total reads (nRPKM) as described previously (Srinivasan et al., 2016). Graphs were generated in the R environment or using Partek Software.

TCGA gene expression was obtained as raw counts that were estimated with the PanCanAtlas pipeline (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/ Expression_mRNA_Pipeline) using the human genome reference (NCBI Build 38) and Gencode (V22) gene models. Subsequently, read counts were scaled by library size, quantile normalized and precision weights calculated in the same fashion as for our in-house data (above). PCA was performed on the normalized data and UMAPs were produced using the first 1000 principal components.

Histology, immunohistochemistry and in situ RNA hybridization

Tumors were isolated, fixed in 10% neutral buffered formalin overnight, washed in two changes of ddH₂O and dehydrated through a series of increasing concentrations of ethanol and xylenes for paraffin embedding. Tissue sections were prepared at 5μm thickness and slides were incubated at 60C for 10 minutes, de-paraffinized using Xylenes and rehydrated through an ethanol gradient prior to any staining procedure.

For evaluation of morphology, sections were stained with Hematoxylin and Eosin (Fisher Scientific, # 23245637 and #6766010) using standard methods.

For immunohistochemical (IHC) and immunofluorescence (IF) staining, sections were boiled for 15 min in Target Retrieval solution (DAKO-S2369) using a pressure cooker, treated with 3% hydrogen peroxide (VWR-BDH7690-1) for 10 min, washed in PBS, and subjected to a 1 h blocking step in 5% bovine serum albumin (Sigma-A9647) in PBT (PBS with 0.1% Triton-X, Sigma-T9284). For IHC, incubation with primary antibodies from Cell Signaling Technologies against N1-ICD (Val1744, #4147), NOTCH2 (D76A6, #5732), or JAGGED1 (D47Y1R, #70109) diluted 1:2000 in blocking solution was performed overnight in a humidified chamber. Three 10 min washes in PBS were followed by incubation with Poly-HRP anti Rabbit IgG secondary antibody (PowerVision, PV6119) for one hour, two washes in PBS and one wash in PBT, each for 10 min and application of DAB substrate (25μl/2ml buffer, DAB+ Substrate Chromogen system, DAKO-K3468). Staining was performed in the dark for 2 minutes after which the slides were immediately immersed into tap water. To visualize the nuclei, the sections were counterstained for 40 seconds using Gill’s hematoxylin III diluted 1:5 in water. A brief bluing step in Scotts tap water substitute (Eng Scientific, #8442) was included prior to dehydration through ethanols and xylenes and slides were mounted in Permount (Fisher Scientific, #SP15).

For IF staining, sections were incubated in a 1:100 dilution of anti-BrdU antibody (Biorad, MCA2060), 1:100 dilution of anti-HNF4A (Perseus Proteomics, PP-H1415-00), or 1:1000 dilution of anti-Ki67 primary antibody (Thermo Scientific, RM-9106-50). Alexa Fluor 555 goat anti-rat or goat anti-rabbit, or Alexa Fluor 488 goat anti mouse secondary antibodies were used at a 1:250 dilution (Thermo Scientific A-21428, A-21434, and A-11001). Following nuclear staining with DAPI (Invitrogen, D1306), sections were mounted in Pro Long Gold antifade reagent (Invitrogen, P36934).

For in situ RNA hybridization, mRNAs were visualized by RNAscope® technology (Wang et al., 2012). Manual assays were performed using the RNAscope® 2.5 HD Reagent Kit-BROWN (Advanced Cell Diagnostics, #322300) according to the manufacturer’s instructions. Probes against FGF9 (Cat. No.439431) or GAS7 (custom order) as well as positive and negative control probes (310451, 310043) were designed by Advanced Cell diagnostics.

All tissue stainings were performed for at least 3 animals per treatment group.

Images were obtained on a Zeiss Axioscope 40 or Nanozoomer slide scanner.

Quantification of Immunofluorescence staining data

For quantification of BrdU+ or Ki67+ cells, 3 randomly chosen 20x frames were imaged per tumor. Image J software was used to count DAPI stained or marker-stained nuclei. First, images were converted to 8-bit format, a threshold was set and the “watershed” option was used to separate clustered nuclei. The “Analyze Particle” command was used with ‘‘size = 0.005-∝ inch 2’’ and ‘‘circularity = 0.00-∝ ”. Cell counts reported by ImageJ were expressed as % of marker-positive nuclei/DAPI positive nuclei. Graphs were generated using PRISM software. Each data point graphed is the average value from one mouse; data are presented as mean ± SEM.

Analysis of cell cycle phase by flow cytometry

PDX tumors were harvested, dissociated into single cells and subjected to depletion of CD45 and Ter19-postive populations as described above. Following two washes in ice-cold PBS, each sample was resuspended in 1ml ice-cold PBS and 3ml of ice-cold 100% ethanol were added to the tube, drop-wise while vortexing on slow speed to avoid clumping while fixing the cells. For proprio iodine staining, samples were pelleted, the ethanol removed, and following 2 washes with ice-cold PBS the cells were resuspended in 0.5ml PI/RNAse staining solution (BD550825) per 1x10^6 cells. After passing the cells through a 35um nylon mesh-strainer cap into FACS tubes (Corning Catalog #352235) they were incubated for 15 minutes at RT prior to analyzing on a BD Symphony Analyzer.

Generation, processing and analysis of single-cell RNAseq data

NCR nude mice bearing similarly sized LIV78 tumors (volumes between 300-900mm³) were injected intravenously with 30mg/kg bodyweight of either anti-Notch2 or anti-Jag1 or anti-gD control antibody 72hrs prior to tumor harvest. 2 samples (one with Notch inhibition and one with control antibody treatment) were processed in the same day.

The tumor cell isolation procedure was adapted from a liver cell isolation technique by Mederacke et al (Mederacke et al., 2015). Tumor-bearing mice were sacrificed and the tumors isolated and placed in a petri-dish containing 5mls of a total of 50ml of enzyme buffer solution (8000mg/L NaCl, 400mg/L KCl, 88.17mg/L NaH₂PO₄*H₂O, 2380mg/L

HEPES, 350mg/L NaHCO₃, 560mg/L CaCl₂ *2H₂O) with 25mg of Pronase (Sigma-Aldrich #P5147), 4.4 U of Collagenase D (Roche #11088882001), 32μM Actinomycin D (Chemodex #A0043) and 1% DNAse I (Roche #10104159001). Tumors were minced to small pieces with a razor blade, added to the remaining 45mls of Pronase/Collagenase/ DNAse/Actinomycin D solution in a 50 ml conical tube and digested for ∼1hr rotating in a hybridization oven at 42C. The digested tumors were filtered through a 70μm cell strainer into a new 50ml tube and centrifuged at 450xg for 10min at 4C. Following aspiration of all but 10ml of the supernatant, 120ul of DNAse I stock solution (100mg DNAse I in 50 ml Gey’s Balanced Salt Solution B (GBSS/B; Sigma-Aldrich #G9779)) were added and the cells were resuspended. The tube was filled up to 50ml with GBSS/B to wash the cells prior to centrifugation at 450xg for 10 min at 4C. Tumor cells were resuspended in 3% heat inactivated fetal bovine serum (Sigma-Aldrich #19G057-AH1) in PBS and subjected to removal of host immune and red blood cells using the EasySep^TM Mouse Mesenchymal Stem/Progenitor Cell Enrichment kit (Stem Cell Technologies #19771) according to the manufacturer’s instructions.

Following the depletion of CD45+ and Ter119+ cells, the samples were processed for single-cell RNA-seq (scRNAseq) as described previously (Long et al., 2019) using the Chromium Single Cell 3’ Library and Gel bead kit v2, following the manufacturer’s manual (CG00052 Chromium Single Cell 3 Reagent Kits v2 User Guide RevA; 10x Genomics). Cell density and viability of the single-cell suspensions were determined by Vi-CELL XR cell counter (Beckman Coulter). All of the processed samples had more than 95% viable cells. Cell density was used to impute the volume of single-cell suspension needed in the RT master mix. cDNAs and libraries were prepared following the manufacturer’s manual (10X Genomics). Libraries were profiled by Bioanalyzer High Sensitivity DNA kit (Agilent Technologies) and quantified using Kapa Library Quantification Kit (Kapa Biosystems, Wilmington, MA). Each library was sequenced in one lane of HiSeq 2500 (Illumina) following the manufacturer’s sequencing specification (10X Genomics).

Sequencing data were processed as described previously (Long et al., 2019). In short, reads were tallied and demultiplexed based on their association with cell-specific barcodes. Only barcodes with >10K associated reads were considered for further processing. For each cell, reads were mapped to the human and mouse reference genomes using the same GSNAP settings as for bulk RNA, respectively. The number of transcripts per gene was quantified based on unique molecular identifiers (UMI) using reads that overlapped exonic regions in sense direction. Subsequently, cells were classified as human (PDX) or mouse (stroma) if their mapping rate was higher than 42% to the respective genome. Only cells with a mitochondrial UMI fraction <0.25 were considered for analysis. Per-gene UMI counts were normalized by the total number of transcripts per cell and scaled by the median number of transcripts across all cells. Dataset alignment (between multiple processing dates), cell clustering, visualization, calculation of gene set scores and differential expression analysis were performed according to best practices using Seurat (v3.0.0.9000) (Stuart et al., 2019). Differential expression between treatment conditions was based on the likelihood-ratio test for single cell gene expression (McDavid et al., 2013).

To characterize cells based on the activity of transcriptional programs, we performed topic modeling using functions from the CountClust R package (Dey et al., 2017). We evaluated a range between 2 to 15 topics for their fit using the FitGoM function. A total of 9 topics were selected to represent the data after inspecting the AIC and BIC statistics. Top genes associated with each topic were identified by applying the ExtractTopFeatures function, which uses the relative gene expression profile of the GoM clusters and applies a KL-divergence based method to obtain a list of top features that drive each of the clusters.

To identify transcription factors (TFs) that drive specific transcriptional programs (characterized by individual topics) we quantified the transcription factor activities in individual cells and correlated them with the previously obtained topic weights. We used a total of 289 high-confidence regulons (comprised of TF downstream target genes) from the non-academic version of the Dorothea database (Garcia-Alonso et al., 2019; Holland, Szalai, et al., 2020; Holland, Tanevski, et al., 2020) and quantified their activity using AUcell (Aibar et al., 2017).

To perform differential abundance analysis of cell groups between treatment conditions, we used the MiloR package (Dann et al., 2022). In brief, milo uses a k-nearest neighbor (KNN) graph to define a set of representative neighborhoods, followed by counting of cells in each neighborhood by experimental condition. These counts are used for differential abundance testing using a negative binomial GLM framework followed by a weighted FDR procedure that takes the spatial overlap of neighborhoods into account.

Graphs were generated in the R environment or using Partek Software.

Gene set enrichment

We tested for the enrichment of particular gene categories to identify relevant biological processes associated with (1) differential expression or (2) transcriptional programs (topics), by using functions from the "clusterProfiler" R package (T. Wu et al., 2021) and the MSigDB gene set collection (Liberzon et al., 2011). In the case of differential expression, treated samples were contrasted with control samples, genes were ranked based on the π statistic (π=logFC * -log10(p-val)) (Xiao et al., 2014) and subsequently gene set enrichment was performed using the clusterProfiler::GSEA function or the GSEA4.1.0 software. In the case of transcriptional programs, we extracted the top 50 genes associated with specific topics and performed a hypergeometric test using the clusterProfiler::enricher function.

Generation and analysis of bulk ATAC-seq data

Tumors were isolated, snap frozen in liquid nitrogen and stored at -80C. Tumor tissue was pulverized using a Covaris CP02 CryoPrep Pulverizer. Flash frozen pulverized tumor samples were sent to Active Motif to perform the ATAC-seq assay.

The tissues were manually disassociated, isolated nuclei were quantified using a hemocytometer, and 100,000 nuclei were tagmented as previously described (Buenrostro et al. 2013), with some modifications based on (Corces et al. 2017) using the enzyme and buffer provided in the ATAC-Seq Kit (Active Motif). Tagmented DNA was then purified using the MinElute PCR purification kit (Qiagen), amplified with 10 cycles of PCR, and purified using Agencourt AMPure SPRI beads (Beckman Coulter). Resulting material was quantified using the KAPA Library Quantification Kit for Illumina platforms (KAPA Biosystems), and sequenced with PE42 sequencing on the NextSeq 500 sequencer (Illumina).

The ATAC-seq data were aligned to human reference genome (NCBI Build 38) using the ENCODE ATAC-seq pipeline (https://github.com/ENCODE-DCC/atac-seq-pipeline) with minimum p-value cutoff for peak calling set to 1e-6. The resulting peak lists were sorted and overlapping peaks were merged using bedtools (v2.26.0). The differential accessibility analysis was performed using the bioconductor package DiffBind (v3.2) in R (v4.1.0) (Ross-Innes et al., 2012; Stark & Brown, 2011). For differential accessibility analysis, peaks called in only one sample were filtered out. The read counts in the peak matrix were normalized using the dba.normalize() function of DiffBind with the parameter library = DBA_LIBSIZE_PEAKREADS. The statistical significance of differential accessible regions was estimated using the dba.analyze() function of DiffBind in its default setting. The motif enrichment analysis was performed using the AME (v5.4.1) tool of the MEME suite using the non-differential ATAC-seq peaks as background and HOCOMOCOv11_core_HUMAN motif database (Kulakovskiy et al., 2018; McLeay & Bailey, 2010). The motif localization analysis was performed using the CENTRIMO (v5.3.0) tool from the MEME suite (Bailey & MacHanick, 2012).

Supplemental Figures

Supplemental Figure 1 (goes with Main Figure2)

Supplemental Figure 2 (goes with Main Figure3)

Supplemental Figure 3 (goes with Main Figure4)

Supplemental Figure 4 (goes with Main Figure 5)

Acknowledgements

We wish to thank Nayab Adibi, Simona Hankeova, and Gloria Hernandez for helpful discussions. We are grateful to Heinrich Jasper, Louis Vermeulen, and Xin Ye for critical reading of the manuscript. We would like to thank the Genentech Necropsy, Histology and Pathology and Legal teams for their contributions and advise.

Additional information

Author contributions

K.S., R.P. and C.W.S. designed the study; C.W.S. supervised and directed the project. K.S. performed the experiments and data analysis; R.P., T.T.T. N. and M.A. processed, analyzed and interpreted sequencing data; C.C., A.Sh. and J.C. contributed to experiments and data analysis; J.G. performed in house PDX model studies; C.D.L.C, U.S., and M.M facilitated PDX model studies in house and at CROs and performed data analysis; Y.A.Y., J.C. and Z.M. performed RNA sequencing experiments; A.Sc. provided histopathological analysis and discussion; R.V. and Y.W. generated and validated antibodies used in this study; L.M. provided input for experiments, discussion and edits for the manuscript. C.M. performed experiments and data analysis and provided edits and critical input throughout this study; K.S., R.P., and C.W.S. cowrote the paper. All authors contributed to the manuscript and approve its content.