Eriobotrya Belongs to Rhaphiolepis (Maleae, Rosaceae): Evidence From Chloroplast Genome and Nuclear Ribosomal DNA Data

The Eriobotrya-Rhaphiolepis (ER) clade consists of about 46 species distributed in East and Southeast Asia. Although Eriobotrya and Rhaphiolepis have been supported to form a clade, the monophyly of Eriobotrya and Rhaphiolepis at the genus level has never been well tested and their phylogenetic positions in Maleae still remain uncertain. This study aims to reconstruct a robust phylogeny of the ER clade in the framework of Maleae with a broad taxon sampling and clarify the phylogenetic relationship between Eriobotrya and Rhaphiolepis. This study employed sequences of the whole plastome (WP) and entire nuclear ribosomal DNA (nrDNA) repeats assembled from the genome skimming approach and included 83 samples representing 76 species in 32 genera of Rosaceae, especially Maleae. The Maximum Likelihood (ML) and Bayesian Analysis (BI) based on three datasets, i.e., WP, coding sequences of plastome (CDS), and nrDNA, strongly supported the paraphyly of Eriobotrya, within which Rhaphiolepis was nested. Our plastid tree supported the sister relationship between the ER clade and Heteromeles, and the nrDNA tree, however, did not resolve the phylogenetic placement of the ER clade in Maleae. Strong incongruence between the plastid and the nuclear trees is most likely explained by hybridization events, which may have played an important role in the evolutionary history of the ER clade. Molecular, morphological, and geographic evidence all supports the merge of Eriobotrya with Rhaphiolepis, which has the nomenclatural priority. We herein transferred 36 taxa of Eriobotrya to Rhaphiolepis. We also proposed a new name, Rhaphiolepis loquata B.B.Liu & J.Wen, for the economically important loquat, as the specific epithet “japonica” was pre-occupied in Rhaphiolepis.

The generic delimitation in Maleae has been notoriously difficult, which may be due to the low sequence divergence resulted from ancient, rapid radiations (Wolfe and Wehr, 1988;Campbell et al., 2007) and/or frequent hybridizations Lo and Donoghue, 2012;Liu et al., 2019). As the one of the few genera of Maleae largely distributed in subtropical and tropical regions , Eriobotrya consists of ca. 15-20 species ranging from the Himalayas throughout continental southeast Asia to Japan and the islands of western Malesia (Kalkman, 2004). The well-known fruit loquat, Eriobotrya japonica, has been widely cultivated all over the world, making Eriobotrya of great economical importance. Characterized by persistent sepals and craspedodromous lateral veins of leaves, Eriobotrya has been treated to be distinct from Rhaphiolepis in all taxonomic literature and floras (Vidal, 1968;Vidal, 1970;Robertson et al., 1991;Kalkman, 1993;Lu et al., 2003;Kalkman, 2004). The monophyly of the Eriobotrya-Rhaphiolepis complex has been strongly supported by recent molecular phylogenetic studies based on limited taxon sampling, however, its phylogenetic placement was still controversial (Campbell et al., 2007;Li et al., 2009;Li et al., 2012;Lo and Donoghue, 2012;Liu et al., 2019;Xiang et al., 2017;Yang et al., 2012;Zhang et al., 2017). In addition, with the synapomorphies of the proportionally large seed with rounded or wide-elliptic cross-section and the lack of endosperm (Aldasoro et al., 2005;Rohrer et al., 1991), Eriobotrya and Rhaphiolepis were morphologically distinct from other members of Maleae. Although the nuclear Adh topology indicated the nonmonophyly of Eriobotrya with Rhaphiolepis nested within it, Yang et al. (2012) did not provide any discussions on the taxonomic implications.
As a rapid and cost-effective method in phylogenomics, genome skimming (i.e. low-coverage genome shotgun sequencing) could generate a large number of phylogenetically informative sites from the whole chloroplast genome, partial mitochondrial genome, and entire nuclear ribosomal DNA (nrDNA) repeats (Straub et al., 2012;Zhang et al, 2015;Zhang et al, 2017;Zimmer and Wen, 2015). Strongly supported incongruence between the maternally inherited plastomes and nuclear markers has been utilized to address polyploidy and reticulate evolution (Bock et al., 2014;McKain et al., 2018;Liu et al., 2019). Recently Liu et al. (2019) have successfully clarified the phylogenetic relationships between Photinia and its morphologically similar allies in the framework of Maleae using whole chloroplast genomes and nuclear ribosomal DNA (nrDNA) repeats assembled from genome skimming approach. The case study showed great potential in resolving the phylogenetic relationships in Maleae using genome skimming data. These datasets have been successfully applied to reconstruct the phylogenetic history of various plant lineages (e.g., in Bock et al., 2014;Zhang et al., 2015;Dillenberger et al., 2018;Mabry and Simpson, 2018;Thomson et al, 2018;Wen et al., 2018;Liu et al., 2019).
The main objectives of the present study were to clarify the phylogenetic placements of Eriobotrya and Rhaphiolepis with a broader sampling using sequences of the whole chloroplast genomes and entire nrDNA repeats via genome skimming. We aim to (1) explore the phylogenetic position of the Eriobotrya-Rhaphiolepis complex; (2) test the monophyly of Eriobotrya and clarify its phylogenetic relationship with Rhaphiolepis;and (3) discuss the taxonomic implications of the new phylogenetic results.

