Dissecting genetic architecture of grape proanthocyanidin composition through quantitative trait locus mapping

Plant material

The two grapevine populations used in this study have been previously described [25–27]. Briefly, the QTL mapping population (S × G) consisted of a pseudo-F1 progeny of 191 individuals from a cross between two wine grape cultivars, Syrah (S) and Grenache (G) and was maintained under classical local training system (3300 plants/ha plant density) at Montpellier SupAgro Domaine du Chapitre (Hérault, France). The S × G population was planted in two blocks. Each individual from the progeny was planted in two elementary plots (one per block) comprising five plants each. Parental cultivars were also planted in each block with nine and 43 elementary plants for Syrah and Grenache, respectively. The association mapping population (CC) consisted of a core-collection of 141 cultivars maximising agro-morphological diversity for 50 quantitative and qualitative traits, and maintained at INRA Domaine de Vassal (Hérault, France) [26].

DNA extraction, marker genotyping and gene sequencing

DNA extraction and marker genotyping were already described in Adam-Blondon et al. [28]. In order to densify a previous 97-SSR linkage map of the S × G population [25], additional SSR markers, heterozygous both in Syrah and Grenache, were chosen from recent grapevine reference maps [29]. Based on the 12X grapevine reference sequence (Grape Genome Browser http://www.genoscope.cns.fr), we designed primers for candidate gene amplification using Primer3 (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) with default parameters. Primers used in this study are listed in Additional file 1. For SNP analysis, gene fragments were amplified, sequenced and analysed as described in [30].

Linkage map construction

Framework maps were constructed based on the 97-marker linkage map of [25] with 56 additional SSRs. All 153 markers had a genotypic error rate lower than 1.5%, using Tmap as check [31]. Linkage maps were constructed using CarthaGène 0.999R [32] as described in [33] with Haldane mapping function. The "Syrah" and "Grenache" framework maps were composed of 121 SSRs (total length 1118.8 cM), and 133 SSRs (1349.4 cM total length) respectively. The "Consensus" framework map spanned 1256.4 cM based on 153 SSRs among which over 70% allowed segregation in four genotypic classes in the pseudo F1 progeny (ab × ac and ab × cd). Marker order reliability was ensured at LOD 2 threshold. Segregation distortion on genotypic classes was verified by a χ² test according to the segregating type of each marker for the different maps (e.g. for markers segregating as ab × cd in the consensus map, the H₀ hypothesis was ac:bc:ad:bd = 1:1:1:1). Twenty-five markers out of 153 exhibited distorted segregation (P < 0.05) and were mainly grouped on chromosomes 3 (4 markers), 4 (6 markers) and 10 (5 markers). Markers on chromosome 4 exhibited the most significant allelic deviation (P < 0.001) due to segregating distortion between Syrah alleles (aa:ab ~2/3:1/3).

Phenotyping and PA variable construction

Grapes were harvested at maturity (20° Brix). For each genotype, eight representative berry clusters were harvested from the five plants of the elementary plot. Sample homogenisation was based on the accumulation of total solutes (principally sucrose), a major marker of berry development during ripening. Berry density was assessed by floatation in salt solutions [34]. Twenty-five berries with a density between 130 and 160 g NaCl/L were randomly selected. In the present study, we focused on the analysis of berry skin and seeds since PA concentration is quite low in flesh, flesh PAs accounting for only 2-6% of the total berry PA content [35]. Berry skin and seeds were separated, ground in liquid nitrogen and stored at-80°C until analysis. PAs were extracted and analysed by high performance liquid chromatography (HPLC) after acid-catalysed cleavage in the presence of phloroglucinol according to [36]. For the S × G population, both skin and seeds were analysed for 2 consecutive years: skin was analysed in 2005 (1 block) and 2006 (2 blocks) while seeds were analysed in 2006 (2 blocks) and 2007 (2 blocks). For the CC population, skin was analysed in 2005 and 2006 and seeds in 2006.

Phenotypic data analysis

All statistical analyses were performed with R software [37]. We identified the best-fit mixed model for each PA variable through Bayesian information criterion (BIC) in order to extract the best linear unbiased predictors (BLUPs) for genotypic values and to estimate the broad sense heritability (H ² ). Mixed model fit was performed with lme4 package [38]. The mixed model assumption of normality of residual and BLUPs was checked after model fitting by quantile-quantile plot comparing the distribution of residual and random effect predictors to a theoretical normal distribution (Additional file 2). No data transformation was performed for PA variables measured in the two populations. More information regarding phenotypic data analysis and best-fit model for each PA variable is in Additional file 3.

QTL analysis

QTL analysis was performed on the genotypic BLUPs with R/qtl package [39]. Multiple QTL regression was carried out with "stepwiseqtl" function. This approach uses forward/backward selection to identify a multiple-QTL model with inclusion of both main effect QTLs and pairwise interactions. Maximum QTL number was set to 10 for forward selection (max.qtl = 10). Model choice was made via a penalized LOD score (pLOD) which is the LOD score for the model (the log likelihood ratio comparing the full model to the null model without QTL) with penalties on the number of QTLs and pairwise QTL × QTL interactions [40]. For each PA variable, specific penalties for main effect and digenic pairwise interaction terms were derived from 1000 permutations of two-dimensional scan (the "scantwo" function, method = "hk", n.perm = 1000) and penalties at genome-wide error rate of 0.05 were used for multiple-QTL model fitting. The QTL model with the largest pLOD was identified as the most probable one. Once determined the multiple QTL model, we refined QTL position ("refineqtl" function) and estimate R ² for the whole model and each term of the model, the individual LOD score of each term and the genotypic effect ("fitqtl" function). The "lodint" function was used to derive LOD-1 QTL location confidence interval. Allelic effects for consensus QTLs were estimated as described by Segura et al. [41]. Genome scan was performed with a 1 cM step.

Association analysis

Nine candidate genes were selected for association test according to their function and co-localisation with QTLs. Prior to association test, we used R kinship package [42] to perform model comparison among different nested models according to [43] in order to select the best fitted model for association test for each PA variable. Ancestry structure and kinship matrix were estimated based on 20 SSR markers located throughout the whole genome as described in [25].

After model comparisons, we used TASSEL package to perform association tests [44]. Two models were used: one accounting for ancestry structure effect (with General Linear Model, or GLM in TASSEL) the other for both ancestry structure effect and random genetic background effect (with Mixed Linear Model, or MLM in TASSEL). Association tests were performed on BLUPs for skin variables and raw data for seed variables since seed data were available for 2006 only. For GLM analyses, tests were run with 1000 permutations allowing the determination of site-wise P value, whi ch is the probability of a greater F value under the null hypothesis that polymorphism was independent of the phenotype. The adjusted P value (called p_adj_Marker in TASSEL), is the site-wise P value adjusted for multiple tests which takes into account the dependence between SNPs due to linkage disequilibrium. Because each gene was tested independently, we used an additional Bonferroni correction to correct for the number of studied genes (nine) which led to a threshold of 0.0056 for the adjusted P value. As the permutation method is not available for MLM, we used the threshold proposed by Benjamini and Hochberg [45] with q equal to 0.05 which led to a threshold of 0.0039. The effect of minor genotypic frequency and non-normality of observed trait distribution was checked (details in Additional file 4).

Phenotype analysis

PA variable distribution and heritability

For the S × G population, all PA variables showed continuous distribution and transgressive segregation and variation extent was equivalent in the S × G and CC populations (Figure 2, Additional file 5). In agreement with previous studies [15, 16, 46], samples taken in 2006, for which both berry skin and seeds were analysed, displayed different PA composition between tissues as illustrated in Figure 2A with the mean values in the S × G population. For both S × G and CC populations, (-)-epicatechin (epiEx) was the predominant extension subunit in all tissues while (-)-epigallocatechin (egcEx) was only detected in skin. (+)-catechin (catT) was the predominant terminal subunit/monomer in both skin and seed while galloylated units (galEx and galT) were more abundant in seed PA. Each subunit exhibited large variation according to genotype. For instance in skin, (-)-epigallocatechin (instead of (-)-epicatechin) could be the predominant subunit in the extension position (67.9%) (Additional file 5).

PA content variables (concP, concB, concK) reached higher values in seeds than in skin regardless of the unit, as illustrated for concP (Figure 2B), which exhibited the largest difference between tissues. Comparison of composite PA ratio variables showed different range of variation between skin and seed. PAs were on average 8-times shorter in seeds than in skin with wide variation in skin (Figure 2C). All three ratio variables assessing the flux between 2,3-trans and 2,3-cis forms (Ftrancis series) pointed to different kinetics for extension and terminal positions: trans subunits were more abundant in skin for terminal units/monomers (Ftranscis_T) while they were much reduced in seeds when considering extension positions alone (Ftranscis_Ex) or extension plus terminal subunits/monomers (Ftranscis_all, Figure 2D-F). Since the major extension blocks were in cis-configuration, i.e. (-)-epi(gallo)catechin, Ftranscis_Ex was always less than 1 both in skin and seeds. The higher Ftranscis_T in skin conformed to the fact that (+)-catechin was the predominant terminal subunits/monomer in skin. (Figure 2B). Means of each PA variable measured in both skin and seeds were systematically different (paired t-test, P < 0.001, data not shown).

For S × G population, average H ² of PA variables was 0.56 (from 0.24 to 0.82) and 0.44 (from 0.26 to 0.54) in skin and seeds, respectively. No significant difference in H ² magnitude was detected between skin and seeds (t-test, P = 0.053). Nevertheless, higher H ² were observed for skin variables, especially catT and mDP (0.76 and 0.82, respectively, Additional file 5). A high H ² value was also found for these two traits in CC (0.86 and 0.72 for catT and mDP, respectively, Additional file 5).

PA variable correlation

We performed PA variable correlation on genotypic values (BLUPs) from S × G population because the two-year data available both in skin and seeds allowed us to work with a much reduced environmental effect. All three total content variables were highly correlated within a given tissue while significant correlations between tissues were only observed for concB and concK (Figure 3). Among subunit percentage variables in skin (Figure 3), the most noticeable features were the significant negative correlation between egcEx (B-ring tri-hydroxylated subunit) and all other units (B-ring di-hydoxylated subunits) and the positive correlation between (-)-epicatechin and (+)-catechin either in extension position (Ex) or in terminal/monomer position (T). In seeds, the most noticeable feature was the significant correlation of epiT with all other variables (negative with extension units and positive with terminal units) while epiEx was negatively correlated with other subunits. For the same subunit in a given tissue, no highly significant correlation was observed between extension and terminal position, except the negative correlation between epiEx and epiT in seeds (P < 0.001). Significant correlation between subunit percentage variable and composite variables inside a tissue reflected the variable construction (Table 1). Between tissues, significant positive correlations were observed for concK, most of terminal subunits/monomers pairs, galEx and also mDP and Ftranscis_all.

QTL analysis

Global features of PA QTLs

We performed QTL detection with genotypic BLUPs both on consensus and parental maps. In total, 103 vs 43 QTLs and 24 vs 2 digenic epistatic interactions were identified on the consensus map for seed and skin PA variables, respectively (Figure 4). QTLs detected on parental maps were generally also detected on the consensus map except for skin concP, concK, catEx, galEx, epiT, F3pr35 and seed concP where additional QTLs were identified through parental detection (Additional file 6). More QTLs and digenic pairwise interactions were identified on the consensus map than on parental maps, allowing some QTL models to explain more than 80% of the BLUP variance in consensus mapping (Additional file 6), as illustrated in the case of epiT in seeds (Figure 5D). Some loci were involved in phenotypic variation almost exclusively through digenic epistasis such as locus 10@32 for seed concB or locus 14@16.0 for seed epiT (Figure 5). Loci were mainly of additive effect while dominance was predominant at some loci for concK, epiEx, mDP and Ftranscis_T in skin and galEx, epiT, Ftranscis_T and Ftranscis_all in seeds (Figure 4). Among all detected QTLs, only 10 main effect loci overlapped for the same variable in both tissues: 1 for concB, 2 for concK, 1 for epiEx, 1 for galEx, 2 for catT, 1 for epiT, 1 for mDP and 1 for Ftranscis_T (see Figure 5 for some examples). Parental alleles contribution to these common loci was not always consistent across tissues (Additional file 6), which could be an indication of tissue-specific genetic mechanisms. Different genetic architectures were observed for the same PA variable between berry skin and seeds as illustrated in Figure 5: few moderate QTLs (< 3) or no QTL in skin vs several (> 5) small to moderate QTLs with possible involvement of epistasis in seeds (for concP, concB, catEx, galEx, Ftranscis_Ex and Ftranscis_T, illustrated by concB in Figure 5A); many QTLs with involvement of epistasis in skin vs a small number (2) of main effect QTLs in seeds (epiEx in Figure 5B); a major QTL (R ² > 50%) and some QTLs of moderate effect in skin v.s. many QTLs of small to moderate effect in seeds (for catT and Ftranscis_all, illustrated by catT in Figure 5C), or few moderate QTLs in skin vs many moderate QTLs in the presence of a QTL of large effect and epistasis in seeds (epiT, Figure 5D). Conversely, similar genetic architecture between skin and seeds was observed for concK (only moderate main effect QTLs, Figure 5E) and mDP (a major QTL and a few QTLs of moderate effect, Figure 5F). Details regarding position, major allelic effect, LOD score, LOD-1 confidence interval and percentage of explained variation (R ² ) for each QTL are given in Additional file 6.

Association analyses on candidate genes

PA variation extent compared to previous studies

A first characterisation of PA composition in a grapevine pseudo-F1 population was provided by Hernandez-Jimenez and co-workers [46]. Their population was composed of 42 offsprings, derived from a cross between Syrah and Monastrell. In all tissues, the subunit percentage in extension position and Ftranscis-series variables of the Syrah × Monastrell population were of a magnitude and extent equivalent to those of the present study. More divergent results were observed for 1) epiT, which is more abundant in the Syrah × Monastrell population; 2) mDP, which is higher in our study and 3) total content variables, for which the population mean and variation extent was three-to two-fold larger in the present study than in [46]. Syrah, the common parent, behaved similarly in both studies although we observed 10-fold and two-fold higher total PA contents in our study for skin and seeds, respectively. The difference observed in offsprings may result from the fact that the two populations differed by one parent but environmental differences as well as PA extraction and quantification methods might also have affected PA variables, as suggested by the different PA composition of the same parental cultivar. Indeed, mixed model fit suggested that year had a major effect on PA-related variables for both tissues, which was consistent with a previous study where PA content and composition were measured in two cultivars for two consecutive years [53]. The large quantitative variation in PA variables in S × G, of equivalent extent in the CC diversity panel, underlines the interest to implement a quantitative genetic approach on a F1 population for grape PA studies.

