Reduced local mutation density in regulatory DNA of cancer genomes is linked to DNA repair

Reduced local density of somatic mutations in DHSs

To examine mutation frequencies in regulatory DNA, we mapped DHSs genome-wide in 12 cancer cell lines, as well as normal cellular counterparts of major malignancies (see Methods). We then analyzed whole-genome sequencing data from 34 tumor/normal pairs from seven distinct datasets: small-cell lung cancer³, melanoma², 23 multiple myeloma⁵ (MM) samples, and 9 colon cancers¹³. We used published mutation data for small-cell lung cancer³ and melanoma cell lines² (http://icgc.org) and re-analyzed primary tumor data on multiple myeloma and colon cancer using MuTect ¹⁴ (http://www.broadinstitute.org/cancer/cga/mutect). These 34 cancer genomes contained 364,226 somatic point mutations in about 2.6 Gbp of sequence that could be uniquely mapped in the DHSs assay, i.e. density of 0.000139 per base-pair (bp).

We observed a substantial reduction in the frequency of somatic nucleotide substitutions in DHSs compared to the genome average (Fig. 1 and Supplementary Fig. 1). This reduction is highly significant and consistent across all tumors (P <10^-36, chi-square test). The reduction was most prominent in the core TF binding regions of DHSs marked by the maxima of DNase I cleavage intensity (Fig. 1).

We next confirmed that the reduction of frequency of somatic mutations in DHSs was not the result of confounding factors influencing local variation in cancer mutation density, nor the result of sequencing and mapping biases¹⁵. Confounding factors may include differences between intergenic regions and genes (including both exons and introns), distance from transcription start sites² (Supplementary Fig. 2), time of DNA replication during the S-phase¹⁶, distances to telomeres and centromeres, and local G+C content¹⁵. Relative density of somatic mutations also depends on sequence context, especially flanking nucleotides, and different tumors exhibit different context dependencies^2,3,13 (Supplementary Fig. 3). The relative density of mutations expected from the sequence context is higher in DHSs, magnifying our observation (P <5*10^-181).

We observed significant relative reduction of somatic mutations in DHSs located in both intronic and intergenic regions, and in DHSs proximal (<1 kb) to transcription start sites versus more distal DHSs (Supplementary Fig. 4). More notably, the reduction was evident in comparison to immediately flanking 1kb regions (P <2.2*10^-16, chi-square test in lung, MM and colon; P=7.21*10^-13 in melanoma). As such, the observed reduction in the density of somatic mutations cannot be explained by a regional factor acting over long ranges¹⁷ such as transcription or DNA replication timing.

To rule out biases related to sequencing and mapping, for two of the cancer types (colon¹³ and MM⁵) with available raw sequencing data we repeated the analysis restricting it to nucleotide positions with above 80% detection power based on sequencing coverage. This analysis confirmed that the density of somatic sequence alterations is significantly reduced in DHSs compared to 1kb flanking sequences.

To account collectively for all of the above potential confounding factors, we applied Poisson regression model¹⁸. DHSs remained a significant and substantial contributor to the local somatic mutation frequency on top of other factors, including DNA replication timing¹⁹, distance from transcription start sites, distance from the DHS itself, CpG islands status, G+C content, and region type (exonic, intergenic, intronic) (Supplementary Tables 1 and 2). Because our regression analysis included neighboring windows, short-range regional dependencies could potentially inflate statistical significance. We repeated the analysis using a small subset (20%) of spatially separated windows and confirmed that DHSs remain a highly significant contributor even if only 20% of data are used (Supplementary Tables 1 and 2).

Notably, the effect of chromatin accessibility is monotonic and continuous and thus does not depend on the specific thresholds used to define DHSs (Supplementary Fig. 5). Finally, DHSs mapped in potential cells or tissues of origin (e.g., lung tissue for lung carcinoma, etc; see Methods) substantially contribute to the regression model that already includes pooled DHSs from multiple cell types. This demonstrates that cell-selective chromatin architecture and not simply genomic location is the driving feature.

The observed reduction in the frequency of somatic sequence changes within DHSs might be explained by either reduced occurrence of somatic mutation or by the action of purifying selection. At present, purifying selection in cancers has not been carefully studied, so we lack information that would support or contradict the action of purifying selection. In general, population genetics and comparative genomics studies in a variety of organisms suggest that purifying selection is usually stronger in coding regions than in regulatory regions^20,21. To investigate the possible action of purifying selection, we compared relative mutation densities in regulatory and protein-coding regions. The average reduction (across cancers) in frequency of somatic mutations in coding sequences relative to flanking sequences is smaller than analogous reduction in DHSs (Supplementary Fig. 6). Furthermore, the observed reduction of mutation frequency in exons may not necessarily represent the action of purifying selection. The frequency of missense mutations is not lower (and is even apparently higher) than frequency of synonymous changes (Supplementary Fig. 6). Thus, although it is possible that cancers may differ from evolving populations and cannot rule out the action of purifying selection in regulatory regions, we suggest that mutation attenuation plays a more important role.

Association with nucleotide excision repair

Relative mutation density depends on replication fidelity, levels of DNA damage or efficiency of DNA repair. The fact that the observed relative reduction of mutation density was highly limited to DHSs, makes it difficult to explain it by an increase in global replication fidelity. It is also unlikely that more accessible DNA at DHSs would be less prone to damage than less accessible DNA elsewhere. In fact, mutation frequencies observed in model organisms are reduced by positioned nucleosomes^22,23, while the effects of nucleosome positioning on somatic mutations in cancer are relatively small and differ in directions between various cancer types¹⁵.

Chromatin accessibility plays a major role in targeting nuclear proteins to regulatory DNA, and may provide a mechanism for preferential access by the repair machinery. Preferential activity of DNA repair proteins in accessible regulatory DNA may thus offer an explanation for the observed effect, analogous to the action of transcription-coupled repair in protein coding genes. We hypothesized specifically that potentiation of nucleotide excision repair (NER) and base excision repair (BER) by chromatin accessibility could be responsible for the observed relative reduction of mutation density at DHSs. The level of oxidative stress and the subsequent accumulation of lesions targeted by BER²⁴ is higher in malignant vs. normal cells²⁵. BER is an evolutionary conserved DNA repair pathway, which starts from the recognition and excision of various base lesions by specific DNA glycosylases, followed by the processing of the resulting AP sites and then DNA repair synthesis and ligation. Direct access of glycosylases to DNA lesions is pivotal for this repair process²⁶. Not surprisingly, therefore, BER complexes preferentially assemble in non-nucleosomal regions in response to oxidative stress²⁷, which naturally targets them to DHSs.

The NER pathway consists of two converging branches: global genome repair and transcription-coupled repair (TCR). DNA damage is first recognized by the XPC complex, DNA duplex is opened by XPD and XPB helicases, followed by the incision of the damaged strand XPF-ERCC1 and XPG nucleases, then gap filling by the replication polymerases and, finally DNA ligation²⁸. A priori, NER machinery shall be able to correct damaged DNA regardless of its chromatin state. A fully assembled NER complex, however, has a footprint of ~100 bps in DNA, which is significantly longer than the length of an internucleosomal linker. As a result, chromatin structure inhibits functional NER complex assembly and function^29-31. In case of TCR, this problem is circumvented by the fact that DNA damage is first sensed by the RNA polymerase followed by NER recruitment to an already unraveled chromatin requires by the CSB and CSA proteins (reviewed in ref. 6). This is not the case for global genome repair, and the problem is additionally acerbated by the fact that the damage sensor, the XPC complex, cannot bind to DNA adducts embedded in nucleosomes³². Thus, both lesion recognition and repair complex assembly may be potentiated in accessible chromatin^29,33. The fact that roughly half of DHSs lie in intergenic regions suggests GGR as the more likely candidate. Notably, nucleosomal chromatin appears to inhibit global genome repair function^29-31. Moreover, the XPC complex involved in damage recognition is inhibited from binding to lesions in nucleosomal DNA³². Thus, both lesion recognition and repair complex assembly may be potentiated in accessible chromatin^29,33.

Failure of NER predisposes to cancer. This is best illustrated by the extreme frequency of cancers caused by exposure to sunlight in Xeroderma pigmentosum (XP) patients, evidently due to their inability to repair UV photoproducts in DNA. Mutations in NER are commonly detected in melanomas. Similarly to XP, most somatic mutations observed in melanoma cells are C:G→T:A transitions caused by ultraviolet (UV)³⁴ damage, which are primarily repaired by NER.

BER and NER are high fidelity pathways as suggested by studies showing that deactivation of these pathways leads to increased chance of mutation from cryptic lesions or exogenous DNA damage^34-37. As was shown recently, NER defects lead to increase in density of C:G→T:A mutations under chronic low-dose UV radiation, conditions specifically relevant to melanoma³⁸.

We, therefore, reasoned that the relative reduction of mutation density in DHSs in individual melanoma genomes should parallel the integrity of NER pathway components. To test this, we analyzed 29 individual melanoma genomes sequenced at high coverage¹. Figure 2 and Supplementary Table 3 show the continuous dependence of C:G→T:A mutation densities on chromatin accessibility in melanocytes in both intergenic regions, introns and exons. We note that the observed effect is inconsistent with the recently discovered action of APOBEC proteins^39-41. First, APOBEC acts on single strand DNA and it is unlikely that inaccessible and untranscribed DNA would more readily adopt single strand conformation than DNA in DHSs and transcribed regions. Second, action of APOBEC would preferentially increase rates of C→T transitions and C→G transversions within TCA/TCT motif. Our observation is not confined to this motif, and we do not detect a parallel effect for C→G changes (Supplementary Fig. 7). Quantitatively, dependence of mutation density within TCA/TCT motif on number of DNase I cleavages is slightly lower (rather than higher) than the dependence of mutations within TCG/TCC motif (p-value for the interaction term in the regression analysis in Supplementary Figure 7 is 2.11×10^-5). This is inconsistent with the hypothesis that APOBEC action induces the observed dependency on chromatin accessibility.

This effect was far more pronounced in melanocyte chromatin vs. that of other cell. Furthermore, the signature of TCR activity is also observed by the higher density of C→T over G→A in the non-template strand² of exons and introns (Supplementary Fig. 8). This demonstrates that NER provide genes a multi layer protection against UV-light damage due to activity of TCR and the accessibility of chromatin to global genome repair.

Overall, 9 out of 29 melanoma genomes harbored non-synonymous mutations in NER genes. Notably, four melanoma genomes with the lowest levels of mutation reduction at DHSs all harbor mutations in NER genes^42-44 (P < 0.0237, Wilcoxon-Mann-Whitney test; Fig. 3). In three of these samples, mutation frequencies in DHSs return close to the genomic baseline. The presence of genomes with mutations in NER genes and reduced mutation frequencies in DHSs is not surprising because NER mutations may appear late in cancer development and some of the missense mutations may be functionally benign. This result implicates NER into observed reduction of mutation frequency in DHSs. It also provides an additional argument against selection explanation because purifying selection would not be expected to differentially impact melanoma samples.

Eight out of the nine samples with mutations in NER pathway genes harbored mutations in major components of the NER machinery (XPG/ERCC5, XPF/ERCC4 and LIG1). These lesions would be expected to compromise both NER and TCR, and therefore should affect mutation density in both DHSs and transcribed regions. Intriguingly, one of these samples also harbored a mutation in CETN2, which recognizes DNA distortions and therefore preferentially impacts GGR function over TCR⁴⁵. In agreement with this reasoning, three genomes carrying mutations in core subunits had markedly reduced or negligible reduction of mutation density in transcribed regions compatible with defective TCR (Fig. 4). Concordantly, in the genome carrying a mutation in the CETN2 gene, strong suppression of mutations in transcribed regions remained, implying that TCR function was not significantly compromised.

DNaseI hypersensitivity mapping

DNaseI mapping was conducted on cultured cancer cell lines, primary ex vivo hematopoietic cells, cultured primary cells, and isolated fetal tissues using appropriate nuclei isolation protocols (below), followed by a standard processing pipeline. Data from lines A549, HepG2, LNCap, CACO2, PANC1, CLL, K562, CMK, NB4, and MCF7 derive from Reference 8. Generation of new data from M059J (Glioblastoma), RPMI_7951 (melanoma), CD19 (B-cell), CD20 (B-cell), melanocytes, fetal lung and fetal intestine are described below.

Isolation of nuclei from cultured cancer cell lines

Cells were cultured in accordance with the detailed protocols provided at http://www.uwencode.org/protocols. To prepare nuclei, freshly grown cells were centrifuged at 500g for 5 minutes (4°C) in an Eppendorf Centrifuge 5810R, and washed in cold PBS (Cellgro/Mediatech Inc.). Cell pellets were resuspended in Buffer A (15 mM Tris-Cl pH 8.0, 15 mM NaCl, 60 mM KCl, 1 mM EDTA (Ambion/Life Technologies Corp) pH 8.0, 0.5 mM EGTA (Boston BioProducts) pH 8.0, 0.5 mM spermidine (MP Biomedicals, LLC) and 0.15 mM spermine (MP Biomedicals, LLC) to a final concentration of 2 × 106 cells/mL. Nuclei were obtained by dropwise addition of an equal volume of Buffer A containing 0.04% IGEPAL CA-630 (Sigma-Aldrich) to the cells, followed by incubation on ice for 10 min. Nuclei were centrifuged at 1,000g for 5 min and then resuspended and washed with 25 mL of cold Buffer A. Nuclei were resuspended in 2 mL of Buffer A at a final concentration of 1 × 107 nuclei/mL.

Isolation of nuclei from hematopoietic cells

CD19+ and CD20+ cells (separately) were isolated by immunomagnetic separation by the Large Scale Cell Processing Facility at the Fred Hutchinson Cancer Research Center from normal volunteer donors under an IRB-approved protocol. Cells were pelleted by centrifugation for 5 minutes at 500g at 4 °C. Cells were washed in ice-cold PBS, then resuspended to 5 million cells per mL in Buffer A. An equal volume of ice-cold 2X IGEPAL CA-630 solution (ranging from 0.02%-0.06%) was added and the tube was incubated for 5 - 6 minutes on ice to lyse the cells. Nuclei were pelleted by centrifugation for 5 minutes at 500g at 4 °C, resuspended in Buffer A and counted with a hemocytometer.

Isolation of nuclei from fetal tissues

Fetal lung and intestine tissues were obtained from morphologically normal fetuses by the Birth Defects Research Laboratory in the Dept. of Pediatrics at University of Washington, collected under an IRB-approved protocol. Tissue was minced, resuspended in cold 250 mM sucrose, 1 mM MgCl2, 10 mM Tris-Cl pH 7.5, with added EDTA Protease Inhibitor Cocktail (Roche Applied Science Corp.). Resuspended tissue from fetal brain, fetal lung, fetal kidney, and fetal adrenal was dissociated by slowly homogenizing with a Dounce homogenizer. Resuspended tissue from fetal heart or fetal intestine was dissociated in a gentleMACS Dissociator (Miltenyi Biotech Inc.). Following dissociation, all fetal tissues were filtered through a 100 uM filter, and nuclei pelleted by centrifugation 600g for 10 minutes. Pelleted nuclei were washed with Buffer A, resuspended in Buffer A and counted in a hemocytometer.

DNaseI mapping from isolated nuclei

Briefly, DNaseI digestion was performed as described in Reference 8, with minor modifications. Isolated nuclei (2 × 106) from suspension cells or dissociated tissue were washed with 15 mM Tris-Cl pH 8.0, 15 mM NaCl, 60 mM KCl, 1 mM EDTA pH 8.0, 0.5 mM EGTA pH 8.0, 0.5 mM spermidine and 0.15 mM spermine then subjected to DNaseI digestion for 3 min at 37 °C in 13.5 mM Tris-HCl pH 8.0, 87 mM NaCl, 54 mM KCl, 6 mM CaCl2, 0.9 mM EDTA, 0.45 mM EGTA, 0.45 mM Spermidine. Digestion was stopped by addition of 50 mM Tris-HCl pH 8.0, 100 mM NaCl, 0.1% SDS, 100 mM EDTA pH 8.0, 1 mM spermidine, 0.3 mM spermine. A range of DNaseI (Sigma-Aldrich), 10–80 U/mL) concentrations was used for each preparation of nuclei and the sample giving the optimum difference between DNaseI treated and untreated was used for sequencing library construction. DNaseI double-hit fragments were collected by ultra-centrifugation and gel-purified. Adaptors were ligated to the ends of purified fragments, and the resulting libraries sequenced on an Illumina Genome Analyzer IIx according to a standard protocol.Processing of DNaseI-seq data. 36-base reads with up to two mismatches were mapped to the human genome (GRCh37/hg19) using the sequence aligner BOWTIE. DHSs were identified using the Hotspot algorithm (8) at a false discovery rate (FDR) threshold of 5%. Genomic feature overlaps and distance calculations were performed using the BEDOPS suite of software tools available at http://code.google.com/p/bedops/.

