Increasing Taxa Sampling Provides New Insights on the Phylogenetic Relationship Between Eriobotrya and Rhaphiolepis

Eriobotrya (Rosaceae) is an economically important genus with around 30 species. It is widely distributed in tropical and warm temperate regions of Asia, with most of its species in China, Myanmar, and Vietnam. However, Eriobotrya is often confused with the smaller genus Rhaphiolepis, and the phylogenetic relationships between the two genera are controversial. Here we present phylogenetic analyses of 38 newly generated Eriobotrya and Rhaphiolepis nrDNA together with 16 sequences of nrDNA and 28 sequences of ITS obtained from GenBank, representing 28 species of Eriobotrya and 12 species of Rhaphiolepis, in order to reconstruct highly supported relationships for the two genera. Contrary to previous research based on limited sampling, our results highlight the monophyly of Eriobotrya as well as Rhaphiolepis. The topology recovered here is consistent with key morphological synapomorphies such as the persistent sepals in Eriobotrya. Our findings show that increased sampling of taxa can provide a more robust phylogeny through reducing phylogenetic error and increasing overall phylogenetic accuracy.


INTRODUCTION
Eriobotrya Lindl. and Rhaphiolepis Lindl., two genera of the tribe Maleae in the family Rosaceae (Kalkman, 2004), include about 30 and 15 species respectively, which are distributed throughout tropical and warm temperate regions from East Asia to tropical Southeast Asia. Loquat [Eriobotrya japonica (Thunb.) Lindl.] was endemic and originally domesticated in China and has been widely cultivated throughout the world. The nutritious and fleshy fruits have an attracted increasing number of consumers worldwide (Badenes et al., 2009;Blasco et al., 2014). Rhaphiolepis indica (L.) Lindl. also has nutritious fruits, and the red pigment in the pericarp can be used as a colorant (Huang et al., 2008).
Eriobotrya japonica (Thunb.) Lindl. was first described in Flora Japonica by Thunberg (1784) under the genus name Mespilus L.. John Lindley revised the genus Mespilus and established Eriobotrya as a new genus in 1882 (Lindley, 1822). Rhaphiolepis indica (L.) Lindl. was first described in Species Plantarum (1753) under the genus name Crataegus L.. In 1820, Lindley separated Rhaphiolepis indica from the genus Crataegus because its fruits had a papery endocarp, and he later published the genus Rhaphiolepis (Lindley et al., 1820). Eriobotrya is characterized by its fruits that have persistent sepals and leaves with excurrent lateral veins, whereas the calyx of Rhaphiolepis is quickly deciduous as a unit, leaving an annular ring and leaves have curved lateral veins (Vidal, 1968;Vidal, 1970;Kalkman, 1993;Lu et al., 2003;Kalkman, 2004).
It has long been difficult to classify the genera of the Maleae tribe, which may be due to polyploidy events, rapid radiations, frequent hybridizations, and/or ancient diversification among some clades (Wolfe and Wehr, 1988;Robertson et al., 1991;Vamosi and Dickinson, 2006;Campbell et al., 2007;Dickinson et al., 2007;Li et al., 2012;Lo and Donoghue, 2012;Xiang et al., 2016;Liu et al., 2019). The latest research also shows that multiple ancient hybridization and chloroplast capture events within Eriobotrya in the Yunnan-Guizhou Plateau (Chen et al., 2021). This makes the taxonomic study of Eriobotrya more complicated. In addition, the phylogenetic status of Eriobotrya and Rhaphiolepis in the Maleae tribe has always been uncertain, and the phylogenetic relationship between the two genera is also controversial. A sister relationship between Eriobotrya and Rhaphiolepis was reported for the first time by Campbell et al. (1995) using the nuclear ribosomal internal transcribed spacer (nrITS) (Campbell et al., 1995). Six chloroplast DNA (cpDNA) regions, two GBSSI genes (1A and 2B), and nrITS sequences supported a sister relationship between Eriobotrya and Rhaphiolepis (Campbell et al., 2007). Phylogenetic relationships among 88 genera of Rosaceae were investigated using nucleotide sequence data from six nuclear and four chloroplast regions, and the results showed that Eriobotrya and Rhaphiolepis were sister groups (Li et al., 2012). In addition, nrITS data supported the monophyly of Eriobotrya, with Rhaphiolepis indica as a sister to the Eriobotrya clade (Li et al., 2012). Further studies in Maleae, also showed this sister relationship was supported using chloroplast, nrITS, and even whole plastome sequences (Lo and Donoghue, 2012;Xiang et al., 2016;Zhang et al., 2017;Sun et al., 2018;Liu et al., 2019;Idrees et al., 2020b). In a preliminarily phylogenic study of the Eriobotrya genus based on the nuclear ribosomal DNA (nrDNA) Adh sequences, R. indica was shown to be sister to a subclade of Eriobotrya and they suggested paraphyly (Yang et al., 2012). Close morphological and genetic relationships have been found in almost all studies involving Eriobotrya and Rhaphiolepis, there have been very few cases showing paraphyly and there was no case for merging them, despite the existence of intergeneric hybrids (Aldasoro et al., 2005;Sun et al., 2018).
A present genomic study, however supports incorporating Eriobotrya into the Rhaphiolepis genus and renaming all species within the Eriobotrya genus (Liu et al., 2020). In that research analyzed 16 nrDNA sequences and 21 complete plastomes indicating that the Rhaphiolepis species were nested among the Eriobotrya taxa (Liu et al., 2020). Morphologically, the persistent sepals and the excurrent lateral veins of the leaves are used to distinguish the two genera (Vidal, 1968;Vidal, 1970;Kalkman, 1993;Lu et al., 2003;Kalkman, 2004), but researchers have found that the sepals of Eriobotrya henryi Nakai are obviously caducous, the lateral veins of the leaves in E. henryi and Eriobotrya seguinii J. E. Vidal are curved, and the lateral veins of Rhaphiolepis ferruginea F. P. Metcalf sometimes terminate at the leaf margins (Liu et al., 2020). Furthermore, the seeds of the two genera are rounded or widely elliptic with an absence of endosperm (Liu et al., 2020). Geographically, these two genera overlap broadly in East and Southeast Asia (Liu et al., 2020). However, the latest research (Kang et al., 2021) does not support the results of Liu et al. These latter results show that Bayesian inference (BI) and maximum likelihood (ML) trees exhibit a well-supported monophyly for Eriobotrya, which is separate and distinct from Rhaphiolepis (Kang et al., 2021). To better estimate the phylogenetic relationship between Eriobotrya and Rhaphiolepis, it is necessary to sample more taxa to reliably reconstruct the phylogenetics of the Eriobotrya and Rhaphiolepis genera.

