A large-scale genome and transcriptome sequencing analysis reveals the mutation landscapes induced by high-activity adenine base editors in plants

A few homozygous SNVs and indels were detected in all sequenced plants (Additional file 2: Fig. S5a, b). We counted the number of plants with the same mutation sites and found that the homozygous mutations tended to be present in more than one plant, while the heterozygous mutations tended to be present in a single plant (Additional file 2: Fig. S5c, d). These homozygous mutations could be the remaining background mutations or mutations induced by tissue culture, Agrobacterium infection, or ABEs. The induced mutations in the two alleles are two independent events following binomial distribution, so the probability of the homozygous mutations is p², the probability of being wild type (WT) is (1-p)², and the probability of the heterozygous mutations is 2 * p * (1-p), assuming that the induced mutation ratio for each allele was p and the ratio of the WT allele was 1-p. A binomial test for all loci of homozygous SNVs or indels revealed that these loci did not follow a binomial distribution (Additional file 2: Fig. S5e and Additional file 1: Tables S4 and S5), indicating that these homozygous mutations remain background SNVs and indels. These data suggests that ABEs induce sgRNA-independent heterozygous DNA mutations.

We next examined whether Cas proteins or TadA variants play distinct roles in inducing off-target DNA mutations by comparing plants harboring rBE46b with those harboring rBE49b as well as rBE50 to rBE53 to characterize TadA8e and TadA9, and compared plants with rBE46b to rBE50 and rBE49b to rBE53 to characterize the role of SpCas9n and SpCas9n-NG in off-target effects. Although there was no significant difference between TadA8e and TadA9 when the total number of SNVs was considered, plants harboring TadA9 had a higher number and a higher percentage of A>G SNVs (Fig. 2c), indicating that TadA9-based ABEs lead to a higher number of A>G SNVs. Plants harboring SpCas9n-NG had a higher number of SNVs as well as a higher number of A>G SNVs, but not a higher percentage of A>G SNVs (Fig. 2d), indicating that SpCas9n-NG-based ABEs lead to a higher number of SNVs.

We classified all SNVs into six types and calculated the percentage of each type of SNV versus the total number of SNVs. We observed a higher percentage of C>A/G>T SNVs in plants harboring TadA8e (Additional file 2: Fig. S7). We further mapped all SNVs and A>G SNVs to different genic and intergenic regions and calculated the ratio of SNVs in given regions versus in the whole genome. As a result, the number of A>G SNVs and the total number of SNVs were higher at all genic and intergenic regions in plants for all four types of ABEs, while A>G SNVs were enriched in genic regions and depleted in intergenic regions (Fig. 2e and Additional file 2: Fig. S8). In addition, we mapped total SNVs as well as A>G SNVs to the 12 rice chromosomes and established that they were distributed throughout the rice genome (Additional file 2: Fig. S9).

We next examined the integrity of T-DNA regions containing both a complete left border (LB) and right border (RB) and identified four plants with a partial T-DNA insertion characterized by the missing TadA8e, TadA9, or SpCas9n-NG fragment (Additional file 2: Fig. S11a). However, desired on-target mutations were detected in three out of four plants (Additional file 2: Fig. S3), suggesting that sgRNA-dependent on-target A>G editing could occur before T-DNA integration into the rice genome. We further checked the off-target SNVs between plants with or without complete T-DNA insertion and found that plants with a complete T-DNA insertion had a higher number of total SNVs, a higher number of A>G SNVs, and a higher percentage of A>G SNVs when compared to those with partial T-DNA insertion (Fig. 3c).

It was known that T-DNAs can be integrated in rice genome in more than one copy [28], so we divided the plants into two groups based on whether one copy or multiple copies of T-DNAs were integrated. We examined the number of total SNVs, the number of A>G SNVs, and the percentage of A>G SNVs in plants with rBE46b, rBE49b, rBE50 and rBE53 separately and did not observe a consistent influence of the copy number of T-DNA insertion (Additional file 2: Fig. S12).

We performed similar studies on DNA off-target SNVs but did not observe general patterns of flanking A>G editing. However, we did identify 25 loci with more than one A>G SNV from 12 plants (Additional file 1: Table S10); some loci contained 5–10 A>G SNVs, and others contained 2–3 A>G SNVs (Fig. 5d and Additional file 2: Fig. S18). Overall, 45% of these SNVs were located in the genic region, which is higher than the 30% observed for all A>G SNVs in the genic region, consistent with the tendency of off-target A>G SNVs to occur in the genic region (Fig. 5e). We classified these 12 plants into group 1, and the remaining 36 plants carrying ABEs into group 2. The number of SNVs and A>G SNVs and the percentage of A>G SNVs were significantly higher for plants in group 1 compared to plants in group 2 (Fig. 5f).