Data availability

All DNaseI data used in this study have been released to the ENCODE Project repository or to the Roadmap Epigenomics Mapping Consortium data coordination center. These data have been deposited in GEO under accession numbers GSE29692 and GSE18927. Data are also available for download through www.uwencode.org/data and through the data links at www.epigenomebrowser.org.

Cancer datasets

We restricted our study to cancer genomes sequenced by Sanger Institute and by Broad Institute Whole genomes of COLO-829 and NCI-H209 cell lines that have been sequenced by Sanger Institute^2,3. COLO-829 cell line has been derived from metastatic tissue. NCI-H209, an immortal cell line of a small cell lung cancer derived from a bone marrow metastasis. Nine colon cancer genome¹³, 23 multiple myelomas⁵ and 29 melanoma genomes have been sequenced by Broad Institute. For the 4 Broad institute datasets we identified sites in the genome were we have over 80% power to detect mutations (this defined by at least 14 reads the covered his position in the tumor sample and 8 reads that cover this position in normal). About 81% of the bases in colon and MM genomes are sufficiently covered and about 86% of the melanoma genomes are sufficiently covered.

Annotations

gene and exon coordinates were retrieved for hg19 from UCSC genome browser. For the flank analysis we took 1000 bps windows around the DHSs, since some of flank regions overlapped with DHSs we removed the DHS defined sites from the original set of flanks.

Density of mutations

We calculated the number of mutations per bps that can be mapped uniquely by the DHS essays.

Poisson regression

In order to test whether DHS regions have an additional impact on rates on top of other cofounding factors we used multivariate Poisson regression¹⁸. We divided the genomes into non-overlapping 400 bp windows. Every bin was classified as intergenic, intronic and exonic regions. The windows were also classified as DHS regions when overlap at least 80bp of any DHSs. For each window we calculated the following quantities: GC content, CpG content, distance from the nearest transcription start site, distance from nearest CpG islands, distance from nearest DHSs, and mutation counts and coverage (i.e. the number of bases which we have 80% power to detect mutation) for Colon, multiple myeloma, and similarly the number of GC bases in a window that are covered in Broad Institute melanoma samples. Then using glm function in R we calculated the estimated rates¹⁸. We repeated the analysis for spatially separated windows comprising 20% of the data to ensure that possible short-range interdependence between neighboring windows does not artificially inflate statistical significance.