Genome-wide binding of the CRISPR endonuclease Cas9 in mammalian cells

Genome-wide binding of dCas9-sgRNA

To map dCas9 in vivo binding sites, we generated mESCs with a stably integrated vector encoding HA-tagged dCas9 (Fig. 1a), and performed chromatin immunoprecipitation followed by sequencing (ChIP-seq) with cells transfected with either no sgRNA or one of each of 4 sgRNAs (Phc1-sg1, Phc1-sg2, Nanog-sg2 and Nanog-sg3) _targeting the promoters of Phc1 or Nanog, respectively. For each sgRNA, we observed ∼100 fold enrichment for dCas9 at the on-_target site compared to flanking regions, and the spatial resolution is sufficient to distinguish between two binding sites separated by 22 base pairs (bps) (Nanog-sg2 and Nanog-sg3) (Fig. 1b).

Using the standard ChIP-seq peak-calling procedure MACS³⁴ – comparing immunoprecipitated material and input (whole cell extract) DNA – we identified between 2,000 and 20,000 peaks in each sequencing library (Supplementary Fig. 1a). Cells expressing dCas9 but not transfected with sgRNAs (dCas9-only ChIP) exhibited 2,115 peaks. Most (77%) of the peaks detected in the dCas9-only ChIP were also detected in libraries prepared from dCas9-sgRNA immunoprecipitations (Supplementary Fig. 1b). The peaks in dCas9-only ChIP were enriched in open chromatin regions (Supplementary Fig. 2a) and 41% contained GG/CC-rich motifs that closely resemble CTCF binding motifs (Supplementary Fig. 2b-d). Such peaks could either represent ‘sampling’ by dCas9 of accessible sites containing NGG³³, or transcription-dependent artifacts as previously reported for GFP ChIP in yeast³⁵.

To identify sgRNA-dependent dCas9 binding sites, we matched sequencing depth by randomly sampling an equal number of reads from all six libraries and then performed pair-wise peak calling with MACS using each of the other five libraries as the control; we only retained peaks that were enriched over all five controls (Fig. 1c). Only 3 background peaks were called using this approach for dCas9-only ChIP. The number of sgRNA-specific peaks varied substantially; for example there were nearly 6,000 peaks for Nanog-sg3 but only 26 peaks for Nanog-sg2 (Fig. 1b). Many of the off-_target peaks showed high binding levels, as defined by the peak height relative to on-_target peaks after subtracting dCas9-only reads. For example, there were 91 off-_target peaks with more than 50% of the binding level of the on-_target site for Nanog-sg3 (Supplementary Table 1). These results suggest that there are substantial numbers of off-_target binding sites and the majority of the dCas9-sgRNA complex binds outside the designed _target site.

A 5-nucleotide seed for dCas9 binding

Sequence motifs enriched within 50 bps of peak summits were identified using MEME-ChIP³⁶. The top motif found for each ChIP library matched the PAM-proximal region of the transfected sgRNA plus the PAM NGG (Fig. 2a and Supplementary Fig. 3). For 3 of the 4 sgRNAs, only PAM-proximal positions 1 to 5 in the _target DNA showed preference of base match to the guide (Fig. 2a). We therefore define position 1-5 as the ‘seed’ region of the sgRNA. For Nanog-sg2, the guide match extends to about 10-12 bases to the 5′ end, possibly due to the presence of multiple Us in the seed that lowers the thermodynamic stability of the sgRNA-DNA interaction. For Nanog-sg3 and Phc1-sg2, exact match to the 5-nucleotide seed followed by NGG (seed+NGG) within 50 bps of peak summits explained 96% and 97% of the peaks, respectively. When the 50 nucleotides flanking peak summits were shuffled preserving dinucleotide frequency, less than 5.7% of the shuffled sequences contained seed+NGG (Fig. 2a) for all four sgRNAs. Moreover, the seed+NGG sites were highly enriched at the center of the peak (Fig. 2a, right), suggesting the _target sites identified are directly bound by sgRNA-guided dCas9.

