Genetic structure in cultivated grapevines is linked to geography and human selection

Diversity

Our sample of 2,096 cultivated genotypes (Additional file 1: Table S1) displayed from 5 to 34 alleles per locus, with a total of 324 alleles over the 20 loci, an unbiased expected heterozygosity (He) of 77% and a mean polymorphism information content (PIC) of 0.740 (Additional file 2: Table S2). Owing to the large number of loci with moderate allele frequencies, the single parent exclusion probability was quite high (7 × 10^-6).

Population structure

The similarity pattern among the 10 STRUCTURE analysis replicates (Figure 1) and Evanno’s ΔK_s statistics (Additional file 3: Figure S1) indicated K_s = 3 and K_s = 5 as the most pertinent levels of population subdivision. No converging solutions were identified for the subsequent Ks levels (6 to 12), which were therefore not further considered with STRUCTURE.

Using a threshold of >85% for group assignation, 1,001 genotypes (out of the 2,096) were assigned to a cluster at K_s = 3 and 817 at K_s = 5. The proportion of admixed genotypes was thus large, i.e. 52% and 61% of the total number of cultivars, at K_s = 3 and K_s = 5 respectively.

Geographic origin and viticultural traits of the identified subgroups

STRUCTURE clustering at K_s = 3 highlighted three well-distinct groups (Figure 2, Additional file 4: Table S3):

a Western Europe group (S-3.1) of wine cultivars (93%) containing 55% of the Western and Central Europe genotypes, without any Asian, Balkans or Maghreb cultivars;

a East group (S-3.2) mostly composed of table cultivars (71% of table grapes and 9% of ‘double-use’ cultivars), including 96% of the Far- and Middle-East genotypes, notably all genotypes from Uzbekistan (n = 33), Afghanistan (8), Tajikistan (4), Turkmenistan (4) and Iran (23), as well as 66% of the Eastern Mediterranean and Caucasus cultivars, and almost no Western and Central Europe cultivars (less than 2%). Interestingly, 43% of the cultivars from Maghreb were positioned within this group;

a Balkan and Eastern Europe group (S-3.3) of mostly wine cultivars (71%), with 34% and 60% of the varieties from the Balkans and East Europe, and less than 4% Western Europe cultivars.

The relationship between the different clusters and cultivar geographical origin was analyzed (Figure 3). The Eastern Mediterranean & Caucasus (EMCA), Middle & Far East (MFEAS) and Maghreb (MAGH) geographic groups (Table 1) were composed almost exclusively of genotypes clustered in the S-3.2 STRUCTURE genetic group; the Western & Central Europe (WCEUR) and the Balkans (BALK) geographic groups were also mainly composed of genotypes from their corresponding STRUCTURE group (S-3.1 and S-3.3, respectively). By contrast, the Russian, Iberian and New World cultivars were distributed in several STRUCTURE genetic groups. For K_s = 3, it is worth noting that 100% of the Italian cultivars were assigned to the “admixed” class, while the Middle and Far East group displayed a very low level of admixture (3.6%). Apart from its meaningful geographic distribution, the admixed group was composed of even proportions of wine or table cultivars, of black, red or white color grapes, and of aromatic or non-flavored grapes.

The clustering at K_s = 5 (Figure 2, Additional file 4: Table S3) identified in addition an Iberian and Maghreb group (S-5.1), and a group comprising mostly table grapes (80%) of recent origin, also called “obtentions”, from Italy and Central Europe (S-5.4). The group S-5.1 derived partially from the S-3.2 group (41% of the varieties), with Iberian varieties composing 69% of the group. The group S-5.4 mostly derived from the admixed Ks = 3 group (78% of the varieties).

Axes 1 and 2 of a PCA on SSR data of the genotypes belonging to Ks = 5, explaining 30.3% and 21.4% of the total variance respectively (Figure 4), clearly separated the WCEUR, BALK and East groups. The two additional groups at (S-5.4 and S-5.1) were separated by PCA only on axes 3 and 4 respectively (not shown), which explained 14.3% and 8.8% of the total variance. Thus genetic clustering at K_s = 3 appeared more structuring than the one at K_s = 5.

