Targeted gene correction minimally impacts whole-genome mutational load in human disease-specific induced pluripotent stem cell clones

The combination of stem cells and targeted genome editing provides a powerful tool to model human disease as well as to develop curative cellular therapies for genetic disorders. Custom-designed nucleases, including Zinc Finger Nucleases (ZFNs), Transcription Activator-Like Effector Nucleases (TALENs) and Clustered Regularly Interspaced Short Palindromic Repeat (CRISPR)/CAS9 nucleases, specifically induce double-strand breaks (DSBs) in the target genomic loci, which facilitates genome editing by HR (Ding et al., 2013a; Ding et al., 2013b; Hockemeyer et al., 2009; Hou et al., 2013; Lombardo et al., 2007; Soldner et al., 2011; Zou et al., 2009). Several strategies have also been developed to allow for efficient HR without inducing DSBs, which could be genotoxic (Li et al., 2013). Recently, both synthetic nucleases and nuclease-independent methods (HDAdV and Bacteria artificial chromosome) have been independently used for targeted correction of pathogenic mutations of multiple genetic diseases, which results in effective rescue of disease phenotypes (Corti et al., 2012; Fong et al., 2013; Liu et al., 2012; Liu et al., 2011b; Reinhardt et al., 2013; Sanders et al., 2014; Yusa et al., 2011). These and similar studies provide a rationale for applying genome-editing technologies towards developing novel cellular therapies for a variety of debilitating genetic disorders.

An important concern that needs to be addressed before clinical translation of the current targeted gene-correction approaches is the possibility of unwanted genetic variations introduced by the gene targeting procedure. HDAdVs are highly efficient in targeted gene correction in hiPSCs (Li et al., 2011; Liu et al., 2011a; Liu et al., 2012; Liu et al., 2011b). To determine the mutation frequency associated with this method, we first performed deep and pairwise whole-genome sequencing (WGS) to assess DNA sequence variation in disease-specific hiPSCs derived from Hutchinson-Gilford progeria syndrome (HGPS), sickle cell disease (SCD) and Parkinson’s disease (PD) following gene correction by HDAdV (Groups 1–3, Fig. 1A).

We generated on average 60× coverage WGS data (Table S1A). At this depth, greater than 99% (Table S1A) of the bases were sufficiently covered to pass our thresholds for variant calling (Cheng et al., 2012). With stringent criteria to eliminate bias from the sequencing process, we discovered 452, 440 and 665 single nucleotide variants (SNVs) when comparing cHGPS-iPSC, cSCD-iPSC and cPD-iPSC to their corresponding reference lines, respectively (Table 1). We also observed on average 471 loss-of-heterozygosity (LOH) variants per genome in gene corrected cells (Table 1).

Our analysis intially detected a large number of small insertions and deletions (indels) (average 892, range: 892–1052) in all cases. However, Sanger sequencing of indel candidates in exonic regions and additional 35 random loci showed that all of them were false positive. The high false positive rate of indel calling has been reported by others and likely represents a common technical difficulty associated with current WGS technology (Goldstein et al., 2013; Young et al., 2012). Subsequent analysis allowed us to filter out most of the false positive candidates and improve the estimate of the number of real indels (Table 1, see below and Supplemental Experimental Procedures).

The results from the aforementioned cell lines showed a relatively high incidence of sequence variants per sample and big sample-to-sample variation. These observations could be due to the fact that the group 1–3 samples were bulk-passaged and the pairs of lines being compared were not passaged in parallel. In vitro passaging induces mutations and contributes to genetic heterogeneity of hiPSCs (Ji et al., 2012), but heterogeneity in bulk culture iPSCs can be reduced by single-cell subcloning (Sugiura et al., 2014). Our own evaluations have also indicated greater overall variation relative to an hg19 reference genome in bulk passaged lines than in single cell subclones (data not shown). We therefore evaluated whether our analysis could be further improved by: (1) maintaining the reference cell lines in parallel with cells undergoing gene correction; (2) using single-cell derived clones to minimize pre-existing genetic heterogeneity; (3) including additional gene corrected clones generated using alternative methods to evaluate the contribution to the mutational load by these factors.

To this end we performed gene-correction experiments de novo in the SCD context using nuclease-independent (HDAdV) and nuclease-based (TALEN) methods. All samples in this group (group 4) started from a single-cell derived and integration-free SCD-iPSC clone (Li et al., 2011), which was also maintained without genetic manipulation in parallel with gene-targeted lines as an unmodified control (cHBB-ctrl). We designed two independent TALEN expression vectors that target exon 2 (TALEN1) and a downstream region of the HBB locus (TALEN2) (Fig. S1A). The efficacy of the TALENs was confirmed by Surveyor assay in HEK293 cells (Fig. S1B). For TALEN-mediated gene correction, an HR donor vector was also used (Fig. S1C and D). Additionally, we performed gene targeting using two independent HDAdV-based gene-correction vectors (HDAdV1 and HDAdV2), one of which was described previously (Li et al., 2011) (Fig. S1E and F). The TALENs- and HDAdVs mediated gene-correction with comparable efficiencies per integration at the HBB locus (Fig. 1B). Validated gene-corrected clones (cHBB-TALEN1, cHBB-TALEN2, cHBB-HDAdV1 and cHBB-HDAdV2) were subjected to deep WGS (coverage: 56.3–66.4x, 99.15–99.20% of the bases, Fig. 1A, Fig. S1G and H, and Table S1A). Since group 4 TALEN- and HDAdV-mediated gene corrected clones and unmodified clone showed substantially fewer SNVs (range 74–129) compared to groups 1–3 (Table 1) and less variation in SNV numbers between samples, we focused on group 4 for subsequent analysis.