Sampling, DNA Extraction, and Sequencing
We used 83 chloroplast genomes and 68 nrDNA repeats for this analysis, in which 37 samples were sequenced for this study, 31 samples were from our previous study (Liu et al., 2019), and 15 accessions of plastomes were from GenBank. These 83 samples of plastomes represented 76 species and 32 genera which represented nearly all the genera recognized currently in Maleae expect Chamaemeles Lindl. As the core ingroup, nine species of Rhaphiolepis and 12 species of Eriobotrya were included in the phylogenetic analysis, representing all the morphological and geographic diversity ( Figure 1). In order to resolve the phylogenetic positions of the Eriobotrya-Rhaphiolepis complex, we extensively sampled its closely related genera in the plastid tree, i.e., Cotoneaster Medik., Heteromeles M.Roem., and Photinia, as well as its closely related genus Stranvaesia Lindl. in the nrDNA tree (Liu et al., 2019). We also sampled 67 accessions of plastomes from 30 genera of Maleae as the ingroup, i.e., Amelanchier, Aronia, Chaenomeles Lindl., Cormus Spach, The total genomic DNAs sequenced for this study were extracted either from silica-gel dried leaves or from the specimens in herbarium PE and US (Supplementary Table 1). The DNAs of the samples with species names marked with a (Supplementary Table 1) were extracted with a modified CTAB (mCTAB) method (Li et al., 2013b) in the lab of the State Key Laboratory of Systematic and Evolutionary Botany, CAS, China. Library preparations were done at Novogene, Beijing, China using NEBNext ® Ultra ™ II DNA Library Prep Kit and the libraries were sequenced on the NovaSeq 6000 Sequencing System. The DNAs of the samples with b (Supplementary Table 1) were extracted using Qiagen DNeasy ® plant mini-kit (Qiagen Gmbh, Hilden, Germany) following the manufacturer's protocol and the libraries were prepared with the same kit in the Laboratories of Analytical Biology (LAB), National Museum of Natural History (NMNH), Smithsonian Institution (SI), USA. The pooled libraries were shipped to Novogene in UC Davis Sequencing Centre, Davis, California, USA on the Illumina HiSeq 2500 instrument. Paired-end reads of 2 × 150 bp were generated on both sequencing platforms.