To support the STRUCTURE analysis, we performed a Ward clustering (Figure 5; the full dendrogram for the 2,096 cultivars is given in Additional file 5: Figure S2). The Ward and STRUCTURE clustering were found consistent (Additional file 6: Table S4), with a correspondence among clusters composition of 90% and 87% for K = 3 and K = 5 respectively. In addition to the main partitions already explored at Ks = 3 and 5, the Ward clustering level at K_w = 12 (Additional file 7: Table S5) identified local germplasms, groups of cultivars with a particular characteristic (white, seedless or muscat flavor grapes), or parentages linked to human selection and breeding (next paragraph).

Genetic diversity and family structure within and among clusters

The average genetic differentiation among STRUCTURE groups for K_s = 3 and K_s = 5 was D _est  = 0.166 and 0.213 respectively (harmonic means; in both cases, standard deviation D _est  = 0.005). The largest differentiation between pairs of subpopulations was found between the Western-Central Europe and the Eastern groups, for both K_s = 3 and 5, with D _est  = 0.217 and 0.256 respectively). The Eastern group and the group from the Iberian Peninsula and Maghreb displayed a D _est of 0.139, the lowest of all comparisons (Additional file 8: Table S6).

Since the presence of family groups may affect population genetic structure, we explored the distribution of cultivar family relationships within and among STRUCTURE groups. By comparing all possible genotype pairs, the ML-relatedness software outlined in total 1,069 likely parent pairs involving 1,099 putative parents.

All the analyses above allowed us to characterize the different groups (Table 4 and Additional file 9: Table S7) for their genetic diversity, family relationships and phenotypic characteristics, and finally propose an interpretation of grapevine genetic history and geographic partitioning.

Geography and history

The three main clusters revealed by our study, both with STRUCTURE and Ward’s methods, confirmed previously obtained molecular results [5, 9] and the eco-geographic grouping proposed by Negrul [38], in particular the correspondences between the “proles” occidentalis and S-3.1/W-3.1 groups, the pontica and S-3.3/W-3.3 groups, and the orientalis and S-3.2/W-3.2 groups. Our results allow us to subdivide these clusters according to cultivar putative geographical origins: i) West and Central Europe (S-3.1), ii) East Mediterranean, Caucasus, Middle and Far East (S-3.2), and iii) Balkans and East Europe (S-3.3). Clustering at K = 5 identified two new groups, an Iberian Peninsula group and a group of table grape obtentions with Italian Peninsula and Central Europe origins.

Genetic characterization of the groups clearly showed the East table grape group (S-3.2 and S-5.2 for K = 3 and 5 respectively) as the most diverse in terms of mean number of alleles, number of private alleles, and non-biased heterozygosity. This is consistent with the hypothesis that grapevine domestication initially occurred in Eastern regions (Caucasus and Fertile Crescent) as suggested earlier [2–4, 9], repeatedly introducing genes from the wild. The high frequency of private alleles in S-3.2 and S-5.2 could also be explained by a history of limited exchanges from East to West, as attested by the high differentiation values (D_est) between these regions, and a slower development of grape breeding in the East, as indicated by the low frequency of family-related genotypes in that region as compared to other regions, revealing a weaker selection bottleneck effect there. However, given the high genetic diversity of grapevine at all subdivision levels, the selection and breeding bottlenecks seem in general weak for this crop.

The second most diverse group was the West and Central Europe wine grape group, probably as a result of this area’s long history of grapevine cultivation and development, in combination, as already stated by other authors, with gene flow from local wild or primo-domesticated grapevines [9, 10, 18]. The Balkans and East Europe cluster also formed a well identified STRUCTURE group with an intermediate diversity. The two additional groups at Ks = 5 (the Iberian Peninsula group and the group of table grape obtentions), appeared as secondary groups with a lesser global diversity.

