Ninety-five accessions from European and Asiatic Pyrus spp. were included in this study (Table 1 and Table S1). Many of the 95 accessions were well documented by historically reliable sources. For others, lacking some information, the original name, like that given by the donor, was maintained. In all other cases they were named Unknown. Therefore, based on the initial information, the 95 accessions used in this study were grouped into reference species (RS), commercial varieties (CV), local varieties (LV) and unknown accessions (UA), where CV and RS were used as controls. Moreover, 75 accessions were divided into three geographic groups, following the classification of Fideghelli (2007): Mediterranean (2), East-Asian (4), and European (69 in total: 24 CV, 2 RS, 43 LV). Because the remaining 20 accessions could not be assigned to any geographic group, they were grouped as Unknown (Table S1).

Total genomic DNA was isolated from young leaves using the DNeasy Plant Kit (Qiagen) following the protocol provided by the manufacturer. Nine pear SSR primer combinations (Yamamoto et al., 2002a,b) were used (Table S2). PCRs were carried out with the Type-it Microsatellite PCR Kit (Qiagen) containing 1X Type-it master mix with 0.2 μM of each forward and reverse primer and 20 ng of DNA and H₂O to a final volume of 20 μl. Amplification was performed as follow: an initial step at 95°C for 5 min followed by 30 cycles at 95°C for 30 s, 54–62°C for 30 s, and 72°C for 30 s, and a final extension at 72°C for 10 min.

PCR products were separated and analyzed on a 3130 XL DNA Analyzer (Applied Biosystems). The size of the amplified products was determined on internal standard DNA (GeneScan 500 Liz, Thermo Fischer Scientific) and the scorable peaks were assigned by GeneMapper software (Applied Biosystems).

The hypervariable accD-psaI intergenic spacer was amplified by PCR using accD F2 (Zheng et al., 2014) and psaI 75R (Small et al., 1998) primers. PCR reactions contained: 1X Phusion HF Buffer, 200 μM dNTPs each, 0.5 μM each primer, 0.2 U Phusion Taq DNA polymerase (Thermo Fischer Scientific), 30 ng genomic DNA and H₂O to a final volume of 50 μl. PCR amplification was carried out by GeneAmp PCR system 9700 (Applied Biosystems) programmed as follow: 98°C for 30 s, followed by 30 cycles of 98°C for 10 s, 66°C for 15 s, 72°C for 30 s, and then 72°C for 10 min.

The purified PCR products of about 1,000 bp were sequenced at the Polo di Innovazione Genomica, Genetica e Biologia (Perugia, Italy). The new cpDNA sequences were recorded in GenBank with accession numbers from KY606436-KY606530.

The statistical analysis of the SSR data matrix included (i) the estimation of observed (Ho) and expected (He) heterozygosity (Nei, 1978), (ii) the F statistics (F_is and F_st) (Weir and Cockerham, 1984), and (iii) the analysis of molecular variance (AMOVA) by estimating the fraction of the genetic variation among and within populations (Excoffier et al., 1992; Michalakis and Excoffier, 1996). The software packages GENODIVE (Meirmans and Van Tienderen, 2004) and SPAGeDi1.2 (Hardy and Vekemans, 2002) were used for these purposes, being able to analyze data files containing diploid and triploid accessions.

SSR data were also converted to a binary data matrix by assigning “0” to the absence of a defined allele and “1” to its presence, and were used to estimate a similarity matrix using the coefficient of Dice (Dice, 1945), and the individuals were clustered by the unweighted pair group method with arithmetic mean (UPGMA) and validated by 1,000 bootstrap replicates using PAST software (Hammer et al., 2001).

The SSR profiles of the 95 accessions were used to investigate the population structure through the Bayesian model-based clustering procedure of STRUCTURE ver. 2.2.3 (Pritchard et al., 2000). The analyses were based on an admixture ancestral model with correlated allele frequencies, and the number of K clusters was determined by simulating a range of K-values starting from one to ten. A burn-in and a run length of the Monte Carlo Markov Chain (MCMC) of 200,000 and 500,000 iterations for data collection with 10 runs per K-value were used. The best K-value was determined through the ΔK method (Evanno et al., 2005) by using Structure Harvester ver. 0.6.193 application (Earl and vonHoldt, 2012). A second step analysis of STRUCTURE was then performed separately at these K-values with 600,000 burning period and 1,000,000 MCMC repeats after burning. The 95 individuals were assigned to the groups according to their highest membership coefficient, considering a strong affinity when the assigning probability (qI) was ≥0.80 (Breton et al., 2008; Pereira-Lorenzo et al., 2008; Miranda et al., 2010; Urrestarazu et al., 2012).