In this study, five rice (Oryza sativa) genomic loci (OsACC, OsGS1, OsMPK13, OsGSK3, and OsGSK4) and four rice genomic loci (OsACC, OsGS1, OsMPK13, and OsTms9) were targeted by rBE46b and rBE49b, respectively. Three genes (OsSERK2, OsDEP2, and OsGSK4) were targeted by both rBE50 and rBE53. Plant IDs and their corresponding information are described in Additional file 1: Table S1. The rBE46b, rBE49b, rBE50, and rBE53 expression plasmids were constructed as previously reported [16]. The empty entry vector without any spacer was cloned into pUbi:rBE46b, pUbi:rBE49b, pUbi:rBE50, and pUbi:rBE53 using Gateway technology to yield ABEs without sgRNAs (Additional file 1: Table S1).

The genome editing constructs were individually introduced into the Agrobacterium tumefaciens strain EHA105 via the freeze-thaw transformation method, and 2-week-old calli derived from immature seeds of the Geng rice variety Kitaake were infected by each Agrobacterium strain. After 4 weeks of culture on MSD medium supplemented with 50 mg/L hygromycin (Roche, Germany), the resistant callus lines were transferred onto RM plates to generate transgenic rice seedlings. All information on target gene mutations of each seedling examined in this study is given in Additional file 1: Table S1.

To eliminate background mutations, 10 individual Kitaake plants grown from seeds were used directly. Seedlings were regenerated from rice calli without Agrobacterium infection (namely C1) and regenerated from calli co-cultured with the empty EH105 strain (namely C2). Also, seedlings were regenerated from calli infected with EH105 strains harboring SpCas9 only (namely Cas). All rice materials were grown in the greenhouse under a 16-h-light/8-h-dark photoperiod, 28/25 °C temperature cycle, and 75% humidity.

Genomic DNA of 4-week-old rice plants was extracted using the CTAB method (Li et al., 2016). Approximately 200 mg of fresh rice leaves was collected in a 2-ml centrifuge tube containing disposable metal balls. After being quickly frozen in liquid nitrogen, samples were ground to a fine powder using a tissue grinding apparatus (Jingxin, China). Following chloroform extraction, isopropanol precipitation, and 70% EtOH washing, genomic DNAs were eluted with 50 μL of double-distilled water supplemented with 1 μL of 10 U/μL RNase I (Thermo Fisher Scientific, USA) and stored at − 80 °C for later experiments.

RNA was extracted with TRIzol reagent (Takara, Japan) according to the manufacturer's instructions. Briefly, 100 mg of fresh rice leaves was sampled, quickly frozen in liquid nitrogen, and ground to a powder with a tissue grinding apparatus. Then, 1 ml of TRIzol reagent was added to the sample followed by chloroform and isopropanol treatment. Finally, RNA pellets were dissolved in 50 μL of RNase-free water (0.1% DEPC-treated) and stored at − 80 °C for later experiments.

The on-target genomic regions were amplified using Phanta Max Super-Fidelity DNA Polymerase (Vazyme, China) and locus-specific primers (Additional file: Table S1, Table S11, Table S12) with genomic DNAs and cDNAs used as the template. PCR amplicons were subjected to Sanger sequencing, and Bioedit software was used for sequence data analysis. 1

RNA-free genomic DNAs (0.2 μg) from each sample were used to construct the DNA libraries using a NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's instructions. DNA libraries were sequenced on the Illumina platform in the 150-nt paired-end mode with an average coverage depth of 40× (Additional file: Table S2). 1

