Frequent loss of heterozygosity in CRISPR-Cas9–edited early human embryos

Segmental Losses and Gains at a CRISPR-Cas9 On-_target Site Identified by Cytogenetics Analysis.

In our previous study (17), in vitro fertilized zygotes donated as surplus to infertility treatment were microinjected with either an sgRNA-Cas9 ribonucleoprotein complex to _target POU5F1 or Cas9 protein alone as a control and cultured for up to 6 d (_targeted and control samples, respectively). We collected a single cell or a cluster of 2–5 cells from these embryos for cytogenetic, genotyping, or transcriptomic analysis (SI Appendix, Fig. S1).

To determine whether CRISPR-Cas9 genome editing leads to complex on-_target DNA damage that would have been missed by our previous _targeted amplicon sequencing, we reanalyzed low-pass WGS data following whole-genome amplification (WGA) from 23 OCT4-_targeted and 8 Cas9 control samples (SI Appendix, Table S1). Given the small sample size, we microinjected additional human embryos with a ribonucleoprotein complex to _target POU5F1, or the Cas9 enzyme as a control, followed by single-cell WGA and low-pass WGS, as before (17). Here and below, the prefix that distinguishes the processing steps is followed by an embryo number and a cell number. The samples used for low-pass WGS were identified with prefix L_ (SI Appendix, Fig. S1). The letter C precedes the embryo number to distinguish CRISPR-Cas9 _targeted from control samples (SI Appendix, Fig. S1). Low-pass WGS data were used to generate copy number profiles for each sample to investigate the presence of abnormalities with a focus on chromosome 6 (Fig. 1A). As an additional comparison, we performed single-cell WGA and low-pass WGS of uninjected control embryos and distinguish these samples with a letter U preceding the embryo number (SI Appendix, Fig. S1)

After preprocessing and quality control, we examined the profiles of 65 samples (25 CRISPR-Cas9 _targeted, 16 Cas9 controls, and 24 uninjected controls; SI Appendix, Fig. S2 A and B). Fifty-six samples exhibited two copies of chromosome 6 with no obvious cytogenetic abnormalities (Fig. 1 C and D and SI Appendix, Figs. S3–S5). Seventeen of the CRISPR-Cas9–_targeted samples, or 68%, had no evidence of abnormalities on chromosome 6. By contrast, we observed that 8 of the 25 _targeted samples had evidence of abnormalities on chromosome 6. Four _targeted samples presented a segmental loss or gain that was directly adjacent to or within the POU5F1 locus on the p-arm of chromosome 6 (Fig. 1 B and D and SI Appendix, Fig. S5). Interestingly, this included two cells from the same embryo where one exhibited a segmental gain and the other a reciprocal loss extending from 6p21.3 to the end of 6p (Fig. 1B). Altogether, segmental abnormalities were detected in 16% of the total number of CRISPR-Cas9–_targeted samples that were evaluated. We also observed that four _targeted samples had evidence of a whole gain of chromosome 6 (Fig. 1 B and D and SI Appendix, Fig. S5), which also represents 16% of the _targeted samples examined. Conversely, a single Cas9 control sample (6.25%) had evidence of a segmental gain on the q-arm of chromosome 6, which was at a site distinct from the POU5F1 locus (SI Appendix, Fig. S4). The uninjected controls did not display any chromosomal abnormalities (Fig. 1D and SI Appendix, Fig. S3).

The number of segmental and whole-chromosome abnormalities observed in the CRISPR-Cas9–_targeted human cells was significantly different from that in the Cas9 (P = 0.0144, two-tailed Fisher’s test) and uninjected control (P = 0.0040, two-tailed Fisher’s test) samples (Fig. 1D). Moreover, this significant difference can be attributed to the observed segmental abnormalities on 6p, because excluding them from the comparison results in a negligible difference in whole-chromosome abnormalities between _targeted and Cas9 control samples (P = 0.1429, two-tailed Fisher’s test). This conclusion is further supported by the fact that none of the _targeted samples show segmental losses or gains on the p-arm of chromosomes 5 and 7, the closest in overall size to chromosome 6, but the frequency of whole chromosome abnormalities is similar to that observed for chromosome 6, suggesting that genome editing does not exacerbate the rates of whole chromosome errors (SI Appendix, Fig. S2C). The comparison we performed between Cas9 control and CRISPR-Cas9 genome edited samples includes a combination of both cleavage and blastocyst stage samples (SI Appendix, Table S1). Because rates of aneuploidy are known to be significantly higher at the cleavage stage compared to the blastocyst (18), we wondered whether excluding the samples at the earlier cleavage stage would alter the conclusions drawn about the rates of aneuploidy in CRISPR-Cas9–_targeted cells. Here, we found that in comparison to uninjected controls there remained a significantly higher proportion of chromosome 6 aneuploidies in OCT4-_targeted cells collected at the blastocyst stage (SI Appendix, Fig. S2D). Altogether, low-pass WGS analysis suggests that a significant proportion of unexpected on-_target events leads to segmental abnormalities following CRISPR-Cas9 genome editing in human preimplantation embryos.