Chloroplast Genome and nrDNA Assembly, Annotation, Gene Map, and Alignment
We removed the adaptors introduced by Illumina sequencing using cutadapt 2.4 (Martin, 2011) with AGATCGGAAGAGC as the forward and the reverse adaptor. The results were checked for quality control with FastQC 0.11.8 (Andrews, 2018). NOVOPlasty 3.6 (Dierckxsens et al., 2016) was then used to assemble the plastomes de novo with the ribulose-1,5bisphosphate carboxylase/oxygenase large subunit (rbcL) gene sequence as the seed. Twenty-two circular plastomes (59.5%) of the 37 samples were assembled with NOVOPlasty. The remaining 15 samples were assembled using Zhang et al. (2015)'s successive method combining mapping-based and de novo assembly. For these 15 samples, the raw data generated from the Illumina HiSeq runs were trimmed using Trimmomatic v.0.39 (Bolger et al., 2014), removing bases below PHRED 15 within a sliding window of four bases, keeping only reads of 36 bases or longer. The results were also checked by FastQC for quality control. Zhang's successive assembly used a multistep approach, including reference-guided assembly and de novo assembly to obtain a high-quality plastome. The referenceguided assemblies were performed by Bowtie2 2.3.5.1 (Langmead and Salzberg, 2012) using only the forward and reverse paired reads generated by Trimmomatic v.0.39 (Bolger et al., 2014). We used the closely related plastomes generated from the previous step of this study or the GenBank as the reference genome. After obtaining a complete plastome for one species, we used it as the reference to assemble its closely related samples (see Zhang et al., 2015). The resulting plastomes were used as subsequent references. De novo assemblies were constructed by SPAdes 3.13.1 with careful error correction and  K-mer length of 21, 33, 55, 77 (Nurk et al., 2013). To correct errors and ambiguities resulting from each approach, the scaffolds obtained by de novo assembly, in combination with the contigs generated from NOVOPlasty were mapped to the plastome from the reference-guided assembly. Through the combined effort using NOVOPlasty and Zhang et al. (2015)'s method, we assembled high-quality plastomes for all samples. We also used Zhang et al. (2015)'s method to assemble the entire nrDNA repeats. The raw data were trimmed by Trimmomatic v.0.39 (Bolger et al., 2014). We used the previously assembled nrDNA sequences (Eriobotrya cavaleriei (H.Lév.) Rehder [GenBank accession no. MN215982]and E. deflexa (Hemsl.) Nakai [GenBank accession no. MN215978]) to capture and map the reads to reference by Bowtie2 2.3.5.1 (Langmead and Salzberg, 2012). The scaffolds generated by SPAdes 3.13.1 (Nurk et al., 2013) in assembling the plastome were mapped to the consensus sequences, and the gaps and ambiguous sites generated by the reference method were corrected. If the quality of nrDNA sequences from de novo assembly was not good, we used a different reference to assemble the sequences again with the above procedure. Generally, we obtained high-quality nrDNA sequences using this successive approach that combined the strengths of reference-guided and de novo approaches.
The assembled plastid genomes and nrDNA repeats were annotated using Geneious Prime (Kearse et al., 2012) with a closely related and well-annotated sequence downloaded from NCBI as a reference (KT633951), and the results of automated annotation were checked manually. The coding sequences of plastome were translated into proteins for checking the start and stop codons manually in Geneious Prime. With the two reverse complementary repeats in the plastome of Rosaceae, the boundary of the large-single copy (LSC), small-single copy (SSC), and inverted-repeats (IRs) for each plastome were verified by the Find Repeats in Geneious Prime. Eighty-three plastomes were aligned by MAFFT v.7.409 (Nakamura et al., 2018) with default parameters. To reduce the systematic errors produced by poor alignment, we used the trimAL v.1.2 (Capella-Gutiérrez et al., 2009) with a heuristic method to decide on the best-automated method to trim the alignment of the whole plastome (WP). All 79 coding sequences (CDSs) of each plastome were extracted separately by Geneious Prime (Kearse et al., 2012), each coding sequence was aligned by MAFFT as specified above, and the alignment of each gene was concatenated by AMAS (Borowiec, 2016). Sixty-eighty nrDNA sequences were also aligned by MAFFT as specified above, and each region of nrDNA was extracted separately by Geneious Prime (Kearse et al., 2012), and then concatenated by AMAS (Borowiec, 2016). The gene map of Eriobotrya japonica and Rhaphiolepis indica (L.) Lindl. ex Ker Gawl. were generated by OrganellarGenomeDRAW (OGDRAW) version 1.3.1 (Greiner et al., 2019).

