A total of 124 plastomes were preliminary selected, among which 91 were newly assembled and 33 were downloaded from the National Center for Biotechnology Information (NCBI) database (Supplementary Table 2). Our samples covered 13 genera of subfamilies Amygdaloideae (tribes Amygdaleae, Exochordeae, Spiraeeae, and Maleae) and Rosoideae of family Rosaceae. Three species, Morus mongolica (Moraceae), Ziziphus jujuba (Rhamnaceae), and Elaeagnus macrophylla (Elaeagnaceae), were used as out-groups. Based on our previous studies (Huang et al., 2013; Chen et al., 2015, 2016a,b, 2020; Liu et al., 2016; Zhang et al., 2018), 90 representative true cherry (Cerasus) and dwarf cherry (Microcerasus) accessions were selected for whole-genome re-sequencing, consisting of 34 P. pseudocerasus accessions (11 landraces and 23 wild types) representing diverse genotypes, phenotypes and geographical distributions, 6 accessions referring to 4 European cherry taxa [P. avium, Prunus fruticosa, P. cerasus × Prunus canescens (Gisela 5), and Prunus mahaleb], 46 accessions covering 20 other Cerasus taxa, and 4 dwarf cherry accessions (P. tomentosa and P. tianshanica). We also obtained the genomic pair-end reads of Prunus cerasifera (SRR4036106) from the GenBank database to assemble its plastome.

Total genomic DNA was extracted from silica-gel dried leaf tissues following the modified cetyltrimethyl ammonium bromide (CTAB) method used by Chen et al. (2013). Ninety-one Illumina paired-end (PE) libraries with 500-bp insert size were constructed and sequenced using an Illumina HiSeq 2000 (Illumina, San Diego, CA, United States) instrument by BGI-Shenzhen (Shenzhen, China). Taking Prunus persica (Jansen et al., 2011) and P. pseudocerasus (Feng et al., 2017) as reference plastomes, we obtained plastid reads for these 91 accessions. These reads were then assembled into contigs and scaffolds using SPAdes v.3.9.0 (Bankevich et al., 2012). The scaffolds were aligned to the reference plastomes of P. persica in Geneious v.8.1 (Kearse et al., 2012) and then were manually ordered as the genomes in the SnapGene v.2.3.2 software². The newly assembled plastomes were deposited in the GenBank database under the following accession numbers: MT576845-MT576934 (Supplementary Table 2).

Gene annotation was conducted in GeSeq³ (Tillich et al., 2017). All ambiguous annotations, such as the absence of start/stop codons, were manually corrected in SnapGene, referring to the downloaded 33 Rosaceae plastomes. Genome structures were drawn with Circos v.0.69.6 (Krzywinski et al., 2009).

The plastome sequences were aligned in MAFFT v.7.037b (Katoh and Standley, 2013). Nucleotide contents and coefficients of sequence similarities were calculated with BioEdit v.7.0.5 (Hall, 1999). Genetic distances were calculated with Tajima-Nei model in MEGA v.5.1 (Tamura et al., 2011). Genetic differentiation coefficient (F_ST value) was estimated using DnaSP v.6 (Rozas et al., 2017).

Microsatellites in each plastome were screened with MISA perl script⁴ with the following parameters: mononucleotide SSR repeat units ≥ 10, dinucleotide repeat units ≥ 6, trinucleotide repeat units ≥ 5, tetranucleotide repeat units ≥ 4, and pentanucleotide and hexanucleotide repeat units ≥ 3.

Taking P. pseudocerasus (NC030599.1) as the reference plastome, we employed MUMmer v.3.3 (Kurtz et al., 2004) to identify SNP and InDels for each Rosaceae plastome. Independent SNP and InDel files from different individuals were transformed and combined into one SNP and InDel vcf file with BCFTools v.1.7. BEDTools v.2.26.0 (Quinlan and Hall, 2010) and SnpEff_latest_core⁵ were used to detect the distributions of SSRs, SNPs, and InDels across the plastomes and estimate the effects of SNPs and InDels on gene functions.

To calculate the proportion of mutational events (Gielly and Taberlet, 1994; Huang et al., 2014), we also detected SNPs and InDels between pairwise plastomes in MUMmer v.3.3. The proportion of mutation events (PME) was calculated as [(NS + NI)/LA] × 100, where NS represented the number of nucleotide substitution between plastomes, NI represented the number of InDels, and LA denoted the length of the aligned plastome sequences.

