Evaluation of CRISPR gene-editing tools in zebrafish

Background Zebrafish have practical features that make them a useful model for higher-throughput tests of gene function using CRISPR/Cas9 editing to create ‘knockout’ models. In particular, the use of G0 mosaic mutants has potential to increase throughput of functional studies significantly but may suffer from transient effects of introducing Cas9 via microinjection. Further, a large number of computational and empirical tools exist to design CRISPR assays but often produce varied predictions across methods leaving uncertainty in choosing an optimal approach for zebrafish studies. Methods To systematically assess accuracy of tool predictions of on- and off-target gene editing, we subjected zebrafish embryos to CRISPR/Cas9 with 50 different guide RNAs (gRNAs) targeting 14 genes. We also investigate potential confounders of G0-based CRISPR screens by assaying control embryos for spurious mutations and altered gene expression. Results We compared our experimental in vivo editing efficiencies in mosaic G0 embryos with those predicted by eight commonly used gRNA design tools and found large discrepancies between methods. Assessing off-target mutations (predicted in silico and in vitro) found that the majority of tested loci had low in vivo frequencies (< 1%). To characterize if commonly used ‘mock’ CRISPR controls (larvae injected with Cas9 enzyme or mRNA with no gRNA) exhibited spurious molecular features that might exacerbate studies of G0 mosaic CRISPR knockout fish, we generated an RNA-seq dataset of various control larvae at 5 days post fertilization. While we found no evidence of spontaneous somatic mutations of injected larvae, we did identify several hundred differentially-expressed genes with high variability between injection types. Network analyses of shared differentially-expressed genes in the ‘mock’ injected larvae implicated a number of key regulators of common metabolic pathways, and gene-ontology analysis revealed connections with response to wounding and cytoskeleton organization, highlighting a potentially lasting effect from the microinjection process that requires further investigation. Conclusion Overall, our results provide a valuable resource for the zebrafish community for the design and execution of CRISPR/Cas9 experiments. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08238-1.

