The Plant Journal (2020) 101, 455–472 doi: 10.1111/tpj.14538 RESOURCE Transposons played a major role in the diversification between the closely related almond and peach genomes: results from the almond genome sequence Tyler Alioto1,2, Konstantinos G. Alexiou3,4, Amelie Bardil4, Fabio Barteri4, Rau  l Castanera4, Fernando Cruz1,2, Amit Dhingra5, 6  Henri Duval , Angel Fernandez i Martı , Leonor Frias , Beatriz Galan , Jose L. Garcıa9, Werner Howad3,4, 7,8 1,2 9 Jessica Go  mez-Garrido1,2, Marta Gut1,2, Irene Julca2,10, Jordi Morata4, Pere Puigdomenech4, Paolo Ribeca1,2,11, Marıa J. Rubio Cabetas12,13, Anna Vlasova10, Michelle Wirthensohn14, Jordi Garcia-Mas3,4, Toni Gabaldo  n2,10,15, 4, Josep M. Casacuberta * and Pere Aru s 3,4 1 CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, 08028 Barcelona, Spain, 2 Universitat Pompeu Fabra (UPF), 08005 Barcelona, Spain, 3 IRTA, Campus UAB, Edifici CRAG, Cerdanyola del Valle s (Bellaterra), 08193 Barcelona, Spain, 4 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Campus UAB, Edifici CRAG, Cerdanyola del Valle s (Bellaterra), 08193 Barcelona, Spain, 5 Department of Horticulture, Washington State University, 99164-6414 Pullman, WA, USA, 6 INRA, UR1052, Unite de Ge ne tique et Ame lioration des Fruits et Le gumes (GAFL), Domaine St. Maurice CS 60094, 84143, Montfavet Cedex, France, 7 Department of Environmental Science Policy and Management, University of California, Berkeley, 94720, CA, USA, 8 Innovative Genomics Institute (IGI), 94720 Berkeley, CA, USA, 9 Department of Environmental Biology, Center for Biological Research (CIB-CSIC), Spanish National Research Council (CSIC), Ramiro de Maeztu 9, 28040 Madrid, Spain, 10 Bioinformatics and Genomics Programme, Centre for Genomic Regulation (CRG), Dr Aiguader, 88, 08003 Barcelona, Spain, 11 The Pirbright Institute, Woking, Surrey, GU24 0NF, UK, 12 Centro de Investigacio n y Tecnologıa Agroalimentaria de Arago n (CITA), Unidad de Hortofruticultura, Gobierno de Arago n, Avda. Montan~ ana 930, 50059 Zaragoza, Spain, 13 Instituto Agroalimentario de Arago n – IA2 (CITA-Universidad de Zaragoza), Calle Miguel Servet 177, 50013 Zaragoza, Spain, 14 University of Adelaide, Waite Research Institute, School of Agriculture, Food and Wine, PMB 1, Glen Osmond, SA, 5064, Australia, and 15 Institucio Catalana de Recerca i Estudis Avancß ats (ICREA), Pg Lluıs Companys 23, 08010 Barcelona, Spain Received 11 June 2019; revised 29 August 2019; accepted 2 September 2019; published online 16 September 2019. *For correspondence (e-mail pere.arus@irta.cat and josepcasacuberta@cragenomica.es). SUMMARY We sequenced the genome of the highly heterozygous almond Prunus dulcis cv. Texas combining short- and long-read sequencing. We obtained a genome assembly totaling 227.6 Mb of the estimated almond genome size of 238 Mb, of which 91% is anchored to eight pseudomolecules corresponding to its haploid chromosome complement, and annotated 27 969 protein-coding genes and 6747 non-coding transcripts. By phylogenomic comparison with the genomes of 16 additional close and distant species we estimated that almond and peach (Prunus persica) diverged around 5.88 million years ago. These two genomes are highly syntenic and show a high degree of sequence conservation (20 nucleotide substitutions per kb). However, they also exhibit a high number of presence/absence variants, many attributable to the movement of trans- posable elements (TEs). Transposable elements have generated an important number of presence/absence variants between almond and peach, and we show that the recent history of TE movement seems markedly different between them. Transposable elements may also be at the origin of important phenotypic differ- ences between both species, and in particular for the sweet kernel phenotype, a key agronomic and domes- tication character for almond. Here we show that in sweet almond cultivars, highly methylated TE © 2019 The Authors. 455 The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. 456 Tyler Alioto et al. insertions surround a gene involved in the biosynthesis of amygdalin, whose reduced expression has been correlated with the sweet almond phenotype. Altogether, our results suggest a key role of TEs in the recent history and diversification of almond and its close relative peach. Keywords: Prunus dulcis, Prunus persica, genome sequence, variability, divergence, indels, transposable elements, crop evolution, seed bitterness. INTRODUCTION The genome sequences of some Prunus species are available, including the high-quality genome of peach Almond, Prunus dulcis (Miller) D.A. Webb (syn. Prunus (Verde et al., 2017) and those of sweet cherry (Prunus amygdalus Batsch), is a rosaceous tree species cultivated avium L.) (Shirasawa et al., 2017), mume (Prunus mume for its seeds; it has a diploid (2n = 2x = 16) and compact L.), a relative of apricot (Zhang et al., 2017), and Prunus genome (about 300 Mbp) (Baird et al., 1994). The genus yedoensis, a wild cherry tree (Baek et al., 2018), the latter Prunus comprises a group of approximately 200 species, two used for ornamental purposes. The genome sequence some of which, such as the stone fruits (peach, apricot, of almond cv. Lauranne has been very recently added cherry and plum) and almond, have high economic value (Sanchez-Perez et al., 2019). In this paper, we present the (Aranzana et al., 2019). The high level of genomic resem- whole genome sequence of almond cv. Texas, a self-in- blance and synteny among the species of this genus (Dir- compatible and highly heterozygous genotype that was lewanger et al., 2004) enables production of hybrids that obtained in the USA from materials imported from western are sometimes fertile. Europe. Texas (also called Texas Prolific and Mission) was Humans used almond as a food long before the advent bred at Houston, Texas, USA, as a seedling of French culti- of agriculture, and the oldest remains have been found in var Languedoc (Wickson, 2001), and became one of the Israel, dating from 19 000 years ago (Kislev et al., 1992), leading cultivars in California in the twentieth century although its domestication probably occurred 14 000 years along with ‘Nonpareil’ (Kester et al., 2015). Texas was also later (Spiegel-Roy, 1976). The origin of the almond tree is one of the parents, the other was peach cv. Earlygold, of not well established; its closest wild relatives live in central the interspecific progeny used for the construction of the and western Asia, stretching from the Himalayas to the reference linkage map of Prunus (Joobeur et al., 1998). We eastern Mediterranean Basin (Yazbek and Al-Zein, 2009). also compared the Texas almond genome sequence with Based on the distribution of the cultivated species, two other sequenced genomes, including that of its close rela- alternative hypotheses place the domestication site of tive peach, and found that, in addition to other aspects of almond in the Levant (Browicz and Zohary, 1996) or in cen- diversity between these genomes already reported (Yu tral Asia (Ladizinsky, 1999). Diamond (1997) proposed et al., 2001; Velasco et al., 2014), transposable elements almond as an example of simple domestication, where a (TEs) played a key role in their recent diversification. dominant mutation at a single gene conferring a sweet taste to the otherwise bitter and toxic kernel would result in an edible and cultivable crop. This gene, sweet kernel RESULTS Sk/sk, was initially described by Heppner (1923) and later Texas almond sequence assembly, annotation and mapped to the central region of chromosome 5 (Sanchez- comparison with the linkage map and the peach sequence Perez et al., 2007); it has recently been proposed to corre- spond to a basic helix–loop–helix (bHLH) transcription fac- A total of 138.6 Gb of Illumina (>5009 coverage) and tor (Sanchez-Perez et al., 2019). The closest relatives of 10.2 Gb (379) of Oxford Nanopore Technologies (ONT; almond are within the subgenus Amygdalus, encompass- score ≥ 7.0 and a read N50 length of 7.3 kb) sequence were ing peach [Prunus persica (L.) Batsch] and a group of 25 produced (Table S1 in the online Supporting Information). wild species (Yazbek and Al-Zein, 2009). Peach and almond By analyzing k-mer frequency, the lower bound for gen- hybrids are fertile. In fact, peach was proposed by Darwin ome size was estimated to be 238 Mb (Figure S1). We col- (1968) as a possible direct derivative of almond with a fle- lapsed the assembly into a haploid representation and shy, non-dehiscent and juicy mesocarp. However, molecu- anchored it to eight pseudomolecules, the number of the lar phylogenetics has identified a clear separation between haploid almond chromosome complement. The final peach and almond consistent with their geographical ori- assembly, P. dulcis Texas v.2.0 (also known as pdulcis26) gin and distribution: peach and its closest relatives are totals 227.6 Mb (91.5% of which is anchored to the eight native to China and eastern Asia, whereas almond and its pseudomolecules) and has a contig and scaffold N50s of wild relatives are native to central and western Asia (Del- 103.9 and 381.5 kb, respectively (Table 1). The complete- plancke et al., 2016). ness of the assembly as determined by BUSCO analysis is © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 457 96.4%: 95.4% complete (89.4% unique and 6.0% dupli- usually coincide with pericentromeric regions and occurred cated), 1.0% fragmented and 3.6% missing BUSCOs. k-mer at similar regions to those in the peach genome (Verde analysis confirmed the BUSCO results (Figure S1). et al., 2017). We annotated a total of 27 969 protein-coding genes that Phylogenomic analysis produce 34 039 transcripts (1.22 transcripts per gene) and encode for 32 559 unique protein products (Table 1). We To shed light on the evolutionary history of the genome of were able to assign some type of functional annotation to P. dulcis in the context of 16 other sequenced plant species 92% of them. In addition, we annotated 6747 non-coding (Table S4), we generated the phylomes of almond and transcripts, of which 3590 and 3153 are long and short peach, that is their complete collections of gene phyloge- non-coding RNA genes, respectively. Most of the main nies (see Experimental Procedures). These phylomes were assembly and annotation parameters of the Texas genome filtered to remove the gene trees containing proteins with presented here were similar to those obtained by Sanchez- domains associated with transposons. After filtering, a Perez et al. (2019) on cv. Lauranne (see Table S2). total of 18 475 and 20 812 gene trees were kept for almond Out of the 1833 single nucleotide polymorphisms (SNPs) and peach, respectively. These filtered phylomes were that comprise the Texas 9 Earlygold (T9E) linkage map scanned to infer duplications and speciation events and (Donoso et al., 2015), 1609 (87.8%) mapped onto the derive orthology and paralogy relationships from individ- almond assembly with single high-quality hits (percentage ual gene trees (Gabaldo  n, 2008). These analyses produced identity ≥ 90% and marker coverage ≥ 90%) (Table S3), a catalog of gene duplication events and phylogeny-based with 1597 (93.4%) mapping onto the anchored assembly. homology relationships for genes in the 17 considered Of the anchored SNPs, 1578 aligned with the pseudo- species, which were used in subsequent analyses. molecules and were syntenic and collinear to the eight We concatenated the protein alignments of 262 genes chromosomes of peach. Only 19 SNPs had a different that had single-copy orthologs in all the 17 species consid- order on the assembly, two of which mapped to a pseudo- ered to reconstruct a phylogeny of these species. The molecule that was different from the linkage group of the resulting highly supported topology (Figure 1a) was con- T9E map and 10 on an unassigned scaffold (Table S3, Fig- gruent with current views on plant phylogeny (Shaw and ure S2), which we attributed to contig reordering or minor Small, 2004) and results in P. dulcis and P. persica forming misassemblies. Similarly, we observed high synteny and a clade, to the exclusion of P. mume and P. avium (Bade- collinearity between the genome sequence of peach v.2.0 nes and Parfitt, 1995; Scholz et al., 2013). The same topol- a1 (Verde et al., 2002) and that of the Texas almond gener- ogy was obtained when all individual gene trees were ated here (Figure S3). A comparison of physical versus combined into a single species phylogeny by using a gene genetic distances of the eight pseudomolecules is pre- tree parsimony approach. We estimated the divergence sented in Figure S4. Regions of low recombination rates times on this topology using a Bayesian relaxed molecular clock approach. According to our results, P. dulcis diverged from P. persica approximately 5.88 million years Table 1 Texas genome assembly and annotation statistics ago (Mya), from P. mume 20.84 Mya and from P. avium Assembly length 227.6 Mb 62.04 Mya (Figure 1a, Table S5). Contig N50 103.9 kb We then calculated the duplication frequency (i.e. the Scaffold N50 381.5 kb Pseudomolecule N50 24.8 Mb average number of duplications per gene) for each node of Per cent anchored to pseudomolecules 91.47% the species tree, and observed a slightly high duplication BUSCO complete genes 95.4% frequency (about 0.20; Figure 1b) at the common ancestor BUSCO fragmented genes 1.0% of all Prunus species, which is supported by both almond BUSCO missing genes 3.6% and peach phylomes. A functional analysis of protein fami- Genomic GC content 37.65% Number of protein-coding genes 27 969 lies duplicated at this branch shows enrichment of some Median gene length (bp) 2288 molecular functions such as methyltransferase activity, Number of transcripts 34 039 ionotropic glutamate receptor activity, terpene synthase Number of unique protein products 32 559 activity, oxidoreductase activity and transferase activity. In Number of exons 184 149 addition, some biological processes were enriched: Number of unique exons 148 374 Number of coding exons 140 538 response to auxin, metabolic process and oxidation–reduc- Coding GC content 44.12% tion process (Table S6). Median intron length (bp) 171 Then we focused on duplications specific to almond or Exons/transcript 5.41 peach, including large expansions. A total of 1175 (4.4%) Transcripts/gene 1.22 almond proteins and 831 (3.1%) peach proteins have an in- Multi-exonic transcripts 81% Gene density (genes Mb–1) 123 paralog (a recent paralog resulting from a duplication that specifically occurred in the almond and peach lineage, © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 458 Tyler Alioto et al. they were present in the same or different scaffolds, and if present in the same scaffold we counted the number of genes that are present between the resulting in-paralogs. From the total number of species-specific paralogous pairs or in-paralogs (732), the majority (384, 52%) were located in the same scaffold, and 285 pairs (39%) were located in close proximity, with fewer than 10 intervening genes. A group of 87 (12%) paralogs were directly neighboring each other, and most of them had functions associated with stress responses. We next analyzed protein gains and losses in the lin- eages leading to almond and peach, as inferred from the analysis of the gene trees in the phylome. When we ana- lyzed the almond phylome, a total of 1471 proteins were gained in almond and 1146 were lost in peach. For the peach phylome we found that 1984 proteins were gained in peach and 3157 were lost in almond. Functional analysis shows that the proteins lost in almond with respect to peach are enriched in functions related to serine-type endopeptidase inhibitor activity, nutrient reservoir activity, lipid transport, response to auxin, oxidation–reduction pro- cess and ion transport. Conversely, genes lost in peach with respect to almond are mainly enriched in functions related to transferase activity, transcription and ATP syn- thesis-coupled proton transport (Table S7). Variability of almond cultivars Read mapping rate, depth and genome coverage. Align- Figure 1. Species tree obtained from the concatenation of 262 widespread ment of the 919 019 814 trimmed reads from the 10 re-se- single-gene families. quenced almond cultivars (Table S8) to our reference (a) Full species tree. All Prunus species are highlighted in pink. All bootstrap assembly resulted in mapping of 825 914 441 ‘clean’ reads values that are not maximal (bootstrap 100%) are indicated in red. Green numbers correspond to the nodes in Table S4. Bars at the nodes indicate (after removal of unmapped reads, PCR duplicates and the uncertainty around mean age estimates based on 95% credibility inter- reads with a mapping quality < 10) which corresponds to vals. Scale at the bottom shows the divergence time in Mya (million years an average mapping rate of 89.6% (Table S9). Nonpareil ago). Green dots represent selected calibration points. (b) Zoom-in of the Prunus group. Numbers indicate the duplication ratio for and Vivot, and Falsa Barese and Genco showed the highest each branch calculated with the phylome of almond (red) and peach (94%) and lowest (79%) mapping rate, respectively (blue).[Colour figure can be viewed at wileyonlinelibrary.com] (Table S9). Regarding sample depth and genome coverage, an average read depth of 39.7 was detected for the 10 culti- respectively). These paralogs could be assigned to 542 vars, whereas 96.7% of the assembly was covered by the almond-specific gene expansions and 367 peach-specific re-sequencing data on average. Marcona and Falsa Barese gene expansions. In both almond and peach most expan- had the highest (51.49) and lowest (27.79) read depth, sions [540 (99.6%) for almond; 363 (98.9%) for peach] have respectively (Table S9). a moderate size (two to five in-paralogs; Figure S5). Some almond expansions of size two encode putative members Variant calling and phylogenetic analysis. Genetic vari- of the lignin biosynthesis pathway (Vanholme et al., 2010) ability analysis resulted in the detection of 2 253 377 vari- such as caffeic acid 3-O-methyltransferase (COMT; Prudu- ants, of which 2 203 582 (87%) were SNPs and 330 795 l26A009858P1-Prudul26A011895P1) and shikimate O-hy- (13%) were insertions/deletions (indels). Genome-wide dis- droxycinnamoyl transferase (HCT; Prudul26A003924P1- tribution of SNPs and indels can be seen in Table S10. Prudul26A028947P1, Prudul26A000852P1-Prudul26A0228 Nonpareil had the highest number of SNPs and indels, 43P2). Interestingly, these genes have undergone parallel 1 072 759 and 142 142, respectively, whereas Ripon had duplications in P. persica, P. mume and P. avium that the lowest number of SNPs and indels, 827 397 and occurred independently in the three species. In order to 94 070, respectively. Average SNP density was calculated check whether the almond species-specific duplications are as 6.2 SNPs per kb, whereas the average heterozygosity tandem or dispersed duplications we assessed whether for the 10 cultivars was 0.44% (Table S11). Our calculations © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 459 of SNP density and heterozygosity for almond are lower accounting for 18 Mb of sequence, equivalent to 8% of the than those recently published (a SNP density of 19.1 and almond genome (Table S13). heterozygosity of 0.69; Yu et al., 2001). This discrepancy Re-sequencing data from the peach cv. Earlygold were could be attributable to the use of the peach genome as a compared with the almond reference genome and dele- reference sequence during the variant calling in the pub- tions were identified, as previously done with almond culti- lished study, a varietal panel of larger genetic diversity or var re-sequencing data (Table S12). The average number different filters and tools used for variant calling. of deletions in almond cultivars compared with the almond A graphical representation of SNP and indel distribution reference sequence was 62 238, whereas in Earlygold this in windows of 100 kbp showed a similar profile for most of figure was more than double (126 137) (Table S12). How- the almond cultivars analyzed (Figure S6). Nevertheless, ever, when considering only deletions larger than 50 bp, lower overall variant density was observed in certain culti- peach had almost 12 times more (1436 versus 120) than vars such as Genco, Falsa Barese and Ripon, which we almond. attributed to the lower number of reads mapped in these Transposable element landscape genotypes. An analysis of deletions in the collection of almond culti- Using the REPET pipeline we annotated 38.21% of the vars is presented in Table S12. For small deletions (1–50 bp), almond genome as TE-related sequences (Tables S14 and numbers were about half of those estimated for indels in S15). The distribution of TEs along almond pseudochromo- Table S10 (about 60 000 versus about 120 000 per cultivar), somes shows an inverse correlation with respect to the as expected considering that only deletions and not inser- gene density, with TE-rich regions showing low gene den- tions are considered. On average, 27% of these deletions sity per chromosome, coinciding with pericentromeric overlapped with TEs. Considering large deletions (>50 bp), regions, and lower TE densities in the gene-rich chromoso- we detected 1219 unique events, the majority of which fell mal arms (Figure 2a). The almond TE landscape was com- within the range 51–500 bp, with an average number of 88 pared with that of peach. For that purpose, we annotated deletions per cultivar. Marcona had the largest number of peach TEs with the same strategy and found very similar deletions (219) whereas Ripon had the lowest number (12) results: 37.60% of TE content (Tables S14 and S15) and a (Table S12). Five hundred and eighty-eight large deletions comparable TE and gene distribution to that of almond (48.2% of the total) overlapped with TEs, with Nonpareil and chromosomes (Figure 2b). Ripon showing nearly 60% of overlap with TEs and Cristo- In addition to the general TE annotation, we performed morto with the lowest percentage of overlap (38.8%). For a dedicated annotation of the long terminal repeat (LTR) deletions larger than 500 bp, almost all the events were retrotransposons and miniature inverted-repeat transpos- found to overlap with TEs (92.6% in the range 501–10 000 able elements (MITEs) in the almond and peach genomes. and 100% in the range 10 001–50 000; Table S12, Figure S7). A conservative search for LTR retrotransposons with a A SNP-based phylogenetic analysis grouped the almond well-preserved structure (i.e. presence of LTRs and coding cultivars into two main clades where the first clade con- capacity for retrotransposon-related proteins) resulted in tained Cristomorto, Falsa Barese and Genco while the sec- the annotation of approximately 2200 elements in both ond clade is split into two subclades, the first containing A€ı, almond and peach (Table S16). Whenever possible, these Belle d’Aurons, Nonpareil and Ripon and the second Des- elements were classified as Copia or Gypsy (i.e. when cod- mayo largueta, Marcona and Vivot (Figure S8). This phy- ing regions for both integrase and reverse transcriptase logeny is in agreement with the geographical origin of the were detected, which allowed us to classify them) or analyzed cultivars, grouping the Italian cultivars, the French remained as unclassified LTR retrotransposons. Although and US cultivars and finally the Spanish cultivars in the same the content of these elements was similar in both gen- clade. The fact that French and US cultivars are clustered omes, the number of LTR retrotransposons that remained together agrees with the known origin of US materials com- unclassified in almond was slightly higher. The distribution ing from French imported accessions (Kester et al., 2015). of the almond and peach retrotransposons along chromo- somes is also highly similar, with Gypsy elements showing Indel variants between peach and almond and their rela- a tendency to concentrate in a central region of chromo- tionship to TE sequences. To assess the structural vari- somes, probably coinciding with the centromeric regions, ability between almond and peach genomes we aligned whereas Copia elements are more evenly distributed (Fig- almond genome contigs to the peach reference genome. A ure 2). A conservative search for MITEs with well-pre- total of 92.96% of the Texas almond assembly could be served Toll-interleukin like regions (TIRs) rendered 10 460 aligned to the Lovell reference peach genome sequence MITEs in almond and 8738 MITEs in peach (Table S16). with an average identity of 95.59%, which increased to The distribution of these elements along chromosomes is 97.99% (20 SNPs per kb) when only regions that align 1:1 similar in peach and almond and follows that of Copia LTR are considered. We detected a total of 20 418 indel variants retrotransposons (Figure 2). © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 460 Tyler Alioto et al. Figure 2. Distribution of gene and transposable element (TE) abundance. Distribution of gene and TE abundance along Prunus dulcis (a) and Prunus persica (b) chromosomes. Outer to inner tracks represent the coverage per 100 kb of genes, TEs, Copia long terminal repeat (LTR) retrotransposons, Gypsy LTR retrotransposons, and miniature inverted-repeat transposable elements. The chromo- some scale is in Mbp. [Colour figure can be viewed at wileyonlinelibrary.com] The LTR retrotransposon dynamics in almond and peach LTR retrotransposon insertions older than the estimated speciation time are polymorphic in almond, peach contains To gain insight into the evolution of almond and peach very few old polymorphic insertions, suggesting that they LTR retrotransposons, we grouped all almond and peach were lost in this species (Figure S10). elements into clusters showing sequence identity higher As all these analyses may be somewhat biased by a dif- than 80% along more than 80% of their length. Most of the ferent degree of assembly of peach and almond genomes, elements (66.2%) were grouped into clusters of at least we performed a detailed analysis of the LTR retrotrans- two elements. Two hundred and fifty-nine clusters (81%) poson insertions comparing orthologous loci in both spe- were mixed clusters and contained 93% of the almond and cies. We were able to unambiguously identify the peach clustered elements. An analysis of the insertion orthologous locus for 1155 full-length LTR retrotrans- times of these LTR retrotransposons shows that the num- poson peach insertions and 1134 almond insertions ber of recent (≤5 Mya) LTR retrotransposons is clearly (around 51% of the insertions in both species). These cor- higher in peach than in almond (Figure 3a). An analysis of respond to 142 insertions found intact in both species the insertion time distribution of individual clusters within (conserved insertions), 592 and 440 specific insertions in LTR retrotransposon families (belonging to Gypsy and peach and almond, respectively, 422 peach insertions par- Copia superfamilies, or that are unclassified), shows that tially deleted or rearranged in almond and 562 almond many of them contain insertions that are younger in peach insertions partially deleted or rearranged in peach. An that in almond (Figure S9), suggesting that peach has analysis of the ages of these LTR retrotransposon inser- experienced higher LTR retrotransposon activity after the tions belonging to the different categories showed that, evolutionary split of these two species. as expected, the majority of the specific insertions in both To further understand LTR retrotransposon dynamics we genomes were younger than the estimated speciation analyzed the prevalence in the species of the LTR retro- date. The almond genome contains a larger fraction of transposon insertions found in peach and almond refer- specific insertions that are older, which probably corre- ence genomes by analyzing re-sequencing data from 10 sponds to elements that were polymorphic in the ancestor peach and 10 almond cultivars. This analysis shows that and that were subsequently lost in peach. Also as the LTR retrotransposon insertions are frequently polymor- expected, the vast majority of the conserved insertions phic among almond cultivars whereas they often appear were older than the estimated speciation time of peach fixed in peach. An analysis of the insertion time distribu- and almond (Figure S11). In addition, analysis of the pres- tion of fixed and polymorphic LTR retrotransposons in ence of these insertions in peach and almond cultivars both species shows that whereas an important number of using Pindel (see Experimental Procedures) showed that © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 461 Figure 3. Dynamics of long terminal repeat (LTR) retrotransposons in peach and almond. (a) Insertion time of complete LTR retrotransposons in Prunus dulcis and Prunus persica. (b) Insertion time (MYA, million years ago) of polymorphic and fixed orthologous LTR retrotransposons in almond (left) and peach (right). [Colour figure can be viewed at wileyonlinelibrary.com] while all insertions are fixed in peach an important frac- peach and other Prunus species bitter and toxic. In order tion is polymorphic in almond (Figure 3b). to shed light on the genetic differences underlying these phenotypic differences we compared the genomic regions Transposon-induced variability in peach and almond traits containing the genes known to determine the expression The almond fruit resembles that of peach and other Prunus of these characters in both species. species, the major differences being that in almond the It has been recently shown that the sweet almond phe- mesocarp does not develop to produce the fleshy tissue notype is due to the reduced expression of the genes typical of other Prunus fruit crops, and that the almond encoding two cytochrome P450 enzymes catalyzing the seed does not accumulate the high levels of the cyano- first steps of amygdalin biosynthesis in sweet almond vari- genic diglucoside amygdalin that renders the seeds of eties compared with bitter almond varieties (Thodberg © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 462 Tyler Alioto et al. et al., 2011). It has also been shown that this reduced et al., 2011). We have analyzed the structure of these genes expression is not related to differences in the gene in peach and almond and found that six of them present sequence, which points to a difference in the regulation of differential TE insertions within the genes in their proximal the expression of those genes (Thodberg et al., 2011). We (less than 1000 nucleotides) upstream region that probably compared the sequence of one of these almond genes, contains their promoter. These species-specific TE inser- CYP71AN24, with its homologs in peach, sweet cherry and tions are all highly methylated and may have altered their P. mume and found that it is highly conserved. However, expression (Table S18). In addition to the potential muta- CYP71AN24 is flanked by several almond-specific TE inser- tion of transcriptional regulatory elements, MITE insertions tions, and in particular two MITE insertions in its proximal could have provided other transcription factor-binding upstream region (Figure 4a). A preliminary analysis of the sites (TFBS), as MITEs have been shown to frequently methylation of this region shows that the almond-specific amplify and mobilize TFBS in plants (Morata et al., 2018). TE insertions flanking the CYP71AN24 gene are highly DISCUSSION methylated. The insertion of TEs in the proximal upstream region of the CYP71AN24 gene may have affected, directly Using a hybrid strategy based on short and long-read DNA or indirectly, its expression due to its high methylation. In sequences and the information provided by existing link- addition, the presence of the TEs also correlates with a age maps we have assembled the highly heterozygous much higher methylation of this gene in almond compared genome of almond cv. Texas into a rather complete, con- with peach (11.6% methylation at CG and 0.1% methylation tiguous and low-redundant assembly with eight pseudo- at CHG in almond versus 0.1% methylation at CG and molecules corresponding to the eight chromosomes and 0.02% at CHG in peach), which may be the result of methy- comprising 227.6 Mbp of sequence. Annotation of this lation spreading from the TEs. An analysis of the structure genome has resulted in the identification of 27 969 pro- of this locus in Prunus webbii, a wild species closely tein-coding genes and 6747 non-coding transcripts. The related to almond which produces bitter seeds, and in two assembled sequence is highly syntenic with the genome sweet almond cultivars (Cristomorto and Marcona) and sequence of peach (Verde et al., 2017), as was expected two almond cultivars producing bitter seeds (D05-187 and considering previous information on the close genetic sim- S3067) performed by mapping the re-sequencing data of ilarity between these two species (Dirlewanger et al., 2004). these genomes to the almond reference genome, shows Based on molecular data and the use of fossil records of that the presence of the TE insertions (and in particular a diverse sample of 17 plant species we estimated the that of the MITE named TIR2), correlates with the sweet divergence times of P. dulcis with respect to other versus bitter seed phenotype (Figure 4b). Moreover, an sequenced Prunus. Our estimate of 5.88 Mya for the diver- article published after the submission of our research gence of peach and almond from a common ancestor is reports a strong association between the sweet almond similar to that of recent molecular evidence that places this phenotype and a point mutation in a gene (bHLH2) encod- figure between 5.0 and 8.0 Mya (Yu et al., 2001; Velasco ing for a bHLH transcription factor that renders it unable to et al., 2014; Delplancke et al., 2016). This is compatible with bind to the promoter of one of these genes, CYP71AN24 the separation of the ancestral species by the uplift of the (Sanchez-Perez et al., 2019). An analysis of the sequence of Central Asian Massif in two subpopulations that faced the genome produced here, the sweet almond Texas, as completely different environments: one (almond and its well as the mentioned sweet and bitter accessions and close relatives) in the arid steppes of central and western P. webbii shows that this mutation is absent in Texas and, Asia and the other (peach) in the subtropical climate of on the contrary, it is present in P. webbii, being as southwestern China, close to where the first fossil endo- expected in the rest of the cases. Additionally, Sanchez- carps were found dated at 2.6 Mya (Su et al., 2015). Perez et al. (2019) reported that one of the sweet almond In agreement with results from earlier studies of other cultivars they analyzed (Atocha), did not have the same sequenced diploid Prunus genomes (Shirasawa et al., mutation in bHLH2 but had another one in the neighbor- 2017; Verde et al., 2017; Zhang et al., 2017; Baek et al., hood. 2018) we have not found evidence for any recent whole- The analysis of model plant species, such as Arabidopsis genome duplication of almond. Analysis of duplicated and tomato, has shown that the development of the fruit is gene sequences indicates a parallel gene expansion for all the result of the combined action of genes involved in sequenced Prunus species genomes for genes involved in meristem organization, floral development and fruit cell lignin biosynthesis, such as COMT and HCT. One of the proliferation and expansion. We selected 97 genes distinctive aspects of Prunus is that its fruit is a drupe, (Table S17) that belong to gene families involved in these characterized by the formation of a strongly lignified meso- processes, including WUSCHEL (WUS) and CLAVATA carp (the ‘stone’), unlike most of its closest taxa that have (CLV) genes whose mutants lead to larger fruits in tomato follicetum and nuculanium as fruit types (Xiang et al., and genes known to determine fruit shape (Rodriguez 1990). These duplications may be at the origin of the © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 463 Figure 4. Analysis of the locus of the CYP71AN24 gene in almond varieties and related Prunus species. (a) Nucleotide conservation of the CYP71AN24 region between Prunus avium, Prunus mume, Prunus persica and Prunus dulcis based on a Mauve multiple alignment (physical distance scale is in bp). White boxes represent inserted regions in P. dulcis. The Integrative Genomics Viewer (IGV) tracks of the gene and transposable element (TE) annotations of P. dulcis and P. persica along with their DNA methylation levels in the three different contexts (CG, CHG and CHH) are shown below. (b) The IGV spnapshot of the region containing the CYP71AN24 gene and the polymorphic TE insertions displaying the coverages of mapped DNA-seq reads from re-se- quencing data of sweet- and bitter-kernel P. dulcis varieties, as well as from that of the closely related Prunus webbii. [Colour figure can be viewed at wileyonlinelibrary.com] formation of the stone, determine its characteristics and be differences, particularly those of larger sizes, were related crucial to understanding its evolution and possible modifi- to TE sequences. The LTR retrotransposons and MITEs cation, with important consequences for plant breeding, constitute the two most prevalent superfamilies of TEs in including the production of stoneless cultivars in Prunus plants (Casacuberta and Santiago, 2003). While the propor- fruit (Callahan et al., 2015). tion of the almond genome consisting of TEs (38%) was As we have shown from the comparison of the peach similar to that of other sequenced genomes of similar size, and almond reference genomes, as well as the analysis of that is from 30% in peach (Verde et al., 2013a,b) and 43– the structural variants between both genomes, indel events 47% in other Prunus species (Shirasawa et al., 2017; Zhang seem to explain a substantial part of the divergence of et al., 2017; Baek et al., 2018) and very similar chromoso- peach and almond genomes from their common ancestor. mal distributions of TEs were observed for almond and In this study we show that most of such structural peach, detailed analysis of the dynamics of LTR © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 464 Tyler Alioto et al. retrotransposon evolution has revealed key aspects of the P. webbii) with unexpected genotypes. This shows that, divergence of almond and peach genomes after their spe- although there is a good correlation of this mutation in the ciation. In short, almond has maintained more ancestral bHLH transcription factor gene and the sweet almond phe- LTR retrotransposons, which are still in some cases poly- notype (Sanchez-Perez et al., 2019), this association may morphic within the species, whereas peach has lost most not be perfect, suggesting that other mechanisms to avoid polymorphic ancestral insertions but seems to have wit- the activation of CYP71AN24 could contribute to the sweet nessed a higher level of recent retrotransposition activity. kernel phenotype. We propose here that the presence of The contrasting data on the polymorphism of ancestral highly methylated TE insertions in the close proximity of TEs between peach and almond may be explained by (i) CYP71AN24 could help to ensure the low expression of this the mating types (selfing in peach and outcrossing in gene in the seed tegument of sweet almonds. Additionally, almond) and (ii) differences in their recent history, with a TEs may be at the origin of other important traits in almond. strong reduction of population size in peach prior to its We identified six genes related to fruit flesh formation that recent expansion as a cultivated species 2000 years ago also contain highly methylated TEs in their upstream region (Velasco et al., 2014), while almond population sizes in either almond or peach. Although we have not been able remained higher (Yu et al., 2001). Further analyses should to associate any of these genes with the positions of major help to clarify this in the future. In any case, the results pre- genes or quantitative trait loci described so far in these spe- sented here suggest that LTR retrotransposons, and in gen- cies, they are clear candidates for further studies. eral TEs, may explain an important fraction of the In summary, peach and almond diverged less than 6 interspecific variability between peach and almond, as well Mya, which for a perennial species with a long intergenera- as the intraspecific variability of both species. tional period of around 10 years, is a short evolutionary Long terminal repeat retrotransposons are at the origin of time for sequence divergence. Indeed, our results show somatic mutations in plant species, some of which have that the genomes of peach and almond are highly similar high agricultural value (Foster and Aranzana, 2018). Only in and only show a mean of 20 nucleotide substitutions per peach, white versus yellow fruit color (Falchi et al., 2013), kb. Our results here show that this relatively low sequence hairy versus glabrous fruit (Vendramin et al., 2013) and divergence is accompanied by an important number of stonyhard versus melting flesh texture (Tatsuki et al., 2018) indels frequently resulting from specific TE insertions, are caused by the action of transposon movement. The which in some cases could be at the origin of important MITE insertions have also been linked to crop traits, such as almond characteristic traits. The activity of TEs is not con- sex determination in melon (Martin et al., 2009) or a stant through evolution, alternating between quiescent drought tolerance phenotype in maize (Mao et al., 2015). periods and transposition bursts, which, as our results sug- The analysis of almond 9 peach interspecific progeny iden- gest here for almond, could allow for a rapid phenotypic tified 11 Mendelian genes explaining the inheritance of diversification between closely related species. some key agronomic characters: one of them, responsible The genome sequence of the almond will accelerate for the formation of the thick mesocarp that constitutes the genetic research and facilitate breeding of this species by peach flesh, and another, conferring juiciness to the fleshy providing useful information on genes and markers with mesocarp, represent major contributions to the difference an unprecedented level of detail, as is the case in other between almond and peach fruits (Donoso et al., 2016). Our Prunus species (Aranzana et al., 2019). It can also help fur- results suggest that TEs could be responsible for some of ther our understanding of the evolution and domestication the genomic changes at the origin of the agronomic traits of closely related crop species that share a slower rate of that distinguish peach from almond, such as mesocarp evolution due to their long intergeneration period, and development and bitterness of the kernel. For one of them, may enable the detection of modes and aspects of evolu- sweet versus bitter kernel, which is essential for the domes- tion that could be different or otherwise difficult to identify tication of the almond, we show here that the sweet almond in herbaceous crops (Gaut, 2015). phenotype correlates with the presence of TE insertions sur- rounding the gene CYP71AN24. This gene is involved in the EXPERIMENTAL PROCEDURES synthesis of one of the key enzymes of the amygdalin path- Extraction of DNA, sequencing and k-mer analysis way (cytochrome P450), and it has been proposed that its reduced expression, together with the lack of expression of Fresh young leaves of the Texas almond were ground in liquid CYP79D16, results in the sweet kernel trait (Thodberg et al., nitrogen to a fine powder, and 100 mg of ground leaves was used 2011). It has been recently shown that the sweet kernel phe- for DNA extraction using the DNeasyâ Plant Mini Kit (Qiagen, https://www.qiagen.com/) according to the manufacturer’s instruc- notype is closely associated with a point mutation in a bHLH tions. The same method was employed for the extraction of 10 transcription factor (Sanchez-Perez et al., 2019). These additional almond cultivars (A€ı, Belle d’Aurons, Cristomorto, Des- authors identified a sweet almond cultivar that did not have mayo largueta, Falsa Barese, Genco, Marcona, Nonpareil, Ripon the mutation, and we found two more examples (Texas and and Vivot), eleven peach cultivars (Armking, Belbinette, BigTop, © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 465 Blanvio, Catherine, Earlygold, Flatmoon, Nectalady, Platurno, genome size (512 Mbp). Heterozygosity and repeats were clearly Sweetdream and Tiffany) and one accession of the wild almond going to be a problem. Fosmid-pool sequencing (150/250 nucleo- P. webbii (R755). For sequencing with an Oxford Nanopore Tech- tide PE reads) and assembly combining whole-genome sequenc- nologies (ONT; https://nanoporetech.com/) MinION sequencer, ing PE and MP reads as previously described (Abascal et al., 2016; high-molecular-weight DNA from Texas almond was extracted Cruz et al., 2016) increased contiguity (N50 = 142 kb) and reduced with the method described by Mayjonade et al. (2016). the assembly size (238 Mbp); however, the resulting assembly Whole-genome shotgun sequencing was performed for Texas exhibited high levels of discordance with the available genetic DNA using the Illumina HiSeq2000 sequencing instrument. The maps, as well as with the peach assembly (v.2.0.1). The final strat- standard Illumina protocol with minor modifications was followed egy, combining long-read ONT sequencing and whole-genome for the creation of short-insert paired-end (PE) libraries (Illumina sequencing PE reads, resulted in the best balance of contiguity Inc., cat. no. PE-930-1001, https://www.illumina.com/). In brief, three and concordance with the genetic map. libraries were generated from > 2.0 lg of genomic DNA each. For The hybrid assembler MaSuRCA v.3.2.3 (Zimin et al., 2013) was one library the DNA was amplified by PCR while for the other two run with default parameters (no linking mates; Celera assembly of libraries (and the majority of the sequence) the DNA was not ampli- super-reads). The input was the two PE Illumina libraries of insert fied in order to reduce GC bias. Then the DNA was sheared on a Cov- sizes 317 and 354 bp for a total of 2859 coverage (Table S1) and arisTM E220, the fragmented DNA was size selected on an agarose the self-corrected ONT reads for a total of about 209 coverage. gel to obtain three PE libraries with incremental insert sizes of 263, Redundans v.0.13 (Pryszcz and Gabaldo n, 2016) was run on the 317 and 354 bp. The fragments were end-repaired, adenylated and non-deduplicated output of MaSuRCA (the step 9-terminator ligated to Illumina indexed PE adaptors. The PE libraries were run genome.ctg.fasta file totaling 470 Mb with N50 = 53 kb) using the on the of Illumina HiSeq2000 platform in 2 9 101 PE mode accord- –longreads option for scaffolding, reducing the size of the assem- ing to standard Illumina operation procedures. Primary data analy- bly to 248 211 Mb of which was scaffolded further with the ONT sis was carried out with the standard Illumina pipeline (HCS reads. 2.0.12.0, RTA 1.17.21.3). A total of 97 Gb of raw sequence A first round of corrections to the assembly was carried out (>3509 coverage) was produced. Post-processing of sequence using consistency with nanopore data as the main criteria. ONT reads involved detection and trimming of Illumina adapter reads were mapped to the assembly with NGM-LR v.0.2.6 (Sed- sequences with cutadapt, quality trimming with trim_galore and PE lazeck et al., 2018), the assembly was broken at regions of zero overlap detection and merging with FLASH. Mate pair (MP) libraries coverage and then re-scaffolded with SSPACE-LongRead v.1.1 of Texas DNA (fragment sizes 3.1 and 5.2 kb) were constructed (Boetzer and Pirovano, 2014). A second round of corrections was according to the Nextera Mate Pair Preparation protocol, which made utilizing collinearity with the T9E genetic F2 and BC1 (to leaves a linker of known sequence at the junction. The resulting Earlygold) linkage maps (Donoso et al., 2015) as the main crite- libraries were run on the HiSeq2000 platform in 2 9 101 bp read rion, with break points guided by synteny with the peach genome length runs. Post-processing of sequence reads involved detection and coverage of ONT reads. Peach transcripts (annotation and trimming of the Nextera linker sequence with cutadapt, and Pp2.01a) were mapped to the almond assembly with GMAP quality trimming with trim_galore. v.2014-12-23 (Wu and Watanabe, 2017). Marker sequences were Fosmid pools from Texas almond were used to prepare 2D mapped with BWA mem, keeping those mappings with mapping genomic libraries using the Ligation Sequencing Kit SQK-MAP005/ quality ≥ 20 and identity ≥ 90%. The broken assembly was again 006 or SQK-NSK007. The sequencing run was performed on an scaffolded with SSPACE-LR. The assembly at this stage had a con- R7.3 chemistry FLO-MAP103 or R.9 FLO-MIN104 flow cell (ONT) tig N50 of 99 kb and scaffold N50 of 151 kb. according to the manufacturer’s recommendations. Genomic DNA A third round of corrections was performed with improvements from Texas almond was also used to prepare 1D and 1D2 genomic in mapping and break detection. First, peach transcripts were libraries using the Ligation Sequencing Kits SQK-LSK108 and mapped only in the sense direction, and second, discrepant mar- SQK-LSK308, respectively. The sequencing run time on FLO- ker mappings were screened for mapping artifacts. Moreover, MIN106 or FLO-MIN107 flow cells using the MinION/MKI Pk.1 Sniffles v.1.0.11 (Sedlazeck et al., 2018) was used for detection of instruments (ONT) was 48 h. The quality parameters of the structural variants. Additional breaks were made and duplicate sequencing runs were further monitored by the MinKNOW plat- sequences were also detected. In the end, we were able to merge form, while the run was base-called using the Metrichor agent 30 Mb with minimus2 from AMOS v.3.1.0 (Sommer et al., 2007; (https://metrichor.com) in real time. Treangen et al., 2010) and any remaining duplicate sequence For the almond cultivars other than Texas and the peach and (>99% identical, >5 kb) was manually reviewed. Overlapping P. webbii accessions, we developed PE libraries of fragment size regions were joined into new longer contigs using nanopore read 300 bp and sequenced them with Illumina as described in the pre- mappings to confirm new joins. Also at this stage, the putative vious paragraph. chloroplast sequence was identified by coverage and homology Jellyfish v.2.2.6 (Marcßais and Kingsford, 2011) was run on the and set aside. PE300 library (insert size 317 bp, 2 9 100-nucleotide reads) with Finally, the assembly was anchored to pseudomolecules using the canonical k-mer (-C) option and a mer-size of 21. Genome both the T9E genetic map and synteny with the peach genome Scope (Vurture et al., 2006) was then used to analyze the resulting using ALLMAPS (jcvi-0.7.3) (Tang et al., 2018), with more weight k-mer distribution. The ONT reads were error-corrected using given to the map marker order. Remaining conflicts were resolved Canu (v1.5) (Koren et al., 2017). Corrected trimmed reads were manually. ALLMAPS uses a genetic algorithm for placing and ori- used for hybrid assembly, scaffolding and assembly correction. enting scaffolds, and sometimes it does not converge completely on the optimal solution, even with a large number of generations. Genome assembly Thus, we had to manually review and fix the order and orientation of some scaffolds which still exhibited discordance with either the Short-read whole-genome sequence (WGS) assembly of 100 bp genetic map or synteny with peach. Further improvement to the PE read Illumina libraries with AbySS v.1.3.6 (Simpson et al. 2009) assembly was made by joining adjacent scaffolds if they could be resulted in a fragmented assembly (N50 = 4867 bp) with inflated linked together with split nanopore read mappings. A few © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 466 Tyler Alioto et al. additional overlaps were also detected in this fashion and longer present in the genome assembly. To annotate long non-coding contigs were constructed. RNAs (lncRNAs) we first selected PASA assemblies that had not Assembly completeness was estimated in two ways. First, gene been included in the annotation of protein-coding genes. Those completeness was determined by running BUSCO v.3.0.2 (Sima~o longer than 200 bp and whose length was not covered to at least et al., 2015) using the embryophyta_odb9 database comprising 80% by a small ncRNA were incorporated into the ncRNA annota- 1440 single-copy plant orthologous groups (BUSCOs). Second, a tion as lncRNAs. The resulting transcripts were clustered into pairwise comparison of k-mers present in both input reads and genes using shared splice sites or significant sequence overlap as the assembly was performed using KAT (Mapleson et al., 2017) criteria for designation as the same gene. using all whole-genome sequencing PE Illumina reads and a k- mer length of 27 (Figure S12). Functional annotation. Functional annotation was performed by integrating several data sources to infer protein function based Comparison of the P. dulcis anchored assembly to the on sequence similarity to annotated sequences or/and the pres- linkage map and peach v.2.0 a1 genome sequence ence of particular domains and sequence motifs. We used the InterPro (Hunter et al., 2012), KEGG (Kanehisa et al., 2012), signalP The almond assembly was compared with the T9E linkage map (Petersen et al., 2011) and NCBI CDsearch (Marchler-Bauer et al., that contains 1833 SNP markers (Donoso et al., 2015). Markers 2011) databases. InterProScan v.5.19-58 (Zdobnov and Apweiler, were mapped onto the almond pseudomolecule-based assembly 2012) was used to scan though all available InterPro databases, using BLAST and coordinate data of both almond and peach were including PANTHER, Pfam, TIGRFAM, HAMAP and SUPERFAMILY. used as input in MapChart software (Voorrips, 2017) for represent- Initial sequence similarity search was determined using BLASTP ing graphically the comparison between the two species. The v.2.6.0 + against the NCBI non-redundant (NR) collection of pro- genetic and physical distances of SNP markers from the T9E pop- tein sequences (release 2018-08). KEGG orthology (KO) groups ulation were used for calculating the recombination rate across were assigned by the KEGG Automatic Annotation Server (KAAS) the pseudomolecules of the almond assembly. (Moriya et al., 2007) using the bi-directional best hit (BBH) method The peach genome sequence and annotation data were down- against a representative gene set from 27 different species, which loaded from the GDR (ftp://ftp.bioinfo.wsu.edu/species/Prunus_per includes a core set of species for gene annotation and additional sica/Prunus_persica-genome.v2.0.a1/). Synteny between almond plant species from the Rosaceae family. The KO identifiers were with peach genomes was assessed using SyMap software v.4.2 then used to retrieve the relevant KEGG functional annotation (Soderlund et al., 2007) with default parameters, except that the using the KEGG REST-based API service, KEGG release v.87.1. ‘min dot’ parameter was set to 25. To predict plant disease resistance genes, each protein was searched against a manually curated list of ‘reference’ R-genes Annotation with the DRAGO pipeline (Sanseverino et al., 2013). For each hit, The P. dulcis genome assembly was annotated by combining classes were assigned based on combination of specific domains, such as TIR, nucleotide-binding site (NBS), leucine-rich region transcript alignments, protein alignments and ab initio gene pre- dictions. A flowchart of the annotation process is shown in Fig- (LRR) or coiled-coil domain. Putative transcription factor genes ure S13. Scripts are available at https://github.com/jesgomez/ were predicted using the Plant Transcriptional Factor database annotation_pipeline. (Jin et al., 2017) v.4.0 First, almond RNA-seq reads were downloaded from NCBI with the accession number SRR1251980 and aligned to the genome Annotation and analysis of transposable elements in gen- with STAR (v.2.5.3a) (Dobin et al., 2013). Transcript models were ome assemblies. The Ilumina PE reads corresponding to the subsequently generated using Stringtie (v.1.0.4) (Pertea et al., re-sequencing of the almond and peach cultivars described in sec- 2015) and, along with the P. persica transcriptome (annotation tion "DNA extraction, sequencing and K-mer analysis" were Pp2.0a) and 4509 almond expressed sequence tags downloaded trimmed with SKEWER (v.0.2.2, –mean-quality 25, –min 35) and from NCBI on July 2015, were assembled into a non-redundant aligned to their respective reference genome with BWA-aln/ set by PASA (v.2.3.3) (Haas et al., 2008). The TransDecoder pro- sampe1 (v.0.7.5, parameters: -t 6, -n 5, -o 1, -e 3) (Li and Durbin, gram, which is part of the PASA package, was run on the PASA 2009) and SAMTOOLS (v.0.1.18) (Li, 2011). Bam files were later assemblies to detect coding regions in the transcripts. Second, the submitted to the package PINDEL (v.0.2.5, parameters: -T 4, -x 5, -r complete Rosaceae proteome was downloaded from Uniprot on false, -t false, -A 35) (Ye et al., 2018) to identify deletions in sam- July 2015 and aligned to the genome using Exonerate (v.2.4.7) ples. The TEs were annotated in P. dulcis and P. persica assem- (Slater and Birney, 2005). Third, ab initio gene predictions were blies using TEdenovo and TEannot pipelines of the REPET performed on the repeat masked pdulcis26 assembly with three package (Quesneville et al., 2005; Flutre et al., 2011) installed in different programs: GeneID v.1.4 (Alioto et al., 2018), Augustus the PiRATE virtual machine (Berthelier et al., 2018). Classification v.3.2.3 (Stanke et al., 2015) and GeneMark-ES v.2.3e (Lomsadze of TEdenovo consensus sequences at the order level was done et al., 2014) with and without incorporating evidence from the with PASTEC (Hoede et al., 2014). RNA-seq data. Finally, all the data were combined into consensus Annotation of TEs for masking purposes was done using coding sequence models using EvidenceModeler-1.1.1 (EVM) RepeatMasker (http://www.repeatmasker.org/) with a reduced TE (Haas et al., 2008). Additionally, untranslated regions and alterna- representative library. The TE representatives obtained using the tive splicing forms were annotated through two rounds of PASA TEdenovo pipeline library were screened for coding domains with annotation updates. hmmscan (HMMER 3.1b1, http://hmmer.org/) against the PFAM Non-coding RNAs were annotated as follows: first, the program database (Finn et al., 2016). The TE representatives containing cmsearch v.1.1 (Cui et al., 2016) from the INFERNAL package regions potentially coding for known domains of non-TE proteins (Nawrocki and Eddy, 2013) was run against the RFAM (Nawrocki usually found in multigene families (kinases, NB-ARC, LRR, TIR) or et al., 2015) database of RNA families (v.12.0). Also, tRNAscan-SE with an N content higher than 30% or with length shorter than 200 v.1.23 (Lowe, 1997) was run to detect the transfer RNA genes nucleotides were discarded. Moreover, all TE representatives not © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 467 categorized in one of the classical TE superfamilies (defined as AMPure XP beads (Beckmann Coulter, Brea CA, USA). The ‘noCat’ by TEdenovo) were also removed. A total of 661 represen- libraries were constructed using the Kapa Library Preparation kit tatives were removed from the library. The final library contained (Roche Kapa Biociences, Pleasanton CA, USA) for short-insert 6898 TE representatives. paired-end libraries for Illumina with some minor modifications. MITE-hunter (Han and Wessler, 2010) was run to detect poten- After ligation of the NEXTFLEXâ Bisulfite-Seq Barcodes (Perkin tial MITE families. In order to complete the annotation, the poten- Elmer, https://www.perkinelmer.com/) the library was treated with tial almond MITE families were combined with the P. persica sodium bisulfite using the EpiTect Bisulfite Kit (Qiagen), following family annotation available in the PMITE database (families carry- the manufacturer’s instructions for formalin-fixed, paraffin-embed- ing target site duplication) (Chen et al., 2014). These sequences ded tissue samples. Two rounds of bisulfite conversion were per- were grouped in clusters of 90% identity with cd-hit (Fu et al., formed to ensure a conversion rate of over 99%. Enrichment for 2012) to remove redundancy and produce a final library of family adaptor-ligated DNA was carried out through seven PCR cycles representatives. RepeatMasker (http://www.repeatmasker.org/) using KAPA HiFi Uracil + DNA Polymerase (Kapa Biosystems, was run to annotate all regions having significant similarity to https://www.kapabiosystems.com/). Library quality was monitored MITE families, and the results were filtered to retain only full- using the Agilent 2100 Bioanalyzer, and the library concentration length elements (consensus length  20%). The same pipeline was estimated using quantitative PCR with the library quantifica- was used to identify MITEs in the P. persica assembly. tion kit from Roche Kapa Biosystems. Paired-end DNA sequencing Candidate LTR retrotransposons were predicted by running (2 9 101 + 8 bp) was then performed using the HiSeq2500 (Illu- LTRharvest (Ellinghaus et al., 2008) with default parameters. The mina) following the manufacturer’s protocol. Image analysis, base internal conserved domains of these elements were identified calling and quality scoring of the run were processed using the using HMMER hmmscan (Johnson et al., 2010) and only coding manufacturer’s software Real Time Analysis (RTA 1.18.66.3) and elements were retained for further analyses. Elements displaying followed by generation of FASTQ sequence files. Raw reads were either a single hit on the genome, more than 10% of gaps or more trimmed with TrimGalore! v.0.4.5 (http://www.bioinformatics.bab than 50% of tandem repeats were filtered out. Classification of the raham.ac.uk/projects/trim_galore/). Low-quality bases (Phred remaining elements (hereafter referred to as ‘coding LTR retro- score < 20) were trimmed before adapter removal and reads with transposons’) into Copia and Gypsy superfamilies was performed a length less than 20 were discarded. The total of trimmed reads based on the order of the internal coding domains, as defined by was 82 073 678 and 94 040 426 in almond and peach, respec- Xiong and Eickbush (2006). Elements lacking one or more tively. Trimmed reads of each species were mapped to their domains were tagged as ‘incomplete’. respective reference genome and methylation was analyzed using Bismark v.0.19.1 (Krueger and Andrews, 2011). The gene/TE The LTR regions of every coding element were extracted and methylation was analyzed with SeqMonk v.1.41 (http://www.bioinf aligned with MUSCLE (Edgar, 2004). The Kimura two-parameter ormatics.babraham.ac.uk/projects/seqmonk/). Only cytosine posi- distance of every aligned LTR pair was calculated and used to esti- tions that had been sequenced at least three times were included. mate insertion ages following the approach described in (SanMi- guel et al., 1998), using a substitution rate of 108 nucleotides per Analysis of the CYP71AN24 locus in almond cultivars and site per year and a generation time of 10 years (Velasco et al., Prunus-related species 2014). The flanking sequences (500 bp) of every coding LTR retrotrans- A 2-Mb region of the P. dulcis genome containing the CYP71AN24 poson were extracted from P. dulcis and used as query for a gene was compared with the corresponding genomic regions of BLASTn (Altschul et al., 1990) search (cutoff E-value < 10–10) P. avium, P. mume, P. persica using Mauve (Darling et al., 2004). against the P. persica assembly, and vice versa. Concordant Re-sequencing data from the P. dulcis cultivars Texas, Marcona, mapped flanks were defined when both flanks of an element Cristomorto, D05-187 (SRX245830) and S3067 (SRX245832), as mapped in the same scaffold at a distance smaller than 25 kb. well as from P. webbii (R755), were mapped to the almond refer- Every internal region between two concordant mapped flanks was ence genome using BWA aln/sampe (Li and Durbin, 2009). aligned to the putative orthologous element using EMBOSS Nee- dle (Rice et al., 2000). The two elements were considered ortholo- Prunus dulcis phylome reconstruction gous if they could be aligned over 80% of their length with at least 80% identity. To assess the orthology of MITE insertions the The P. dulcis and P. persica phylomes, that is the complete collec- approach followed was as the one described for LTR retrotrans- tion of evolutionary histories of all encoded genes, were recon- posons except that the sequences flanking the insertions were structed using the PhylomeDB pipeline (Huerta-Cepas et al., 2011). mapped to the corresponding genome using BBmap (https://sour In brief, for each protein-coding gene in the almond and peach ceforge.net/projects/bbmap/) instead of Blast. genome we searched for homologs (Smith–Waterman Blast search, E-value cutoff < 1 9 10–5, minimum contiguous overlap In order to search for polymorphic LTR retrotransposon and over the query sequence cut-off ≥ 50%) in a database containing MITE insertions within or close to genes we used BEDTools the proteomes of 17 species with sequenced genomes represent- (v.2.27.0) (Quinlan and Hall, 2010). Only those TEs located within ing most of the important plant families (Table S3). The most sim- genes or at less than 1000 nucleotides upstream of a gene were ilar 150 homologs were aligned using three different programs, kept. MUSCLE (Edgar, 2004), MAFFT (Katoh et al., 2009) and KALIGN Analysis of DNA methylation in almond and peach (Lassmann and Sonnhammer, 2005), in forward and reverse orien- tation. These six alignments were combined using M-COFFEE Genomic DNA (1.5–2 lg) from young leaves of P. dulcis (cv. Tex- (Wallace et al., 2008) and trimmed with trimAl v.1.3 (Capella- as) and P. persica (cv. Earlygold) was spiked with unmethylated Gutierrez et al., 2009), using a consistency cut-off of 0.16667 and a bacteriophage k DNA (5 ng of k DNA/lg of gDNA; Promega, gap threshold of 0.1. Phylogenetic trees were built using a maxi- https://www.promega.com/) and with methylated T7 phage DNA mum likelihood approach as implemented in PhyML v.3.0 (Guin- (5 ng of T7 DNA lg–1 of gDNA). The gDNA was sheared on a Cov- don and Gascuel, 2003) using the best fitting model among seven arisTM E220 and fragments of 150–300 bp were size-selected using different ones (JTT, LG, WAG, Blosum62, MtREV, VT and Dayhoff). © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 468 Tyler Alioto et al. The model best fitting the data was determined by comparing the death prior, and a prior on the root of the tree (183 Mya; Bell likelihoods estimated on an initial neighbor joining tree topology et al., 2010). Two independent Markov chain Monte Carlo chains and using the Akaike information criterion. In all cases we used were run for 20 000 cycles, sampling posterior rates and dates four rate categories and inferred the fraction of invariant positions every 10 cycles. The initial 25% were discarded as burn-in. Poste- and rate parameters from the data. Then, these phylomes were fil- rior estimates of divergence dates and associated 95% credibility tered to remove the gene trees that contain proteins associated intervals were then computed from the remaining samples of each with transposon-related functional terms. All alignments and trees chain. are available for browsing or download at PhylomeDB with the PhylomeID 406 (almond phylome) and 407 (peach phylome) Re-sequencing of almond cultivars and comparison of (Huerta-Cepas et al., 2014) (http://www.phylomedb.org/). almond–peach structural genome variability Genetic variability analysis was performed on 10 traditional Prediction of orthology and paralogy and detection of gene almond cultivars and one peach cultivar, Earlygold, that was used duplications. Orthology and paralogy relationships were pre- as an outgroup. These accessions were re-sequenced using PE dicted based on phylogenetic evidence from the almond and Illumina sequencing as described before. Selection of the almond peach phylomes. We used ETE v.3 (Huerta-Cepas et al., 2010) to lines was based on their origin as representing a range of the infer duplication and speciation relationships using a species major areas of production in Europe (France, Spain, Italy) and the overlap approach and a species overlap score of 0. In brief, the USA, and morphological characteristics (shell hardness, bloom algorithm traverses the tree from the tip to the root and, for time and self-incompatibility) (Table S8). each node, evaluates whether the two daughter branches con- The PE Illumina sequencing data from the almond cultivars were tain genes from the same species, in which case a duplication is trimmed (length ≥ 35 bp, mean sliding window of 4 bp, phred qual- inferred, and the genes in each of the two splitting branches are ity score ≥ 20) using Trimmomatic (Bolger et al., 2014). and the out- considered paralogous to each other (Gabaldo  n, 2008). The rela- put was quality checked using FastQC (https://www.bioinformatic tive age of detected duplications was estimated using a phy- s.babraham.ac.uk/projects/fastqc/). Trimmed data were aligned lostratigraphic approach that uses the information on which against the almond assembly using the BWA-MEM algorithm species diverged before and after the duplication node (Huerta- v.0.7.16a-r1181 (http://bio-bwa.sourceforge.net/bwa.shtml) with Cepas and Gabaldo  n, 2011). Duplication frequencies at each default parameters. After removal of unmapped reads, PCR dupli- node in the species tree were calculated by dividing the number cates and reads with mapping quality < 10 we obtained the subset of duplications mapped to a given node in the species tree by of ‘clean’ reads used for variant calling, which was performed with all the gene trees that contain that node. To calculate duplication Samtools v.1.5 (Li, 2011) with default parameters, except from the frequencies we excluded gene trees that contained large (more following: -q 10 -Q 20. Commands for trimming, alignment, bam fil- than five paralogs) species-specific expansions (expansions that tering and variant calling can be found in the github repository contained more than five members). This filter is applied to (https://github.com/kostasgalexiou/sample-processing.git). Variant avoid the contribution of transposon-related gene families or calling format (VCF) files were filtered by applying the following cri- pseudogenes present in the other analyzed genomes. For the teria: global quality ≥ 30, genotype quality ≥ 30, 8 ≤ depth ≤300, rest of the analyses all duplications were considered. All orthol- biallelic sites, minor allele frequency (MAF) ≥ 0.1. Graphical repre- ogy and paralogy relationships are available through PhylomeDB sentation of variant distributions was done with Circos (Krzywinski (Huerta-Cepas et al., 2014). et al., 2009) in non-overlapping windows. Gene Ontology term enrichment analysis was performed using Large deletions between almond and peach re-sequencing data FatiGO (Al-Shahrour et al., 2007). We compared three lists of pro- were identified using Pindel (Ye et al., 2018) using default parame- teins against all the other proteins encoded in the genome. The ters and an insert size of 300 bp for all samples. Peach- and three lists were composed of the proteins involved in a duplication almond-specific deletions were obtained by selecting positions at the ancestral node of all Prunus species, the proteins specifically with at least 20 reads/cultivar supporting the event. We also lost in almond and the proteins specifically lost in peach. removed deletions that overlapped with N-regions (1000 bp) in The trimmed alignments of 262 genes that had single-copy the almond genome. For detecting variants that overlap with TEs, orthologs in the 17 species considered were selected and concate- we considered that a position overlaps with a TE if at least one of nated. The final alignment containing 141 911 amino acid posi- the two elements, the deletion or the TE, had at least 80% of its tions was used to reconstruct the maximum likelihood species sequence overlapping with the other element. tree with PhyML v.3.1 (Guindon et al., 2010) using the LG amino Prunus dulcis contigs were aligned to the P. persica reference acid substitution model and 100 bootstrap replicates. Additionally, genome using Nucmer from the Mummer3 package (Delcher a super-tree was reconstructed using all trees in the phylome and et al., 2003). Assemblytics (Nattestad and Schatz, 2016) was used a gene tree parsimony approach as implemented in duptree to filter the alignment and detect genome-wide variants with the (Wehe et al., 1889). following cut-offs: unique sequence length required for consider- Divergence dates were estimated on the topology derived from ing an alignment = 10 000 bp, minimum variant size = 20 bp, the maximum likelihood approach by using the Bayesian relaxed maximum variant size = 25 000 bp. Structural variants were inter- molecular clock approach as implemented in PhyloBayes v.4.1c sected with TE annotations of P. dulcis (insertions and repeat (Lartillot et al., 2013). An uncorrelated relaxed clock model was expansions) and P. persica (deletions and repeat contractions). A applied, and four fossil constraints specified to the most recent variant was considered to be TE-associated when at least 50% of common ancestor: Prunus (47.8 Mya; Li et al., 2011), Rosaceae its sequence was spanned by a TE. (98.25 Mya; Crepet and Nixon, 1996; Zhang et al., 2013), the split between Fagales and Cucurbitales (84 Mya; Herendeen et al., DATA AVAILABILITY 2002; Sims et al., 2002; Wikstro€ m et al., 2005), Eudicots (124 Mya; Hughes and McDougall, 1990). These calibration constraints were Raw Illumina reads for the four peach cultivars (Bigtop, used with soft bounds (Yang and Rannala, 2014) under a birth– Earlygold, Platurno and Sweetdream), the 10 almond © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 469 cultivars (A€ı, Belle d’Aurons, Cristomorto, Desmayo lar- Figure S3. Synteny analysis of almond versus peach genomes. gueta, Falsa Barese, Genco, Marcona, Nonpareil, Ripon, Figure S4. Distribution of recombination along chromosomes in Vivot) and one P. webbii are available at the European almond. Nucleotide Archive (ENA) under the study with a primary Figure S5. Distribution of size of in-paralog groups resulting from species-specific duplications. accession PRJEB32985 and corresponding experiment IDs Figure S6. Circos graphical representation of single nucleotide ERX3390856-ERX3390868 (10 almond cultivars; there are polymorphism and insertion/deletion distribution across the two entries for Ai, Belle d’Aurons and Desmayo largueta almond genome. due to double Illumina runs), ERX3391776-ERX3391779 Figure S7. Percentages of non-transposable element (TE) and TE (four peach cultivars) and ERX3391780 for P. webbii. Raw events for the different deletions in 10 almond varieties and one Illumina data for the seven peach cultivars (Armking, Belbi- peach variety. nette, Blanvio, Catherine, Flatmoon, Nectalady and Tiffany) Figure S8. Single nucleotide polymorphism-based phylogenetic analysis of 10 almond and cultivars one peach (EG2) cultivar. were downloaded from SRA, corresponding to the SRA Figure S9. Insertion time distribution of individual long terminal accessions ERS1801609-ERS1801614 and ERS1801617. repeat retrotransposon families of the Copia and Gypsy superfam- Raw Illumina data for almond cultivars D05-187 and S3067 ilies or that remained unclassified. were downloaded from SRA with the corresponding acces- Figure S10. Insertion time distribution of fixed and polymorphic sion IDs of SRX245830 and SRX245832. Sequencing reads long terminal repeat retrotransposon insertions in peach and and assembly data of P. dulcis cv. Texas are available via almond. the ENA (PRJEB32994). The assembly and annotation are Figure S11. Insertion time distribution of new and orthologous additionally accessible via the Genome Database for Rosa- long terminal repeat retrotransposon insertions in peach and almond. ceae (https://www.rosaceae.org/analysis/295) and the Figure S12. A stacked histogram based on the 27-mer matrix of CNAG-CRG (http://denovo.cnag.cat/almond). the assembly and the paired-end Illumina libraries. ACKNOWLEDGEMENTS Figure S13. Protein-coding gene annotation pipeline. Table S1. Summary of sequence data used for Texas almond gen- This research was supported in part by grants from the Spanish ome sequencing. Ministry of Economy and Competitiveness (MINECO/FEDER pro- jects AGL2012-40228-C02-01, AGL2015-68329-R, AGL2016-78992-R Table S2. Comparison between the Prunus dulcis cv. Texas gen- ome sequence assembly and annotation statistics and that of cv. and RTA2015-00050-00-00), Severo Ochoa Program for Centres of Excellence in R&D 201-2019 SEV-2015-0533 and CERCA Pro- Lauranne obtained by Sanchez-Perez et al. (2019). gramme-Generalitat de Catalunya. MW acknowledges grant support Table S3. Mapping of single nucleotide polymorphism markers from the Waite Research Institute of the University of Adelaide and from the Texas almond 9 Earlygold peach linkage map onto the MJR from the Spanish Instituto Nacional de Investigacio  n y Tec- almond assembly. nologıa Agraria y Alimentaria (INIA) project no. RTA-2014-00062. Table S4. List of species used in the phylome reconstruction. Table S5. Estimated dates (million years ago) and 95% highest AUTHOR CONTRIBUTIONS posterior density. PA, TA, JMC and TG designed the project and supervised Table S6. List of the Gene Ontology terms enriched in protein fam- ilies of almond and peach that duplicated at the last common research. TA, FC, LF, JGG, MG, PR and AV performed the ancestor of Prunus species. genome assembly and annotation. JMC, AB, FB, RC, JM Table S7. List of the Gene Ontology terms enriched in the protein provided the analysis of transposable elements. TG and IJ families lost specifically in peach and almond. performed the phylome analysis. KA and PA performed the Table S8. Almond cultivars selected and their main characteristics. variability and synteny analysis. AD, HD, AFM, MJRC, MW Table S9. Mapping statistics for the re-sequenced almond cultivars. and WH collected materials, extracted DNA and provided Table S10. Variant distribution across the almond pseudo- re-sequencing data. JLG and BG generated the fosmid molecules. library. PP and TA provided long-range DNA sequencing. Table S11. Comparison of single nucleotide polymorphism vari- PA, TA, JMC, JGM, TG and KA wrote the paper. ability parameters in Prunus species with whole genome sequences available. CONFLICT OF INTEREST Table S12. Deletions in 10 almond cultivars and one peach cultivar compared with the almond reference sequence and deletions that The authors declare no conflicts of interest. contain transposable element sequences. SUPPORTING INFORMATION Table S13. Summary of variants detected between Prunus dulcis and Prunus persica assemblies. Additional Supporting Information may be found in the online ver- Table S14. General statistics of transposable element annotation sion of this article. in Prunus dulcis and Prunus persica. Table S15. Percentage of transposable element coverage at the Figure S1. Genomescope k-mer coverage model fit using 21-mers order level in Prunus dulcis and Prunus persica. of the Texas almond genome. Table S16. Detailed annotation of long terminal repeat retrotrans- Figure S2. Synteny between the almond genome and the Texas posons and miniature inverted-repeat transposable elements in almond 9 Earlygold peach linkage map. Prunus dulcis and Prunus persica. © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 470 Tyler Alioto et al. Table S17. List of the 97 genes potentially involved in mesocarp Delcher, A.L., Salzberg, S.L. and Phillippy, A.M. (2003) Using MUMmer to development. identify similar regions in large sequence sets. Curr. Protoc. Bioinformat- ics, 10–13. Table S18. Methylation status on genes potentially involved in Delplancke, M., Yazbek, M., Arrigo, N., Espındola, A., Joly, H. and Alvarez, mesocarp development and presenting transposable element N. (2016) Combining conservative and variable markers to infer the evo- insertions in peach or almond. lutionary history of Prunus subgen. Amygdalus s.l. under domestication. Genet. Resour. Crop Evol. 63, 221–234. REFERENCES Diamond, J. (1997) Guns, germs, and steel: the fates of human societies. New York: W.W. Norton & Company, p. 14. Abascal, F., Corvelo, A., Cruz, F. et al. (2016) Extreme genomic erosion after Dirlewanger, E., Graziano, E., Joobeur, T., Garriga-Caldere, F., Cosson, P., recurrent demographic bottlenecks in the highly endangered Iberian Howad, W. and Arus, P. (2004) Comparative mapping and marker-as- lynx. Genome Biol. 17, 251. sisted selection in Rosaceae fruit crops. Proc. Natl Acad. Sci. USA, 101, Alioto, T., Blanco, E., Parra, G. and Guigo , R. (2018) Using geneid to identify 9891–9896. genes. Curr. Protoc. Bioinformatics, 64, e56. Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Al-Shahrour, F., Minguez, P., Tarraga, J., Medina, I., Alloza, E., Montaner, D. Batut, P., Chaisson, M. and Gingeras, T.R. (2013) STAR: ultrafast univer- and Dopazo, J. (2007) FatiGO +: a functional profiling tool for genomic sal RNA-seq aligner. Bioinformatics, 29, 15–21. data. Integration of functional annotation, regulatory motifs and interac- Donoso, J.M., Eduardo, I., Pican ~ ol, R., Batlle, I., Howad, W., Aranzana, M.J. tion data with microarray experiments. Nucleic Acids Res. 35, W91–6. and Aru  s, P. (2015) High-density mapping suggests cytoplasmic male Altschul, S.F., Gish, W., Miller, W., Myers, E.W. and Lipman, D.J. (1990) sterility with two restorer genes in almond 9 peach progenies. Hortic. Basic local alignment search tool. J. Mol. Biol. 215, 403–410. Res. 2, 15016. Aranzana, M.J., Decroocq, V., Dirlewanger, E. et al. (2019) Prunus genetics Donoso, J.M., Pican~ ol, R., Serra, O., Howad, W., Alegre, S., Ar us, P. and and applications after de novo genome sequencing: achievements and Eduardo, I. (2016) Exploring almond genetic variability useful for peach prospects. Hortic. Res. 6, 58. improvement: mapping major genes and QTLs in two interspecific Badenes, M.L. and Parfitt, D.E. (1995) Phylogenetic relationships of culti- almond 9 peach populations. Mol. Breed. 36, 1–17. vated Prunus species from an analysis of chloroplast DNA variation. The- Edgar, R.C. (2004) MUSCLE: multiple sequence alignment with high accu- or. Appl. Genet. 90, 1035–1041. racy and high throughput. Nucleic Acids Res. 32, 1792. Baek, S., Choi, K., Kim, G.B. et al. (2018) Draft genome sequence of wild Ellinghaus, D., Kurtz, S. and Willhoeft, U. (2008) LTRharvest, an efficient Prunus yedoensis reveals massive inter-specific hybridization between and flexible software for de novo detection of LTR retrotransposons. sympatric flowering cherries. Genome Biol. 19, 127. BMC Bioinformatics, 9, 18. Baird, W.V., Estager, A.S. and Wells, J.K. (1994) Estimating nuclear DNA Falchi, R., Vendramin, E., Zanon, L., Scalabrin, S., Cipriani, G., Verde, I., Viz- content in peach and related diploid species using laser flow cytometry zotto, G. and Morgante, M. (2013) Three distinct mutational mechanisms and DNA hybridization. J. Am. Soc. Hortic. Sci. 119, 1312–1316. acting on a single gene underpin the origin of yellow flesh in peach. Bell, C.D., Soltis, D.E. and Soltis, P.S. (2010) The age and diversification of Plant J. 76, 175–187. the angiosperms re-revisited. Am. J. Bot. 97, 1296–1303. Finn, R.D., Coggill, P., Eberhardt, R.Y. et al. (2016) The Pfam protein families Berthelier, J., Casse, N., Daccord, N., Jamilloux, V., Saint-Jean, B. and Car- database: towards a more sustainable future. Nucleic Acids Res. 44, rier, G. (2018) A transposable element annotation pipeline and expres- D279–D285. sion analysis reveal potentially active elements in the microalga Flutre, T., Duprat, E., Feuillet, C. and Quesneville, H. (2011) Considering Tisochrysis lutea. BMC Genomics, 19, 378. transposable element diversification in de novo annotation approaches. Boetzer, M. and Pirovano, W. (2014) SSPACE-LongRead: scaffolding bacte- PLoS ONE, 6, e16526. rial draft genomes using long read sequence information. BMC Bioinfor- Foster, T.M. and Aranzana, M.J. (2018) Attention sports fans! The far-reach- matics, 15, 211. ing contributions of bud sport mutants to horticulture and plant biology. Bolger, A.M., Lohse, M. and Usadel, B. (2014) Trimmomatic: a flexible trim- Hortic. Res. 5, 44. mer for Illumina sequence data. Bioinformatics, 30(15), 2114–2120. Fu, L., Niu, B., Zhu, Z., Wu, S. and Li, W. (2012) CD-HIT: accelerated for clus- Browicz, K. and Zohary, D. (1996) The genus Amygdalus L. (Rosaceae): spe- tering the next-generation sequencing data. Bioinformatics, 28, 3150– cies relationships, distribution and evolution under domestication. 3152. Genet. Resour. Crop Evol. 43, 229–247. Gabaldo  n, T. (2008) Large-scale assignment of orthology: back to phyloge- Callahan, A., Dardick, C., Tosetti, R., Lalli, D. and Scorza, R. (2015) 21st cen- netics? Genome Biol. 9, 235. tury approach to improving Burbank’s ‘stoneless’ plum. HortScience, 50, Gaut, B.S. (2015) Evolution is an experiment: assessing parallelism in crop 195–200. domestication and experimental evolution. Mol. Biol. Evol. 32, 1661–1671. Capella-Gutierrez, S., Silla-Martınez, J.M. and Gabaldo  n, T. (2009) trimAl: a Guindon, S. and Gascuel, O. (2003) A simple, fast, and accurate algorithm tool for automated alignment trimming in large-scale phylogenetic anal- to estimate large phylogenies by maximum likelihood. Syst. Biol. 52, yses. Bioinformatics, 25, 1972–1973. 696–704. Casacuberta, J.M. and Santiago, N. (2003) Plant LTR-retrotransposons and Guindon, S., Dufayard, J.F., Lefort, V., Anisimova, M., Hordijk, W. and Gas- MITEs: control of transposition and impact on the evolution of plant cuel, O. (2010) New algorithms and methods to estimate maximum-likeli- genes and genomes. Gene, 311, 1–11. hood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. Chen, J., Hu, Q., Zhang, Y., Lu, C. and Kuang, H. (2014) P-MITE: a database 59, 307–321. for plant miniature inverted-repeat transposable elements. Nucleic Acids Haas, B.J., Salzberg, S.L., Zhu, W., Pertea, M., Allen, J.E., Orvis, J., White, Res. 42, D1176–D1181. O., Robin, C.R. and Wortman, J.R. (2008) Automated eukaryotic gene Crepet, W. and Nixon, K. (1996)The fossil history of stamens. In D’Arcy, W. structure annotation using EVidenceModeler and the program to assem- and Keating, R., eds. The anther: form, function and phylogeny. Cam- ble spliced alignments. Genome Biol. 9, R7. bridge: Cambridge University Press, pp. 25–28. Han, Y. and Wessler, S.R. (2010) MITE-Hunter: a program for discovering  mez-Garrido, J. et al. (2016) Genome sequence of the Cruz, F., Julca, I., Go miniature inverted-repeat transposable elements from genomic olive tree, Olea europaea. Gigascience, 5, 29. sequences. Nucleic Acids Res. 38, e199. Cui, X., Lu, Z., Wang, S., Jing-Yan Wang, J. and Gao, X. (2016) CMsearch: Heppner, M.J. (1923) The factor for bitterness in the sweet almond. Genet- simultaneous exploration of protein sequence space and structure space ics, 8, 390–392. improves not only protein homology detection but also protein structure Herendeen, P.S., Crane, P.R. and Drinnan, A.N. (2002) Fagaceous flowers, prediction. Bioinformatics, 32, i332–i340. fruits, and cupules from the Campanian (Late Cretaceous) of Central Darling, A.C.E., Mau, B., Blattner, F.R. and Perna, N.T. (2004) Mauve: multi- Georgia, USA. Int. J. Plant Sci. 156, 93–116. ple alignment of conserved genomic sequence with rearrangements. Hoede, C., Arnoux, S., Moisset, M., Chaumier, T., Inizan, O., Jamilloux, V. Genome Res. 14, 1394–1403. and Quesneville, H. (2014) PASTEC: an automatic transposable element Darwin, C. (1968) The variation of animals and plants under domestication. classification tool. PLoS One, 9, e91929. New York: John Murray. © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 Almond genome 471 Huerta-Cepas, J. and Gabaldo n, T. (2011) Assigning duplication events to Mapleson, D., Accinelli, G.G., Kettleborough, G., Wright, J. and Clavijo, B.J. relative temporal scales in genome-wide studies. Bioinformatics, 27, 38– (2017) KAT: a K-mer analysis toolkit to quality control NGS datasets and 45. genome assemblies. Bioinformatics, 33, 574–576. Huerta-Cepas, J., Dopazo, J. and Gabaldo n, T. (2010) ETE: a Python environ- Marcßais, G. and Kingsford, C. (2011) A fast, lock-free approach for efficient ment for tree exploration. BMC Bioinformatics, 11, 24. parallel counting of occurrences of k-mers. Bioinformatics, 27, 764–770. Huerta-Cepas, J., Capella-Gutierrez, S., Pryszcz, L.P., Denisov, I., Kormes, Marchler-Bauer, A., Lu, S., Anderson, J.B. et al. (2011) CDD: a conserved D., Marcet-Houben, M. and Gabaldo  n, T. (2011) PhylomeDB v3.0: an domain database for the functional annotation of proteins. Nucleic Acids expanding repository of genome-wide collections of trees, alignments Res. 39, D225–D229. and phylogeny-based orthology and paralogy predictions. Nucleic Acids Martin, A., Troadec, C., Boualem, A., Rajab, M., Fernandez, R., Morin, H., Res. 39, D556–60. Pitrat, M., Dogimont, C. and Bendahmane, A. (2009) A transposon-in- Huerta-Cepas, J., Capella-Gutierrez, S., Pryszcz, L.P., Marcet-Houben, M. duced epigenetic change leads to sex determination in melon. Nature, and Gabaldo  n, T. (2014) PhylomeDB v4: zooming into the plurality of 461, 1135–1138. evolutionary histories of a genome. Nucleic Acids Res. 42, D897–D902. Mayjonade, B., Gouzy, J., Donnadieu, C., Pouilly, N., Marande, W., Callot, Hughes, N.F. and McDougall, A.B. (1990) Barremian–Aptian angiospermid C., Langlade, N. and Mun ~ os, S. (2016) Extraction of high-molecular- pollen records from southern England. Rev. Palaeobot. Palynol. 65, 145– weight genomic DNA for long-read sequencing of single molecules. 151. Biotechniques, 61, 203–205. Hunter, S., Jones, P., Mitchell, A. et al. (2012) InterPro in 2011: new develop- Morata, J., Marın, F., Payet, J. and Casacuberta, J.M. (2018) Plant lineage- ments in the family and domain prediction database. Nucleic Acids Res. specific amplification of transcription factor binding motifs by miniature 40, D306–D12. inverted-repeat transposable elements (MITEs). Genome Biol. Evol. 10, Jin, J., Tian, F., Yang, D.C., Meng, Y.Q., Kong, L., Luo, J. and Gao, G. (2017) 1210–1220. PlantTFDB 4.0: toward a central hub for transcription factors and regula- Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A.C. and Kanehisa, M. (2007) tory interactions in plants. Nucleic Acids Res. 45, D1040–D1045. KAAS: an automatic genome annotation and pathway reconstruction ser- Johnson, L.S., Eddy, S.R. and Portugaly, E. (2010) Hidden Markov model ver. Nucleic Acids Res. 35, W182–W185. speed heuristic and iterative HMM search procedure. BMC Bioinformat- Nattestad, M. and Schatz, M.C. (2016) Assemblytics: a web analytics tool ics, 11, 431. for the detection of variants from an assembly. Bioinformatics, 32, 3021– Joobeur, T., Viruel, M.A., de Vicente, L.M. et al. (1998) Construction of a sat- 3023. urated linkage map for Prunus using an almond 9 peach F2 progeny. Nawrocki, E.P. and Eddy, S.R. (2013) Infernal 1.1: 100-fold faster RNA Theor. Appl. Genet. 97, 1034–1041. homology searches. Bioinformatics, 29, 2933–2935. Kanehisa, M., Goto, S., Sato, Y., Furumichi, M. and Tanabe, M. (2012) KEGG Nawrocki, E.P., Burge, S.W., Bateman, A. et al. (2015) Rfam 12.0: updates to for integration and interpretation of large-scale molecular data sets. the RNA families database. Nucleic Acids Res. 43, D130–D137. Nucleic Acids Res. 40, D109–D114. Pertea, M., Pertea, G.M., Antonescu, C.M., Chang, T.C., Mendell, J.T. and Katoh, K., Asimenos, G. and Toh, H. (2009) Multiple alignment of DNA Salzberg, S.L. (2015) StringTie enables improved reconstruction of a sequences with MAFFT bioinformatics for DNA sequence analysis. Meth- transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. ods Mol. Biol. 537, 39–64. Petersen, T.N., Brunak, S., von Heijne, G. and Nielsen, H. (2011) SignalP 4.0: Kester, D.E., Gradziel, T.M. and Grasselly, C. (2015) Almonds (Prunus). Acta discriminating signal peptides from transmembrane regions. Nat. Meth- Hortic. 701–760. ods, 8, 785–786. Kislev, M.E., Nadel, D. and Carmi, I. (1992) Epipalaeolithic (19,000 BP) cereal Pryszcz, L.P. and Gabaldo  n, T. (2016) Redundans: an assembly pipeline for and fruit diet at Ohalo II, Sea of Galilee, Israel. Rev. Palaeobot. Palynol. highly heterozygous genomes. Nucleic Acids Res. 44, e113. 73, 161–166. Quesneville, H., Bergman, C.M., Andrieu, O., Autard, D., Nouaud, D., Ash- Koren, S., Walenz, B.P., Berlin, K., Miller, J.R., Bergman, N.H. and Phillippy, burner, M. and Anxolabehere, D. (2005) Combined evidence annotation A.M. (2017) Canu: scalable and accurate long-read assembly via adaptive of transposable elements in genome sequences. PLoS Comput. Biol. 1, k-mer weighting and repeat separation. Genome Res. 27, 722–736. 0166–0175. Krueger, F. and Andrews, S.R. (2011) Bismark: a flexible aligner and methy- Quinlan, A.R. and Hall, I.M. (2010) BEDTools: a flexible suite of utilities for lation caller for Bisulfite-Seq applications. Bioinformatics, 27(11), 1571– comparing genomic features. Bioinformatics, 26, 841–842. 1572. Rice, P., Longden, L. and Bleasby, A. (2000) EMBOSS: the European molecu- Krzywinski, M., Schein, J., Birol, I., Connors, J., Gascoyne, R., Horsman, D., lar biology open software suite. Trends Genet. 16, 276–277. Jones, S.J. and Marra, M.A. (2009) Circos: an information aesthetic for Rodriguez, G.R., Munos, S., Anderson, C., Sim, S.-C., Michel, A., Causse, comparative genomics. Genome Res. 19, 1639–1645. M., Gardener, B.B.M., Francis, D. and van der Knaap, E. (2011) Distribu- Ladizinsky, G. (1999) On the origin of almond. Genet. Resour. Crop Evol. 46, tion of SUN, OVATE, LC, and FAS in the tomato germplasm and the rela- 143–147. tionship to fruit shape diversity. Plant Physiol. 156, 275–285. Lartillot, N., Rodrigue, N., Stubbs, D. and Richer, J. (2013) Phylobayes mpi: Sanchez-Perez, R., Howad, W., Dicenta, F., Aru s, P. and Martınez-Go  mez, P. phylogenetic reconstruction with infinite mixtures of profiles in a parallel (2007) Mapping major genes and quantitative trait loci controlling agro- environment. Syst. Biol. 62, 611–615. nomic traits in almond. Plant Breed. 126, 310–318. Lassmann, T. and Sonnhammer, E.L.L. (2005) Kalign – an accurate and fast Sanchez-Perez, R., Pavan, S., Mazzeo, R. et al. (2019) Mutation of a bHLH multiple sequence alignment algorithm. BMC Bioinformatics, 6, 298. transcription factor allowed almond domestication. Science, 364, 1095– Li, H. (2011) A statistical framework for SNP calling, mutation discovery, 1098. association mapping and population genetical parameter estimation SanMiguel, P., Gaut, B.S., Tikhonov, A., Nakajima, Y. and Bennetzen, J.L. from sequencing data. Bioinformatics, 27, 2987–2993. (1998) The paleontology of intergene retrotransposons of maize. Nat. Li, H. and Durbin, R. (2009) Fast and accurate short read alignment with Bur- Genet. 20, 43–45. rows–Wheeler transform. Bioinformatics, 25, 1754–1760. Sanseverino, W., Hermoso, A., D’Alessandro, R., Vlasova, A., Andolfo, G., Li, Y., Smith, T., Liu, C.J., Awasthi, N., Yang, J., Wang, Y.F. and Li, C.Sen Frusciante, L., Lowy, E., Roma, G. and Ercolano, M.R. (2013) PRGdb 2.0: (2011) Endocarps of Prunus (Rosaceae: Prunoideae) from the early towards a community-based database model for the analysis of R-genes Eocene of Wutu, Shandong Province, China. Taxon 60, 555–564. in plants. Nucleic Acids Res. 41, D1167–D1171. Lomsadze, A., Burns, P.D. and Borodovsky, M. (2014) Integration of mapped Scholz, A., Rabaey, D., Stein, A., Cochard, H., Smets, E. and Jansen, S. RNA-Seq reads into automatic training of eukaryotic gene finding algo- (2013) The evolution and function of vessel and pit characters with rithm. Nucleic Acids Res. 42, e119. respect to cavitation resistance across 10 Prunus species. Tree Physiol. Lowe, T. (1997) tRNAscan-SE: a program for improved detection of transfer 33, 684–694. RNA genes in genomic sequence. Nucleic Acids Res. 25, 955–964. Sedlazeck, F.J., Rescheneder, P., Smolka, M., Fang, H., Nattestad, M., von Mao, H., Wang, H., Liu, S., Li, Z., Yang, X., Yan, J., Li, J., Tran, L.S.P. and Haeseler, A. and Schatz, M.C. (2018) Accurate detection of complex Qin, F. (2015) A transposable element in a NAC gene is associated with structural variations using single-molecule sequencing. Nat. Methods, drought tolerance in maize seedlings. Nat. Commun. 6, 8326. 15, 461–468. © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472 472 Tyler Alioto et al. Shaw, J. and Small, R.L. (2004) Addressing the ‘hardest puzzle in American Verde, I., Jenkins, J., Dondini, L. et al. (2017) The peach v2.0 release: high- pomology:’ phylogeny of Prunus sect. Prunocerasus (Rosaceae) based resolution linkage mapping and deep resequencing improve chromo- on seven noncoding chloroplast DNA regions. Am. J. Bot. 91, 985–996. some-scale assembly and contiguity. BMC Genomics, 18, 225. Shirasawa, K., Isuzugawa, K., Ikenaga, M., Saito, Y., Yamamoto, T., Hira- Voorrips, R.E. (2002) MapChart: software for the graphical presentation of kawa, H. and Isobe, S. (2017) The genome sequence of sweet cherry linkage maps and QTLs. J. Hered. 93, 77–78. (Prunus avium) for use in genomics-assisted breeding. DNA Res. 24, Vurture, G.W., Sedlazeck, F.J., Nattestad, M., Underwood, C.J., Fang, H., 499–508. Gurtowski, J. and Schatz, M.C. (2017) GenomeScope: fast reference-free Sima~o, F.A., Waterhouse, R.M., Ioannidis, P., Kriventseva, E.V. and Zdob- genome profiling from short reads. Bioinformatics, 2202–2204. nov, E.M. (2015) BUSCO: assessing genome assembly and annotation Wallace, I.M., O’Sullivan, O., Higgins, D.G. and Notredame, C. (2006) M-Cof- completeness with single-copy orthologs. Bioinformatics, 31, 3210–3212. fee: combining multiple sequence alignment methods with T-Coffee. Simpson, J.T., Wong, K., Jackman, S.D., Schein, J.E. Jones, S.J.M. and Nucleic Acids Res. 34, 1692–1699. Birol, I. (2009) ABySS: A parallel assembler for short read sequence data. Wehe, A., Bansal, M.S., Burleigh, J.G. and Eulenstein, O. (2008) DupTree: a Genome Res. 19(6), 1117–1123. program for large-scale phylogenetic analyses using gene tree parsi- Slater, G.S.C. and Birney, E. (2005) Automated generation of heuristics for mony. Bioinformatics, 24, 1540–1541. biological sequence comparison. BMC Bioinformatics, 6, 31. Wickson, E.J. (1889) The California fruits and how to grow them. San Fran- Soderlund, C., Nelson, W., Shoemaker, A. and Paterson, A. (2006) SyMAP: a cisco: Dewey & Co. system for discovering and viewing syntenic regions of FPC maps. Gen- Wikstro€ m, N., Savolainen, V. and Chase, M.W. (2001) Evolution of the ome Res. 16, 1159–1168. angiosperms: calibrating the family tree. Proc. R. Soc. B Biol. Sci. 268, Sommer, D.D., Delcher, A.L., Salzberg, S.L. and Pop, M. (2007) Minimus: a 2211–2220. fast, lightweight genome assembler. BMC Bioinformatics, 8, 64. Wu, T.D. and Watanabe, C.K. (2005) GMAP: a genomic mapping and align- Spiegel-Roy, P. (1976) L’amandier en Israel. Options Me ditetrrane ennes, 32, ment program for mRNA and EST sequences. Bioinformatics, 21, 1859– 92–95. 1875. Stanke, M., Keller, O., Gunduz, I., Hayes, A., Waack, S. and Morgenstern, B. Xiang, Y., Huang, C.H., Hu, Y., Wen, J., Li, S., Yi, T., Chen, H., Xiang, J. and (2006) AUGUSTUS: ab initio prediction of alternative transcripts. Nucleic Ma, H. (2017) Evolution of Rosaceae fruit types based on nuclear phy- Acids Res. 34, W435–W439. logeny in the context of geological times and genome duplication. Mol. Su, T., Wilf, P., Huang, Y., Zhang, S. and Zhou, Z. (2015) Peaches preceded Biol. Evol. 34, 262–281. humans: fossil evidence from SW China. Sci. Rep. 5, 16794. Xiong, Y. and Eickbush, T.H. (1990) Origin and evolution of retroelements Tang, H., Zhang, X., Miao, C., Zhang, J., Ming, R., Schnable, J.C., Schnable, based upon their reverse transcriptase sequences. EMBO J. 9, 3353–62. P.S., Lyons, E. and Lu, J. (2015) ALLMAPS: robust scaffold ordering Yang, Z. and Rannala, B. (2006) Bayesian estimation of species divergence based on multiple maps. Genome Biol. 16, 3. times under a molecular clock using multiple fossil calibrations with soft Tatsuki, M., Soeno, K., Shimada, Y. et al. (2018) Insertion of a transposon- bounds. Mol. Biol. Evol. 23, 212–226. like sequence in the 50 -flanking region of the YUCCA gene causes the Yazbek, M.M. and Al-Zein, M.S. (2014) Wild almonds gone wild: revisiting stony hard phenotype. Plant J. 96, 815–827. Darwin’s statement on the origin of peaches. Genet. Resour. Crop Evol. Thodberg, S., Cueto, J.Del, Mazzeo, R., Pavan, S., Lotti, C., Dicenta, F., 61, 1319–1328. Jakobsen Neilson, E.H., Møller, B.L. and Sanchez-Perez, R. (2018) Eluci- Ye, K., Schulz, M.H., Long, Q., Apweiler, R. and Ning, Z. (2009) Pindel: a pat- dation of the amygdalin pathway reveals the metabolic basis of bitter tern growth approach to detect break points of large deletions and med- and sweet almonds (Prunus dulcis). Plant Physiol. 178, 1096–1111. ium sized insertions from paired-end short reads. Bioinformatics, 25, Treangen, T.J., Sommer, D.D., Angly, F.E., Koren, S. and Pop, M. (2011) 2865–2871. Next generation sequence assembly with AMOS. Curr. Protoc. Bioinfor- Yu, Y., Fu, J., Xu, Y. et al. (2018) Genome re-sequencing reveals the evolu- matics, 1–18. tionary history of peach fruit edibility. Nat. Commun. 9, 5404. Vanholme, R., Demedts, B., Morreel, K., Ralph, J. and Boerjan, W. (2010) Zdobnov, E.M. and Apweiler, R. (2001) InterProScan – an integration plat- Lignin biosynthesis and structure. Plant Physiology, 153(3), 895–905. form for the signature-recognition methods in InterPro. Bioinformatics, Velasco, D., Hough, J., Aradhya, M. and Ross-Ibarra, J. (2016) Evolutionary 17, 847–848. genomics of peach and almond domestication. G3: Genes Genomes Zhang, Q., Chen, W., Sun, L. et al. (2012) The genome of Prunus mume. Genetics, 6, 3985–3993. Nat. Commun. 3, 1318. Vendramin, E., Pea, G., Dondini, L. et al. (2014) A unique mutation in a MYB Zhang, S.D., Jin, J.J., Chen, S.Y., Chase, M.W., Soltis, D.E., Li, H.T., Yang, gene cosegregates with the nectarine phenotype in peach. PLoS ONE, 9, J.B., Li, D.Z. and Yi, T.S. (2017) Diversification of Rosaceae since the Late e90574. Cretaceous based on plastid phylogenomics. New Phytol. 214, 1355– Verde, I., Abbott, A.G., Scalabrin, S. et al. (2013) The high-quality draft 1367. genome of peach (Prunus persica) identifies unique patterns of Zimin, A.V., Marcßais, G., Puiu, D., Roberts, M., Salzberg, S.L. and Yorke, genetic diversity, domestication and genome evolution. Nat. Genet. J.A. (2013) The MaSuRCA genome assembler. Bioinformatics, 29, 2669– 45, 487–494. 2677. © 2019 The Authors. The Plant Journal published by Society for Experimental Biology and John Wiley & Sons Ltd, The Plant Journal, (2020), 101, 455–472