We employed site models and likelihood ratio test (LRT) implemented in PAML v.4.9h (Yang, 2007) to detect the signatures of positive selection for 81 unique plastid protein-coding genes. At the Rosaceae level, we removed rps19-fragment (not completely assembled in some taxa), rps12 (the special gene structure), and infA (pseudogene) from this analysis. Except for the three genes, 43 highly conserved plastid protein-coding genes have also been excluded for true cherries. Seventy-eight and 35 protein-coding genes remained for Rosaceae taxa and true cherries, respectively. First, selective pressures were computed in CodeML (included in PAML package) with three site models: M0 (model = 0, NS sites = 0), M1a (model = 0, NS sites = 1), and M2a (model = 0, NS sites = 1) (Nielsen and Yang, 1998). Then, likelihood ratio test was conducted to compare M1a against M2a by calculating the χ² critical value and P value. Finally, when the log likelihoods between the two models were statistically different (P < 0.05 in LRT), positively selected sites of genes were identified by Bayes empirical Bayes (BEB) analysis (posterior probabilities for site class > 0.95 and ϖ > 1) (Zhang et al., 2005) in the CodeML program. In addition, a branch-site model (Yang et al., 2005; Yang, 2007) was also used to investigate branch-specific selection for true cherries. Likelihood ratio test for positive selection on each examined branch was compared model A (model = 2, NS sites = 2, fix_omega = 0, omega = 1.5) against null model (model = 2, NS sites = 2, fix_omega = 1, omega = 1). Also, positively selected sites were determined by LRT and BEB analyses.

Systematic errors are thought to mainly result from the inaccurate alignment caused by rapidly evolving sites and may lead to an incorrect tree with strong supports, while the removal of problematic regions is an effective method for improving the robustness of phylogenomic reconstruction (Rodríguez-Ezpeleta et al., 2007). To reduce potential systematic errors, we constructed 12 different datasets to carry out the phylogenomic analyses at the Rosaceae family level and at the tribe Amygdaleae level. Since gene order and content were highly conserved in the studied plastomes of Rosaceae and the three outgroups, the alignment could be straightforward. The 12 datasets were as follows: (i) WCGD (whole plastomes dataset, n = 124) and PCGD (Amygdaleae plastomes dataset, n = 107) were constructed with complete plastome sequences, both removing all missing data (N) and long insertion (> 50 bp) sequences that were only detected in an individual; (ii) WOID (whole one inverted-repeat dataset, n = 124) and POID (Amygdaleae one inverted-repeat dataset, n = 107) were generated using large single-copy (LSC), short single-copy (SSC), and one IR sequences, also both removing all missing data and long insertion (> 50 bp) sequences; (iii) VSWD (variant sites of whole plastomes dataset, n = 124) and VSPD (variant sites of Amygdaleae plastomes dataset, n = 107) were constructed with the variant sites (SNPs) of WCGD and PCGD using the custom bash script, respectively; (iv) WGSD (whole gene sequence dataset, n = 124) and PGSD (Amygdaleae gene sequence dataset, n = 107) contained sequences of 102 unique genes; (v) PCWD (protein-coding sequence of whole plastomes dataset, n = 124) and PCPD (protein-coding sequence of Amygdaleae plastomes dataset, n = 107) consisted of exon sequences of 72 unique protein-coding genes; and (vi) PWGD (pruned whole plastomes dataset, n = 124) and PPGD (pruned Amygdaleae plastomes dataset, n = 107) were generated by removing rapidly evolving sites, and large InDels and sequences with rich structural variation of WCGD and PCGD in GBlocks v.0.91b (Castresana, 2000) (parameters: minimum sequences per conserved position, 65; minimum sequences per flank position, 110 (PWGD)/100 (PPGD); maximum number of contiguous non-conserved positions, 8; minimum block length, 10; allowed gap positions, none).

We employed Maximum likelihood (ML) methods to generate phylogenetic trees for each dataset mentioned above. For each dataset, the best-fit model of nucleotide substitution was selected with jModelTest v.2.1.7 (Darriba et al., 2012) using the Akaike Information Criterion (AIC). Bayesian inference (BI) analyses were conducted with MrBayes v.3.2.6 (Ronquist et al., 2012). Two independent Markov chain Monte Carlo (MCMC) algorithm chains were carried out, and each of them ran with one cold and three heated chains for 12,000,000 generations and started with a random tree and sampling one tree every 100 generations. When the average standard deviation of split frequencies was below0.01 between the two runs, analyses were considered to reach stationarity. The first 25% of generations were treated as burn-in. ML analyses were performed using IQ-TREE v.1.5.5 (Nguyen et al., 2015) with 1,000 regular bootstrap replicates (-b 1,000). All the phylogenetic trees were edited and presented using iTOL v.5 (Letunic and Bork, 2021) and FigTree v.1.4.4.⁶