Phylogenetic Analyses
We first performed independent phylogenetic analyses for each dataset (WP and nrDNA) obtained via genome skimming using the Maximum Likelihood (ML) and Bayesian Inference (BI). As the sequence of two IR copies was completely or nearly identical, only one copy of inverted repeat (IR) region was included for the whole plastome (WP) analyses. We used all the 79 coding sequences (CDSs) of each plastome extracted separately by Geneious Prime for phylogenetic analyses. For nrDNA sequences, the intergenic spacer (IGS) and external transcribed spacer (ETS) regions were difficult to align, only part of ETS and the complete sequences of 18S, ITS1, 5.8S, ITS2, and 26S in nrDNA were used for the phylogenetic analyses.
The best-fit partitioning schemes and nucleotide substitution models for the coding sequences of whole plastome and nrDNA dataset were estimated using PartitionFinder2 (Stamatakis, 2006;Lanfear et al., 2016) for unpartitioned whole plastome or partitioned coding sequences (CDS) of plastomes and nrDNA sequences. Under the corrected Akaike information criterion (AICc) and linked branch lengths, the PartitionFinder2 were performed by greedy (Lanfear et al., 2012) and rcluster (Lanfear et al., 2014) algorithm option for WP, plastid CDS, and nrDNA datasets, respectively, with prior defined data blocks by codon positions of each protein-coding genes and all models. The partitioning schemes and evolutionary model for each subset were used for the downstream ML and BI analyses. The ML tree was inferred by IQ-TREE v.1.6.9 (Nguyen et al., 2015) with 1000 bootstrap replicates using UFBoot2 (Hoang et al., 2017) and collapsing near zero branches option. The BI was performed with MrBayes 3.2.7 (Ronquist et al., 2012). The Markov Chain Monte Carlo (MCMC) analyses were run for 10,000,000 generations. The stationarity was regarded to be reached when the average standard deviation of split frequencies remained below 0.01. Trees were sampled at every 1,000 generations, and the first 25% of samples were discarded as burn-in. The remaining trees were used to build a 50% majority-rule consensus tree. The ML and BI trees were visualized using Geneious Prime (Kearse et al., 2012).

RESULTS
The chloroplast genome of Rhaphiolepis indica (the type species of Rhaphiolepis) was 159,466 bp in length with a classic quadripartite structure that comprised of inverted repeat's pairs (Figure 2), and that of Eriobotrya japonica (type species of Eriobotrya) had the same structure as R. indica with a length of 159,156 bp (Supplementary Figure 1). They contained the same number of coding sequences (79), tRNAs (37), and rRNAs (8). No significant rearrangements or gene losses were found in the other species of Eriobotrya and Rhaphiolepis.
The aligned matrix of the 83 whole chloroplast genomes was 129,168 bp in length with the pool alignment trimmed by trimAL (Capella-Gutiérrez et al., 2009). The best-fit model of nucleotide substitutions for the ML analysis was TVM+I+G calculated by PartitionFinder, and that for the BI analysis was GTR+G+I. The 79 concatenated CDS sequences from 83 plastomes were generated an aligned matrix of 67,961 bp in length, which was split into 25 sets of sites (aka data blocks) with the same best scheme for each data block for the subsequent ML and BI analysis. The ML and BI analysis of WP dataset resulted in five strongly supported clades within Maleae ( Figure 3): clade A (Kageneckia and Lindleya), clade B (only Vauquelinia), clade C (only Pyracantha), clade D (Amelanchier, Crataegus, Hesperomeles, Malacomeles, and Peraphyllum), and two large clades (E and F). These five major clades were also recovered by the phylogenetic analysis of CDS (Supplementary Figure 2). Two tropical American genera (Kageneckia and Lindleya) formed clade A, which constituted the first diverged major clade of Maleae. Clade A was sister to the North American genus Vauquelinia (clade B), and then together they were sister to the core Maleae ( Figure 3 and Supplementary Figure 2). As the basalmost group of the core Maleae, Pyracantha (clade C) was sister to a large clade, including clades D, E, and F. Most of the members of clade D were from the New World, and those of clades E and F were from Eurasia except that Aronia and Heteromeles are from North America. The two major Eurasian clades (E and F) were sister to each other, and then together sister to the New World clade D (Figure 3 and Supplementary  Figure 2).
The concatenated nrDNA data (six regions, i.e. partial ETS, 18S, ITS1, 5.8S, ITS2, and 26S) generated an aligned matrix of  6,361 bp in length. Each region of nrDNA was treated as a data block, and these six sets of sites were partitioned into five subsets with each corresponding models of molecular evolution for the subsequent ML and BI analyses. The ML and BI analyses had nearly the same topology ( Figure 4). The ER clade was wellsupported in the nrDNA tree (BS = 100, PP = 1), however, its phylogenetic position was not resolved (Figure 4). Two species of Eriobotrya (E. henryi and E. seguinii) were strongly supported to be sister to all Rhaphiolepis species sampled here (BS = 90, PP = 0.96), and then together they were sister to a clade of the other species of Eriobotrya (E. cavaleriei, E. deflexa, E. japonica, E. malipoensis, E. obovata, and E. salwinensis) (BS = 100, PP = 1).
The topology of the ER clade generated from the nrDNA dataset showed strong conflicts with that from the plastome dataset ( Figure 5).