DNA Extraction and Sequencing
Thirty-eight taxa from Eriobotrya and Rhaphiolepis were collected, with additional 16 nrDNA and 28 nrITS sequences from GenBank added ( Table 1, and Supplementary Table S1). The species of Eriobotrya and Rhaphiolepis collected covered all main plant-distribution regions ( Figure 1). Genomic DNA was extracted from 3 g of fresh leaves or from silica-dried leaf materials using the modified cetyltrimethylammonium bromide (CTAB) method (Liu et al., 2005), in which 4% CTAB was used, and we added~1% polyvinyl polypyrrolidone (PVP) and 0.2% DL-dithiothreitol (DTT). The nrDNA was sequenced following Zhang et al. (2016), and the long-range PCR was used for next-generation sequencing with a primer pair for nuclear ribosomal DNA (rRNA_2F: TAAGCCATGCATGTG TAAGTATGAAC; rRNA_2R: CGTATTTAAGTCGTCTGC AAAGGATT). Sequencing was performed at Annoroad Gene Technology Co., Ltd., Beijing, China.

Mutation Events Analysis
To identify the microstructural mutations between Eriobotrya and Rhaphiolepis, the two genera-aligned sequences were further analyzed with MAFFT version 7 software (Katoh and Standley, 2013) and then manually aligned with MEGA X (Kumar et al., 2018). Indel and single-nucleotide polymorphism (SNP) events were counted and positioned in the nrDNA using manual statistics. We also conducted a sliding-window analysis to evaluate the nucleotide variability (Pi) throughout the nrDNA in DnaSP version 6 software (Rozas et al., 2017). The window length was set to 100 bp and the step size to 50 bp.  Table S3). To evaluate potential conflict among regions, we divided matrix I into nine subsets: ETS, ITS1, ITS2, ETS-ITS1, ETS-ITS2, ITS1-ITS2, 26S, and 18S-5.8S-26S. The sequence matrix was aligned with MAFFT version 7 software (Katoh and Standley, 2013) and then manually edited with MEGA X (Kumar et al., 2018).

