Assessment of genomic changes in a CRISPR/Cas9 Phaeodactylum tricornutum mutant through whole genome resequencing

Target sequence design

Putative target sequences found in the locus chr9: 533409–537647 of the P. tricornutum genome were BLASTed against the same genome with the following criteria on the NCBI (https://www.ncbi.nlm.nih.gov/) website: RefSeq at Representative Genome Database, optimized for somewhat similar sequences, as general parameters: expect threshold 10,000, word size 7, as scoring parameters: match/mismatch 1, −3, without filters and without mask. The sequence showing the lowest identity percentage within the rest of the genome was chosen with particular attention at the 3′ end (seed region). The oligonucleotide pair chosen was sgRNAfw 5′-tcgaatactatgcttggatggcgg-3′ and sgRNArv 5′-aaacccgccatccaagcatagtat-3′ (GC% = 50).

Preparation of the constructs

The target oligonucleotides (two μg each in 50 μl, 40 ng/μl) were annealed 10′ at 100 °C, cooled at RT, then diluted 20×. PtPuc3_diaCas9_sgRNA (https://www.addgene.org/109219/) has been digested with BsaI and 50 ng of the linearized construct were ligated by T4 (2 h at RT) with four ng of insert. The ligation reaction was microdialized and two μl of it were electroporated in one shot top 10 cells (Invitrogen, Carlsbad, CA, USA). The clones were verified by digestion with BsaI and HindIII (in positive clones the BsaI restriction site is lost) and sequenced with pdiaCasfw 5′-CGGTGAGCTGGAAATTGGTT-3′and pdiaCasrv 5′-CCAAGACATGCTACCCGCA-3′ to verify the presence of the target sequence. The obtained plasmid pPtPuc3_diaCas9_sgRNA containing the target sequence (cargo vector) was chemically transformed in DH10B cells containing the conjugation vector pTA-Mob (Strand et al., 2014). Aliquots of transformed DH10B containing conjugation and cargo vectors were cryopreserved and inoculated the day before the conjugation for overnight growth.

P. tricornutum conjugation with E. coli

Axenic CCMP632 strain of P. tricornutum Bohlin was obtained from the Provasoli-Guillard National Center for Culture of Marine Phytoplankton. Culture were grown in f/2 medium (Guillard, 1975) at 18 °C under white fluorescent lights (70 μmol m⁻² s⁻¹.), 12 h:12 h dark–light cycle. The DH10B containing pTA-Mob and pPtPuc3_diaCas9_sgRNA were used for the conjugation of P. tricornutum performed as in Karas et al. (2015) using f/2 as diatom medium as the only modification.

Selection of resistant clones

Resistant clones were selected on ½f/2, 1% agarose, phleomycin 20 μg/ml containing plates. Colony PCR was performed as in De Riso et al. (2009) and PCR was performed on cell lysates with ShBlefw (5′-acgacgtgaccctgttcatc-3′)-ShBlerv (5′-gtcggtcagtcctgctcct-3′) and CEN6ARS4fw (5′-gaccacacacgaaaatcctg-3′)-HIS3rv (5′-gctctggaaagtgcctcatc-3′) primer pairs to ascertain the presence of the episome (pPtPuc3_diaCas9_sgRNA) in the cells.

Sanger sequence analysis of the resistant clones

Resistant clones, resulting positive to the colony PCR amplification, were amplified by PCR with the locusfw (5′-atcctgtcgaaatggcaaaa-3′)-locusrv (5′-ttgtatcgctctacggcttg-3′) primer pair to produce and sequence the 403 bp amplicon encompassing the target locus.

RT-PCR

Total RNA was extracted from P. tricornutum axenic cultures of wild type and mutant strains and cDNA retrotranscription was performed as in Russo et al. (2015). PCR was performed on cDNA with locusfw-locusrv primer pair to produce and sequence the 403 bp amplicon encompassing the target locus. PCR amplification was performed with the locusfw2 (5′-attggagaccattcgtgaag-3′)-locusrv (5′-ttgtatcgctctacggcttg-3′) primer pair to produce and sequence the 984 bp amplicon encompassing the target locus and containing single nucleotide polymorphisms (SNPs).

P. tricornutum genomic DNA preparation and WGS

Genomic DNA was extracted from 400 ml wild type and mutant P. tricornutum axenic cultures in exponential growth phase by centrifugation, followed by freezing the cells in liquid nitrogen. The frozen pellet was then resuspended in two ml of lysis buffer, incubated at 37 °C for 15′ and centrifuged at 10,000×g at 4 °C for 10′. An equal volume of a mix of phenol:chloroform:isoamyl alcohol (25:24:1) was added to the supernatant and centrifuged for 15′. An equal volume of chloroform:isoamyl alcohol (24:1) was added to the supernatant and centrifuged in the same conditions. Two volumes of NaAc 3M pH 5.4 and ethanol 100% were added to the supernatant and precipitated at −20 °C. Following a centrifugation of 30′, five ml of ethanol 70% was added to the pellet twice to remove salts. The resulting pellet was resuspended in water and then treated with RNAse.

Genomic DNA libraries were prepared using KAPA Hyper Prep kit (PCR-free) following manufacturer’s standard protocol. Libraries have been sequenced by Illumina HiSeq 4000 as paired-end (2 × 150 bp) following manufacturer’s standard protocol. Data have been deposited in NCBI with the accession number PRJNA453101.

WGS analyses

Raw sequencing reads were processed with FastQC (v0.11.3, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and BBDuk (version 35.82, https://jgi.doe.gov/data-and-tools/bbtools/) to perform a quality check and to remove Illumina adapters and bases with a Phred-Like score less than 30. At the end of the quality check all the reads shorter than 35 bp were removed. High-quality reads were then mapped against the P. tricornutum reference genome (Ensembl Protists accession number ASM15095v2) with BWA MEM (version 0.7.12-r1039) (Li & Durbin, 2009) with default options. The obtained SAM files were then converted to BAM with samtools (v1.2) (Li et al., 2009) and then picard-tools (v1.131, https://broadinstitute.github.io/picard/index.html) was used to add read group and to remove optical duplicates. SNPs and small insertions or deletions (indels) were detected by processing the obtained BAM files with FreeBayes (v1.1.0-3-g961e5f3) (Garrison & Marth, 2012) with the following options: -m 30 -q 20 --min-coverage 5 --max-coverage 200 --genotype-qualities. The obtained VCF was filtered to keep only the variants with a quality score higher than 30 and with coverage in both the strains. Bigger structural variations were identified with the software Breakdancer (v1.4.5) (Chen et al., 2009). A manual curation was performed on the detected variants to focus only on those with a clearer difference between wild type and mutant. The heterozygosity rate across the genome was calculated as follows: the genome was divided in 50 kb windows with bedtools (Quinlan & Hall, 2010), then the windows were intersected with the Freebayes filtered VCF and the number of homozygous (1/1 or 0/0) and heterozygous (0/1) variants were calculated. Plots were generated with ggplot2 (Wickham, 2009) applying the loess smooth algorithm. The same procedure was followed to process public available Pt1 strain reads (SRX3566562).

Potential off-targets were predicted with the online Cas-OFFinder tool (http://www.rgenome.net/cas-offinder/new) (Bae, Park & Kim, 2014).

Design of the target sequence

The coding sequence corresponding to the Phatr3_J46193 gene in P. tricornutum chr9: 533409–537647 locus, was scanned for the presence of NGG PAM sequences. Among the putative target sequences, the ones that occurred just upstream of the active sites of the protein were chosen. These sequences were submitted to a nucleotide BLAST analysis on the NCBI website against the P. tricornutum genome, and target sequences showing higher identity with sequences in other regions of the P. tricornutum genome were excluded. Particular attention, in the identity analysis, was paid to the 3′-end, because it has been proposed that the 8–12 PAM-proximal bases, known as the seed sequence, determine targeting specificity of the Cas9 protein (Nishimasu et al., 2014). No software for off-target prediction in diatoms were available at the beginning of this study. A target sequence with a GC content of 50% was chosen considering that it has been reported that a GC content higher than 70% may increase the probability of off-target effects (Tsai et al., 2015).

Selection of mutant clones

The transformation efficiency was 2.0 × 10⁻⁵. One month after transformation we performed colony PCR on 14 clones to verify the presence of the episome in the cells, checking for the presence of two different regions of the episome, a fragment of the bleomycin antibiotic coding sequence and a fragment of the yeast centromeric CEN6-ARSH4-HIS3 sequence that enables episome maintenance in P. tricornutum (Karas et al., 2015) (Fig. S1). All the analyzed clones were successfully transformed. On a first group of 20 clones, resistant to the antibiotic selection, we performed colony PCR and sequenced 403 bp, encompassing the target sequence, obtaining two mutant clones. These displayed overlapping traces starting at the mutation site, an indication that two different sequences are present on the two alleles. After 1 month, we analyzed 20 additional clones with the same procedure and we found 16 mutated clones showing overlapping traces in chromatograms, observing an increase in the percentage of the mutated clones. The overlapping traces obtained in sequencing resulted from the presence of a mixture of PCR products that can be due to: (i) the occurrence of a monoallelic mutation (one allele is mutated and the other one is wild type), (ii) biallelic heterozygous mutations (both the alleles are mutated in different ways) or to (iii) a phenomenon of colony mosaicism, that is, a colony consists of a combination of cells with wild type and mutant alleles due to mutations occurring after transformed cells have started dividing (Daboussi et al., 2014).

We chose two clones and subcloned the PCR products in TOPO vector, obtaining, respectively, for clone#1: wild type, del_8, del_2, ins_7; for clone#2: del_13.1, del_23, del_13.2, del_33 (Fig. S2). We decided to focus on clone#2 by manually isolating 24 cells to obtain monoclonal cultures deriving from a single cell. From these 24 clones we obtained: one wild type, two clones with a del_1 with single peaks in the sequence (Fig. 1A), while the other 21 clones showed overlapping peaks indicating the occurrence of two different deletions on the two alleles or the presence of a monoallelic deletion together with a wild type allele (data not shown). We sequenced the same amplicon from the cDNA of the #2del_1 mutant clone and found the same mutation, confirming the presence of the single nucleotide deletion and excluding any repair events at the mRNA level (Fig. 1A).

Isolation of a single clone with specific biallelic mutations

We decided to focus on the mutant clone#2 del_1 with the clear deletion of a single nucleotide. The single nucleotide deletion affected a guanine at position Chr9_535251, corresponding to position -4 immediately adjacent to the PAM sequence. This change in the transcript generated a frame-shift mutation causing the creation of a premature stop codon, and resulted in a truncated protein product lacking the active sites. The mutant strain was subcloned on plate and all of the resulting colonies displayed the same mutation, reinforcing the evidence that they derived from a single clone. Since all the sequences showed single peaks, we imagined two scenarios: the one nucleotide deletion was present on both the alleles (biallelic and homozygous mutation) or only this small mutation was visible by PCR and Sanger analysis and a deletion larger than the obtained amplicon was present (biallelic and heterozygous mutation). To distinguish between the two possibilities of a homozygous or a heterozygous mutation, we performed a PCR amplification on a larger cDNA fragment of 984 bp encompassing the locus of interest and containing SNPs. We obtained fragments of identical length from the wild type and from the mutant strain. Since SNPs were observed in the wild type sequence but not in the mutant (Fig. 1B), the presence of a large deletion on one of the alleles in the mutant could be hypothesized.

Exclusion of Cas9 activity

Once we had identified the desirable mutation in the target locus, we eliminated the selective pressure of the antibiotic in the medium and repeated the colony PCR, verifying that the episome had been lost in the cells grown for 3 or 5 months in the absence of the antibiotic, while the episome was still present in the cells grown for 5 months under antibiotic pressure (Fig. S3).

Analysis of off-target mutations in the genome by WGS

A wild type strain and a mutant strain in which the episome had been removed by relieving the selective pressure were sequenced by WGS. About 21 and 14 million of paired-end 150 nt reads were obtained from the wild type and the mutant strain, respectively. After trimming and duplicate removal, 10,237,553 and 16,194,883 of properly paired reads were mapped against the P. tricornutum genome from wild type and mutant, respectively, corresponding to a predicted coverage of 44× and 70×.

Mutant and wild type strains were compared to the reference genome and between themselves in order to identify both small variants (SNPs, indels) and larger structural variations, such as large insertions and deletions. Overall, 279,279 (81.31% SNPs) and 274,730 (81.32% SNPs) small variants were identified in mutant and wild type respect to the reference genome, respectively. Of those, 265,209 variants were identified both in the mutant and in the wild type, whereas 12,390 and 7,843 were found only in mutant and in the wild type, respectively. A total of 713 larger structural variations were identified comparing mutant and wild type with the reference genome, with a support of at least 10 reads. Those variants included 303 inter-chromosomal translocations, 294 deletions, 37 inversions and 79 intra-chromosomal translocations, with a mean variant size of 14,297 bp. When comparing wild type and mutant, we found seven variants specific for the mutant strain (one intra-chromosomal translocation and six deletions) and one deletion specific for the wild type (Table S1).

The results obtained from the variant calling successfully detected the mutant specific deletion of the G nucleotide (Chr9:535251) on one allele of the target Phatr3_J46193 locus (Figs. 2A and 2B). On the other allele in the same locus we detected a mutant-specific 3,639 nucleotides deletion (Chr9:535239–538878) (Fig. 2C) which was supported by the coverage reduction, the distance of the paired-ends and the presence of split-reads on the breaking point. This result confirmed our hypothesis of a large deletion on one of the alleles in the mutant (Fig. 1B).

By using the Cas-OFFinder tool we predicted 98 possible off-targets of Cas9, distributed on 15 chromosomes (Table S2). Interestingly none of the observed structural variations nor any specific SNP or indel were found in any of the potential off-targets in the mutant strain.

On the other hand, we found long stretches of loss of heterozygosity (LOH) in the mutant clone respect to the wild type but also vice versa (Figs. 3A and 3B). More specifically, LOHs were observed for chromosomes 3, 6, 15, 21, 24 and 25 in the mutant, and for chromosomes 5, 10, 19 and 22 in the wild type. LOHs were very broad in size, ranging from a minimum of 20 kb (i.e. LOH on chromosome 15) to nearly an entire chromosome (LOH on chromosome 25). Some of the observed LOH included regions predicted as putative off-targets such as those on chromosomes 3, 6 and 25 (Table 1), and might be due to a break in one of the alleles successively repaired by homology-directed repair. However, no off-targets were predicted on the other chromosomes affected by LOH. It is important to point out that when LOH was observed this was not due to a deletion on one of the two chromosomes since the number of reads in that region remained as high as in the other regions (when a deletion occurs the number of mapped reads to that region is halved), therefore all the observed LOHs might be due to gene conversion.

In order to get more insights into the observed variation in heterozygosity, a public dataset of WGS of the Pt1 strain, hereafter named “public Pt1” (Rastogi et al., 2017) was processed with the same bioinformatics pipeline. The results in Figs. 3A and 3B highlight that the three strains show very similar profiles of heterozygosity across the genome and the only variations affect either the mutant or the wild type strains. Interestingly, the Pearson correlation of the allele frequencies of the observed variants of wild type and mutant with the public Pt1 was higher (0.78) then the correlation between wild type and mutant (0.61). Besides, the similarity of Pt1, wild type and mutant strain was evaluated by calculating Pearson Correlation and Euclidean distance scores using all the detected variants (Table 2). Both the measures highlighted that the WT strain is more similar to Pt1 than the mutant strain.