1. Introduction

Angiosperms exhibit a great diversity of inflorescence types, such as racemes, panicles, corymbs, and cymes, depending on the number and arrangement of flowers on shoots, and their branching patterns (Benlloch et al., 2007; Endress, 2010). The evolution of inflorescences has attracted the interest of botanists for a long time (e.g., Parkin, 1914; Stebbins, 1973; Takhtajan, 1991; Endress, 2010). The rapid development of DNA sequencing technology, computational resources, and phylogenetic inference methods provide an opportunity for inferring accurate phylogenies and conducting rigorous tests of hypotheses concerning inflorescence evolution for some angiosperm groups based on genomic and phylogenetic data (Endress, 2010; Gerrath et al., 2017; Ma et al., 2017).

Prunus L. consists of 250–400 species of trees and shrubs widely distributed in the northern temperate zone and subtropical and tropical regions (Rehder, 1956; Yü et al., 1986; Lu et al., 2003), with eastern Asia being its center of diversity (Wen et al., 2008; Chin et al., 2014). This genus is characterized by simple leaves with stipules and leaf glands, a superior single ovary, and drupe fruits (Lu et al., 2003; Zhao et al., 2016). Many species of Prunus are economically significant as food crops, such as peach, plum, almond, and sweet cherry (Fig. 1; Bortiri et al., 2001; Lee and Wen, 2001), and many others have been used as ornamentals, timber, and medicine (Andro and Riffaud, 1995; Lee and Wen, 2001; Wen et al., 2008). Species of Prunus show a diversity of floral displays with several distinct inflorescence types, including racemes, corymbs, and solitary flowers, making it an ideal group to investigate the evolutionary transitions among different inflorescence organizations.

