Fiol et al. Plant Methods (2022) 18:105 https://doi.org/10.1186/s13007-022-00937-4 Plant Methods METHODOLOGY Open Access An efficient CRISPR‑Cas9 enrichment sequencing strategy for characterizing complex and highly duplicated genomic regions. A case study in the Prunus salicina LG3‑MYB10 genes cluster Arnau Fiol1, Federico Jurado‑Ruiz1, Elena López‑Girona2 and Maria José Aranzana1,3*    Abstract Background: Genome complexity is largely linked to diversification and crop innovation. Examples of regions with duplicated genes with relevant roles in agricultural traits are found in many crops. In both duplicated and non-dupli‑ cated genes, much of the variability in agronomic traits is caused by large as well as small and middle scale structural variants (SVs), which highlights the relevance of the identification and characterization of complex variability between genomes for plant breeding. Results: Here we improve and demonstrate the use of CRISPR-Cas9 enrichment combined with long-read sequenc‑ ing technology to resolve the MYB10 region in the linkage group 3 (LG3) of Japanese plum (Prunus salicina). This region, which has a length from 90 to 271 kb according to the P. salicina genomes available, is associated with fruit color variability in Prunus species. We demonstrate the high complexity of this region, with homology levels between Japanese plum varieties comparable to those between Prunus species. We cleaved MYB10 genes in five plum varie‑ ties using the Cas9 enzyme guided by a pool of crRNAs. The barcoded fragments were then pooled and sequenced in a single MinION Oxford Nanopore Technologies (ONT) run, yielding 194 Mb of sequence. The enrichment was confirmed by aligning the long reads to the plum reference genomes, with a mean read on-target value of 4.5% and a depth per sample of 11.9x. From the alignment, 3261 SNPs and 287 SVs were called and phased. A de novo assembly was constructed for each variety, which also allowed detection, at the haplotype level, of the variability in this region. Conclusions: CRISPR-Cas9 enrichment is a versatile and powerful tool for long-read targeted sequencing even on highly duplicated and/or polymorphic genomic regions, being especially useful when a reference genome is not available. Potential uses of this methodology as well as its limitations are further discussed. Keywords: Cas9 enrichment, Complex regions, Long-read, Targeted sequencing, Variability, Gene duplication, MYB10 Background The origins of agriculture and the development of new crops are tightly linked to whole genome duplications *Correspondence: mariajose.aranzana@irta.cat (WGD), and/or to chromosome rearrangements, which 1 Centre for Research in Agricultural Genomics, CSIC-IRTA-UAB-UB, Campus are among the major determinants of useful trait diver- UAB, Barcelona, Spain sity [1–3] as well as of genome complexity [4]. In addition Full list of author information is available at the end of the article © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​ mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Fiol et al. Plant Methods (2022) 18:105 Page 2 of 16 to WGD, small-scale duplications (SSD), usually pro- 22]; however, although the cost of long-range sequencing duced by segmental duplications and gene expansion, methods is decreasing, the availability of large amounts also play an important role in crop innovation [5, 6]. of high quality DNA, whole genome data storage, align- Reduced selective pressure acting on duplicated genes ment and computation costs may be drawbacks when may, in some cases, result in them acquiring novel func- the objective is to resolve and explore the variability in tionalities that contribute to diversification and lead to a large panel of varieties of a unique complex loci for an important traits for agriculture [7–10]. interesting agricultural trait. Gene diversification results in high levels of diversity Over the last few years, several skimming and targeted- within gene families. For example, a substantial number enrichment sequencing methods have been released and of copies of celiac disease related α-gliadin genes are esti- successfully used to pinpoint small variants in single mated in wheat (ranging from 25–35 to 100–150 cop- regions, especially in phylogenomics. However, a cost- ies per haploid genome), and are probably generated by effective method to scan the variability (including SV) duplication, deletion events and retrotransposon inser- in a complex or highly diverse region in a panel of geno- tion [11, 12]. These genes show high variability between types (for example varieties or seedlings) using LR-HTS genomes [13], with at least half of the copies being inac- strategies is still lacking. tive or pseudogenes [14]. Recently, a methodology has been reported for long- In both duplicated and non-duplicated genes, a large read targeted sequencing using CRISPR-Cas9 technology fraction of the variability is caused by large as well as to direct the sequencing adapters to the regions of inter- small and middle scale structural variants (SVs). For est [23, 24]. In this procedure, selected regions of high example, Chawla et al. [15] found that up to 10% of all molecular weight dephosphorylated DNA are cut by the genes in Brassica napus were affected by small- to mid- Cas9 enzyme directed by specific guide RNAs (gRNAs). scale SV events. In tomato, Alonge et al. [16] found The digested fragments are enriched by long-read multiple SVs, in many cases mediated by transposable sequencing thanks to the preferential ligation of sequenc- elements (TEs); many of the SVs responsible for changes ing adapters towards the 5’ phosphate groups generated in gene dosage and expression levels modified agro- at the cleaved ends. Gilpatrick et al. [23] used a pool of nomic traits such as fruit flavor, size, and production. gRNAs to cut and sequence several loci associated with Recent studies of SVs between and within plant species human breast cancer on the same run, resolving SNPs have revealed extensive genome content variation [17], and SVs at the haplotype level. This methodology is sim- leading to the pan-genome concept, which is centered ple to perform, as it does not require physical separa- on the study of the entire gene repertoire of a species by tion of the DNA for sequencing, cloning or amplification sequencing multiple individuals [18]. steps, thus allowing native DNA strands to be read, and The major effect of duplicated genes and SVs, especially it is not affected by amplification bias and allows visu- in crop agronomic traits, highlights the importance of alization of allele-specific methylation patterns [23, 25]. studying the variability between genomes in these genes CRISPR-Cas9 enrichment has recently been used on a for plant breeding. plant genome (as a proof of concept) showing that, at the Whole-genome sequencing (WGS) methods have con- haplotype level, it can resolve the SV that causes red flesh tributed enormously to broadening our understanding of color in apples [24]. Although this methodology was used the genetic basis of useful traits. However, genome com- to sequence a short (7.8 kb) and low complex region in plexity still represents a challenge for genome assembly one apple variety, the successful results augur well for its and annotation. In general, large divergent duplications use for fine-mapping in other situations. with less than 97% homology can be resolved by WGS The objective of this study was to evaluate the potential assembly, whereas reads of duplications with higher of CRISPR-Cas9 enrichment to identify the variability homology are frequently collapsed, producing assembly (in particular SV) in highly complex regions using a pool errors [19]. of DNA of highly diverse varieties. The region selected A large number of genomes have been obtained, and was the Japanese plum (Prunus salicina) MYB10 region SNPs and small InDels discovered using short-read high located on linkage group 3 (LG3), found recently to con- throughput sequencing (SR-HTS), but long-read high tain at least three MYB10.1 (MYB10.1a, MYB10.1b and throughput sequencing (LR-HTS) is required to resolve MYB10.c) genes, one MYB10.2 gene and one MYB10.3 highly homologous duplicated regions and SVs, includ- gene, while additional copies could not be discarded ing repetitive regions and different forms of copy num- [26]. A marker system designed in this region was able ber variations, such as presence-absence variants [20]. to characterize the LG3-PsMYB10 allelic organization Multiple reference genomes for crop species are cur- in Japanese plum varieties and progenies into six haplo- rently being generated using LR-HTS methods [16, 21, types (H1 to H6), each with a different combination of Fiol et al. Plant Methods (2022) 18:105 Page 3 of 16 the MYB10.1–3 alleles [26]. For example, H1 and H3 hap- of five extra exons along the next 6 kb of sequence. lotypes contained three PsMYB10.1 alleles, one of them By BLAST we identified three additional genes: one (MYB10.1a-a356) associated with the anthocyanin skin was homologous to MYB10.1 and two homologous to color; polymorphisms in the promoter and intron of this MYB10.2. In the Zhongli-2 region, four MYB10 genes allele in H1 and H3 were observed, confirming the high were annotated: two MYB10.1 copies, one MYB10.2 and variability of this region. one MYB10.3. The second exon was also mis-annotated We designed seven gRNAs in conserved sites of the in most of the gene copies of this region. BLAST analysis PsMYB10 genes and pooled them, together with two did not identify additional MYB10 genes. guides flanking the region, to allow for multiple Cas9- To identify and visualize the differences in size and mediated cuts in the DNA of five Japanese plum com- gene copy number and organization, we produced dot mercial varieties, all with a distinct and heterozygous plots comparing the three LG3-MYB10 region assem- genotype for the MYB10 region, while sharing one haplo- blies (Fig. 1). These reflected the high variability and poor type by pairs. The fragments were labeled with a barcode homology along the region, with inverted assemblies, system and pooled and sequenced in a single MinION mis-assemblies and missing fragments. The Zhongli-1 Oxford Nanopore Technology (ONT) run. Reads align- MYB10 region was assembled in an inverted orientation, ment with reference genomes or de novo alignments suc- when compared with ‘Sanyueli’ and Zhongli-2. Despite cessfully identified SNPs and SVs, most of them newly this, the Zhongli-1 (of 271 kb) and the Zhongli-2 (90 kb) discovered and others observed in our previous work, so alignments were collinear and with high homology (69% validating the methodology for sequencing and variant of sequence hits), with an evident gap of ~ 150 kb due to detection of regions containing complex duplications. a missing fragment in Zhongli-2 explaining most of their difference in size (271 kb vs 90 kb). The percentage of hits Results between ‘Sanyueli’ and Zhongli-1 was lower (51%), with The Japanese plum LG3‑MYB10 region is highly complex multiple misalignments along the region that explained We used the Japanese plum LG3-MYB10 region as a the smaller size in ‘Sanyueli’. Similarly, a discontinuous model of high complexity. Japanese plum is a self-incom- alignment and an inversion was also observed between patible fruit tree with a complex history of interspecific Zhongli-2 and ‘Sanyueli’ (51% of sequence hits). crosses between diploid plums: the LG3-MYB10 region Similar comparisons were constructed for other Pru- in Prunus genomes, particularly in Japanese plum, is nus species with more than one available reference highly diverse and contains a cluster of at least three par- genome (Additional Files 1, 2). The MYB10 regions of the alogous MYB10 gene copies. two sweet cherry (P. avium) genomes were highly col- To test the complexity of this region and validate its linear with 75% of sequence hits. In peach (P. persica), suitability for analysis, we identified and compared the the region was collinear and highly similar in the ‘Lovell’ assembly of the LG3-MYB10 region in the two Japanese and ‘Chinese Cling’ genomes (85% of sequence hits), and plum genomes currently available, one for ‘Sanyueli’[27] three times larger than in the wild relative P. mira, with a (https://​www.​rosac​eae.​org/​Analy​sis/​94507​78) and one mean of 55.5% of sequence hits. In apricot, we compared for ‘Zhongli No. 6’ [28] (https://​www.​rosac​eae.​org/​Analy​ five P. armeniaca accessions and three wild apricot spe- sis/​90196​55). In the ‘Sanyueli’ genome (v1.0) the region cies (two P. sibirica and one P. mandshurica). Homology spanned 135 kb in chromosome 4 (Chr4:12,192,580– between each P. armeniaca pairwise comparison ranged 12,327,184). Although the genome annotation contained from 44 to 91% (mean of 67.7% of hits). When these only one gene, homologous to PsMYB10.2, BLAST iden- were compared to the wild apricot genomes, the homol- tified an additional PsMYB10.2, two PsMYB10.1 and one ogy ranged from 58 to 93% (mean of 74.8% of hits) with PsMYB10.3 gene copies in the region. P. sibirica, and from 47 to 59% (mean of 53.2% of hits) In the ‘Zhongli No. 6’ genome we identified two MYB10 with the P. mandshurica genome. The mean percentage regions 2 Mb apart in chromosome 3: one 271 kb long of sequence hits between the Prunus sections was 19.4% (from now on, Zhongli-1; LG03:30,838,572–31,109,359) with values ranging from 8.6% to 37%. and the other 90 kb long (Zhongli-2; LG03:28,592,935– 28,682,529). The Zhongli-1 region had five MYB10 genes The designed gRNAs successfully cut LG3‑MYB10 genes annotated, corresponding to two copies of MYB10.1, For selective high-throughput sequencing of the desired two copies of MYB10.2 and one copy of a MYB10.3 gene. LG3-MYB10 region in Japanese plum varieties, we These genes were mis-annotated: the second exon was designed seven CRISPR RNAs (crRNAs) in genic regions omitted in all but one of the PsMYB10.1 copies, while (exons and introns) of the PsMYB10.1, PsMYB10.2 and in this copy the STOP codon in the last exon was not PsMYB10.3 genes, including regions conserved between detected, causing the spurious automatic identification the genes to minimize the number of required crRNAs. Fiol et al. Plant Methods (2022) 18:105 Page 4 of 16 Fig. 1 a Dot plot comparison of the three LG3-MYB10 assemblies in the two Japanese plum genomes published: ‘Zhongli No. 6’ (with two assemblies: Zhongli-1 and Zhongli-2) and ‘Sanyueli’. Top right of the diagonal: homologous hits; inversions colored in red. Bottom left of the diagonal: homologous blocks colored by their score, from green to colder colors while the homology value decreases. Coordinates of each region are written in base pairs for both axes. b Distribution of the homologous sequences in the three compared genome regions. Grey areas correspond to sequences not found in the Zhongli-1 MYB10 region. Homologous regions are represented in the same yellow to purple colors The crRNAs were designed at the two DNA strands to possible SNP in the second exon of the MYB10.2 cleav- ensure forward and reverse sequences and to direct the age site of crRNA-s1 (although only identified in P. ceras- sequencing from the inside to the outside of the gene, ifera, P. domestica and P. avium) we designed crRNA-s2 covering the whole region with concatenated sequences for the alternative nucleotide. As reported by Fiol et al., (Fig. 2). The PsMYB10.1, PsMYB10.2 and PsMYB10.3 variability in the introns was higher than in exons, so sequences obtained in Fiol et al. (2021), the peach ‘Lovell’ the crRNAs were designed taking into consideration the reference genome, and other Prunus genomes were used variants known. Introns of PsMYB10.1, PsMYB10.2 and to design the crRNAs with higher on-target capacity and PsMYB10.3 were targeted with crRNA-a1, crRNA-a3 and lower off-target risk (see Methods). Therefore, crRNA- crRNA-a4, respectively, while crRNA-a2 was designed to s1 was designed to cut within the exon 2 of MYB10.1 include an SNP variant in one of the PsMYB10.1 alleles and MYB10.2 genes in the plus strand, while crRNA-s3 present in the sample (PsMYB10.1-H1 in [26]). In addi- was directed to the same exon and strand in MYB10.3 tion, to cut the MYB10 genomic region at each flank, (Additional File 3). To avoid a putative mis-cut due to a two crRNAs (crRNA-f1 and crRNA-f2) were designed in Fiol et al. Plant Methods (2022) 18:105 Page 5 of 16 Fig. 2 Schematic representation of the CRISPR RNA (crRNA) design to specifically cleave and sequence a hypothetic LG3-MYB10 genomic region in Japanese plum. The guide RNAs (gRNAs) are formed by the interaction of trans-activating RNAs with crRNAs, then with the Cas9 enzyme are assembled into the Cas9- Ribonucleoprotein (Cas9-RNP) complex. The Cas9-RNPs allows the specific cleavage of the LG3-MYB10 region, generating cuts at the two DNA strands for further sequencing in both directions from the inside to the outside of the genes. Gene number, order and size of the region, as well as the expected number of fragments were unknown at the time of the crRNA design conserved exon sequences of the nearest flanking genes MYB10.2 and MYB10.3 complete genes in the plum vari- upstream (Prupe.3G162900 in peach, Pav_sc0000464.1_ ety ‘Angeleno’ and cleaved the products with a pool of g320.1.mk in sweet cherry genomes) and downstream the seven gRNAs targeting the MYB10 genes only. As (Prupe.3G163400 in peach, Pav_sc0000464.1_g090.1.mk intron 2 in one of the PsMYB10.1 alleles in ‘Angeleno’ is in sweet cherry genomes). (Fig. 2, Additional File 3). about 1.5 kb larger [26], the MYB10.1 primers yielded We first validated the cleavage ability of the crRNAs in fragments visualized in two bands in Fig. 3a (the one PCR amplified DNA. For that we amplified the MYB10.1, containing the large intron is very faint in the gel, due Fig. 3 a MYB10 amplicon bands before and after cleavage with a pool of the guide RNAs (gRNAs) targeting the genes, with their size annotated in base pairs (bp). Bands (a, b) correspond to the faint bands from the allele with a large intron 1, (c) to the remaining undigested product. Wells M1: Lambda DNA HindIII-EcoRI digested molecular weight marker; M2: DNA Ladder 50 bp ready-to-use (GeneON). b Expected fragment sizes after MYB10.1 and MYB10.2 digestion with the Cas9-Ribonucleoprotein (RNP) complex. Guide-RNAs assembled with the trans-activating CRISPR RNA (tracrRNA) and crRNA-s (s1, s2 or s3) cut the exon 2 of the gene, gRNAs with crRNA-a (a1, a2, a3 or a4) cut by intron 1. Fragment 6 was not obtained after cleavage, fragment 3 and 5 bands overlapped in the gel Fiol et al. Plant Methods (2022) 18:105 Page 6 of 16 to the preferential amplification of the shorter one). The the output reactions was sequenced with the ONT Min- fragments were cut and visualized in five bands. Three ION device, which produced 194 Mb of sequences in of them corresponded to the fragments 2 (526 bp), 3 33,794 reads of average mean length of 6,097 bp (Table 1). (1078 bp), 4 (613 bp) and 5 (991 bp) in the in-silico design After demultiplexing the reads into varieties through shown in Fig. 3b, with fragments 3 and 5 overlapping in their barcode, we obtained sequence yields ranging from the gel. The two faint additional bands corresponded also 29.32 Mb to 45.64 Mb in reads of mean length from 4,110 to fragments 3 and 5 (of about 2.5 kb) but of the allele to 7,778 bp; the average N50 was 14,510 bp. with the large intron 2, and to an undigested fraction of The reads of each of the five varieties were aligned the product of the preferentially amplified copy. Similar against the ‘Sanyueli’ and ‘Zhongli No. 6’ genomes. The results were observed for the MYB10.2 amplicon, which coverage and depth of the alignments within the MYB10 was cut into four fragments. The band corresponding to region varied for each sample and reference genome the small gene fragment (fragment 6 in Fig. 3b) spanning used (Table 2, Additional File 4). Higher coverage was the exon and intron crRNA sites was not produced in observed for the alignments against the shorter region either the MYB10.1 or the MYB10.2 cleavages. The pool of Zhongli-2 (90 kb long, 63.4% of coverage in average) was not able to cut the MYB10.3 amplicon. A posterior while the coverage of the alignments against Zhongli-1 sequencing of that amplicon revealed an unexpected 210- (54.2% of coverage, 271 kb long) was only slightly lower bp deletion affecting part of intron 1 and all of exon 2, than that for ‘Sanyueli’ (57% of coverage, 135 kb long). and so lacking the two scission target sites. The on-target alignment did not correlate with the num- ber of reads yielded, which was uneven for each variety. Sequencing, alignment to the reference genomes In particular, ‘Golden Japan’ was the variety with the and variant calling lowest number of on-target aligned reads (168, 67 and The genomic DNA of five plum varieties (Table 1) was 45 in ‘Sanyueli’, Zhongli-1 and Zhongli-2, respectively). cleaved with the pool of gRNAs and barcoded. A pool of The mean sequence depth was also unequal for each Table 1 Sequencing statistics for each plum variety after read demultiplexing Variety Barcode Total reads Yield (Mb) Mean length (bp) N50 Mean Q ‘Angeleno’ NB01 4995 38.85 7777.9 19,049 11.0 ‘Black Gold’ NB02 11,106 45.64 4109.7 10,288 10.9 ‘Fortune’ NB03 5417 29.32 5412.2 11,786 10.8 ‘Golden Japan’ NB04 5063 34.98 6908.0 16,740 11.0 ‘TC Sun’ NB05 7213 45.27 6275.8 14,686 10.9 Table 2 Alignment of each cultivar to the three Japanese plum MYB10 reference regions Sample MYB10 region ‘Sanyueli’ Zhongli-1 Zhongli-2 (135 kb) (271 kb) (90 kb) On-target Cov. bases Mean depth On-target Cov. bases Mean depth On-target Cov. bases Mean depth reads (%) reads (%) reads (%) (%) (%) (%) ‘Angeleno’ 1,036 (5.63%) 73.3 kb 12.2x 489 (4.17%) 96.7 kb 3.8x 486 (4.15%) 48.7 kb 9.1x (54.49%) (35.7%) (54.34%) ‘Black Gold’ 1,300 (5.37%) 76.8 kb 19.1x 572 (2.58%) 123.0 kb 5.5x 611 (2.75%) 77.7 kb 14.1x (57.07%) (45.42%) (86.31%) ‘Fortune’ 1,692 78.1 kb 23.2x 841 (8.60%) 165.1 kb 9.1x 616 (6.30%) 72.7 kb 18.5x (10.47%) (58.05%) (60.97%) (81.18%) ‘Golden 168 71.1 kb 2.0x 67 (0.63%) 133.8 kb 1.2x 45 (0.43%) 31.6 kb 1.0x Japan’ (0.99%) (52.85%) (49.42%) (35.31%) ‘TC Sun’ 1,824 (7.04%) 84.5 kb 29.1x 677 (6.23%) 214.8 kb 19.3x 280 (2.58%) 53.6 kb 11.1x (62.75%) (79.33%) (59.84%) On-target reads: reads that aligned to the region and the percentage from the total reads aligned to the genome; Cov. bases: the bases covered with at least one read and the percentage of the region size; Mean depth: the mean depth of coverage for the whole region Fiol et al. Plant Methods (2022) 18:105 Page 7 of 16 variety-reference genome combination. In general, for all ‘Angeleno’ and ‘Black Gold’; H3 by ‘Angeleno’ and ‘For- varieties, a higher depth and higher number of on-target tune’; and H4 by ‘TC Sun’ and ‘Golden Japan’). This was aligned reads was observed with ‘Sanyueli’ (Table 2). The used to assign the phase of some of the polymorphisms: variety with the lowest depth was ‘Golden Japan’ due to ‘Angeleno’ and ‘Black Gold’ had 259 private SNPs which the poor alignment, while ‘TC Sun’ performed better. The may correspond to the shared H1; ‘Angeleno’ and ‘For- off-target regions in ‘Golden Japan’ corresponded mostly tune’ had 62 that could correspond to the shared H3. to small regions (< 1 kb) scattered all over the genome, ‘Fortune’ and ‘Black Gold’ had 87 private SNPs which homologous to the chloroplast DNA sequence and align- should be in H2 and H6. Unfortunately, due to the poor ing reads clipped from both sides. The same off-target quality of the ‘Golden Japan’ sequencing the phase of reads were observed aligned to these regions in the other the SNPs in H4, H5 and H9 could not be inferred. sample sequences, indicating that it was not specific of Sniffles [29] and NanoVar [30] software were used ‘Golden Japan’ and the low-depth in its MYB10 region to call SVs and polymorphisms affecting more than 30 was caused by a poor Cas9 digestion. bps, identifying 287 variants: 72 in ‘Angeleno’, 63 in In general, the sequence depth was unevenly distrib- ‘Black Gold’, 68 in ‘Fortune’, 30 in ‘Golden Japan’ and 54 uted along the MYB10 region, with scattered regions in ‘TC Sun’. Most corresponded to sequence breakends of high depth (Additional File 4). This was caused by (i.e., specific positions where aligned reads were clipped sequences that aligned partially and were clipped from because they lost the homology with the reference one or both sides, indicating low homology with the sequence) (47.39%), followed by deletions (27.87%), genomes used as reference. Sequence clipping was not insertions (12.89%), tandem duplications (8.36%) and observed in ‘TC Sun’ aligned to Zhongli-1. inversions (3.48%) (Fig. 4b). Variants were called for the alignment with ‘Sanyueli’, The accuracy of the pipeline was tested by search- the one with the highest depth. A total of 3,261 SNPs in ing a previously Sanger-validated insertion of 44 bp in the MYB10 region were found, 62.34% of them (2,033) the PsMYB10.1a promoter in the haplotypes associ- being present in at least two varieties (Fig. 4a), 98.5% ated with red skin color (H1, H3) and in H9. Sniffles with the same alternative allele (Additional File 5). All detected the polymorphism in homozygosity on ‘Ange- six observed MYB10 haplotypes (H1–H6) defined by leno’ (H1/H3) and in heterozygosity on ‘Black Gold’ Fiol et al. (2021) plus one inferred (H9) were present in (H1/H2); NanoVar could detect it even in low depth on the five varieties, three of them shared by pairs (H1 by ‘Golden Japan’ (H4/H9). The insertion was not called in Fig. 4 a Venn diagram showing the number of SNPs identified in the MYB10 region from each of the five Japanese plum samples sequenced and the overlap between them. The number next to the cultivar stands for their total SNP count. b Classification and percentage of all the structural variants identified: breakends (BND), deletions (DEL), duplications (DUP) and inversions (INV) Fiol et al. Plant Methods (2022) 18:105 Page 8 of 16 ‘Fortune’ (H3/H6), although some reads were found in only. We used as an example a 250 bp contig upstream low frequency. PsMYB10.1 (previously sequenced with Sanger by Fiol et al. [26]), which included the 44 bp sequence previ- de novo assembly and homology visualization ously used to validate the SV call, plus a 41 bp insertion Despite being able to identify new polymorphisms with and several SNPs differing between haplotypes. We suc- this method, we carried out a de novo assembly to over- cessfully identified the SNPs and InDels in heterozygo- come the low coverage caused by sequence clipping in sis, among them the 41 bp deletion in H3 and the 44 bp the regions with low homology. First we constructed an insertion in the H1, H3 and H9 haplotypes. The former assembly for each variety, then we used the ‘Zhongli No. (position 66 to 76 in the contig, Fig. 5) was shared by 6’ and ‘Sanyueli’ genomes to scaffold the contigs along ‘Angeleno’ (H1/H3) and ‘Fortune’ (H3/H6). The latter the region. As expected, the degree of homology and col- (position 200 to 243 in the alignment, Fig. 5) was pre- linearity varied for each assembly (Additional File 6). In sent in the varieties with H1, H3 or H9, in agreement general, higher homology was observed in comparisons with the Sanger sequencing results of Fiol et al. [26]. of the assemblies of the five varieties sequenced with the Most of the SNPs detected were consistent with those CRISPR-Cas9 selective cleavage strategy, particularly in previously found, however we identified a few sequenc- those sharing a haplotype. In particular, the assemblies of ing errors. This was the case of the SNPs in positions ‘Angeleno’ and ‘Black Gold’ were always the most similar, 52, 137, 175, and 244–245 in the alignment (Fig. 5). with hits ranging from 41 to 53%, depending on the refer- Despite these few sequencing errors, intrinsic to the ence region used for assembly. ONT sequencing technology, the results validate this methodology for the resolution and identification of Identification of polymorphisms in the PsMYB10.1a variability in highly complex regions. promoter Given the low homology between the samples and the ‘Zhongli No. 6’ and ‘Sanyueli’ genomes, we searched for polymorphisms among the de novo assemblies Fig. 5 Visualization of the de novo sequences in the PsMYB10.1a promoter region containing a G-box motif (underlined) associated with the red skin color, SNPs and two large InDels. The de novo sequences had resolution to the SNP level when compared between them and to the Sanger sequences provided in Fiol et al. Heterozygous positions (shadowed bases) on the de novo sequences matched the nucleotide on the same haplotype in a different variety (red rectangles). Sequencing mistakes are marked with asterisks Fiol et al. Plant Methods (2022) 18:105 Page 9 of 16 Discussion genome level within the Prunus genus [42] and confirms In this work we confirm sequence enrichment CRISPR- the high complexity of this region in Prunus in general Cas9 selective cleavage as a cheap and efficient high and in Japanese plums in particular. throughput strategy to identify variants and validate its Previously in Fiol et al. [26], we observed six haplo- use in highly complex regions in heterozygous genomes. types from the segregation of PsMYB10 alleles in six F1 Complex regions frequently include small and large vari- Japanese plum progenies and advanced breeding lines, ants, transposable elements (TE) and clusters of tandem finding high allelic variability and the triplication of the or segmental duplicated multi-copy genes. Such genes MYB10.1 gene in some haplotypes. Using the limited are more likely to carry higher levels of variation [31] and homology between Prunus genomes and by designing be affected by presence/absence variation (PAV) than numerous PCR primer combinations, we were able to singleton genes [32], and often intervene in processes clone some of the alleles and phase the variants. This was related to the evolution of agronomic traits [31, 33–36]. a slow and costly strategy, and we were not able to obtain This makes characterizing such regions highly relevant the full sequence of the haplotypes. In contrast, here we for evolutionary studies as well as for plant breeding [37]. complemented the Cas9 enrichment strategy with a mul- It is difficult, or even impossible, to completely resolve tiplex approach to decrease the cost of sequencing six plant genome complexity by short read sequencing. of the observed haplotypes (named H1 to H6), which Long-read approaches have improved the contiguity of allowed to quickly obtain high-quality sequences along genome assemblies, while special effort is currently being the MYB10 region for each of them. made to resolve particularly complicated regions of rela- tively small sizes in selected genotypes. In humans, sev- Design of the crRNAs eral authors have used a strategy based on enrichment The MYB10 proteins belong to one of the largest tran- and sequencing using targeted cleavage of chromosomal scription factors (TF) groups, the R2R3-MYB class, one DNA with Cas9 combined with Nanopore sequencing of the four MYB gene classes in seed plants [43]. As to resolve genomic SVs as well as mobile elements [23, occurs with most multi-copy gene families, they share 38–40]. In plants, this strategy has been validated to fine conserved motifs within the gene sequence. Such regions map short regions [24] and to characterize the transpo- can be successfully used to direct crRNAs. We designed son insertion landscape [41]. Here we extend this strat- seven non-overlapping forward and reverse crRNAs able egy for the characterization of much longer regions in a to provide sequences of the genes and their upstream pool of genotypes in species without a reference genome and downstream regions. This work was started several or in regions with large variability in gene content. months before the two Japanese plum genomes were As the region of interest (ROI), we selected the Prunus released. Although their availability could have helped in chromosome 3 MYB10 region, which consists of a clus- designing the crRNAs, our results indicate that this strat- ter of MYB10 gene copies with a gene content variable egy can be successfully applied to those species lacking a among Prunus species. reference genome. To reflect the complexity within the ROI we ana- Due to the small length of the MYB10 sequences avail- lyzed the synteny within and between Prunus sections, able and the high variability in the region, the on-target including two recent P. salicina genomes (one for the activity of the crRNAs was far from optimal (51 to 68, variety ‘Sanyueli’ and one for ‘Zhongli No. 6’) that were on a scale of 1–100). Their cleavage activity was tested released when our study was already at a very advanced in vitro with MYB10 amplicons of the variety ‘Angeleno’. stage. Despite belonging to the same species, the level The results revealed that while cleavage would likely not of homology in this region of these two Japanese plum perform well for in vivo gene editing, they were suffi- varieties was low, comparable to that observed between cient for in vitro cleavage. In addition, in vitro digestion cultivated and wild members from the same section of revealed that Cas9 cannot simultaneously excise two tar- apricot and peach. By contrast, the homology within each geted regions when these are only a few base pairs apart of the Prunus sections studied (apricot, sweet cherry and (about 100 bp in our design). This may be due to the Cas9 peach) was higher or moderately higher. Such low level remaining bound in the cleaved DNA sequence on the of homology in the ROI within P. salicina could be attrib- 5’-side of the gRNA, allowing preferential ligation of the utable to the interspecific origin of the species together sequencing adapters to the available end [23, 44, 45]. with a putative mis-assembly of the region in the ‘Zhongli The pool of crRNAs did not cleave the MYB10.3 ampli- No. 6’ genome, in which two LG3-MYB10 regions were con tested because one of the gene copies missed both assembled 2 Mb apart instead of the single one expected cleavage target sites, revealing the existence of a new for Prunus. The low identity between the Prunus ROI allele, probably an unknown duplicated MYB10.3 copy. contrasts with the high synteny described at the whole When we later visualized the genome alignments, reads Fiol et al. Plant Methods (2022) 18:105 Page 10 of 16 were observed starting from the MYB10.3 annotated reads, we would not have obtained much higher coverage position in all the samples, indicating that the Cas9 could or depth. The complexity of the region, its size, and the indeed cleave this gene. number of cutting sites need to be considered in deciding the number of pooled varieties. Sequencing and alignment to reference genomes The crRNAs were pooled to cleave the DNA of each vari- Variant identification in the alignments to a reference ety in a single reaction. This strategy made it possible to genome enclose the ROI and cleave the MYB10 genes along the The strategy described here can be successfully used to whole region, regardless of their copy number and organ- rapidly extract useful information for marker-assisted ization, which was unknown. selection. The generation of long reads offers a higher Pooling the barcoded cleaved DNAs lowered the mapping certainty [46], which is extremely useful for sequencing costs, although it reduced the number of regions with segmental duplications to pinpoint sequence reads obtained per variety, which could have risked the variants with high confidence levels. depth of coverage. We obtained 38.8 Mb per variety, on We identified SNPs and SVs along the ROI involved average, which represents 288 times the number of bases in the regulation of fruit skin color. For SNP calling we in the ROI of ‘Sanyueli’, 433 in that of Zhongli-2 and 143 used Longshot software [47], which not only considered times in Zhongli-1. The depth and coverage obtained the error tendency of the ONT long-reads to avoid false after alignment were lower than expected from these positives but has also been reported to outperform the figures. The lowest depth values were for the alignments SNP calling method using only short-read sequences with the ‘Zhongli No. 6’ genome, which was reason- in regions with duplications. A total of 129 SNPs were able considering that sequences were aligned to the two shared by the three cultivars with red-to-purple skin and assembled MYB10 regions, therefore reducing the depth absent in the two with yellow hues; since their phases are at each of them. However, the depth in the alignments known, they could be selected, after validation, as alter- with ‘Sanyueli’ was only slightly better (from 12.2 × to natives to the haplotype (PCR) marker provided in Fiol 29.1x, if we dismiss the extremely low values of ‘Golden et al. [26]. Japan’). The depth is highly determined by the percent- For SV calling we used Sniffles software [29] comple- age of coverage and of the on-target reads. We did not mented with NanoVar [30] to overcome the low depth of observe a correlation between the number of reads and ‘Golden Japan’. The results indicate that the tools comple- the coverage nor with the on-target ratios. Indeed, the mented each other well in overcoming their limitations, sequencing statistics for ‘Golden Japan’ were similar and that SVs could still be called even in low depth. In to those obtained for the other varieties, while all the addition, we were able to phase the variants, finding the alignment values were extremely low, indicating a poor 8 bp insertion in exon 1 of the PsMYB10.1a gene in the enrichment probably caused by an inefficient Cas9 cleav- H9 of ‘Golden Japan’ in phase with the 44 bp insertion age reaction. associated with red skin. This 8 bp insertion produces The ratio of on-target alignments (mean of 7.13% in an early STOP codon, which explains the yellow skin of ‘Sanyueli’ after discarding ‘Golden Japan’ value) was this cultivar despite it carrying the 44 bp insertion (Addi- higher than in other studies (4.61% for multiple cuts in tional File 7). Gilpatrick et al. [23]; 2.08 and 3.04% in López-Girona et al. [24]. This could be attributable to the higher num- De novo assemblies ber of cleavage sites designed, generating more Cas9- We have shown that reads can be aligned to a reference digested DNA fragments with compatible ends for their genome or assembled de novo in the case of its absence posterior ligation with the sequencing adapters. This or of bad homology. Here, to recover the positions lost conclusion was also reached by McDonald et al. [40], who due to the lack of homology with the reference genomes, obtained a mean on-target value of 44% in the CRISPR we used the pipeline for read correction and subsequent Cas9-enrichment of TEs. de novo assembly provided in López-Girona et al. [24]. Since the complexity of the ROI is out of our con- One of the main limitations in our study was that reads trol, higher depths would be more easily obtained by produced from two very distant adjacent genes had lit- increasing the number of cleavage sites. As the crRNAs tle chance to overlap. This was the case in the Zhongli-1 are pooled in a single cleavage reaction, the cost of the region, where two adjacent genes 124 kb apart left 55 kb method is only increased by the cost of the crRNAs of the region uncovered. The other limitation for obtain- employed, which is a small proportion of the total cost. ing a whole-region assembly was that reads could not In addition, we consider that by reducing the number of overlap the cleavage positions because these were pro- pooled varieties, and therefore increasing the number of duced in a single reaction. Therefore, despite reducing Fiol et al. Plant Methods (2022) 18:105 Table 3 Strengths, weakness, improvements and new opportunities of the CRISPR-Cas9 sequencing enrichment strategy described in this study Strengths Weakness Improvements New opportunities The design of crRNAs does not necessarily The search for polymorphisms in the de novo Digesting the DNAs with sub-pools of crRNA Identification of methylation variants between require a reference genome sequences might require the manual isolation might improve the region assembly, at the the pooled samples in the target region [23, 25] Pooling the crRNAs and barcoded DNAs results and comparison of contigs expenses of simplicity and cost-effectiveness Enrichment of genes scattered in more than one in a simple and cost-effective method ONT technology is prone to sequence errors in The use of samples homozygous for the target region of interest, such as the enrichment of Polymorphisms can be extracted and phased homopolymer regions [48] region might increase the efficiency of polymor‑ specific gene or transposon families[41] with high efficiency from genome alignments The Cas9 digestion with a single crRNA pool phism discovery and sequence scaffolding and/or de novo sequences cannot produce overlapping fragments The method is computationally inexpensive Page 11 of 16 Fiol et al. Plant Methods (2022) 18:105 Page 12 of 16 the cost, cleaving in a single reaction reduced the accu- The quality of the extracted DNA was evaluated using racy of the assembly, meaning that there should be a Nanodrop and low voltage 0.8% agarose TAE gels, and compromise between cost and coverage of the assembly quantified using a Qubit fluorometer. when designing the experiment. Design of guide RNAs targeting conserved sites of MYB10 The Alt-R CRISPR-Cas9 System (IDT®) was selected to genes Variant identification in the de novo assemblies The pipeline used to identify variants between de novo design and assemble the CRISPR RNAs (crRNA) and assemblies succeeded in the identification and assign- the trans-activating crRNAs (tracrRNA) into functional ment of polymorphisms into haplotypes in three of the guide RNAs (gRNAs). samples rather than compiling them, improving the positions with the Alt-R® Custom Cas9 crRNA Design The crRNA sequences were designed in the desired results obtained when using a reference genome (the 44 bp insertion was detected in ‘Fortune’). However, in Tool, which listed the potential crRNAs by their calcu- cases of high similarity between haplotypes or low depth lated on-target activity and sequence position. Those coverage, the use of a reference genome performs bet- with the highest on-target capacity were BLAST to the ter. Therefore, where a reference genome is available, the Peach genome v2.0.a1 [52] and the Sweet cherry genome alignment to it and de novo assemblies should be com- v1.0 [53] to discard gRNAs with off-target risk regions. bined to ensure all polymorphisms are identified. SNP variability identified in other genomes was also con- sidered in the crRNA design. Those gRNAs directed at Strengths, weaknesses, improvements, and opportunities the most highly conserved LG3-MYB10 sequences were We present here an economic and powerful strategy to finally selected. extract sequences and polymorphisms from highly com- plex regions which may contain duplicated and presence/ Assembly of the Cas9 RNA‑Ribonucleoprotein complex absence variants (PAVs). The main strengths and weak- (RNP) ness of this strategy is summarized in Table 3. The designed crRNAs and the tracrRNA were resus- pended to 100 µM in Nuclease-Free Duplex Buffer (IDT). To form the gRNA duplex, all the crRNAs were pooled in Conclusions an equimolar mix and combined with the same volume of The CRISPR-Cas9 enrichment strategy enabled long-read tracrRNA to a final concentration of 10 µM. The mixture sequences from a highly polymorphic and duplicated was denaturalized for five minutes at 95 ºC and allowed region to be produced and polymorphisms related to to cool at room temperature. To assemble the Cas9-RNP fruit skin color to be deduced. The methodology was sim- duplex and 0.08 µL of Alt-R® S.p. HiFi Cas9 Nuclease V3 complex for a single cleavage reaction, 1 µL of the gRNA plified through using a pool of crRNAs in duplicated con- served regions and the sequencing cost was optimized (IDT) were mixed with 1X of CutSmart Buffer (NEB) and by multiplex sequencing while keeping the protocol nuclease-free water to a total of 10 µL then incubated for computationally inexpensive. The tool could be further 30 min at room temperature and kept on ice for immedi- improved by considering methylation patterns and also ate use. adapted for the full assembly of complex regions, being also appropriate for those crops with a lack of a reference Cleavage assay of the designed gRNAs with PCR amplicons genome or with one unrepresentative for the genomic region of interest. The assembled Cas9-RNPs activity was tested using PCR amplicons of PsMYB10.1a, PsMYB10.2 and PsMYB10.3 full gene sequences from ‘Angeleno’. The primer Materials sequences and their annealing temperature (Ta) are Plant material and nucleic acid isolation listed in Additional File 8. Each PCR reaction contained Young leaves from five commercial Japanese plum vari- 1.5 mM ­MgCl2 and 1X ­NH4 buffers, 1U of BioTaq poly- eties (Table 4) were collected and kept at − 80 ºC. Prior merase (Bioline), 0.2 mM dNTPs, 0.2 µM of each primer, to DNA extraction, nuclei were separated as described in 40 ng of DNA from ‘Angeleno’ and MilliQ water to a Naim et al. [49] and detailed online [50], but the vegetal total of 25 µL. Thermocycler conditions were 94 ºC for material was homogenized by grinding in liquid nitrogen 1 min, 35 cycles of 94 ºC for 30 s, the Ta described for using a pestle and mortar. DNA was extracted from the each primer pair for 20 s and 72 ºC for 2 min, followed by isolated nuclei using the Doyle CTAB method [51], intro- a final extension step of 5 min at 72 ºC. Reactions were ducing an RNAse treatment before the chloroform cen- purified using ExoSap-IT (Thermo Fisher) and then 3 µL trifugation step. of each amplicon was diluted with 32 µL of nuclease-free Fiol et al. Plant Methods (2022) 18:105 Page 13 of 16 Table 4 Commercial varieties selected for CRISPR-Cas9 targeted Expansion 1–12 PCR-free kit (EXP-NBD104, ONT) and sequencing of their MYB10 region, their fruit color and the purified using Agencourt AMPure XP beads (Beckman MYB10 haplotype combination described in Fiol et al. [26] Coulter). The barcoded samples were quantified using Variety Fruit color MYB10 a Qubit fluorometer and pooled equimolarly to reach Haplotypes 700 ng in 65 µL. The Cas9 Sequencing Kit (SQK-CS9109, Skin Flesh ONT) was used to ligate the AMX adapters before puri- ‘Angeleno’ Black Yellow H1/H3 fication of the fragments with AMPure XP Beads. The ‘Black Gold’ Black Red H1/H2 beads were washed twice with Long Fragment Buffer ‘Fortune’ Red Yellow H3/H6 (LFB) and then resuspended in Elution Buffer (EB) for ‘Golden Japan’ Yellow Yellow H4/H9 30 min at room temperature. The sequencing library was ‘TC Sun’ Yellow Yellow H4/H5 prepared by mixing 12 µL of the DNA library, 37.5 µL of Sequencing Buffer (SQB) and 25.5 µL of Loading Beads (LB) (SQK-CS9109, ONT). The sample of pooled frag- ments was sequenced in a MinION (v9.4.1) flowcell on water before adding 10 µL of assembled Cas9-RNP. a GridION mK1 device operated by MinKNOW version Cleavage reactions were incubated at 37 ºC for 20 min 3.6.5 software. and five minutes at 72 ºC, then loaded in a 2% agarose TAE gel to visualize the cleaved bands compared to the uncut amplicons. Basecalling, quality filtering, adapter removal, alignment to reference and variant calling Nanopore read sequences were basecalled and demulti- Comparison of the MYB10 region in several Prunus plexed using Guppy (v3.2.10 + aabd4ec) high accuracy whole‑genome assemblies model (ONT). MinIONQC v1.4.2 [57] was used to evalu- A total of fifteen genomes were added to the analysis: two ate the read quality scores before and after trimming bar- Japanese plums, eight apricots, two sweet cherries, two codes and adapters with Porechop v0.2.4 (https://​github.​ peaches and a wild peach relative. For all the genomes, com/​rrwick/​Porec​hop). Reads were aligned to Prunus the MYB10 region was considered from the BLAST salicina ‘Sanyueli’ v1.0 [27] and P. salicina ‘Zhongli No. [54] regions of crRNA-f1 and crRNA-f2 sequences. All 6’ v1.0 [28] whole genome assemblies by using Minimap2 genomes and the regions considered in this study are [58], and were visualized in IGV [59] and represented in listed in Additional File 2. The homology between the bwtool [60]. Alignment information was obtained with regions was calculated and visualized pairwise using the Samtools version 1.11 [61]. Variant calling was done DECIPHER v2.0 [55] package for R version 4.0.3 [56]. with software specifically scripted for third-generation The MYB10 genes in the two Japanese plum genomes sequencing reads: Longshot v0.4.3 [47] was run for SNP were identified by BLAST using the PsMYB10 sequences calling, while both Sniffles v1.0.12 [29] and NanoVar from Fiol et al. [26]. v1.3.9 [30] were run, and their output files merged with VCFtools v0.1.16 [62], to obtain the list of SV calls. VFC- CRISPR‑Cas9 enrichment, library preparation tools was further used to compare the lists of SNP calls and sequencing between the five ‘Sanyueli’ genome alignments, and the The high-quality nuclei DNA from the five selected results represented in a Venn diagram using pyvenn commercial plum varieties and the pooled crRNAs (https://​github.​com/​tctia​nchi/​pyvenn). (100 µM) were sent to CNAG (Centre Nacional d’Anàlisi Genòmica, Barcelona) for CRISPR-Cas9 enrichment tar- Reads correction, Flye assembly and reference‑guided geted sequencing using Oxford Nanopore Technology scaffolding (ONT). For each sample, 1000 ng of DNA diluted in 27 Canu v 2.1.1 [63] was used to correct the trimmed reads µL of nuclease-free water was dephosphorylated by the and perform an initial de novo assembly. Nanopolish addition of 3 µL of 10X CutSmart Buffer and 3 µL of v0.11.1 [64] was used to improve the assembly results. Quick Calf Intestinal Phosphatase (NEB) and then incu- For each sequenced cultivar, their Canu corrected reads bated at 37 ºC for 10 min, 80 ºC for two minutes and and polished contigs were used as input for Flye assem- 20 ºC for two minutes. To each reaction, 10 µL of previ- bler v2.8.3 [65] enabling the option to keep haplotypes. ously assembled Cas9-RNP complex was added together The resulting contigs were genome-guided scaffolded with 1 µL of dATP (10 mM) and 1 µL of Taq polymer- with the RagTag software [66], using separately the ase (NEB) and incubated at 37 ºC for 20 min and at 72 ºC three MYB10 regions from the two P. salicina available for five minutes for A-tailing. Samples were tagged genomes as template reference sequences. DECIPHER using barcodes NB01 to NB05 of the Native Barcoding v2.0 [55] was used to display the homology between the Fiol et al. Plant Methods (2022) 18:105 Page 14 of 16 references and reference-guided sequences in a pair- Acknowledgements We thank Joan París, Ferran Contreras (ADV de Fruita del Baix Llobregat) and wise fashion and obtain their similarity values. Jorge Naranjo (Tany Nature) for facilitating the plant material used in the experiment. We would like to thank David Chagne for thoughtful comments and helpful advice while reviewing our manuscript. Identification of SNPs and InDels in the PsMYB10.1a Author contributions promoter The experiments were conceived and designed by AF and MJA. ELG provided A region of 250 bp in the PsMYB10.1a promoter con- support in the design of the experiment. Experiments were conducted by AF. Bioinformatics analysis were performed by AF and FJR. The paper was taining several SNPs and large InDels was used for written by AF and MJA. All authors have critically revised the manuscript and comparison. Sanger sequences of this region obtained approved the final document. for different haplotypes [26] were used in a BLAST Funding analysis to recover the haplotypes in the de novo (Flye) AF is recipient of grant BES-2016–079060 funded by MCIN/ contigs. In the case of the H4 sequence, which could not AEI/10.13039/501100011033 and by “ESF Investing in your future”. FJR is recipi‑ be recovered, reads were isolated from their alignment ent of grant PRE2019-087427 funded by MCIN/AEI/10.13039/501100011033 and by “ESF Investing in your future. This research was supported by project to the ‘Zhongli No. 6’ genome PsMYB10.1a promoter RTI2018-100795-B-I00 and by project PID2021-128885OB-I00 both funded by region and assembled manually. No Sanger sequences MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. were available for H5 and H6, and their contigs on ‘TC We acknowledge support from the CERCA Programme (“Generalitat de Cata‑ lunya”), and the “Severo Ochoa Programme for Centres of Excellence in R&D” Sun’ and ‘Fortune’ were isolated as alternative to H4 2016–2019 (SEV‐2015‐0533) and 2020–2023 (CEX2019-000902-S) both funded and H3 sequences. The selected contigs from each cul- by MCIN/AEI/10.13039/501100011033. tivar were aligned and compared using Sequencher 5.0 Availability of data and materials (Gene Codes Corporation, Ann Arbor, MI, USA). The raw sequencing data presented in this study has been deposited at the EBI European Nucleotide Archive repository under the accession number Supplementary Information PRJEB48338. The online version contains supplementary material available at https://​doi.​ org/​10.​1186/​s13007-​022-​00937-4. Declarations Additional file 1. Dot plot comparing, pair-wise, the MYB10 regions Ethics approval and consent to participate identified in 15 Prunus genomes, represented as in Figure 1a. The colored Not applicable. squares border the Prunus sections considered: purple for Japanese plums, orange for apricots, red for sweet cherries, and green for peaches Consent for publication and its wild relative. Not applicable. Additional file 2. Details of the MYB10 region isolated from each of the Competing interests fifteen Prunus genomes considered for pair-wise comparison. The results The authors declare that they have no competing interests. of the percentage of sequence hits and the number of shared blocks for every comparison is included. Author details 1 Additional file 3. Details of the crRNAs designed for the Japanese plum Centre for Research in Agricultural Genomics, CSIC-IRTA-UAB-UB, Campus LG3-MYB10 region enrichment. The SNPs identified are underlined and UAB, Barcelona, Spain. 2 The New Zealand Institute for Plant and Food Research crRNAs were designed including each variant. Limited (Plant & Food Research), Private Bag 11600, Palmerston North 4442, New Zealand. 3 Institut de Recerca I Tecnologia Agroalimentàries, Barcelona, Additional file 4. Visualization of the depth of the sequences aligned to Spain. the ‘Sanyueli’, Zhongli-1 and Zhongli-2 regions. The pink bars highlight the coordinates with previously identified MYB10 gene sequences, where Received: 25 November 2021 Accepted: 17 August 2022 the Cas9 enzyme has two cutting points enabling sequencing in both directions. Additional file 5. The total number of SNPs called from each sequenced sample and the count of the shared SNP positions between them. Additional file 6. Pair-wise visualization of the homologous hits (above References diagonal) and homologous blocks (below) between each reference region 1. Taylor JS, Raes J. Duplication and divergence: the evolution of new genes and the de novo contigs scaffolded using the ‘Sanyueli’, Zhongli-1 or and old ideas. Annu Rev Genet. 2004;38:615–43. Zhongli-2 regions, represented in the same colors as in Figure 1a. The light 2. Van de Peer Y, Maere S, Meyer A. The evolutionary significance of ancient pink lines show the position of the MYB10 genes in each of the reference genome duplications. Nat Rev Genet. 2009;10(10):725–32. regions used for the scaffolding. 3. Gabur I, Chawla HS, Snowdon RJ, Parkin IA. Connecting genome structural variation with complex traits in crop plants. Theor Appl Genet. Additional file 7. Visualization of two phased variants 1 kb apart on H9. 2019;132(3):733–50. The 44 bp insertion is present in H1, H3 and H9 and was associated to the 4. Soltis DE, Visger CJ, Soltis PS. The polyploidy revolution then…and now: red skin color. The polymorphism is phased with an 8 bp insertion at the Stebbins revisited. Am J Bot. 2014;101(7):1057–78. start of exon 1 of the PsMYB10.1 gene, which explains its lack of function 5. Giannuzzi G, D’Addabbo P, Gasparro M, Martinelli M, Carelli FN, Antonacci on H9. The reads from H4 do not show either of the two polymorphisms. D, et al. Analysis of high-identity segmental duplications in the grapevine Additional file 8. Primer sequences and their annealing temperatures (Ta) genome. BMC Genom. 2011;12(1):1–14. used to PCR amplify the MYB10 gene sequences. 6. Fares MA, Keane OM, Toft C, Carretero-Paulet L, Jones GW. The roles of whole-genome and small-scale duplications in the functional specializa‑ tion of Saccharomyces cerevisiae genes. Plos Genet. 2013;9(1): e1003176. Fiol et al. Plant Methods (2022) 18:105 Page 15 of 16 7. Zmienko A, Marszalek-Zenczak M, Wojciechowski P, Samelak-Czajka A, variants using low-depth nanopore sequencing. Genome Biol. Luczak M, Kozlowski P, et al. AthCNV: a map of DNA copy number varia‑ 2020;21(1):1–15. tions in the Arabidopsis genome. Plant Cell. 2020;32(6):1797–819. 31. Panchy N, Lehti-Shiu M, Shiu S-H. Evolution of gene duplication in 8. Vandepoele K, Simillion C, Van de Peer Y. Evidence that rice and other plants. Plant Physiol. 2016;171(4):2294–316. cereals are ancient aneuploids. Plant Cell. 2003;15(9):2192–202. 32. González VM, Aventín N, Centeno E, Puigdomènech P. High presence/ 9. Lye ZN, Purugganan MD. Copy number variation in domestication. Trends absence gene variability in defense-related gene clusters of Cucumis Plant Sci. 2019;24(4):352–65. melo. BMC Genom. 2013;14(1):1–13. 10. Xu C, Nadon BD, Kim KD, Jackson SA. Genetic and epigenetic diver‑ 33. Qiao X, Yin H, Li L, Wang R, Wu J, Wu J, et al. Different modes of gene gence of duplicate genes in two legume species. Plant Cell Environ. duplication show divergent evolutionary patterns and contribute 2018;41(9):2033–44. differently to the expansion of gene families involved in important fruit 11. Gu YQ, Crossman C, Kong X, Luo M, You FM, Coleman-Derr D, et al. traits in pear (Pyrus bretschneideri). Front Plant Sci. 2018;9:161. Genomic organization of the complex α-gliadin gene loci in wheat. Theor 34. Gu C, Wang L, Wang W, Zhou H, Ma B, Zheng H, et al. Copy num‑ Appl Genet. 2004;109(3):648–57. ber variation of a gene cluster encoding endopolygalacturonase 12. Camerlengo F, Sestili F, Silvestri M, Colaprico G, Margiotta B, Ruggeri R, mediates flesh texture and stone adhesion in peach. J Exp Bot. et al. Production and molecular characterization of bread wheat lines 2016;67(6):1993–2005. with reduced amount of α-type gliadins. BMC Plant Biol. 2017;17(1):1–11. 35. Barragan AC, Weigel D. Plant NLR diversity: the known unknowns of 13. Ruiz M, Giraldo P, Royo C, Villegas D, Aranzana MJ, Carrillo JM. Diversity Pan-NLRomes. Plant Cell. 2021;33(4):814–31. and genetic structure of a collection of Spanish durum wheat landraces. 36. Badet T, Croll D. The rise and fall of genes: origins and functions of Crop Sci. 2012;52(5):2262–75. plant pathogen pangenomes. Curr Opin Plant Biol. 2020;56:65–73. 14. Van Herpen TW, Goryunova SV, Van Der Schoot J, Mitreva M, Salentijn 37. Lallemand T, Leduc M, Landès C, Rizzon C, Lerat E. An overview of E, Vorst O, et al. Alpha-gliadin genes from the A, B, and D genomes of duplicated gene detection methods: why the duplication mechanism wheat contain different sets of celiac disease epitopes. BMC Genom. has to be accounted for in their choice. Genes. 2020;11(9):1046. 2006;7(1):1–13. 38. Gabrieli T, Sharim H, Fridman D, Arbib N, Michaeli Y, Ebenstein Y. 15. Chawla HS, Lee H, Gabur I, Vollrath P, Tamilselvan-Nattar-Amutha S, Selective nanopore sequencing of human BRCA1 by Cas9-assisted Obermeier C, et al. Long-read sequencing reveals widespread intragenic targeting of chromosome segments (CATCH). Nucleic Acids Res. structural variants in a recent allopolyploid crop plant. Plant Biotechnol J. 2018;46(14):e87-e. 2021;19(2):240–50. 39. Bruijnesteijn J, van der Wiel M, de Groot NG, Bontrop RE. Rapid char‑ 16. Alonge M, Wang X, Benoit M, Soyk S, Pereira L, Zhang L, et al. Major acterization of complex genomic regions using Cas9 enrichment and impacts of widespread structural variation on gene expression and crop Nanopore sequencing. bioRxiv. 2021. improvement in tomato. Cell. 2020;182(1):145-61. e23. 40. McDonald TL, Zhou W, Castro CP, Mumm C, Switzenberg JA, Mills RE, et al. 17. Della Coletta R, Qiu Y, Ou S, Hufford MB, Hirsch CN. How the pan- Cas9 targeted enrichment of mobile elements using nanopore sequenc‑ genome is changing crop genomics and improvement. Genome Biol. ing. Nat Commun. 2021;12(1):1–13. 2021;22(1):1–19. 41. Kirov I, Merkulov P, Gvaramiya S, Komakhin R, Omarov M, Dudnikov M, 18. Tao Y, Zhao X, Mace E, Henry R, Jordan D. Exploring and exploiting pan- et al. Illuminating the transposon insertion landscape in plants using genomics for crop improvement. Mol Plant. 2019;12(2):156–69. Cas9-targeted Nanopore sequencing and a novel pipeline. bioRxiv. 2021. 19. She X, Jiang Z, Clark RA, Liu G, Cheng Z, Tuzun E, et al. Shotgun sequence 42. Jung S, Jiwan D, Cho I, Lee T, Abbot AG, Sosinski B, et al. Synteny of Prunus assembly and recent segmental duplications within the human genome. and other model plant species. BMC Genom. 2009;10:76. Cah Rev The. 2004;431(7011):927–30. 43. Stracke R, Werber M, Weisshaar B. The R2R3-MYB gene family in Arabidop‑ 20. Amarasinghe SL, Su S, Dong X, Zappia L, Ritchie ME, Gouil Q. Opportuni‑ sis thaliana. Curr Opin Plant Biol. 2001;4(5):447–56. ties and challenges in long-read sequencing data analysis. Genome Biol. 44. Sternberg SH, Redding S, Jinek M, Greene EC, Doudna JA. DNA inter‑ 2020;21(1):1–16. rogation by the CRISPR RNA-guided endonuclease Cas9. Cah Rev The. 21. Liu Y, Du H, Li P, Shen Y, Peng H, Liu S, et al. Pan-genome of wild and 2014;507(7490):62–7. cultivated soybeans. Cell. 2020;182(1):162-76. e13. 45. Brinkman EK, Chen T, de Haas M, Holland HA, Akhtar W, van Steensel B. 22. Song J-M, Guan Z, Hu J, Guo C, Yang Z, Wang S, et al. Eight high-quality Kinetics and fidelity of the repair of Cas9-induced double-strand DNA genomes reveal pan-genome architecture and ecotype differentiation of breaks. Mol Cell. 2018;70(5):801-13.e6. Brassica napus. Nat Plants. 2020;6(1):34–45. 46. Li W, Freudenberg J. Mappability and read length. Front Genet. 23. Gilpatrick T, Lee I, Graham JE, Raimondeau E, Bowen R, Heron A, et al. 2014;5:381. Targeted nanopore sequencing with Cas9-guided adapter ligation. Nat 47. Edge P, Bansal V. Longshot enables accurate variant calling in diploid Biotechnol. 2020;38(4):433–8. genomes from single-molecule long read sequencing. Nat Commun. 24. López-Girona E, Davy MW, Albert NW, Hilario E, Smart ME, Kirk C, et al. 2019;10(1):1–10. CRISPR-Cas9 enrichment and long read sequencing for fine mapping in 48. Wick RR, Judd LM, Holt KE. Performance of neural network basecalling plants. Plant Methods. 2020;16(1):1–13. tools for Oxford Nanopore sequencing. Genome Biol. 2019;20(1):1–10. 25. Giesselmann P, Brändl B, Raimondeau E, Bowen R, Rohrandt C, Tandon R, 49. Naim F, Nakasugi K, Crowhurst RN, Hilario E, Zwart AB, Hellens RP, et al. et al. Analysis of short tandem repeat expansions and their methylation Advanced engineering of lipid metabolism in Nicotiana benthamiana state with nanopore sequencing. Nat Biotechnol. 2019;37(12):1478–81. using a draft genome and the V2 viral silencing-suppressor protein. PLoS 26. Fiol A, García-Gómez BE, Jurado-Ruiz F, Alexiou K, Howad W, Aranzana ONE. 2012;7(12): e52717. MJ. Characterization of Japanese Plum (Prunus salicina) PsMYB10 alleles 50. Hilario E. Plant nuclear genomic DNA preps 2018 [Available from: https://​ reveals structural variation and polymorphisms correlating with fruit skin www.​proto​cols.​io/​view/​plant-​nucle​ar-​genom​ic-​dna-​preps-​rncd5​aw. color. Front Plant Sci. 2021;12:1057. 51. Doyle JJ, Doyle JL. A rapid DNA isolation procedure for small quantities of 27. Liu C, Feng C, Peng W, Hao J, Wang J, Pan J, et al. Chromosome- fresh leaf tissue. 1987. level draft genome of a diploid plum (Prunus salicina). GigaScience. 52. Verde I, Jenkins J, Dondini L, Micali S, Pagliarani G, Vendramin E, et al. 2020;9(12):giaa130. The Peach v2.0 release: high-resolution linkage mapping and deep 28. Huang Z, Shen F, Chen Y, Cao K, Wang L. Chromosome-scale genome resequencing improve chromosome-scale assembly and contiguity. BMC assembly and population genomics provide insights into the adapta‑ Genom. 2017;18(1):225. tion, domestication, and flavonoid metabolism of Chinese plum. Plant J. 53. Shirasawa K, Isuzugawa K, Ikenaga M, Saito Y, Yamamoto T, Hirakawa H, 2021;108(4):1174–92. et al. The genome sequence of sweet cherry (Prunus avium) for use in 29. Sedlazeck FJ, Rescheneder P, Smolka M, Fang H, Nattestad M, Von Haese‑ genomics-assisted breeding. DNA Res. 2017;24(5):499–508. ler A, et al. Accurate detection of complex structural variations using 54. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment single-molecule sequencing. Nat Methods. 2018;15(6):461–8. search tool. J Mol Biol. 1990;215(3):403–10. 30. Tham CY, Tirado-Magallanes R, Goh Y, Fullwood MJ, Koh BT, Wang W, 55. Wright ES. Using DECIPHER v2.0 to analyze big biological sequence data et al. NanoVar: accurate characterization of patients’ genomic structural in R. R J. 2016;8(1):352. Fiol et al. Plant Methods (2022) 18:105 Page 16 of 16 56. R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing.Vienna, Austria. 2020. http://​www.R-​ proje​ct.​org/. 57. Lanfear R, Schalamun M, Kainer D, Wang W, Schwessinger B. MinIONQC: fast and simple quality control for MinION sequencing data. Bioinformat‑ ics. 2019;35(3):523–5. 58. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinfor‑ matics. 2018;34(18):3094–100. 59. Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6. 60. Pohl A, Beato M. bwtool: a tool for bigWig files. Bioinformatics. 2014;30(11):1618–9. 61. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9. 62. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8. 63. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Biotechfor. 2017;27(5):722–36. 64. Loman NJ, Quick J, Simpson JT. A complete bacterial genome assem‑ bled de novo using only nanopore sequencing data. Nat Methods. 2015;12(8):733–5. 65. Kolmogorov M, Yuan J, Lin Y, Pevzner PA. Assembly of long, error-prone reads using repeat graphs. Nat Biotechnol. 2019;37(5):540–6. 66. Alonge M, Soyk S, Ramakrishnan S, Wang X, Goodwin S, Sedlazeck FJ, et al. RaGOO: fast and accurate reference-guided scaffolding of draft genomes. Genome Biol. 2019;20(1):1–17. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub‑ lished maps and institutional affiliations. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions