Chromatin accessibility changes at intergenic regions are associated with ovarian cancer drug resistance

Abstract

Background

Resistance to DNA damaging chemotherapies leads to cancer treatment failure and poor patient prognosis. We investigated how genomic distribution of accessible chromatin sites is altered during acquisition of cisplatin resistance using matched ovarian cell lines from high grade serous ovarian cancer (HGSOC) patients before and after becoming clinically resistant to platinum-based chemotherapy.

Results

Resistant lines show altered chromatin accessibility at intergenic regions, but less so at gene promoters. Clusters of cis-regulatory elements at these intergenic regions show chromatin changes that are associated with altered expression of linked genes, with enrichment for genes involved in the Fanconi anemia/BRCA DNA damage response pathway. Further, genome-wide distribution of platinum adducts associates with the chromatin changes observed and distinguish sensitive from resistant lines. In the resistant line, we observe fewer adducts around gene promoters and more adducts at intergenic regions.

Conclusions

Chromatin changes at intergenic regulators of gene expression are associated with in vivo derived drug resistance and Pt-adduct distribution in patient-derived HGSOC drug resistance models.

Background

Platinum-based chemotherapeutics, such as cisplatin and carboplatin, are clinically important first line therapies in the treatment of a wide variety of solid cancers [1, 2]. These drugs exert their DNA damaging, cytotoxic effect by the formation of platinum-DNA adducts, inter- and intra-strand cross-links, which induce cell death through a number of pathways if the adduct is not repaired [3]. While many patients initially respond to platinum-based chemotherapy, they will eventually relapse with disease that fails to respond to treatment leading to poor survival [3, 4]. Understanding the changes that occur in platinum-resistant tumors is essential for rational approaches to circumvent resistance and developing molecularly targeted agents for use in recurrent cancers.

Epigenetic mechanisms play a key role in the development of platinum-resistance and drug tolerance [5]. Cells surviving cisplatin exposure, such as transient drug tolerant persisters, or drug-resistant cells surviving cisplatin selection have genome wide epigenetic alterations [5,6,7,8]. Drug-tolerant persisters exhibit a repressed chromatin state and can serve as founders for further genetic and epigenetic change leading to resistance [7, 9,10,11]. The promoters of genes susceptible to hypermethylation in ovarian tumors during the emergence of resistance to platinum resistance are marked by H3K27 and H3K4 bivalent methylation domains present in tumor cells pre-chemotherapy, further emphasizing the relationship between chromatin states and epigenetic adaptation to platinum treatment [12]. Epigenome profiling of ovarian cell line models selected in vitro for platinum resistance has suggested a role for distal super-enhancers (SEs) and their gene targets in maintain the transcriptional program of the platinum-resistant state [13, 14]. These studies highlight the multi-factorial nature of the transcriptional states driving platinum resistance, while identifying SEs that play a critical role in transcriptional regulation of platinum resistance during in vitro selection. However, it is still poorly understood whether such changes at SEs occur following in vivo selection following platinum-based chemotherapy of patients’ tumors as they develop clinical resistance.

Using matched chemosensitive and chemoresistant ovarian cancer cell lines isolated from high-grade serous ovarian cancer patients before and following treatment, we examined the relationship between chromatin accessibility, platinum-DNA adduct distribution, and chemotherapy resistance. To examine further the potential functional significance of any changes observed, we have correlated chromatin changes, identified using ATAC-seq, between the lines at gene promoters, CpG islands, clusters of cis-regulatory elements (COREs) and other genomic regions with changes in gene expression.

The platinum atom of cisplatin can form covalent bonds to the N7 positions of purine bases and the density of GG dinucleotides is considered the main factor influencing distribution of platinum adducts [15]. However, chromatin states could also affect the pattern of cisplatin cross-linking [16, 17]. Indeed, epigenetic therapies, such as histone deacetylase (HDAC) inhibitors, can increase the accessibility of chromatin to cisplatin through global chromatin decompaction, leading to enhanced cell death [18]. Given the known roles for chromatin in mediating the emergence of drug resistance and since chromatin changes might influence the formation and effect of platinum adducts, we aimed to examine platinum-adduct distribution in the genome of ovarian cancer cells and relate this to chromatin conformation and gene expression.

We developed Pt-exo-seq to map platinum adducts, genome-wide by adapting methods for exonuclease mapping of transcription factor binding sites [19,20,21]. Pt-exo-seq is primarily based on the inhibition of exonuclease activity by the bulky platinum adduct and the subsequent removal of the adduct by cyanide treatment, followed by next generation sequencing. A schematic of the Pt-exo-seq approach is shown in Additional file 1: Fig. S1. This complements existing methods for examining the distribution of platinum adducts genome wide such as Damage-seq [22, 23]. In Damage-seq, cisplatin-damaged DNA fragments are immunoprecipitated using an antibody specific to platinum-bound DNA prior to sequencing using a library preparation protocol dependent on the stalling of DNA polymerase by DNA adducts. Fragments without adducts which are non-specifically bound by the immunoprecipitation step are removed by subtractive hybridization prior to sequencing. This method provides accurate and strand-specific mapping of DNA damage, but is dependent on immunoprecipitation and therefore heavily influenced by the sensitivity and specificity of the antibody used and the conditions for subtractive hybridization. Pt-exo-seq avoids these potential complications while still providing accurate and strand-specific mapping at high resolution.

Results

Chromatin conformation in matched resistant lines following platinum-based treatment