Divergence time was estimated with BEAST v.2.6.6 (Bouckaert et al., 2019). Since 8 (WCGD, WOID, PWGD, VSWD, PCGD, POID, PPGD, and VSPD) of the 12 datasets produced congruent topologies with high statistical supports, we used the PWGD in molecular dating following the method of Zhang et al. (2017). To further decrease computation power requirement, with the exception of P. pseudocerasus and P. avium, we kept only one accession for each taxon within inter nodes to construct a pruned PWGD. Software parameters were finally set as the GTR substitution model and exponential uncorrelated relaxed model with Yule process. Two independent MCMC runs were conducted, each with 300,000,000 generations and sampling every 1,000 generations. The first 12,000,000 generations in each run were removed as burn-in. The fossil P. wutuensis from Wutu Formation, Shandong province, China, has been dated to 47.8–55 Mya in Early Eocene (Li et al., 2011), and age estimates of Prunus was from 60.7 to 62.4 Mya (Chin et al., 2014; Table 1). Therefore, the age of crown Prunus (N1) was constrained by a log-normal distribution with a mean of 55 Mya and a standard deviation of 0.09 in our study. The divergence time of the tribe Maleae (only containing Malus, Pyrus, Chaenomeles, Cydonia, Docynia, Eriobotrya, and Sorbus) was estimated at approximately 42 Mya in recent molecular study (Xiang et al., 2017). Leaf fossils distinguishing Malus from Pyrus has been dated to 45 Mya (Wehr and Hopkins, 1994). Here, we set the age of crown Maleae (N2) as a log-normal distribution with a mean value of 45 Mya and a standard deviation of 0.01. Based on mesofossil (Edelman, 1975), ages of crown Rosoideae (only containing Rosa, Fragaria, and Potentilla) (N3) were constrained by the log-normal distribution with a mean value of 48.6 and standard deviation of 0.05. The recent age estimate of the divergence between Rosaceae and other Rosales taxa was at 106.5 Mya (Zhang et al., 2017). Thus, the crown Rosales (N4) was constrained by log-normal distribution with a mean value of 106.5 Mya and standard deviation of 0.05. Tree files and log files from the two independent runs were combined with LogCombiner v.2.6.6 (part of the BEAST package). The effective sample size (ESS) for each logged statistic was estimated in Tracer v.1.7.2 (Rambaut et al., 2018), and most of the ESSs were above 200. Finally, the consensus tree and divergence time were calculated and annotated in the TreeAnnotator v.2.6.6 (part of the BEAST package).

Ninety-one new plastomes of 27 Cerasus and Microcerasus taxa (25 species and two varieties) and one closely related species (P. cerasifera) were assembled (Supplementary Table 3). Mean coverage of these newly assembled plastomes ranged from 171 (P. pseudocerasus) to 9,065 × (P. cerasifera) (Supplementary Table 3). Thirty previously published plastomes of Rosaceae were also downloaded. The overall Guanine-Cytosine (GC) content of the 121 Rosaceae plastomes was ∼37% (Supplementary Table 4).

All the 121 examined Rosaceae plastomes exhibited typical genomic structures, consisting of one LSC and SSC, and two conserved inverted repeats (IRa and IRb) (Figure 1). They also possessed conserved gene orders and gene contents with 132 identified genes (115 unique genes), namely, 87 (81 unique) protein-coding genes, and 37 tRNA and 8 rRNA coding genes (Figure 1, Table 2, and Supplementary Table 4). Among these identified genes, nine protein-coding and five tDNA genes contained one intron, and three protein-coding and one tDNA genes contained more than one intron (Table 2). Most of the plastid genes were linearly concentrated on the plastomes (Figure 1), while the overlapping genomic regions that were detected between matK and trnU-UUU, between ycf1-fragment and ndhF, and between psbC and psbD. MatK were completely nested in the intron region of trnU-UUU in all the examined plastomes. The overlapping genomic region between ycf1-fragment and ndhF was observed in most plastomes of subfamily Amygdaloideae, with 21-bp length in the tribe Maleae and 4–173-bp length in the tribe Amygdaleae (Supplementary Table 4). A conserved overlapping genomic region (53-bp size) was detected between psbC and psbD in all the Rosaceae plastomes and three other Rosales plastomes (Supplementary Table 4).