Two studies only have characterised PA composition in different grape cultivars [54, 55] with at most a 37-cultivar sample [54]. Different biochemical analyses did not allow for result comparison between this latter study and the present work. Nevertheless, CC in our work was composed of 141 grapevine cultivars of broad geographical origin (from East to West Europe) and was initially defined to maximise the diversity of 50 agro-morphological traits [26]. The PA composition variation in the diversity panel provides thus the potential to refine PA QTLs in a population of larger genetic background.

Multiple QTL mapping in a pseudo-F1 population for grape PA composition

To our knowledge, this study presents the first QTL analysis on grape PA composition with comparisons between skin and seeds of grape berry. This is also the first work on grape using multiple QTL models taking into account both main effects and digenic epistasis during the mapping procedure. QTL mapping in animals has shown that epistasis effects are often large enough to be detected and thus merit a systematic scan regardless of population size, although larger populations (> 500 individuals) allow a more powerful epistasis detection [56]. By employing the multiple QTL mapping approach, we actually showed the important involvement of epistatic interaction in shaping PA composition variation; indeed, some loci were involved in phenotypic variation almost exclusively through pairwise interaction. Our mapping population is of sufficient size (191 individuals) to allow identification of small effect QTLs. However, one should keep in mind that the R ² estimate of individual QTLs is usually overestimated [57] and may have a wide confidence interval [58]. Some of the identified QTLs may therefore be of smaller effect in reality. One should thus be cautious in result interpretation and further identification of causal polymorphism although we did check initially the genome-wide first type error rate.