We studied chromatin conformation by ATAC-seq in three isogenic pairs of sensitive and resistant human ovarian cell lines. PEO1 and PEO4 were both isolated from the same patient diagnosed with high-grade serous ovarian cancer (HGSOC) following platinum-based chemotherapy, but before and after the clinical development of resistance [25]. PEA1 and PEA2 were isolated from ascites of a patient with HGSOC prior to platinum chemotherapy and after relapse with resistant disease [25]. A2780cp70 (abbreviated to CP70) is an in vitro derived platinum-resistant derivative of the ovarian cell line A2780 [24]. The characteristics and sensitivities of the cell line to cisplatin are summarized in Additional file 1: Table S1.

Chromatin accessibility profiles were compared between the sensitive and resistant lines in each pair, based on ATAC-seq data produced in triplicate. ATAC-seq reads were assigned to 1 Kb windows tiling across the reference genome (hg19). FPKM values (Fragments Per Kilobase of DNA, per Million mapped reads of sequencing library) were calculated for each window, in each replicate. Individual windows were called differentially accessible based on a log2 fold change in coverage in the sensitive-to-resistant lines > ±2 and moderated t test FDR < 0.001. Substantial alterations to the landscape of chromatin accessibility were found in all three resistant lines compared to their sensitive counterparts, affecting 297 windows (0.009%) of the genome in the A2780 pair (Fig. 1a), 8950 windows (0.27%) of the genome in the PEA pair (Fig. 1b) and 7298 windows (0.22%) in the PEO pair (Fig. 1c) out of 3,485,781 windows examined in the genome. A matrix of pair-wise correlations of chromatin accessibility across consistently differentially accessible windows (adjusting for cell line) showed the resistant PEA2 and PEO4 lines clustering together, separately from their respective sensitive PEA1 and PEO1 counterparts (Fig. 1d).

Fig. 1

Correlation of ATAC-seq chromatin accessibility profiles across ovarian cancer cell lines. Heatmaps of ATAC-seq coverage in regions showing differential accessibility between the cisplatin-sensitive and cisplatin-resistant line in each pair. Each row is a 1-Kb window showing log2FC > ± 2 and FDR < 0.001. Dendrograms show similarity between samples as defined by unsupervised hierarchical clustering a A2780/CP70, b PEA1/PEA2, c PEO1/PEO4. d Heatmap of Spearman’s correlation coefficients for each experimental replicate for each cell line analyzed, based on all windows passing filter

Thus, the resistant HGSOC tumor cells isolated at time of recurrence with non-responsive disease show similarity, suggesting common chromatin changes occurring during the acquisition of resistance in vivo. The HGSOC lines were in turn different in their chromatin accessibility to the non-HGSOC, in vitro resistance derived, A2780-CP70 pair, potentially reflecting their different origin and in vitro method of drug selection.

Platinum resistance associated with differential chromatin accessibility in intergenic regions

Each of the windows showing differential accessibility between sensitive and resistant lines by ATAC-seq was annotated by genomic class using HOMER [35]. Those that were associated with genes were classed as either CpG island, promoter-TSS (− 1 Kb to 100 bp 3′ of Transcription Start Site, TSS), exon, intron, 5′ untranslated region (UTR), 3′ UTR, or TTS (− 100 bp to 1 Kb 3′ of TTS). The number of windows falling within each genomic class are shown in Additional file 1: Table S2. Those falling outside these classes were classed as intergenic. We then calculated odds ratios (OR) of enrichment for each class of genomic element in windows showing either increased or decreased accessibility in the resistant line compared with all windows (shown as log(OR) in Fig. 2).

Fig. 2

Enrichment for genomic elements in regions showing differential accessibility in cisplatin-resistant ovarian cancer cell lines. Log(OR) from Fisher’s exact test for enrichment for each class of genomic element in sets of windows showing increasing or decreasing chromatin, between the sensitive and resistant lines in each pair. a A2780/CP70, b PEA1/PEA2, c PEO1/PEO4. OR > 1 or < 1 corresponds to enrichment for a class. Error bars represent 95% confidence interval

In all three sensitive/resistant pairs, there was under-representation of promoter-TSS, CpG island and exonic windows in the sets of differentially accessible windows, suggesting fewer chromatin changes at these genomic regions during the acquisition of resistance (Fig. 2). In all cases, there was enrichment for intronic windows in the sets of differentially accessible windows, demonstrating chromatin changes at intronic regions as strongly associated with resistance. Intergenic elements were also highly over-represented in differentially accessible windows in all three pairs.

Intersection of expression and promoter chromatin accessibility

Analysis of differential gene expression between sensitive and resistant lines was performed by RNA-seq to detect alterations to gene expression in the resistant lines. Genes were defined as differentially expressed based on a log2FC in normalized read count of > ± 2 and moderated t test FDR < 0.05. The summary of differentially expressed (DE) genes and their direction of change for each pair of cell lines is shown in Additional file 1: Table S3. Although many genes show changes in expression, few show universal differential expression in the same direction in all three pairs. By overlapping lists of genes showing significantly altered expression in all three resistant lines (DE = log2FC > |2|, FDR < 0.05), we identified only four genes which were consistently upregulated, PARP9, SPHK1, DDX60L and BCAM, while only 5 genes were consistently downregulated CCDC80, TLE4, THBS1, MAP1A and VEPH1. Using alternative methods of analysis (a sensitive–resistant linear model incorporating all three pairs, and a rank product analysis) only THBS1 was consistently defined as downregulated using all methods. We confirmed the differential expression of THBS1 and PARP9 using RT-QPCR (p < 0.05, t-test, Additional file 1: Fig. S2A). We further explored the similarities in transcriptional dysregulation in the resistant line in the three pairs using gene set enrichment analysis considering the Hallmark gene sets. While this showed a trend towards similar pathways being altered in all three lines, only genes involved in regulation of the interferon gamma response and regulation of apoptosis were significantly dysregulated in all 3 lines (Additional file 1: Fig. S2B). These results suggest there is minimal overlap between alterations of individual gene expression across these three models, although we cannot exclude more subtle changes in gene expression or gene networks being in common.