We found that seed+NGG alone is sufficient for Cas9 binding in vivo and in vitro. For example, there were 92 peaks in the Nanog-sg3 sample containing only seed+NGG matches, i.e. mismatches at all the other 15 positions. The strongest peak containing only seed+NGG showed 52% binding activity relative to the on-_target (Fig. 2b). In vitro gel shift assays confirmed specific binding to seed+NGG only substrates but with lower affinity than the on-_target site (Fig. 2c).

The peak motif analysis (Supplementary Fig. 3) revealed no enrichment of binding at seed sites followed by NAG, an alternative PAM previously reported to function in Cas9-mediated cleavage^9,16,27. For example, of all 996 (33%) Phc1-sg1 ChIP peaks without seed+NGG sites, only 18 had seed+NAG within 50-bp of the peak summit, even less than expected by chance (Supplementary Fig. 4). ChIP-seq in human HEK293FT cells transfected with dCas9 and the same sgRNAs used in a previous study⁹, where NAG cleavage was reported, also failed to detect binding at those NAG off-_target sites (Supplementary Fig. 5). In vitro we observed >10 fold decrease in affinity when NGG was mutated to NAG in the on-_target substrate (Supplementary Fig. 6). Our in vivo and in vitro binding data are consistent with previous in vitro cleavage data showing that NAG or other variants rarely function as PAMs under enzyme-limiting conditions²⁷.

Chromatin accessibility is the major determinant of in vivo binding

There are hundreds of thousands of seed+NGG sites in the genome for each sgRNA, for example, 621,651 for Nanog-sg3. To understand why only a small fraction of sites (<1%) were bound, we first looked for a correlation between the number of base match to the 20-nucleotide guide region and binding levels of ChIP peaks. Overall the correlation was very weak (Pearson correlation coefficient r=0.03, 0.12, 0.15 and 0.55 for Nanog-sg3, Phc1-sg2, Phc1-sg1 and Nanog-sg2, respectively (Fig. 3a and Supplementary Fig. 7)).

To identify determinants influencing in vivo binding, we applied a linear regression model of a set of sequence (mono- and di-nucleotide frequency), structural (melting temperature, DNA energy and flexibility³⁷) and epigenetic (chromatin accessibility as assayed by DNase I hypersensitivity (DHS)³⁸ and DNA CpG methylation³⁹) features around the seed+NGG sites for each sgRNA (Online Methods). We found that chromatin accessibility (DHS) is the strongest indicator of binding in vivo, explaining up to 19% of the variation in binding when considering all individual seed+NGG sites in the genome (Fig. 3b). The difference in the number of seed+NGG sites in DHS peaks (i.e. accessible seed+NGG sites) explained 92% of the variation in the number of dCas9 peaks among the four sgRNAs (Fig. 3c, n = 4, p<0.05, F-test). Although this is based on a limited set of sgRNAs, it suggests that it might be possible to predict the approximate number of off-_target peaks based on the seed sequence in cell types where chromatin accessibility data are available.

Previous data suggested that Cas9 cleavage activity is not affected by DNA CpG methylation⁹. However, for the 17% of seed+NGG sites in the genome that contain CpG dinucleotides within the 20mer guide match and NGG, CpG methylation became the strongest predictor of dCas9 binding and negatively correlated with binding (Fig. 3d, Supplementary Fig. 8a-b). In a regression model, adding CpG methylation to DHS for sites containing CpGs almost doubled the amount of variation explained (Supplementary Fig. 8c). Our data suggests that CpG methylation likely reflects an aspect of chromatin accessibility not fully captured by DHS or that when combined with extensive mismatches, CpG methylation may impede binding.

