Nucleotide-resolution DNA double-strand breaks mapping by next-generation sequencing

Method workflow

To avoid specificity issues associated with ChIP or TdT, we devised a protocol for specific and direct in situ single-nucleotide resolution labeling and capture of individual DSBs in nuclei purified from mammalian cells (Online Methods). Direct in situ labeling prevents labeling of DSBs artificially formed during gDNA extraction, thus minimizing the risk of false positives. After fixation to stabilize chromatin and prevent artificial DSBs, cells are lysed and shortly incubated with proteinase K to purify intact nuclei. As a proof of concept, we optimized this step using both human cell lines as well as mouse tissues from which single-cell suspensions can be easily obtained (Supplementary Fig. 1a). After purification of intact nuclei, DSBs are blunted, 5′-phosphorylated, and finally ligated to a biotinylated linker (proximal) using the highly specific T4 ligase enzyme, which can only ligate double-stranded but not single-stranded breaks. The linker forms a hairpin-like structure, and thus can be ligated to either a blunt DSB or to an identical linker molecule, but cannot form concatemers. The ligatable end of the linker consists of a barcode sequence marking the exact position of ligation followed by the XhoI recognition site. Genomic DNA (gDNA) is extracted and fragmented, and labeled fragments are captured by streptavidin. A second linker (distal) also containing the XhoI site is attached to the free extremity of captured genomic fragments, enabling PCR amplification and sequencing (Fig. 1a-b). The use of barcoded linkers is a powerful strategy to unequivocally mark DSBs, avoiding background subtraction procedures as in ChIP data analysis.

Method implementation and validation

To implement our DSBs direct labeling and capture procedure, we performed pilot BLESS experiments in HeLa cells and mouse B-lymphocytes, followed by Sanger sequencing and next-generation sequencing on the Roche 454 platform. We performed various controls to exclude substantial false positive labeling due to incomplete washout of proximal linkers, unspecific binding of gDNA to streptavidin beads or mispriming. Only by following the complete BLESS protocol, DNA fragments could be amplified and subjected to sequencing (Supplementary Fig. 1b-d). 88% ± 6.5% (mean ± s.d., n = 2) of Roche 454 barcodes reads contained both proximal and distal barcodes, whereas obly 1.5% ± 2% (mean ± s.d., n = 2) contained the proximal barcode on both ends (Supplementary Table 1). As a first proof of specificity, we searched for reads mapping in the immunoglobulin heavy chain (IgH) locus among sequences derived from activated mouse B-lymphocytes. Upon B-lymphocyte activation, DSBs are formed in the IgH donor Su region and the downstream acceptor S region, enabling antigen class switch ²⁰. Accordingly, the density of correctly barcoded reads within these regions was significantly higher than the average read density in the genome (2-fold enrichment, P = 0.02, hypergeometric test), even with the relatively modest throughput achievable with the Roche 454 platform.

To increase data throughput, we performed deeper sequencing of BLESS samples using Illumina GAII and HiSeq 2000 platforms (Supplementary Table 1). All sequencing data, including Sanger and Roche 454 sequences, can be accessed at http://www.breakome.eu. In single-end sequencing experiments, the proportions of proximal and distal barcodes among barcoded reads were greatly similar (proximal 52.3% ± 9.8%, and distal 47.7% ± 9.8%, mean ± s.d., n = 9). Pair-end sequencing confirmed that over 99% (99.3% ± 0.2%, mean ± s.d., n = 2) of BLESS barcoded fragments contained both proximal and distal barcodes, whereas less than 0.8% contained the same barcode on both ends (Fig. 1c and Supplementary Table 1). This result demonstrates that the false positive DSBs labeling rate in BLESS is lower than 1%. We initially deep sequenced HeLa cells – a model system for which a large amount of genome-wide data is available and in which telomeric ends have been well characterized ²¹. During BLESS, the 3′ G-overhang of unprotected telomeres – which resembles a DSB repair intermediate ²² – is trimmed down to the first nucleotide of the complementary C-strand, where the biotinylated linker is ligated (Supplementary Fig. 2a). Therefore, we expected accessible telomeric ends to be detected by BLESS. Accordingly, we retrieved telomeric reads derived from the C-strand, with CTAACC being the most frequent (73%) C-strand end, as previously reported ²¹ (Supplementary Fig. 2b). We also deep sequenced U2OS cells carrying an I-SceI transgenic cassette, after transfecting them with HA-I-SceI (SCE) to induce a single DSB per cell within this cassette (Supplementary Fig. 2c). The density of barcoded reads inside the I-SceI cassette was almost 13,000 fold higher than the average read density in the genome (P = 6 × 10⁻³⁰⁰, hypergeometric test), further validating our labeling method. These results demonstrate that direct in situ labeling of DSBs followed by next-generation sequencing is an effective strategy to identify DSBs at various genomic locations.