In order to relate gene expression with chromatin accessibility, ATAC-seq coverage was measured in the 2 Kb region spanning the transcription start site (TSS) of the genes defined as differentially expressed using rank product analysis in each cell line pair (Additional file 1: Fig. S3). Genes which showed reduced expression in the resistant lines had slightly reduced accessibility compared to the resistant lines, although this did not reach statistical significance (t test p > 0.05) (Fig. 3a, b). Genes which were upregulated in the resistant lines showed the inverse, whereby there was a stronger ATAC-seq signal around the TSS of these genes in the resistant lines in which they were upregulated, compared to the sensitive lines (t test p < 0.01). ATAC-seq data at example genes upregulated or downregulated are shown in Additional file 1: Fig. S3.

Fig. 3

ATAC-seq coverage around TSS of differentially expressed genes. a ATAC-seq coverage around TSS of differentially expressed genes for genes with decreased expression and those with increased expression defined by rank product analysis. b Mean ATAC-seq coverage in 2 Kb around TSS of differentially expressed genes (t-test, NS, not significant; **p < 0.01).

Thus, these trends are consistent with an accessible promoter being associated with a permissive state for gene transcription and an inaccessible promoter with transcriptional silencing. However, relatively few changes in chromatin accessibility were observed in the resistant lines at promoter regions (Fig. 2) and those that did occur only weakly correlated with changes in gene expression (Fig. 3).

Chromatin accessibility change at clusters of cis-regulatory elements associate with gene expression in resistant lines