The clean reads were mapped to the Kitaake genome V3 from Phytozome (https://data.jgi.doe.gov/refine-download/phytozome↗) via BWA [41] and sorted using samtools (v1.9) [42]. The Genome Analysis Toolkit (GATK v4.2) was used to mark duplicated reads and recalibrate base qualities [25]. To identify high-quality genetic changes at the genomic scale, we applied three independent germline variant-calling methods: GATK, LoFreq [23], and Strelka2 [23]. We documented SNVs identified by all three methods and indels identified by GATK and Strelka. All genetic changes identified by the three methods in the 10 Kitaake plants were combined and used as background mutations. Sanger sequencing has been performed to validate the overlapping set of SNVs called by the three methods (Additional file 2: Fig. S19). The genetic mutation ratios were calculated using an in-house R program and 'AC' value from GATK's results. Both background mutations and homozygous mutations were removed from the SNVs as well as indels. The IGV browser was used to demonstrate sgRNA-directed on-target mutations [43]. Then, the on-target mutations were removed for off-target analysis. sgRNA-dependent off-target mutations were discovered using Crisflash [26], and the genetic on-target mutations were assessed using the IGV browser. A gene annotation file (OsativaKitaake_499_v3.1.gene_exons.gtf) from the Phytozome website was used to define different genomic regions, such as gene regions, exon regions, and intergenic regions. The ggpubr, ggbio, and VennDiagram R libraries were used to draw the graphs.

The clean reads were mapped to T-DNA sequences using BWA and sorted using samtools. The T-DNA insertion sites were located through T-LOC (Li et al. in preparation). The coverage of T-DNAs between the left border (LB) and right border (RB) was assessed using the R library ShortRead. The expression of ABEs was quantified as the average raw read number of Cas proteins and TadA variants normalized by the total read number in millions. Since we used T0 plants, the copy number of T-DNA integration was calculated as the relative T-DNA coverage versus half coverage of the rice genome.

DNA-free RNAs (0.2 μg) were used to construct the RNA-seq libraries using a NEB Next Ultra RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's instructions. RNA-seq libraries were sequenced on the Illumina platform in the 150-nt paired-end mode (Additional file: Table S8). 1

The clean reads were mapped to the Kitaake V3 genome and annotation from Phytozome via STAR aligner with a maximum of eight mismatches per paired-end read [44]. GATK was used to mark duplicate reads and split reads that contained Ns in their cigar string and to recalibrate base qualities. SNVs were called by GATK, LoFreq, and Strelka2 for each transcriptome dataset and corresponding genome dataset. The SNVs identified by three methods in the transcriptome data but not in the genome data were kept for later analysis. Sanger sequencing has been performed to validate the overlapping set of SNVs called by the three methods (Additional file 2: Fig. S20). All the genetic changes identified by the three methods in three Agrobacterium-infected plants were combined and used as background mutations and were removed from the SNVs identified in plants transformed with SpCas9, rBE46b, and rBE49b. The A>G mutation ratios of off-target RNA loci were calculated through in-house Python programs. The 30- and 3-bp flanking sequences of the off-target RNA SNVs were extracted from the Kitaake reference genome and subjected to motif prediction using WebLogo3 (http://weblogo. threeplusone.com↗ /) [45].

We combined all A>G off-target SNVs obtained from plants with RNA off-target activities. For each A>G SNV, we calculated the number of reads with nucleotide A, T, G, and C separately in the 5′ and 3′ 30-bp region with a read coverage larger than 10. The genetic change ratio was calculated as the number of Gs divided by the total number of As and Gs if the reference is A. The genetic change ratio was calculated as the number of Cs divided by the total number of Cs and Ts if the reference is T. Positions with an A>G mutation ratio of higher than 0.05 were used as the numerator, while positions of A/T with a read coverage larger than 10 were used as the denominator. Similarly, we combined all A>G off-target SNVs obtained from plants through WGS and calculated the percentage of A>G mutations at the 5′ and 3′ 30-bp flanking regions.

The horizontal line in the box represents the median value, and the bottom and top of the box are the lower (Q1) and upper quartiles↗ (Q3), respectively. The upper whisker is min(max(x), Q3 + 1.5 × IQR), and the lower whisker is max(min(x), Q1 − 1.5 × IQR). IQR (interquartile range) = Q3 − Q1. Black dots located outsides the whiskers are outliers.

Additional file 1: Supplementary tables. Table S1, Summary of plants with ABEs. Table S2, Mapping statistics of whole-genome sequencing. Table S3, Summary of sgRNA-dependent on-target and off-target loci. Table S4, Summary of all the homozygous SNVs. Table S5, Summary of all the homozygous indels. Table S6, Summary of genomic SNVs detected through WGS. Table S7, Summary of the overlapping SNVs between each of the plants with whole-genome sequencing. Table S8, Mapping statistics of whole-transcriptome sequencing. Table S9, Summary of all the transcriptomic SNVs. Table S10, Summary of clustered A>G DNA SNVs. Table S11, Primers used to verify DNA SNVs by Sanger sequencing. Table S12, Primers used to verify RNA SNVs by Sanger sequencing.Additional file 2: Supplementary figures. Fig. S1. Sanger sequencing chromatograms of on-target mutations in plants harboring rBE46b and rBE49b. Fig. S2. Sanger sequencing chromatograms of on-target mutations in plants harboring rBE50 and rBE53. Fig. S3. IGV browser views showing the on-target mutations for 36 plants harboring ABEs. Fig. S4. Analysis of SNVs and indels identified by whole-genome sequencing. Fig. S5. Analysis of the remaining background homozygous DNA mutations. Fig. S6. Characterization of ABE-induced genomic mutations. Fig. S7. Distribution of six types of SNVs. Fig. S8. Distribution of SNVs at given regions of the genome. Fig. S9. Chromosomal distribution of SNVs. Fig. S10. On-target and off-target mutations in plants from the same calli. Fig. S11. Off-target SNVs in plants with incomplete T-DNA insertions. Fig. S12. Distribution of SNVs with different copy numbers of T-DNA insertions. Fig. S13. Transcriptome-wide distribution of ABE-induced off-target mutations. Fig. S14. Heatmap demonstrating A>G mutations in transcriptomes with more than 5 A>G SNVs detected. Fig. S15. The 5′ and 3′ flanking A>G mutations in transcriptomes with ABEs containing A>G RNA SNVs and in transcriptomes with SpCas9 only lacking A>G RNA SNVs. Fig. S16. The 5′ and 3′ flanking A>G mutations in transcriptomes with ABEs but without A>G RNA SNVs. Fig. S17. IGV genome browser views showing the off-target RNA mutations. Fig. S18. IGV genome browser views showing A>G mutations with flanking A>G SNVs in genome sequencing data. Fig. S19. Sanger sequencing chromatograms of off-target A>G DNA mutations. Fig. S20. Sanger sequencing chromatograms of off-target A>G RNA mutations.Additional file 3: Review history