The correlation with chromatin accessibility suggested that dCas9 off-_target binding might preferentially occur at active genes. For Nanog-sg3, 70% of the off-_target sites were associated with genes, including 18% in promoter region (< 2 kb upstream of gene TSS), 6% near enhancer regions and 46% within genes (Fig. 3e). For example, an off-_target peak that co-localized with the Dusp19 gene TSS and a DHS peak showed 74% binding relative to the on-_target with only 7 base matches to NanogE sg3 (Fig. 3f).

Seed sequences influence sgRNA abundance and specificity

The Nanog-sg2 sgRNA had substantially fewer off-_target binding sites than predicted by accessible seed+NGG sites (Fig. 3c). Although the same amount of sgRNA plasmids were transfected, the abundance of Nanog-sg2 was more than seven fold lower than the other three sgRNAs as determined by Northern blot (Fig. 4a). The same pattern of sgRNA abundance was observed when cells were transfected with sgRNA expression plasmids without co-transfecting dCas9, although all four sgRNAs showed substantially decreased levels of abundance, consistent with previous reports that Cas9 stabilizes sgRNA in cells¹³.

To test if sgRNA abundance influences the number of off-_target sites bound, we repeated the ChIP experiments after transfection with various amounts of sgRNA plasmids. Northern blots confirmed the decrease in sgRNA when less plasmid was transfected (Fig. 4a) and we identified decreased numbers of peaks with decreased amounts of plasmid (Fig. 4b). When the level of Nanog-sg3 was reduced to a similar level as Nanog-sg2 (Fig 4a, comparing lane 13 to lanes 16 and 17), the number of peaks for Nanog-sg3 was still much higher than for Nanog-sg2, presumably due to the presence of more accessible Nanog-sg3 seed+NGG sites in the genome (Fig. 3c). When 0.02 μg plasmid was transfected, Nanog-sg3 RNA was barely detected (lane 14); the 122 peaks identified in this library showed low overlap (9%) with our previous Nanog-sg3 ChIP, suggesting these were mostly non-specific signals.

A comparison of the seed regions of the four sgRNAs suggested that UUU in the seed of Nanog-sg2 might be responsible for decreased sgRNA abundance and increased specificity, consistent with recent observation that U in PAM-proximal position 1-4 leads to low gene knockout efficacy¹⁴. Indeed, two mutations (U to G and U to A) in the Nanog-sg2 seed region that converted the seed (GUUUC) to the same sequence as the Phc1-sg2 seed (GGUAC), led to higher levels of sgRNA (sgRNA N2b in Fig. 4c). Considering the presence of GUUUUA adjacent to the seed and because sgRNAs are transcribed by RNA polymerase III which is terminated by U-rich sequences^40,41, we speculate that together with the downstream U-rich region, multiple Us in the seed might induce termination of sgRNA transcription. Consistent with this, three sgRNAs with seeds UUAUU, ACUUU and UUUUU also showed very low abundance (Fig. 4c, sgRNA P3, N5 and N6). When GUUUC was placed upstream of the seed thus away from GUUUUA in the sgRNA, the sgRNA was well expressed (sgRNA C4 in Fig. 4c).

One of 295 off-_target sites is mutated above background

To test if binding correlates with Cas9 nuclease-induced mutation, we examined the indel frequencies of the four on-_target sites and 295 selected off-_target sites by _targeted PCR and sequencing⁹. These sites were selected to cover a broad range of binding levels and numbers of mismatches to the sgRNA. We ranked all peaks by binding (background subtracted read counts) and for each binding level selected a peak with the fewest mismatches and another peak with most mismatches to the guide.

We determined the indel frequency of the 299 selected binding sites in wild type mESCs transfected with active Cas9 and each of the four sgRNAs, for three independent biological replicates (Supplementary Table 3). The level of Cas9 protein transiently expressed in the cells was 2.6 fold higher than in cells with stably integrated dCas9 used for ChIP (Supplementary Fig. 9a, comparing lane 1 to lane 8). The same ChIP and peak calling procedures in cells transiently transfected with dCas9 identified 2.7 times more Nanog-sg3 peaks (16,119 versus 5,957 in dCas9 stable cell lines), including 96% (85) of the 89 peaks selected for indel analysis. The amount of Cas9/dCas9 plasmids used for transfection is similar to levels used for genome editing applications by the field (Supplementary Fig. 9b).