More generally, the full hierarchical partitioning obtained with the STRUCTURE and Ward methods as well as the D_est differentiation statistics appeared consistent with historical data, such as the diffusion of viticulture around the Mediterranean Sea, with one route connecting Eastern (W-3.2) to Western Europe through the Balkans and Central Europe (W-3.3, W.3.1) [2, 9], and a Southern route to the Maghreb and Iberian peninsula (W-3.2 /W-5.1 / W-12-4).

The Balkans and Eastern Europe group and the Western and Central Europe group were both characterized by a large proportion of genotypes belonging to one STRUCTURE group only, probably corresponding to separate regional grapevine cultivar development and selection. In contrast, other regions as Russia and Ukraine, the Iberian Peninsula, and the New World countries, contain a mix of two or three STRUCTURE groups, in relation to their regional position. In particular, varieties found in Russia and Ukraine appear to have either East (S-3.2), Balkans and East Europe (S-3.3) origins, consistently with what we know of the centralizing impact that Russian agricultural research had during the Soviet period [39]. Similarly, the Iberian peninsula group include cultivars from West Europe (S-3.1), East (S-3.2) and Maghreb (S-5.1) as well as a high proportion of admixed genotypes, in coherence with the long historical exchange relationships this region had both with Europe and North Africa. Based on maternally inherited chloroplast markers, Arroyo-Garcia et al. suggested that the Iberian Peninsula could be a secondary center of domestication [9]. Our results add a new view of Spain and Portugal as platforms of centralization, intermixing and exchange of varieties throughout history.

Finally, at K_w = 12, the genotypes from the eastern regions (proles orientalis[38]) further subdivided into two sub-groups, one mainly composed of wine cultivars of Caucasian origin (including Georgia, Armenia, Azerbaijan and Turkey, W-12.12), and the other comprising table cultivars from Central Asia (Tajikistan, Uzbekistan, Turkmenistan) together with Iran and Afghanistan (W-12.11). The separation of these two groups may be a trace of divergent selection for the main local use for grapevine (table vs. wine). On the other hand, the absence of admixture in the Middle and Far East group, in particular for the 72 cultivars from Uzbekistan, Afghanistan, Tajikistan, Turkmenistan and Iran, and the high K scores of its genotypes may be an indication that the corresponding center of domestication was larger than formerly believed (several authors indeed placed it in a geographic region between the Black Sea and Iran [2, 3, 40, 41]), an hypothesis already proposed in 1976 by Olmo [12], but not confirmed by later studies. It is difficult to decide between these two scenarii since the information available on grapevine crop development is quite limited for Central Asian countries.

A large proportion of admixed genotypes was found by STRUCTURE, both at K_s = 3 and K_s = 5. A previous study on maize indicated that, in crops, STRUCTURE grouping is generally coherent for first cycle inbreds with simple parentage relationships, while the presence of multiple levels of family relationships and cohort overlapping in more advanced breeding systems leads to different grouping possibilities and low STRUCTURE stability [42]. We can infer that our sample contains both types of material, with a number of ancient varieties anchoring the main clusters (founders), and recent breeds complicating structure resolution. The stability of K_s = 3 and K_s = 5 groupings and the individual percentage of cluster ancestry allowed us to discriminate among these two types of materials. The geographic distribution of the admixed genotypes is not “random” (Table 1): the Middle-Far East is the region displaying the lowest level of admixture, while Italy in particular and secondly the Iberian Peninsula, display the larger proportion of admixed genotypes. We were unable to find other traits characterizing the admixed group: it is composed of even proportions of phenotypic classes of grape use, berry color, flavor, berry seed number, or sex.

While confirming and reinforcing the observation of geographic structure of the cultivated gene pool already described by other authors [5, 9, 16, 18, 38], our results are also coherent with the study of Cipriani et al. [17] suggesting that Italian varieties present weak or no structure: indeed in our study the Italian cultivars appear to be admixed, probably as a result of the inter-regional exchange role that Roman culture has certainly played.

Phenotypic traits