LOH Identified by _targeted Deep Sequencing.

The copy-number profiles described above with low-pass WGS data can only provide a coarse-grained karyotype analysis. To independently investigate the prevalence of LOH events at finer resolution and increased sequencing depth, we designed PCR primer pairs to amplify 15 fragments spanning a ∼20-kb region containing the POU5F1 locus. We also included a control PCR amplification in the ARGFX locus located on chromosome 3 (SI Appendix, Table S4). The PCR amplicons were used to perform deep sequencing by Illumina MiSeq using the gDNA isolated and amplified from 137 single cells or a cluster of 2–5 microdissected cells (111 CRISPR-Cas9 _targeted and 26 Cas9 controls) (SI Appendix, Fig. S1 and Table S2). The prefix W_ distinguished samples whose gDNA was isolated solely for WGA and the prefix G_ was used to demarcate samples that underwent WGA via the genome and transcriptome-sequencing (G&T-seq) protocol (19). All of these samples were different from the samples used for the cytogenetic analyses above.

We then took advantage of the high coverage obtained at each of the sequenced fragments to call single-nucleotide polymorphisms (SNPs), which allowed us to identify samples with putative LOH events: Cases in which heterozygous variants, indicative of contribution from both parental alleles, cannot be confidently called in the amplicons flanking the CRISPR-Cas9 on-_target site directly. Since we do not have the parental genotype from any of the samples that we analyzed, we cannot exclude the possibility that they inherited a homozygous genotype. Therefore, we required the presence of heterozygous SNPs in at least one additional cell from the same embryo to call putative LOH events.

The variant-calling pipeline that we implemented was specifically adjusted for MiSeq data from single cell amplified DNA and includes stringent preprocessing and filtering of the MiSeq reads (Methods). To have sufficient depth of coverage and to construct reliable SNP profiles, we only considered samples with ≥5× coverage in at least two-thirds of the amplicons across the POU5F1 locus (Methods and SI Appendix, Fig. S6A). This threshold allowed us to retain as many samples as possible and still be confident in SNP calling (20). In addition, we implemented a step in our SNP calling pipeline to control for allele overamplification bias, which is a common issue with single cell-amplified DNA (21). This step changes homozygous calls to heterozygous if the fraction of reads supporting the reference allele is above the median value across samples (SI Appendix, Fig. S6 B and C and Methods). Thus, we proceeded with 42 CRISPR-Cas9–_targeted and 10 Cas9 control samples with reliable SNP profiles for subsequent analysis. These data led to the identification of four different patterns: samples without clear evidence of LOH, samples with LOH at the on-_target site, bookended, and open-ended LOH events (Fig. 2A and SI Appendix, Figs. S7–S12).