Phylogenetic Analyses
Two different data matrices were aligned and analyzed using BI, ML, and Parsimony (P) methods. The BI was performed with BEAST version 2.6.3 (Bouckaert et al., 2019), using the best-fit DNA replacement model selected by jModelTest 2.1.10 for the phylogenetic reconstruction (Darriba et al., 2012). The Markov chain Monte Carlo (MCMC) algorithm was run for 100,000,000 generations, and the BI analysis started with a random tree and was sampled every 1,000 generations. The first 20% of the tree was discarded as burn-in, and the remaining tree was used to produce a majority-rule consensus tree. The ML analysis was carried out with IQ-TREE version 1.6.7 (Nguyen et al., 2015) with 1,000 bootstrap(BS) replicates using UFBoot2 (Hoang et al., 2018) and collapsing the near zero branches option. The P analysis was performed with MEGA X (Kumar et al., 2018) and 1,000 BS replicates.

Size and Organization of the nrDNA Sequences
The size of the 38 newly determined Eriobotrya and Rhaphiolepis nrDNA sequences ranged from 6,775 bp (Eriobotrya glabrescens J.E. Vidal) to 6,803 bp (E. henryi) ( Table 2). These nrDNA sequences included three rRNA genes and three transcribed spacers. In the 26S large-subunit rRNA (26S) region, the length varied from 3,358 bp (E. glabrescens) to 3,387 bp (E. henryi); in the 18S small-subunit rRNA (18S) region, 1,808 bp; in the 5.8S rRNA (5.8S) region, 159 bp; in the external transcribed spacer (ETS) region, from 1,016 to 1,019 bp, in the ITS1 region from 217 to 219 bp; and in the ITS2 region from 207 to 217 bp. The overall G + C content was 56.0-56.5%. The G + C content of the rRNA region varied from 55.0 to 55.3% (18S, 49.8-49.9%; 5.8S, 56.6-57.9%; 26S, 57.7-58.1%), and that of gene spacer region varied from 59.5 to 60.8% (ETS, 56.3-57.6%; ITS1, 64.1-66.4%; ITS2, 67.4-72.1%). Correlations were tested between the size of the nrDNA sequence and each of the six regions and the GC content of the nrDNA sequence and each of the six regions. The r 2 and Pearson results (r 2 > 0.5 and p < 0.05) were considered significant. There were significant correlations in the sequence size between nrDNA and 26S and in the GC content between nrDNA and ETS, ITS2, and 26S (Table 3).

Numbers and Pattern of Indel and Single-Nucleotide Polymorphism Mutations in nrDNA Sequences
To detect variable sites in the nrDNA of Eriobotrya and Rhaphiolepis, indel mutations among the 38 sequences were identified. A total of 21 indels were detected in the 38 sequences, including seven indels in the 26S region, six indels in the ETS region, four indels in the ITS2 region, and four indels in the ITS1 region ( Table 4). All 21 indels occurred in the nrDNA sequences of the Eriobotrya taxa, rather than that of the Rhaphiolepis taxa.
The SNP markers were also counted in the nrDNA sequences of Eriobotrya and Rhaphiolepis species. The nrDNA of E. japonica (RHA10014) was used as a reference. We detected a total of 348   Figure S1). The Ts-to-Tv ratio was 1:0.35. Among the Tv, 10 were Tv between the T and the A, 23 were Tv between the C and the G, and the other 315 were related to GC content changes.  Figure S2).
To elucidate the level of sequence divergence, the Pi values within 100 bp in the nrDNA of both Eriobotrya and Rhaphiolepis, were calculated with DnaSP 6.0 software (Figure 2). Within the combined Eriobotrya and Rhaphiolepis genera, those values varied from 0 to 0.08315, with a mean of 0.01988 ( Figure 2A). Within Eriobotrya, those values varied from 0 to 0.05959, with a mean of 0.00665 ( Figure 2B). Within the Rhaphiolepis, those values varied from 0 to 0.01709, with a mean of 0.001812 ( Figure 2C). The results show that the differences between the two genera were larger than those among congeneric species. Three regions including ETS, ITS1, and ITS2 were particularly highly variable between Eriobotrya and Rhaphiolepis and among the congeneric species.
To better understand the phylogenetic relationships among the sequenced taxa from Eriobotrya and Rhaphiolepis, we