Genomic landscape of sensitivity to replication stress

Deep sequencing experiments revealed DSBs sparse throughout the genome, which from now on we refer to as “breakome”. Cells grown in cell culture carry a non negligible amount of DNA breaks caused by a combination of replication stress, physiological apoptosis, and damage induced by reactive oxygen species ¹². Indeed, even in the absence of any exogenous treatment HeLa cells carried a substantial burden of γH2A.X foci (mean = 3.8, s.d. = 6.7 foci per nucleus, n = 300) (Supplementary Fig. 3a-b). In order to single out breaks caused by replication stress, which are believed to be a major source of genome instability ³, we exposed HeLa cells to a dose of aphidicolin – a DNA polymerase inhibitor – that induces replication fork stalling without arresting progression in S-phase ²³. This treatment resulted in significant accumulation of breaks (P = 10⁻³⁴, Kolmogorov-Smirnov test), and increased the amount of labeled DSBs captured by BLESS (Supplementary Fig. 3a-d).

We analyzed aphidicolin sensitivity at various resolutions by comparing deep sequencing data from replicates of samples treated (A) or not (C) with aphidicolin (Supplementary Fig. 4a). We analyzed the data by comparing read numbers within genomic windows with constant mappable length in A vs. C samples, and calculated enrichment P values based on the hypergeometric probability distribution. We computed final Q values based on the Benjamini-Hochberg approach for multiple hypotheses testing ²⁴ (Online Methods). We mapped aphidicolin-sensitive regions (ASRs) at a resolution of 48 kilobases (kb) without and with correction for copy number variation effects due to karyotype and aphidicolin treatment (Supplementary Fig. 4b) which yielded 2,307 and 2,429 significantly correlated ASRs, respectively (P = 10⁻³²³, hypergeometric test). To calculate the false discovery rate (FDR) related to our approach, we analyzed reads that were mapped to the chromosome Y through sequencing errors (in HeLa cells no reads from chromosome Y are expected). At 48 mappable kb resolution, the calculated FDR was 0.3%. The full list of ASRs and Q values is available at http://www.breakome.eu.