Enigmatic Phylogenetic Position of the Eriobotrya-Rhaphiolepis Clade
Our phylogenetic results based on the whole chloroplast genome and entire nrDNA repeats strongly supported the monophyly of the Eriobotrya-Rhaphiolepis complex (Figure 3 and Supplementary Figure 2). The ER clade was first reported by Campbell et al. (2007) based on phylogenetic analyses using six cpDNA regions, and this relationship was supported by subsequent studies using chloroplast regions or even plastomes, but with limited taxon sampling (Lo and Donoghue, 2012;Zhang et al.,2017;Sun et al., 2018;Liu et al., 2019). In addition, Eriobotrya and Rhaphiolepis have also been supported to form a clade based on nuclear data, e.g., GBSSI-1A and nrITS plus GBSSI-2B sequences (Campbell et al., 2007), as well as nuclear Adh sequences (Yang et al., 2012). Although the topology of [Eriobotrya, (Rhaphiolepis, Vauquelinia)] was reported by Campbell et al. (1995) using nrITS and a small portion of the 5.8S sequences, all the subsequent studies based on the nrITS sequences supported the close relationship between Eriobotrya and Rhaphiolepis (Campbell et al., 2007;Li et al., 2009;Li et al., 2012;Lo and Donoghue, 2012) instead of between Rhaphiolepis and Vauquelinia (Campbell et al., 1995). Our nrDNA tree showed a distant relationship between Vauquelinia and the ER clade. In addition, Vauquelinia can be easily distinguished from the members of the ER clade morphologically by the former's dry capsular fruits, winged seeds with endosperms, and different wood ray anatomy (Zhang, 1992). Morphologically the ER clade is supported by two synapomorphies: the rounded or widely elliptic cross-section of seeds and the absence of endosperm (Aldasoro et al., 2005). The ML and BI analyses based on the WP and CDS dataset strongly supported the sister relationship between the ER clade and Heteromeles, then together they were sister to Photinia.   Liu et al., 2019). Such uncertainties may be due to the limited markers and/ or samples used in the previous studies, or extinctions of the closest relatives of the ER clade. The ER clade was found to be closely related to the small North American genus Heteromeles based on our plastome data ( Figure 3 and Supplementary Figure 2). But we will test this relationship in our further analyses of Maleae. By contrast, the ML and BI analyses based on the nrDNA dataset did not resolve the phylogenetic position of the ER clade. The ML result supported that the ER clade was sister to a large clade that includes Chaenomeles, Cydonia, Dichotomanthes, Docynia, Heteromeles, Malus, Phippsiomeles, Photinia, Pourthiaea, Pseudocydonia, Pyracantha, Pyrus, and Stranvaesia ( Figure 4). The BI result did not resolve the relationships among the ER clade, the Central American Phippsiomeles, and the remaining genera (Chaenomeles, Cydonia, Dichotomanthes, Docynia, Heteromeles, Malus, Photinia, Pourthiaea, Pseudocydonia, Pyracantha, Pyrus, and Stranvaesia). Campbell et al. (1995) recovered the ER clade as well as Vauquelinia as forming the basalmost clade of Maleae using ITS sequences, and this result was supported by the cladistic analysis using 16 morphological and anatomical characters (Aldasoro et al., 2005). While the ER clade was placed either as sister to Pyrus using GBSSI-1A sequences or as a weakly supported sister group with Chaenomeles using GBSSI-2B and nrITS plus GBSSI-2B sequences (Campbell et al., 2007). The ER clade was moderately supported to be sister to Osteomeles and Stranvaesia using entire nrDNA repeats (Liu et al., 2019). Based on the transcriptome data, Xiang et al. (2017) placed the ER clade as sister to a large clade that includes Chaenomeles, Cydonia, Docynia, Eriolobus, Malus, Photinia, Pseudocydonia, Pyrus, Sorbus, and Stranvaesia, and this relationship does not conflict with our result. The ER clade is largely distributed in the warm temperate to subtropical and tropical regions of China, Indochina and the Malesian region, whereas most other Asian Maleae including its closely related genera are well developed in cool temperate and warm temperate regions, extending to subtropical regions of the Northern Hemisphere. Fossils of the ER clade were reported from northern China and northeastern Siberia in the Miocene when the earth was warmer (Hsu and Chaney, 1940;Baikovskaja, 1974).
The ER clade can be easily distinguished from other genera of Maleae. They have evergreen, coriaceous leaves with unlobed, entire to finely or coarsely serrate margins (Robertson et al., 1992), and semi-inferior to inferior ovaries  and its seed morphology is unique with the lack of endosperm and proportionally larger in size with a rounded or wide-elliptic cross-section Aldasoro et al., 2005). The clade was placed as sister to all pome-bearing members of Maleae based on morphological characters (Aldasoro et al., 2005), supporting its isolated position. The incongruent positions of the ER clade based on our plastome and the nrDNA data may also point to the potential impact of hybridization. So the ER clade remains an enigmatic position within Maleae.

