Abstract
Oesophageal adenocarcinoma (OAC) patients show poor survival rates and there are few targeted molecular therapies available. However, components of the receptor tyrosine kinase (RTK) driven pathways are commonly mutated in OAC, typified by high frequency amplifications of the RTK ERBB2. ERBB2 can be therapeutically targeted, but this has limited clinical benefit due to the acquisition of drug resistance. Here we examined how OAC cells adapt to ERBB2 inhibition as they transition to a drug resistant state. ERBB2 inhibition triggers widespread remodelling of the accessible chromatin landscape and the underlying gene regulatory networks. The transcriptional regulators HNF4A and PPARGC1A play a key role in this network rewiring. Initially, inhibition of cell cycle associated gene expression programmes is observed, with compensatory increases in the programmes driving changes in metabolic activity. Both PPARGC1A and HNF4A are required for the acquisition of resistance to ERBB2 inhibition and PPARGC1A is instrumental in promoting a switch to dependency on oxidative phosphorylation. Our work therefore reveals the molecular pathways that support the acquisition of a resistant state and points to potential new therapeutic strategies to combat cellular adaptation and ensuing drug resistance.
Similar content being viewed by others
Introduction
Cancer is predominantly caused by DNA mutations and genomic rearrangements. However, it is becoming increasingly clear that rewiring of the epigenetic landscape also plays a pivotal role in tumourigenesis [1,2,3,4]. This epigenetic rewiring leads to changes in gene expression programmes and the molecular pathways that are operational in the cell [5]. The molecular changes manifested by a cancer cell provide the opportunity for personalised treatment which is exemplified by the use of inhibitors like trastuzumab and lapatinib to treat cancer patients with amplifications of the RTK ERBB2 [6, 7]. However, drug resistance often arises, limiting the effectiveness of treatment. This is especially the case for OAC where trastuzumab and lapatinib have been shown to have relatively limited effects on patient survival [8, 9]. Drug resistance in OAC often arises due to the selection of cells containing compensatory mutations [10,11,12]. However, it is becoming increasingly recognised that changes to the epigenetic landscape can also play an important role in drug resistance, particularly in enabling the survival of “persistor” cells, which ultimately gather additional mutations to adopt a stable resistant form [13,14,15,16,17].
The incidence of OAC is rapidly expanding in the Western world while the survival rates remain poor [18]. A pre-malignant state known as Barrett’s oesophagus (BO) is thought to be the precursor to OAC [19]. Genome sequencing studies have revealed numerous mutational changes in the transition from BO to OAC but there are few high frequency recurrent oncogenic driver events [20,21,22]. However, at the pathway level, components of the receptor tyrosine kinase (RTK) driven pathways are frequently mutated in OAC (60–76%; [22, 23]), typified by relatively high frequency amplifications and mutations of the RTK ERRB2 (ranging from 18–32% tumours; [22, 23]). At the epigenetic level, the accessible chromatin landscape of BO is vastly different to the surrounding normal oesophageal tissue and this landscape is further altered during the transition to OAC [24, 25]. These chromatin changes are accompanied by alterations to the transcriptional regulatory networks, with many changes to transcription factor activity being common to BO and OAC, including transcription factors like HNF4A, GATA6, FOXA, HNF1B and PPARG [24, 26, 27]. Conversely, other transcription factors appear more dominant in OAC such as AP1 [28] or their regulatory activity is repurposed and directed to alternative transcriptional programmes as exemplified by KLF5 [25]. However, it is unclear how the deregulated RTK pathways impact on these gene regulatory networks.
Here, we have investigated how inhibition of the RTK ERBB2 influences gene regulatory networks as OAC cells acquire drug resistance. We demonstrate rapid and widespread remodelling of the accessible chromatin landscape in response to ERBB2 inhibition. This revealed changes in transcriptional regulatory activities, which converged on HNF4A and PPARGC1A and their influence on metabolic programmes. Both HNF4A and PPARGC1A are required for the acquisition of a drug resistant state and their target gene networks represent potential therapeutic vulnerabilities to prevent the emergence of drug resistance.
Results
ERRB2 amplified OAC cells develop resistance to lapatinib
Previous studies indicated that several gastro-oesophageal adenocarcinoma cancer cell lines harbour amplifications of the gene encoding ERBB2. By using ATAC-seq data, we validated this amplification event which encompasses the ERRB2 locus in all of the cancer lines but is absent in the HET1A and CPA cell lines, derived from normal and Barrett’s oesophageal tissue respectively (Supplementary Fig. S1A). These amplifications lead to variable levels of ERBB2 protein expression with the highest levels in the OE19 and NCI-N87 cell lines (Fig. 1A). To select cell lines which most closely resembled patient derived samples, we analysed the open chromatin landscapes of the cell line panel compared to three OAC patient biopsies which harbour ERBB2 amplifications. Principal component analysis demonstrated tight clustering of samples from normal or Barrett’s oesophagus, and a distinct looser cluster of the OAC samples (Fig. 1B). The OE19, ESO26, KYAE1 and OAC organoid WTSI-OESO_009 (CAM408; [29]) cluster together with these tumour samples whereas OE33 and NCI-N87 are more distantly associated. Furthermore, clustering based on Pearson’s correlations of the same data gave the same broad conclusions with OE33 cells being a clear outlier. Importantly, we also examined the expression of a set of transcription factors we previously associated with Barrett’s and OAC [24] and two markers of squamous epithelium, TP63 and PAX9, in a range of patient samples [30] and cell lines (Supplementary Fig. S1B). OE19, ESO26, and KYAE1 again clustered with OAC patient samples and exhibited expression of GATA6, FOXA2 and HNF4A. We therefore took OE19, ESO26, and KYAE1 cells forward for further analysis as representative examples of ERBB2-amplified OAC.
Next, we studied the response of OAC cells to treatment with the ERBB2/EGFR inhibitor lapatinib [31]. In all cell lines, lapatinib caused a reduction in the activating phosphorylation (Y1196) of ERBB2 and downstream ERK1/2 (Supplementary Fig. S1C). In addition to signalling via ERK, ERBB2 also signals via a parallel pathway involving AKT and similar reductions in AKT phosphorylation levels were observed with the exception of ESO26 cells where co-treatment with an AKT inhibitor was required to gain maximal inhibition due to the presence of activating mutations in PIK3CA in this cell line [10]. Genetic depletion of ERBB2 caused similar reductions in ERK and AKT phosphorylation in OE19 cells (Supplementary Fig. S1D). Lapatinib treatment caused dose-dependent reductions in growth of all the OAC lines (Fig. 1D; Supplementary Fig. S1E). Part of this reduced growth was due to increased levels of apoptosis as exemplified by OE19 cells (Fig. 1E). Having established an optimal concentration of lapatinib, we examined the response of OE19 cells over a 5 week period. ERBB2 phosphorylation remained inhibited for the duration of treatment but both AKT and ERK phosphorylation levels partially recovered after 5 weeks treatment (Fig. 1F). This recovery was accompanied by a resumption of proliferation as exemplified by the numbers of cells in S phase (Fig. 1G; Supplementary Fig. S1F) and increased cell numbers following an initial decrease following lapatnib treatment (Supplementary Fig. S1G). Thus by 5 weeks, the OE19 cells can be considered as resistant to ERBB2 inhibition with lapatinib and represent a good model for further investigation of the underlying resistance mechanisms.
ERBB2 inhibition causes widespread changes to the transcriptional regulatory landscape in OAC cells
To understand how ERRB2 inhibition rewires OAC cells along the path to resistance, we focussed on OE19 cells and examined changes to their gene expression profile and associated open chromatin landscape following treatment with lapatinib for 1, 7 and 35 days (Fig. 2A). RNA-seq analysis revealed over a thousand genes changing expression (>2 fold change, FDR < 0.05) after 1 day and 7 days (roughly equally split as up- and down-regulated), with a slightly lower number at 35 days (Fig. 2B; Supplementary Fig. 2A; Supplementary Table 1). Clustering based on PCA analysis and Pearson’s correlations indicated that by day 35, the transcriptome had begun to resemble the starting population (Supplementary Fig. S2B, C) which is also reflected in the recovery in expression levels of the differentially expressed genes (Fig. 2C, D; Supplementary Fig. S2D). Importantly, there is a substantial overlap in gene expression changes observed after 1 day of lapatinib treatment and RNAi-mediated ERBB2 depletion, which demonstrates the specificity of response to the inhibitor (Supplementary Fig. S2E). Gene ontology analysis of gene expression changes between sequential timepoints revealed a decrease in cell cycle gene expression which was further reduced at 7 days but reversed at 35 days as cells become resistant to ERBB2 inhibition (Fig. 2E) as exemplified by FOXM1, MYBL2 and MYC expression (Fig. 2F). These genes are generally typical of cluster 5 (Fig. 2C, D). Conversely, different processes are upregulated at day 1, including metabolic and signalling categories (Fig. 2E) which tend to be maintained in day 35 resistant cells (Fig. 2G) and is emphasised by the expression profiles of genes like DGAT2, EPHX2, IDH1 and PRKAB2 (Supplementary Fig. S2F). The changes in gene expression categories as cells become resistant are further emphasised by comparing GO terms relative to the parental OE19 cells, although this comparison demonstrates that the cell cycle is still suppressed at day 35 whereas various lipid metabolic processes are rapidly altered and maintained at higher levels throughout the time course (Supplementary Fig. S2G). Ingenuity pathway analysis (IPA) revealed the expected inhibition of the upstream regulator ERBB2 but also other cell cycle regulators like MYC and FOXM1 in resistant cells (Supplementary Fig. S2H). Collectively these data reveal that a substantial rewiring of the transcriptome occurs as cells progress to becoming resistant with initial decreases in cell cycle gene expression, accompanied by the sustained activation of metabolic processes.
To further probe the gene regulatory changes and uncover potential regulatory pathways, we next performed ATAC-seq to assess changes to the accessible chromatin landscape across the same time points towards a lapatinib resistant sate. Widespread changes were observed in the accessible chromatin landscape, which are initiated after 1 day and result in thousands of differentially accessible regions after 7 and 35 days (Fig. 2H; Supplementary Fig. S2I, J; Supplementary Table 2). The majority of the differentially accessible regions are found in intra and intergenic regions that represent potential enhancers (Supplementary Fig. S2K, L). There is a notable increase in promoter accessibility after 7 days whereas there is a tendency for putative enhancer regions to close (Supplementary Fig. S2L). Two broad groups of chromatin regions can be identified, which show general opening or closing across the time course (Fig. 2I) as exemplified by the MYC and the EPHX2 loci which mirrors changes in their gene expression (Fig. 2J). Indeed, more broadly, both promoters and enhancer regions show correlations between changes in accessibility and associated putative target gene expression (Supplementary Fig. S2M). Consistent with transcriptional changes, each set of regions is associated with genes linked to different biological processes, chiefly lipid metabolism and other metabolic processes for regions with increased accessibility and regulation of apoptosis and MAPK cascades for those closing (Fig. 2I; Supplementary S2N). The latter observation is consistent with the expected effects of ERBB2 inhibition.
Collectively, these data demonstrate widespread dynamic changes to the accessible chromatin landscape and associated transcriptome as cells adapt to inhibition of ERBB2 signalling through lapatinib treatment.
The remodelling of the accessible chromatin state is reversible
Having established the accessible chromatin landscape in lapatinib resistant cells, we next asked whether this is a stable state potentially promoted through the selection of rare mutational events or through the establishment of fixed epigenetic states in the cell population. We therefore profiled the open chromatin landscape of OE19 cells which had been cultured in lapatinib for 35 days (hence adopting a resistant state) followed by drug withdrawal at several timepoints between 1 and 14 days (Fig. 3A). Phosphorylation levels of ERBB2 and downstream ERK and AKT were re-established after 1–2 days of drug withdrawal, demonstrating reinstatement of ERBB2 signalling activity (Fig. 3B). Consistent with this, after 14 days drug withdrawal, the cells again became sensitive to lapatinib re-addition (Fig. 3C). PCA analysis showed a gradual drift of the open chromatin state of resistant cells back towards the parental OE19 cells (Fig. 3D; Supplementary Fig. S3A) which was further supported by Pearson’s correlation analysis which clustered the parental OE19 cells with the resistant cells which have been released from drug treatment for 14 days (Fig. S3B). Furthermore, when we examined the fate of the peaks which changed accessibility in the resistant cells, there was a rapid decay in number after drug withdrawal until barely any differences were detected between the parental cells and the resistant cells after 14 days (Fig. 3E, F). This is exemplified at the locus-specific level where closed chromatin is re-established around the CEBPA and PPARA loci whereas chromatin reopening is reinstated at the CCND3 and DUSP6 loci (Fig. 3G). A number of transcription factor motifs are enriched in peaks showing differentially increased or decreased accessibility after 35 days of drug treatment (see Fig. 4B, C). Interrogation of these enriched transcription factor motifs showed a gradual loss of binding sites for HNF4, GATA and CEBP factors in peaks exhibiting increased accessibility relative to parental cells following drug withdrawal from d35 cells (Fig. 3H, left). Conversely, there was reduced enrichment of the motifs for SOX, AP1 and ETS (ELF3 and ETV1) transcription factors in regions exhibiting decreased accessibility relative to parental cells following drug withdrawal (Fig. 3H, right). Thus, the underlying transcriptional regulatory networks of OE19 cells are also rapidly re-established following cessation of drug treatment.
Taken together these data show a rapid reactivation of the ERBB2 signalling pathway and almost complete reversion of the accessible chromatin landscape to the parental form following drug withdrawal. This argues against selection of a cellular subclone with a unique chromatin landscape and instead demonstrates that the underlying chromatin landscape is reprogrammed based on the intrinsic regulatory pathways available to the OAC cells, and that this is fully reversible upon re-activation of ERBB2 signalling.
Different transcriptional regulatory pathways are associated with chromatin opening and closing events
Having established the changes in chromatin accessibility during acquisition of resistance to ERBB2 inhibition, we next harnessed this data to uncover the transcription factor networks affected. To focus on the resistance mechanisms, we concentrated on putative enhancer elements (i.e., non-promoter regions) which showed increased or decreased accessibility in cells treated with lapatinib for 35 days (a time point at which cells are proliferating again and hence fully resistant to growth inhibition). More accessible regions already showed evidence of opening after 1 day of treatment whereas less accessible regions only closed substantially after 7 days (Fig. 4A; Supplementary Fig. S4A). Opening regions showed enrichment of motifs for a range of transcription factors which drive a gastro-intestinal programme in OAC including HNF4, GATA, FOXA, KLF and CEBP factors after 35 days (Fig. 4B; Supplementary Table 2; [24,25,26,27, 32]). These events were already initiated after 1 day treatment (Supplementary Fig. S4F). In contrast, closing regions contained over represented motifs for SOX transcription factors and members of the AP1 and ETS transcription factor families that we previously linked to OAC (Fig. 4C; Supplementary Table 2; [28, 33]). Closer examination of chromatin accessibility at HNF4A and AP1 motifs further emphasised the dynamic nature of the accessibility around these transcription factor binding sites (Fig. 4D). To examine whether any transcription factors that bind to these motifs exhibit changes in expression, we examined our RNA-seq data and found congruency between directionality of expression and motif presence (Fig. 4E, F). HNF4A, PPARA and several other nuclear hormone receptors show upregulation following lapatinib treatment, alongside CEBPA and several FOX family members (Fig. 4E). Similar increases in expression of the majority of these genes were observed upon depletion of ERBB2 (Supplementary Fig. S4B). Conversely, several SOX and ETS transcription factors (including ETV4/5), and AP1 proteins from the FOS subfamily are downregulated (Fig. 4F). PPARA upregulation is associated with the opening of multiple potential regulatory regions throughout its locus (Fig. 4G) whereas SOX13 downregulation coincides with the closing of several regions (Supplementary Fig. S4C).
We also asked whether these transcription factor networks were similarly affected in other OAC/GAC cell lines containing ERBB2 amplifications and focussed on the initial responses to lapatinib treatment for 1 day because an increase in accessibility centred on the TF motifs like HNF4 was already apparent after 1 day of treatment (Fig. 4D; Supplementary Fig. S4F). There was a moderate overlap in their overall accessibility profiles (Supplementary Fig. S4D) and in the identities of the specific peaks showing differential accessibility following lapatinib treatment (Supplementary Fig. S4E). In all cell lines, opening regions showed enrichment of binding motifs for HNF4A and FOX transcription factors whereas closing regions exhibited enrichment of AP1 motifs (Supplementary Fig. S4F). Thus, the changes to the transcriptional regulatory pathways that we observe in response to drug treatment are generally detected in a range of OAC cell lines.
To further examine which upstream pathways are operational in resistant cells, we performed IPA analysis on genes which were differentially activated after 35 days of lapatinib treatment. Consistent with our transcription factor motif analysis, HNF4A was identified as the most significant upstream regulator alongside a range of other nuclear hormone receptors and the co-activator PPARGC1A/PGC1α (Fig. 4H). Due to its known role as a co-activator for HNF4A [34,35,36] and other nuclear hormone receptors like ESRRA and PPARA/G [37, 38], we further investigated PPARGC1A and found that its expression was rapidly induced upon ERRB2 inhibition in OE19 cells and was subsequently maintained well above basal levels in resistant cells (Fig. 4I, J). We also detected a lapatinib-induced decrease in PPARGC1A phosphorylation at S571 (an inhibitory phosphorylation event; [39]), providing further activation potential. Similarly, lapatinib treatment caused increased PPARGC1A expression in a range of other OAC/GAC lines (Fig. 4K), consistent with increases we observed in PPARGC1A expression following ERBB2 depletion in OE19 cells (Supplementary Fig. S4B). Moreover, PPARGC1A levels are already elevated in samples from patient OAC tumours containing high ERBB2 levels compared to the precursor Barrett’s oesophagus (Fig. 4L). Several nuclear hormone receptors showed changed expression in tumours with high levels of ERBB2 with some being elevated (ESRRA) whereas others were either unchanged (HNF4A) or showed reduced levels (PPARA) (Supplementary Fig. S4G). There is therefore no uniform response to ERBB2 amplification in patients. Thus, a simple model whereby ERBB2 signalling directly suppresses nuclear hormone receptor expression and activity relative to the precursor Barrett’s state does not appear to be operational.
Collectively these data demonstrate that the lapatinib resistant state is characterised by the activation of a transcription factor network that comprises a set of gastro-intestinal transcription factors including HNF4A, and additional nuclear hormone receptors alongside the coactivator protein PPARGC1A.
HNF4A activity is enhanced and required for lapatinib resistance
To provide further mechanistic insights into how transcriptional regulatory networks influence the lapatinib resistant cell state, we first focussed on the nuclear hormone receptor HNF4A due to the enrichment of its binding motifs in accessible regions appearing in resistant cells after 5 weeks lapatinib treatment (Fig. 4B). We performed ChIP-seq analysis to ask whether HNF4A gained additional activities through changing its targets following lapatinib treatment (Supplementary Fig. S5A, B; Supplementary Table 3). This analysis was performed after 2 days lapatinib treatment as HNF4A motif enrichment in accessible chromatin (Supplementary Fig. S4D) and increased HNF4A expression (Fig. 4E) was already apparent after 1 day lapatinib treatment. We first examined the HNF4A signal at chromatin loci that increased in accessibility one day after lapatinib treatment and found a general increase in HNF4A binding occurred after lapatinib treatment (Fig. 5A, top; Supplementary Fig. S5C). This is consistent with the increases in HNF4A expression we observe (Fig. 4E). Furthermore, this increased binding was accompanied by increased levels of the active histone mark H3K27ac, consistent with an activating event. In contrast, closing regions showed only low levels of HNF4A binding in parental cells with little change in binding triggered by lapatinib treatment, although inactivation was suggested by the decreases in H3K27ac (Fig. 5A, bottom; Supplementary Fig. S5C). We also partitioned the HNF4A binding regions into those with increased (ERBB2 suppressed) or decreased (ERBB2 activated) occupancy in lapatinib treated cells (Supplementary Table 3). HNF4A peaks which show increased occupancy also exhibited increases in both H3K27ac levels and chromatin opening following lapatinib treatment (Fig. 5B; top) whereas HNF4A peaks showing reduced occupancy also demonstrated reduced levels of acetylation but little change in chromatin opening (Fig. 5B; bottom). These results indicate that elevated levels of HNF4A binding in response to lapatinib treatment are associated with chromatin activation events and these accompany gene expression changes as exemplified by the IMPA1 (Fig. 5C), KDM5B (Supplementary Fig. S5D), loci. Indeed, there is a more widespread strong correlation between the changes in HNF4A binding activity and the changes of H3K27ac observed following lapatinib treatment (Supplementary Fig. S5E). Further analysis of genes neighbouring HNF4A binding regions that exhibit increased occupancy in response to lapatinib treatment, showed that these activation events correlate with increased gene expression under the same conditions (Fig. 5D). These HNF4A binding regions are rich in HNF4A binding motifs (70% of regions), consistent with direct binding but we also observed enrichment for GATA motifs suggesting combinatorial control at a subset of sites (Fig. 5E). In contrast, in regions showing decreased HNF4A occupancy, HNF4A binding is less associated with HNF4A motifs (25% of regions) and instead, AP1 is the highest ranked motif, suggesting potential indirect recruitment of HNF4A to these regions (Fig. 5F). Functionally, regions which both opened (d1), and harboured HNF4A binding regions (d2) following lapatinib treatment are close to genes associated with GO terms encompassing a variety of mitochondrial related and metabolic processes, consistent with a potential role in metabolic reprogramming (Supplementary Fig. S5F). To establish the relevance of these regulatory events in drug resistance, we depleted HNF4A (Supplementary Fig. S5G) and monitored the re-emergence of cell proliferation following lapatinib treatment. Loss of HNF4A blunted the re-establishment of cell growth and resistance in OE19 (Fig. 5G; Supplementary Fig. S5H) cells. In contrast, HNF4A loss did not cause reduced growth of parental OE19 cells over a 5 day period in the absence of lapatinib treatment (Supplementary Fig. S5I) although we have not studied the longer term consequences of this.
These results therefore demonstrate that ERBB2 inhibition triggers enhanced HNF4A binding to chromatin which is associated with chromatin and gene activation events. These events are important for the emergence of resistance, potentially through metabolic reprogramming.
PPARGC1A promotes metabolic changes during acquisition of the resistant state
To further understand the transcription factor networks driving the resistant state, we turned to PPARGC1A, which has previously been shown to act as a potential coactivator protein for HNF4A [34,35,36]. Our results have implicated PPARGC1A in driving the chromatin changes in the resistant state (Fig. 4H), therefore we performed ChIP-seq in OE19 cells to identify the binding sites for PPARGC1A (Supplementary Fig. S6A; Supplementary Table 4). Unexpectedly, we found that PPARGC1A and HNF4A exhibited virtually mutually exclusive binding locations (Fig. 6A; Supplementary Fig. S6B, C). These two proteins therefore likely contribute to lapatinib resistance through different regulatory activities and target genes.
Next, we analysed PPARGC1A in more detail to uncover its regulatory activities. PPARGC1A binding was detected at 1,438 sites in OE19 cells that are largely associated with promoter regions (Fig. 6B). These regions exhibited a general increase in accessibility upon lapatinib treatment, which was most apparent at day 35 (Fig. 6A, C; Supplementary Fig. S6D). In addition, the PPARGC1A binding regions are associated with the active histone mark H3K27ac, in both parental and lapatinib treated cells (Fig. 6A) and exhibit slight increases in acetylation after drug treatment (Supplementary Fig. S6D). We tested a panel of these binding sites in lapatinib treated cells and found little change in PPARGC1A occupancy, although significant increases were found at the ACO2 and ESRRA loci (Fig. 6D). Furthermore, we examined genome-wide binding of PPARGC1A following lapatinib treatment and found very little change in PPARGC1A occupancy, with little evidence for new binding events following ERBB2 inhibition (Supplementary Fig. S6E, F). This suggested that lapatinib induces a change in PPARGC1A activity rather than promoting de novo binding. Indeed, reductions in inhibitory phosphorylation events on PPARGC1A were observed indicating a likely increase in pre-bound PPARGC1A activity (Fig. 4J). To gain insights into the potential transcription factors that might recruit PPARGC1A, we searched for enriched DNA motifs and uncovered the NRF1 motif in both promoters and potential distal regulatory regions in the vast majority of binding regions (Fig. 6E, F). In keeping with the lack of overlap with the HNF4A ChIP-seq data, we did not observe co-occurrence of the HNF4A binding motif. We took advantage of ChIP-seq data for NRF1 from other cell lines and observed a large overlap with PPARGC1A binding regions (Supplementary Fig. S6B). NRF1 has previously been implicated in controlling mitochondrial function and its metabolic activities [40,41,42] and consistent with this, two of the top GO terms for genes associated with PPARGC1A binding regions are “mitochondrial organisation” and “TCA cycle and respiratory electron transport” (Supplementary Fig. S6G).
To uncover the regulatory consequences of PPARGC1A binding, we next depleted PPARGC1A (Supplementary Fig. S7A and 7B) and performed RNA-seq on OE19 cells after lapatinib treatment for 24 h (Supplementary Fig. S7C). PPARGC1A depletion mainly led to downregulation of gene expression with 86 genes exhibiting significantly reduced levels (Fig. 6G; 1.4 fold change, FDR < 0.05; Supplementary Table 4). By analysing all significantly downregulated genes, we observed a large number of direct targets for PPARGC1A that are associated with PPARGC1A binding peaks (Supplementary Fig. S7D; left; Supplementary Table 4) including ESRRA and ACO2 (Fig. 6H). The number of direct target genes increased following lapatinib treatment (Supplementary Fig. S7D; right; Supplementary Table 4). Many of these directly activated targets (36.6%) show upregulation in response to lapatinib treatment as exemplified by ESRRA and ACO2 (Fig. 6H; Supplementary Fig. S7E). GO term analysis of the directly activated PPARGC1A target genes uncovered several terms associated with mitochondrial function and aerobic respiration (Fig. 6I) and many of these genes are key components of the TCA cycle or the electron transport chain (Fig. 6J). We therefore examined mitochondrial activity by measuring the oxygen consumption rate in parental OE19 cells and cells treated with lapatinib in the presence and absence of PPARGC1A depletion (Fig. 6K). Basal and maximal respiration rates were dampened by lapatinib treatment and maximal respiration was further reduced by PPARGC1A depletion. In contrast, the glycolytic rate as assessed by the extracellular acidification rate, was unaffected by PPARGC1A depletion (Supplementary Fig. S7F). However, a lapatinib-mediated decrease in glycolytic rate was observed, suggesting a general decrease in metabolic activity.
Having established that the PPARGC1A-driven pathway is one of the predominant transcriptional pathways activated in OAC cells following lapatinib treatment, we asked whether this transcriptional regulator plays a role in acquired resistance. We created two different shRNA constructs to deplete PPARGC1A (Supplementary Figs. S5G,S7G) and tested the re-emergence of cell growth following extended treatment with lapatinib. Control OE19 cells transduced with scrambled shRNA constructs restarted proliferation after 30 days lapatinib treatment, but cells depleted of PPARGC1A failed to resume growth (Fig. 6L). This reduced growth was also apparent at earlier times of lapatinib treatment (Supplementary Fig. S5H). In contrast, OE19 cells grown in the absence of lapatinib showed no reductions in growth following PPARGC1A depletion (Supplementary Fig. S5H). Reciprocally, we performed a gain of function experiment in the gastric adenocarcinoma cell line NCI-N87 which harbours an ERBB2 amplification. This cell line is sensitive to ERBB2 inhibition (Supplementary Fig. S1E) and has a similar open chromatin profile to several OAC cell lines and patient samples (Fig. 1C) but has relatively low levels of PPARGC1A expression (Supplementary Fig. S7H). We inserted a doxycycline inducible PPARGC1A transgene (Fig. 6M) and tested growth in the presence of lapatinib. High levels of PPARGC1A promoted enhanced cell growth in the presence of lapatinib (Fig. 6N). To see if we could prevent the emergence of resistance by targeting the metabolic pathways regulated by PPARGC1A, we treated OAC cells with lapatinib plus the IDH2 inhibitor enasidenib. IDH2 is an important enzymatic component of the citric acid cycle and enasidenib has been previously suggested as a mechanism to target PPARGC1A-driven mitochondrial respiration [43], and whilst enasidenib is more specific to mutant IDH2, it can still inhibit WT IDH2 [44] that is present in parental OE19 cells [45]. Treatment of OAC cells with enasidenib in the absence of lapatinib only had a slight, generally insignificant effect on cell growth over a two week period (Supplementary Fig. S7I, J). In contrast, co-treatment of OE19, ESO26 and KYAE1 cells with lapatinib and enasidenib impaired the emergence of resistance after 5 weeks of treatment (Fig. 6O), and this inhibitory effect of enasidinib on the growth of OE19 cells was observed consistently throughout a timecourse of lapatinib treatment (Supplementary Fig. S7K). Thus, re-purposing of clinically used mitochondrial inhibitors could potentially prevent the emergence of resistance to RTK inhibitors in OAC.
Collectively, our data demonstrate that a PPARGC1A-driven pathway becomes activated by ERBB2 inhibition with lapatinib. Regions bound by this coactivator protein become more accessible and associated target genes are activated following lapatinib treatment, leading to changes in mitochondrial activity and oxidative phosphorylation. This regulatory activity is essential for the resumption of proliferative activity during the acquisition of resistance.
Discussion
Drug resistance is a major clinical problem for cancer treatment, especially when targeting signalling pathways [46]. This applies to OAC where compensatory mutational changes have been shown to be acquired to permit resistance to the ERBB2 inhibitor trastuzumab [10,11,12]. However, there are many ways for tumours to acquire resistance, and recently it has been shown that changes to the regulatory chromatin landscape can permit resistance in several scenarios, including endocrine resistance in breast cancer [13], BET inhibitor resistance in acute myeloid leukaemia [47], and PDGFRA inhibitor resistance in glioblastoma [48]. Here we have focussed on RTK pathways driven by ERBB2 amplifications in OAC and demonstrate that widespread rewiring of the regulatory chromatin landscape is a major driver towards acquiring drug resistance.
The integration of chromatin changes with transcriptomic changes enabled us to uncover transcriptional regulatory networks that are rewired as cells respond to drug treatment and develop resistance. Shortly after treatment, chromatin closing is instigated which corresponds to transcriptional programmes driven by AP1 and ETS family transcription factors. These families of transcription factors have previously been shown to play important roles in OAC [26, 28, 33], and lapatinib acts to shut these down. Conversely, chromatin opening causes these programmes to be replaced with a network driven by a set of transcription factors that are usually associated with early intestinal development, typified by HNF4, KLF5, FOXA and GATA factors (Fig. 6P). This set of transcription factors has previously been shown to be operational in both Barrett’s and OAC [24,25,26,27, 32] and suggests that ERBB2 activation leads to at least partial suppression of their activity in OAC cells. Alternatively, additional compensatory pathways may be activated in response to drug treatment which further enhance the activity of this set of transcription factors. Further analysis also implicated the transcriptional co-activator PPARGC1A in regulating the gene expression programmes that arise as cells acquire resistance (Fig. 6P). These transcription factors are associated with changes in gene expression changes linked to various metabolic programmes, suggesting that metabolic reprogramming allows cells to persist in the presence of ongoing lapatinib treatment. Although we only measured switching to oxidative phosphorylation after short term lapatinib treatment, it appears likely that this effect persists after longer term treatment as several OXPHOS pathway genes remain upregulated after 5 weeks lapatinib treatment (e.g., ACO2, OGDH, COX15 and UQCRC2; Supplementary Fig. S7E). Our findings in OAC are consistent with work in ERBB2-amplified breast cancer cells where metabolic adaptations to lapatinib treatment occur which shift the cells towards mitochondrial energy metabolism [49]. However, in breast cancer these changes are driven by the transcription factor, ERRα, demonstrating that resistance to the same inhibitor can occur in two different cancers through reprogamming of different pre-existing gene regulatory networks.
By focussing on HNF4A and PPARGC1A, we provide more detailed analysis of their downstream regulatory activities. PPARGC1A has previously been shown to be a co-activator for several nuclear hormone receptors, including HNF4A [34,35,36]. However, unexpectedly we found little evidence for chromatin co-occupancy. This is reflected in the metabolic programmes they are implicated in, with PPARGC1A primarily associated with oxidative phosphorylation whereas HNF4A controls various lipid metabolic processes. Interestingly, both converge on programmes controlling mitochondrial functions, suggesting that mitochondrial changes may be key in promoting resistance. Both HNF4A and PPARGC1A activity is required for acquiring resistance, underlying the importance of their downstream regulatory programmes. Indeed, when we inhibit oxidative phosphorylation by targeting IDH2, resistance is impaired, providing more direct evidence that this metabolic switch is important for cancer development. These findings are consistent with studies in other cancer types where PPARGC1A-driven programmes play an important role in drug resistance in melanomas [50] and glioblastoma [51]. In both cases, PPARGC1A-dependent metabolic reprogramming is elicited following inhibition of either a RTK (MET) or a downstream pathway component (BRaf), demonstrating a common response to RTK pathway inhibition in rewiring metabolic programmes in OAC and these cancers. Furthermore, HNF4A loss in pancreatic cancer also drives metabolic reprogramming and correlates with a metabolic switch to dependence on glycolytic activity [52], suggesting a more widespread role for HNF4A in this context.
After 5 weeks of lapatinib treatment, OAC cells begin to proliferate again despite ERBB2 being inactive (as defined by lack of Y1196 phosphorylation), suggesting that other compensatory signalling pathways must be activated which is reflected by the re-activation of ERK. This may be due to other RTKs being activated, potentially through remodelling their gene regulatory regions. However, this is not a stable state and is fully reversible after drug withdrawal and the open chromatin landscape reverts to its initial state within days. This has important implications for considering treatment regimes, as drug holidays could allow cells to be re-sensitised to ERBB2 inhibition, thereby reducing the selective pressure to select for additional mutational events. Our work has additional therapeutic implications as selective targeting of the persistor cells may prevent transition to a resistant state. One route would be through direct inhibition of HNF4A or PPARGC1A, or alternatively through targeting their regulatory programmes as exemplified by IDH2 inhibition with enasidenib, a clinically approved drug [44, 53]. As depletion of HNF4A and PPARGC1A in the presence of lapatinib reduces cell growth and metabolic reprogramming occurs early in the treatment time course, early intervention appears feasible and leads to a persistent change in cell viability.
In summary, we have demonstrated that OAC cells undergo dramatic remodelling of their accessible chromatin landscape following ERBB2 inhibition. This chromatin remodelling revealed that the HNF4A and PPARGC1A regulated transcriptional regulatory networks play a pivotal role in promoting the emergence of drug resistant cancer cells and represent potential therapeutic targets to enhance the efficacy of ERBB2 inhibition strategies. However prior to embarking on clinical trials, further pre-clinical investigation of the efficacy of combinatorial drug treatment should be undertaken in animal models. More generally, such an approach could be combined with inhibitors of other RTK pathway components given the high frequency of mutational events in this pathway in OAC [22, 23].
Methods
Cell culture and treatments
OE19, CP-A, HET1A, NCI-N87 and HEK293T cells were purchased from ATCC and were authenticated by STR profiling. KYAE1 and ESO26 cells were purchased from ECACC. Cell lines were routinely tested for mycoplasma. OE19, ESO26 and NCI-N87 cells were cultured in RPMI 1640 (ThermoFisher Scientific, 52400) supplemented with 10% foetal bovine serum (ThermoFisher Scientific, 10270) and 1% penicillin/streptomycin (ThermoFisher Scientific, 15140122). KYAE1 cells were cultured in 1:1 RPMI 1640:F12 (Thermo Fisher, 11765054) supplemented with 10% foetal bovine serum and 1% penicillin/streptomycin. HEK293T and HET1A cells were cultured in DMEM (ThermoFisher Scientific, 22320-022) supplemented with 10% foetal bovine serum. CP-A cells were cultured in keratinocyte serum free media (ThermoFisher Scientific, 17005042) supplemented with 10% foetal bovine serum, 5 ng/mL EGF (ThermoFisher Scientific, 17005042), 50 μg/mL bovine pituitary extract (ThermoFisher Scientific, 17005042) and 1% penicillin/streptomycin. Cell lines were cultured at 37 °C, 5% CO2 in a humidified incubator. For inhibitor treatment, cells were seeded at a density of 2 × 104 cells/cm2. 24 h after seeding, cells were treated with inhibitors [ERBB2/EGFR, Lapatinib (Selleckchem, S1028); AKT, MK-2206 (Selleckchem, S1078); IDH2, Enasidenib (Selleckchem, S8205) were reconstituted in dimethyl sulfoxide (DMSO)] or vehicle control. For long term treatments, inhibitor and media was replenished every 72 h. To ensure all timepoints were treated equally, 24 h prior to the endpoint fresh media/inhibitor was added to cells. For drug withdrawal experiments lapatinib resistant cells were first generated by treatment with lapatinib for 35 days. After 35 days of drug treatment cells were trypsinised, counted and seeded at a density of 2 × 104 cells/cm2 in the continued presence of lapatinib. 24 h after seeding, media containing lapatinib was removed, cells were washed 3X with PBS to wash out inhibitor before adding fresh media to the cells.
Organoid culture
WTSI_OESO-009 (CAM408; [29]) was a gift from Rebecca Fitzgerald. Organoids were cultured based on protocols described previously [29], except organoids were cultured in IntestiCult™ Organoid Growth Medium (STEMCELL Technologies, 06010). Organoids were cultured in 6-well plates. To passage organoids, media was removed and organoids were washed with PBS. 1 mL PBS was added and organoids were dissociated from BME-2 (Cultrex, 3533-005-02) by pipetting. The organoid suspension was centrifuged at 1000 RCF for 5 min and the supernatant was discarded. A single cell suspension was then form by the addition of 1 mL TrypLE™ Express (Gibco, 12609-013), and cells were transferred to a 1.5 mL tube and incubated for 15 min, 1000 rpm, 37 °C. Single cells were centrifuged at 1800 RCF for 5 min and the supernatant was removed by pipetting. Single cells were re-suspended in 200 μL 10 mg / mL BME-2 and 20 μL cells / BME-2 suspension were then pipetted into a 6-well plate, before incubating for 15 min at 37 °C, 5% CO2 to enable BME-2 to polymerise. 1.5 mL IntestiCult™ Organoid Growth Medium was then added to culture the cells.
For MTS growth assays, organoids were formed in 3D culture as described and then dissociated from BME-2 by pipetting. Organoids were then counted using a haemocytometer and 5 × 103 organoids were seeded into 96-well plates. For ATAC-seq, organoids were dissociated from BME-2 into single cells using TrypLE. ATAC-seq libraries were then generated.
Cell growth, cell cycle and apoptosis assays
MTS growth assays (Promega, G3580) were performed in 96-well plate format according to the manufacturer’s protocol. Absorbance readings were taken at 490 nm on a SPECTROstar Nano Micoplate Reader (BMG LABTECH).
Crystal violet assays were performed by fixing cells with 4% paraformaldehyde for 10 min. Cells were stained using 0.1% crystal violet (Sigma-Aldrich, HT90132) for 30 min at room temperature. Plates were rinsed with water and left to dry before solubilising dye in 10% acetic acid for 10 min at room temperature with gentle shaking. Absorbance readings were taken at 570 nm on a SPECTROstar Nano Micoplate Reader (BMG LABTECH). Data was uniformly transformed so that the mean of the control sample was represented as 100%.
For cell cycle analysis, media was collected from cells and cells were dissociated using trypsin, using the collected media to quench trypsin. Cells were collected by centrifugation, washed with PBS and then fixed by the addition of 70% ethanol, pre-cooled to −20 °C. Fixed cells were then stored at −20 °C until required. Cells were then pelleted and washed with PBS before the addition of 400 μL 50 μg/mL propidium iodide (Sigma, P4864). Cells were then analysed by the University of Manchester Flow Cytometry Core Facility on a BD Biosciences LSRFortessaTM. Data was analysed using ModFit LTTM software to determine the percentage of cells in G0/G1, S or G2/M phase.
For apoptosis assays, 2 × 105 cells were seeded into 6-well plates. 24 h later, fresh culture media was added, and cells were treated with 500 nM lapatinib or vehicle control. 30 μM propidium iodide (Sigma, P4864) was added to culture medium to measure apoptosis. Cells were imaged every 20 min for 72 h using an Incucyte ZOOM (ESSEN Bioscience), maintained at 37 °C, 5% CO2. The number of apoptotic cells was determined by the number of red fluorescent cells and the data was exported to Prism 8 (GraphPad).
siRNA transfection
4 × 105 cells were reverse transfected with 25 pmol siRNA using LipofectamineTM RNAiMAX transfection reagent (ThermoFisher Scientific, 13778150) according to the manufacturer’s instructions. Cells were seeded into 6-well plates. SMART-pool siRNAs for control non-targeting siRNA (Dharmacon, D-001810-10-0020), siERBB2 (Dharmacon, L-003126-00-0005) and siPPARGC1A (Dharmacon, L-005111-00-0005) were used.
Lentiviral vectors, production and transduction
The PPARGC1A (PGC1α) gene was amplified from the vector pcDNA myc PGC-1 alpha (Addgene, 10974) using primers containing BamHI and EcoRI cloning sites (PPARGC1A_BamHI_F, PPARGC1A_EcoRI_R, see Supplementary Table 5) and sub-cloned into pENTR1A (Invitrogen, A10462). A 3’ 3x FLAG tag was then inserted by site directed mutagenesis using a Q5® Site-Directed Mutagenesis Kit (NEB, E0554S) and the primers PPARGC1A_FLAG_SDM_F and PPARGC1A_FLAG_SDM_R (see Supplementary Table 5). PPARGC1A-3xFLAG was then cloned into pINDUCER20 (Addgene, 44012) using Gateway™ LR Clonase™ II Enzyme mix (Invitrogen, 11791-020), forming the vector pINDUCER20-PPARGC1A-3xFLAG.
shRNA vectors were created by annealing shRNA overlapping oligonucleotides and then cloned into pLKO.1 (Addgene, 10878). shRNA oligonucleotides sequences are detailed in Supplementary Table 5. Scramble shRNA pLKO.1 vector (Addgene, 1864) was used as a control.
Lentivirus was produced as described previously [54]. Briefly, 3 × 106 HEK293T cells were seeded in T75 flasks. The following day, HEK293T cells were transfected with 2.25 μg psPAX2 (Addgene, 12260), 1.5 pMD2.G (Addgene, 12259) and 3 μg target vector using PolyFect transfection reagent (Qiagen, 301107). Media containing virus was collected both 48 and 72 h post-transfection and viral particles were precipitated using PEG-itTM Virus Precipitation Solution (System Biosciences, LV810A-1). To transduce cells with virus cells were treated with both virus (MOI 0.5–1.0) and 5 μg / mL Polybrene infection reaction (EMD Millipore, TR-1003). For pLKO.1 vectors, polyclonal cells were selected using 500 ng/mL puromycin (Sigma P7255) for 2 weeks; for pINDUCER20 vectors, polyclonal cells were selected using 250 μg/mL G418 (ThermoFisher Scientific, 10131027) for 2 weeks.
Western blots
Cells were lysed in RIPA buffer (150 mM NaCl, 1% IGEPAL CA-630, 0.5% sodium deoxycholate, 0.1% SDS, 50 mM Tris pH 8.0, 1 mM EDTA) supplemented with protease inhibitors (Roche, 11836170001). Protein concentration was determined by BCA assay (Pierce, 23227). 5x SDS loading buffer (235 mM SDS, 50% glycerol, 0.005% bromophenol blue, 10% β-mercaptoethanol, 210 mM Tris-HCl pH 6.8) was then added to protein lysates to a 1x concentration and then incubated at 90 °C for 10 min. Samples were then analysed by SDS-PAGE on 8% polyacrylamide gels using a PageRuler™ Prestained Protein Ladder (Thermo Scientific, 26616). Proteins where then transferred onto a nitrocellulose membrane (GE Healthcare, 10600002) and blocked using Odyssey® Blocking Buffer (LI-COR Biosciences, P/N 927-40000). Primary antibodies used: anti-ERBB2 (Thermo Fisher, MA5-14057, 1:1,000), anti-phospho-ERBB2 (Cell Signalling Technologies, 6942, 1:5000), anti-ERK (Cell Signalling Technologies, 4695 S, 1:1,000), anti-phospho-ERK (Cell Signalling Technologies, 9106 S, 1:2,000), anti-AKT (Cell Signalling Technologies, 2920 S, 1:2,000), anti-phospho-AKT (Cell Signalling Technologies, 4060 S, 1:2,000), anti-HNF4A (R&D Systems, PP-H1415-00, 1:1,000), anti-PPARGC1A (Novus Biologicals, NBP1-04676, 1:1,000), anti-phospho-PPARGC1A (Novus Biologicals, AF6650, 1:1,000), anti-GFP (Santa Cruz, sc-8334, 1:2,000). For the anti-PPARGC1A antibody, membranes were blocked and the primary antibody was diluted in 5% (w/v) milk powder in PBS. Secondary antibodies used: anti-rabbit (LI-COR Biosciences, 926-32213, 1:10,000) and anti-mouse (LI-COR Biosciences, 926-32210, 1:10,000). The membranes were visualised using a LI-COR Odyssey® CLx Infrared Imager. Western blots were quantified using Empiria studio v1.1.
RNA extraction and RT-qPCR
Total RNA was extracted from cells using a RNeasy Plus RNA extraction kit (Qiagen, 74136) according to the manufacturer’s protocol. RT-qPCR reactions were run using a QuantiTect SYBR Green RT-qPCR kit (Qiagen, 204243) on a Qiagen Rotorgene Q. Relative copy number of transcripts was determined from a standard curve and then normalised by the expression of RPLP0 control gene. Primer pairs are listed in Supplementary Table 5.
Seahorse metabolism assays
1 × 104 OE19 cells were reverse transfected with siRNAs and seeded into 96-well plates (Agilent, 102601-100). 24 h post-transfection cells were treated with 500 nM lapatinib or vehicle control. 24 h after drug treatment normal culture medium (RPMI 1640) was changed to Seahorse XF RPMI assay medium (Agilent, 103681-100) and a mitochondrial stress test (Agilent, 103015-100) was carried out according to the manufacturer’s protocols using a Seahorse XFe96 analyser (Agilent) to measure oxygen consumption and extracellular acidification rates. Following completion of the mitochondrial stress test, a crystal violet assay was performed and oxygen consumption and extracellular acidification results were normalised by crystal violet 570 nm absorbance readings. Concentrations of inhibitors used: 1.5 μM oligomycin, 1 μM FCCP and 0.5 μM rotenone/antimycin A.
ATAC-seq processing and analysis
Two biological replicates were sequenced per condition. Omni-ATAC-seq was performed following published protocols [55]. Cells were dissociated from plates using trypsin (Gibco, 25300-054) and 2 × 105 cells were collected and centrifuged at 500 RCF for 5 min. Cells were then washed with PBS and centrifugation repeated. The omni-ATAC-seq protocol [55] was then followed until the isolation of nuclei. Nuclei were centrifuged at 500 RCF, 4 °C for 10 min and re-suspended in 10 μL nuclease free water (Ambion, AM9937). Nuclei were then counted using a haemocytometer, and 5 × 104 nuclei were used for the transposition reaction. Libraries were size selected using Ampure XP beads (Beckman Coulter Agencourt, A63881) using a two-sided selection (0.5x reaction volume and 1.25x reacton volume) and eluted in 12 μL nuclease free water. ATAC-seq libraries were then sequenced by the University of Manchester Genomic Technologies Core Facility on a HiSeq 4000 System (Illumina).
Initial processing of ATAC-seq was performed as described previously [28]. Reads were trimmed using Trimmomatic v0.32 [56] and aligned to the human genome (GRCh37, hg19) using Bowtie2 v2.3.0 [57] with the following options: -X 2000 -dovetail. Using SAMtools v1.9 [58], only mapped reads(>q30) were retained. Reads mapping to blacklisted regions were removed using BEDtools v2.27.1 [59]. Duplicates were then marked using Picard (https://broadinstitute.github.io/picard/). Peaks were called using MACS2 v2.1.1 [60] with the following parameters: -q 0.01, -nomodel-shift -75 -extsize 150 -B -SPMR. For the drug withdrawal timecourse all samples were processed as described above but the -SPMR option was not used when calling peaks with MACS2.
To create a union peakset for each experiment, biological replicates were checked for concordance (r > 0.90) and then merged into a single alignment file. Peaks were then re-called using MACS2 and the top 50,000 most significant peaks from each condition were retained. Peak summits were extended + /− 250 bp using BEDtools slop and then a union peakset of all conditions was created using HOMER v4.9 [61] mergePeaks.pl using the -d 250 parameter. This union peakset was then used to identify differentially accessible regions by counting reads mapping to peaks (i.e., accessible regions) in individual replicates using featureCounts v1.6.2 [62]. Read counts were then used in DESeq2 v1.14.1 [63] to call differentially accessible regions. Typically, an FDR (adjusted p-value) < 0.05 and 2 fold linear fold change cut off was used to define differential regions. To identify transcription factor binding motifs enriched in differentially accessible regions, peaks were separated into promoter or non-promoter associated peaks based on whether the peaks were −2.5 kb / + 0.5 kb from the TSS using BEDtools intersectBed. De novo or known motif enrichment was then performed in non-promoter peaks using HOMER v4.9 [61].
ChIP-qPCR and ChIP-seq processing and analysis
ChIP-qPCR and ChIP-seq was performed as described previously [25]. Primers for ChIP-qPCR are listed in Supplementary Table 5. 5 × 106 nuclei and 2.5 μg antibody were used for transcription factor/co-activator immunoprecipitation, and 2 × 106 nuclei and 2 μg antibody were used for histone marks. 50 μL of protein A or G Dynabeads were used for immunoprecipitation (Invitrogen, 10002D and 10004D). Antibodies used: anti-H3K27ac (abcam, ab4729), anti-HNF4A (R&D Systems, PP-H1415-00), anti-PPARGC1A (Novus Biologicals, NBP1-04676). For ChIP-seq, two biological replicates were sequenced per condition and replicates were checked for concordance (r > 0.80). Whilst spike in control chromatin was supplemented to chromatin, analysis of results showed that ‘reads in peaks’ normalisation of HNF4A and H3K27ac ChIP-seq data was more appropriate because global changes to the levels of HNF4A and H3K27ac were not observed.
Reads were trimmed using Trimmomatic v0.32 [56] and aligned to the human genome (GRCh37, hg19) using Bowtie2 v2.3.0 [57]. Mapped reads (>q30) were retained using SAMtools v1.9 [58]. Reads mapping to blacklisted regions were removed using BEDtools v2.27.1 [59]. Duplicates were then marked using Picard (https://broadinstitute.github.io/picard/). Peaks were called using MACS2 v2.1.1, using input DNA as control [60]. For differential binding analysis peak summits were extended + /- 250 bp using BEDtools slop and a union peakset was created by merging peaks from each sample using HOMER mergePeaks.pl (-d 250). Reads mapping to peaks in the union peakset were then counted using featureCounts v1.6.2 and analysed in DESeq2 v1.14.1 to call differentially bound sites (FDR < 0.1). Note that for PPARGC1A ChIP-seq, one of the day 2 lapatinib treatment replicates showed significantly lower quality, so we used a single replicate for comparing to day 0 merged data to maximise the chances of detecting new binding events. Day 2 replicate 1 showed high overall correlation with day 0 merged data (r = 0.8) demonstrating the quality of the data.
RNA-seq processing and analysis
Three biological replicates were sequenced per condition. Total RNA was extracted from cells using a RNeasy Plus RNA extraction kit (Qiagen, 74136). An on-column DNase digest (Qiagen, 79254) was performed according to the manufacturer’s protocol. RNA-seq libraries were generated using a TruSeq stranded mRNA library kit (Illumina, RS-122-2001) and sequenced by the University of Manchester Genomic Technologies Core Facility on a HiSeq 4000 System (Illumina).
Reads were trimmed using Trimmomatic v0.32 [56] and aligned to the human genome (GRCh37, hg19, RefSeq transcript annotation) using STAR v2.3.0 [64]. Gene expression counts were obtained using featureCounts v1.6.2 [62] and differentially expressed genes were identified using DESeq2 v1.14.1 using FDR < 0.05 [63]. For results from the lapatinib treatment timecourse, differentially expressed genes were filtered to remove genes in which no timepoint had an FPKM value > 1. Metascape [65] was used for gene ontology analysis of differentially expressed genes. Ingenuity Pathway Analysis [66] was used to predict upstream regulators.
Patient tumour samples in the OCCAMS dataset expressing ERBB2 at high levels (ERBB2HIGH) were defined by ERBB2 expression levels greater than the median +2 SD.
Bioinformatics and data visualisation
To visualise data, the bigwig files from an experiment were normalised using the reciprocal of scale factors obtained in DESeq2. The bigwigCompare tool from deepTools v3.1.1 [67] was then used to scale individual replicate bigwigs. Bigwigs for each condition were then merged using UCSC bigwigmerge and the resulting bedgraph files were converted to bigwigs for visualisation in IGV v2.7.2 [68]. Heatmaps of epigenomic data were generated using deepTools. Tag density plots were generated in deepTools and the data was then plotted in Microsoft Excel. Heatmaps of expression data and Pearson correlations were generated using Morpheus (https://software.broadinstitute.org/morpheus/). Peaks were annotated to genes using HOMER for the nearest gene model or GREAT [69] for the basal plus extension model. Principal component analysis was performing using the prcomp function in R v3.6.0. Euler diagrams were generated using the R Eulerr library.
Statistical tests
To determine whether an overlap of genes or peaks is statistically significant a hypergeometric test was used, using the dhyper function in R v3.6.0. Other statistical tests were performed using GraphPad Prism v8. All T-tests were two-tailed.
Datasets
All data was obtained from ArrayExpress, unless stated otherwise. Cell line RNA-seq data was obtained from: EBI, EGAD00001001357 (Cancer Genome Project, cancer.sanger.ac.uk[45]); NCI-N87 RNA-seq data, Sequence Read Archive SRP091839[70]; OE33 RNA-seq data, E-MTAB-5175 [28]. HET1A ATAC-seq data, E-MTAB-6931 [24]. CP-A ATAC-seq data, E-MTAB-8994 [25].
Patient tissue RNA-seq data was obtained from: E-MTAB-4054 [30] and the OCCAMS consortium (EGAD00001007496). Human tissue ATAC-seq data was obtained from: E-MTAB-5169 [28], E-MTAB-6751 [24], E-MTAB-8447 [25] and The Cancer Genome Atlas OAC ATAC-seq data were obtained from the GDC data portal (portal.gdc.cancer.gov; [71]).
OE19 HNF4A ChIP-seq data was obtained from E-MTAB-6858 [24]. OE19 siERBB2 RNA-seq data was obtained from E-MTAB-8579 [25]. OE19 H3K27ac ChIP-seq data was obtained from NCBI SRA SRP201335 [26]. NRF1 ChIP-seq data was obtained from ENCODE: HepG2, ENCSR853ADA; K562, ENCSR837EYC; MCF7, ENCSR135ANT [72].
Data availability
Sequencing data have been deposited in ArrayExpress. OE19 lapatinib treatment timecourse ATAC- and RNA-seq: E-MTAB-10302, E-MTAB-10304. Lapatinib treatment of WTSI-OESO_009, ESO26, KYAE1 and NCI-N87 cells ATAC-seq: E-MTAB-10306, E-MTAB-10307, E-MTAB-10310, E-MTAB-10313. Lapatinib withdrawal timecourse ATAC-seq: E-MTAB-10314. OE19 siPPARGC1A RNA-seq: E-MTAB-10317. OE19 HNF4A, PPARGC1A and H3K27ac ChIP-seq data: E-MTAB-10319. OE19 PPARGC1A lapatinib treated ChIP-seq data: E-MTAB-11300.
References
Somerville TDD, Xu Y, Miyabayashi K, Tiriac H, Cleary CR, Maia-Silva D, et al. TP63-mediated enhancer reprogramming drives the squamous subtype of pancreatic ductal adenocarcinoma. Cell Rep. 2018;25:1741–55.
Flavahan WA, Drier Y, Johnstone SE, Hemming ML, Tarjan DR, Hegazi E, et al. Altered chromosomal topology drives oncogenic programs in SDH-deficient GISTs. Nature 2019;575:229–33.
Allis CD, Jenuwein T. The molecular hallmarks of epigenetic control. Nat Rev Genet. 2016;17:487–500.
Jones PA, Issa JP, Baylin S. Targeting the cancer epigenome for therapy. Nat Rev Genet. 2016;17:630–41.
Bradner JE, Hnisz D, Young RA. Transcriptional addiction in. Cancer Cell. 2017;168:629–43.
Blackwell KL, Burstein HJ, Storniolo AM, Rugo H, Sledge G, Koehler M, et al. Randomized study of Lapatinib alone or in combination with trastuzumab in women with ErbB2-positive, trastuzumab-refractory metastatic breast cancer. J Clin Oncol. 2010;28:1124–30.
Cobleigh MA, Vogel CL, Tripathy D, Robert NJ, Scholl S, Fehrenbacher L, et al. Multinational study of the efficacy and safety of humanized anti-HER2 monoclonal antibody in women who have HER2-overexpressing metastatic breast cancer that has progressed after chemotherapy for metastatic disease. J Clin Oncol. 1999;17:2639–48.
Bang YJ, Van Cutsem E, Feyereislova A, Chung HC, Shen L, Sawaki A, et al. Trastuzumab in combination with chemotherapy versus chemotherapy alone for treatment of HER2-positive advanced gastric or gastro-oesophageal junction cancer (ToGA): a phase 3, open-label, randomised controlled trial. Lancet 2010;376:687–97.
Hecht JR, Bang YJ, Qin SK, Chung HC, Xu JM, Park JO, et al. Lapatinib in combination with capecitabine plus oxaliplatin in human epidermal growth factor receptor 2-positive advanced or metastatic gastric, esophageal, or gastroesophageal adenocarcinoma: TRIO-013/LOGiC–A randomized phase III trial. J Clin Oncol. 2016;34:443–51.
Kim J, Fox C, Peng S, Pusung M, Pectasides E, Matthee E, et al. Preexisting oncogenic events impact trastuzumab sensitivity in ERBB2-amplified gastroesophageal adenocarcinoma. J Clin Invest. 2014;124:5145–58.
Janjigian YY, Sanchez-Vega F, Jonsson P, Chatila WK, Hechtman JF, Ku GY, et al. Genetic predictors of response to systemic therapy in esophagogastric cancer. Cancer Disco. 2018;8:49–58.
Wang DS, Liu ZX, Lu YX, Bao H, Wu X, Zeng ZL, et al. Liquid biopsies to track trastuzumab resistance in metastatic HER2-positive gastric cancer. Gut 2019;68:1152–61.
Bi M, Zhang Z, Jiang YZ, Xue P, Wang H, Lai Z, et al. Enhancer reprogramming driven by high-order assemblies of transcription factors promotes phenotypic plasticity and breast cancer endocrine resistance. Nat Cell Biol. 2020;22:701–15.
Sharma SV, Lee DY, Li B, Quinlan MP, Takahashi F, Maheswaran S, et al. A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell 2010;141:69–80.
Ramirez M, Rajaram S, Steininger RJ, Osipchuk D, Roth MA, Morinishi LS, et al. Diverse drug-resistance mechanisms can emerge from drug-tolerant cancer persister cells. Nat Commun. 2016;7:10690.
Hammerlindl H, Schaider H. Tumor cell-intrinsic phenotypic plasticity facilitates adaptive cellular reprogramming driving acquired drug resistance. J Cell Commun Signal. 2018;12:133–41.
Shen S, Vagner S, Robert C. Persistent cancer cells: the deadly survivors. Cell 2020;183:860–74.
Coleman HG, Xie S-H, Lagergren J. The epidemiology of esophageal adenocarcinoma. Gastroenterology 2018;154:390–405.
Peters Y, Al-Kaabi A, Shaheen NJ, Chak A, Blum A, Souza RF, et al. Barrett oesophagus. Nat Rev Dis Prim. 2019;5:35.
Weaver JMJ, Ross-Innes CS, Shannon N, Lynch AG, Forshew T, Barbera M, et al. Ordering of mutations in preinvasive disease stages of esophageal carcinogenesis. Nat Genet. 2014;46:837–43.
Stachler MD, Camarda ND, Deitrick C, Kim A, Agoston AT, Odze RD, et al. Detection of mutations in barrett’s esophagus before progression to high-grade dysplasia or adenocarcinoma. Gastroenterology 2018;155:156–67.
Frankell AM, Jammula S, Li X, Contino G, Killcoyne S, Abbas S, et al. Oesophageal cancer clinical and molecular stratification (OCCAMS) consortium, fitzgerald RC. The landscape of selection in 551 esophageal adenocarcinomas defines genomic biomarkers for the clinic. Nat Genet. 2019;51:506–16.
The Cancer Genome Atlas Research Network, Kim J, Bowlby R, Mungall AJ, Robertson AG, Odze RD, et al. Integrated genomic characterization of oesophageal carcinoma. Nature 2017;541:169.
Rogerson C, Britton E, Withey S, Hanley N, Ang YS, Sharrocks AD. Identification of a primitive intestinal transcription factor network shared between esophageal adenocarcinoma and its precancerous precursor state. Genome Res. 2019;29:723–36.
Rogerson C, Ogden S, Britton E, OCCAMS Consortium, Ang YS, Sharrocks AD. Repurposing of KLF5 activates a cell cycle signature during the progression from a precursor state to oesophageal adenocarcinoma. Elife 2020;9:e57189.
Chen L, Huang M, Plummer J, Pan J, Jiang YY, Yang Q, et al. Master transcription factors form interconnected circuitry and orchestrate transcriptional networks in oesophageal adenocarcinoma. Gut 2020;69:630–40.
Ma S, Zhou B, Yang Q, Pan Y, Yang W, Freedland SJ, et al. A transcriptional regulatory loop of master regulator transcription factors, PPARG, and fatty acid synthesis promotes esophageal adenocarcinoma. Cancer Res. 2021;81:1216–29.
Britton E, Rogerson C, Mehta S, Li Y, Li X, Fitzgerald RC, et al. Open chromatin profiling identifies AP1 as a transcriptional regulator in oesophageal adenocarcinoma. PLoS Genet. 2017;13:e1006879.
Li X, Francies HE, Secrier M, Perner J, Miremadi A, Galeano-Dalmau N, et al. Organoid cultures recapitulate esophageal adenocarcinoma heterogeneity providing a model for clonality studies and precision therapeutics. Nat Commun. 2018;9:2983.
Maag JLV, Fisher OM, Levert-Mignon A, Kaczorowski DC, Thomas ML, Hussey DJ, et al. Novel aberrations uncovered in barrett’s esophagus and esophageal adenocarcinoma using whole transcriptome sequencing. Mol Cancer Res. 2017;15:1558–69.
Xia W, Mullin RJ, Keith BR, Liu LH, Ma H, Rusnak DW, et al. Anti-tumor activity of GW572016: a dual tyrosine kinase inhibitor blocks EGF activation of EGFR/erbB2 and downstream Erk1/2 and AKT pathways. Oncogene 2002;21:6255–63.
Pan J, Silva TC, Gull N, Yang Q, Plummer JT, Chen S, et al. Lineage-specific epigenomic and genomic activation of oncogene HNF4A promotes gastrointestinal adenocarcinomas. Cancer Res. 2020;80:2722–36.
Keld R, Guo B, Downey P, Cummins R, Gulmann C, Ang YS, et al. PEA3/ETV4-related transcription factors coupled with active ERK signalling are associated with poor prognosis in gastric adenocarcinoma. Br J Cancer. 2011;105:124–30.
Yoon JC, Puigserver P, Chen G, Donovan J, Wu Z, Rhee J, et al. Control of hepatic gluconeogenesis through the transcriptional coactivator PGC-1. Nature 2001;413:131–8.
Rhee J, Inoue Y, Yoon JC, Puigserver P, Fan M, Gonzalez FJ, et al. Regulation of hepatic fasting response by PPARgamma coactivator-1alpha (PGC-1): requirement for hepatocyte nuclear factor 4alpha in gluconeogenesis. Proc Natl Acad Sci USA. 2003;100:4012–7.
Charos AE, Reed BD, Raha D, Szekely AM, Weissman SM, Snyder M. A highly integrated and complex PPARGC1A transcription factor binding network in HepG2 cells. Genome Res. 2012;22:1668–79.
Mootha VK, Handschin C, Arlow D, Xie X, St Pierre J, Sihag S, et al. Erralpha and Gabpa/b specify PGC-1alpha-dependent oxidative phosphorylation gene expression that is altered in diabetic muscle. Proc Natl Acad Sci USA. 2004;101:6570–5.
Vega RB, Huss JM, Kelly DP. The coactivator PGC-1 cooperates with peroxisome proliferator-activated receptor alpha in transcriptional control of nuclear genes encoding mitochondrial fatty acid oxidation enzymes. Mol Cell Biol. 2000;20:1868–76.
Li X, Monks B, Ge Q, Birnbaum MJ. Akt/PKB regulates hepatic metabolism by directly inhibiting PGC-1alpha transcription coactivator. Nature. 2007;447:1012–6.
Evans MJ, Scarpulla RC. NRF-1: a trans-activator of nuclear-encoded respiratory genes in animal cells. Genes Dev. 1990;4:1023–34.
Wu Z, Puigserver P, Andersson U, Zhang C, Adelmant G, Mootha V, et al. Mechanisms controlling mitochondrial biogenesis and respiration through the thermogenic coactivator PGC-1. Cell 1999;98:115–24.
Cam H, Balciunaite E, Blais A, Spektor A, Scarpulla RC, Young R, et al. A common set of gene regulatory networks links metabolism and growth inhibition. Mol Cell. 2004;16:399–411.
De Vitto H, Bode AM, Dong Z. The PGC-1/ERR network and its role in precision oncology. NPJ Precis Oncol. 2019;3:9.
Yen K, Travins J, Wang F, David MD, Artin E, Straley K, et al. AG-221, a first-in-class therapy targeting acute myeloid leukemia harboring oncogenic IDH2 mutations. Cancer Disco. 2017;7:478–93.
Tate JG, Bamford S, Jubb HC, Sondka Z, Beare DM, Bindal N, et al. COSMIC: the catalogue of somatic mutations in cancer. Nucleic Acids Res. 2019;47:D941–7.
Boumahdi S, de Sauvage FJ. The great escape: tumour cell plasticity in resistance to targeted therapy. Nat Rev Drug Disco. 2020;19:39–56.
Bell CC, Fennell KA, Chan YC, Rambow F, Yeung MM, Vassiliadis D, et al. Targeting enhancer switching overcomes non-genetic drug resistance in acute myeloid leukaemia. Nat Commun. 2019;10:2723.
Liau BB, Sievers C, Donohue LK, Gillespie SM, Flavahan WA, Miller TE, et al. Adaptive chromatin remodeling drives glioblastoma stem cell plasticity and drug tolerance. Cell Stem Cell. 2017;20:233–46.
Deblois G, Smith HW, Tam IS, Gravel SP, Caron M, Savage P, et al. ERRα mediates metabolic adaptations driving lapatinib resistance in breast cancer. Nat Commun. 2016;7:12156.
Haq R, Shoag J, Andreu-Perez P, Yokoyama S, Edelman H, Rowe GC, et al. Oncogenic BRAF regulates oxidative metabolism via PGC1α and MITF. Cancer Cell. 2013;23:302–15.
Zhang Y, Nguyen TTT, Shang E, Mela A, Humala N, Mahajan A, et al. MET inhibition elicits PGC1α-dependent metabolic reprogramming in glioblastoma. Cancer Res. 2020;80:30–43.
Brunton H, Caligiuri G, Cunningham R, Upstill-Goddard R, Bailey UM, Garner IM, et al. HNF4A and GATA6 loss reveals therapeutically actionable subtypes in pancreatic cancer. Cell Rep. 2020;31:107625.
Stein EM, DiNardo CD, Pollyea DA, Fathi AT, Roboz GJ, Altman JK, et al. Enasidenib in mutant IDH2 relapsed or refractory acute myeloid leukemia. Blood 2017;130:722–31.
Tiscornia G, Singer O, Verma IM. Production and purification of lentiviral vectors. Nat Protoc. 2006;1:141–245.
Corces MR, Trevino AE, Hamilton EG, Greenside PG, Sinnott-Armstrong NA, Vesuna S, et al. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nat Methods. 2017;14:959–62.
Bolger AM, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 2010;26:841–2.
Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.
Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. (2010). Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576–89.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2012;29:15–21.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10:1523.
Krämer A, Green J, Pollard JJ, Tugendreich S. Causal analysis approaches in ingenuity pathway analysis. Bioinformatics. 2014;30:523–30.
Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44:W160–5.
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–26.
McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol 2010;28:495–501.
Zeng Y, Shen Z, Gu W, Wu M. Bioinformatics analysis to identify action targets in NCI-N87 gastric cancer cells exposed to quercetin. Pharm Biol. 2018;56:393–8.
Corces MR, Granja JM, Shams S, Louie BH, Seoane JA, Zhou W, et al. The chromatin accessibility landscape of primary human cancers. Science 2018;362:eaav1898.
ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489:57–74.
Acknowledgements
We thank Guanhua Yan for excellent technical assistance, Connor Rogerson for advice and help with RNA-seq data processing and staff in the Bioinformatics, Genomic Technologies and Flow cytometry core facilities. Rebecca Fitzgerald for providing the OAC organoid. We also thank Nicoletta Bobola, Tim Somervaille and Shen-Hsi Yang for critical appraisal of the manuscript. This work was funded by grants from the Wellcome Trust (103857/Z/14/Z and 102171/Z/13/Z).
Author information
Authors and Affiliations
Contributions
SO, KC and IA performed the experiments and data analysis in this study; JB led the metabolic studies; ADS led the experimental parts of the programme. All authors contributed to manuscript preparation and/or critically appraised manuscript drafts. OCCAMS consortium (see Supplementary File 1) provided OAC sample collections and associated RNA-seq and DNA sequencing data.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Ogden, S., Carys, K., Ahmed, I. et al. Regulatory chromatin rewiring promotes metabolic switching during adaptation to oncogenic receptor tyrosine kinase inhibition. Oncogene 41, 4808–4822 (2022). https://doi.org/10.1038/s41388-022-02465-w
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41388-022-02465-w