Rehder (1956) divided Prunus into five subgenera, namely Amygdalus L., Cerasus Mill., Lauro-cerasus Tourn. ex Duhamel, Padus Mill. and Prunus, and several taxonomists have accepted this classification (Kalkman, 2004; Wen et al., 2008; Chin et al., 2010, 2014; Zhao et al., 2016). Over the last two decades, researchers have tried to reconstruct the phylogenetic relationships of Prunus using different genomic regions, such as plastid markers, nuclear ribosomal ITS, and other nuclear loci (Bortiri et al., 2001, 2002, 2006; Lee and Wen, 2001; Wen et al., 2008; Chin et al., 2010, 2014; Shi et al., 2013; Zhao et al., 2016, 2018; Hodel et al., 2021). These analyses have recovered three main groups within Prunus, each of them presenting a synapomorphic inflorescence structure: (1) the solitary-flower group consisting of Prunus subg. Amygdalus (L.) Foche and P. subg. Prunus (including section Armeniaca (Scop.) Nakai); (2) the corymbose group, which includes P. subg. Cerasus (Mill.) Pers; and (3) the racemose group comprising P. subg. Lauro-cerasus and P. subg. Padus, as well as the Maddenia Hook. f. & Thomson and Pygeum Gaertn. groups (Chin et al., 2014). Most species of the solitary-flower and corymbose groups are diploid, whereas taxa of the racemose group usually have higher ploidy levels (Chin et al., 2014; Zhao et al., 2016; Hodel et al., 2021; CCDB, http://ccdb.tau.ac.il/). Nevertheless, the frequent, ancient hybridization events in the evolutionary history of Prunus have greatly challenged the clarification of the phylogenetic relationships among these three groups. The monophyly of the racemose group has been supported by plastid sequences (e.g., Chin et al., 2010, 2014); however, most analyses of nuclear sequence data have resolved the racemose group as paraphyletic (e.g., Bortiri et al., 2001, 2002; Lee and Wen, 2001; Chin et al., 2014), leading to hypotheses of multiple hybrid origins of this mostly polyploid lineage (Chin et al., 2014; Zhao et al., 2016). The monophyly of the racemose group has been strongly supported by hundreds of single-copy nuclear genes and chloroplast genomes from 21 transcriptomic data points (Hodel et al., 2021); however, only two species of the racemose group were included in their study, questioning the representativeness of the sampling regarding the taxonomic and morphological diversity in this lineage. Therefore, these results need to be further tested using more comprehensive taxon sampling.

Due to limited taxon sampling and informative loci in previous molecular phylogenetic studies of Prunus, the phylogenetic relationships among and within the major lineages still need to be reevaluated. Although phylogenomic data from transcriptomes and plastomes have generated a well-supported backbone for Rosaceae, limited taxon sampling of Prunus in these analyses failed to untangle all relationships (Xiang et al., 2016; Zhang et al., 2017). Moreover, the divergence times of some lineages within Prunus estimated in previous studies need to be further tested. Thus, exploring phylogenetic relationships within Prunus with extensive taxon sampling and molecular characters is essential.

The plastid genome (plastome) is characterized by its conserved structure and typically maternal inheritance, which makes it helpful in reconstructing phylogenetic relationships at both shallow and deep taxonomic levels (e.g., Li et al., 2019; 2021a; Liu et al., 2019, 2020a; 2020b; Li et al., 2020; Thode et al., 2020; Walker et al., 2022; Wang et al., 2020; Bai et al., 2021; Lei et al., 2021; Li et al., 2021b; Wu et al., 2021; Su et al., 2021; Lee et al., 2022; Liu et al., 2022a, 2022b; Xu et al., 2022; Zhang et al., 2022). However, the uniparental inheritance of plastomes limits their power to fully elucidate the evolutionary histories of lineages with pervasive reticulate evolution, such as the family Rosaceae. Therefore, nuclear and plastome data are needed to infer robust phylogenetic reconstruction. Restriction-site associated DNA sequencing (RAD-seq) can readily generate thousands of unlinked nuclear loci that can resolve dense reticulation on shallow systematic levels (e.g., Vargas et al., 2017; Ma et al., 2018; Mu et al., 2020; Zhou et al., 2020; Hodel et al., 2022).

In this study, we used plastomes from genome skimming data and single nucleotide polymorphism (SNPs) data from RAD-seq data to infer the phylogeny of Prunus, and we also reconstructed the evolution of inflorescence types in Prunus. We aim to (1) resolve phylogenetic relationships among the major lineages, (2) detect potential hybridization events, (3) estimate divergence times, and (4) infer the ancestral states and evolutionary transitions among the states of inflorescence types in Prunus.

2. Materials and methods 2.1. Sampling, DNA extraction, library preparation and sequencing

We sampled 38 accessions, including 32 Prunus species and two outgroup species (Table 1), for specific-locus amplified fragment sequencing (SLAF-seq), which is a reduced-representation genomic sequencing approach. Our sampling was taxonomically and morphologically representative, as these 32 Prunus species represented all subgenera and inflorescence types of Prunus. These species are mainly from Asia, which is considered the origin and diversity center of Prunus (Chin et al., 2014). SLAF-seq is a version of RAD-seq that uses a pre-determined reduced representation scheme to optimize and evenly space loci (Sun et al., 2013). For the plastome data, our sampling was identical to the sampling used for nuclear data, with 37 accessions newly sequenced for this study and one (Prunus davidiana (Carrière) Franch., accession number: MH460864) downloaded from GenBank. Total genomic DNA was extracted from the silica-gel-dried leaves using the CTAB method (Doyle and Doyle, 1987).

For RAD-seq data, the Prunus mume (Siebold) Siebold & Zucc. genome (accession number: GCF_000346735.1) was used in an in silico double enzyme digestion to determine the enzyme pair before preparing the library. The library preparation of RAD-seq was completed with the RsaI + HaeⅢ enzyme pair. The paired-end (264–414 bp) sequencing was carried out on the Illumina Hiseq™ 2500 sequencing platform (Illumina, Inc; San Diego, CA, USA) at Beijing Biomarker Technologies Corporation (Beijing, China). For the plastome data, the libraries were prepared in the Molecular Biology Experiment Center, Germplasm Bank of Wild Species in Southwest China using NEBNext® Ultra™ Ⅱ DNA Library Prep Kit. Paried-end sequencing was conducted on a HiSeq 2500™ (Illumina) in BGI (Shenzhen, China).

2.2. SNP calling and plastome assembly

Raw sequencing reads were assembled de novo using ipyrad v.0.9.74 (Eaton, 2014; Eaton and Overcast, 2020) with default parameters, except the clustering threshold set at 0.85 and the maximum alleles per site set at 4. The resulting assembly consisted of 25, 770, 082 nucleotides per species representing 100, 579 loci. This unfiltered 100, 579 locus alignment was used in subsequent coalescent- and concatenation-based phylogenetic inferences. Additionally, the ipyrad assembly was filtered using vcftools (Danacek et al., 2011) to remove loci with minor allele frequency less than 0.05 (--maf = 0.05), and using several filtering levels for missing data to investigate the impact of missing data on the downstream analyses. We used three filtering levels, allowing a maximum of 40%, 60%, or 80% missing data per site (--max-missing = 0.6, 0.4, 0.2). These three SNP matrices were also used in the downstream phylogenetic analyses. We used BLAST searches to filter the RAD-seq data to remove any chloroplast- and mitochondrion-related loci with references from NCBI GenBank (P. mume chloroplast genome (accession number: NC_023798.1) and mitochondrial genome (accession number: NC_060491.1)).

Raw reads generated by genome skimming (Zhang et al., 2015; Liu et al., 2021) were assembled into plastomes using the GetOrganelle pipeline (Jin et al., 2020). The word sizes were set as 125 bp, round numbers were 15, and k-mers were 75, 85, 95 and 105. Complete plastid genomes from the following species were used as seed sequences for the assembly of all sequenced accessions: Prunus kansuensis Rehder (accession number: NC_023956), Prunus maximowiczii Rupr (accession number: NC_026981), P. mume (accession number: NC_023, 798), Prunus padus L. (accession number: NC_026982), Prunus persica (L.) Batsch (accession number: NC_014697), P. pseudocerasus Lindl (accession number: NC_030599). and P. yedoensis Matsum (accession number: NC_026980). All the plastomes were successfully assembled. The seven published sequences were also used as references for annotating the plastomes with PGA (Qu et al., 2019). The annotated sequences were manually revised in Geneious v.11.0.2.

2.3. Phylogenetic analyses

Phylogenetic relationships were inferred for each of these four data sets (full sequence matrix, and those allowing a maximum of 40%, 60%, or 80% missing data per site). Concatenated Maximum Likelihood (ML) phylogenies were inferred using RAxML v.8.2.12 with 100 bootstrap replicates and 20 independent ML searches. The GTRGAMMA model of evolution was used for the full sequence matrix; the ASC_GTRGAMMA model with the Lewis ascertainment correction was used for the three SNP data sets.

Species trees were constructed using the nuclear data in a coalescent framework with SVDQuartets (Chifman and Kubatko, 2014). All quartets were evaluated for each of the four data matrices, and 100 bootstrap replicates were implemented to assess confidence in phylogenetic relationships. The 100, 579 (~250 bp) loci were considered unlinked genes for the full sequence matrix; however, each SNP was considered unlinked in the three SNP data sets.

All the plastome sequences were aligned using MAFFT (Katoh and Standley, 2013) implemented in Geneious v.11.0.2 (Kearse et al., 2012). The plastome data set was used to infer the phylogeny of Prunus based on the maximum likelihood (ML) method. The analyses were performed by RAxML-HPC Black Box 8.2.12 (Stamatakis, 2014) with 1000 bootstrap replicates and the GTR + G model via the CIPRES Science Gateway website (Miller et al., 2010). For the plastome data set, we also reconstructed the phylogeny of Prunus using MrBayes v.3.2 (Ronquist et al., 2012) based on Bayesian inference (BI) methods. Bayesian inference was performed using ten million generations, the first 25% of trees were discarded as burn-in, and trees were sampled every 1000 generations.

2.4. Hybridization tests

We tested the hypothesized allopolyploid origin of the racemose group using the software package Hybrid Detector (HyDe; Blischak et al., 2018). HyDe is an approach similar to ABBA-BABA tests, and uses phylogenetic invariants arising under a coalescent model to identify hybridization among three ingroup taxa polarized by an outgroup. Within each quartet examined, the parameter γ represents the probability of one ingroup species X being sister to species Y, whereas 1-γ represents the probability of species X being sister to species Z. A significant γ value of 0.5 implies that species X is an F1 hybrid resulting from species Y and Z. We ran three separate analyses using this approach. We ran HyDe analyses exhaustively on all combinations of individual species, and then we combined accessions into three groups based on inflorescence type, i.e., solitary, corymbose and racemose. Additionally, we investigated hybridization within each of these three groups. Physocarpus amurensis (Maxim.) Maxim and Prinsepia uniflora Batalin were used as outgroups for all analyses.

2.5. Divergence time estimation

Divergence times were estimated using the plastome matrix of the 38 individuals in BEAST v.2.5.2 (Bouckaert et al., 2019). Based on the published fossil endocarp Prunus wutuensis from Shandong Province in eastern China (Li et al., 2011), a lognormal prior with an offset of 55.0 million years ago (Mya) and a standard deviation of 1 Mya were used to set the stem age of Prunus subg. Cerasus. A normal prior with offset of 35.0 Mya and a standard deviation of 2 Mya were set for the crown age of the racemose group according to the published amber of a staminate flower of Prunus hirsutipetala D.D. Sokoloff, Remizowa et Nuraliev from northwestern Ukraine (Sokoloff et al., 2018). A nodal calibration was set for the crown age of Prunus (Chin et al., 2014), using a normal prior with an offset of 61.5 Mya and a standard deviation of 3.0. The Markov chain Monte Carlo (MCMC) run was conducted for 100, 000, 000 generations. Trees were sampled every 1000 generations. The log files of BEAST analysis were checked with Tracer 1.5 (Rambaut et al., 2018). Results were considered reliable once the effective sampling size (ESS) for all parameters exceeded 200. We used TreeAnnotator to generate the maximum clade credibility tree with 25% of trees discarded as burn-in.

2.6. Ancestral character state reconstruction

We traced the inflorescence evolution of Prunus using the BI tree inferred from the plastome data set. Three inflorescence types were defined (0) simple flowers, (1) corymbs, and (2) racemes. The detailed inflorescence information was obtained from Flora Republicae Popularis Sinicae (http://www.iplant.cn/frps) and observations in the field listed in Table 1. The 'Stochastic Character Mapping' method implemented in Mesquite v.3.51 (Maddison and Maddison, 2018) was employed to reconstruct ancestral character states.

3. Results 3.1. Characteristics of RAD-seq and plastome data sets

RAD-seq yielded 12.9 Gb raw data across the 38 samples included in this study. The total read number ranged from 1, 436, 523 to 5, 977, 937 among all sequenced individuals. The average Q30 percentage was 93.5% and the average GC content was 42.8% (Table S1). The SNP matrices constructed via de novo assembly using ipyrad consisted of a sequence matrix totaling 25, 770, 082 nucleotide sites per species and a SNP matrix of 1, 587, 718 nucleotides per species. Subsequent filtering of the assembly resulted in data matrices of 9, 637, 65, 176, and 448, 083 sites, respectively (Table S2).

The size of Prunus plastomes ranged from 157, 660 bp (P. davidiana) to 159, 000 bp (P. stipulacea Maxim.) in length. Plastomes of all Prunus species we studied had a quadripartite structure, including a large single-copy (LSC, 85, 764–87, 692 bp) region, a small single-copy (SSC, 18, 855–19, 161 bp) region, and two inverted repeat (IR, 26, 198–26, 458 bp) regions (Table S3). The total GC content of all the Prunus plastomes ranged from 36.6% to 36.8%. All Prunus plastomes encoded 113 unique genes, including 79 protein-coding genes (CDS), four ribosomal RNAs (rRNAs), and 30 transfer RNAs (tRNAs). In addition, 17 genes were duplicated in the IRs, of which six, four, and seven were protein-coding genes, rRNAs and tRNAs, respectively (Table S3). The matrix length of the plastome sequence was 169, 232 bp after trimming.

3.2. Phylogenetic relationships of Prunus

Overall, the nuclear and plastome data both showed three strongly supported major clades (Fig. 2 and Figs. S1–S6), each presenting a characteristic inflorescence architecture: the racemose group (clade 1/A), solitary-flower group (clade 2/B) and corymbose group (clade 3/C). Cytonuclear discord was not detected in the backbone of Prunus. Still, there were some conflicts, especially within the major clades, e.g., for section Armeniaca and section Microcerasus M. Roem. of the solitary-flower group, and within the racemose group, as well as within the corymbose group (Fig. 2). The racemose group (clade 1/A) included the subg. Lauro-cerasus, subg. Padus, and Maddenia and Pygeum groups, and its monophyly was strongly supported with a sister relationship to the clade of the remaining two major groups of Prunus (clades 2/B and 3/C). Padus species formed a monophyletic group based on ipyrad assemblies (Figs. 2, S1, S2 and S4). The nuclear phylogenies inferred using the concatenated-based method were broadly congruent with the phylogenies constructed using a coalescent-based approach (i.e., SVDQuartets) (Figs. 2 and S1–S6). However, the species trees and concatenation-based phylogenies showed a different phylogenetic placement of Prunus serotina Ehrh. Whereas P. serotina was always sister to the remaining Padus species in the SVDQuartets trees (Figs. 2, S1 and S2), this was the case in the concatenation-based phylogeny for only one data set—the one with a maximum of 40% missing data (Fig. S4). In the data sets with no limit on missing data (Figs. S3 and S6), P. serotina was sister to a clade containing all species of Padus and Maddenia. In the data set allowing up to 80% missing data, P. serotina was sister to Maddenia (Fig. S5), which matched the chloroplast topology but not the nuclear coalescent topology.

In the rest of this section, the tree generated using the ipyrad assembly allowing a maximum of 60% missing data will be used to discuss the phylogenetic relationships of Prunus because there were only minor differences among assemblies (Figs. 2 and S1–S6). Lauro-cerasus was monophyletic in none of the analyses. One species, Prunus laurocerasus L., was sister to all taxa from Maddenia and Padus in the nuclear tree, in contrast with the sister relationship to the Padus species excluding P. serotine, as shown in the plastome tree. Both nuclear and plastome phylogenies showed that Maddenia and Pygeum each formed a monophyletic group (Fig. 2).

The solitary-flower group (clade 2) contained not only species traditionally assigned to it (subg. Amygdalus, section Armeniaca and subg. Prunus s. str.), but also section Microcerasus species of subg. Cerasus. Amygdalus was the first diverging lineage within the solitary flower group, but there were some conflicts between the nuclear and plastome trees within clade 2. Section Microcerasus was a monophyletic clade in the nuclear tree but not in the plastome tree. Armeniaca and Prunus s. str. formed a sister clade, and together they were sister to Microcerasus in the nuclear tree. Armeniaca (Prunus mandshurica (Maxim.) Koehne and P. mume) formed a clade with Microcerasus, and one Armeniaca species (Prunus sibirica L.) was sister to Prunus s. str. in the plastome tree.

The corymbose group was composed of all subg. Cerasus species (clade 3) and was divided into two subclades, but the position of Prunus cerasoides Buch.-Ham. ex D. Don was not congruent in the plastome and nuclear trees. P. cerasoides and P. serrulata Lindl. formed a clade with P. maximowiczii that was siter to the clade of P. stipulacea and P. discadenia Koehne in the nuclear tree. P. cerasoides and (P. stipulacea + P. discadenia) formed a clade that was sister to the clade of P. serrulata and P. maximowiczii in the plastome tree (Fig. 2).

3.3. Hybridization tests

Of the 13, 485 hypotheses of hybridization tested using HyDe, 4180 contained significant evidence of hybridization (Table S4). There was a wide range of γ values among the significant tests, and over 2/3 of γ values were less than 0.3 or greater than 0.7 (Table 2). Typically, γ values that deviate substantially from 0.5, as is the case here, indicate a more ancient history of hybridization within the group, as opposed to a first-generation hybridization event that a γ value of 0.5 would suggest. Species from the racemose group make up the majority of species included in the significant hybridization comparisons, although racemose species are often found included as either parental species or hybrids (Table 2). The HyDe analyses grouping together accessions into three groups suggested that the corymbose group may have originated via hybridization of the racemose and solitary-flower groups (Table S5). HyDe analyses of hybridization within each of the three inflorescence groups indicated high levels of hybridization (i.e., greater than 30%) (Tables S6–S8).

3.4. Divergence time estimation

The diversification of Prunus began around 67.32 Mya (95% HPD: 55.66–79.77 Mya) in the Late Cretaceous. Clades 2 and 3 split at around 58.38 Mya (95% HPD: 41.23–73.71 Mya), and diversified around 25.92 Mya (95% HPD: 12.83–58.59 Mya) and 11.45 Mya (95% HPD: 3.86–49.58 Mya), respectively (Fig. 3). Clade 1 diversified around 35.41 Mya (95% HPD: 26.51–43.39 Mya). The divergence of Pygeum and Maddenia groups in clade 1 occurred at 12.24 Mya (95% HPD: 2.72–29.29 Mya) and 2.81 Mya (95% HPD: 0.86–12.82 Mya), respectively.

3.5. Ancestral character state reconstruction

Our results indicated that the ancestral inflorescence type of Prunus was the raceme. The ancestral state of clades 2 and 3 in Prunus was a simple flower. As such, the simple flower and corymb states were derived. In the corymbose group, there was an evolutionary reversal to raceme in P. discadenia SN502 (Fig. 4).

4. Discussion

Below, we discuss the main findings of our phylogenetic results and inferences of hybridization events, divergence times, and inflorescence evolution in Prunus.

4.1. Phylogenetic relationships

Our study integrated nuclear and plastome data to infer phylogenetic relationships of Prunus based on species sampling representing all major lineages within the genus. Data sets from both nuclear and plastid genomes supported that Prunus consisted of three main clades—the racemose group, the solitary-flower group, and the corymbose group. The corymbose group was sister to the solitary-flower group, and then together sister to the racemose group with strong support. This topology was consistent with previous findings based on plastid genes (Chin et al., 2010, 2014; Zhao et al., 2016) and those from some single-copy nuclear genes and plastomes (Hodel et al., 2021), but in conflict with results from other analyses of nuclear genes (Chin et al., 2010, 2014). Although our two data sets resulted in a consistently robust backbone, strong phylogenetic conflict was recovered at some shallow-level nodes within each of the three main groups.

The racemose group was previously considered paraphyletic according to phylogenetic trees using a few nuclear gene loci (Chin et al., 2014; Zhao et al., 2016), but it was monophyletic in both nuclear and plastome trees (Fig. 2 and Hodel et al., 2021). The observed monophyly of the racemose group may have resulted from the majority of nuclear loci included in this study tracking the maternal phylogeny of the racemose group. Often, a few outlier genes with outsized influence may skew phylogenetic relationships, so even when relationships appear strongly supported, underlying cytonuclear or gene tree discord may give false confidence in incorrect relationships (Walker et al., 2022). The uncertain phylogenetic placement of the racemose species P. serotina, with one concatenation phylogeny matching the plastid topology as opposed to the nuclear coalescent topology, may indicate that the maternal phylogeny readily dominated some of the SNP data sets when constructing concatenated trees (Fig. S5). In any case, the origin of the racemose group needs to be better explored with single-copy nuclear data that track multiple homologous loci of a specific gene (Zimmer and Wen, 2015) and/or whole nuclear genome data. Both Pygeum and Maddenia groups were strongly supported as monophyletic, which is consistent with previous results (Zhao et al., 2018; Su et al., 2021). The paraphyly of subg. Padus and subg. Lauro-cerasus was suggested in previous studies (Bortiri et al., 2002; Wen et al., 2008; Chin et al., 2010, 2014; Liu et al., 2013; Shi et al., 2013; Zhao et al., 2016, 2018). Our results based on plastome data further confirm these conclusions. Padus and Lauro-cerasus may have each derived from multiple hybridization and allopolyploidy events (Zhao et al., 2016). However, the three trees based on ipyrad de novo assemblies show a monophyletic Padus, and the difference in the position of P. serotina was the most obvious incongruence among the topologies obtained with de novo assembly and those published in previous studies (Liu et al., 2013; Chin et al., 2014; Zhao et al., 2018) in nuclear trees.

Prunus maackii Rupr., a species mainly distributed in northeast China, Korea, and eastern Russia, has been treated as a member of subg. Padus (Rehder, 1956; Yü et al., 1986; Lu et al., 2003). This species has racemose and deciduous inflorescences similar to some species of Padus (Rehder, 1956; Wen et al., 2008). Yet it resembles some members of Cerasus with highly incised stipules (Liu et al., 2013) and a short inflorescence. P. maackii has been supported to be nested within subg. Cerasus based on nuclear and plastid sequences (Lee and Wen, 2001; Bortiri et al., 2006; Wen et al., 2008; Shi et al., 2013; Liu et al., 2013; Chin et al., 2014; Zhao et al., 2016). This species was reported to hybridize with P. maximowiczii of the Cerasus group and formed natural hybrids (Wen et al., 2008; Liu et al., 2013; Shi et al., 2013). Thus, some researchers suggested P. maackii be placed in subg. Cerasus (Wen et al., 2008; Liu et al., 2013). Yet, we found that it formed a clade with Padus species in both nuclear and plastome trees with high support for the first time (Fig. 2), which challenges previous studies. Nevertheless, P. maackii is similar to most members of subg. Padus, and they share two morphological characteristics, i.e., terminal racemes with a few leaves at the base of their peduncle and leaves without glands on the petiole. Our results hence support the placement of P. maackii in subg. Padus.

The solitary-flower group showed several conflicts between plastome and nuclear trees (Fig. 2). Section Armeniaca was monophyletic in the nuclear tree, but not in the plastome tree. Due to maternal inheritance of plastomes and biparental inheritance of nuclear data, section Armeniaca may be derived from one paternal parent and several maternal lineages.

Even though the corymbose group refers to all species of subg. Cerasus, species of Microcerasus were scattered in the solitary-flower group in our analyses. Cerasus, excluding section Microcerasus, was monophyletic and had some cytonuclear discord among Cerasus species (Fig. 2). Microcerasus was regarded as a section of subgenus Cerasus by Rehder (1956) and as a subgenus of the more narrowly defined genus Cerasus by Yü et al. (1986). Microcerasus comprises shrubby and woody species with three axillary buds and a short pedicel, similar to the solitary-flower group members (Lee and Wen, 2001; Wen et al., 2008; Shi et al., 2013). However, true cherries share only one axillary bud at each leaf axil. Moreover, it was reported that Microcerasus species more readily hybridized with species of the solitary-flower group than with those of the Cerasus group (Kataoka et al., 1988).

Two decades ago, Lee and Wen (2001) were the first to systematically investigate the phylogenetic relationships of Prunus based on ITS sequences, and showed that two Microcerasus species had been closely related to the solitary-flower group rather than to the corymbose group. This result was supported by those inferred from different molecular markers in many subsequent studies (Bortiri et al., 2001, 2002, 2006; Wen et al., 2008; Liu et al., 2013; Shi et al., 2013; Chin et al., 2014; Zhao et al., 2016). Our nuclear and plastid results revealed that the Microcerasus species (Prunus tomentosa Thunb. and P. japonica var. nakaii (H. Lév.) Rehder) were nested within the solitary-flower group, which is congruent with previous studies. Our analysis did not support the inclusion of Microcerasus into Cerasus. We found that the ovule developmental characteristics of P. tomentosa were more similar to members of the solitary-flower group than to Cerasus (personal observation by L. Zhao). Shi et al. (2013) also suggested that the Cerasus group only included true cherries, and that Microcerasus did not belong to it. Based on morphological and molecular evidence, it seemed reasonable to assign Microcerasus to the solitary-flower group. Mowrey and Werner (1990) pointed out that due to its paraphyly in all analyses Microcerasus was not a clade, which has been supported by subsequent studies (Lee and Wen, 2001; Bortiri et al., 2002; Wen et al., 2008; Chin et al., 2014). Our results also show that Microcerasus species did not form a clade in the plastome tree. Therefore, further studies should sample additional Microcerasus species and closely related taxa to better understand their origin(s) and systematic placement.

4.2. On the hybrid origins of major Prunus lineages

Zhao et al. (2016) hypothesized that the racemose group derived from multiple hybridization events and allopolyploid speciation events. The maternal parent may have been extinct, and the paternal parents shared common ancestry with various members of the corymbose-solitary flower lineages, leading to the observed conflicts between the plastid and nuclear topologies (Zhao et al., 2016). However, our results support the monophyly of the racemose group (Fig. 2).

Although these results challenge the hypothesis of multiple hybrid origins of the racemose group, the HyDe analysis based on RAD-seq data detected frequent hybridization events in this group. In concert with previous studies, the HyDe results indicate pervasive hybridization and/or allopolyploidy within Prunus, especially within the racemose group (Table 2). Our goal was to test the hypothesized allopolyploid origin of the racemose group, but it is challenging to use hybrid detection methods based on extant taxa to assess ancient hybridization that may have involved extinct lineages. Furthermore, subsequent introgression and/or repeated hybridization can lead to conflicting genomic signals in different species. Nonetheless, a large number of hybridization events detected by HyDe highlights the frequency of hybridization within this genus. Analyses of individual species and within groups indicate that hybridization was pervasive throughout Prunus, and our data do not readily confirm any specific hypotheses of ancient hybridization (Tables S5–S8). Because many species in the racemose group are polyploid, HyDe may pick up signals of genome doubling (i.e., allopolyploidy) instead of homoploid hybridization.

Previous studies with less representation of racemose species (Hodel et al., 2021) found that approximately 1% of hybridization tests were significant, in contrast with the present study in which 31% of the tests were significant. A similar approach in another Rosaceae lineage with frequent hybridization (the apple genus Malus) found that 24% of hybridization tests were significant (Liu et al., 2022a). A study of the Hawaiian genus Myrsine using similar markers as the present study (i.e., RAD-seq) found similarly high levels of hybridization using HyDe (between 27.5% and 30.4% of tests significant, depending on the data set; Appelhans et al., 2020). RAD-seq data have been used successfully for non-polyploid taxa with fewer reticulation events (e.g., Hodel et al., 2022). However, RAD-seq loci cannot be used for accurately estimating ancient hybridization events between the polyploid species and its close diploid relatives (e.g., Wang et al., 2021), especially in the presence of extinction events, genome doubling and/or subsequent recent hybridizations.

Our HyDe analyses with accessions grouped based on inflorescence type showed a potential signal of ancient hybridization of the corymbose group, i.e., between the racemose group and solitary-flower group (Table S5), in contrast to the hypothesis of allopolyploid origins of the racemose group suggested by previous studies (Chin et al., 2014; Zhao et al., 2016). While homoploid hybrid origin of the corymbose group with a maternal 'Solitary flower' lineage and a diploid paternal ancestral 'Racemose' lineage is plausible (Chin et al., 2014), it is not supported by the higher ploidy levels of all extant members of the racemose group, nor is it consistent with patterns of cytonuclear discord concerning the relationships among the three major groups observed in previous studies (Chin et al., 2014; Zhao et al., 2016) or the lack of cytonuclear discord concerning those relationships observed here. Moreover, our HyDe analyses using extant taxa grouped based on inflorescence types could not explicitly test Zhao et al.'s (2016) hypothesis of multiple allopolyploid origins of the racemose group, because that hypothesis invoked members of a now-extinct lineage as the maternal parents in the hybridization events. In fact, if an ancestral member of the corymbose group acted as one or more of paternal parents in the hybridization events, as suggested by Zhao et al. (2016), then one would expect exactly the results seen in our species group HyDe tests here, because some nuclear genes would group the corymbose species with their sister solitary-flower group and others would group them with the hybrid racemose lineages to which they gave rise via hybridization. On the other hand, the corymbose inflorescence seems to represent an intermediate morphological state between the raceme and the simple flowers, which may be consistent with the above hybridization hypothesis of the corymbose group. Future additional studies are needed to rigorously test these competing hypotheses using phylogenomic and developmental morphological data.

4.3. Temporal diversification of main clades in Prunus

A previous study based on four plastid genes and ITS sequences estimated that the ancestor of Prunus emerged between 56.7 and 67.4 Mya (Chin et al., 2014). The age of Prunus in the present study slightly predates their estimate (Fig. 3). In our chronogram, the ancestor of clades 2 and 3 firstly split around 58.38 Mya (95% HPD: 41.23–73.71 Mya). This is congruent with the fossil record of P. wutuensis, which mostly resembled extant P. yedoensis of subgenus Cerasus (Li et al., 2011). The evolution of plants has often been impacted by paleoclimatic and geologic events. During the boundaries between epochs of the Paleogene, key paleoclimatic events caused multiple effects on the evolution and distribution of plants. Around 55.8 Mya, there was an abrupt period of global warming caused by a transient burst of CO₂, known as the Paleocene-Eocene Thermal Maximum (PETM) (Currano et al., 2008; McInerney and Wing, 2011). During this period, floristic changes have been shown to occur in response to climate fluctuations (Wing and Currano, 2013). Clades 2 and 3 split around 58.38 Mya, which corresponds to the PETM. Therefore, our results seem to indicate that climatic change may have influenced the diversification of Prunus during the PETM.

In addition, the Eocene–Oligocene transition (EOT, 30–40 Mya) in the Cenozoic witnessed a global cooling and led to a reorganization of organisms (Prothero, 1994; Sun et al., 2014; Deng et al., 2020). The diversification of clades 1 and 2 was herein estimated to be around 35.41 Mya (95% HPD: 26.51–43.39 Mya) in the Late Eocene and 25.92 Mya (95% HPD: 12.83–58.59 Mya) in the Late Oligocene, respectively (Fig. 3). There was another climatic cooling event after the Middle Miocene Climatic Optimum (MMCO, around 15 Mya) (Flower and Kennett, 1994). The molecular dating analysis indicate that clade 3 of Prunus initially diverged when a climatic cooling event occurred after MMCO. The divergence time of the three clades coincided with the two cooling events, which also indicates that the paleoclimatic events might have impacted Prunus diversification and evolution.

4.4. Ancestral character state reconstruction of inflorescences

By tracking the evolution of inflorescence types on the Prunus phylogeny, we inferred that the raceme was the ancestral state of inflorescence in this genus. The solitary flower and corymbose inflorescence types might be derived from the reduction of the flower number and suppression of the rachis, respectively. This result is consistent with the hypothesis that racemes are a primitive state in angiosperms (Stebbins, 1973). Although terminal solitary flowers (Parkin, 1914) and panicles (Takhtajan, 1991) have also been regarded as an ancestral state of inflorescences, there has been some skepticism about the reliability of any of these ideas because the earliest-diverging angiosperms, the ANITA grade, do not possess any panicles, and single flowers terminal on shoots are only present in several genera, such as Austrobaileya C.T. White (in part), Schisandra Michx (in part), and Illicium L. (Endress, 2010).

Simple flowers were the ancestral state of clades 2 and 3 of Prunus based on mature inflorescence types. In fact, a lateral floral primordium was initiated but it did not develop further, contributing to simple flowers at maturity in P. persica (personal observation by L. Zhao). P. salicina and P. americana Marshall have a typical inflorescence type with 2–3 flowers clustered in a bud (personal observation by L. Zhao and J. Wen), which resemble a corymb more than a simple flower. Such developmental/morphological evidence suggests that there may be a transition between simple flowers and corymbs. Hence, we assume the ancestral state of clades 2 and 3 had a transitional inflorescence type, which shows simple flowers only at maturity, with more than one flower observed even then. Corymb and simple flowers in Prunus were formed by continuing and stopping development of lateral floral primordium from their ancestor, respectively. Furthermore, P. discadenia SN502, a member of subg. Cerasus, had racemose inflorescence. It seems that this represents a case of reverse evolution in Prunus, i.e., racemes likely evolved from corymbs by elongation of the rachis. However, we only analyzed the mature character state of Prunus and did not provide robust evidence for developmental evolutionary histories of inflorescences in this study. Future studies integrating inflorescence development and phylogenetic analyses are needed to understand the evolutionary directionality of Prunus inflorescences better.

5. Conclusions

This study sheds light on the phylogenetic relationships, hybridization events and inflorescence evolution of Prunus, an economically important plant group. RAD-seq data have been successfully used to elucidate the phylogenetic relationships of various polyploid lineages, such as bamboos (Guo et al., 2021) and Salix L. (Wagner et al., 2020). However, RAD-seq data cannot distinguish orthologs from paralogs for the polyploid species, limiting their utility for tracing ancient hybridization events. Thus, additional studies using whole-genome sequencing or genome resequencing including more extensive sampling of the racemose group and their allies are needed to test the competing hypotheses of multiple allopolyploid origins of the racemose group vs. the homoploid hybrid origin of the corymbose group (Zhao et al., 2016).

Acknowledgments

We thank Prof. Zhong-Hu Li and Dr. Miao–Miao Ju of Northwest University and Dr. Li Feng of Xi'an Jiaotong University, Dr. Fu-Zhen Guo and Xiao-Hua He from Instrument Sharing Platform of Northwest A&F University and Molecular Biology Experiment Center, Germplasm Bank of Wild Species in Southwest China for help with data analysis, and Mr Hai-Ning Li, Xiao-Jie Li and Prof. You Zhou for providing plant photographs (Fig. 1L–O). We also thank Dr. Florian Jabbour and Dr. Valéry Malécot for their valuable comments. Computational analyses were partially conducted on the Smithsonian Institution High Performance Computing Cluster (SI/HPC, "Hydra": https://doi.org/10.25572/SIHPC). This work was supported by grants from National Natural Science Foundation of China (No. 32170381 and 31770200).

Authors contributions

Na Su: Formal analysis, Investigation, Methodology, Software, Visualization, Writing - original draft, Writing - review & editing. Richard G. J. Hodel: Formal analysis, Investigation, Methodology, Software, Visualization, Writing - original draft, Writing - review & editing. Xi Wang: Formal analysis, Investigation, Methodology, Software, Visualization, Writing- original draft. Jun-ru Wang: Formal analysis, Investigation, Methodology, Software, Visualization, Writing - original draft. Si-yu Xie: Methodology, Software, Visualization, Writing - original draft. Chao-xia Gui: Writing - original draft. Ling Zhang: Writing - original draft. Zhao-yang Chang: Writing - original draft. Liang Zhao: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing - review & editing. Daniel Potter: Conceptualization, Investigation, Methodology, Project administration, Supervision, Validation, Writing - review & editing. Jun Wen: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing - review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2023.03.013.