Paraphyly of Eriobotrya
Our results strongly supported the paraphyly of Eriobotrya, with Rhaphiolepis nested within it, based on the plastome as well as the nrDNA datasets. The paraphyly of Eriobotrya was never reported in any previous studies. Yang et al. (2012) sampled only one species of Rhaphiolepis (R. indica) in their preliminary phylogenetic analyses of Eriobotrya using the nuclear Adh sequences. Rhaphiolepis indica was shown to be sister to a subclade of Eriobotrya (E. bengalensis Hook.f., E. prinoides Rehder & E.H.Wilson f. angustifolia J.E.Vidal, E. henryi, and E. seguinii), and this clade was then sister to another subclade of nine Eriobotrya species, showing the biphyly of Eriobotrya, with an extremely limited sampling of Rhaphiolepis. Yang et al. (2012) never discussed the nonmonophyly of Eriobotrya nor its taxonomic implications. Due to the general limited taxon sampling of Eriobotrya and Rhaphiolepis, most previous studies only emphasized the close relationships between these two genera (Campbell et al., 2007;Lo and Donoghue, 2012;Xiang et al., 2017;Yang et al., 2017;Zhang et al., 2017), but their data were insufficient in addressing the phylogenetic relationships between the two genera.
Our study represents the first phylogenetic analysis that was designed to test the monophyly of Eriobotrya and Rhaphiolepis with the taxon sampling representing their respective morphological and geographic diversity. The paraphyly of Eribototrya revealed in our plastid and nrDNA supports the merge of Eriobotrya and Rhaphiolepis into one genus. Eriobotrya was thought to be easily distinguished from Rhaphiolepis by the former's persistent sepals (vs. caducous sepals in Rhaphiolepis) and the excurrent lateral veins of leaves (vs. curved lateral veins in Rhaphiolepis). However, these two characters were not always stable to be used to distinguish Eriobotrya and Rhaphiolepis. For example, based on our field observations and herbarium studies in PE and US, the sepals of Eriobotrya henryi are obviously caducous ( Figure  6A) and the lateral veins of leaves in E. henryi and E. seguinii are curved ( Figures 6B, C), furthermore, the lateral veins of Rhaphiolepis sometimes terminate at the leaf margins ( Figure  6D). Treating the ER clade as one genus is also supported by two synapomorphies: the presence of rounded or widely elliptic seeds and the absence of endosperm). Geographically, these two genera overlap broadly in East and Southeast Asia ( Figure  1). The phylogenetic, morphological and geographic evidence all supports merging Eriobotrya into Rhaphiolepis, which has the nomenclatural priority. Nevertheless, describing the clade II of Figures 3 as a new genus may seem a likely alternative solution to ensure the monophyly of each genus, however, the phylogenetic incongruence (e.g., between E. obovata and E. salwinensis) in the chloroplast and nuclear trees is consistent with the extensive gene flow between clades II and III (see below). It thus seems not a good taxonomic solution to segregate clade II of Eriobotrya into a separate genus.
We detected topological incongruence within the expanded Rhaphiolepis between the chloroplast and nrDNA trees ( Figure  5). Processes that might explain the incongruence between chloroplast and nuclear phylogenies include incomplete lineage sorting, allopolyploidy, and hybridization (Rieseberg and Soltis, 1991;Mckinnon et al., 2001;Zou and Ge, 2008). Allopolyploidy has never been reported in Eriobotrya except for the cultivars of E. japonica (Mehra et al., 1973;Singhal et al., 1990;Chen, 1993;Liang et al., 1999;Liang et al., 2001;Yahata et al., 2005;Zhang et al., 2012;Li et al., 2013a), and this mechanism may be excluded for explaining the conflicts between the plastid and nrDNA trees. Lineage sorting could also result in incongruence between chloroplast and nuclear topologies, however, it is very difficult to distinguish hybridization from incomplete lineage sorting (Joly et al., 2009). Hybridization has been shown to be very common in Maleae Lo and Donoghue, 2012;Liu et al., 2019), and artificial hybrids were reported even between Eriobotrya and Rhaphiolepis (Coombes and Robertson, 2008). The conflicts shown in our results likely reflect frequent hybridization events in the evolutionary history of the expanded Rhaphiolepis. The extent and impact of hybridizations in Maleae will need to be further analyzed using next-generation sequencing and genomic tools . We will discuss the evolutionary events, involving hybridization, chloroplast capture, introgression, and/or allopolyploidy that occurred in the expanded Rhaphiolepis with a comprehensive sampling in a follow-up paper.
Below we transfer all taxa of Eriobotrya to Rhaphiolepis and make the nomenclatural changes.