In samples without LOH (20% of control and 11.9% of _targeted samples), we were able to call heterozygous SNPs in multiple amplified fragments (G_8.04, G_C16.05, and W_C16.05, Fig. 2A). Cases with putative LOH at the locus have heterozygous SNPs in the amplicons covering exons 1 and 5 of the POU5F1 gene (fragments E1-2, G1, and E4 in Fig. 2A) and homozygous SNPs in between (50% of control and 2.4% of _targeted samples). These putative LOH samples would have had to have a cell isolated from the same embryo that had a detectable SNP(s) anywhere in between these flanking exons (e.g., see samples G_8.03 versus G_8.04 in SI Appendix, Fig. S7). Interestingly, this was the most prevalent pattern in Cas9 control samples (Fig. 2B and SI Appendix, Fig. S7), which may indicate the possibility of technical issues due to sequencing or overamplification of one parental allele (see below). Bookended samples have two heterozygous SNPs flanking the cut site but in fragments outside the POU5F1 locus (20% of control and 23.8% of _targeted samples). These LOH events could represent deletions of lengths between ∼7 kb (G_C12.03, SI Appendix, Fig. S10) and ∼12 kb (W_C11.04, SI Appendix, Fig. S9). Finally, in open-ended samples (10% of control and 61.9% of _targeted samples), it was not possible to find heterozygous SNPs in any of the amplified fragments (G_C12.07, Fig. 2A) or there was one or a few heterozygous SNPs on only one side of the region of interest (G_C16.02, SI Appendix, Fig. S12). This was the most common pattern in _targeted samples (Fig. 2B and SI Appendix, Figs. S8–S12) and could represent large deletions of ∼20 kb in length (the size of the region explored) or larger.

As mentioned above, the MiSeq data must be interpreted with caution given the presence of LOH in Cas9 controls. The gDNA employed in these experiments was extracted and amplified with a kit based on multiple displacement amplification (MDA, Methods), which is common in single-cell applications but is known to have high allelic dropout and preferential amplification rates (22). Although, as mentioned above, we implemented a step to control for these biases, this estimate likely undercalls samples with heterozygosity. For example, some homozygous SNPs had 5% of reads mapping to the reference allele but remained homozygous because they fall below the threshold that we used. Considering that we lack the parental genotypes as a reference to choose a more informed cutoff, our method to calculate one from the data represents an unbiased means to correct the presumed allele overamplification in the samples. Moreover, we cannot exclude the possibility that the analyzed single cells inherited a homozygous genotype in the explored region. Nevertheless, the fact that there is a significant number of CRISPR-Cas9–_targeted samples with the largest LOH patterns is notable (Fig. 2B).

Unexpected CRISPR-Cas9–Induced On-_target Events Do Not Lead to Preferential Misexpression of Genes Telomeric to POU5F1.

Our low-pass WGS and SNP analysis above indicate mutations at the POU5F1 locus that are larger than discrete indels. We therefore wondered if this on-_target complexity may encompass the mutations of genes adjacent or telomeric to POU5F1 that could complicate the use of CRISPR-Cas9 to understand gene function in human development or other contexts where the analysis of primary cells is required. To address this, we reanalyzed the single-cell RNA sequencing (scRNA-seq) transcriptome datasets (SI Appendix, Table S6) we generated previously (17) and focused on the chromosome location of transcripts (Fig. 3 A–C). This analysis indicated that differentially expressed genes are not biased to a specific chromosome (Fig. 3A). Moreover, differentially expressed genes are not enriched to either chromosome 6 or the region telomeric to the CRISPR-Cas9 on-_target site (Fig. 3D). These results suggest that the transcriptional differences observed as a consequence of POU5F1 _targeting are not confounded by mutations of genes adjacent, or telomeric, to the on-_target locus. This could be due to a number of reasons. For example, given that the proportion of samples that exhibit unintended CRISPR-Cas9–induced mutations (e.g., segmental aneuploidies or LOH events) is low, the sample size used is sufficiently high to mask any transcriptional differences in genes adjacent to the cut site in samples with segmental loss of the p-arm of chromosome 6. It is also possible that the extent of the on-_target complexity is exaggerated using the gDNA-based pipelines we developed. Notably, because we use single-cell samples, as mentioned above, these are prone to allele overamplification and this can confound the interpretation of on-_target mutation complexity.

No Evidence of On-_target Complexity Using Digital Karyotype and LOH Analysis of the Single-Cell Transcriptome Data.

The use of RNA-sequencing (RNA-seq) data to detect chromosomal abnormalities (23) has great potential to complement the informative low-pass WGS or array CGH methods currently used for embryo screening in the context of assisted reproductive technologies (24, 25). In addition to karyotype analysis, transcriptome data may also provide information about embryo competence at the molecular level. Groff et al. have demonstrated that aneuploidy can be estimated based on significant variations in gene expression in the affected chromosome(s) compared to reference control samples (24). In addition, Weissbein et al. developed a pipeline, called eSNP-Karyotyping, for the detection of LOH in chromosome arms (26). eSNP-Karyotyping is based on measuring the ratio of expressed heterozygous to homozygous SNPs. We applied these two approaches, hereinafter referred to as z-score– and eSNP-Karyotyping, to the scRNA-seq samples (distinguished with the prefix T_) obtained using the G&T-seq protocol (14) (SI Appendix, Table S3). This allowed us to investigate whether transcriptome data could be used to determine the frequency of LOH events in CRISPR-Cas9–_targeted embryos.