Phylogenetic Analyses Based on Six Regions of nrDNA Sequences
Incongruence is significant among the topologies obtained from the transcribed spacer regions and rRNA gene sequences. BI analyses of the ETS, ITS1, ITS2, ETS-ITS1, ETS-ITS2, ITS1-ITS2, and ETS-ITS1-ITS2 sequences fully resolved phylogenetic relationships among the major clades and most genera, and the Eriobotrya and Rhaphiolepis groups had high internal support, the exception being the data matrices of 18S-5.8S-26S and 26S ( Figure 5). The phylogenetic analyses with 18S-5.8S-26S and 26S gene sequence found that two Eriobotrya species E. henryi and E. seguinii and all Rhaphiolepis species are located in the same clade, and other Eriobotrya species forming the next sister group, followed by Phippsiomeles.

nrDNA Sequence Variation
Rapidly developing molecular markers, such as allozymes, DNA sequence including single nucleotide polymorphism (SNP), and simple DNA sequence repeated (SSR) loci have great potential in species identification, population structure analysis, and phylogenetic analysis. The standardized DNA regions include plastid rbcL, matK, and trnH-psbA and ribosomal DNA ITS1 or ITS2 (China Plant BOL Group et al., 2011). Among the 38 ITS1 sequences of Eriobotrya and Rhaphiolepis species, we manually identified mutation events including 28 SNPs and four SSR indels, and 61 SNPs, one SSR indel, and four non-SSR indels were accurately located in the ITS2 sequences. In addition, two non-SSR indels, four SSR indel, and 151 SNPs were found in the ETS regions, and four highly variable regions including ITS2, ETS, 26S, and ITS1 among the Eriobotrya and Rhaphiolepis species were identified. Both ITS1 and ITS2 regions were used to elucidate relationships among the taxa of Eriobotrya and    Rhaphiolepis (Idrees et al., 2018;Idrees et al., 2020a;Idrees et al., 2020b;Kang et al., 2021). Here, two rarely reported highly variable loci ETS and 26S were present in Eriobotrya and Rhaphiolepis nrDNA sequence (Figures 2, 5). It was stressed that complementary ETS and 26S markers to the recommended ITS1 and ITS2 should continue to be assessed from nrDNA sequence. Through analysis of nrDNA sequences, additional plant ETS and 26S have been found and, in turn, have become valuable molecular markers for the identification of interspecific germplasm, which is helpful for the phylogeny relationship.

Relationship of Eriobotrya and Rhaphiolepis
All ML, BI, and MP analyses of the nrDNA sequences fully resolved phylogenetic relationships between Eriobotrya and Rhaphiolepis and confirmed that the monophyly of the Eriobotrya and Rhaphiolepis, respectively, in agreement with previously published phylogenetic relationships (Yang et al., 2012;Yang et al., 2017;Idrees et al., 2018;Idrees et al., 2020a;Idrees et al., 2020b;Kang et al., 2021). The topology obtained shows that nrDNA sequence, with appropriate sampling, can provide robust and significantly supported relationship among deep lineages of Eriobotrya and Rhaphiolepis. Seven such phylogenetically meaningful clades were identified among the deep lineages of the Eriobotrya. The backbones of the phylogenomic topologies obtained here are consistent with previously published phylogenetic relationships (Kang et al., 2021), but problems within several major clades in the Eriobotrya were solved. All the species in the genus Rhaphiolepis form a sister group of Eriobotrya, consistent with the study by Chen et al. (2020) and Kang et al. (2021). The Vietnamese species E. condaoensis is located in the earliestdiverging extant lineage within the Eriobotrya, which is in agreement with the previous phylogenetic results by Kang et al. (2021) who defined the relationships among 17 Eriobotrya species, respectively. This species is located in Con Dao National Park in southern Vietnam (Idrees et al., 2018), relatively far away from E. henryi and E. seguinii (Supplementary Figure S8). In the Clade 2, the sisterhood of E. henryi and E. seguinii is clarified, as found in previous studies (Idrees et al., 2020a;Idrees et al., 2020b;Liu et al., 2020;Kang et al., 2021). Previous phylogenetic analyses with plastid genome found that members of E. laoshanica were sister to E. malipoensis (Chen et al., 2020). However, our phylogenetic analyses show both species are located in different clades (Figures 3, 4). E. laoshanica is located in the Clade 4, while E. malipoensis is located in the Clade 3 with E. grandiflora, E. hookeriana, E. japonica, E. petiolate, E. sp2, and E. × daduheensis. E. deflexa, E. cavaleriei and E. fragrans are located in the Clade 5, likewise significant support in the ITS data (Kang et al., 2021) and the nuclear genes data (Chen et al., 2021) rather than the study of Idrees et al. (2020a); Idrees et al. (2020b), Chen et al. (2020), and Liu et al. (2020). In clade 6, E. prinoides is closely related to E. elliptica and E. serrata, but the relationship is not supported in the study of Idrees et al. (2020a); Idrees et al. (2020b); Kang et al. (2021), and Chen et al. (2021). Three Myanmar Eriobotrya species E. glabrescens, E. platyphylla, and E. sp4, were located in Clade 7 with six Chinese Eriobotrya species, E. bengalensis, E. bengalensis var. angustifolia, E. obovate, E. salwinensis, E. sp3, and E. tengyuehensis (Figures 3, 4). We further determined the relationships of 17 Eriobotrya and Rhaphiolepis species, E. elliptica, E. glabrescens, E. laoshanica, E. platyphylla, E. sp1, E. sp2, E. sp3, E. sp4, E. sp5, E. × daduheensis, R. brevipetiolata, R. indica var tashiroi, R. mekongensis, and R. sp.
Our nrDNA sequences of Eriobotrya and Rhaphiolepis yielded a fully resolved tree, consistent with the study of Kang et al. (2021), rather than that of Liu et al. (2020). The phylogenomic analysis showed E. henryi and E. seguinii is not nested among the members of Rhaphiolepis, which is incompatible with the chloroplast and nrDNA data in Liu et al. (2020). Liu et al. (2020) collected 16 taxa from Eriobotrya and Rhaphiolepis, and the molecular, morphological, and geographic evidence supported merging these two genera into one genus. In that research, the sampling proportions of Eriobotrya and Rhaphiolepis accounted for 27.6 and 70%, respectively, that is, only a small proportion of the sample was Eriobotrya. In addition, according to the distribution of Eriobotrya and Rhaphiolepis, no samples were collected from southeast Asia, southern Yunnan, Hainan, Sichuan, and Tibet. Among the 7 clades of Eriobotrya in our phylogenetic tree, Liu et al. (2020) only sampled species in clade 2, 3, 5, and 7. It is known that sample deviation could lead to phylogenetic errors, increasing the sampling of taxa is one of the most important ways to increase overall phylogenetic accuracy (Hillis, 1996;Hillis, 1998;Graybeal, 1998;Poe, 1998;Soltis et al., 1999;Pollock and Bruno, 2000;Pollock et al., 2002;Murphy et al., 2001;Zwickl and Hillis, 2002). Research shows that increased taxon sampling provides new insights into the phylogeny and evolution of the subclass Calcaronea (Porifera, Calcarea) (Alvizu et al., 2018). In our research, adequate sampling of Eriobotrya species required sampling from Myanmar, Vietnam, Yunnan, Sichuan, Tibet, Hainan, and other places, greatly increasing the taxa and reducing the phylogenetic errors ( Figure 1). In addition, we calculated the proportions of SNP mutations in the two genera of nrDNA. Among 23 species of Eriobotrya and 11 species of Rhaphiolepis, shared SNP sites accounted for 19.25% (Supplementary Figure S6A), whereas the shared sites accounted for 22.42% in the 8 species of Eriobotrya and 7 species of Rhaphiolepis (Supplementary Figure S6B). According to the distribution heat map of the two genera, although there are some overlapping areas, Eriobotrya species are mainly distributed in southwestern China and Indo-China (Supplementary Figure S7A), while Rhaphiolepis species are centered in southeastern China, and very scarcely recorded in Yunnan, Sichuan and Myanmar (Supplementary Figure S7B). Yunnan is the diversity center of Eriobotrya species, while the diversity center of Rhaphiolepis species is not.
The phylogenetic tree obtained with the 18S-5.8S-26S and 26S dataset showed a low resolution at all taxonomic levels, rendering most relationships inconclusive, which may be caused by the conservation of the rRNA gene functions. In analyzing the hypervariable regions of the two genera, we found that the mutation frequency was low in 18, 5.8, and 26S ( Figure 2). Because of the conservatism of the rRNA gene sequencing, combined with a low mutation rate and limited information loci. In addition, there is an overlap between E. henryi, E. seguinii and Rhaphiolepis species, and hybridization may have occurred between E. henryi, E. seguinii and Rhaphiolepis species (Supplementary Figure S8). Resulting in the insertion of E. henryi and E. seguinii into the genus Rhaphiolepis. It is the same not only in plants but also in animals. This is due to the small number of informative sites in the18S rRNA. 18S proved to be highly conserved within Calcaronea and does not have sufficient signal to resolve phylogenetic relationships within the subclass (Alvizu et al., 2018).

Morphological Difference Between Eriobotrya and Rhaphiolepis
In the taxonomic literature and flora, the persistence of sepals was used to distinguish between Eriobotrya and Rhaphiolepis (Vidal, 1968;Vidal, 1970;Kalkman, 1993Kalkman, , 2004Lu et al., 2003). However, Liu et al. (2020) found that the sepals of E. henryi fell early in the field, and it was considered that the persistent sepal could not be used to distinguish between Eriobotrya and Rhaphiolepis. In addition, those authors argued that the camptodromous leaf venation in some loquat species of Eriobotrya and Rhaphiolepis lacked stability (Liu et al., 2020). However, taxonomic studies clearly show that both camptodromous and craspedodromous venation can be observed in Eriobotrya (Robertson et al., 1991;Gu and Spongberg, 2003.). We found that the sepals of E. henryi were persistent in the field ( Figure 6D), and there was a significant difference between the two genera ( Figure 6, Supplementary Table S5 and Supplementary Figure S9). Kang et al. (2021) reviewed the same picture and found that, although the E. henryi fruit calyx was caducous, it was intact in some photos. The (circular) annular ring after sepal senescence can only be observed in Rhaphiolepis (Robertson et al., 1991;Gu and Spongberg, 2003;Kang et al., 2021). The two genera can be well separated according to whether the sepals fall off and whether there is an annular ring after sepal senescence or not (Supplementary Table S5 and Supplementary Figure S9). In addition, Shaw (2020) stressed the importance of maintaining nomenclatural stability for Eriobotrya species with horticultural and agricultural value. Because of the conflicting issues we have found, we do not recommend Eriobotrya being incorporated into the genus Rhaphiolepis.

CONCLUSION
Phylogenetic analysis of nrDNA sequences strongly supports Eriobotrya and Rhaphiolepis being monophyletic. In addition, phylogenetic analysis using nrDNA combined with ITS sequences, both the Eriobotrya and Rhaphiolepis were 100% supported monophyletic. Moreover, we speculate that the phylogenetic evidence for Eriobotrya as monophyly is congruent with the morphological characteristics of its leaves and the persistence of its sepals. It is not recommended that Eriobotrya be merged into Rhaphiolepis.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ACKNOWLEDGMENTS
Our thanks go to Wei Kang, Qunfei Yu, Fayu Feng, Linyi Yang, and Yaya Qu for the sample collections. We also thank Wei Sun, Renbin Zhu, Xiao Luo and Big Data Platform for Ex Situ Plant Conservation (https://espc.cubg.cn/) for supplying the photos of Eriobotrya henryi and Eriobotrya bengalensis var. angustifolia, Rhaphiolepis indica, Rhaphiolepis major, and Rhaphiolepis umbellate.