Abstract
A major hurdle to the application of precision oncology in pancreatic cancer is the lack of molecular stratification approaches and targeted therapy for defined molecular subtypes. In this work, we sought to gain further insight and identify molecular and epigenetic signatures of the Basal-like A pancreatic ductal adenocarcinoma (PDAC) subgroup that can be applied to clinical samples for patient stratification and/or therapy monitoring. We generated and integrated global gene expression and epigenome mapping data from patient-derived xenograft models to identify subtype-specific enhancer regions that were validated in patient-derived samples. In addition, complementary nascent transcription and chromatin topology (HiChIP) analyses revealed a Basal-like A subtype-specific transcribed enhancer program in PDAC characterized by enhancer RNA (eRNA) production that is associated with more frequent chromatin interactions and subtype-specific gene activation. Importantly, we successfully confirmed the validity of eRNA detection as a possible histologic approach for PDAC patient stratification by performing RNA-ISH analyses for subtype-specific eRNAs on pathologic tissue samples. Thus, this study provides proof-of-concept that subtype-specific epigenetic changes relevant for PDAC progression can be detected at a single-cell level in complex, heterogeneous, primary tumor material.
Subtype-specific enhancer activity analysis via detection of eRNAs on a single-cell level in patient material can be used as a potential tool for treatment stratification.
Introduction
Pancreatic ductal adenocarcinoma (PDAC) is the most common form of pancreatic cancer with a complex molecular landscape and highly heterogeneous nature. Despite the decline in overall cancer death rates, the mortality rate of pancreatic cancer remains high due to the lack of proper early detection methods and effective therapeutic options beyond surgery. Consequently, pancreatic cancer is projected to become the second leading cause of cancer-related deaths by 2030 (1). Accordingly, there is an urgent need to elucidate the molecular mechanisms underlying tumor development, therapeutic resistance, and metastasis of this lethal disease.
A number of large-scale studies involving transcriptomic, genetic and epigenetic analyses of PDAC tumor samples (2–6) consistently identified two principally distinct molecular subtypes referred to as “classical-” and “Basal-like”, which display differences in phenotypic behavior, outcome and therapeutic responsiveness. Subsequently, a more detailed transcriptomic analysis including samples from advanced-stage PDAC tumors (7) revealed that the previously identified major subtypes can be further segregated into two subclusters each based on distinct gene expression signatures. Among the identified subtypes, Basal-like A PDAC was shown to associate with particularly aggressive and chemoresistant behavior. While the origin of Basal-like cell identity in PDAC is still poorly understood, it has been closely linked to the expression of the the ΔN isoform of the transcription factor TP63 (ΔNp63) (3, 8, 9).
ΔNp63 is a master regulator that plays an essential role in the development of stratified epithelial tissues such as epidermis, breast and prostate (10). A number of studies have revealed the oncogenic potential of ΔNp63, which was found to be specifically overexpressed in the “Basal” or “squamous” subtypes of diverse epithelial cancers that display a more aggressive clinical course and poor prognosis (11, 12). While ΔNp63 expression was previously linked to Basal-like PDAC (3, 13, 14), its functional relevance in PDAC was unclear until recently. ΔNp63 directs subtype-specific gene expression in PDAC through the dynamic epigenetic regulation of enhancer landscapes (8, 9, 15) and promotes inflammatory changes in the tumor microenvironment (16).
As a lineage-defining TF, ΔNp63 occupies a large number of binding sites at distal enhancers compared with promoters (8, 9). Enhancers are distal regulatory elements which have a prominent role in driving the expression of lineage-specific genes upon binding of TFs and cofactors, which in turn mediate long-range interactions with target gene promoters (17). Growing evidence suggests that tumor-associated molecular aberrations result in a rewiring of the enhancer networks, thereby directing cancer cell identity (18, 19). Notably, in addition to epigenetic marking, some enhancers are transcribed and produce short-lived noncoding RNA molecules referred as enhancer RNAs (eRNA; refs. 20, 21). Although the functional relevance of eRNAs remains to be definitively established, multiple studies have uncovered an activating role of eRNAs in transcriptional regulation (22–26). Spatio-temporal control of gene expression is achieved by the dynamic interaction of distal enhancers with the target promoters. Interestingly, eRNAs have been reported to contribute to the maintenance and stabilization of enhancer-promoter loop structures (18, 27, 28). In fact, eRNAs can directly associate with TFs and stabilize their binding to the cognate enhancers (27–30). Hence, based on the current body of evidence, we hypothesized that eRNAs can serve as a precise readout of cell-type specific enhancer activity. Consistently, enhancer activity profiling across breast cancer cell lines was used to identify the TFs driving subtype-specific biology of breast cancer (31).
In this study, we performed comprehensive multi-omics analyses to characterize the molecular determinants of ΔNp63-dependent transcription in PDAC. For the first time, we report a Basal-like A subtype-specific transcribed enhancer program (B-STEP) which is essential for subtype-specific gene expression in Basal-like A PDAC. Moreover, utilization of a CRISPR-dCas9-KRAB epigenome editing approach confirmed the essential role of several B-STEP enhancers. Importantly, we show that the B-STEP program and in vivo growth of Basal-like patient-derived xenograft (PDX) can be inhibited by pharmacologic targeting of Bromodomain- and Extra-Terminal (BET) domain proteins. Finally, we validated the detectability of B-STEP eRNAs in formalin-fixed paraffin-embedded (FFPE) tissue samples from xenograft and patient samples. Overall, this study uncovers a unique repertoire of transcribed, subtype-specific enhancers in PDAC and displays the potential utility of their detection for clinical stratification and/or accurate monitoring of therapeutic response.
Materials and Methods
PDX
PDX models generated in Göttingen, Germany (GöPDX) were prepared as previously described (32). Informed and written consent was obtained from all patients. Most PDX samples used for mRNA sequencing (mRNA-seq) were F1 generation with exception of GöPDX5 (F4), GöPDX12 (F2), and GöPDX13 (F3). Animal experiments were performed following the protocols approved by the Institutional Animal Care and Use Committee (IACUC; 33.9–42502–04–17/2407). The generation and utilization of PDX models was approved by the ethical review board of the University Medical Center Göttingen (UMG; 70112108).
Mouse xenografts
L3.6pl xenografts were generated during a previous study (33). HT29 xenografts were generated by injecting one million cells under sterile conditions subcutaneously into the right abdominal flank of SCID mice. Tumors were allowed to grow for 21 days. Mouse weight and tumor size were measured daily. Animal experiments were performed following the protocols approved by the IACUC (33.9–42502–04–15/2039).
Cell culture
L3.6pl (RRID:CVCL_0384; ref. 34) were kindly provided by Elisabeth Hessmann (UMG, Germany) and AsPC1 cells (RRID:CVCL_0152) – by Frauke Alves (35) [Max-Planck-Institute for Multidisciplinary Scienes Göttingen, Germany]. L3.6pl cells were maintained and cultured in Eagle Minimum Essential Medium (Gibco, Life Technologies) and AsPC1 cells - in RPMI medium (Gibco, Life Technologies). Cells were regularly tested for Mycoplasma using the MycoAlert Mycoplasma Detection Kit (Lonza, #LT07–318). Cells were maintained in culture for a maximum of 10 to 15 passages on average before thawing new stock. Although cell line authentication was not performed, cell lines were regularly monitored for phenotypic appearance and behavior. Transfection of siRNAs was achieved using Lipofectamine RNAiMAX (Invitrogen) according to the manufacturer's instructions. Control non-targeting siRNA #5 (NT5) (siGENOME, #D-001210–05–20, target sequence: UGGUUUACAUGUCGACUAA) and TP63 siRNA pool (siGENOME; #M-003330–01–0005 SMARTpool consisting of D-003330–05 CAUCAUGUCUGGACUAUUU; D-003330–06 CAAACAAGAUUGAGAUUAG; D-003330–07 GCACACAGACAAAUGAAUU; D-003330–08 CGACAGUCUUGUACAAUUU) were purchased from Dharmacon, ThermoScientific. For JQ1 experiments, 24 hours after seeding, cells were treated with 500 nmol/L JQ1 (dissolved in DMSO) or an equivalent volume of DMSO for 24 hours.
For proliferation assay, 6,000 cells were seeded in each well of a 48-well plate and left to grow overnight. Next day, cells were treated with different concentrations of pan BET inhibitor CC-90010 (Trotabresib, BMS-986378) (MedChemExpress) and transferred to Incucyte S3 Live Cell Analysis Instrument (Sartorius; RRID:SCR_023147) for live cell imaging. Cell growth inhibition was assessed for total duration of 96 hours while treatment was replenished every 48 hours. Confluence was measured by Incucyte 2022B Rev2 software. IC50 was calculated from dose-response curves (three parameters) using GraphPad Prism 9.5.1 (RRID:SCR_002798).
Establishment of stable dCas9-KRAB expressing cell lines
HEK293T (RRID:CVCL_0063) cells were transfected with pMD2.G (Addgene, #12259, RRID:Addgene_12259), psPAX2 (Addgene, #12260, RRID:Addgene_12260), and Lenti-dCas9-KRAB-blast (Addgene, #89567, RRID:Addgene_89567) using Lipofectamine 3000 (Invitrogen) according to the manufacturer's instructions. Viral supernatant was pooled by collecting the media at 24, 48, and 72 hours post-transfection and centrifuged at 500 × g for 10 minutes to remove any cellular material. 500,000 L3.6pl cells were transduced by spin-oculation at 1,000 × g for 2 hours at 30°C using 8 μg/mL polybrene in a dilution of viral supernatant and normal culture media. After 48 hours, the cells were resuspended in appropriate culture medium supplemented with 10 μg/mL Blasticidin S (EMD Millipore Corp). Selection was maintained throughout all experiments.
For epigenome editing, 400,000 L3.6pl-dCas9-KRAB cells per well (in a 6-well plate) were reverse transfected using Lipofectamine RNAiMax (Invitrogen) with custom designed synthetic guide-RNAs (Integrated DNA Technologies) targeting individual super-enhancers (SE) or the SNAI2 transcriptional start site (Supplementary Table S1) following the manufacturer's instructions. RNA was harvested 48 hours post-transfection and gene expression was assessed via qRT-PCR.
IHC staining
IHC staining was performed in the Clinical Pathology Core of the Mayo Clinic (Rochester, MN) with the Leica Bond platform (Leica) using antibodies against p40 (ΔNp63, Biocare, #ACI 3066 A, C; clone: BC28; RRID:AB_2858274), HNF1B (Sigma, #HPA002083, Polyclonal; RRID:AB_1080232), KRT5 (Leica, #NCL-L-CK5; clone: XM26) and KRT7 (Dako, # M7018; clone: OV-TL 12/30) after heat-induced antigen retrieval. Studies on human samples were approved by the Institutional Review Board of the Mayo Clinic (#19–005710).
RNAScope
RNA-ISH staining was performed at the Mayo Clinic Pathology Research Core (Rochester, MN) using the Leica Bond RX stainer (Leica). RNA-ISH staining was performed on-line. Slides were dewaxed and baked prior to being exposed to the retrieval solution, Epitope Retrieval 2 (EDTA; Leica) for 15 minutes at 95°C and the Protease III pretreatment (ACD) for 15 minutes at 40°C. The slides were incubated for 2 hours at 40°C with either a positive control probe (Hs-PPIB probe, ACD), negative control (DapB, ACD) or custom designed probes against the eRNA-producing regions contained within the SEs associated with SNAI2 (Hs-eSNAI2–2, ACD) or TP63 (Hs-eTPRG1, ACD). After hybridization, the slides were incubated in a series of amplification reagents (RNAscope 2.5 LSx Reagent Kit-Brown, ACD) according to the manufacturer's protocol. Slides were incubated for 20 minutes in DAB for visualization and counterstained in hematoxylin for 5 minutes, followed by a bluing agent for 2 minutes. Slides were removed from the stainer and rinsed in tap water, dehydrated in increasing concentrations of ethanol, and cleared 3 times in xylene prior to mounting with coverslips in xylene-based medium.
Chromatin immunoprecipitation
Chromatin immunoprecipitation (ChIP) was performed as previously described (32). PDX tissue samples were homogenized on dry ice using the Sciences BioMascher II Micro Tissue homogenizer. Powdered tissue and cell pellets from L3.6pl cells were fixed in 1% formaldehyde in PBS for 20 minutes and quenched with 125 mmol/L glycine for 5 minutes at room temperature. Sonicated chromatin extract was immunoprecipitated overnight with 1 μg H3K27ac antibody (Diagenode, #C15410196, RRID:AB_2637079) or 1 μg control rabbit IgG (Diagenode, #C15410206, RRID:AB_2722554). Protein A sepharose was used to capture immune complexes. The ChIP efficiency was verified via qRT-PCR and normalized to input DNA.
Length extension chromatin run-on followed by sequencing
A detailed protocol for length extension chromatin run-on followed by sequencing (leChRO-seq) was kindly provided by Charles G. Danko and was performed as previously described (36).
H3K27ac HiChIP
HiChIP was performed as previously described (37). Cells were crosslinked with 1% formaldehyde. For IP, 6 μg of H3K27ac antibody Diagenode, #C15410206, RRID:AB_2722554) were used. Sequencing libraries were prepared using the KAPA HyperPrep (Roche) kit. The libraries were sequenced by the Genome Analysis Core at Mayo Clinic, Minnesota.
Statistical analysis
For gene expression analyses, statistical comparisons were performed using Student t test or one-way ANOVA followed by Dunnett post test in GraphPad Prism (GraphPad Software, RRID:SCR_002798). Individual P values and type of statistical test are indicated in the figure legends. Data is represented by the mean and standard deviation.
Data availability
The sequencing data generated in this study have been deposited into the Gene Expression Omnibus (RRID:SCR_005012) under GSE213397 superseries. The additional data used in this study were obtained from ArrayExpress (http://www.ebi.ac.uk/arrayexpress) under accession codes E-MTAB-5632 (PDX H3K27ac ChIP-seq), E-MTAB-5639 (PDX mRNA-seq; ref. 38); E-MTAB-7034 (L3.6pl TP63 & BRD4 ChIP-seq; ref. 9).
Further methodologic details including tumor growth analysis following JQ1 treatments of PDX-bearing mice, library preparation and next-generation sequencing and data analysis can be found in Supplementary Data.
Results
PDX models faithfully recapitulate PDAC subtype-specific gene expression signatures
A recent study showed that the two major PDAC molecular subtypes “Basal-like” and “classical” PDAC can be further classified into further subgroups based on distinct gene expression patterns (7). To determine whether these molecular subtypes are faithfully reflected in experimental model systems, we established an independent cohort of PDX samples from tumor specimens resected at the University Medical Center Göttingen, Germany (GöPDX). Application of signature gene analysis with subsequent clustering of these specimens identified two PDX samples (GöPDX4 and GöPDX15) clearly classifiable as Basal-like A, expressing higher levels of “Signature 2” genes (7), which define the most aggressive Basal-like A PDAC subtype in patient samples (Fig. 1A, red box). Of the 13 samples analyzed, eight (GöPDX5, 8, 9, 11, 12, 19, 21, 23) clearly expressed “Signature 1 and 6” genes and were hence classified as “classical A/B” (Fig. 1A, blue box). We verified subtype identity by examining the expression of several signature genes in PDX samples classified as Basal-like A (GöPDX4 and GöPDX15) and “classical A/B” (GöPDX5 and GöPDX8; Supplementary Fig. S1A). Importantly, IHC analyses of the Basal, Signature 2 marker KRT5 and the classical, Signature 6 marker HNF1B conducted in primary PDAC tissue from the donor patients matching these PDX tumors supported the respective molecular subtype (Fig. 1B). Hence, these PDX models indeed recapitulate molecular features of the original tumors and thus represent a powerful preclinical tool to understand the underlying molecular biology and targeting thereof.
ΔNp63 defines the Basal-like A specific enhancer program
Subtype-specific gene expression patterns in cancer are controlled by the epigenetic enhancer landscape and the tissue-specific TFs that define them (39). The oncogenic function of ΔNp63 in PDAC has been linked with the activation of a “squamous” gene expression program that defines the aggressive clinical course of Basal-like PDAC (7–9). To further explore the direct role of ΔNp63 in preclinical PDX models, we performed ChIP-seq analysis for ΔNp63 in a Basal-like A PDX (GöPDX4). Consistent with a role in gene activation, ΔNp63-occupied genomic regions in GöPDX4 were marked by the active epigenetic mark H3K27ac (Fig. 1C). As a key driver of cell fate specification, ΔNp63 exerts its function via establishment of squamous-specific enhancer clusters, i.e., SEs, which have been strongly implicated in oncogenic transcriptional pattern (40). In fact, ΔNp63 displayed around 18% promoter-specific binding pattern, whereas around 40% of ΔNp63 occupancy occurred at distal intergenic regions (Supplementary Fig. S1B). Consistently, of 551 SEs identified in GöPDX4 (Fig. 1D), ca. 60% were occupied by ΔNp63 (Fig. 1E; Supplementary Table S2). Moreover, among the genes associated with ΔNp63-bound SEs (Supplementary Table S3) 23 belonged to Basal-like A “Signature 2” as indicated in (Fig. 1E). To verify the “Basal-like” nature of identified SEs we used data from an independent published cohort of PDX PDAC samples (38). Consistent with previous classification, three PDX samples (1.033; AO.IP; foie8) that were characterized as squamous based on ΔNp63-driven epigenetic and gene expression patterns (9), displayed a particularly high expression of genes from “Signature 2” genes (Supplementary Fig. S1C, red box). The remaining three PDX tumors classified as Basal in the initial publication (37) expressed genes of both “Basal” and “classical” signatures to differing degrees, probably more closely representing the “hybrid” subtype of PDAC. Nine PDX specimens consistently expressed “classical A” and “classical B” genes, with no clear separation into either of the two classical subclusters. Thus, we classified these as “classical A/B” (Supplementary Fig. S1C, blue box). Consistent with a key role in cellular identity specification, the ΔNp63-enriched SEs from GöPDX4 displayed a subtype-specific activation pattern that was not marked by H3K27ac in the classical PDX samples (Fig.1F; Supplementary Fig. S1D). In summary, the identification of Basal-like A-specific enhancers supports a role for ΔNp63 in shaping the enhancer architecture of the most aggressive PDAC subtype.
ΔNp63 activates a transcribed enhancer program to maintain oncogenic transcription
To further delineate the underlying mechanisms of ΔNp63-dependent enhancer activation in shaping the Basal-like A specific disease phenotype, we sought to use a cell culture model where ΔNp63 expression could be modulated. Consistent with our previous identification of ΔNp63 as a central component defining L3.6pl cellular identity (9), an unbiased principal component analysis of all H3K27ac-enriched regions revealed that L3.6pl cells and GöPDX4 were both similar to the Basal-like A PDX samples from Lomberk, and colleagues (Fig. 2A). Interestingly, upon loss of ΔNp63 in L3.6pl cells (Fig. 2B; Supplementary Fig. S2A) whereas some of “Signature 2” genes were downregulated, the expression of the three other signatures was mainly enriched (Fig. 2C; Supplementary Fig. S2B), thereby confirming the central role of ΔNp63 in driving the squamous gene expression program and suppressing other subtype identities.
After validating the Basal-like A identity of L3.6pl cells and their high degree of similarity to patient-derived tumor samples, we sought to identify the ΔNp63-depedent enhancers that control the identity of these cells. To specifically elucidate the function of ΔNp63-bound distal regulatory elements, we focused on distal ΔNp63 peaks not contained within the transcribed region of any active genes (Fig. 2D) and performed differential H3K27ac occupancy analysis on these upon ΔNp63 depletion. In total, 718 ΔNp63-bound H3K27ac regions displayed significant change, with the 83% of these sites losing H3K27ac occupancy upon ΔNp63 depletion (Fig. 2D and E), supporting the key function of ΔNp63 in establishing the epigenetic status of the distal regulatory sites. Moreover, these ΔNp63-dependent distal sites (Supplementary Table S4) displayed typical active enhancer features characterized by BRD4 and H3K4me1 occupancy (Fig. 2F).
A number of recent studies have revealed that enhancer activity strongly correlates with noncoding eRNA production (41). To identify to what degree the ΔNp63-dependent enhancer repertoire was associated with active transcription, we performed leChRO-seq (36) that allows mapping of RNAPII-driven nascent transcription. Whereas nascent transcription remained unaffected on ΔNp63-independent H3K27ac bound regions (Supplementary Fig. S2C and S2D), bidirectional RNA transcription was substantially reduced on ΔNp63-dependent enhancer regions following loss of ΔNp63 (Fig. 2G), thereby further supporting a central regulatory role of ΔNp63 in controlling a Basal-like A transcribed enhancer program.
Transcribed enhancers display frequent higher order chromatin interactions with target genes
Enhancers regulate gene expression via long-range chromosomal interactions with target gene promoters, and eRNAs have been shown to contribute to stabilization of these contacts (27, 28). To identify genes directly associated with ΔNp63-dependent enhancers we performed HiChIP in L3.6pl cells, using H3K27ac to identify genomic interactions with active enhancer regions. In total, 626 genes were found to engage with ΔNp63-dependent enhancers (Fig. 3A; Supplementary Table S5) including 76 ΔNp63-dependent genes. In a recent work we identified a subgroup of highly interacting enhancers associated with chemoresistance in PDAC that displayed particularly high rates of eRNA production (42). In addition, high transcriptional activity was found to be associated with subset of enhancers forming dense interactome with genes defining the cell identity (43). Thus, we hypothesized that ΔNp63-controlled eRNA transcription may be associated with enhancer-gene interaction frequency as well. We classified ΔNp63-dependent enhancers into 5 groups based on the number of interactions originating from them (Fig. 3B) and observed that eRNA production positively correlated with interaction frequency (Fig. 3C). This finding is consistent with a potential role of eRNAs and/or the active RNA Polymerase II in mediating and/or stabilizing the enhancer-promoter looping. Consistently, the 338 GöPDX4 ΔNp63-bound SEs identified in (Fig. 1D) displayed similar ΔNp63-dependent active enhancer features in L3.6pl cells based on H3K27ac and nascent transcription profile changes (Fig. 3D-F; Supplementary Fig. S3A and S3B).
One example of a “Signature 2” gene associated with a ΔNp63-dependent eRNA-producing SE is the SNAI2 gene (Fig. 3E). SNAI2 was not previously recognized as a direct target of ΔNp63-dependent SE due to the limitations of the 2D genome-association rules typically used (nearest genes), but was identified via our HiChIP analyses and thereby revealed a mechanistic link between ΔNp63 expression and SNAI2-driven epithelial-to-mesenchymal (EMT) transition (44–46). To assess the specific role that ΔNp63-dependent, eRNA-producing enhancers play in controlling target gene expression, we performed CRISPR/dCas9-KRAB–mediated epigenome editing to inactivate the SE associated with the FAT2, SNAI2 and TP63 genes (Supplementary Fig. S3C; Fig. 3G). In case of the latter, TPRG1 was shown to harbor a SE that directly regulates ΔNp63 expression (15). Notably, dCas9-KRAB-mediated inactivation of SE activity resulted in the significant and specific downregulation of all 3 target genes. Furthermore, consistent with the master regulatory function of ΔNp63, silencing of the TPRG1 SE upstream of the TP63 gene not only attenuated expression of ΔNp63, but also resulted in downregulation of the ΔNp63-dependent, SE-driven genes SNAI2 and FAT2. In contrast, inactivation of the SE associated with SNAI2 or FAT2 specifically downregulated the expression of the associated target gene, but not the other SE-associated genes investigated (Fig. 3G). Altogether, these results indicate that the ΔNp63-established B-STEP is a defining feature of the most aggressive subtype of pancreatic cancer.
Subtype-specific eRNA detection
The highly specific expression pattern of eRNAs suggests that they might serve as specific and sensitive markers of enhancer activity and therefore aid in the selection of therapy and monitoring of therapeutic efficacy. To test the potential clinical utility of detecting B-STEP eRNAs, we evaluated the detectability of TPRG1 and SNAI2 eRNAs in FFPE xenograft samples from the Basal-like PDAC cell line L3.6pl and a negative control colorectal cancer cell line (HT29) by using the RNAScope RNA-ISH approach (Fig. 4A and B). Remarkably, one to two individual foci could be observed for TPRG1 and SNAI2 in many nuclei of the L3.6pl xenograft sample, but not in the colorectal cancer xenograft. To confirm the feasibility of this approach on clinical material, we used histologic tumor material obtained from a Mayo Clinic patient following resection, which we confirmed to be of the Basal subtype based on the expression of both ΔNp63 (“p40”) and its downstream target KRT5 (Fig. 4C). Notably, positively stained nuclear foci could be observed for both TPRG1 and SNAI2 eRNAs in individual cells (Fig. 4D).
Our previous work revealed an essential role of the Bromodomain and Extraterminal (BET) histone reader protein BRD4 in controlling eRNA expression in breast cancer (47). Moreover, BET family member bromodomain testis-specific protein (BRDT) was found to co-localize and interact with ΔNp63 to drive a unique transcriptional program and modulate cell phenotype in esophageal squamous cell carcinoma (37). Given the central role of BET proteins as a defining feature of SE activity (48, 49) in cancer, we hypothesized that small molecule inhibition of BET proteins may inhibit B-STEP-dependent transcription of the Basal-like A-specific gene expression pattern. Consistently, the BET inhibitor JQ1 downregulated ΔNp63-dependent (Fig. 4E) and Basal-like A-defining “Signature 2” genes (Fig. 4F). Moreover, similar to the effects of ΔNp63 loss (Supplementary Fig. S3D and S3E), JQ1 treatment abrogated expression of B-STEP eRNAs (Fig. 4G). Importantly, based on cell growth inhibition in L3.6pl cells ΔNp63-program conferred higher sensitivity to trotabresib, a reversible and orally active pan-BET inhibitor CC-90010 (50) currently under phase-I clinical trials (NCT03220347), when compared with cells of classical origin AsPC1 (ref. 8; Fig. 4H). Finally, in line with the anti-cancer effect of BETi, JQ1 treatment of mice implanted with a Basal-like A PDX (Supplementary Fig. S3F) substantially impaired tumor growth (Fig. 4I).
In summary, we conclude that B-STEP eRNAs can be detected on a single-cell level in histologic and primary patient material, are highly sensitive to BETi treatment and can therefore likely serve as predictive markers for therapy stratification.
Discussion
A number of recent studies have examined the transcriptomic and genomic diversity within PDAC and uncovered distinct subtypes of PDAC with different clinical courses and therapeutic responsiveness (2–5, 51, 52). In this work, we have uncovered important epigenetic principles controlling the most aggressive Basal-like A subtype of PDAC. Growing evidence indicates that epigenetic reprogramming of the enhancer landscape is a key regulatory event controlling tumor cell identity. However, while technological advances have enabled the detailed interrogation of the mechanistic and epigenetic underpinnings of tumor subtype-specific gene expression in cell culture, the application of these findings to the clinical realm has been halted by the impracticality and difficulty of applying complex high throughput procedures to clinical material. In this study, we used a multi-omics approach to characterize the epigenetic and transcriptional landscape of clinically relevant model systems and identified a novel B-STEP we have termed B-STEP. Our findings demonstrate that eRNAs can function as a robust proxy of subtype-specific enhancer activity in PDAC.
To date, the application of genomic methodologies for the identification of subtype-specific enhancers in patients has been fraught with difficulties including technical limitations due to low quantities of primary biopsy material, extensive desmoplasia and high tumor heterogeneity. For the first time we not only identified tumor subtype-specific eRNAs, but also demonstrated the ability to detect these in histologic tumor material. We propose that detection of eRNAs in histologic material may serve as a novel diagnostic approach for patient stratification on a single-cell level. eRNA expression reflects mRNA expression kinetics more closely than H3K27ac occupancy at enhancers (53) and alterations in enhancer transcription precede the successive wave of changes in transcription (54), suggesting that analysis of eRNAs underlying complex cell identity programs in patient tumors and/or circulating tumor cells can be alternatively used as a tool for patient stratification and the assessment of therapeutic responses.
A major hurdle to identifying mechanism-based, clinically meaningful biomarkers is the availability of experimental model systems that faithfully recapitulate patient tumor characteristics. Established cancer cell lines often significantly diverge from primary tumor samples due to the high selective pressure of 2D cell culture conditions and lack microenviromental signals. Our results show that PDX models largely reflect the diversity in gene expression programs observed in clinical PDAC samples. Moreover, consistent with our previous work (9), our results demonstrate that the highly metastatic L3.6pl cell line closely mirrors the transcriptional and epigenetic enhancer program present in Basal-like A PDX models.
Functionally, eRNAs likely promote target gene activation and/or promoter/enhancer looping by interacting with transcriptional activators and co-activators (23, 27–29). Given their short length and unstable nature, discovery of eRNAs has proven challenging, especially in patient-derived clinical samples with substantial RNA degradation. However, this limitation was recently overcome by the development of leChRO-seq (36). Application of this technique in L3.6pl cells and a Basal-like A PDX enabled us to identify the B-STEP. Moreover, consistent with the recent data by Hamdan and colleagues (42) and Briend and colleagues (43), we similarly observed a positive correlation between contact frequency and the transcriptional activity of B-STEP enhancers, which also supports a role of eRNAs in promoting enhancer-promoter interactions.
Notably, while epigenetic activation of enhancers has frequently been correlated with cell state-specific gene expression, the direct interaction of enhancers with target genes and their necessity for the expression of specific genes has rarely been investigated. Our HiChIP analyses enabled us to identify direct targets of B-STEP enhancers. Moreover, we confirmed the essential role of these enhancers via epigenome editing. Intriguingly, one B-STEP ΔNp63-dependent target gene is the EMT TF SNAI2, which is one of the “Signature 2” genes (7). Our results reveal that SNAI2 expression in Basal-like A PDAC is under the direct control of an upstream ΔNp63-dependent SE. The aggressive and highly invasive nature of the Basal-like A PDAC subtype has been partially explained by its higher enrichment of EMT gene signatures (7, 55) and recently SNAI2 was shown to contribute to a more aggressive behavior of PDAC by inducing EMT in response to nutrient deprivation (46). Although the role of the ΔNp63-SNAI2 axis in other cancers has been reported before (44, 45), we describe the first mechanistic link explaining ΔNp63 and SNAI2 co-expression in Basal-like A PDAC. Moreover, we confirm that ΔNp63 plays a key master regulatory function in establishing the transcribed enhancer landscape of Basal-like A PDAC subtype. Furthermore, perturbation of the TPRG1 SE that controls ΔNp63 expression (15) is sufficient to block downstream ΔNp63-driven enhancer activity and target gene expression.
Another important finding from our work is the discovery that the B-STEP program is particularly amenable to modulation by small molecule BET inhibitor treatment both in vitro and in vivo. This finding is particularly relevant for the clinical application of BETi for PDAC treatment. A major limitation to the clinical utility of BETi is the occurrence of treatment-emergent adverse events (TEAE) such as thrombocytopenia. To date, BETi have only been marginally effective as a monotherapy in clinical trials for solid tumors, probably due to the limited tumor cell cytotoxicity that they invoke. On the basis of our findings we propose a new subtype-specific mechanism-based strategy to BETi application and dosage in PDAC. Specifically, we propose that B-STEP eRNA transcription may be a better indicator of on-target effects and can serve as a highly sensitive and accurate readout for clinical dosing and monitoring of BETi treatment in Basal-like A PDAC. Importantly, we believe this approach may enable lower, more accurate BETi dosing that will help to minimize TEAE. On the basis of the findings that Basal-like A PDAC is generally refractory to standard of care chemotherapy, we anticipate that transient inhibition of the B-STEP program and hence of the Basal-like A phenotype by BETi treatment can serve as an ideal inductive therapy to sensitize Basal-like A PDAC tumors to chemotherapy treatment. According to our model, verification of B-STEP eRNA expression may serve as an ideal biomarker for early phase clinical trial dosing studies and later for patient stratification (i.e., inclusion criteria for BETi clinical trials) as well as for monitoring of on-target therapeutic responsiveness. We envision that therapeutic strategies targeting cell identity programs will significantly benefit from the utilization of improved mechanism-based biomarkers and enable a better, more accurate identification and evaluation of therapy-induced effects.
Limitations of the study
To validate the potential utility of subtype-specific eRNA expression as a marker for patient/treatment stratification, more patient and PDX samples need to be included in the study. However, given the low level of eRNA expression and their high instability it is technically challenging to obtain suitable samples for analysis. Very often eRNA integrity is highly compromised in archived FFPE patient material with varying fixation conditions. Moreover, given the fact that most Basal A tumors are diagnosed at an advanced stage (i.e., using biopsies) fresh material is not readily available and not in sufficient quantity. While we agree that large numbers would strengthen our findings, we do believe that the data presented in the study demonstrate that eRNAs can be detected under suitable conditions, both in xenograft and primary patient material.
Authors' Disclosures
P. Hermann reports grants from Deutsche Krebshilfe (Max Eder Program), Hector Foundation; and grants from Deutsche Forschungsgemeinschaft (SFB1279) during the conduct of the study. J.T. Siveke reports grants and personal fees from Bristol-Myers Squibb/Celgene, Roche/Genentech, AstraZeneca,; grants from Eisbach Bio, Abalos Therapeutics; personal fees from Immunocore, Bayer, MSD Sharp Dohme, Novartis, Servier; and other support from Pharma15 outside the submitted work. No disclosures were reported by the other authors.
Authors' Contributions
X. Wang: Data curation, software, formal analysis, investigation, visualization, writing–original draft, writing–review and editing. A.P. Kutschat: Formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. J. Aggrey-Fynn: Investigation, visualization, writing–review and editing. F.H. Hamdan: Investigation, writing–review and editing. R.P. Graham: Methodology, writing–review and editing. A.Q. Wixom: Investigation, writing–review and editing. Y. Souto: Validation, investigation, writing–review and editing. S. Ladigan-Badura: Resources, methodology, writing–review and editing. J.A. Yonkus: Methodology, writing–review and editing. A.M. Abdelrahman: Resources, methodology, writing–review and editing. R. Alva-Ruiz: Resources, methodology, writing–review and editing. J. Gaedcke: Resources, methodology, writing–review and editing. P. Ströbel: Resources, methodology, writing–review and editing. R.L. Kosinsky: Resources, methodology, writing–review and editing. F. Wegwitz: Resources, methodology, writing–review and editing. P. Hermann: Funding acquisition, writing–review and editing. M.J. Truty: Resources, methodology, writing–review and editing. J.T. Siveke: Funding acquisition, writing–review and editing. S.A. Hahn: Resources, funding acquisition, methodology, writing–review and editing. E. Hessmann: Resources, funding acquisition, methodology, writing–review and editing. S.A. Johnsen: Conceptualization, formal analysis, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. Z. Najafova: Conceptualization, supervision, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.
Acknowledgments
This project was funded by the Robert Bosch Stiftung (to S.A. Johnsen, Z. Najafova, R.L. Kosinsky), the Deutsche Krebshilfe (70112505; PIPAC Consortium; to S.A. Johnsen, E. Hessmann, S.A. Hahn, J.T. Siveke, P. Hermann, F. Wegwitz, J. Gaedcke, S. Ladigan-Badura, and Z. Najafova), a grant from NCI (CA 102701; to S.A. Johnsen), a grant from National Institute of Diabetes and Digestive and Kidney Diseases (T32 DK07198; to A.Q. Wixom), grants from the Deutsche Forschungsgemeinschaft (DFG, JO 815/3–2 and KFO 5002, to S.A. Johnsen and E. Hessmann, P. Strobel, J. Gaedcke, respectively). Z. Najafova was partially supported by the Dorothea Schlözer Program (University of Göttingen). This work was supported by the Optical Microscopy Core of the Mayo Clinic Center for Cell Signaling in Gastroenterology (P30DK084567). Work in the lab of J.T. Siveke is supported by the German Cancer Consortium (DKTK), Deutsche Forschungsgemeinschaft (DFG; 405344257/SI1549/3–2 and SI1549/4–1); the German Federal Ministry of Education and Research (BMBF; 01KD2206A/SATURN3). The authors thank the Genome Analysis Core at the Mayo Clinic for performing the sequencing and the Pathology Research Core at Mayo Clinic for performing IHC analysis and RNAScope.
The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.
Note: Supplementary data for this article are available at Molecular Cancer Research Online (http://mcr.aacrjournals.org/).