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).
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).
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.
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).
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).
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.
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).
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.