Since eSNP-Karyotyping relies on SNP calls from gene expression data, it is very sensitive to depth and breadth of sequencing (26). Therefore, we used results from this method as a reference to select high quality samples for our transcriptome-based analyses (SI Appendix, Fig. S13 A–C). After these filtering steps, we retained 38 samples (22 CRISPR-Cas9 _targeted and 16 Cas9 controls) to analyze further.

In general, we found good agreement between the chromosomal losses detected by z-score–karyotyping and the LOH events identified by eSNP-Karyotyping (SI Appendix, Fig. S14 A and B). For example, the digital karyotype of SI Appendix, Fig. S14A shows the loss of chromosome 4, the p-arm of chromosome 7, and the q-arm of chromosome 14 in sample T_7.01, as well as the loss of chromosome 3 and the p-arm of chromosome 16 in sample T_C16.06. These abnormalities are identified as LOH events in the eSNP-Karyotyping profiles of the same samples (SI Appendix, Fig. S14B). Moreover, the copy number profiles built from low-pass WGS data for different cells from the same embryos also corroborates these chromosomal abnormalities (SI Appendix, Fig. S13 D and E). In terms of events that could be associated with CRISPR-Cas9 on-_target damage, z-score–karyotyping identified the loss of chromosome 6 in sample T_C12.07 (Fig. 4A), which is consistent with the open-ended LOH pattern observed in the gDNA extracted from the same cell G_C12.07 (SI Appendix, Fig. S10) and the segmental loss detected in sample L_C12.01 from the same embryo (Fig. 1B). Also, the gain of the p-arm of chromosome 6 was detected in sample T_C12.15 (Fig. 4A), which is consistent with the segmental gain observed in sample L_C12.02 from the same embryo (Fig. 1B). The gains and losses of chromosome 6 in samples T_2.02, T_2.03, T_2.14, T_7.02, and T_C16.06 (Fig. 4A) are difficult to interpret due to the low quality of their MiSeq data or the lack of amplicon information for the q-arm (SI Appendix, Figs. S7 and S12). Interestingly, eSNP-Karyotyping did not detect any LOH events in chromosome 6 (SI Appendix, Fig. S15), suggesting that this approach is not sensitive enough to detect segmental abnormalities in single-cell samples. Overall, the transcriptome-based karyotypes did not confirm the trends observed in the gDNA-derived data (Fig. 4B).

Ethics Statement.