To independently validate the SNVs, we obtained high quality Sanger sequencing data for the 105 SNVs in unmodified cells. Eighty-two percent of the sites contained the annotated mutation (Table 1). We calculated the number of predicted SNVs in all of the analyzed clones in group 4 (range 67–118) as well as in group 1–3 (range 363–549). Importantly, regardless of targeting vector design and methodology, there was no dramatic increase of SNVs in gene corrected clones compared with the control line (Table 1). Additionally, HDAdV and TALEN approaches resulted in comparable numbers of SNVs (Table 1). On average, ~83 SNVs were introduced per genome following gene correction. This number is significantly smaller than the number of SNVs (over 1,000) reported to accumulate during somatic cell reprogramming (Cheng et al., 2012). There was no specific chromosomal region or functional pathway that showed significant enrichment and none of the variants were common to all samples (data not shown).

The majority of the samples in group 4 did not contain any exonic SNVs. Two gene corrected clones – cSCD-HDAdV1 and cSCD-TALEN2 – contained 3 and 1 validated exonic SNVs, respectively (Table S1B). Among them, two were nonsynonymous changes. Based on Polyphen2 (Adzhubei et al., 2010), the missense mutation in the CBX2 gene is predicted to be benign, while the one in the AMBP gene is potentially damaging. However, given the heterozygous nature of these SNVs, the potential harm caused by the mutations could be mitigated by the presence of a wild-type allele.

We next sought to understand the functional relevance of the non-coding SNVs found in our samples (Table S1C). We compared all 489 non-coding SNVs in group 4 with a catalog of 12,345 SNPs identified in 1,791 GWAS studies published to date and found zero overlap. Similarly, none of the SNVs was found in the catalogue of somatic mutations in cancer (Forbes et al., 2008). A low numbers of SNVs were found to overlap with potential regulatory elements such as promoter CpG islands, ENCODE DNaseI hypersensitive sites and conserved regions in vertebrates. Together, although the functional impact of non-coding SNVs is difficult to ascertain, these results indicate that targeted gene correction by HDAdVs or TALENs did not introduce more SNVs above the background level of the unmodified control, and that it should be possible to avoid clones that bear harmful mutations for downstream applications, especially considering the heterozygous nature of the SNVs.

LOH has been commonly reported in genomic sequencing of tumor samples (Ha et al., 2012). We noted big variations in LOH number in group 1–3, which might be attributed to their distinct genetic/disease background. In contrast, LOH number in group 4 was consistently and significant lower than that of group 2, which shared the same genetic background as group 4 (Table 1). All of the LOH variants observed in protein-coding genes matched one of the parental chromosome DNA sequences, indicating that at least these LOH variants are not harmful to gene encoding proteins.

To estimate the number of true indels in the genome, we performed Sanger sequencing for all of the 327 indels detected in SCD-ctrl. Out of 208 successful assays, two (0.96%) showed concordant indels as predicted. Based on this analysis, we generated a new filter that removed most false positive indel calls and therefore enriched the true indels by about 10-fold (see Supplemental Experimental Procedures). We then performed Sanger sequencing for all of the indel candidates, and concluded that the occurrence of novel indels in gene-corrected clones was low (1–4 per genome) and not above the background level observed in the unmodified control line (Table 1 and Table S1D). Additionally, we compared all indels and LOHs with known disease-associated SNPs and found no overlap (data not shown).

We identified 40 SNV candidates that appeared in at least two of the eight lines of groups 1–4. After PCR and Sanger sequencing, no common SNV was found among all of the gene-corrected and uncorrected clones. Similarly, there were no common indels among groups 1–4. This result suggests that mutations are induced in a random manner.

To understand the extent of large insertions and deletions that may not be detected by the indel-calling algorithm, we performed structure variations (SVs) analysis by CREST (Wang et al., 2011). Six SVs ranging from 187 bp to 1,581 kbp were detected and further validated by PCR/Sanger sequencing. Importantly, five of 7 gene-corrected clones did not show any structure variance (Table S1E). We also found that there is only one SV in the case of TALEN1, and that it is far from any predicted off-target sites, suggesting no relationship between SVs and TALEN off-target sites (Table S1D and S1E, and see below).