We generated a total of 24 phylogenetic trees with 12 strictly and carefully proceeded datasets using both the ML and BI methods. Aligned sequences for each dataset ranged from 8,515 to 178,338 bp in length (Supplementary Table 17). The best-fit model GTR was set for VSWD and VSPD, and the best-fit model GTR + I + G for the remaining 10 datasets (Supplementary Table 17). The BI and ML trees generated with the same dataset had the same or highly congruent topologies; thus, we only presented the ML tree for each dataset (Figure 3 and Supplementary Figures 7–13). The topology of major clades and their ML bootstrap (BS) and BI posterior probabilities (PPs) are shown in Supplementary Figure 7. Except for the phylogenetic trees generating with exon sequences (PCWD), the remaining phylogenomic analyses suggested that Rosaceae formed four well-supported (78–100% BS, 1 PP) clades (A, B, C, and D) (Figure 3 and Supplementary Figure 7). Clade A included all examined tribe Amygdaleae taxa. Clade B only contained tribe Exochordeae species. Clade C contained taxa from tribes Spiraeeae and Maleae. Clade D corresponded to the three subfamily Rosoideae species.

In all of the ML and BI trees, clade A was further divided into three subclades (A_I, A_II, and A_III) with high support values (100% BS, 1 PP) (Figure 3 and Supplementary Figures 7–13). Prunus hypoleuca, P. serotina, and Prunus padus were assigned as basal subclade A_III. Microcerasus species were grouped with Armeniaca, Prunus, and the three Amygdalus taxa in subclade A_II. All Cerasus taxa formed a distinct monophyletic group (A_I). Within subclade A_I, ML and BI trees under eight datasets (WCGD, WOID, PWGD, VSWD, PCGD, POID, PPGD, and VSPD) showed six major groups (A_I1, A_I2, A_I3, A_I4, A_I5, and A_I6) with different support values (37–100% BS, 0.96–1 PP) (Figure 3 and Supplementary Figure 7). In these trees, A_I1 contained all cultivated P. pseudocerasus accessions, 11 wild P. pseudocerasus accessions mainly from Longmenshan Fault Zones, and 11 accessions of six Cerasus taxa. A_I2 and A_I4 consisted of 26 accessions of 9 Cerasus taxa and 26 accessions of 16 Cerasus taxa, respectively. A_I3 was composed of five European cherry accessions of P. avium, P. cerasus × P. canescens, and P. fruticosa. The groups A_I5 and A_I6 corresponded to P. cerasoides and P. mahaleb, respectively. In the eight datasets, difference was only observed in the BI tree under the VSWD (Supplementary Figure 7A), where A_I3 was further divided into two non-sister subgroups, A_{I3_1} (P. avium) and A_{I3_2} (P. cerasus × P. canescens and P. fruticosa). Further subdivision of A_I2, A_I3, and A_I4 was also observed in the phylogenetic trees generated with gene sequence (WGSD and PGSD) and exon sequence (PCWD and PCPD), while these subdivisions were weakly supported (9–73% BS) (Supplementary Figures 7B,C). Overall, despite the aforementioned differences among phylogenetic trees, resemblance between tree topology and taxonomy was completely lost, and no phylogeographic subdivision was detected among true cherries in all the ML and BI trees (Figure 3 and Supplementary Figures 7–13).

The divergence time between Rosaceae and other Rosales taxa was estimated to be 113.8 Mya (95% highest posterior density (HPD): 103.66–124.46 Mya), and the estimated divergence time for Rosaceae was 92.18 Mya (95% HPD: 81.74–108.71 Mya) (Table 4 and Figure 4). Maleae-Spiraeeae divergence occurred at 67.31 Mya (95% HPD: 58.28–78.36 Mya). The estimated time of divergence of the crown Maleae was 44.93 Mya (95% HPD: 44.07–45.83 Mya). The first divergence of the three subclades A_I (Cerasus), A_II (Microcerasus- Armeniaca- Prunus -Amygdalus), and A_III (Padus-Maddenia) occurred at 49.84 Mya (95% HPD: 42.39–57.82 Mya). Cerasus (true cherry) and relatives (Microcerasus, Armeniaca, Amygdalus, and Prunus) separated from each other at 28.21 Mya (95% HPD: 16.19–42.17 Mya). In subclade A_I, A_I6 (P. mahaleb) diverged from the common ancestor of true cherries at 15.28 Mya (95% HPD: 8.58–24.91 Mya), followed by A_I5 and A_I4 with an estimated divergence time of 11.16 (95% HPD: 6.42–17.87 Mya) and 9.51 Mya (95% HPD: 5.64–14.53 Mya), respectively. European cherry taxa (A_I3) and the cherry taxa of A_I1 and A_I2 diverged at 8.48 Mya (95% HPD: 4.95–13.01 Mya). A_I1 and A_I2 diverged from each other at 5.96 Mya (95% HPD: 3.06–9.73 Mya). In A_I1, the divergence of cultivated Chinese cherry, wild Chinese cherry, and relatives occurred between0.05 (95% HPD: 0–0.32 Mya) and0.83 Mya (95% HPD:0.26–1.96 Mya) (Figure 4).

Rosaceae plastomes exhibited conserved gene order and gene contents and showed a typical quadripartite structure (LSC, SSC, Ira, and IRb) that has been widely reported in green plants (Howe, 2016; Gao et al., 2019). The overlapping between matK and trnU-UUU, between ycf1-fragment and ndhF, and between psbC and psbD was also detected in most of the Rosaceae plastomes. Ultraconserved psbC-psbD overlapping regions (53 bp) were identified at the Rosales level. Combining with the results of Theaceae (52 bp) (Huang et al., 2014) and Malvaceae (53 bp) (Xu et al., 2012), we infer that the psbC-psbD overlapping region may have already been existing before the divergence of flowering plants.

Contraction or expansion of the IR region has been widely proposed to be the main reason for the variation in plastome size (Gao et al., 2019; Xue et al., 2019). Partial losses of plastid genes have been predicted as another main reason for the decreases of some Rosaceae plastomes (Xue et al., 2019). In this study, our data confirm that IR contraction, sequence losses in SSC region, and gene losses (atpF and rps19-fragment) in LSC regions are predominantly responsible for the relatively small sizes of subfamily Rosoideae plastomes. Our study also verifies that the plastome sizes of subfamily Amygdaloideae are particularly subject to the size increases/decreases in LSC and SSC regions.

The construction of the first accurate map of genomic variation (InDels and SNPs) has widely presented the hotspots across plastomes in Rosaceae. Most of the genomic variations were distributed in intergenic and intronic regions, which was consistent with those reported in rice (Gao et al., 2019), citrus (Bausher et al., 2006), tea (Huang et al., 2014), and other land plants (Howe, 2016; Dobrogojsk et al., 2020). High proportions of synonymous mutations and InDels with three multiple sizes exhibited slight or moderate influences on gene function and expression, suggesting strong constrains to maintain gene functions in Rosaceae plastomes. In addition, since the GC content can greatly influence gene expression (Barahimipour et al., 2016), the universally low GC contents in plastid protein-coding genes and intergenic regions of the Rosaceae plastomes may contribute to maintain efficient biological activities of chloroplasts while responding to diverse and extreme environment and climate changes.

We also detected strong signatures of positive selection in 11 plastid genes (rpoA, rps16, rps18, psaA, psbL, rbcL, ndhD, ndhF, accD, ycf1, and ycf2), especially in the rbcL gene that encodes the large subunit of Rubisco. The evolution of RuBisco large subunit has been thought to be driven by environmental pressures (Hermida-Carrera et al., 2017). Therefore, for Rosaceae crops, the positive selection in this gene probably contribute to their adaptation to various environmental stress (such as low CO₂ partial pressure) and climate shifts (Hermida-Carrera et al., 2017; Jiang et al., 2018; Yao et al., 2019; Shen et al., 2020). The NdhF gene has been suggested to poise the level of redox and consequently maintain or improve the photosynthetic performance of plant crops under extreme temperature and changing light intensity (Martín et al., 2009). The remaining nine genes also play crucial roles in chloroplast protein synthesis, energy transformation and regulation, and photosynthesis. These results indicate the diverse adaptive evolution in plastid genes of Rosaceae crops. Among the 11 plastid genes, 5 photosynthesis-related genes (rpoA, rps18, ndhF, rbcL, and ycf1) contained positively selected sites unique in most woody fruit trees. The positive selection in the five genes probably help Rosaceae woody fruit trees efficiently capture light energy to produce adequate nutrition to adapt to their growth and development under extreme and variable environmental conditions.

Tribes Amygdaleae and Maleae of the Rosaceae family contain many economically important crops that exhibit wide adaption to various environments and remarkably diversified phenotypes and genotypes; therefore, the origin, evolution, and domestication of the taxa of the two tribes have been widely explored by botanists and horticulturalists (Yamamoto and Terakami, 2016; Duan et al., 2017; Baek et al., 2018; Wu et al., 2018; Aranzana et al., 2019). Recent genomic and transcriptomic studies (Xiang et al., 2017; Yi et al., 2020) suggested that majority of tribe Amygdaleae members (x = 8) only underwent the ancient WGD shared by all Rosaceae members, while tribe Maleae members (x = 17) might have experienced two additional WGD events. These different evolution histories possibly contribute to the evolution of distinct fruit types of the two tribes (Xiang et al., 2017). Here, our plastome data also indicate completely different evolutionary patterns between the plastomes of tribes Amygdaleae and Maleae. In comparison with tribe Amygdaleae plastomes, tribe Maleae plastomes exhibited trivial genomic structural variation. This may suggest that the selective pressure from climatic fluctuations, environment changes, and human activities (e.g., domestication) probably results in slighter plastome variation in tribe Maleae than in tribe Amygdaleae.

The plastomes of the Amygdaleae tribe are highly varied in nucleotide substitution and genome structures, implying their highly divergent evolution. Our phylogenomic analyses demonstrated a well-supported relationship of (A_I (Cerasus) + A_II (Microcerasus, Armeniaca, Prunus, and Amygdalus)) + A_III (Maddenia and Padus) (Figure 3 and Supplementary Figures 7–13). This topology verifies the three lineages (Cerasus-Prunus-Padus) in previous molecular studies (Wen et al., 2008; Liu et al., 2013; Shi et al., 2013), where the Prunus in Cerasus-Prunus-Padus also only contained Microcerasus, Armeniaca, Prunus, and Amygdalus but excluded Cerasus of the broad-sensed genus Prunus. This study suggests that Cerasus forms a distinct group and that Microcerasus (dwarf cherry) is genetically closer to Armeniaca, Prunus, and Amygdalus than to Cerasus. These results can be supported by a recent whole-genome analysis (Wang et al., 2020).

It has been long controversial that Cerasus should be treated as a separate genus (de Tournefort, 1700; Linnaeus, 1754; Bate-Smith, 1961; Komarov, 1971; Shishkin and Yuzepchuk, 1971; Yü et al., 1986; Takhtajan, 1997), or as one of the subgenera of the broad-sensed genus Prunus (Bentham and Hooker, 1865; Focke, 1894; Schneider, 1905; Koehne, 1911; Rehder, 1940; Ingram, 1948; Hutchinson, 1964; Krüssmann, 1978; Ghora and Panigrahi, 1995). Here, the comparative genomic analyses have revealed that Cerasus (true cherry) contained diverse genomic variations and a unique nucleotide substitution pattern with transversion from G to T. The phylogenomic study showed that Cerasus was monophyletic and genetically distinct from relatives (Microcerasus, Armeniaca, Prunus, and Amygdalus). Molecular dating indicated that Cerasus and relatives diverged around the late Oligocene (28.21 (95% HPD: 16.19–42.17) Mya), a period before 66% angiosperm genera in China originated (∼23 Mya) (Lu et al., 2018). The level of plastome-wide divergence between Cerasus and relatives (Microcerasus, Armeniaca, Prunus, and Amygdalus) was even higher than those among genera of tribe Maleae. Moreover, Cerasus taxa show significant morphological differences in lenticels, axillary winter buds, petiole, and endocarp from relatives (Microcerasus, Armeniaca, Prunus, and Amygdalus) (Supplementary Table 26 and Supplementary Figure 14), and, generally, Cerasus taxa bear inflorescences umbellate or corymbose-racemose with moderate pedicels and conspicuous bracts, while relatives show solitary or two to three sessile flowers with absent bracts (Supplementary Table 26 and Supplementary Figure 14). Therefore, our results strongly support that Cerasus (true cherry) is treated as a separate genus, and this will be more convenient for horticulturist. As for the classification of groups Microcerasus, Armeniaca, Prunus, and Amygdalus, based on our present results, we prefer to treat them as subgenera or sections of the genus Prunus as in the previous study (Supplementary Table 1), while a further analysis of nuclear genome data is necessary.