We reprocessed the DNA and reanalyzed the data generated in our previous study (17). This corresponds to 168 samples (134 OCT4-_targeted and 34 Cas9 controls) across 32 early human embryos (24 OCT4-_targeted and 8 Cas9 controls). For the present work, we used 56 additional single-cell samples (19 OCT4-_targeted, 12 Cas9 controls, and 25 uninjected controls) across 22 early human embryos (1 OCT4-_targeted, 1 Cas9 control, and 20 uninjected controls). This study was approved by the UK Human Fertilization and Embryology Authority (HFEA) (research license no. R0162) and the Health Research Authority’s Research Ethics Committee (Cambridge Central reference no. 19/EE/0297). Our research is compliant with the HFEA code of practice and has undergone inspections by the HFEA since the license was granted. Before giving consent, donors were provided with all of the necessary information about the research project, an opportunity to receive counseling, and the conditions that apply within the license and the HFEA Code of Practice. Specifically, patients signed a consent form authorizing the use of their embryos for research including genetic tests and for the results of these studies to be published in scientific journals. No financial inducements were offered for donation. Patient information sheets and the consent documents provided to patients are publicly available (https://www.crick.ac.uk/research/labs/kathy-niakan/human-embryo-genome-editing-licence). Embryos surplus to the in vitro fertilization treatment of the patient were donated, cryopreserved, and transferred to the Francis Crick Institute, where they were thawed and used in the research project.

CRISPR-Cas9 _targeting of POU5F1.

We analyzed single cells or trophectoderm biopsies from human preimplantation embryos that were CRISPR-Cas9 edited in our previous study (17) plus an additional 56 samples used in the present work. In vitro fertilized zygotes donated as surplus to infertility treatment were microinjected with either a sgRNA–Cas9 ribonucleoprotein complex or with Cas9 protein alone and cultured for 5–6 d (_targeted and control samples, respectively). The sgRNA was designed to _target exon 2 of the POU5F1 gene, and experiments were performed as previously described (17). Genomic DNA from Cas9 control and OCT4-_targeted samples was isolated using the REPLI-g Single Cell Kit (QIAGEN, 150343). DNA samples isolated for cytogenetic analysis were amplified with the SurePlex Kit (Rubicon Genomics). See SI Appendix for more details.

Cytogenetic Analysis.

Low-pass whole genome sequencing (depth of sequencing <0.1×) libraries were prepared using the VeriSeq PGS Kit (Illumina) or the NEB Ultra II FS Kit and sequenced with the MiSeq platform as previously described (17) or with Illumina HiSeq 4000, respectively. Reads were aligned to the human genome hg19 using BWA v0.7.17 (44) and the copy number profiles generated with QDNaseq v1.24.0 (45). See SI Appendix for more details.

PCR Primer Design and Testing.

PCR primer pairs were designed with the Primer3 webtool (https://bioinfo.ut.ee/primer3/, SI Appendix, Table S4). We restricted the product size to 150–500 bp and used the following primer temperature settings: Min = 56, Opt = 58, Max = 60. We tested all primers using 1 μL of genomic DNA from H9 human ES cells in a PCR containing 12.5 μL of Phusion High Fidelity PCR Master Mix (NEB, M0531L), 1.25 μL of 5 μM forward primer, 1.25 μL of 5 μM reverse primer, and 9 μL of nuclease-free water. Thermocycling settings were: 95 °C at 5 min, 35 cycles of 95 °C at 30 s, 58 °C at 30 s, 72 °C at 1 min, and a final extension of 72 °C at 5 min. We confirmed the size of the PCR products by gel electrophoresis. See SI Appendix for more details.

PCR Amplification and _targeted Deep Sequencing.

Isolated DNA was diluted 1:100 in nuclease-free water. We used the QIAgility robot (QIAGEN, 9001531) for master mix preparation (see above) and distribution to 96-well plates (SI Appendix, Table S5). Then, the Biomek FX liquid handling robot (Beckman Coulter, 717013) was used to transfer 1 μL of DNA to the master mix plates and to mix the reagents. The PCR was run with the settings described above. PCR products were cleaned with the Biomek FX robot using the chemagic SEQ Pure20 Kit (PerkinElmer, CMG-458). Clean PCR amplicons from the same DNA sample were pooled to generate 137 libraries that were sequenced by Illumina MiSeq v3. See SI Appendix for more details.

SNP Typing.

We trimmed the MiSeq paired-end reads with DADA2 (46), corrected substitution errors in the trimmed reads with RACER (47), and mapped the corrected reads to the human genome hg38 with BWA v0.7.17 (44). Subsequently, SAM files were converted to the BAM format and postprocessed using Samtools v1.3.1 (48). SNP calling was performed with BCFtools v1.8 (49) using mpileup and call. SNPs supported by less than 10 reads and with mapping quality below 50 were filtered out. To control for allele overamplification, homozygous SNPs were changed to heterozygous if the fraction of reads supporting the reference allele was at least 6% of the total (21). This threshold corresponds to the median of the distribution of the fraction of reads supporting the reference allele across samples. See SI Appendix for more details.

scRNA-Seq Data Analysis.

scRNA-seq reads from G&T-seq samples were processed as previously described (17). Samples with a breadth of sequencing below 0.05 were not considered for any downstream analysis (SI Appendix, Fig. S13 A–C). Differential gene expression analysis was carried out with DESeq2 v1.10.1 (50). For digital karyotyping based on gene expression, we adapted the method described in ref. 24 to identify gains or losses of chromosomal arms (z-score-karyotyping). For digital karyotyping based on SNP expression, we applied the eSNP-Karyotyping pipeline with default parameters (26). See SI Appendix for more details.