We looked for evidence of TALEN-induced off-target mutations near potential off-target sites (20,542 for TALEN1 and 29,026 for TALEN2, see Supplemental Experimental Procedures) predicted by TALENoffer in cHBB-TALEN1 and cHBB-TALEN2 iPSCs. Of note, the performance of TALENoffer has been validated recently (Guilinger et al., 2014; Liu et al., 2014; Miyaoka et al., 2014). However, we found no incidence of the typical >2bp indels near predicted spacer regions in either TALEN corrected line (Table S1D). DSB-induced recombination in eukaryotes is also known to result in a 100- to 1000-fold increase in the mutation rate of genes adjacent to the DSB (Holbeck and Strathern, 1997; Rattray et al., 2002; Strathern et al., 1995). Among all of the SNVs and indels, three SNVs were found within 1,000 base pairs of a predicted off-target site only in cSCD-TALEN2 (Table S1D). Although it is possible that these SNVs were caused by a TALEN-mediated DSB associated point mutation, further evidence is necessary to support this hypothesis.

Given that WGS analysis showed that both HDAdV and TALEN approaches did not significantly elevate mutational burden in gene corrected hiPSCs, we hypothesized that the combination of these two methods could result in a more robust gene-editing tool without sacrificing genome integrity. Accordingly, we created a novel HDAdV vector (HBB-c-telHDAdV) harboring both HBB-specific long homology arms (13.3 kb) and TALEN expression cassettes (Fig. S1I and J). The use of this hybrid vector (hereafter called telHDAdV) resulted in significantly higher genome-editing efficiency when compared with either single vector (up to 12.5 fold if calculated per drug resistant colony or up to 30 fold if calculated per cell, Fig. 1B and Fig. S1K). Recently, CRISPR/CAS9 nucleases has been shown to outperform TALENs in editing the same loci in human pluripotent stem cells (Ding et al., 2013b). To compare the performance of CRISPR/CAS9 nucleases and telHDAdVs, we designed two CRISPR/CAS9 nucleases targeting the same loci as the TALENs encoded by the telHDAdVs (Fig. S1A). The gene targeting efficiency of these CRISPR/CAS9 nucleases was on par with the corresponding TALENs at the HBB locus (Fig. 1B and Fig. S1K). Interestingly, a recent report also showed that although CRISPR/CAS9 nucleases are more efficient in NHEJ mutagenesis than TALENs at the HBB locus, these methods are comparable in terms of homology directed repair (Hendel et al., 2014). Importantly, telHDAdVs showed clear superiority over both CRISPR/CAS9 nucleases and TALENs in gene targeting and gene correction of the HBB locus (Fig. 1B and Fig. S1K). These results show that the combination of the nuclease-dependent and –independent techniques may constitute a novel and productive strategy to optimize genome-editing approaches.

Finally, we performed Sanger sequencing to determine if the hybrid vector-corrected lines share any common nucleotide variations with the lines corrected with the corresponding TALENs or HDAdVs (Table S1D). No common SNV/indel was found between the HDAdV- and telHDAdV-corrected lines among 51 validated variants. Two SNVs out of 35 validated SNVs /indels in the TALEN1-corrected line were found in tel-HDAdV1-corrected cells. Similarly, one of the 32 validated SNVs /indels was shared between TALEN2- and tel-HDVAdV2-corrected lines. Because these common SNVs between TALEN- and telHDAdV-corrected clones were not found elsewhere, they are likely to be due to off-target effects of the TALENs. However, all of the common mutations are far away from any predicted TALEN off-target sites (range 17 kb – 249 kb; Table S1D). It is possible, as it has been previously reported with ZFN (Gabriel et al., 2011; Pattanayak et al., 2011), that these mutations are near cryptic off-target sites that are not predicted in silico. Thus, while we did not find mutations specifically associated with HDAdV, our results suggest that TALEN might induce rare SNVs in a non-random manner.

In summary, our study provides evidence that HDAdV- or TALEN-mediated genome editing results in low levels of sequence variants in single-cell clones, which are not significantly different from what is observed during routine hiPSC culture. Our WGS analysis of TALEN-edited hiPSCs shows no clear evidence of widespread off-target modification. Interestingly, analysis of the telHDAdV-corrected clones revealed rare SNVs that were also in TALEN-corrected clones and were far from any predicted off-target sites, suggesting a need for improvement in off-target site prediction. Our findings only cover TALENs of a 20+20 configuration, and the specificity of TALENs is sensitive to the length of the TALE recognition sites (Guilinger et al., 2014; Liu et al., 2014; Miyaoka et al., 2014). Additionally, our analysis of off-target effects is broad but limited in depth, so it is very likely and in fact expected that other low-frequency off-target effects are not detected due to the small sample number.

In sum, therefore, it seems reasonable to conclude that the mutational load attributable to genome editing technologies in single-cell clones can be minimal and should not constitute an inherent obstacle for future application of single-cell clone-based gene editing approaches in human regenerative medicine. In addition, the development of combined strategies such as the hybrid gene-editing vector technology presented here may further accelerate clinical translation.