AUTHOR CONTRIBUTIONS
B-BL, G-NL, D-YH, and JW designed the study. B-BL and G-NL conducted the experiments, analyzed the data, and drafted the manuscript. JW and D-YH provided suggestions on structuring the paper and the main points of the discussion, and revised the manuscript. All authors approved the final manuscript.

ACKNOWLEDGMENTS
We thank John Wiersema (US) for his advice on nomenclature and Chao Xu (PE) for assistance with experiments. Our thanks also go to Xin-Tang Ma (PE), Xiao-Hua Jin (PE), You-Sheng Chen (SCBG), Hui-Min Li (NAS), and Ting Wang (HHBG) for the sample collections and Jian Huang (HITBC) for supplying the photo of Eriobotrya henryi. The computational analyses were conducted on the Smithsonian Institution High Performance Cluster (SI/HPC, "Hydra").

SUPPLEMENTARY MATERIALS
The Supplementary Material for this article can be found online a t : h t t p s : / / w w w . f r o n t i e r s i n . o r g / a r t i c l e s / 1 0 . 3 3 8 9 / fpls.2019.01731/full#supplementary-material SUPPLEMENTARY FIGURE 1 | Gene map of the chloroplast genome of Rhaphiolepis loquata (Eriobotrya japonica). The genes inside and outside of the circle are transcribed in the clockwise and counterclockwise directions, respectively. Genes belonging to the different functional group are shown in different colors. The thick lines indicate the extent of the inverted repeats (IRa and IRb) that separate the genomes into small single-copy (SSC) and large single-copy (LSC) regions.

SUPPLEMENTARY FIGURE 2 | The phylogenetic relationships between
Eriobotrya and Rhaphiolepis in the framework of Maleae resolved by Bayesian inference of the coding sequences (CDS) of chloroplast genome. Numbers associated with the branches are ML bootstrap value (BS) and BI posterior probabilities (PP), and asterisks (*) indicate 100/1 support.
SUPPLEMENTARY TABLE 1 | Accessions of Maleae used in this study. † and ‡ indicate that the samples were extracted from silica-gel dried leaves and herbarium specimens, respectively.