Our results also provide information about the effect human selection on morphological traits had on shaping the genetic diversity of cultivated grapevine. Table and wine grapes have different berry size and bunch shapes, both important traits used for cultivar classification [22]. Table and wine grapes are clearly separated by STRUCTURE at K _s  = 3. At K _s  = 5; only the group including Iberian and Maghreb cultivars (S-5.1) is composed of a mix of table and wine cultivars, which is likely the result of artificial selection and intimate cultivars intermixing in this area.

The black color of berries is considered as an ancestral trait compared with the other colors, both at phenotypic [1] and molecular level. The molecular basis of the apparition of red, rose, grey and white berry colors has been previously documented [43–45] and the diffusion of the major causal mutations – Gret1 insertion and K980 mutation – within the cultivated compartment was described by Fournier-Level et al. [19]. In the present STRUCTURE analysis, the Central and West Europe subgroup (S-5.3) is composed of a majority of black cultivars. This can be explained by the isolation of these regions from the Eastern cultivars, by local domestication and gene flow from endemic black-berried V. v. sylvestris, or human selection. All other subgroups include a large number of white cultivars, reinforcing the idea of a wide and strong diffusion of Gret1 over the whole geographic range of grapevine [19]. Most of the intermediary phenotypes (red, rose and grey) are concentrated within two groups: Balkans and Central Europe (S-5.5), and East (S-5.2), confirming these regions as putative sources of color variation [19].

The geographical origin of Muscat flavor is assumed to be Greece or the Balkan Peninsula [46, 47]. Thereafter, human selection aimed to spread this desirable trait in both table and wine grapes [20]. With STRUCTURE, we found the majority of Muscat founders within the Central Europe table group (S-5.4). Only a small number of them were involved in breeding, essentially in the Balkans, forming kingroups with other known parents such as Chasselas.

Seedless cultivars clustered essentially with cultivars of Turkish, Caucasian and Asian origins, belonging to the proles orientalis[38], coherently with available historical data about their origins from Turkey and Near-East [48].

Plant material

The plant material was composed of cultivated grapevine varieties belonging to V. vinifera subsp. vinifera held in the INRA grape repository at Vassal (France). This collection includes 3,727 accessions available as field-grown plants and genotyped with 20 microsatellite markers [8].

Geographic assignation

Geographic origin of referenced cultivars was derived from general bibliography on ampelography and viticulture [53, 54]. For non-referenced cultivars, the origins were estimated on the basis of the accession origin. Recently bred cultivars (e.g., Tarrango) were assigned to the breeder’s country (in this example, Australia) and not to the countries of origin of their progenitors (in this case, Portugal and Turkey). Countries of origin were grouped as shown in Table 1.

Cultivar characterization

DNA extraction and genotyping

DNA extraction, PCR amplification and genotyping of microsatellites were carried out according to This et al. [4] and Laucou et al. [8]. The 20 nuclear microsatellite loci were chosen according to their polymorphism level and their position on the linkage groups [23]. Differences of one base pair between alleles at one given locus were double-checked by re-amplification and re-analysis; a test for the presence of null alleles was also carried out [57].

Since genotypes with only one or two allele differences represent closely related material, such as clones or recent mutants, and provide very little additional information to the analysis, these were considered redundant and not taken into account. Thus, out of the 2,323 single genotypes identified in Laucou et al. [8], 2,096 genotypes presenting no missing data and at least three allele differences were analyzed in our study for further structure and clustering analyses (full list and characteristics are given in Additional file 1: Table S1).

Statistical analysis

Main diversity statistics for the 20 microsatellite markers, such as the total number of alleles, expected heterozygosity and total gene diversity [58] were calculated using Genetix [59]. A private allele index adjusted for sample bias was estimated with ADZE [60], following a generalized rarefaction approach. Genetic differentiation D_est was computed using the SMOGD software [61, 62], based on the method of Jost [63]. Confidence intervals were calculated using 1000 bootstraps, Polymorphism Information Content (PIC) according to Botstein et al. [64], and the single parent exclusion probabilities according to Jamieson & Taylor [65].