transparency, small size, and the availability of effective gene-editing tools [3,[5][6][7][8][9][10][11][12]. Continuous improvements of CRISPR editing in zebrafish have allowed efficient targeting of multiple genes simultaneously leading to rapid generation of either mosaic (G 0 ) or stable mutant lines and subsequent characterizations of phenotypes [8,[13][14][15][16][17][18][19]. Such mutants have subsequently been used to test candidate genes associated with human diseases and developmental features [20]. The trend towards more affordable higher-throughput protocols using zebrafish requires a careful evaluation of methods used for the design of CRISPR-based genetic screens and potential confounders that may arise from the microinjection process that could artificially impact phenotypes.
New and creative CRISPR-based approaches in zebrafish address biological questions related to developmental processes (e.g., cell-lineage tracing) as well as gene functions (e.g., epigenome editing and targeted mutagenesis, reviewed in [21]). In the latter application, important factors in generating CRISPR gene knockouts include predicting/maximizing 'on target' Cas9 cleavage activity, predicting/minimizing unintended 'off-target' editing events, and rapidly detecting small insertions or deletions (indels). Presence of indels at candidate loci can be determined in an affordable manner via a number of approaches (reviewed in [22]), ranging from simple identification of heteroduplexes-arising from multiple alleles coexisting in the sampled DNA-visualized using polyacrylamide gel electrophoresis (PAGE) [23] to more sophisticated sequencing approaches that precisely identify and quantify mutant alleles [14,24]. On-target activity of a particular guide RNA (gRNA) can be predicted using tools that provide efficiency scores, often defined by information gathered across empirical assays [25]. One relevant example is CRISPRScan, a predictive-scoring system built from experimental zebrafish gene-editing data based on multiple factors such as nucleotide GC and AT content and nucleosome positioning [9,26]. Bioinformatic tools also exist that define potential regions prone to off-target edits mainly based on sequence similarity and the type/amount of mismatches relative to the on-target region [26]. More recently, several methods have been devised to experimentally identify off-target cleavage sites (reviewed in [27]), including CIRCLE-Seq [28] and GUIDE-seq [29], that do not depend on prior sequence similarity information. These approaches are meant to provide a blind assessment of editing sites but do not necessarily reflect the in vivo activity of on-target activity of the CRISPR/Cas9 complex.
Previous studies have shown CRISPR off-target activity in vivo to be relatively low in zebrafish [8,12,18]. A cross-generational study identified no inflation of transmitted de novo single-nucleotide mutations due to CRISPR-editing using exome sequencing and a stringent bioinformatic pipeline [30] in a similar approach used to identify off-target mutations in mouse trios [31,32]. Other studies have observed off-target mutation rates ranging from 0.07 to 3.17% in zebrafish by sequencing the top three to four predicted off-target regions based on sequence homology [11,12,18]. Although offtarget mutations should not significantly impact studies of stable mutants, since unwanted mutations can be outcrossed out of studied lines relatively easily [14,21], they may be problematic in rapid genetic screens using G 0 mosaics that quickly test gene functions in a single generation.
The increasing number of tools available for the design and execution of CRISPR screens provide an important resource to the zebrafish community. Here, we assayed different available CRISPR on-and off-target prediction methods using empirical data from Cas9-edited zebrafish embryos. We quantified CRISPR cleavage efficiencies in vivo employing a variety of experimental approaches and used these results to compare the accuracy of in silico and in vitro tools for predicting Cas9 on-and/or offtarget activity. Finally, to examine potential confounders that may arise from microinjection of Cas9 into embryos on resulting phenotypes, we assayed G 0 'mock' negative control embryos injected with a buffer containing either Cas9 enzyme or mRNA in the absence of gRNAs by performing RNA-seq and obtained a list of genes with significant differential expression versus uninjected wild-type siblings. In all, these results will serve as a useful resource to the research community as larger-scale G 0 CRISPR screens become more common in assaying gene functions in zebrafish.

Identification of CRISPR-induced indels in zebrafish
We generated a dataset of experimentally confirmed indels within 14 protein-coding genes from injected NHGRI-1 wild-type zebrafish larvae targeted by 50 gRNAs (2-4 different gRNAs/target gene, assembled through the annealing of crRNA:tracrRNA) (Fig. 1A, Supplementary Tables 1 and 2). These 50 gRNAs were designed using CRISPRScan [26] and include a range of predicted editing efficiencies (mean 57.6, range 23-83). To obtain experimental in vivo editing efficiency values for each gRNA, DNA extracted from a pool of 20 G 0 mutant embryos-generated via microinjections of individual gRNAs at the one-cell stage and harvested at 5 days post-fertilization (dpf )-and ~200 bp regions surrounding predicted cut sites for all gRNAs were amplified and Illumina sequenced. To extract the proportion of reads carrying indel alleles, we used Cris-pRVariants [33] with uninjected batch siblings DNA as reference (Supplementary Table 2 for all scores obtained). From this, we inferred an in vivo 'efficiency score' , calculated as the percentage of DNA from injected embryos harboring indels compared to uninjected batch siblings (Fig. 1B).
To compare our efficiency scores with those produced from Sanger-based tools, we also amplified and sequenced ~500 bp fragments surrounding the targeted sites from the same DNA. We extracted the percentage of indels using two different tools that deconvolve major mutations Fig. 1 Workflow for the evaluation of CRISPR cleavages in NHGRI-1 zebrafish embryos. A The cartoon depicts our experiment, which included 50 gRNAs individually microinjected into one-cell stage embryos, DNA extracted from 20 pooled G 0 larvae, and genomic regions targeted by the gRNA amplified. Lightning symbols represent a cleavage event. B An in vivo score was obtained from the Sanger sequencing traces using the ICE and TIDE tools, with an example output from ICE pictured. Scores for the two tools were plotted with values below the median in orange and above the median in purple. C Scores from ICE and TIDE tools were compared to mosaicism percentages from Illumina sequencing of the same regions. D From PAGE, an empirical intensity ratio was obtained and compared to the in vivo scores from Illumina and Sanger sequencing (ICE). Spearman correlation results are shown in the scatter plots with the line of best fit included and their frequencies within Sanger traces-Tracking of Indels by DEcomposition (TIDE) [34] and Inference from CRISPR Edits (ICE) [35] (Fig. 1B). Briefly, these tools use the gRNA sequence to predict the cutting site in the control sequence trace, map the sample trace to this reference, identify indels by deconvolving all base reads at each position, and provide a frequency of the indel spectrum [34,35]. As previously reported [35], both tools provided positively correlated in vivo scores across all gRNAs (Spearman ρ = 0.87, p = 6.78 × 10 − 15 ) with an average score difference of 8.8 ± 12.1 between tools (Fig. 1B). We noted a higher correlation between tools in scores below the median (Spearman ρ = 0.96, p = 1.23 × 10 − 13 ) than above the median (Spearman ρ = 0.65, p = 0.00072; Fig. 1B), suggesting that the deconvolution process of both tools is more accurate when fewer molecules from the pool carry indels. Both ICE and TIDE efficiency scores were correlated with our Illumina-based editing scores (ICE: Spearman ρ = 0.88, p = 9.14 × 10 − 16 ; TIDE: Spearman ρ = 0.59, p = 7.33 × 10 − 6 ; Fig. 1C), though they significantly underestimated editing efficiencies with, for example, Illumina estimates 19.4 ± 16.3 higher than ICE estimated scores (Fig. 1C). Based on its higher correlation, we reported Sanger-based ICE in vivo scores for the rest of this study.
To ascertain consistency of editing efficiencies across embryos, we also repeated microinjections for four gRNAs targeting a single gene (srgap2) and assessed in vivo efficiencies scores of 20 individual larvae using the ICE tool. This resulted in low variance across injections and relative parity of efficiencies versus results from our pooled-larvae DNA preparations (e.g., low efficiency gRNA targeting exon 2; average ± standard error ICE score in individual larvae 18.9 ± 3.3 versus an ICE score of 13 in pooled larvae and a high efficiency gRNA targeting exon 4; average ± standard error in individual larvae 69.2 ± 5.3 ICE score versus an ICE score of 68 in pooled larvae).
A quicker and more affordable approach to quantify CRISPR cleavage efficiency is via PAGE, which takes advantage of the heteroduplexes produced from DNA harboring a mosaic mix of different types of indel mutations [23]. We performed PAGE on ~200 bp regions surrounding the predicted target site for each gRNA and quantified the PCR 'smear' intensity ratio of injected versus uninjected controls (see Methods). These intensity ratios were weakly correlated with our Illumina-(Spearman ρ = 0.37, p = 0.016) and ICE-estimated scores (Spearman ρ = 0.38, Fig. 1D), indicating that accurate quantitative efficiencies cannot be directly deduced from PAGE but that the intensity of PCR 'smear' does qualitatively convey CRISPR-cleavage efficiency.

Accuracy of CRISPR on-target predictions by in silico methods
We next compared the accuracy of CRISPR on-target predictions computed by several published algorithms, including our chosen design tool CRISPRScan [26] (among the most popular tools used by the zebrafish community), CHOPCHOP [36][37][38] using two different scoring methods [39,40] (among the most widely-used tool generally), E-CRISP [41], CRISP-GE [42], CCTop [43], CRIS-PRon [44], DeepSpCas9 [45], as well as the design tool from Integrated DNA Technologies (IDT, www. idtdna. com). Additionally, to assess if strain variability may have impacted our analysis-since all prediction tools used the Tübingen-derived reference genome (GRCz11) [4] whereas our study was performed in the NHGRI-1 strain (a cross between wild-type strains AB and Tübingen [46])-we obtained re-calculated CRISPRScan scores for our gRNAs using a modified zebrafish reference that included known NHGRI-1 variants [46] (now available as an additional reference in the tool browser at www. crisp rscan. org). We then compared all in silico predicted efficiency scores to our in vivo mutagenesis from Illumina sequencing and ICE and observed a generalized underestimation of editing efficiency ( Fig. 2A). Strikingly, only scores predicted by CRIS-PRScan exhibited significant, albeit weak, correlation with in vivo scores. Further, the correlation with Illumina-based scores was significant only when using our NHGRIzed reference (Spearman ρ = 0.31, p = 0.028; Fig. 2B), though CRISPRScan scores were highly concordant with those obtained using the Tübingen-derived reference (Spearman ρ = 0.88, p = 5.02 × 10 − 17 ; Fig. 2B), with an average difference between scores of 4.2 ± 4.6 (range 0-31) (Supplementary Table 2). CHOPCHOP values exhibited correlations with scores from four other in silico tools (E-CRISP, CRIS-PRon, DeepSpCas9, IDT) but none were concordant with our in vivo results (Fig. 2B). Additionally, two tools that utilize deep learning methods (CRISPRon and DeepSpCas9) were significantly correlated with each other but failed to predict in vivo editing efficiencies in our assay (Fig. 2B). Thus, despite the research community broadly adapting all methods for designing gRNAs, there is little consensus in Fig. 2 Correlation of on-target efficiencies calculated using different methods. A Heatmap of the efficiency scores obtained from the design tool (CRISPRScan), in silico prediction tools, and cutting cleavages obtained in vivo using Illumina sequencing and a deconvolution tool from Sanger sequencing [35] for 50 gRNAs. Each box represents a gRNA and the efficiency scores range from 0 (blue) to 100 (red). B Spearman correlations between all efficiency scores from in silico predictions, an in vitro protocol [28], and in vivo cutting assays. Each box includes the correlation result with the p-value in parenthesis. The color of the boxes represent the correlation values, ranging between − 1 (blue) and 1 (red). CHOPCHOP scores were obtained using two different scoring methods, CHOPCHOP (D) (based on [39]) and CHOPCHOP (X) (based on [40]) predicting activity of a particular gRNA among these tools. Based on our results, we recommend using CRISPRScan for choosing gRNAs in zebrafish experiments.

Accuracy of CRISPR on-target predictions by an in vitro method
Next, we evaluated the possibility of using the in vitro protocol CIRCLE-seq [28], an approach designed to identify target sites of a given gRNA by subjecting naked genomic DNA to Cas9 enzyme/gRNA cleavage followed by Illumina sequencing, to obtain an editing efficiency score. It is important to emphasize that such in vitro assays are not designed for predicting on-target editing efficiencies. Nevertheless, we sought to understand if such an approach could be used for this purpose. We tested individually the 50 gRNAs described above using the CIRCLE-seq protocol [47], following the standard recommendations, and computed a log enrichment score normalized by the sequence library size, termed reads per million normalized (RPMN) (see Methods). We found that in vitro-obtained enrichment scores were not correlated with in vivo efficiencies (Illumina: Spearman ρ = 0.10, p = 0.494; ICE: Spearman ρ = − 0.02, p = 0.911, Fig. 2B) or with in silico predictions, with the exception of CHOPCHOP using the scoring method by [40] ( Fig. 2). This indicates that the CIRCLE-seq assay does not necessarily predict on-target CRISPR cleavage activity, at least quantitatively. Previous work from in vivo CRISPR studies of zebrafish suggests that increased GCcontent predicts increased activity of gRNAs [26]. Examining GC content of our tested gRNAs, ranging from 31.8 to 77.3%, we observed a positive correlation with CRISPRScan in silico scores (linear model: beta = 68.18, p = 0.003, adjusted-r 2 = 0.16) and CIRCLE-seq in vitro RPMN scores (linear model: beta = 6.4, p = 0.006, adjusted-r 2 = 0.14) (Supplementary Fig. 1A and B); however, our experimentally determined in vivo scores were not correlated with GC content (Illumina linear model: beta = − 18.57, p = 0.689, adjusted-r 2 = − 0.02; ICE linear model: beta = 12.36, p = 0.817, adjusted-r 2 = − 0.02; Supplementary Fig. 1C), suggesting that additional variables should also be considered (e.g., depletion of A nucleotide bases, nucleosome positioning or DNA accessibility [26,48]).

CRISPR off-target mutation prediction methods
To avoid spurious phenotypes, off-target mutations should be minimized when choosing gRNAs in CRISPR experiments. To characterize off-target mutations for our set of 50 gRNAs, we queried predictions from in silico (CRISPRScan) and in vitro (CIRCLE-seq) methods. CRISPRScan provides a list of predicted off-target sites (between 55 and 1350, median 206.5; Supplementary Table 3) for each gRNA within the zebrafish NHGRIzed reference genome (GRCz11/danRer11) based on a cutting frequency determination (CFD) score that primarily takes into account sequence similarity, location, and type of sequence mismatches [26,38]. The CIRCLE-seq empirical approach also produced variable numbers of sites (between 18 and 874, median 113.5; Supplementary Table 3) per gRNA (defined as 'CIRCLE-seq sites') relative to the control library digested solely with Cas9 enzyme. The number of off-target sites predicted by CRISPRScan exhibited a significant, albeit weak, correlation with the number of CIRCLE-seq sites per gRNA (Spearman ρ = 0.33, p = 0.022, Fig. 3A). Focusing on putatively impactful off-target predictions, an average of 20 ± 13% CRISPRScan-predicted and 64 ± 7% CIRCLEseq sites per gRNA intersected at least one gene (Supplementary Table 3). The sites predicted in silico or in vitro intersecting genes predominantly did not overlap, with an average of 1.6 ± 1.8 (range 0-7) genes per gRNA overlapping between the two approaches for the same gRNA.
To verify if predicted off-target sites were subjected to in vivo Cas9 cleavage, we performed Sanger sequencing of sites within genes identified in silico (n = 17) and in vitro (n = 20) for eight gRNAs, an average of six regions per gRNA (see Supplementary Table 1 for description of sites). Using the ICE tool, we found mosaic mutations at frequencies between 0 and 11%, with 23 out of the 37 sites evidencing indel frequencies below 1% (Fig. 3B), and no differences observed between offtarget sites predicted by CRISPRScan or CIRCLE-seq (Mann-Whitney U = 175.5, p = 0.873; Fig. 3B). To validate the accuracy of ICE at these low indel frequencies, we again performed Illumina sequencing of predicted off-target sites for six of the eight evaluated gRNAs (see Supplementary Table 1 for description) and found significant concordance in results (0.29-7.62% of mosaicism; Spearman ρ = 0.83, p = 0.039; Fig. 3C). The average difference in mosaicism between ICE and Illumina was low (1.6 ± 2.0), with ICE tending to slightly underestimate indel frequencies, highlighting its utility to quickly and economically assess predicted off-targets regions.
We also tested if sites predicted with higher likelihoods of off-target cutting events resulted in higher mutation rates by comparing the indel frequencies among the different levels of prediction (top 1, 2, or 3 prediction scores by CRISPRScan or CIRCLE-seq). No differences were found between prediction groups (Kruskal-Wallis: H (2) = 2.26, p = 0.320; Fig. 3B), suggesting that the information used by the tools to assign probabilities of off-target activity (e.g., CFD scores in CRISPRScan or normalized read counts in CIRCLE-seq) do not necessarily predict the efficiency of cutting at off-target sites in vivo. Thus, off-target cutting mutations at the assessed sites exhibited low frequencies with no clear method performing best. Moreover, none of the on-target scores previously obtained (in silico, in vitro, or in vivo) correlated with the number of predicted off-target sites per gRNA (using either CRISPRScan or CIRCLE-seq), nor the frequency of indels at validated off-target sites (Spearman ρ = 0.27, p = 0.111, Fig. 3D), suggesting that higher on-target efficiencies do not necessarily translate into increased frequencies of spurious off-target mutations.

Evaluating CRISPR Cas9-injection controls
A commonly used 'mock' injection control for phenotypic screens of CRISPR-generated G 0 mosaic lines are embryos injected with buffer and Cas9 in the absence of a gRNA. We sought to determine if such control treatments could significantly impact the genome or transcriptome of our zebrafish larvae. To characterize its impact on genes, we performed RNA-seq of wild-type NHGRI-1 embryos injected with either Cas9 enzyme or Cas9 mRNA (three biological replicates of a pool of five injected larvae), uninjected batch siblings (two biological replicates of a pool of five larvae), and uninjected siblings from another batch (three biological replicates with a pool of five larvae) as controls.

Potential genomic mutations in controls
Recently, Sundaresan and colleagues [49] found that Cas9 in the presence of Mn + 2 ions can result in double-strand cleavage of genomic DNA in the absence of a gRNA. Although their study did not show this same off-target cleavage activity in the presence of Mg + 2 , we hypothesized that aberrant genomic mutations could be incurred by Cas9 due to the presence of MgCl 2 in our injection buffer since Mg + 2 has been shown to compete with Mn + 2 in activating common enzymes [50]. Using our RNA-seq data, we used an optimized pipeline [51] to identify somatic mosaic mutations with uninjected wild-type controls as a reference for common polymorphisms. Focusing only on high-confidence variants (minimum sequence read depth of 20), we filtered already-reported variants in the NHGRI-1 zebrafish line [46], and used the Variant Effect Predictor tool from ENSEMBL to obtain a list of frameshift mutations in protein-coding genes present in our Cas9-injected larvae. A total of 48 and 38 genes were identified with frameshifting variants in larvae injected with Cas9-enzyme and Cas9-mRNA, respectively, with 14 of these genes shared across both injection types (Fig. 4A, Supplementary Table 4 Table 4). Additionally, frameshift variants were positioned closer to a potential Cas9 PAM site (NGG) than by random chance (4 bp median observed distance to closest PAM site; empirical p = 0.0016 using the whole-genome and p = 0.006 using protein-coding regions only, from 10,000 permutations). Therefore, we  (see Supplementary Tables 1 and 7 for the description of the genes) decided to evaluate if indels would consistently arise in these genes in an additional set of microinjections. We performed a new set of microinjections in NHGRI-1 larvae using these same controls (Cas9 enzyme and Cas9 mRNA) and two additional ones commonly used in CRISPR experiments (catalytically dead Cas9 (dCas9) enzyme and a scrambled gRNA coupled with Cas9 enzyme, sequence published in [19]) and evaluated the presence of mutations in 21 genes, including 14 genes with identified frameshift mutations in our RNA-seq data and seven controls with no mutations observed (Fig. 4B, see Supplementary Tables 1 and 5 for the description of all sites). Briefly, genomic DNA was harvested from (1) three pools of five larvae from each group injected at the one-cell stage (Cas9 enzyme, Cas9 mRNA, dCas9, scrambled gRNA); (2) three pools of five uninjected batch siblings larvae; and (3) finclips of the crossing parents as controls. Subsequently, ~200 bp regions surrounding the closest Cas9 PAM site to the previously RNA-seq-identified variants were Illumina sequenced and the alleles extracted using CrispRVariants [33]. We did not observe evidence of inflation of indels in any of the injected groups relative to the uninjected batch siblings or the parental fish, with an overall average mosaicism of 3.1 ± 0.8% per site (below the expected 10% allele ratio for a heterozygous variant in a single individual from a pool of five; Fig. 4B, Supplementary Table 5). Our NHGRI-1 zebrafish carried common single nucleotide variants in the targeted regions, particularly in gene si:ch1073-110a20 where two variants were present in close to 50 and 20% of the reads (Supplementary Fig. 2). Interestingly, we did observe a subtly higher mosaicism in the genes previously detected with variants in our RNA-seq data relative to the regions used as controls (Mann-Whitney U = 2251.5, p = 0.00074, median mosaicism in tested genes 3.4%, median mosaicism in control genes 2.88%; Fig. 4B, Supplementary Table 5). Thus, it is possible that the genes we identified with variants in our RNA-seq data may be naturally prone to carry variants. In summary, these results suggest that currently used CRISPR controls do not suffer systematic DNA cleavages in the absence of a gRNA.

Differential gene expression in controls
We also characterized the impact of injecting Cas9 enzyme or mRNA on the transcriptomes of our zebrafish larvae. Comparisons of transcripts abundances show significant variance across biological replicates when quantifying in both Cas9 treatments, particularly evident in samples injected with the Cas9 enzyme, versus wild-type uninjected larvae (Fig. 5A). This suggests that considerable stochasticity may exist regarding the effects of Cas9 injections in these controls. Examining the genes impacted, we identified hundreds of differentiallyexpressed (DE) genes in our Cas9-injected versus uninjected controls, with a greater number of upregulated genes than downregulated genes (Fig. 5B, Supplementary  Table 6). Specifically, Cas9-enzyme injections resulted in a total of 1100 DE genes (3.6% of the genes assayed), with 756 genes (68.7%) upregulated (fold change > 1) and 344 (31.3%) downregulated (fold change < − 1). Cas9-mRNA injected larvae exhibited 548 DE genes (1.8% of the genes assayed), 376 (68.6%) of these upregulated and 172 (31.4%) downregulated (Fig. 5B). We observed 248 (197 upregulated and 51 downregulated) common DE genes between the two treatments (Fig. 5C), which could be part of a common response to the microinjection process. Network analyses identified commonalities in the shared DE genes enriched in key regulators of different KEGG pathways [52], including spliceosome and ribosome (including genes eif4g2b, eif4g1a, hnrnpd, magoh, hnrnpa0a), hedgehog signaling (shha), glutathione metabolism (gsto2, gsr), GnRH signaling (dusp6), aminoacyl-tRNA biosynthesis (yars), cell cycle (kif2c), glycolysis (aldoca), and cellular senescence (ppp3cca) (Supplementary Fig. 3). Furthermore, while we observed no enrichment in gene ontology terms for downregulated genes, common upregulated genes from both treatments were related to response to wounding (GO:0009611, adjusted p-value = 0.009) and cytoskeleton organization (GO:0045104, adjusted p-value = 0.009) (Supplementary Table 7), revealing molecular consequences of the microinjection process that were still detectable five days later.

Discussion
Our study presents a comprehensive evaluation of empirical and predictive tools currently used for CRISPR editing in zebrafish. Cleavage scores obtained by an in vivo assessment of 50 gRNAs via Sanger sequencing and deconvolution tools (ICE and TIDE) were concordant with Illumina sequencing, the gold standard in predicting efficiencies, as previously reported [35]. Both tools underestimated the presence of non-edited alleles by ~20%, contrary to previous comparisons of TIDE and Illumina sequencing in cell lines, where TIDE showed a ~10-20% overestimation of non-edited alleles [53]. For sites with lower indel frequencies, as we observed for predicted off-target mutations, ICE scores were more concordant with Illumina results (~1-2% difference, again mostly underestimates). Therefore, we suggest that Sanger sequencing deconvolution tools are valuable for establishing relative gRNAs efficiencies but do not necessarily accurately predict absolute cleavage efficiencies in zebrafish in vivo, except at sites with low indel frequencies. In addition, we formalized an empirical 'intensity ratio' score from the commonly-used PAGE approach to assay CRISPR indels and verified its utility in approximating cleavage efficiencies, making it a more affordable and rapid approach to assay editing efficiencies versus sequencing.
On-target efficiency prediction tools showed large differences using the same set of gRNAs sequences, highlighting the importance of understanding features accounted for by each tool. A recent review [25] provides a comprehensive overview of different design tools available and the source of experimental data used to train each one. CRISPRScan [26] was the only tool that could predict on-target efficiency in our set of gRNAs, while no other method provided scores that were correlated with cleavage activities observed in vivo. One limitation of our study was the skew in higher efficiency gRNAs (mean predicted CRISPRScan score of 57.6), which could feasibly impact correlations. Notably, we did obtain more accurate CRISPRScan predictions when we utilized our NHGRIzed reference [46] compared to the current Tübingen-derived reference [4], highlighting the importance of accounting for known genetic variation when designing suitable gRNAs [54,55]. Considering CRISPR-Scan was the only tool that incorporated empirical data from zebrafish, with most methods tested using in vitroderived data, our results emphasize the importance of utilizing a tool trained using in vivo experimental data specific to the study's target species.
An in silico (CRISPRScan) and in vitro (CIRCLEseq) method predicted ~ 20 and 65% potential off-target regions impacting genes, respectively. Notably, we did not evaluate if other predicted sites included cisregulatory elements that could also potentially alter gene expression. Future assessments should include tests targeting a diversity of loci for a more thorough understanding of the potential off-target indels caused by unwanted CRISPR cleavage sites. We observed low off-target mutation frequencies (most < 1%), similar to those previously reported from using single [11,12] or multiple gRNAs [18], although did observe off-target indel frequencies as high as 11% for certain gRNAs. Notably, neither predictive method (CRISPRScan or CIRCLE-seq) nor their likelihood score (using CFD or  Table 6. C Differentially-expressed genes across samples injected with Cas9 enzyme or Cas9 mRNA relative to uninjected batch-siblings show significant correlations. Plots include the numbers and percentages (in parentheses) of genes downregulated (blue) and upregulated (red) in both Cas9 treatments from the total amount of genes assayed (n = 30,258) normalized read count) could accurately predict indel frequencies at off-target sites. Typically, such low mutation frequencies should not be of high impact when generating stable knockout zebrafish lines as these could be easily outcrossed. However, such mutations could have significant impacts on phenotypic outcomes when injected G 0 mosaic populations are analyzed directly.
The adequate selection of controls is a fundamental process in evaluating gene function using G 0 knockout crispant zebrafish, as these larvae serve as baselines from which inferences will be made from. Currently, no consensus exists for preferred controls used in high-throughput CRISPR workflows of zebrafish larvae, which can include targeting a known gene as a positive control (e.g., tyr) [14], uninjected larvae [17,18], sham injections with a Cas9:tracrRNA complex [15], and injections of a scrambled gRNA [16,19], among others. Our RNA-seq assay identified several genes carrying frameshift mutations using uninjected clutch siblings as reference. A follow-up analysis of a second set of injections showed existence of mosaic variants in all injected controls (e.g., Cas9 mRNA, enzyme, and scrambled gRNA), in addition to uninjected siblings and crossed parents at low allelic frequencies (~ 3%). Nevertheless, even though we were limited to our targeted regions, we did observe a higher mosaicism in genes identified as carrying frameshift mutations from our RNA-seq assay compared to control genes, suggesting that these genes could be naturally prone to exhibit mutations in the NHGRI-1 zebrafish line. We also observed high variability in gene expression in larvae solely injected with Cas9 enzyme or mRNA, with several of these DE genes involved in response to wounding processes. Notably, these DE genes were retrieved from 5 dpf larvae suggesting that damage incurred during the microinjection process has a lasting effect. These results suggest that caution should be taken in using G 0 mosaic mutants in investigating phenotypes related to pathways found to be significantly skewed in injection controls, including those involving the function of the spliceosome, ribosomes, and cytoskeleton dynamics.

Conclusions
Overall, we performed a simultaneous assessment of gRNA activities predicted by several commonly used in silico and in vitro methods with those determined experimentally in vivo in injected zebrafish embryos. These results provide valuable information that can be incorporated into the design and execution of CRISPR/Cas9 assays in zebrafish using available workflows [8,13,14,17,18]. Namely, we make the following conclusions and recommendations: • Sanger-based efficiency estimates (TIDE and ICE) tend to underestimate indel mosaicism in zebrafish, though they are more accurate when lower mutational mosaicism exists (such as those observed at off-target sites). • Quantifying heterodimers via PAGE gels represents an affordable method to qualitatively assay CRISPR cutting efficiencies. • Of the existing tools, we recommend CRISPRScan for predicting gRNA on-target efficiency, preferably matched to the zebrafish strain being used. • Off-target mutations occur at relatively low rates with neither in silico nor in vitro prediction methods performing significantly better. • Microinjection of Cas9 (enzymes or mRNA) into embryos does not result in spurious genomic mutations but does impact expression of certain genes and pathways. Caution should be exercised if studying phenotypes related to these genes when performing G 0 mosaic zebrafish screens.
Our aim was to provide information to aid in the decision-making process for future projects using affordable and reliable gene-editing tools in zebrafish. As higherthroughput methods continue to be developed for assaying multiple genes simultaneously, it will be important to use optimal tools for predicting and assessing on-and off-target activity in zebrafish larvae for accurate interpretation of phenotypic outcomes.

Zebrafish husbandry
NHGRI-1 wild-type zebrafish [46] were maintained through standard protocols [56] and their use was approved by the Institutional Animal Care and Use Committee (IACUC) from the Office of Animal Welfare Assurance, University of California, Davis. Animals were kept in a temperature (28 ± 0.5 °C) and light (10 h dark/14 h light cycle) controlled modular system with UV-sterilized filtered water (Aquaneering, San Diego, CA), with a density of 25 adult fish per tank. Feeding and general monitoring of all zebrafish was performed twice a day (9 am and 4 pm). Food included rotifers (Rotigrow Nanno, Reed Mariculture, Campbell, CA), brine shrimp (Artemia Brine Shrimp 90% hatch, Aquaneering, San Diego, CA), and flakes (Zebrafish Select Diet, Aquaneering, San Diego, CA). For all experimental procedures, eggs were collected via natural spawning of randomly selected adult NHGRI-1 zebrafish in 1 l crossing tanks (Aquaneering, San Diego, CA), using a minimum of five breeding pairs (1 male, 1 female) unless otherwise specified. Embryos were grown in standard Petri dishes with E3 media (0.03% Instant Ocean salt in deionized water) and incubated at 28 ± 0.5 °C, using a dissecting microscope (Leica, Buffalo Grove, IL) for developmental staging and daily monitoring until their use for molecular procedures.

Microinjections to generate CRISPR G 0 mosaic mutants
All gRNAs were individually injected into NHGRI-1 embryos to estimate the frequency of indels. gRNAs were prepared following the manufacturer's protocol (IDT). Briefly, 2.5 μl of 100 μM crRNA, 2.5 μl of 100 μM tracr-RNA, and 5 μl of Nuclease-free Duplex Buffer using an annealing program consisting of 5 min at 95 °C, a ramp from 95 °C to 50 °C with a − 0.1 °C/s change, 10 min (min) at 50 °C, and a ramp from 50 °C to 4 °C with a − 1 °C/s change. Ribonucleoprotein injection mix was prepared with 1.30 μl of Cas9 enzyme (20 μM, New England Bio-Labs), 1.60 μl of prepared gRNAs, 2.5 μl of 4x Injection Buffer (containing 0.2% phenol red, 800 mM KCl, 4 mM MgCl 2 , 4 mM TCEP, 120 mM HEPES, pH 7.0), and 4.6 μl of Nuclease-free water. Microinjections directly into the yolks of NHGRI-1 embryos at the one-cell stage were performed as described previously [59], using needles from a micropipette puller (Model P-97, Sutter Instruments) and an air injector (Pneumatic MPPI-2 Pressure Injector). Embryos were collected and ~1 nl of ribonucleoprotein mix was injected per embryo, after previous calibration with a microruler. Twenty injected embryos per Petri dish were grown up to 5 dpf at 28 °C.

Illumina and Sanger amplicon sequencing
DNA extractions were performed on 20 pooled embryos by adding 100 μl of 50 mM NaOH, incubation at 95 °C for 20 min, ramp from 95 °C to 4 °C at a 0.7 °C/s decrease, followed by an addition of 10 μl of 1 M Tris-HCl and a 15 min spin at 4680 rpm. We amplified a ~200 bp region surrounding the targeted site of each gRNA (see Supplementary Table 1 for description of primers). PCR amplifications were performed using 12.5 μl of 2X DreamTaq Green PCR Master Mix (Thermo Fisher), 9.5 μl of Nuclease-Free water, 1 μl of 10 μM primers, and 1 μl extracted DNA. Thermocycler program included 3 min at 95 °C, followed by 35 cycles of 15 s at 95 °C, 30 s at 60 °C, and 20 s at 72 °C, and a final 5 min incubation at 72 °C. Reactions were purified using Ampure XP magnetic beads (Beckman Coulter) and Illumina sequenced (Genewiz, San Diego, CA). To obtain percent mosaicism of mutants by mapping paired-end fastq reads to the zebrafish reference genome (GRCz11/danRer11) using bwa [60] and the R package CrispRVariants [33]. Additionally, we amplified a ~500 bp region surrounding the targeted site of each gRNA from the same extracted DNA for six gRNAs and performed and performed Sanger sequencing (Genewiz, San Diego, CA). Raw trace files were used in the TIDE [34] and ICE [35] tools to predict the percentage of indels, which we used as our in vivo editing score for each gRNA. For both Sanger and Illumina sequencing, we used uninjected batch-sibling embryos as a control reference.

PAGE and intensity-ratio estimation
An empirical cleavage analysis from each gRNA was performed using PAGE. Briefly, we amplified a ~200 bp region in DNA around the targeted site from gRNAinjected and uninjected embryos, as described above. Reactions of the uninjected and injected samples from the same amplicon were run on a 7.5% polyacrylamide gel together for 75 min at 110 V and revealed using Gel-Red (VWR International). Gel images were processed in the software Fiji [61]. For each sample, we defined areas A and B as follows: For each gRNA, the mean-intensity value was obtained for the A and B areas in both the injected and uninjected samples. The A and B areas were exactly the same size between samples. The intensity ratio was calculated as: Log-normalized intensity ratios followed a normal distribution (Shapiro-Wilk test: W = 0.96, p = 0.167) with an average value of 1.21 ± 0.70.

CIRCLE-seq
CIRCLE-seq libraries were prepared for each gRNA (IDT) using genomic DNA extracted from NHGRI-1 (DNA Blood & Tissue kit, Qiagen) following the described protocol [47]. Libraries were sequenced using one HiSeq XTen lane (Novogene, Sacramento, CA), providing an average of 7.3 million reads (range: 4.0 -13.3 million reads) and > Q30 for 92% of reads per gRNA library. Raw reads were processed using the bioinformatic pipeline described [47] (mapping rate > 99% in all samples) to identify regions with cutting events relative to a control sample (treated with Cas9 enzyme and no gRNA). In an attempt to obtain an on-target efficiency estimation from in vitro digestions, we calculated the reads per million normalized (RPMN). For this purpose, we used samtools [61] to extract read coverage from aligned bam files. For each gRNA, coverage was obtained for the third and fourth base upstream of the PAM site as it is the region expected to be cut by Cas9 [62]. RPMN for each gRNA was calculated as the sum of coverage at these two sites divided by the total mapped reads per sample and multiplied by 1 million to scale the values. RPMN scores ranged from 4.42 to 881 (median 99.3) so we decided to use a log normalization to reduce this range.

Variant identification from RNA-seq data
We followed a previously described pipeline to identify somatic variants from RNA-seq data [51]. Briefly, we mapped reads with STAR [63] using the 2-pass mode and a genomic reference created with GRCz11/danRer11 assembly and gtf files (release version 100). Variant calling was performed with MuTect2 as part of GATK [58] using the tumor versus normal mode. 'Normal' was defined by the two uninjected samples to identify all somatic mutations in our Cas9 injected embryos. Variants were annotated using the Variant Effect Predictor tool [64]. High confidence variants (minimum sequencing depth of 20) previously reported for the NHGRI-1 line [46] were removed. Only frameshift loss-of-function variants with a minimum read depth of 20 in canonical protein-coding genes were considered. We extracted the median distance between the identified variants and the nearest Cas9 PAM site (NGG sequence) using the coordinates in the CRISPRScan UCSC track. This median observed distance was compared to the result of median distances of 10,000 permutations of random sampling across the genome and their nearest PAM site. One-tailed empirical p values from this comparison were calculated as (M + N)/(N + 1), where M is the number of iterations with a median distance below the observed value and N is the total number of iterations. We orthogonally investigated the presence of variants in 23 genes via Illumina sequencing of a ~200 bp region surrounding the identified variant location and the R package CrispRVariants [33] (Supplementary Table 1 for primers description). For this purpose, we extracted DNA from three pools of five embryos injected with Cas9 enzyme, Cas9 mRNA, dCas9 (Alt-R S. p. dCas9 protein V3 from IDT), a scrambled gRNA (see Supplementary Table 1 for sequence description), or uninjected. In addition, we extracted DNA from a finclip of the crossing parents of the embryos used for the injections (both female and male). In all of these groups, we quantified the percentage of mutations as all alleles different from the reference.

Statistical analyses
All analyses were performed in R version 4.0.2 [72]. Normality of variables was checked using the Shapiro-Wilk test and parametric or nonparametric comparisons made accordingly. Spearman correlation tests (denoted as ρ) and linear regression models were used to determine the relationship between variables. All analyses compared across different experimental batches included batch as a factor in the model to prevent biases caused by inter-batch differences. Averages include the standard deviation unless otherwise specified. Alpha to determine significance across the different tests was set at 0.05 unless otherwise specified. Additional R packages used for making figures included eulerr [73] and pheatmap [74].

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s12864-021-08238-1.  Table 2. On-target efficiency scores for the 50 gRNAs obtained from different in silico prediction tools, empirical cleavage scores in vivo using Sanger sequencing and sequence deconvolution tools (TIDE and ICE), from the in vitro assay CIRCLE-seq (RPMN), and Illumina sequencing using the CrispRVariants software. Supplementary  Table 3. Description of the off-targets predicted by CRISPRScan and CIRCLE-seq for the 50 evaluated gRNAs. The percentage of sites intersecting genes from the CRISPRScan tool was obtained from the top 30 predicted off-targets since the tool only provides detailed information from these. Supplementary Table 4. List of the frameshift variants identified in protein-coding genes in samples injected with either Cas9 enzyme or Cas9 mRNA using uninjected batch-siblings as reference. Supplementary  Table 5. Illumina mosaicism percentages observed in larvae injected with Cas9 (enzyme or mRNA), catalytically dead Cas9 (dCas9), scrambled gRNA, uninjected batch siblings, and a finclip from their crossing parents. Supplementary Table 6. List of DE genes observed in Cas9-injected samples relative to uninjected batch-siblings and siblings from another batch. Supplementary Table 7. List of gene ontology terms enriched in upregulated genes found in both injection treatments (Cas9 enzyme and Cas9 mRNA).