Allele contribution to individual QTL was mainly due to additive effects between Syrah and/or Grenache alleles (Additional file 6). Parental detection allowed the identification of smaller additive QTLs because of a greater power due to more individuals in each genotypic class compared to consensus detection (i.e., 2 and 4 genotypic classes for parental and consensus detection, respectively). On the other hand, QTL detection on the consensus map allowed us to estimate QTL dominance effect s, i.e. the interaction between allelic classes, but not necessarily with the assumption of a dominant-recessive relationship [59]. In the present work, 9.8% and 30% of QTLs detected in skin and seeds, respectively, had dominance as the major allelic effect (D in Figure 4 and Additional file 6). For example, the locus 8@69 of concB in seed was involved in phenotypic variation almost exclusively through dominance (Additional file 6) and this information would have been overlooked if parental mapping only had been performed.

Genetic architecture of grape PA composition

Tissue-specific genetic architecture for PA composition

In accordance with a previous study which demonstrated the tissue specificity of transcriptional profiles in grape berry [74], the present work illustrates different genetic mechanisms for grape PA composition between skin and seeds: QTLs differed in terms of number, position, R ² and allelic effects. For total content variables, the major QTLs differed in skin and seeds. For subunit percentage and composite variables, important differences between tissues were also observed. Globally, for the same PA subunit percentage variables in both tissues, only 5 QTLs among 74 had overlapping intervals between skin and seeds. Another contrasting feature was observed: in skin, 54% of all QTLs accounted for Syrah additive effect and 78% for Grenache effect whereas in seeds, 74% of all QTLs accounted for Syrah additive effect and 53% for Grenache effect (Additional file 6 and Figure 4). Even for loci identified in both tissues for a given variable, major allelic effect and R ² differed (e.g. the QTL in chromosome 17 for catT, Additional file 6). Our results suggest that seed PA variation is controlled by QTLs of moderate and equivalent magnitude with involvement of epistasis. On the other hand, skin variables are mainly under the genetic control of a few large effect loci with a fluctuating variance unexplained by QTLs.