Using our previously validated model⁹, the background indel frequencies due to sequencing noise were determined for each individual _target using two biological replicates transfected with only Cas9 but no sgRNA (“control”). Importantly the control samples showed no evidence of _targeted mutations by Cas9 (note that background indels in the absence of Cas9 might also occur). We manually reviewed sequencing alignments of all loci with indel frequencies above 0.03%. We found 12% to 37% sequencing reads from the on-_target sites contained indels, yet only one off-_target, which was from Nanog-sg2, was mutated at a frequency of 0.7% (Fig. 5). There was no detectable correlation between binding and indel frequency (sites in Fig. 5 are ranked by decreasing binding from left to right for each sgRNA). The selected sites include 7 of the top 10 (including all the top 6) and 36 of the top 50 Nanog-sg3 binding sites with the strongest ChIP signals, and 4 of the 8 Nanog-sg3 off-_target binding sites that have fewer than four mismatches to the sgRNA; none of these off-_target sites showed cleavage significantly above the background level.

Oligonucleotides

All oligonucleotides used in this study were purchased from Integrated DNA Technologies. Sequences are listed in Supplementary Table 2.

Cloning

A two-step fusion PCR was used to amplify Cas9 nickase ORF from pX335 vector (Addgene: 42335) and incorporate H840A mutation to create a nuclease deficient Cas9 (dCas9). This PCR product was inserted into the Gateway donor backbone pCR8/GW/TOPO to create pAC84 (Addgene: 48218). The dCas9 ORF in pAC84 was then transferred to a piggyBac-based destination vector pAC150 (PB-Lox-HygroR- Lox-4×HSInsulators-EF1a-DEST) by LR Clonase reaction (Invitrogen) to create pAC159 (PB-LHL-4×HS-EF1a-dCas9). The sgRNA expression cassette was amplified by PCR from pX335 vector and cloned into a piggyBac vector pAC158 (PB-neo-4×HSInsulators) to create pAC103 (Pbneo-sgExpression). sgRNA was then cloned into BbsI-digested pAC103 by oligo cloning method as described previously⁸. Cas9 transient transfection constructs consisted CBh-driven WT-Cas9 or Cas9-D10AH840A (dCas9) containing a C-terminal HA-tag.

Cell culture

V6.5 mouse embryonic stem cells (mESCs) were cultured in DMEM supplemented with 10% FBS, pen/strep, L-glutamine, nonessential amino acids and LIF. For generation of cells stably integrating dCas9, cells were transfected in a 6-well and selected using Hygromycin B at 100 ug/mL 24 hours post transfection, then raised to 150 ug/mL 48 hours post transfection. Cells were split onto 10 cm plates, and single clones were isolated, expanded, and used for all experiments described. HEK293FT cells were cultured as previously described⁹. All transfection were done with Lipfectamine 2000 (Invitrogen).

ChIP