The differences detected in chromatin conformation between sensitive and resistant lines frequently occurred at intergenic regions. Given the minimal association between chromatin accessibility changes at gene promoters and altered gene expression, we asked how changes to chromatin conformation at intergenic regions associated with transcriptional changes. The algorithm CREAM (Clustering of genomic REgions Analysis Method) was used to define clusters of cis-regulatory elements (COREs) in each of the lines, based on peaks called by MACS2 from the merged alignment file produced from the 3 ATAC-seq replicates from each cell line (https://www.biorxiv.org/content/10.1101/222562v1). The results of this are summarized in Additional file 1: Table S4, along with the numbers of COREs unique to each line. We identified substantial changes to CORE landscapes in all three resistant lines, compared to the sensitive. In both the A2780 and PEO pairs there was a loss of COREs (1030 to 358 and 715 to 537, respectively), while a gain in COREs was found in PEA2 compared to PEA1 (744–811), reflecting the trends in the direction of altered accessibility in each pair. Illustrative examples of COREs with differential accessibility (representative ATAC-seq mapping tracks) for all three cell lines are shown in Additional file 1: Fig. S4.

Having defined COREs that were gained or lost in the sensitive/resistant pairs we examined if this change in chromatin conformation at CORE regions was associated with changes in expression of linked genes. For each CORE gained or lost, genes within 100 Kb were extracted from the Ensembl database and those that were differentially expressed between the sensitive and resistant lines identified from the RNA-seq data. The distribution of t-statistics was then plotted, for the comparison of expression values between RNA-seq replicates from the sensitive and resistant line in each pair. In all three pairs, loss of COREs in the resistant line was associated with a significant shift towards reduced expression of genes proximal to the enhancers. When chromatin accessibility changes led to the gain of a CORE in the resistant line, there was a shift toward increased expression of these genes in the resistant line compared to the sensitive (Fig. 4a).

Fig. 4

Relationship between chromatin changes at COREs and gene expression. a Distribution of t statistics of difference in expression for differentially expressed genes near COREs lost or gained in the resistant line in each pair. Differentially expressed genes show log2FC > ± 2 and FDR < 0.05. Red lines show distribution of t statistics for differentially expressed genes near COREs gained in resistant lines, blue lines show t statistics for differentially expressed genes near COREs lost in resistant lines. b Odds ratios of differentially expressed genes being found near COREs lost or gained in resistant lines compared to their sensitive counterpart. Bars show odds ratios and error bars confidence intervals from Fisher’s exact test. *p < 0.05; ***p < 0.001.

We further tested the association between loss or gain of COREs by testing for enrichment of differentially expressed genes in the 100 Kb surrounding lost/gained COREs. Figure 4b shows the odds ratios of finding differentially expressed genes near COREs which were gained or lost in the resistant line compared to the rest of the genome. In the PEA pair, there was enrichment for differentially expressed genes near COREs lost (p < 0.05, OR 1.3, CI 1.04–1.6, Fisher’s exact test) and a strong enrichment for differentially expressed genes near COREs gained (p < 0.001, OR 1.62, CI 1.32–1.97, Fisher’s exact test). In contrast, there was significant enrichment for differentially expressed genes near COREs lost in PEO4 (p < 0.001, OR 1.98, CI 1.48–2.62, Fisher’s exact test) and to those gained (p < 0.001, OR 1.88, CI 1.31–2.65. Fisher’s exact test). While there was no significant enrichment for differentially expressed genes near COREs which were gained in CP70 from A2780, similar to PEO4, there was a 2.9-fold enrichment for differentially expressed genes near COREs lost in CP70 from A2780 (p < 0.001, OR 2.87, CI 2.28–3.58, Fisher’s exact test).

We further explored the chromatin accessibility data in the PEO pair as a representative case to confirm that the lost/gained COREs identified in this analysis could be associated with known enhancers. We obtained published data on chromatin states in the PEO1 cell line [26] calculated using ChromHMM [27], and compared these with our data on COREs lost/gained between PEO1 and PEO4. COREs lost in PEO4 were strongly enriched for active enhancers and promoters in PEO1, while COREs gained in PEO4 were less strongly enriched for these chromatin states, as these reflect the chromatin organization specific to the PEO1 line (Additional file 1: Fig. S5). Similarly, heterochromatin and latent/inactive chromatin was more strongly underrepresented in COREs lost in PEO4, compared to those gained in PEO4. This demonstrates the link between identified COREs, lost and gained in the resistant PEO4 line, and cell line-specific chromatin states determined from histone modification ChIP-seq data.

Thus, the resistant line in all three pairs show enrichment of differential expression of genes near COREs with altered chromatin accessibility, and that COREs are associated with cell line specific enhancer elements defined by chromatin states. Importantly, within the set of genes near COREs gained in PEO4, we observe enrichment for genes involved in the Fanconi anemia/BRCA DNA damage response (DDR) pathway [28] that are associated with clinical response to DNA damaging cytotoxics in breast cancer (p< 0.05, OR 4.8, CI 1.02–16.0, Fisher’s exact test), highlighting a potential mechanism through which altered CORE landscapes may drive the acquisition of resistance to DNA damaging agents. Furthermore, SOX9 has been recently identified as maintaining chemoresistance in ovarian cancer cells through altered accessibility at its associated enhancer [14], was among the genes near a CORE which was gained in the resistant PEO4 line (see Additional file 1: Fig. S6A).

We examined in more detail the expression of the genes included in the BIOCARTA_ATRBRCA curated list of genes involved in BRCA related cancer susceptibility. As shown in Fig. 5, expression of the 19/22 genes for which we had RNA-seq data was sufficient to separate PEO1 from PEO4 and PEA1 from PEA2 by unsupervised hierarchical clustering. In the PEO pair, 12 of these genes showed small (log2FC < |1|) but significant (FDR < 0.05) changes in expression in PEO4 compared to PEO1. These included RAD51, ATM, FANCA, HUS1, BRCA1 and CHEK1 which were downregulated in PEO4. In the PEA pair, 14 of these showed a similar change in expression, of which NBN, TREX1, ATM, CHEK1, FANCC, FANCG, BRCA1 and BRCA2 were downregulated (FDR < 0.05).

Fig. 5

Expression of HR genes in platinum-sensitive and resistant cell lines. Heatmap of Z scaled normalized counts for 19 genes in the BIOCARTA_ATRBRCA gene set, from RNA-seq from three replicates of PEO1, PEO4, PEA1 and PEA2. Dendrogram shows hierarchical clustering using the ward.D2 algorithm, based on the expression of these genes.

Genome mapping of Pt-adducts in sensitive and resistant lines

To analyse sites of platinum adducts at high resolution, we have adapted methods for exonuclease mapping of transcription factor binding sites and developed Pt-exo-seq [20, 21]. A schematic of the Pt-exo-seq approach is shown in Additional file 1: Fig. S1. Massively parallel DNA sequencing is used to identify positions where 5′ to 3′ exonuclease digestion is blocked by platinum-DNA adducts. Levels of platinum uptake and DNA adducts formed can vary between cell lines, with many platinum resistant cell lines showing lower levels of drug uptake than their sensitive counterparts [29]. While the Pt-sensitive PEA1 and Pt-resistant PEA2 lines show similar levels of Pt-adducts, as measured by inductively coupled plasma mass spectrometry (ICP-MS), when cells are treated over a range of platinum doses, the PEO4 and A2780/cp70 resistant cell line shows markedly fewer total platinum adducts than the sensitive PEO1 and A2780 cell lines (Fig. 6a). Therefore, in order to compare distribution of Pt-adducts across the genomes of PEO1 and PEO4 we have used doses of cisplatin treatment that induce approximately equivalent levels of total Pt-adducts: 16 µM and 32 µM, respectively.

Fig. 6

Platinum adducts acquisition in PEO1 and PEO4. a PEO1 incurs more adducts than PEO4 at the same dose of cisplatin, as measured by ICP-MS. b The difference in adduct formation is reflected in levels of enrichment for platinum target purine dinucleotides at the 5’ end of Pt-exo-seq reads. PEO1 and PEO4 were treated with concentrations of cisplatin that ICP-MS indicated would induce equal amounts of damage.

Cisplatin forms covalent bonds to the N7 positions of purine bases leading to DNA cross-links. The majority of intrastrand cross-links are 1,2-d(GpG) cross-links, followed by a small number of 1,2-d(ApG) cross-links [30]. We would predict therefore that the majority of fragments detected by Pt-exo-seq will have enrichment for platinum target purine dinucleotides at the 5′ end of Pt-exo-seq reads. Indeed at doses of 100 µM cisplatin in A2780, we observe a dose dependent increase in the exonuclease resistant fraction of DNA and a 2.5-fold increase in purine dinucleotides at the terminal position (Additional file 1: Fig. S7), which compares favourably to the twofold enrichment in dinucleotides previously observed using the Damage-seq assay at 200 µM cisplatin. [22, 23]. As shown in Fig. 6b,at less highly toxic doses of cisplatin in the HGSOC lines, the 5’ end of Pt-exo-seq reads from DNA isolated from PEO1 and PEO4 cells treated with cisplatin are also enriched for purine dinucleotides compared to reads from untreated cells. At doses of 16 µM cisplatin, there was less enrichment for purine dinucleotides detected following Pt-exo-seq in the resistant PEO4 line compared to the sensitive PEO1 line, while at 32 µM cisplatin treatment of PEO4 there were approximately equivalent levels of enrichment for purine dinucleotides. Together these data support Pt-exo-seq as being able to detect the location of platinum adducts in genomic DNA.

Distribution of Pt adducts in sensitive and resistant lines

Genome-wide Pt-exo-seq-coverage was calculated in 1 Kb windows, as a log2FC over the untreated signal. This window size falls in the same order of magnitude as the median size of MACS2 called ATAC-seq peaks and provides sufficient depth for the signal-to-noise ratio to be calculated. These profiles were used to compare platinum-DNA adduct formation in PEO1 and PEO4 cells, at equal doses (16 µM vs 16 µM), and doses inducing equal amounts of damage (16 µM vs 32 µM respectively). Mean Pt-exo-seq signal was calculated for each window, across three replicates for each line and treatment condition, and used to calculate a log2FC in damage, after taking account of the untreated controls, between PEO1 and the 32 µM treated PEO4 cells. We found 1170 1 Kb windows showing a log2FC > ± 2 and moderated t test FDR < 0.05. These were then separated into two groups based on direction of change in damage in PEO4 compared to PEO1; increased or decreased frequency of DNA adducts. Hierarchical clustering was then performed on the normalized read counts, using these windows, which separated the PEO1 16 µM replicates from the PEO4 replicates treated at 16 µM and 32 µM, with the PEO4 replicates clustering by dose (Fig. 7a). The PEO1 and PEO4 also separated from each other when clustering on all windows analysed (Additional file 1: Fig. S8).

Fig. 7

Platinum adduct distribution is altered in the PEO4 platinum resistant cell line. a Heatmap of Pt-exo-seq coverage in differentially damaged windows, for PEO1 and PEO4. 1170 1 Kb windows showing mean log2FC in Pt-exo-seq signal > ± 2 and FDR < 0.05, between 16 µM cisplatin treated PEO1 replicates and 32 µM treated PEO4 replicates (5 h treatment). Dendrogram shows result of unsupervised hierarchical clustering of samples by Spearman’s correlation coefficients. Blocks on left show separation of windows into those showing either increased or decreased damaged in PEO4. b Chromatin accessibility in windows showing more or less platinum adduct formation in PEO4 compared to PEO1. PEO1 cells were treated with 16 µM cisplatin and PEO4 treated with 32 µM for 5 h. Difference in means: Decreased damage = 0.5, Increased damage = 0.90 (***p < 0.001, Wilcoxon test). c Enrichment for classes of genomic elements in windows showing differential damage between PEO1 cells treated with 16 µM cisplatin and PEO4 treated with 32 µM for 5 h. Odds ratios for only genomic classes showing significant enrichment (p < 0.05) in sets of windows showing increased or decreased damage in PEO4 compare to PEO1. Error bars show 95% confidence interval (Fisher’s exact test).

Thus, even though cisplatin is a relatively non-specific DNA damaging agent and might be expected to induce adducts at purine dinucleotides throughout the genome, there are clear and consistent differences in distribution of Pt-adducts between these matched sensitive and resistant ovarian cell lines and this is independent of the overall level of adducts formed.

Effect of genomic context on Pt-adduct distribution

Chromatin accessibility was investigated in relation to cisplatin-DNA adduct formation. ATAC-seq coverage was calculated for each of the 1 Kb windows previously defined by Pt-exo-seq as having different levels of Pt-DNA adducts between PEO1 16 µM and PEO4 32 µM replicates. Thus, we aimed to assess changes in chromatin accessibility between PEO1 and PEO4 at regions that showed differential DNA Pt-adduct levels. Regions incurring fewer Pt-adducts in PEO4 showed a significant reduction in accessibility as measured by ATAC-seq in PEO4 compared to PEO1 (p < 0.001, t-test) (Fig. 7b). Surprisingly, when carrying out the same analysis for the smaller set of windows showing increased damage in PEO4, an even greater reduction in accessibility was observed (p < 0.001, t-test).

Odds ratios were then calculated for enrichment of each class of genomic element in the sets of regions showing increased or decreased damage in the 32 µM treated PEO4 samples compared to PEO1 (Fig. 7c). Significant underrepresentation of regions annotated as CpG island, promoter-TSS and exon was detected in the set of regions showing increased damage in PEO4 (CpG OR 0.40, CI 0.16–0.83, promoter-TSS OR 0.56, CI 0.36–0.83, exon OR 0.24, CI 0.03–0.86, Fisher’s exact test). Conversely, regions showing decreased damage in PEO4 were over-represented for these same annotations. Intergenic annotated regions were enriched in the set showing increased damage in PEO4 (OR 1.80, CI 1.41–2.29). The decreased damage around promoters in the resistant line and the increased damage at intergenic regions suggest differential adduct formation in resistant cells is related to different genomic contexts which influence the rate of occurrence, or repair of adducts, or have implications on how they are tolerated by the cell.

To examine this further, we compared Pt-exo-seq data at COREs shared between PEO1 and PEO4, those only in PEO1 (i.e., lost in PEO4) and those only in PEO4 (i.e. gained in PEO4) (Additional file 1: Fig. S9). Significantly more platinum adducts were detected in the COREs in the resistant PEO4 line compared to COREs in the sensitive PEO1 line. Interestingly, the COREs lost in PEO4 incurred more adducts than COREs that are maintained in PEO4 (i.e. the shared COREs). This is consistent with a hypothesis whereby DNA damage at chromatin in less accessible COREs are better tolerated and could be acting as DNA damage ‘sinks’ leading to drug resistance.

Effect of modulating chromatin accessibility on cisplatin-induced DNA adduct levels

Treating ovarian tumor cells with the HDAC inhibitor Vorinostat causes a global increase in chromatin accessibility (Fig. 8a, b). 20 µM Vorinostat treatment for 24 h before chromatin extraction caused a significant increase in fragments of a length associated with mono-, di- and tri-nucleosomes following MNase digestion, along with a corresponding decrease in larger fragments, consistent with a more open chromatin conformation (mono-, p = 0.007, di-, p = 0.006, tri-p = 0.012, larger, p = 0.007, t test) (Fig. 8c). The same Vorinostat treatment followed by cisplatin treatment of cells, leads to a higher level of Pt-induced adducts, as measured by ICP-MS, compared to vehicle pretreatment (p < 0.001, t-test) (Fig. 8d). These data are consistent with a more open chromatin conformation increasing platinum DNA adduct formation. While Vorinostat treatment enhances the sensitivity of cells to cisplatin, this enhanced platinum sensitivity is observed in both platinum-sensitive PEO1 and platinum-resistant PEO4 HGSOC cells (Fig. 8e), suggesting a lack of specificity to drug-resistant cells.

Fig. 8

Effect of Vorinostat treatment on global chromatin accessibility and Pt-adduct formation. a TapeStation trace and b gel images showing, respectively, peaks and bands corresponding to mono-, di- and tri-nucleosome-associated DNA fragments following MNase digestion. Samples were treated with 20 µM Vorinostat or vehicle for 24 h. c Quantification of percent integrated area associated with mono-, di-, trinucleosomal and larger fragments following 24 h treatment with Vorinostat or vehicle, from b. There was a significant difference in integrated area between the vehicle and Vorinostat treated cells. Mono-, p = 0.007, di-, p = 0.006, tri-p = 0.012, larger, p = 0.007, t-test. c Levels of DNA platination as measured by ICP-MS, when PEO4 cells were treated with 0. 16 or 32 µM cisplatin for 5 h following 24 h pre-treatment with 20 µM Vorinostat. Error bars show SEM. ***= p < 0.001, t-test. n = 3. d Cell survival of PEO1 and PEO4 cells on exposure to increasing concentrations of cisplatin for 24 h, following pre-treatment with vehicle, or 20 µM Vorinostat for 24 h. Error bars show SEM. n = 3.

Discussion

Cytotoxic, DNA targeting, chemotherapies, such as platinum based chemotherapy, have major impact on improving patient survival and are still the first line of treatment for many cancers. However, drug resistance is a major clinical problem leading to poor survival for patients. A diverse range of drug resistance mechanisms have been observed experimentally [31]. However, with the exception of THBS1 and depending on the statistical parameters used, we observe relatively few consistent single gene changes in expression or chromatin conformation between the sensitive and resistant matched pairs. We also detected few consistent changes at the gene set level, but those we did detect involved genes from the apoptosis gene set such as EREG, implicated in stimulation of cell proliferation [32] and CCL5 in the interferon gamma gene set, which has been shown to promote invasion and migration of ovarian cancer cells [33]. However, the mechanisms through which these might contribute to drug resistance in ovarian cancer remains to be determined. In contrast to this, the ATAC-seq data of the HGSOC resistant lines in Fig. 1 shows that the resistant lines clusters together rather than with their respective sensitive line, suggests there are common underlying epigenetic mechanisms involving changes in chromatin. The changes in chromatin conformation observed occur more frequently at intergenic regions rather than gene promoters and result in chromatin accessibility changes at COREs which, given their association with active chromatin states serve as a good proxy for enhancers, and hence change the enhancer epigenomic landscape.

Our data argue that epigenetic changes at COREs during development of drug resistance is a common underlying mechanism that can drive expression changes in pathways such as those involved in homologous recombination repair leading to drug resistance. These findings reflect the genome wide chromatin reorganization that underlies resistance to endocrine therapies in breast cancer, or the emergence of repressed chromatin states at intergenic, repetitive regions in non-small lung cancer cells exposed to lethal concentrations of EGFR kinase inhibitors, as well as cisplatin [34, 35]. Such epigenetic changes may alter different resistance pathways and gene expression in a stochastic manner which can be selected on during tumor evolution. We observe enrichment for genes involved in the Fanconi anemia/BRCA DNA damage response pathway within the set of genes near COREs gained during acquisition of drug resistance [28] and observe down regulation of their expression in the resistant PEO4 and PEA2 lines. This is consistent with loss of HR repair leading to platinum sensitivity.

Integrative genome-wide epigenomic and transcriptomic analyses of in vitro derived platinum-sensitive and platinum–resistant ovarian lines identified key distal enhancers associated with platinum resistance and identified SOX9 as a critical CORE regulated transcription factor that plays a critical role in chemoresistance in vitro of ovarian cancer cell lines [14]. We also observe differences in COREs in the in vitro derived A2780 and A2780/cp70 pair of lines, although these cluster separately from the HGSOC, in vivo, derived resistant lines (Fig. 1). In our data, SOX9 is in the set of genes near COREs gained in the resistant PEO4 line, which is consistent with the data from Shang et al. showing increased H3K27ac around SOX9 in chemoresistant ovarian cancer cells [14].

THBS1 was the only gene consistently defined as downregulated in resistant lines using all methods of analysis of the RNAseq data. THBS1 is an important component of the extra-cellular matrix and has an important role in cancer development and regulating tumor cell behaviour [36]. Although no evidence has been shown for THBS1 having a direct role in resistance to platinum drugs, it does have functions in tumor vasculaturization, modulation of immune responses and is a pro-apoptotic factor, which may be important during tumor evolution during treatment. Further, ovarian cancer patients with high THBS1 expression and treated with platinum-based chemotherapy have longer survival [37] and is an independent prognostic factor in multivariate analysis.

In order to explore further the influence of chromatin organization in cisplatin resistance we developed and applied Pt-exo-seq, a method for mapping cisplatin adducts genome-wide at base pair resolution. Genome-wide adduct formation distinguished the resistant PEO4 line from its isogenic-sensitive counterpart, and was associated with the altered chromatin landscape we detected, highlighting altered adduct formation and distribution as a defining characteristic of these resistant ovarian cancer cells. Regions of the genome with fewer observed adducts in the resistant line were primarily around gene promoters, while regions with increased adducts were enriched for intergenic loci, and included COREs which were reprogrammed in the resistant lines.

A limitation of the Pt-exo-seq study is that cells were treated for 5 h to allow platinum cross-links to form before DNA isolation for Pt-exo-seq and some DNA repair of the mono- or cross-linked adduct is likely to have occurred in this time. Transcription-coupled nucleotide excision repair may be responsible for the reduced persistence of adducts near gene promoters; a view supported by the difference in adduct formation at transcription start sites (TSS) of genes showing equally high expression in PEO1 and PEO4 (Additional file 1: Fig. S10). Thus, variability in location of platinum in adducts detected in the genome may be due to differences in rate of DNA repair across the genome [22, 23], although may also be due to differences in chromatin accessibility for platinum damage. The current study has focussed on analysis of patient-derived platinum resistance, but the use of cell lines deficient in specific DNA repair pathways would help address some of these limitations.

It is well established that there is a difference in platinum DNA damage induction and repair in naked versus DNA bound to nucleosomes. For instance, studies on platinum damage at reconstituted chromatin show increased adduct formation predominantly in the nucleosome core. However, what the present study suggests is that it is not only the presence of nucleosomes that is important, but also whether they are in an open or closed conformation. Thus, regions of more compact chromatin had more adduct formation in resistant than in the sensitive HGSOC line, and this was particularly pronounced in intergenic regions. Furthermore, when examining COREs which were shared between PEO1 and PEO4, significantly more platinum adducts were detected in the resistant PEO4 line compared to COREs in the sensitive PEO1 line. Similarly, COREs which are lost in PEO4 incurred more adducts that the COREs which are maintained. This suggests that more platinum damage occurs at COREs in chromatin regions which become more nucleosome dense with the development of drug resistance, and that adducts at these loci are more tolerated by the cell, for instance due to increased error-prone translesion synthesis. While epigenetic approaches to overcome drug resistance have been proposed and are undergoing clinical evaluation, they have generally failed to show clear patient benefit or in some cases have increased toxicity [38, 39]. The use of non-specific epigenetic therapies that affect epigenetic marks throughout the genome will affect normal and tumor cells. This may explain the difficulties in improving therapeutic index in clinical trials of epigenetic therapies to overcome drug resistance. More targeted approaches based on reactivation of specific genes by CRISPR-epigenetic modulator fusions are now being explored [40] but are unlikely to overcome the diverse range of mechanisms involved in drug resistance if they only target specific genes. Targeting chromatin conformation in specific genomic contexts, rather than genome-wide, will be an important future strategy for epigenetic therapies to reverse drug resistance.

Conclusions

We have shown that chromatin changes, particularly at intergenic regulators of gene expression, are associated with in vivo derived drug resistance and Pt-adduct distribution in HGSOC. This has important implications for understanding the mechanisms of how tumor cells can adapt to treatment leading to drug resistance and treatment failure. These data suggest that resistance is driven by epigenomic changes which alter patterns of gene expression primarily through chromatin change at intergenic, CORE regions rather than gene promoters. Therefore, new approaches are required that target the epigenome in terms of chromatin conformation in specific genomic contexts, rather than either the entire genome or single gene targeting.

Materials and methods

Cell lines and cell culture

The derivation of the cell lines has been previously described [24, 25]. The identity of the lines was confirmed by STR profiling (Genetica). All cell lines were grown in RPMI-1640 (Sigma Aldrich) supplemented with glutamine (2 mM) (Gibco) and 10% fetal bovine serum (First Link). Cells were tested for mycoplasma monthly using the MycoAlertTM Mycoplasma Detection Kit (Lonza). Cell viability was assessed using MTS assay: CellTiter 96® AQueous One Solution Cell Proliferation Assay (Promega).

ATAC-seq

Chromatin was extracted and digested as described [41]. Libraries were sequenced on HiSeq 2500 (Illumina) using 100 bp paired-end reads. Quality assessment of libraries was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Duplicate reads were removed using the MarkDuplicates tool from the Picard tools suite (https://github.com/broadinstitute/picard). Reads were then aligned to the genome using Bowtie2 [42] with the parameters -X2000 and -m1. This ensured that fragments up to 2 kb were allowed to align (− X2000) and that only uniquely aligning reads were collected (− m1) as described [43]. Fragment size data were collected using the CollectInsertSizeMetrics tool from the Picard suite. Reads were filtered against genomic blacklist regions, including mitochondrial sequences, defined by The ENCODE Project Consortium [44].

RNA-seq

Libraries were prepared using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB), according to the manufacturer’s instructions. Using the featurecounts function from the RSubread package in R [45], aligned reads were assigned to genes as counts, using a reference file of all known genes on the hg19 genome assembly obtained from the University of California Santa Cruz Genome Browser [46]. Normalization was carried out and differential gene expression was calculated using the limma package [47]. RPKM values were calculated using the rpkm function from the limma package.

Pt-exo-seq

Libraries for Pt-exo-seq were produced using a novel protocol for detecting platinum damage genome wide adapted from [20, 21]. Following treatment with cisplatin (Hammersmith Hospital Pharmacy), DNA was isolated from cells using the Gentra PureGene kit (Qiagen). 1 µg purified DNA was sheared using Bioruptor Pico (Diagenode) to sizes of ~ 150 bp. Fragment sizes were confirmed using TapeStation (Agilent). Fragmented DNA was then incubated at 30 °C for 30 mins, with dNTPs, ATP, T4 DNA polymerase, Polynucleotide Kinase and DNA Polymerase 1 (Klenow fragment) to generate blunt ends before the reaction was stopped by incubation at 75 °C for 20 mins. Fragments were then incubated with the annealed P7 adapters along with ATP, and the T4 DNA ligase. Adapter ligated fragments were purified and incubated with the phi29 DNA polymerase, along with dNTPs, to allow the overhang between the fragment end and the adapter to be filled in. Adapter ligated DNA was incubated with Lambda exonuclease for 1 h at 37 °C. RecJF and single stranded DNA binding protein (SSB) were added with Lambda exonuclease to maximize digestion efficiency. The platinum molecules which were still present in the fragment, adapter ligated, digested DNA were then displaced by incubation with sodium cyanide at 65 °C for 2 h. Complementary strand synthesis was carried out using adapters complementary to the P7 adapter sequence, along with the phi29 DNA polymerase and dNTPs. Finally, the P5 adapter was ligated to the purified, double stranded DNA by incubation for 60 min at 25 °C and then 10 min at 65 °C to inactivate the enzymes. Purified libraries were then amplified by PCR. The enzymes and oligonucleotides used for preparation of Pt-exo-seq libraries are listed in Additional file 1: Tables S5 and S6. Libraries were sequenced on HiSeq 2500 using 100bp single read sequencing. Library quality was assessed by FastQC before adapter trimming with BBDuK, and alignment to the hg19 genome assembly using bowtie2. Reads were filtered against the DAC Encode Blacklisted regions [44] and deduplicated as described for the ATAC-seq libraries. Aligned, deduplicated reads were assigned to genomic loci using a 1 Kb sliding window, moved 900 bp along the genome. Coverage in these windows was then calculated for each platinum treated sample, against the mean coverage of three mock-treated replicates, as a log2 ratio. Further details on assay optimization have been reported [19].

Statistical analysis

ATAC-seq

Counts of ATAC-seq reads in genomic windows were calculated using the featurecounts function from the RSubread package [45]. RPKM values were calculated using the rpkm function in the edgeR package [48]. Reads were further filtered to remove PCR duplicates using the MarkDuplicates tool from the Picard Tools Suite (http://broadinstitute.github.io/picard/). The filtered, deduplicated reads were used in analyses which involved clustering the samples based on chromatin accessibility profiles.

ATAC-seq reads were assigned to 1 Kb windows and these were filtered by requiring a depth > 1 RPKM and < 1 SD above the median depth for all windows analysed. Coverage was normalized using the voom method to estimate the mean-variance relationship of the log of these counts. This generated a precision weight for each window which could be analyzed using the limma Empirical Bayes analysis pipeline [47] in order to calculate log2FC and p values for changes in accessibility for each window. Windows were annotated according to their nearest genomic feature using HOMER [49]. Enrichment odds ratios and confidence intervals were calculated for each class of genomic element in sets of windows showing significantly increased or decreased accessibility, compared to all analysed windows, using Fisher’s exact test.

Coverage around the TSS’ of genes was calculated in the same way as for the sliding window approach described above, except windows were 2 Kb centered around the TSS of genes of interest. Once FPKM values had been calculated, they were compared between cell lines using Students t test.

ATAC-seq peaks were called in each ATAC-seq replicate using MACS2. These peaks were then used as input for the CREAM package (https://www.biorxiv.org/content/early/2017/11/21/222562), which calls clusters of cis-regulatory elements based on the genomic distribution of peaks. Genes within 100 Kb of called COREs were extracted from the Ensembl database using biomaRt. T statistics and odds ratios for differential expression of these genes were calculated from the RNA-seq data produced in this study.

Pt-exo-seq

Odds ratios for enrichment of platinum-target dinucleotides at the 5’ end of reads were calculated using Fisher’s exact test to compare platinum-target dinucleotide frequency at each position (− 5 to + 10 bp) in the read. Windows were filtered by intersection with the regions passing filters in the ATAC-seq data for the PEO pair (depth > 1 RPKM and < 1 SD above the median depth for all windows analysed). Differentially damaged windows were defined as having a Log2 fold change in Pt-exo-seq signal:background ratio > ± 2 and an FDR corrected p value < 0.05, based on a t test between three replicates for each condition. Windows were annotated according to their nearest genomic feature using HOMER [49]. Enrichment odds ratios and confidence intervals were calculated for each class of genomic element in sets of windows showing significantly increased or decreased accessibility, compared to all analysed windows, using Fisher’s exact test. Pt-exo-seq coverage around the TSS’ of genes was calculated in the same way as for the sliding window approach described above, except windows were 2 Kb centered around the TSS of genes of interest. Once RPKM values had been calculated, they were compared between cell lines using Students t test.

MNase laddering following Vorinostat treatment

PEO4 cells were exposed to 20 µM Vorinostat for 24 h, nuclei isolated and incubated in micrococcal nuclease at a final concentration of 0.02 units/µl, in MNase buffer (New England Biolabs) at 37 °C for 10 min before purification using the MinElute PCR cleanup kit (Qiagen). Isolated DNA was run on TapeStation High Sensitivity DNA tape (Agilent). The relative area of each peak was then used to calculate enrichment for mono-, di- and tri-nucleosome associated fragments. Differences in peak area were tested using students t test.

Quantification of DNA platination by inductively coupled plasma mass spectrometry (ICP-MS)

Cells were treated with cisplatin for 5 h, before DNA was extracted. Platinum adduct levels were measured by ICP-MS using an Agilent 7900 mass spectrometer equipped with a micromist nebulizer and a cooled spray chamber (Agilent Technologies, USA). The assay was calibrated using platinum standards (Specure®, Alfa Aesar, Ward Hill, MA) and quality control samples with known platinum concentrations (ClinChek® Serum Controls, RECIPE, Munich, Germany). Platinum adducts were measured as previously described [50] in no-gas mode. Each measurement was made at 6 points on the mass peak in triplicate.