In order to analyze the cpDNA sequences, forward and reverse sequences from each sample were assembled and aligned using Sequencher™ 5.10 (Gene Codes Corporation). Accession gr4_077 P. communis cultivar William was selected as a reference sequence (GenBank accession number KY606501). It is one of the most common and suitable exemplars of European pear cultivars.

Different cpDNA sequence variation parameters were estimated by using DnaSP 5.1 software (Librado and Rozas, 2009). AMOVA was carried out as reported above for the SSR data. The 95 sequences were also compared to all accD-psaI data recorded in GenBank: 203 Pyrus sequences from 26 species and 4 putative inter-specific hybrids. The final dataset included 298 Pyrus accessions from three major geographic areas: Europe, the Mediterranean Area and Asia; the latter included Eastern and Middle Eastern subgroups (Tables S3, S4). A sequence from Malus domestica was included and used as an outgroup (Katayama et al., 2012). Indels of different lengths at the same position were separately treated and coded with the number of inserted or deleted nucleotides (Table S5). The final alignment was compared to the reference sequence, thus allowing sequence classification in different haplotypes. The evolutionary relationships among haplotypes were visualized through the construction of a Median-Joining network using the software Network 4.6 (http://www.fluxus-engineering.com/), each indel was considered as a single mutational event and partitioned by using different indel codes (Table S6).

As for principal component analysis (PCA), performed by XLSTAT (2011) statistical software, each species was considered as a discrete variable, the initial dataset was converted into principal components (PCs) and it was possible to graphically display the relationships among the intergenic regions accD-psaI of all sequences.

Scorable amplicons were produced for all nine nuclear SSRs, with a total of 179 alleles. The average number of alleles per locus was 20, ranging from five (NH030a) to 29 (NH023a), but the number of effective alleles per locus was significantly lower (NAe = 6.8). At locus NH030a, out of five alleles, allele 167 showed a frequency of 0.97; thus loci polymorphism was P = 0.89 (Cavalli-Sforza and Bodmer, 1971).

As many as 24 individuals out of 95 (23%) showed only one locus with a third allele. Ten genotypes showed 2 loci with a third allele (9.5%), while nine individuals showed a third allele at more than 2 loci (Table S7). In particular, except for locus NH030a, all other loci showed at least one genotype with three alleles. NH026a and NH023a identified 16 and 18 individuals with three alleles, respectively; NH029a had only two, while the other loci identified between seven and ten individuals with three alleles.

The AMOVA carried out with all loci showed that most of the existing variability was within groups (97%) rather than among groups (Table S8), which is in agreement with many outbreeding species and similar studies in pear (Jiang et al., 2009; Miranda et al., 2010; Wolko et al., 2015; Wuyun et al., 2015). Nonetheless, differences among the 4 groups (CV, LV, RS, and UA) were found in terms of the effective number of alleles and levels of heterozygosity.

The 95 individuals in the present study were not a panmictic population, as classically defined in population genetics. Therefore, the indices used here (number of effective alleles, observed and expected heterozygosity, F_IS and F_ST) as estimates of existing genetic variability should be considered with caution. It is worth noting that the effective number of alleles (NAe) is a measure of the genetic variability (Zouros, 1979), also valid in the presence of small samples (Nielsen et al., 2003). In our data the highest average number of effective alleles was found in RS (NAe = 12.4), while the lowest was in CV (4.8) (Table 2); it is interesting to note that the values of LV and UA were intermediate and similar to one another (6.9).

A similar trend was observed in terms of heterozygosity. The average expected heterozygosity, considering all of the pear accessions, was He = 0.74, ranging from 0.87 in RS to 0.70 in CV. The highest mean expected value of heterozygosity for all accessions was found at loci Nb109a and NH027a (He = 0.91), while the lowest was at locus NH030a (He = 0.063). Noticeably, at this locus, He in the reference species was 0.61, while in all other groups (LV, CV, and UA) He = 0, monomorphic and due to the fixation of allele 167. In addition, in RS, NH026a, NH027a, and Nb103a loci showed He-values as high as 0.96 (Table 2).

The mean F_ST-value equal to 0.014 (Table 3) indicates a moderate differentiation among the four groups (P < 0.0005), thus confirming that most of the variation is within groups; at locus NH030a the F_ST value of 0.278 indicates a significantly high genetic differentiation (P < 0.0002), suggesting that at this locus the selection for local adaptation could have been so strong as to restrict its variation to a single allele in all cultivated forms.

The nine SSRs were also used to determine the genetic structure among the 95 accessions of Pyrus spp. The plot of the average log-likelihood values for Ks ranging from 1 to 10 and the distribution of ΔK-values (Evanno et al., 2005) according to K-values is shown in Figure S1. Two peaks were found, corresponding to K = 3 and K = 6, and the hierarchical genetic structure was investigated at K = 3. A threshold value qI ≥ 0.80 was used to assign individuals to the clusters (Figure 1). Structure at K = 3 was able to define 3 clusters (from 1A to 3A). Cluster 1A included mostly commercial cultivars and P. syriaca, one of the reference species able to intercross with P. communis. Cluster 2A included 27 accessions, almost all the rest of the reference species (six, except P. cossonii), most of LV provided by the “Archeologia Arborea” collection and only one commercial cultivar (Carmen). Cluster 3A included 22 accessions, mostly LV (12) and 7 CV, four of which were of French origin. Some accessions (11 LV, 6 CV, 5 UA, and P. cossonii) were clustered in the admixture group. Although this was not the object of the present study, as physiological and phenotypic information was scarce and incomplete, Cluster 1A included the majority of accessions characterized by summer ripening, Cluster 2A all accessions with an autumn-winter ripening period, while accessions in Cluster 3A were not clear cut in terms of ripening period.

The structure at K = 6 allowed six Clusters (from 1B to 6B, Figure 1) to be distinguished. Cluster 1B included some unknown accessions and some local varieties provided by “Archeologia Arborea” (gr5). Cluster 2B included two commercial cultivars. Pera Briaca, Pera di Montelupone, Pera Monteleone and Pera Ruzza were assigned to Cluster 3B, characterized by winter ripening and a round fruit shape (Table S1). Cluster 4B included only commercial cultivars. Cluster 5B was represented by two local varieties (Pera Broccolina and Pera Volpina) and some cultivars, all characterized by autumn-winter ripening. On the contrary, Cluster 3B included accessions characterized by winter ripening and a round fruit shape (Supplementary Table S1). Cluster 6B included four local varieties and five reference species (P. betulifolia, P. calleryana, P. caucasica, P. pyrifolia, and P. ussuriensis).

The accessions of Pyrus spp. were clustered by UPGMA following the similarity estimates of Dice's coefficient (Figure 2). Reference species were found in a different subcluster, quite distant from most accessions, with the exception of P. cossonii, P. pyraster, and P. syriaca. In particular, P. syriaca is clustered with Coscia, Coscia Tardiva, Decana del Comizio, Santa Maria Morettini and William, while P. cossonii and P. pyraster with some local varieties, such as Pera di Tiberio and Pera Rubbia.

Several entries clustered at a similarity coefficient equal to 1. For example, the two accessions of Pera Agostina (gr1_007 and gr1_055) clustered with Coscia (gr1_020), Pera della Battitura (gr1_056) and one unknown accession (gr1_024), suggesting cases of synonymies with Coscia. As a matter of fact, “Agostina” in Italian means August, and “battitura” means “threshing,” most likely referring to the harvesting of cereals, normally occurring in summer, which is the ripening period of Coscia.

Similarly, Pera Grassana (gr1_013), Passa Crassana (gr4_073) and the unknown accession gr1_002 are all the same and likely to be Passa Crassana. Moreover, the present study confirmed that Spadona d'Inverno (gr3_050) and Curato (gr4_081) are synonyms and that gr1_022 (Unknown), clustering at a similarity of 1, is also likely to be the same accession. The same holds for Pera Tonda Roggia (gr3_047) and Pera Vernia (gr3_048), probably synonyms for the same genotype.

Furthermore, entry gr1_010 (Unknown) could be named Scipiona (gr2_036) and entry gr1_011 could be named Pera Ruzza (gr1_014). Entries gr1_016 and gr1_043, known as Pera Monteleone, are confirmed to have the same genotype and clustered together with Pera di Montelupone.

In addition to the clusters that joined at a similarity equal to 1, the analysis confirmed that the two accessions of Angelica clustered together, as well as the two accessions of Bergamotte Esperen and all the accessions of Pera Lardaia and Pera Volpina. The analysis summarized in Table S9 was also able to give insights for other accessions, such as Pera Grossa d'Autunno, similar to Curato and Spadona d'Inverno, and some unknown accessions that could be named as the closest commercial or local variety.

The 95 novel chloroplast DNA sequences ranged from 622 to 899 bps depending on the presence of ten indels (four deletions at nps 210, 322, 545, and 601, and six insertions at nps 305, 544, 600, 620, 636, and 665); gr4_088 was the only accession without indels and differed from the reference sequence William only by the transversion at np 567 (567T). After grouping the accessions into four categories (reference species, RS; commercial varieties, CV; local varieties, LV; and unknown accessions, UA), AMOVA was carried out on all haplotypes and the results obtained through the SSR were confirmed. In fact, most of the observed variance was attributable to differences among samples within groups (96.69%), rather than the variability among groups (3.31%).

The analyses was then extended to the entire dataset of 298 cpDNA sequences (95 from the present study and 203 from GenBank). A total of 75 haplotypes were identified and named from HT01 to HT75 (Table S3). The 95 accessions of the present study were represented by 14 haplotypes, the most frequent of which was haplotype HT06 (47%). Out of these 14 haplotypes, 8 were shared among different accessions of P. communis (HT02, HT06, HT07, HT08, HT09, HT11, HT12, HT13), P. caucasica (HT09), P. syriaca (HT06) and some unknown accessions (HT02, HT06, HT08, HT09, HT12, HT13), while 6 were novel (Table S10): one was found in P. ussuriensis (HT18), one in P. cossonii (HT19) and the other four (HT08, HT11, HT12, and HT13) were found in P. communis. As evident in Figure S2, the new haplotypes (circled in red), whose sequences have never been reported before, were present in about 28% of the whole P. communis species and this value rises up to 50% if the most common haplotype HT06 is excluded from the analysis. This means that half of the P. communis samples have sequences never reported before. It is worth noting that four novel haplotypes (HT08, HT11, HT12, and HT13) were well represented in many accessions from Central Italy, most of which were provided by “Archeologia Arborea” (13 out of the 18 accessions, Table S1).

The reconstructed network of the cpDNA intergenic region accD-psaI clearly defined the distribution of the 95 accessions in different branches. All sequences clustered into two main groups, hereafter called Western and Eastern haplogroups (W and E, respectively, Figure S3) and were discriminated by the deletion of 20 bp at np 601. The Eastern group included only the four RS from Eastern Asia (P. betulifolia, P. calleryana, P. pyrifolia, and P. ussuriensis). The Western group included prevalently European and Mediterranean species, and encompassed all the unknown accessions (N = 20). These latter samples showed the two most common Western haplotypes, namely HT06 and HT09, encompassing different Pyrus spp., and five other haplotypes, four of them shared with accessions of P. communis (Figure S3).

A clear-cut geographic subdivision was also observed when all available sequences from GenBank were included in the network (Figure 3). The dominant haplotypes for the Western and Eastern geographic groups were HT06 (including the reference sequence) and HT04, respectively. Out of a total of 298 samples (95 from this study and 203 from GenBank), 169 accessions belonged to the Western haplogroup (57%), mostly derived from European species (61%), Mediterranean species (16%) and only a small amount from Eastern Asia (1.2% including only P. pyrifolia) and from the Middle East (P. spinosa, P. salicifolia, and P. regelii; 8.3%). These two latter species were absent in the Eastern haplogroup which presented only P. pashia (15.5%) and one sequence from P. spinosa, as representative of Middle Eastern species (Figure 3). In addition, the Eastern haplogroup is prevalently characterized by species from Eastern Asia (68.2%), no Mediterranean species, and only three sequences from European species (P. caucasica, P. communis, and P. korshinskyi).

The unknown accessions of the present study were included in the Western branch, and shared the same haplotypes of P. communis (HT08, HT12, and HT13), P. pyraster (HT15) and of other Pyrus spp. (HT02, HT06, and HT09).

The haplotype relationships were then summarized through principal component analysis (PCA) of the entire dataset. A preliminary PCA (not shown) was uninformative because of the distortion due to two evident outliers (P. betulifolia and P. pashia), both exclusively from the Eastern cluster in Figure 3, which moved far because of their unique haplotypes (HT33, HT35, HT39, HT40, HT41, HT42, HT65, HT66, HT67 for P. betulifolia, and HT10, HT34, HT48, HT49, HT53, HT57, HT58, HT74 for P. pashia). In order to better resolve the graph, a further analysis was performed without those two outliers (Figure 4). In the resulting PCA, some species split up into different quadrants: P. communis in the first, P. pyrifolia in the second, and P. calleryana and P. ussuriensis in the third one. The unknown accessions (grouped together here) were placed between P. communis and a mixed Western Eurasian group located in the center. Magnification of this group highlighted that the unknown accessions were nearby to a subgroup of four closely-related Middle Eastern (P. salicifolia and P. spinosa) and European (P. cordata and P. pyraster) species, which in turn were closer to P. caucasica and P. regelii, than to those from Eastern Asia and the Mediterranean area.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.