Good in vitro activity of sgRNA does not assure in vivo efficacy

Analyzing the efficacy of different sgRNAs in vivo is laborious. To improve the screening for efficient sgRNAs, in vitro digestion of the target sequence can be used. We analyzed the mutagenesis activity of six sgRNAs first in vitro and then selected three for analysis in vivo. As shown in Fig 1, some sgRNAs with good in vitro efficiency presented low or no in vivo activity. This suggests that factors present in vivo prevent Cas9 from acting on its target site. Importantly, cxcr2 had neither detectable gene expression nor mutagenesis efficiency, whereas the genes permissive for mutagenesis (pycard, ca6) showed early expression (Fig 1, S1 Fig). This led us to hypothesize that the onset and the level of gene expression could influence the CRISPR-Cas9 mutagenesis. The corresponding results using the T7 Endonuclease I assay are displayed in S2 Fig. In our hands the T7 Endonuclease I assay has a lower resolution compared to the heteroduplex mobility assay, especially with sgRNAs of lower efficiency. On the other hand, the T7 Endonuclease I assay can be readily used for quantitation of mutagenesis efficiency, especially with sgRNAs of higher efficiency.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. In vitro and in vivo CRISPR-Cas9 mutagenesis efficiencies do not correlate for all genes.

a) An in vitro digestion assay shows that sgRNAs differ in their efficiencies. Below the gene name, + and - indicate the presence or absence of Cas9 protein in the reaction. On the right the wild type (wt) and the mutant products are indicated. b) The in vivo CRISPR-Cas9 mutagenesis visualized for ca6, cxcr2 and pycard with a heteroduplex mobility assay, with the wild type (wt) and the mutant products indicated. 5 embryos were collected per sample at 8hpf.

https://doi.org/10.1371/journal.pone.0196238.g001

The onset of mutagenesis differs between sgRNAs in vivo

As we saw a discrepancy between in vivo and in vitro mutagenesis efficiencies for some sgRNAs, we next analyzed whether the onset of mutagenesis correlates with the onset of gene expression. To avoid the delay of mRNA transcription for Cas9 activity, we used a ready Cas9 protein in our experiments with appropriate preincubation step to allow the sgRNA to complex with Cas9. Three of our functional sgRNAs were chosen for the analysis. The sgRNAs targeting ca10a, sema4gb, or ca6 were co-injected with the Cas9 protein into the 1-cell stage embryo and the onset of mutagenesis was analyzed using both a heteroduplex mobility assay and a T7 Endonuclease I mutation detection assays. As shown in Fig 2 using the heteroduplex mobility assay, the first mutations become detectable as soon as 1hpf for ca10a and sema4gb, whereas the first mutations for ca6 appeared at 3hpf (Fig 2). These results indicate that the onset of mutagenesis differs depending on the sgRNAs in zebrafish embryos. Based on these results, we analyzed the relationship of early gene expression and mutagenesis efficiency in more detail with all our sgRNAs. We were able to detect mutagenesis activity at 1hpf (roughly corresponding to 4-cell stage).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Onset of mutagenesis differs between sgRNAs.

Heteroduplex mobility assay to demonstrate the onset of mutagenesis using high efficiency guide RNAs targeting three different genes with different gene expression patterns in early development. Embryos were collected at timepoints 1, 2, 3, 4, 6hpf (15–20 embryos per group). The gene name above the gel image indicates CRISPR-Cas9 injected embryos and control indicates uninjected controls. The legend on the side indicates the positions of wt (wild type) and mutant bands in the gel. Red arrows indicate the point at which first mutations can be detected.

https://doi.org/10.1371/journal.pone.0196238.g002

Likelihood of successful mutagenesis in relation to the expression level of the target gene in zebrafish embryos

Altogether, we have designed 86 sgRNAs using the crispr.mit.edu, ChopChop (V1 and V2) and CRISPRscan softwares [9,14,15]. Of these sgRNAs, 30% showed detectable in vivo activity (S1 Table). As GC-content (%) has been suggested to influence the effectiveness of CRISPR-Cas9 mutagenesis, we analyzed the GC-content of our sgRNAs (S1 Table) [18]. The GC-content of our functional sgRNAs was found to be similar (Mann-Whitney U-test; p-value 0.452) to that of the non-functional sgRNAs.