ASRs were non-uniformly distributed along the genome, with an average of 3% of 48 mappable kb regions per chromosome sensitive to aphidicolin, except for chromosomes 5 and 7 which were significantly more sensitive (5%, P = 10⁻⁵, hypergeometric test) (http://www.breakome.eu). As a comparison, we applied BLESS to HeLa cells treated with neocarzinostatin – a DSBs-inducing drug presumed to yield a more “random” pattern of breaks. Neocarzinostatin-sensitive regions were significantly more uniformly spread along the genome (distance between consecutive sensitive regions: 0.54 ± 1.81 megabases (Mb) and 0.21 ± 0.77 Mb for aphidicolin and neocarzinostatin, respectively, median ± s.d. P = 10⁻¹¹⁸, Kolmogorov-Smirnov test) (Fig. 2a and b). To validate our findings, we compared several top-ASRs with regions displaying no appreciable aphidicolin effect using γH2A.X ChIP. We observed a strong concordance between the aphidicolin effects measured in targeted regions and BLESS results (Fig. 3a-b). Importantly, the genomic locations of ASRs mapped in different experimental replicates were significantly correlated, demonstrating the reproducibility of our method (Supplementary Table 2).

Characterization of aphidicolin-sensitive regions

It has been suggested that repetitive DNA sequences may favor fork stalling and collapse upon replication stress, causing DSBs to appear more frequently at certain genomic regions ^3,4. In particular, repeats prone to form hairpin-like secondary structures might cause the collapse of slowly moving replication forks by directly hindering their progression. To investigate the association between repeats and aphidicolin sensitivity, we applied RepeatMasker ²⁵ to compute the abundance of various DNA repeat families inside ASRs as compared to the rest of the genome (Online Methods). We detected a reproducible strong and significant enrichment in satellites (P = 5 × 10⁻¹³⁷) in particular of alpha-type (P = 8 × 10⁻¹⁹⁸, see Online Methods for derivation of P values), a class of repeats that forms hairpin-like secondary structures and is abundant peri/centromeric regions (Fig. 4a, Supplementary Fig. 5, and Supplementary Table 3). Accordingly, ASRs mapped at high resolution (250 mappable nucleotides) were concentrated in peri/centromeric regions (Supplementary Table 4). Another class of repeats, AT dinucleotides, has been associated with a particular group of genomic regions sensitive to replication stress induced by aphidicolin – common fragile sites (CFSs) ^26,27. While many CFSs were scored as sensitive to aphidicolin following our approach, AT repeats were significantly depleted in ASRs (P = 10⁻¹⁶, hypergeometric test).

During replication stress, slowly moving replication forks will have a higher chance of colliding with transcriptional forks ²⁸, resulting in accumulation of DSBs which may be detected by BLESS. Accordingly, upon aphidicolin treatment, we detected a prominent enrichment of DSBs in transcribed regions with the highest enrichment in coding regions (P = 10⁻¹⁰, hypergeometric test). We next analyzed genes and ranked them based on the computed probability to develop DSBs anywhere along their length. Top 20% aphidicolin-sensitive genes showed significant enrichment of many gene ontology (GO) terms, particularly related to cell death (P = 10⁻³, hypergeometric test). Gene sensitivity to aphidicolin was significantly associated with gene length (P = 10⁻¹⁸, hypergeometric test), in line with the observation that the probability of collisions between replication and transcription forks in active genes seems to increase with gene length ^28,29 (Supplementary Fig. 6).

Finally, we investigated if ASRs mapped by BLESS are also associated with genomic regions or genes frequently rearranged in human cancers. Replication stress-driven genomic instability has been observed in many tumors, where it is thought to be an important cause of cancer rearrangements ^30,31. We used data from a cohort of over 2,700 human cancers ³², and found a modest, but significant enrichment of regions displaying amplifications or deletions inside ASRs as compared to the rest of the genome (P = 0.005, derived as described in Online Methods) (Fig. 4b). We next compared aphidicolin-sensitive genes with the Cancer Gene Census ^33,34, a collection of over 400 well annotated cancer genes, the majority of which is involved in translocations. Cancer genes were more likely to overlap with ASRs than non-cancer genes (P = 0.04, hypergeometric test), and the fraction of genes with 5′ end in 2 Mb vicinity of the center of a 48 mappable kb ASR, and containing that ASR center inside, was higher for cancer genes than others (P = 0.02, hypergeometric test) (Fig. 4c). Among most aphidicolin-sensitive genes, cancer genes were overrepresented (P = 0.04, hypergeometric test), including prominent oncogenes like EGFR, MET, ABL1, and MLL, which are typically mutated by translocation or amplification (Fig. 4d). The full list of genes ranked by aphidicolin sensitivity, and GO analysis results are available at http://www.breakome.eu.

Cells, reagents, and immunocytofluorescence

To obtain primary mouse single-cell suspensions, we squeezed testes and spleens from C57BL6/J mice between two microscope slides in a Petri dish filled with trypsin. We flushed bone marrow out of femurs and tibias from the same animals using a syringe filled with trypsin. We purified and activated B-lymphocytes as previously described ³⁷. Prior to fixation for BLESS, we removed dead cells using a Ficoll gradient. We filtered cell suspensions through MACS Pre-Separation filters (Milteniybiotec), and then fixed them according to the BLESS protocol. We obtained IMR90 primary fibroblasts and HeLa cells from ATCC, and U2OS_DRH-1 cells from Y. Shiloh (Tel Aviv University). Andrew J. Pierce (University of Kentucky) kindly communicated details on the construction of U2OS_DRH-1 cells. We transfected pcBAS-I-SceI and pCAGGs plasmids (kindly donated by Y. Shiloh) into U2OS_DRH-1 cells using Fugene (Roche) following the manufacturer’s instructions. After BLESS, we cloned gDNA fragments into pEGFP-C1 (BD Biosciences). We applied aphidicolin (Sigma) onto cells at 0.4 μM for 18 h, and neocarzinostatin (Sigma) at 200 ng/ml for 45 min. We obtained oligonucleotide linkers from Sigma and annealed them in 1× T4 ligase buffer (NEB). Linkers and primers used are listed in Supplementary Table 5. We visualized γH2A.X foci by immunocytofluorescence (Millipore # 05-636), and counted them as previously described ³⁸.

Breaks Labeling, Enrichment on Streptavidin, and Sequencing (BLESS)

A detailed, step-by-step protocol to perform BLESS can be found in the supporting webpage http://www.breakome.eu. Briefly, to prepare purified nuclei for in situ ligation, we fixed five million cells as single-cell suspensions in growth medium with 2% formaldehyde for 30 min at room temperature, and then washed them one time in ice-cold 1× PBS. To prepare single-nucleus suspensions, we first lysed fixed cells in a buffer containing 10 mM Tris-HCl, 10 mM NaCl, 1 mM EDTA, 1 mM EGTA, 0.2% NP-40 pH 8 for 90 min at 4 °C, and then in a buffer containing 10 mM Tris-HCl, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 0.3% SDS pH 8 for 45 min at 37 °C. We resuspended lysed cells in 1× NEBuffer 2 (NEB) supplemented with 0.1% Triton X-100 and Proteinase K at 100 μg/ml final concentration. We mildly rotated cells for a short time at 37 °C (8 min for HeLa, 4 min for mouse B-lymphocytes), after which we transferred them onto ice. We quenched proteinase K by adding an equal volume of buffer supplemented with PMSF.

We washed purified nuclei twice in 1× NEBuffer 2 supplemented with 0.1% Triton X-100, and then once in blunting buffer (NEB) supplemented with 100 μg/ml BSA. We performed blunting using the Quick Blunting kit (NEB) according to the manufacturer’s instructions in a final volume of 100 μl, for 45 min at room temperature. Afterwards, we washed nuclei twice in 1× NEBuffer 2 supplemented with 0.1% Triton X-100, once in 1× T4 ligase buffer supplemented with 0.1% Triton X-100, and once in 1× T4 ligase buffer. We performed in situ ligation for 18-20 h at 16 °C in 25 μl final volume using 1.5 μl of T4 ligase (NEB) and 5 μl of 10 μM proximal linker previously annealed in 1× ligase buffer. After ligation, we washed nuclei three times ina high-salt buffer (W&B) containing 5 mM Tris-HCl, 1 mM EDTA, 1 M NaCl pH 7.5 supplemented with 0.1% Triton X-100. Afterwards, we extracted gDNA by incubating nuclei in 1× NEBuffer 2 with 0.5% Triton X-100 and Proteinase K at 200 μg/ml final concentration for one hour at 65 °C shaking, followed by isopropanol-ethanol purification. We fragmented purified gDNA for 18-20 h at 16 °C using 6 U of HaeIII (NEB) every million cells fixed.

To capture labeled DSBs, we rotated 20 μg of purified gDNA with 5μl of Dynabeads MyOne C1 (Invitrogen) in W&B buffer supplemented with 0.1% Triton X-100 for 30 min at +4 °C. Afterwards, we washed beads three times in W&B buffer supplemented with 0.1% Triton X-100, and them resuspended in 37 μl of 1× T4 ligase buffer. We added 10 μl of distal linker previously annealed at 10 μM in 1× ligase buffer, and 3 μl of T4 ligase to the beads, and the reaction was carried on for 16-18 h at 16 °C. After distal linker ligation, we washed beads twice in W&B buffer supplemented with 0.1% Triton X-100 at room temperature, and then digested captured fragments with I-Sce (NEB) in 25 μl final volume for 4 h at 37 °C. Afterwards, we centrifuged beads, and stored the supernatant at 20 °C until PCR.

We used all the volume of supernatant after I-Sce digestion to prepare multiple PCR amplification reactions (5 μl per reaction) using Phusion polymerase (NEB) and the appropriate primers pair depending on the downstream sequencing platform (Supplementary Table 5). We performed eighteen amplification cycles using conditions recommended by the manufacturer and T_a = 55 °C. To remove unused primers, we purified PCR products in gel using the DNA Gel Extraction kit (Qiagen). Before Illumina library preparation, we digested purified PCR products with XhoI (NEB) to cleave terminal I-SceI sequences derived from linkers, and again gel-purified them.

Next-generation sequencing

Sequencing was either outsourced (imaGenes GmbH, Berlin, Germany, and ServiceXS, Leiden, Netherlands) or performed in-house as summarized in Supplementary Table 1. We prepared samples for Roche 454 sequencing using indexing barcodes-containing primers during the PCR step in BLESS (Supplementary Table 5). We purified PCR products of size comprised between 300 and 800 nt in gel, and analyzed them on the 2100 Bioanalyzer (Agilent) prior to sequencing. For BLESS Illumina library preparation, we used the TruSeq DNA sample preparation kit v2 (Illumina) without DNA fragmentation and library size selection. For gDNA sequencing, we sheared gDNA with Covaris S220 AFA (Covaris) according to the manufacturer’s instructions prior to Illumina library preparation. We assessed library quality and quantity on the 2100 Bioanalyzer (Agilent) using the High Sensitivity DNA Kit (Agilent), and by qPCR with the Kapa Library quantification Kit (KapaBiosystems). We generated clusters on the Illumina flow cell using the automatic cBot station and the TruSeq PE Cluster Kit v3-cBot-HS. We carried out sequencing by synthesis on Illumina HiSeq 2000 system using the TruSeq SBS Kit v3-HS chemistry.

ChIP and qPCR

We performed ChIP assays as previously described ³⁹, with minor modifications. We purified immuno-precipitated and input DNA with phenol-chloroform, and analyzed it by real-time qPCR using primers listed in Supplementary Table 6. We compared the amount of DNA captured in untreated vs. aphidicolin-treated HeLa cells by qPCR followed by data analysis according to the ΔC_t method ⁴⁰, using the C_t values obtained for each primer pair in sample C1 as a reference for the C_t values obtained for the same primer pair in sample A1. We used three technical replicates for each sample. Primers are listed in Supplementary Table 7.

Computational analyses