The different genetic architectures between tissues could result from divergent functional evolution of PAs in these two berry compartments. For fruits in general, ensuring protection of the embryos is essential. Because of their abundance and their ability to protect plants against biotic stresses, PAs and flavan-3-ols might be the major molecules involved in grape embryo protection. In fact, their influence in maintaining seed dormancy has been demonstrated in Arabidopsis [75] and their interaction with phytohormones has also been reported [76, 77]. Therefore, to prevent biological fluctuation due to a single polymorphism mutation, a network with multiple cross-talking actors as a product of evolution without human selection could be postulated in the case of seed PAs, as suggested by the identification of numerous small effect QTLs and the involvement of epistasis. Conversely, skin is the first protective barrier of the grape berry against its environment. In plants, PAs are thought to be involved in self-defence mechanisms [1]. For berry consumers, skin PAs confer flavour to berries and are also responsible for major organoleptic qualities of wine. Consumers in turn could help the plant in seed dispersion or vegetative propagation. However, high quantities of PAs would confer to berries too much astringency and bitterness, which would lead consumers (specially humans) to reject grapevines producing such berries for direct consumption (or wine-making). Human selection in particular could therefore have narrowed down the genetic basis of skin PAs over time and consequently led to a specific genetic architecture with a few large effect QTLs for skin PA variables.

QTL mapping and association analysis as complementary approaches for candidate locus identification

In this work, we used both QTL mapping and association analyses to identify phenotype-marker associations. For identification of phenotype-associated markers, grapevine segregating populations have a greater diversity than populations derived from inbred lines due to heterozygous parental cultivars. However, their genetic background remains relatively narrow compared to diversity panels. On the other hand, QTL mapping may reveal associations undetected in diversity panels due to low allelic frequency. The inconsistency between both approaches sometimes encountered in this work could therefore result from the fact that the available genetic polymorphisms were different between the two populations: causal polymorphisms in one population might be monomorphic in the other. In addition, our analysis focused on genes of known function co-locating with QTL while other candidates could underlie QTL intervals. Besides time-consuming fine mapping, candidate genes can also be selected by combining QTL results with other data such as transcriptomics. Nevertheless, LAR1 gene evoked a particular interest through association test. Complementing functional studies performed on a single cultivar [14], we provide here additional confirmation of LAR1 gene involvement in grape PA composition through a diversity panel study.