When we compared the expression of the genes that we were able to mutate to those we were not, the genes resistant for mutagenesis more often had a very low expression level (Fig 3). However, the difference did not reach statistical significance (Fischer’s exact test; not significant). Moreover, a majority of genes (79%) permissive for mutagenesis underwent an increase in the number of transcripts around the MZT (identified here as a positive change in the number of transcripts between the oblong sphere stage and 50% epiboly). This occurred more often than in the genes resistant to mutagenesis (50%). However, this observation was not statistically significant (Fischer’s exact test) (Fig 3). To examine whether the lack of statistical significance was due to a type two error, we decided to determine whether there is a correlation between target gene expression and mutagenesis efficiency using larger datasets.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Relationship of transcript levels in early development and low mutagenesis efficiency.

Pie charts of the RNA-seq data corresponding to graphs in S1 Fig. a) Number of transcripts for the genes resistant to (left) or permissive (right) for mutagenesis between the oblong sphere and the 15-somite stage (Fischer’s exact test; not significant). 0.5 RPKM (Reads per Kilobase of transcript per Million mapped reads) was used as a limit for low expression. b) The number of genes resistant (left) or permissive (right) for CRISPR-Cas9 mutagenesis in which the number of transcripts is increased or decreased between the oblong sphere-stage and 50% epiboly (around the MZT) (Fischer’s exact test; not significant).

https://doi.org/10.1371/journal.pone.0196238.g003

Mutagenesis efficiency correlates with gene expression and chromatin accessibility in zebrafish embryos