The dataset of 2,096 unique genotypes was used to run a series of tests, with the Bayesian method implemented in STRUCTURE [27], in order to find the best model to infer population structure (with or without admixture, correlated allele frequencies, or prior information about sampling locations) and the best K_s level of population subdivision, with K_s varying from 2 to 12. Within STRUCTURE, we allowed an iterative process with a burn-in phase of 5 × 10⁴ iterations, and a sampling phase of 5 × 10⁴ replicates. Ten replicates of each assumed K_s-level subdivision were compared to estimate group assignation stability.

We first evaluated the different models of admixture, allele frequencies and prior population information available in STRUCTURE. The most appropriate model to interpret our data appeared to be the uncorrelated allele frequencies and prior geographic information model, which showed a better stability between runs for K_s = 3 and 5, and a lower variance for K_s = 5, as compared to other possible STRUCTURE models (Additional file 3: Figure S1).

Since the geographic groups are not equally represented in Vassal, a second analysis was run to measure a possible sampling effects, in particular the bias that could arise from oversampling one region or one family group. This was tested by running STRUCTURE on two different set of genotypes (Table 1), the full set of 2096 genotypes (set1) and a sub-set of 888 genotypes randomly drawn to constitute equally-sized, geographic origin-based groups (set2).

Finally, the most probable uppermost level of structure subdivision between the successive K_s values was estimated with two methods: 1) the calculation of Evanno’s delta-K as the second order change in the likelihood function divided by the standard deviation of the likelihood [66, 67], and 2) the similarity coefficient between each pair of runs, which provides an evaluation of the stability of the solutions between runs.

Genotypes were assigned to a cluster when 85% or more of their inferred genome belonged to the cluster, the genotypes with a lower score being considered as “admixed”. The chosen clusters for each Ks level were then labeled according to a three digit code (e.g. S-2.1, S-2.2,…, S-5.5) for further geographic and phenotypic characterization. A graphical display of the individual and group distances was obtained with a Principal Component Analysis (PCA) using the package adegenet implemented in R [68, 69].

To validate the STRUCTURE clustering, we compared its output with that obtained using a less constrained method of clustering. Odong et al. [30] highlighted that STRUCTURE and Ward’s method [28] are convergent and complementary. Thus we used Ward’s method to evaluate the distances between clusters minimizing the sum of squares of any two clusters at each step. Using Ward dissimilarity matrix, we built a dendrogram with DARwin software [29]. The advantage of Ward clustering is to provide details of the relationships at any level, as close as family levels. One disadvantage is that it does not deal with admixed genotypes. We indexed the Ward subdivision levels as K_w, and labeled the subgroups accordingly (W-2.1, W-2.2,…, W-12.12).

Genetic structure partitioning between and within groups - Family relationships

To estimate the part of the population genetic structure due to parentage, we first calculated the most probable family relationship among each pair of genotypes using the ML-relatedness software [70]. Genotype pairs (half of a 2,096 x 2,096 matrix minus the diagonal [(n²-n)/2)] = 2,195,569 couples) were declared either unrelated or family-related, this latter category grouping full sibs, half-sibs, and parent-offsprings, in order to lower the chance of false attribution. Only non-ambiguous relationship assignations (according to 99.9% confidence intervals calculated in 100 mating population simulations) with an experimentally determined LOD score > 9 were taken into account. For each of the above subdivisions we also calculated a weighted average relatedness r ² .

The results of this analysis were then assigned to categories of STRUCTURE subdivision (e.g. within or among subgroups), according to the group of each parent. Within-group average relatedness was estimated with the formula of Queller and Goodnight [71] and its standard deviation with 100 jackknifes over loci, using the RERAT software [72].

Phenotypic evaluation of the different K levels

To interpret the population structure in terms of cultivar utilization, movement and history, each subpopulation was finally characterized for its flower and fruit traits and for its geographic origin. Group names were ultimately based on their main characteristics.