Combined Analyses of Chloroplast DNA Haplotypes and Microsatellite Markers Reveal New Insights Into the Origin and Dissemination Route of Cultivated Pears Native to East Asia

Asian pear plays an important role in the world pear industry, accounting for over 70% of world total production volume. Commercial Asian pear production relies on four major pear cultivar groups, Japanese pear (JP), Chinese white pear (CWP), Chinese sand pear (CSP), and Ussurian pear (UP), but their origins remain controversial. We estimated the genetic diversity levels and structures in a large sample of existing local cultivars to investigate the origins of Asian pears using twenty-five genome-covering nuclear microsatellite (simple sequence repeats, nSSR) markers and two non-coding chloroplast DNA (cpDNA) regions (trnL-trnF and accD-psaI). High levels of genetic diversity were detected for both nSSRs (HE = 0.744) and cpDNAs (Hd = 0.792). The major variation was found within geographic populations of cultivated pear groups, demonstrating a close relationship among cultivar groups. CSPs showed a greater genetic diversity than CWPs and JPs, and lowest levels of genetic differentiation were detected among them. Phylogeographical analyses indicated that the CSP, CWP, and JP were derived from the same progenitor of Pyrus pyrifolia in China. A dissemination route of cultivated P. pyrifolia estimated by approximate Bayesian computation suggested that cultivated P. pyrifolia from the Middle Yangtze River Valley area contributed the major genetic resources to the cultivars, excluding those of southwestern China. Three major genetic groups of cultivated Pyrus pyrifolia were revealed using nSSRs and a Bayesian statistical inference: (a) JPs; (b) cultivars from South-Central China northward to northeastern China, covering the main pear production area in China; (c) cultivars from southwestern China to southeastern China, including Yunnan, Guizhou, Guangdong, Guangxi, and Fujian Provinces. This reflected the synergistic effects of ecogeographical factors and human selection during cultivar spread and improvement. The analyses indicated that UP cultivars might be originated from the interspecific hybridization of wild Pyrus ussuriensis with cultivated Pyrus pyrifolia. The combination of uniparental DNA sequences and nuclear markers give us a better understanding of origins and genetic relationships for Asian pear groups and will be beneficial for the future improvement of Asian pear cultivars.


INTRODUCTION
Pear (Pyrus L.), as an important temperate fruit tree worldwide, has been cultivated for more than 3,000 years (Pu and Wang, 1963;Janick, 2000). Pyrus species are naturally distributed in the Eurasian continent and are geographically divided into two native groups: occidental and oriental pears (Bailey, 1917). The occidental pears, including over 20 species, occur in Europe, Northern Africa, Asia Minor, Iran, Central Asia, and Afghanistan, and the majority of cultivars in these areas originated from Pyrus communis L. (Rubtsov, 1944). The oriental pears, comprising 12-15 species, are distributed from the Tian Shan and Hindu Kush Mountains eastward to Japan and are largely concentrated in East Asia, including China, Japan, and Korea (Rubtsov, 1944;Challice and Westwood, 1973). The commercial pear cultivars native to East Asia can be mostly divided into four groups, Ussurian pear (UP; P. ussuriensis Maxim.), Chinese white pear (CWP; Pyrus pyrifolia White Pear Group, generally known as P. x bretschneideri Rehder), Chinese sand pear (CSP; P. pyrifolia Nakai), and Japanese pear (JP; P. pyrifolia Nakai), based on their geographic distributions, and morphological and physiological traits (Teng and Tanabe, 2004). In the limited areas of Gansu, Qinghai, and Xinjiang Provinces, Xinjiang pear (P. x sinkiangensis Yu) is cultivated, which is a hybrid of P. communis and P. pyrifolia (Teng et al., 2001). Asian pears account for over 70% of the worldwide total in terms of both production volume and cultivated area (FAO).
The four main cultivated pear groups native to East Asia, except JP, are native to China, and the total number of cultivars and landraces has been estimated as ∼3,000 in China (Pu and Wang, 1963). UP is the most cold-tolerant cultivar, mainly cultivated in northeast China, and is generally considered to be directly domesticated from the wild species P. ussuriensis. Its fruit are usually small globose or oblate with a persistent calyx, and they need an after-ripening process to become soft and edible like the European pear (Pu and Wang, 1963;Teng et al., 2002). CSP is distributed in southern China and was derived from wild P. pyrifolia in the Changjiang (Yangtze) River Valley of China. Its fruit, with smooth or russet skin, vary greatly in shape, size, and shelf life (Pu and Wang, 1963). CWP grows only between the areas that produce UP and CSP. Therefore, CWP was once proposed to have originated from the hybridization of P. ussuriensis and P. pyrifolia (Rubtsov, 1944;Kikuchi, 1946). In China, this group of cultivars has generally been assigned to P.
x bretschneideri Rehder (Teng and Tanabe, 2004). The latter was named by Rehder (1915) based on seeds sent by Dr. Bretschneider from China. Some researchers speculated that P. x bretschneideri might have originated from the natural hybridization between Pyrus betulifolia Bunge and the large-fruited cultivars grown in North Hebei Province (Kikuchi, 1946;Challice and Westwood, 1973). Thus, P. x bretschneideri is different in the morphological characteristics from the cultivars of CWP (Kikuchi, 1946(Kikuchi, , 1948. Therefore, CWP cultivars should not have originated from P. x bretschneideri (Kikuchi, 1946(Kikuchi, , 1948. The cultivated pears native to Japan are referred to as JP and they have been considered to be from the same germplasm as P. pyrifolia (Teng and Tanabe, 2004). However, the origin of JP and its relationship with CSP are long standing controversies. Some Japanese horticulturists believe that JP cultivars were domesticated from wild P. pyrifolia that once grew in southern Japan (Kikuchi, 1948). However, there has been no solid evidence for the existence of wild P. pyrifolia in Japan (Kajiura and Suzuki, 1980). Others proposed that the progenitors of JP cultivars might have been introduced from China during ancient times (Shimura, 1988).
Over the last few decades, genetic analyses based on DNA markers, such as RAPD (Random Amplified Polymorphic DNA; Teng et al., 2001Teng et al., , 2002, AFLP (Amplified Fragment Length Polymorphism; Bao et al., 2008), and SSR (Simple Sequence Repeats; Bao et al., 2007;Iketani et al., 2012;Nishio et al., 2016) have advanced our understanding of the genetic diversity and genetic relationships among cultivated pears native to East Asia. CWP is closely related to CSP (Teng et al., 2002), and may be treated as an ecotype of P. pyrifolia: P. pyrifolia White Pear Group (Bao et al., 2007). Also, JP shares some identity with pears from China (Iketani et al., 2012), especially cultivars from Zhejiang Province (ZJ) and Fujian Province (FJ) (Teng et al., 2002;Bao et al., 2007;Song et al., 2014), which was supported by our recent studies based on retrotransposons markers (Jiang et al., 2015. Thus, the three major pear cultivar groups (CSP, CWP, and JP) native to East Asia appear to have evolved from the same germplasm: wild P. pyrifolia in China (Jiang et al., 2015. Unfortunately, no population of wild P. pyrifolia still exists. Therefore, we cannot elucidate the origins and levels of geographic differentiation of CSP, CWP, and JP by comparing their genetic makeup with those of wild progenitors, as achieved with apple (Cornille et al., 2012;Nikiforova et al., 2013) and grape (Myles et al., 2011). In addition, most of the aforementioned studies were based on either relatively small or geographically limited samples that under-represented the rich diversity of cultivated pears native to East Asia. There are still some questions that need to be addressed: Where were the P. pyrifolia initially cultivated? How did pear cultivars genetically diverge during the spread of the germplasm? Are there varied genetic differentiation levels among different cultivar groups, especially the three cultivated P. pyrifolia groups? Thus, additional studies are required to answer these questions to better understand the origins and dissemination route of cultivated Asian pears. A better understanding of genetic diversity of existing local cultivars will be helpful for developing new cultivars with the combination of good fruit quality and resistance to abiotic and biotic stresses, and providing information for efficient conservation of pear germplasm resources.
In this study, 441 pear accessions from different geographical regions were collected across China and Japan, covering the major distribution areas of pear cultivars or landraces. Two hypervariable non-coding chloroplast DNA regions, trnL-trnF and accD-psaI, which have been widely applied in molecular phylogeny and genetic diversity studies of the genus Pyrus (Kimura et al., 2003;Liu et al., 2013;Wuyun et al., 2013;Zheng et al., 2014;Zong et al., 2014), were adopted to elucidate the evolutionary history of pear cultivars using a phylogeographic approach. Meanwhile, markers of genome-covering nuclear loci containing simple sequence repeats (nSSRs; microsatellites) were employed to investigate the genetic structure of East Asian cultivated pears using a Bayesian statistical inference based on the theory of population genetics that used no prior information on sample origins. This has been widely used for analyzing genetic structures of fruit crops, such as pear (Iketani et al., 2012), apple (Robinson et al., 2001;Cornille et al., 2012), and grape (Bacilieri et al., 2013). The combination of data from maternally inherited genetic markers (cpDNA sequences) and parentally inherited nuclear markers (nSSRs), will help us achieving a comprehensive understanding of genetic diversity of Asian pears and provide important clues on the origin and dissemination route of Asian pear cultivars.

Plant Materials and DNA Extractions
A total of 441 pear accessions collected from different geographic areas in East Asia were analyzed in this study (Supplementary  Table S1). All cultivated P. pyrifolia accessions (n = 401), including 221 CSPs and 125 CWPs, and 55 JPs were assigned to 14 geographic populations based on their geographical origins, of which CSPs included seven geographical populations of cultivated pears from south of Yangtze river in China, CWPs included five geographical populations from north of Yangtze river in China, and JPs contained two geographical populations from Japan, and the population code of geographical populations were described in Supplementary Table S1. Forty P. ussuriensis accessions, including 22 cultivated P. ussuriensis (UPs) and 18 wild P. ussuriensis (WUPs) sampled from North and Northeast China, were also used for the genetic analysis. Most of the local cultivars were sampled from different pear germplasm repositories in China and Japan, and a few were collected in situ (Supplementary Table S1). DNAs of all pear accessions were extracted from young leaves using a modified cetyl trimethyl ammonium bromide protocol (Doyle, 1987).

Nuclear Microsatellite Amplification
The pear accessions were genotyped using 25 polymorphic nSSR loci distributed across the genome  Supplementary Table S2). The 25 nSSR loci were selected based on the degree of applicability and number of polymorphisms in different Pyrus species. PCR amplification of microsatellite primers, including BGT23b, KU10, KA14, 28f4, 02b1 and "MES", were performed as described in previous protocols (Guilford et al., 1997;Yamamoto et al., 2002a,b;Yao et al., 2010), with fluorescent dye-labeled forward primers. The other recently developed nSSR primers ("NAU, " "CTG, " and "TXY") were used in an economic method for the fluorescent labeling of PCR products (Schuelke, 2000), and the amplification programs were performed using protocols previously described by Yue et al. (2014). Genotyping was performed on an ABI 3700XL Genetic Analyzer (Applied Biosystems, United States), with an internal size standard (GeneScan TM 500 LIZ, Applied Biosystems). Alleles were scored using GENEMAPPER 4.0 software (Applied Biosystems).

Genetic Analysis of Microsatellites
Genotyping errors and null alleles were confirmed by MICROCHECKER (Van Oosterhout et al., 2004). The number of alleles (N), number of observed alleles (Na), number of effective alleles (Ne), Shannon's Information Index (I), expected heterozygosity (H E ), and observed heterozygosity (H O ) were calculated for each nSSR locus and geographic populations by GENEPOP4.0 (Raymond and Rousset, 1995). The number of unique alleles was estimated by GenAIEx 6.5 (Peakall and Smouse, 2012). Analyses of molecular variance (AMOVA) were conducted by Arlequin 3.5 (Excoffier and Lischer, 2010) with standard AMOVA computations and 2,000 permutations using nSSR data to assess the levels of genetic differentiation among cultivar groups and among populations with different groupings. Allele frequency-based pairwise F ST index with 1,000 permutations was calculated by use of Arlequin 3.5 (Excoffier and Lischer, 2010), to describe the nuclear genetic differentiation among different populations. The software STRUCTURE2.3.3 (Pritchard et al., 2000) was applied to the multilocus nSSR data of East Asian pear accessions to infer the different gene pools (K) in the data set. Genotypic clustering was implemented in STRUCTURE using the admixture model with independent allele frequencies without using sites or populations of individuals as priors as described by Richards et al. (2009), and 10 independent runs were performed for each K in the range of 2-16 with 200,000 MCMC (Markov Chain Monte Carlo) iterations after a burn-in of 200,000 steps. The results were uploaded to the STRUCTURE HARVESTER web site (Earl, 2012) to estimate the optimal value of K based on Evanno et al.'s (2005) criterion of K, which used the rate of change in the log probability of data between successive K-values. In addition, CLUMPP (Jakobsson and Rosenberg, 2007) software was employed to average the 10 independent simulations and DISTRUCT (Rosenberg, 2004) was applied to illustrate the results graphically. To further clarify the genetic diversity of pear accessions in each geographic population, pear genotypes were assigned to a certain group when 85% or more of their inferred genome belonged to that group. The others, with lower scores, were considered admixed genotypes as described by Bacilieri et al. (2013).

Amplification and Alignment of Chloroplast Non-coding DNA Regions
The trnL-trnF and accD-psaI intergenic regions were amplified and sequenced for all pear accessions listed in Supplementary  Table S1. Amplification and sequencing conditions of the two cpDNA fragments were as described in a previous study (Zheng et al., 2014). Two cpDNA regions were combined for haplotype identification and phylogeographic analysis. The haplotypes of cpDNA fragments were identified and named following a previous phylogenetic study in Pyrus (Zheng et al., 2014). Chloroplast DNA sequences were aligned by ClustalX 1.81 (Thompson et al., 1997) and refined manually. Each indel (insertion or deletion) was treated as a single mutation event, regardless of its size, and coded as a substitution (A or T) owing to their equal probability levels .

Genetic Analyses of Chloroplast DNA
Haplotype number (H), haplotype diversity (Hd), number of segregating sites, nucleotide diversity, and the average number of nucleotide differences (K) were calculated by DnaSP 5 (Librado and Rozas, 2009). AMOVA were conducted to assess the levels of genetic differentiation among and within different groupings using 20,000 permutations of standard AMOVA computations by Arlequin3.5 (Excoffier and Lischer, 2010). Chloroplast DNA genetic differentiation was examined for geographical populations by a classical analysis of variance calculating haplotype frequency-based F ST , and the existence of a significant genetic differentiation among populations was determined by permutation analyses using 1000 randomly permuted data sets with Arlequin3.5 (Excoffier and Lischer, 2010). A Bayesian clustering of cpDNA data was conducted by BAPS6.0 (Corander et al., 2007). The maximum number of populations was set as 10 replicates for 2-16 during mixture analyses. The other parameters for admixture analyses based on the mixture analyses were set as default.

Phylogeographic Analyses of cpDNA Haplotypes
Two differentiation indices (G ST and N ST ) were computed for haplotypes in geographic populations and compared by a test with 1,000 permutations using the program Permut 2.0 (Pons and Petit, 1996). G ST is measured by simply making use of the allelic frequency to explain the genetic variation among all populations. N ST is considered to represent the similarity levels between the haplotypes, explaining the genetic differentiation influenced by both haplotype frequencies and genetic distances between haplotypes (Pons and Petit, 1996). There would be a phylogeographic structure among geographic populations, if the value of N ST is significantly greater than the value of G ST (Pons and Petit, 1996). Haplotype networks were constructed for the combined cpDNA haplotypes of cultivated pears using the computer program TCS 1.21 (Clement et al., 2000) with the 95% statistical parsimony criterion. The frequency of each haplotype in the TCS network was used to estimate haplotype outgroup probabilities, which correlated with the haplotype age (Clement et al., 2000). Neighbor-Net splits graphs (NN graphs) of cpDNA haplotypes were constructed using SplitsTree (Huson and Bryant, 2006) with the criterion set to the uncorrected p distance. Haplotype sequences of wild species Pyrus pashia 3 (tH4aH12), P. pashia 7 (tH5aH14), and P. betulifolia 1 (tH3aH7), as well as haplotypes of outgroup Malus domestica 'Ralls' (tH19aH34) and M. rockii (tH17aH33) (Zheng et al., 2014) were included in the NN graph analysis.

Estimation of the Dissemination Route of Cultivated P. pyrifolia
To further explore the dissemination route of cultivated P. pyrifolia, pear accessions from 14 geographic populations were assigned into 10 Pops, according to their geographic origins and genetic similarities. The dissemination route of cultivated P. pyrifolia was estimated from cpDNA haplotypes by approximate Bayesian computation using the software DIYABC (Cornuet et al., 2014). Two estimates, including a direct estimate and a logistic regression estimate, were computed to determine the more supported scenario of the evolutionary route of cultivated P. pyrifolia in East Asia. The former estimate simply represents the number of times that a certain scenario was detected in the data sets that were simulated from the historical models. The latter was a polychotomic weighted logistic regression, which was performed on the simulate data sets with dependent variables from a proportion of certain scenarios in the simulated data sets and independent variables from the differences of summary statistics between observed and simulated data sets (Cornuet et al., 2014).

Diversity of Asian Pears Inferred From nSSR Loci
Twenty-five nSSR loci displayed a high level of polymorphism in Asian pears with an average of 24 alleles per locus and ranging from 13 to 43 alleles (Supplementary Table S2). Genome-derived nSSRs generally displayed higher polymorphism levels than ESTderived or transcriptome-derived nSSR markers. The H E and H O values of geographical cultivar populations ranged from 0.660 (JPa) and 0.551 (JPb), respectively, to 0.792 (UP) and 0.665 [Sichuan Province (SC)], respectively, and the overall means of East Asian pears were 0.744 and 0.603, respectively ( Table 1). Geographic populations of Fujian Province (FJ; 10) and Yunnan Province (YN; 11) displayed the largest numbers of unique alleles. The inbreeding coefficient ranged from 0.083 to 0.285, which was close to zero and consistent with the self-incompatibility system of pears.
Analysis of molecular variance based on nSSR data was analyzed for pear cultivars and wild P. ussuriensis accessions native to East Asia (Supplementary Table S3). For different groupings of pear accessions, most of the total variation (88.91-93.18%) was partitioned within geographic populations. The highest variation (10.31%) among populations was observed among P. ussuriensis accessions (UP and WUP). Comparatively, the least of total variation (4.32%) was identified among populations within cultivar groups. Based on the F ST index value (Wright, 1978), levels of nuclear genetic differentiation were generally low for comparisons among adjacent geographical populations of cultivated pears (0.015-0.0621), but fairly high between populations of cultivar groups and the group of wild P. ussuriensis accessions (0.10306-0.18228) (Supplementary Table S9).

Genetic Structure of Asian Pears
STRUCTURE results for the nSSR dataset (n = 441), based on the likelihood output of the genetic admixture analysis and Evanno's K statistics indicated that K = 4 was the most pertinent levels of population subdivision (Supplementary Figure S1). At K = 4 ( Figure 1A), two JP populations were dominated by a gene pool, whereas CSP populations were mainly comprised of two gene pools, and CWP populations were also dominated by two gene pools. The cultivated UPs were composed of one main gene pool, and two minor gene pools, whereas WUP accessions were predominantly comprised of a single gene pool.
Using a threshold of >85% for genetic group assignation, a total of 195 non-admixed genotypes (out of 441 pear accessions) were assigned to four genetic groups (Supplementary Table  S4). The other 246 cultivars were indicated as the admixed genotypes, which accounted for 58.05% of the total number of pear accessions. As seen in Figure 1B, the proportion of non-admixed genotypes displayed great differences among the geographic-based populations and was generally higher in northern populations [i.e., Shandong, Henan, and Hebei Provinces (SHH), and Jilin and Liaoning Provinces (JL)], and JP populations [western Japan (JPa) and northern Japan (JPb)]. Comparatively, admixed genotypes accounted for the largest proportions of geographic-based cultivar populations of southwestern China, northwestern China and the southernmost China, such as SC (89.66%), GQ (89.66%), and Guangdong and Guangxi Province (GG; 82.86%). The proportion of admixed genotypes generally decreased as the distance increased from these western and mid-southern regions, especially in the JP population that had a very low admixed genotype level (3.45%), which was generally in accordance with the lower genetic diversity detected in these areas ( Table 1).
The identified four genetic groups of Asian pear accessions displayed a distinct geographic pattern for their genetic structures ( Figure 1B). The first group, K-4.1 (Supplementary Figure  S2), contained almost all JP genotypes (45 of 55 cultivars), including 27 cultivars of JPa and 18 of JPb, as well as three cultivars from Zhejiang province (ZJ) and GQ, and one from FJ. The second group, K-4.2, was mainly composed of genotypes from the Guizhou province (GZ; 22), and southernmost population in GG (7) and FJ (11). The third group, K-4.3, included pears belonging to CWP and CSP mainly from mid-southern, eastern and northern populations, such as HHJ, Anhui and Jiangsu Provinces (AJ), Shaanxi and Shanxi Provinces (SS), SHH, and JL. The fourth group, K-4.4 included 36.36% of the UP cultivars and all the WUP accessions.

Classification and Diversity of cpDNA Haplotypes
Sequence lengths of the combined alignments of trnL-trnF and accD-psaI ranged from 1,486 to 1,778 base pairs. In total, 10 nucleotide substitutions (5 in trnL-trnF and 5 in accD-psaI; Supplementary Table S6) and 8 insertion and deletion (indels; 2 in trnL-trnF and 6 in accD-psaI) were identified in combined cpDNA alignments, leading to 23 unique haplotypes (H1-H23) in 441 East Asian pears (Supplementary Table S1), including 7 tH haplotypes detected by sequences of trnL-trnF and 11 aH haplotypes discovered by alignments of accD-psaI. A high level of Hd (0.792) was detected in total pear accessions ( Table 2). The number of segregating sites ranged from 6 to 11. The nucleotide diversity ranged from 1.2 (FJ) to 1.76 (SS) and the K-value among whole sequences is 3.097, ranging from 1.634 (WUP) to 3.122 (SS). Of these 23 haplotypes, 19 haplotypes were detected in cultivated pears and four haplotypes (H20, H21, H22, and H23) only in the WUP accessions (Supplementary Table S7). In total, 16 haplotypes were presented in CSP, 9 in CWP, 6 in WUP, 5 in JP and 5 in UP. Seven (H4, H5, H8, H11, H15, H17, and 18) were discovered only in CSP, and H19 and H1 were discovered only in CWP and UP, respectively. Comparatively, no unique haplotype was detected in JP. Among the 19 cpDNA haplotypes from cultivated pears, haplotypes H3, H6, H7, and H14 had global frequencies in cultivated pears, covering all four Asian cultivar groups and accounting for 96.36% in JP, 91.2% in CWP, 82.8% in CSP, and 81.73% in UP of the total haplotypes, respectively (Supplementary Table S8). CSP and CWP shared seven haplotypes H3, H6, H7, H10, H12, H13, and H14; UP shared the common haplotypes H14 and H6 with CSP, CWP, and JP. Contrary to the cultivars of UP, which were domesticated directly from WUP, the most common haplotype H20 (tH1aH6) in WUP has not been detected in any cultivated UP. Cultivated UPs only shared the haplotype H1 with WUP accessions. Additionally, the rare haplotype H2 was shared by CWP and JP, but absent in CSP; H9 was shared by CSP and JP, but absent in CWP; H16 was shared by CSP and UP, but absent in JP and CWP.
Analyses of molecular variance analyses based on cpDNA data indicated there was 42.61-84.31% of total genetic variation partitioned within geographic populations (Supplementary Table  S3). Higher levels of variation were observed among species (19.18%) or among populations (57.39%) in the groupings including wild P. ussuriensis accessions (WUP) than that were identified among cultivar groups [Cultivated P. pyrifolia and cultivated P. ussuriensis, 7.81%; cultivar groups of P. pyrifolia (JP, CSP and CWP), 3.64%]. The value of F ST index for comparisons between populations of cultivated pears and WUP were generally high (0.554 to 0.657; Supplementary Table S9). They were fairly high (0.1251 to 0.4319) between geographical populations from southwestern (i.e., YN, SC, and GZ) and populations from other regions, including mid-southern (i.e., HHJ), eastern (i.e., ZJ) and northern (i.e., JL), as well as from Japan (JPa and JPb). In contrast, the values were small (−0.065 to 0.089) for comparisons between populations from southwestern or between populations from mid-southern, eastern, northern and Japan.

Phylogeography of cpDNA Haplotypes
The number of haplotypes in geographic populations ranged from 3 (JPb) to 8 (YN and SC) (Supplementary Table S7). The two differentiation indexes showed lower values (N ST = 0.202; G ST = 0.162). A permutation test demonstrated that N ST was significantly greater than G ST , indicating the existence of a phylogeographic structure in East Asian pear geographic populations (P < 0.01). To estimate the phylogeographic history of cultivated Asian pears, the TCS-derived network of haplotypes detected in all of the cultivated pears was constructed ( Figure 3B). The results showed that two of the most frequent haplotypes (H14 and H6) corresponded to two major haplotype lineages in the TCA network and that haplotype H13 occupied a central position. The common haplotypes H7 and H3 were closely associated to haplotype H6, and haplotype H14 was surrounded by the rare haplotypes H4, H5, H9, H15, and H17. As shown in Figure 3A, haplotype H14 was very prevalent in southwestern populations (i.e., SC, YN, and GZ), and its frequency was decreased in southern, eastern and northern areas, while H6 was the most frequent haplotype in mid-southern (i.e., HHJ), eastern (i.e., ZJ) and northern populations (i.e., JL), as well as JP populations (i.e., JPa and JPb). Haplotypes H4, H5, H9, H15, and H17 located near H14 in the TCS network were largely confined to southern populations; haplotypes H3, H7, H10, and H11 (located in the same haplotype lineage of H6 in the TCA network) coexisted with H6, occurring across the eastern and northern areas ( Figure 3A). Among these haplotype H3 was largely detected in the northwestern population GQ, and haplotype H7 was primarily detected in southern populations (i.e., FJ and GG). The central haplotype H13 was only detected in a few cultivars from populations GQ, FJ, and GG. Haplotypes H1 and H19 were only detected in northern geographic populations (i.e., GQ) and cultivated UPs. A Neighbor-Net splits graph of haplotypes was constructed for East Asian cultivated pears and related wild species (Supplementary Figure S3). Many parallel edges were observed in the splits of the NN graph, indicating homoplasy among the combined haplotypes. Haplotype H13, in the inner network, was linked to the outside group (cpDNA haplotypes of accD-psaI and trnL-trnF identified in Malus species), and most of the haplotypes were detected in pear accessions, which was consistent with the TCS network ( Figure 3B). Haplotypes H1 and H19 were linked with four unique haplotypes of WUP, forming a separate lineage, and haplotypes H4, H5, H9, H14, H15, and H17 displayed the closest relationships to the two unique haplotypes of P. pashia.

Putative Dissemination Route of Cultivated P. pyrifolia Inferred From cpDNA Data
Populations SC and YN, GQ, GZ, HHJ, AJ, ZJ, JPa and JPb, FJ and GG, SS, and SHH and JL were treated as Pop 1-10, respectively. Populations SC and YN, populations GG and FJ, populations SHH and JL, as well as populations JPa and JPb, were combined separately, since they were not only geographically adjacent and dominated by the same cpDNA haplotype and had the similar composition of gene pool inferred from STRUCTURE analysis based on nSSR data. Three putative historical models of Asian cultivated P. pyrifolia were constructed according to cpDNA haplotype variations among the geographic populations ( Figure 2B). As a result, the estimates indicated that the first scenario was more likely to be the evolutionary route of cultivated P. pyrifolia in East Asia (Supplementary Figure S4

High Levels of Genetic Diversity in Asian Pear Cultivars Across the Whole of China and Japan
In this study, more than 400 local pear cultivars were collected across the main pear distribution areas of China and Japan (Supplementary Table S1), and used in the genetic analysis. The sample size was much larger than any previous related studies (e.g., Bao et al., 2007Bao et al., , 2008Yao et al., 2010;Cao et al., 2012;Iketani et al., 2012;Song et al., 2014;Jiang et al., 2015;Nishio et al., 2016). Twenty-five genome-covered nSSR loci in this study revealed a high level of genetic diversity (H E = 0.744) in Asian pears, which was similar to that (0.73) reported by Liu et al. (2015), and much higher than that (0.64) reported by Song et al. (2014). Genetic analyses of geographic populations revealed that the levels of genetic diversity in southern populations (south of the Yangtze River) were generally greater than those in northern populations (north of the Yangtze River) ( Table 1). Therefore, the conservation priority should be given to local cultivars from southern areas. The cultivars from geographic populations with high levels of genetic diversity are also useful genetic resources for future pear breeding programs (Kumar et al., 2017). A decrease in genetic diversity was also identified in two geographic populations in Japan, in agreement with previous studies (Bao et al., 2007;Iketani et al., 2012). In the future, cultivars with different genetic background, especially from South China should be applied in Japanese pear breeding programs to broaden the genetic diversity of Japanese pears as suggested by Nishio et al. (2016). Additionally, two cpDNA intergenic fragments displayed high degrees of genetic diversity (Hd = 0.792) among Asian pear accessions, which was similar to that in the wild pear species P. betulifolia (Hd = 0.807)  and higher than that in wild P. pashia (Hd = 0.718) (Liu et al., 2013). The generally high levels of genetic diversity observed in Asian pear cultivars may be attributed to the large number of local cultivars and effective nSSR markers used in our studies compared with previous studies.
Wild P. pyrifolia species with large fruit, having five carpels and russet, yellow, or smooth green skin, were supposed to have originated in the Yangtze River Valley and the adjoining southern area (Pu and Wang, 1963). However, no populations of wild P. pyrifolia still exist in either China or Japan, possibly resulting from habitat destruction and land overexploitation (Heywood et al., 2007). Therefore, we cannot investigate the origins of cultivated P. pyrifolia by comparing P. pyrifolia cultivars with their wild ancestors. Instead, we collected as many local pear cultivars across the major areas of China and Japan as possible to cover all genetic variations. In addition, most haplotypes of wild oriental Pyrus species identified by Zheng et al. (2014) were detected in the cultivated pears, including all seven tH (haplotypes of trnL-trnF) and 11 aH (haplotypes of accD-psaI), as well as two newly identified tH haplotypes [tH20 (MH010656) and tH21(MH010657)] and three newly identified aH haplotypes [one in cultivated Asian pears firstly reported in this study [aH36 (MG751777)] and two in WUP accessions [aH37 (MH010658) and aH38 (MH010659)]. The high levels of genetic diversity in cultivated P. pyrifolia based on nSSRs indicated that the enriched level of genetic diversification of cultivated Asian pears was represented by the Asian pear cultivars used in this study.

Poor Differentiation Among Cultivated Asian Pear Groups
Cultivated pears worldwide are traditionally divided into two large geographic types, soft-fleshed European pears (P. communis L.) and the crisp-fleshed Asian pears (Janick, 2000), with a high-level genetic differentiation between them (Kumar et al., 2017). Asian pears contains two or three species that are generally classified into four major cultivar groups, including UP (P. ussuriensis Maxim.), CWP (P. pyrifolia CWP group, also traditionally assigning to P. x bretschneideri), CSP (P. pyrifolia Nakai), and JP (P. pyrifolia Nakai) according to their geographic distribution, and morphological and physiological traits (Teng and Tanabe, 2004).
Chinese white pear has been misassigned to P. x bretschneideri Rehder for a long time (Teng and Tanabe, 2004). Many researchers did not accept this classification because of the morphological differences between CWP and P. x bretschneideri (Kikuchi, 1946;Teng et al., 2002), and the many more similar morphological characteristics, such as leaf morphology and fruit texture (Pu and Wang, 1963), as well as peroxidase isozymic patterns (Lin and Shen, 1983), shared between CWP and CSP. Our previous studies (Teng et al., 2002;Bao et al., 2007;Jiang et al., 2015Jiang et al., , 2016 and other researchers' studies (Iketani et al., 2012;Song et al., 2014;Liu et al., 2015) using different DNA markers suggested that CWP was also genetically close to CSP. Meanwhile, JP cultivars may share the same germplasm as CSP based on similar morphological traits, especially fruit traits (Teng and Tanabe, 2004). Some researchers believed that these cultivars were domesticated from wild P. pyrifolia in Japan (Kikuchi, 1948), which seems to be supported by a few recent studies on the genetic relationships among Asian pears in which the majority of JPs did not cluster together with CSPs in dendrograms (Kimura et al., 2002;Nishio et al., 2016). However, other researchers suggested that JP cultivars might have been introduced from ancient China (Kajiura and Suzuki, 1980;Teng et al., 2002) because a flourishing sea route existed between Japan and the coastal areas of China. Thus, the so-called wild P. pyrifolia individuals in Japan were supposed to be escapees (Teng and Tanabe, 2004). The phenomenon of escapees is complicated in fruit trees owing to the existence of semi-wild or wild populations in nearby cultivated areas and to difficulties in distinguishing true wild types, escapees, and hybrids (Eliezer, 2013). This could be further complicated by gene flow and introgression from cultivated forms into wild relatives, as discovered in grape (Di Vecchi-Staraz et al., 2008), cultivated almond (Mohamed, 2004), and cultivated olive (Lavee and Zohary, 2011). Because of the genetic affinities of JP to pear cultivars from ZJ, based on multiple DNA markers, researchers proposed that this province might be the origin of JP cultivars (Teng et al., 2002;Bao et al., 2007;Jiang et al., 2016).
In the present study, phylogeographic analyses of cpDNA haplotypes suggested that CSP, CWP, and JP shared common cpDNA haplotypes (H14, H6, and H3) and that these haplotypes were distributed throughout geographic populations ( Figure 3A). Even though the two major haplotype lineages (H14 and H6) displayed differential geographic distributions, they accounted for large proportions of all of the cultivar groups. Clusters N-8.1 and N-8.2 contained most of pear cultivars from CSP, CWP, and JP based on the Bayesian clustering of chloroplast DNA (Figure 2B). The results of STRUCTURE clustering inferred from nSSR data indicated that the CWP genotypes shared the major gene pools with CSP genotypes (Figure 1A), AMOVA analyses showed that there were very low level of genetic variations partitioned among cultivated groups of P.pyrifolia (2.49 and 3.61% for nSSR and cpDNA, respectively), which was supported by Iketani et al. (2012), inferring that they had a close genetic relationship, and supporting our previous hypothesis that CWP from northern China might be an ecotype of CSP (Bao et al., 2007). Nevertheless, STRUCTURE clustering clearly showed that the genetic structure of JP cultivars was distinctively different from those of CSP and CWP cultivars in China (Figure 1A), which is consistent with the previous studies where most JPs could not be clustered together with CSP cultivars in the dendrograms (Bao et al., 2007;Iketani et al., 2012;Nishio et al., 2016). The distinctiveness between the JP cultivar group and the CSP/CWP cultivar group, based on nuclear markers, might be explained by the limited primitive genetic backgrounds in the development of JP (Iketani et al., 2012), which was also supported by JP populations having the highest percentages of pure/non-admixed genotypes ( Figure 1B). In addition, the dominant "green" gene pool in JPs was shared with some cultivars from ZJ and FJ in China ( Figure 1B). Higher similarities were detected in maternally inherited cpDNA haplotypes of pear cultivars from both Japan and the coastal areas of China, which indicated that JP seeds could have been introduced from China through human activity. Thus, it was inferred that CSP and JP could be derived from the same progenitor of P. pyrifolia in China.
Because UP cultivars have different morphological traits than the cultivated P. pyrifolia, they have generally been considered as directly derived from the WUP of northeastern China (Pu and Wang, 1963;Teng and Tanabe, 2004), although a few studies revealed levels of genetic differentiation between UP cultivars and WUP (Wuyun et al., 2013;Nishio et al., 2016). In this study, 22 UP cultivars and 18 WUP were collected mainly from northeastern China for the genetic analyses. The cpDNA analysis determined that UPs share the common cpDNA haplotypes (H14 and H6) with cultivated P. pyrifolia: CSP, CWP, and JP ( Figure 3A), and most of UPs and cultivars of CSP and CWP, especially those from southwest populations, were grouped in cluster N-8.1 in the Bayesian clustering of chloroplast DNA ( Figure 2B). The dominant haplotypes of WUP were not detected in UP cultivar accessions, and two common haplotypes of UP did not occur in WUP accessions. Although both UP and WUP accessions were dominated by the same gene pool, as inferred from nSSR markers (Figure 1A), a high proportion of admixed genotypes were detected in UP accessions. Additionally, the two major gene pools (Figure 1A) of CSP and CWP also existed in UP accessions, while all of the WUP individuals were dominated by a single gene pool. And the accessions of WUP were mostly included in cluster N-8.8, separated from most of the cultivated Asian pears (Figure 2B). Thus, UP cultivars could not be directly domesticated from the WUP, and P. pyrifolia cultivars might be involved in their origin. This supports the findings of Jiang et al. (2016) and Yu et al. (2016) in which retrotransposon markers were used to determine that the major gene pools in UP cultivars came from cultivated P. pyrifolia. However, commercial UP cultivars generally have small fruit with poor quality (Kikuchi, 1946). In UP breeding programs, local UP cultivars should be further crossed with cultivated P. pyrifolia to gain a better fruit quality and maintain their cold tolerance. In addition, a higher level of genetic differentiation between UP and WUP was detected in cpDNA (F ST = 0.5739) than in nSSRs (F ST = 0.1031). By comparing Bayesian clustering based on cpDNA variations and nSSR STRUCTURE of the Asian pear accessions, we hypothesized that cultivars of P. pyrifolia might be the maternal parents of UP cultivars and that WUP might be the paternal parents. Thus, the UP cultivars included in this study might be derived from the interspecific hybridization of cultivated P. pyrifolia with local WUP. In addition, based on the sequence divergences of low-copy nuclear genes (such as LFY2) among Pyrus species, Zheng et al. (2014) proposed that occidental species might be the paternal parents of some UPs (i.e., 'Ruanerli'). However, fewer UP cultivars were used in this study. Thus, the origin of UPs in East Asia might be more complex. In the future, different nuclear markers and more UP accessions of different geographic origin will be required to confirm their origins.

Putative Origin Center and Dissemination Routes of Cultivated P. pyrifolia
Compared with annual crops, fewer studies have provided analyses of the origin and domestication of fruit trees (Eliezer, 2013). One of the best known is the study of Zohary and Spiegel-Roy (1975) in which they clarified the evolutionary history of four classical Old World fruits (olive, grape, date, and fig) using "fossil" evidence and morphological clues from living plants and their wild relatives. Recent studies on the origins of apple (M. domestica), based on molecular markers, identified the principal progenitor (M. sieversii) and secondary contributor (M. sylvestris) of modern apples by exploring genetic variations between the domesticated apple and its wild relatives (Velasco et al., 2010;Cornille et al., 2012). Even though the genome of CWP was published by Wu et al. (2013), the origin and evolutionary history of pear cultivars native to East Asia are still obscure and have not been well illuminated. This is mainly because of unavailable wild populations, especially for cultivated P. pyrifolia, which is considered the most important germplasm of Asian pear cultivars owing to its higher commercial value and wide cultivation. For fruit trees, the invention of vegetative propagation based on grafting or cutting revolutionized their domestication processes. These practices were necessary to fix the selected traits discovered in the wild and may have shortened the evolutionary distance between modern fruit cultivars and their wild progenitors (Zohary and Spiegel-Roy, 1975;Eliezer, 2013). The cultivation of Asian pears has been recorded from antiquity (Pu and Wang, 1963), when grafting had already become a well-established technique (Mudge et al., 2009). Thus, although no wild P. pyrifolia individuals were included in this study, the sampled local cultivars of different geographic origin could provide significant clues to the evolutionary history of cultivated Asian pears.
In this study, 401 cultivated P. pyrifolia, including seven geographic populations of CSP from southern China, five geographic populations of CWP from northern China, and two geographic populations of JP from across Japan were used in the genetic analyses (Supplementary Table S1). A greater level of genetic diversity was detected in the geographic populations of CSP from southern China compared with geographic populations of CWP and JP from northern China and Japan based on both nSSRs and cpDNA haplotypes (Tables 1, 2). This suggested that CSP accessions possessed the highest polymorphism level and might have contributed to the development of the CWP and JP. Furthermore, the high level of genetic diversity and the greatest number of cpDNA haplotypes were detected in the southwestern SC and YN population of CSP, suggesting that the southwest might be the origin center of cultivated P. pyrifolia. In contrast, a relative lower genetic diversity level and fewer cpDNA haplotypes were discovered in the two geographic populations of JP, indicating that they might be peripheral populations of cultivated P. pyrifolia.
In contrast to nuclear DNA, which provided information on the current genetic structures, cpDNA, with a lower variation rate, could reveal the evolutionary histories of plants.
The putative dissemination routes of cultivated P. pyrifolia were estimated from cpDNA haplotypes using approximate Bayesian computation (Figure 4), which suggested that cultivated P. pyrifolia might have been initially domesticated from wild P. pyrifolia in the southwest origin center. Then, cultivated pears disseminated across East Asia with the migration of people. Some of the cultivars spread from the original center to its eastern neighbor (i.e., GZ), while others spread to the Middle Yangtze River Valley (i.e., HHJ), then further to southern (ie., FJ and GG) and eastern China (ie., AJ and ZJ), as well as northern regions (i.e., SHH, Shanxi, and Liaoning Provinces). This indicated that cultivars from the Middle Yangtze River Valley areas might play pivotal roles in the southward, eastward and northward spread of cultivated P. pyrifolia, possibly resulting from the long historical commercial and cultural exchanges between Middle-Lower Yangtze Valley and North China. The composition of cpDNA haplotypes in geographical populations has been changed during the dissemination of cultivated P. pyrifolia possibly through conscious or unconscious selection, which may be one of the reasons that an increase of genetic differentiation based on cpDNA haplotype frequency was detected between populations from southwestern and populations from other regions (Supplementary Table S9). However, some rare haplotypes identified in few populations were not detected in cultivars from the putative origin center, such as haplotype H13, identified as the dominant and ancestral haplotype of oriental pears (Zheng et al., 2014), was only detected in a few cultivars scattered in northwestern (GQ) and southern (FJ and GG) populations that had higher numbers of unique nSSR alleles ( Table 1). Cultivars containing these haplotypes might be involved with the interspecific hybridization between cultivated P. pyrifolia and local wild Pyrus species, rather than pollination of cultivated pears by local wild species, since chloroplast DNA was maternally inherited in Rosaceae and could be transmitted only by seeds.
Moreover, geographical population GQ were dominated by the two gene pools that had a relatively high number of unique FIGURE 4 | Three putative dissemination routes of cultivated P. pyrifolia in East Asia. Scenario 1 assumes that cultivated P. pyrifolia from southwestern regions (SC and YN) and Middle-Yangtze valley areas (HHJ) were diverged firstly, and then further spread from HHJ to southern (FJ and GG), eastern (AJ and ZJ) and northern areas (SS, GQ, SHH, and JL); Scenario 2 assumes that cultivated P. pyrifolia in northwestern (GQ) and southern areas (FJ and GG) were directly disseminated from southwestern populations (SC and YN); Scenario 3 assumes that cultivated P. pyrifolia from northwestern population (GQ) were initially disseminated from southwestern populations (SC and YN), from which cultivars were further spread to northern areas (SS, SHH and JL).
alleles (8), might be explained by separate regional cultivar development. By analyzing the putative dissemination route of cultivated P. pyrifolia, the primitive cultivated pears from population GQ might have spread from the neighboring Shanxi Province, as described by Chang et al. (2017), along the ancient Silk Road (Teng et al., 2001). Comparatively, cultivars preserved in the mountain and hill areas of southwestern China might be more ancient types and should have the highest conservation priority.

Comparison of Nuclear DNA-and cpDNA-Based Genetic Structures
Vegetative propagation based on grafting or cutting has undoubtedly revolutionized the domestication and breeding process of perennial trees, allowing the maintenance and spread of desirable lines despite self-incompatibility and a long juvenile phase (Janick, 2005). Because of the spread of grafting, "instant domestication" has been considered in the domestication of apple (Robinson et al., 2001). Under such circumstances, apple individuals with different geographic origins might be grafted and chosen from local wild types. Thus, the parental contributions in the origin of modern apple cultivars could be quite diverse (Robinson et al., 2001). Nuclear DNA is inherited parentally, and chloroplast DNA is inherited maternally in Rosaceae (Hu et al., 2008). The congruency of the nuclear DNA-and cpDNAbased genetic structures in the geographic populations might be attributed to the extensive practice of clonal propagation during the domestication of fruit trees, but the incongruence between them would indicate an intensive influence of cross-pollination (Nikiforova et al., 2013). The combination of nuclear SSR and chloroplast DNA has provided new insight into the origin of European pear cultivars (Ferradini et al., 2017). Therefore, the genetic comparison of parental inherited nuclear DNA markers and maternal inherited cpDNA in this study might generate meaningful results for the domestication of cultivated Asian pears.
In the present study, according to the gene pools of the current pear cultivars detected from nSSRs using Bayesian statistical inferences, three genetic groups were revealed for cultivated P. pyrifolia: (a) JPs; (b) cultivars from South-Central China northward to northeastern China, covering the main pear production area in China; (c) cultivars from southwestern China to southeastern China including Yunnan, Guizhou, Guangdong, Guangxi, and Fujian Provinces (Figure 1B), of which cultivars from South-Central China northward to northeastern China were also largely grouped in one cluster (N-8.2) based on Bayesian clustering results of cpDNA data and were dominated by the same cpDNA haplotype (H6) (Figures 2B, 3A). The congruency of the nuclear DNA-and cpDNA-based genetic structures in the geographic populations from these areas might be attributed to the extensive practice of clonal propagation during the spread of cultivated pears form Middle-Lower Yangtze River areas to the North China Plain. Low levels of genetic differentiation were also detected among populations within these areas based on both nSSRs and cpDNA (Supplementary Table S9). Based on the genetic groups of cultivated P. pyrifolia described above, the following three major groups for cultivated P. pyrifolia were suggested: Japanese, Central South-North China and Southwest to Southeast China.
More genetic clusters were revealed by Bayesian clustering of chloroplast DNA by comparison to genetic groups identified from nSSRs using Bayesian statistical inferences (Figures 1A,  2A). This is possibly due to some rare cpDNA haplotypes specific to certain geographic populations, since most of the cultivated pears were grouped in the first four clusters (Supplementary  Table S5). However, most of cultivars from southernmost populations (i.e., GG and FJ) shared the major gene pools identified from nSSRs with cultivars from southwestern populations (i.e., SC, YN, and GZ), but they were dominated by different haplotypes (Figure 3B) and separately grouped in different clusters base on the Bayesian clustering of cpDNA (Figure 2A). A higher level of genetic differentiation was also detected from cpDNA than that from nSSRs between southwestern populations and southernmost populations. Because pear is a temperate wood species and needs low temperature exposure to break bud dormancy , cultivars from southern China (i.e., GG and FJ) may have become tolerant to warm temperatures in winter through conscious or unconscious selection (Eliezer, 2013), and most of cultivars carrying that haplotype were maintained during the selections. The differentiated genetic structures among geographic populations from different regions used in our study might reflect the synergistic effects of ecogeographical factors and human selection during cultivar improvement and spread. Diversity may have been restored to cultivars through gene introgression from wild types used as stocks when the cultivated clones were grown in new environments (Zohary and Spiegel-Roy, 1975).
Gene introgression from wild relatives has been detected in cultivated pears by some studies, indicating that it may contribute to the genetic diversification of current pear cultivars of different geographic origins (Zheng et al., 2014;Jiang et al., 2016). The introgression could be facilitated by self-incompatibility and cultural practices, such as selection from open-pollinated seeds (Nikiforova et al., 2013). For example, P. pashia may be involved in the development of cultivated P. pyrifolia in southwestern China (Zheng et al., 2014;Jiang et al., 2016) and wild P. ussuriensis may be involved in the evolution of CWP in northern China . Although wild Pyrus species bear fruit with poor quality, they could contribute to disease resistance and to the adaptability of cultivated pears. High proportions of admixture genotypes based on paternal nSSRs detected in geographical populations SC, HHJ and GQ might be due to intensive introgression from local wild relatives. Therefore, samples of local wild species should be collected and genotyped for comparison with local cultivars to determine if wild relatives were involved in the development of local pear cultivars. In summary, a high level of genetic diversity was detected in Asian pear cultivars, suggesting an increase in the level of genetic diversification that was represented by the Asian pears used in this study. AMOVA analyses based on both cpDNA fragments and genome-covering nSSR loci demonstrated close relationships among the four major pear cultivar groups according to relatively lower genetic variations partitioned among groups, but a higher level of genetic differentiation was revealed between UP and WUP. Phylogeographical analyses of cpDNA haplotypes indicated that CSP, CWP, and JP could be derived from the same progenitor of P. pyrifolia in China, and cultivars of P. pyrifolia might be the maternal parent of cultivated UPs. Southwestern populations including SC and YN, having a greater genotypic diversity, might be the origin center of cultivated P. pyrifolia. The putative dissemination route of cultivated P. pyrifolia indicated that cultivars from the Middle Yangtze River Valley might play important roles in the southward, eastward, and northward spread of cultivated P. pyrifolia. Based on the structural clustering of nSSRs, cultivated P. pyrifolia could be divided into the following major groups: Japanese, Central South-North China, and Southwest to Southeast China.