In searching for factors that would explain the poor in vivo activity of some sgRNAs, we investigated available open access datasets. As genes with low expression values tended to be more difficult to mutate in our setting (Fig 3), we analyzed the association between expression levels and mutagenesis efficiency in greater depth, using large datasets in order to avoid type 2 error. We obtained CRISPR-Cas9 sgRNA efficiency data from CRISPRz database for all analyses [33]. We used open access RNA-seq data (E-GEOD-45706) for our primary analysis of the correlation between CRISPR-Cas9 mutagenesis and gene expression [33,34]. We found significant correlations in early development (between 64-cell stage and 36hpf), at the oblong sphere stage (3.66hpf, Spearman correlation 0.227; p-value 0.001) and at 36hpf (Spearman correlation 0.230; p-value 0.001). A strong correlation was observed at the oblong sphere stage which occurs shortly after MBT, around the time of zygotic genome activation. These results suggest that transcriptional activity influences CRISPR-Cas9 mutagenesis at early development (Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Correlation between mutagenesis efficiency and gene expression at different developmental stages during early development.

n = 209.

https://doi.org/10.1371/journal.pone.0196238.t001

As methylation is known to correlate with transcriptional repression, we used zebrafish exon methylation data to analyze whether there is any correlation between exon methylation and the success of CRISPR-Cas9 mutagenesis [33,35]. As is shown in Table 2, there was no significant correlation between exon methylation and CRISPR-Cas9 mutagenesis efficiency at the 1-cell stage or at MBT (Table 2). Similarly, using open access data on embryonic histone methylation, we analyzed whether there is a correlation of H3K4me3 at promoters with CRISPR-Cas9 mutagenesis efficiency. As shown in Table 2, there seemed to be a correlation but this did not reach statistical significance (Spearman correlation 0.263; p-value = 0.074) [33,36].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Correlation between mutagenesis efficiency and chromatin features at different developmental stages and timepoints.

https://doi.org/10.1371/journal.pone.0196238.t002

ATAC-sequencing is a recent next generation sequencing method, which can be used to directly analyze chromatin accessibility. Open access ATAC-seq data for the zebrafish embryo is available at the 4hpf timepoint [37]. We compared mutagenesis efficiency data with ATAC-seq data at transcription start sites for a total of 263 genes. We discovered a significant, albeit rather weak correlation, indicating that chromatin accessibility appears to be one of the factors that explain the efficiency of CRISPR-Cas9 mutagenesis in zebrafish embryos (Table 2).

Zebrafish maintenance

Wild type AB fish were maintained in a flow-through system with a light/dark cycle of 14h/10h according to the standard procedure. Embryos and larvae were grown in an incubator (28.5°C) in embryonic medium/E3 water (5mM NaCl, 0.17mM KCl, 0.33mM CaCl₂, 0.33mM MgSO₄, and 10–15% Methylene Blue).

Ethics statement and data availability

All experiments were carried out in accordance with the EU-directive 2010/ 63/EU on the protection of animals used for scientific purposes, and with the Finnish Act on the Protection of Animals Used for Scientific or Educational Purposes (497/2013) and the Government Decree on the Protection of Animals Used for Scientific or Educational Purposes (564/2013). We have only used zebrafish prior to their independently feeding larval stages in this study, which thus do not require animal permits. Permit for the zebrafish housing and maintenance for the facility at the University of Tampere is ESAVI/10079/04.10.06/2015.

The computational data analysed in this study were collected from open access sources, as detailed in the appropriate sections.

Design and production of sgRNAs for CRISPR/Cas9 mediated genome editing

Target sequences (S1 Table) for sgRNA design were chosen using the online based CRISPR design tool (http://crispr.mit.edu/), ChopChop.V1 or V2 [14,15] or CRISPRscan [9]. Target site uniqueness was verified with the NCBI BLAST analysis against the zebrafish genome (GRCz10). sgRNAs were produced as described previously [46]. Briefly, the sgRNA oligo (Sigma-Aldrich) and the T7 promoter site oligo (S1 and S2 Tables) (Sigma-Aldrich) were annealed and in vitro transcribed using the MEGAshortscript T7 Transcription Kit (Ambion Life Technologies, CA, USA). The integrity and size of the produced sgRNAs were analyzed with gel electrophoresis (1% agarose in Tris-acetate-EDTA, TAE). The concentration of the sgRNAs was measured with the Qubit® RNA BR Assay kit (Thermo Fisher Scientific, MA USA 02451) and Nanodrop 2000 (Thermo Fischer Scientific).

sgRNA and Cas9 microinjection and genomic DNA extraction

The sgRNAs and the Cas9 protein (ToolGen Inc., Seoul, South Korea) were co-injected into one-cell stage zebrafish embryos with a micro injector (PV830 Pneumatic PicoPump, World Precision Instruments) under a Nikon microscope (SMZ645), using borosilicate needles prepared with a Flaming/Brown micropipette puller. Needles were calibrated by injecting solution into a halocarbon oil droplet to achieve a diameter of 12μm (approximately 1nl). The embryos were aligned on 1.2% agarose E3 water plates prior to the injection. An injection solution containing 130ng/μl sgRNA and 250ng/μl of the Cas9 protein in nuclease-free water was incubated 37°C 15min. Rhodamine dextran was added to the solution for the visualization of the injections under a Zeiss Lumar V12 fluorescence microscope. To analyze the onset of the mutagenesis 10–20 CRISPR-Cas9 injected embryos were collected and frozen in liquid nitrogen for DNA extractions at 1, 2, 3, 4, 6hpf (hours post fertilization). To analyze the in vivo mutagenesis efficiency, 5 embryos were collectedat 8hpf and immediately frozen in liquid nitrogen. For DNA extraction, the embryos were lysed 4h 55°C in lysis buffer (10mM Tris pH 8,2, 10mM EDTA, 200mM NaCl, 0.5% SDS, 200μg/ml Proteinase K). DNA was precipitated 1h -20°C using two volumes of ethanol. DNA was then pelleted by centrifuging 16,000g 10min. The pellet was washed with 200μl of 70% ethanol before resuspending in 200μl of water. A purification step with phenol-chloroform was performed after treatment with 15u of RNase A (Thermo Fischer Scientific) per 100μl of sample, 1h 37°C.

Heteroduplex mobility assay

Targeted loci were amplified from the genomic DNA by PCR using the Maxima Hot Start DNA polymerase (Thermo Fischer Scientific) according to the manufacturer’s instructions. The PCR primers (S3 Table) were designed to anneal upstream and downstream of the expected cutting site. The PCR product was purified using Exo I and FastAP (Thermo Fischer Scientific) treatment 15min 37°C, then 15min 85°C. 10μl of the purified PCR product was annealed in a reaction containing 1x NEBuffer 2 (New England Biolabs, MA, USA) and was run on a 10% polyacrylamide gel. The gel was stained with GelRed (Bitium Inc., Fremont, CA).

T7 Endonuclease I mutation detection assay

After purifying and annealing the PCR amplified locus, 10μl of this product was incubated 30min 37°C with 6 units of T7 Endonuclease I (New England Biolabs). The obtained products were separated on a 2.0% agarose TAE gel. The gel was stained with GelRed. The band sizes were compared to control samples.

In vitro digestion of DNA with the Cas9-gRNA complex

To test the in vitro cutting potential, equimolar amounts of the Cas9 protein (ToolGen Inc.) and sgRNA were pre-incubated 15min 37°C in NEB 3 Buffer (New England Biolabs) and 1% Bovine serum albumin (Sigma Aldrich). For the template, a 850–1,200bp site around the target was amplified using Maxima Hot Start DNA polymerase according to the manufacturer’s instructions. The template was then purified (GeneJET PCR Purification kit, Thermo Fischer Scientific). The template was then added to a final 10:10:1 ratio (Cas9:sgRNA:template PCR product). The reaction mix was incubated 3h 28°C as this is the temperature at which zebrafish embryos are maintained. After this, we incubated the sample with 300U of Proteinase K 37°C 10min to release the Cas9. Proteinase K was inactivated by incubation 65°C 10min. Samples were run on a 1% agarose TAE gel to analyze the cutting efficiency.

Gene expression analysis of CRISPR targeted genes

The CRISPRz database contains a list of 1,398 validated zebrafish sgRNAs collected from various published resources [33]. In addition to sgRNA sequences, the associated mutagenesis efficiencies have been recorded in 325 unique zebrafish genes. We compared these mutagenesis efficiencies, from somatic cells, with a publicly available RNA-seq expression dataset housed in the ArrayExpress database [47]. The dataset (ArrayExpress E-GEOD-45706: https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-45706) consists of RNA-seq data performed for samples from multiple stages of zebrafish development: 64-cell, oblong-sphere, 50%-epiboly, 15-somite, 36hpf, 48hpf, 60hpf and 72hpf (and 1 week, excluded from this analysis). Using the Stats package of the SciPy library, we performed Spearman rank correlation analyses of expression data for each sample in each ArrayExpress RNA-seq dataset; using the expression values for genes with available mutagenesis data for somatic cells in CRISPRz [48].

Histone modification in zebrafish promoters

The ArrayExpress dataset E-GEOD-4863 (https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-4863/) is based on custom microarrays for the identification of ChIP binding sites of antibodies against the H3K4me3 in the promoters of zebrafish genes [36]. From the microarray datasets, the log of the median values of the 60-mer probes were summed for each gene, averaged, and then paired with CRISPRz mutagenesis values. Subsequently, these paired values were used to perform the Spearman rank correlation analysis.

Exon methylation analysis of zebrafish genes

McGaughey et al. showed that exon methylation was a better indication of mRNA expression than promoter methylation [35]. Their genome-wide ChIP-seq analysis of whole embryo zebrafish DNA methylation is available as an ArrayExpress dataset E-GEOD-52110 (https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-52110/) at the 1-cell stage and at MBT. We first translated all ChIP-seq peaks from Zv9 genome coordinates to GRCz10 coordinates and then mapped them to exons annotated in the GRCz10 genome. A summation of all ChIP-seq peaks which overlapped exons was calculated for each gene, this sum was divided by the total length of the gene’s exons to generate a methylation coefficient. The methylation coefficients were then combined with mutagenesis data to compute Spearman rank correlations for each timepoint.

ATAC-seq analysis of zebrafish transcriptional units

ATAC-seq is a powerful method for identifying regions of accessible chromatin and it can be used to generate nucleotide resolution mapping of the hyperactive Tn5 transposase binding sites in the genome. An ATAC-seq analysis of 4hpf zebrafish has been previously completed and is available as an ArrayExpress dataset E-GEOD-74231 (https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-74231/) [37]. SRR2747531 was downloaded from the NCBI Sequence Read Archive [37]. Reads were inspected using Fastqc version 0.11.5 and deemed to be of good quality and no further quality filtering or trimming was performed [49]. Subsequently, reads were aligned with Bowtie2 version 2.3.2 [50] using the parameter—very-sensitive-local against the Ensembl Zebrafish reference genome GRCz10. Alignments were filtered and sorted using samtools version 1.4 with the parameter -q 20. Duplicates were removed using Picard Markduplicates version 2.6.0 with the parameters REMOVE_DUPLICATES = TRUE VALIDATION_STRINGENCY = LENIENT [51]. As the alignment was performed against recent reference it was necessary to perform peak- calling independently of the original paper [37]. Furthermore, due to advances in peak calling softwares, peak calling was performed with macs2 version 2.1.1 using the parameters—nomodel—shift -100—extsize 200 -q 0.05 –broad [52]. Transcription start sites (TSSs) for all transcripts annotated in the GRCz10 genome were pooled for each gene. TSSs within 500nt were clustered as a single transcriptional unit. Of the total 22,152 zebrafish genes, 18,687 had a single transcript. Subsequent clustering created single transcriptional units in 2,233 of the remaining 3,465 genes with more than one annotated transcription start site. For clustered TSSs, the midpoint was used as the representative TSS. For each TSS, a +/- 1,000nt region was used to associate ATAC-seq peaks from the E-GEOD-74231 dataset. For each of these regions, an ATAC-seq coefficient was generated by summation of the product of ATAC-seq signal value by total overlap with the TSS region, divided by the length of the region (2,000nt). In cases where after clustering a gene still had more than one TSS ATAC-seq peak, a correlation was performed for all TSS regions and then averaged. Subsequently, all zebrafish genes possessed a single ATAC-seq coefficient. These were then combined with the mutagenesis data from the CRISPRz database in order to compute the Spearman rank correlation.