Three million cells were seeded on to 10cm plates on day 1, transfected with sgRNAs plasmids (or together with HA-dCas9 plasmids) on day 2, transferred to 15cm plates on day 3, and crosslinked on day 4 with roughly 50 million cells. Crosslinking is done by adding 2 mL (0.1 volume) 37% formaldehyde to the plate, incubating at room temperature for 15min, and quenched by adding 1 mL 2.5M glycine. Cells were rinsed twice with cold PBS and scraped to collect in cold PBS. Cells were centrifuged at 1,350g for 5min at 4°C and washed again in cold PBS. Cells were flash frozen in a dry ice/ethanol mix and stored at -80°C. Cell pellet was resuspended in 5mL cold Lysis Buffer 1 (50mM HEPES-KOH pH 7.5, 140mM NaCl, 1mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Trition X-100, 1× Roche complete protease inhibitors), rotated at 4°C for 10min followed by centrifugation at 1,350g for 5min at 4°C. Pellet was resuspended in 5mL Lysis Buffer 2 (10mM Tris-Cl pH 8, 200mM NaCl, 1mM EDTA, 0.5mM EGTA, 1× Roche complete protease inhibitors), rotated at 4°C for 10min followed by centrifugation at 1,350g for 5min at 4°C. Nuclear pellet was resuspended in 2mL Sonication Buffer (20mM Tris-Cl pH 8, 150mM NaCl, 2mM EDTA, 0.1% SDS, 1% Triton X-100, 1× Roche complete protease inhibitors) and sonicated (60min total time, 30sec on, 30sec off, in 6 rounds of 10min) in a Bioruptor (Diagenode). Lysate was centrifuged in eppendorf tubes in a microfuge at 4°C a max speed for 20min. Supernatant was collected and 50μL of this was saved as input. Protein G Dynabeads were conjugated to 5 ug rabbit anti-rat antibody (Thermo) in 0.1M Na-Phosphate pH 8 buffer at 4°C with rotation followed by conjugation to 5 ug HA antibody (Roche 3F10, #11867431001). Beads were resuspended in 50 ul Sonication Buffer and added to samples to immunoprecipitate overnight. The next day, beads were washed twice in sonication buffer, once in sonication supplemented with 500mM NaCl, once in LiCl Buffer (10mM Tris-Cl pH 8, 250mM LiCl, 1mM EDTA, 1% NP-40) and once in TE + 50mM NaCl. Each wash was accomplished with rotation at 4°C for 5min. Chromatin was eluted at 65°C for 15min in Elution Buffer (50mM Tris-Cl pH 8, 10mM EDTA, 1% SDS). Input was combined with elution buffer and both input and IP crosslinks were reversed at 65°C overnight. RNA was digested with RNAse A at 0.2mg/ml final concentration (Sigma) at 37 °C for 2hrs and protein was digested with proteinase K at 0.2mg/mL final concentration (Life Technologies) at 55 °C for 45min. DNA was phenol:chloroform:isoamyl alcohol (Life Technologies) extracted and ethanol precipitated. Barcoded libraries were prepared and sequenced on Illumina HiSeq2000.

ChIP-seq data analysis

Reads were de-multiplexed and mapped to mouse genome mm9 using bowtie⁴⁷, requiring unique mapping with at most two mismatches (-n 2 -m 1 --best --strata). Mapped reads were collapsed and the same number of reads (about 9 million) was randomly sampled from each library to match sequencing depth. Peaks were called using MACS⁵⁰ with default settings. For each sample, the other samples are each used as a control and only peaks called over all five controls are defined as _target sites. To quantify relative binding strength, reads were first extended at the 3′ end to the average fragment length (d) estimated by MACS, and then the number of fragments (extended reads) overlapping with the seed+NGG region is counted and normalized by subtracting counts from dCas9-only control. If multiple seed+NGG match sites were found, the one with the highest relative binding was assigned to the peak.

Analysis on determinants of binding

Mouse ES cell DNase Hypersensitivity data (bigwig file and narrow peak file) were downloaded from UCSC genome browser hosting the mouse ENCODE project³⁸. DNA CpG methylation data was downloaded from GEO dataset GSE30202. Melting temperature (Tm) was calculated using the oligotm program in primer3 version 2.3.6. DNA stability and flexibility were calculated using a table of tetranucleotide scores derived from X-ray crystal structures in a previous study³⁷. The linear regression is performed by using the lm function in R, one feature a time to calculate the R² value for each feature.

Northern blot

Total RNA was isolated using TRIzol (Life Technologies) and 5 ug of total RNA was loaded on 8% denaturing PAGE. Northern blot was done as previously described⁹, using a probe _targeting the scaffold shared by all sgRNAs.

Protein Purification

Human codon-optimized Cas9 (Addgene plasmid 42230) was subcloned into a custom pET-based expression vector with an N-terminal hexahistidine (6×His) tag followed by a SUMO protease cleavage site. The fusion construct was transformed into E. coli Rosetta 2(DE3) competent cells (Millipore), grown in LB media to OD600 0.6, and induced with 0.2 mM IPTG for 16h at room temperature. Cells were pelleted, resuspended and washed with Milli-Q H2O supplemented with 0.2 mM PMSF, and lysed with lysis buffer (20 mM Trizma base, 500 mM NaCl, 0.1% NP-40, 2 mM DTT, 10 mM imidazole). The lysis buffer was supplemented with protease inhibitor cocktail (Roche) immediately prior to use. Whole lysate was sonicated at 40% amplitude (Biologics Inc., 2s on, 4s off) prior to ultracentrifugation (30,000 rpm for 45m). The clarified lysate was applied to cOmplete His-tag purification columns (Roche), washed with wash buffer 1 (20 mM Trizma base, 500 mM NaCl, 0.1% NP- 40, 2 mM DTT, 10% glycerol, 10 mM imidazole) and wash buffer 2 (20 mM Trizma base, 250 mM NaCl, 0.1% NP-40, 2 mM DTT, 10% glycerol, 50 mM imidazole). The 6xHis affinity tag was released via SUMO protease cleavage and bound protein was eluted with a linear gradient of 150 mM – 500 mM imidazole. Eluted protein was concentrated with Amicon centrifugal filter units with Ultracel membrane (Millipore) and stored at -80°C.

In vitro transcription

A T7 promoter forward oligo was annealed to an sgRNA template oligo by heating to 95°C for 3 min in 1× T4 DNA ligase buffer and then cooled at room temperature for 30 min. The annealed product was used as template and transcribed with MEGAshortscript T7 Kit (Life Technologies). RNAs were purified by MEGAclear Kit (Life Technologies) and frozen at -80 °C.

Gel shift assay

Single stranded DNA oligos of 50 nucleotides were purchased from IDT and PAGE purified. Double stranded substrate were generated by mixing 100 pmol each strand in water (10 ul total), heating to 95°C for 3 min and cool to room temperature. The substrates were then 5 end labeled with [γ-32P]-ATP using T4 PNK (New England Biolabs) for 30 min at 37°C, and free ATP removed by G-25 column (GE Healthcare). For each reaction, 100 nM Cas9 was mixed with a 1:4 dilution series of sgRNA (from 0 to 100 nM) in 1× NEBuffer 3 at 37°C for 10 min, and then about 0.5 nM labeled substrate oligos were added and incubated for 5 min at 37°C in a 10 ul reaction. Reactions were stopped on ice and added 1/2 volume of 50% glycerol. Samples were loaded on to 12% native PAGE and run at 300V for 2 hours at room temperature. Gels were visualized by phosphorimaging. Gel quantification is done with ImageJ. The fraction bound shown in Fig. 2c was calculated as the ratio of intensity from the specific binding band to the total intensity of the entire lane.

_targeted sequencing and indel detection

For replicate 1, cells were seeded in 6-well plates (300,000 cells per well), transfected with 2 ug sgRNA plasmid, 2 ug Cas9 plasmid, using 10 ul Lipofectmine 2000 reagent per sample for 3 hours. For replicate 2 and 3, 50% more plasmids were used. DNA was extracted and selected _target sites were PCR amplified, normalized, and pooled in equimolar proportions. Pooled libraries were denatured, diluted to a 14pM concentration and sequenced using the Illumina MiSeq Personal Sequencer (Illumina). Sequencing data was demultiplexed using paired barcodes, mapped to reference amplicons, and analyzed for indels as described previously⁹.