Association mapping reveals the genetic architecture of tomato response to water deficit: focus on major fruit quality traits

Highlight Tomato quality could be improved under deficit irrigation while maintaining yield. The underlying genetic architecture is polygenic and varies with water availability. Candidate genes related to primary metabolism were identified.


Introduction
Global water scarcity will constitute a crucial challenge in the coming years (Jury and Vaux, 2005). Agriculture, which is consuming up to 80% of the worldwide water resources through irrigation, has to move towards a more sustainable use of water (Rost et al., 2008). Utilization of advanced irrigation strategies and development of drought-adapted crops are among the solutions to solve this dilemma (Fereres and Soriano, 2006;Costa et al., 2007).
Beyond these concerns, deficit irrigation practices constitute a way to manage fruit flavor by exploiting the morphological, physiological, and molecular changes (referred to as 'phenotypic plasticity') occurring in water-stressed plants (Ripoll et al., 2014). Under water deficit, plants close their stomata to limit transpiration, impacting resource availability from photosynthetic sources, which may result in a decrease in number and/or size of the fruits. On the other hand, a mild water deficit tends to shift photo-assimilate partitioning towards synthesis of antioxidant compounds (in particular vitamin C) involved in defense against stress-induced reactive oxygen species and compatible solutes (including sugars and acids) involved in osmotic adjustment (Lemoine et al., 2013;Albacete et al., 2014;Osorio et al., 2014). Evidence for the efficiency of deficit irrigation to concentrate the major flavor and nutritional components in fleshy fruits (mainly sugars, acids, and antioxidants), either by a concentration or an accumulation effect, was obtained in many species such as tomato (Kirda et al., 2004;Zheng et al., 2013), grapevine , apple (Leib et al., 2006), and mango (Durán Zuazo et al., 2011). However, these studies focused on a small number of genotypes, while responses to deficit irrigation seem to be highly genotype dependent (Ripoll et al., 2016a, b).
Gene expression studies have revealed hundreds of genes involved in plant survival under severe water limitation, but usually associated with detrimental effects on yield under a realistic drought scenario (Tardieu, 2012;Bac-Molenaar et al., 2016). These studies focused on model species, mainly Arabidopsis thaliana (Seki et al., 2002;Des Marais et al., 2012) and cereals (Langridge, 2006;Barnabas et al., 2007). Up to now, the identification of the genetic determinants of drought response from the natural diversity of fleshy fruit crops remains limited. Quantitative trait locus (QTL) mapping might be particularly valuable to address this question (Des Marais et al., 2013).
Two complementary approaches are commonly applied to dissect genotype by environment interactions into their underlying QTLs (QTL by environment interactions). The first one consists of computing the effects of a given QTL across the environmental conditions using multivariate QTL mapping models (van Eeuwijk et al., 2010;El-Soda et al., 2014b). The second one uses the construction of composite variables measuring phenotypic plasticity and univariate mapping models (El-Soda et al., 2014a). With both approaches, QTLs can be classified according to the prevalence of their effect under the different conditions. A QTL is considered 'constitutive' when its effect is conserved whatever the environment. QTLs whose effect is not significant in every environment are called 'specific', while the effect of 'interactive' QTLs changes direction ('antagonist') or intensity ('differential') according to the environment. With the availability of a high-throughput genotyping assay, this classification can be considered in crop species via conventional linkage mapping (Malosetti et al., 2007;Verbyla et al., 2014) as well as genome-wide association studies (GWASs) Saïdou et al., 2014). A GWAS has the advantage over linkage mapping that it allows exploration of the genetic diversity and the numerous recombination events present in germplasm collections and may lead to higher resolution mapping if the LD (linkage disequilibrium) is low enough in the population (Brachi et al., 2010;Korte and Farlow, 2013;El-Soda et al., 2015;Pascual et al., 2016).
In tomato (Solanum lycopersicum L.), QTLs were mapped for fruit quality traits measured under optimal watering conditions using linkage (Causse et al., 2001;Saliba-Colombani et al., 2001;Tieman et al., 2006;Zanor et al., 2009b;Capel et al., 2015) and association mapping Ruggieri et al., 2014;Sauvage et al., 2014;Sacco et al., 2015). The studies of QTLs by water regime interactions focused on introgression lines between the cultivated tomato and its wild relatives (mainly S. habrochaites and S. pennellii), leading to low mapping resolution (Semel et al., 2007;Gur et al., 2011;Arms et al., 2015). Recently, we analyzed QTLs by watering regime interaction in a segregating population derived from a cross between a small-and a large-fruited S. lycopersicum accession . A total of 56 QTLs were identified for 19 traits, among which 20% were interactive between the control and deficit watering regimes. Nevertheless, these QTLs were limited to the allelic diversity present in the two parental accessions, and the confidence intervals were broad.
The aims of the present study were (i) to explore the pattern of genotype by watering regime interaction in a GWAS panel with a broad genetic basis (including S. pimpinellifolium, S. lycopersicum var. cerasiforme, and admixture genotypes) grown under two different watering regimes in two locations and phenotyped for 27 traits; (ii) to identify with a high resolution QTLs and QTL by watering regime interactions in this collection; (iii) to combine the results with those obtained in the bi-parental progeny to draw an accurate picture of the genetic variability and the genetic determinants of tomato response to water deficit; and (iv) to identify candidate genes related to the variation of major fruit quality traits under water deficit by dissecting some of the QTLs.

Plant material
The population consisted of 141 accessions (2-46 g FW) encompassing the genetic diversity of the cultivated small fruit tomato. Among these, 105 accessions were previously investigated in Blanca et al. (2015). Preliminary genetic analysis of our collection confirmed the genetic structure described by these authors, with clusters reflecting the species and the geographic origin of the accessions (see Supplementary Fig. S1A-D at JXB online). Ten accessions were S. pimpinellifolium (SP; closest wild ancestor of the tomato) originating from Peru and Ecuador. A total of 110 accessions were S. lycopersicum var. cerasiforme (SLC) originating mainly from South America. Finally, 21 accessions belonged to a mixed genetic group mainly including commercial cherry tomatoes and admixed genotypes between SP, SLC, and S. lycopersicum var. lycopersicum. A description of the accessions and their origin is available in Supplementary Table S1. The genetic groups (SLC, SP, and mixture) are used below in the statistical analysis.

Experimental design
The plants were cultivated with the same experimental design as in Albert et al. (2016). Plants were grown in a heated glasshouse in INRA Avignon (Avi, France) from March to July 2014 and in an unheated plastic greenhouse on the experimental site of the seed company GAUTIER Semences in Agadir (Aga, Morocco) from December 2013 to March 2014. Two watering regimes were applied to the plants: control (C) and drought (D). The control treatment was set according to ET (evapotranspiration) and the cultural coefficient for tomato under greenhouse conditions (FAO Water, 2015). A maximal drainage of 25% and a relative humidity of the substrate of 65% were established in the control pots. Drought treatment was applied progressively after flowering of the second truss of the earliest accession. Watering was first reduced by 25% compared with the control for 1 week and then reduced by 60% until the end of the experiments. Relative humidity of the peat substrate was controlled with GRODAN ® moisture probes and monitored between 25% and 30% in drought pots. In both experiments, two plants per watering regime per accession were randomized in the greenhouse.

Plant and fruit phenotyping
A total of 27 traits were assessed in the GWA population as described in Albert et al. (2016). Flowering date (Flw, days after sowing), stem diameter (Diam, mm), leaf length (Leaf, cm), and truss implantation height (Ht, cm) were measured on each plant both in Avignon (sixth truss) and in Agadir (fifth truss). Plant fruit number (Nbfruits, all fruits from the third to sixth truss) was measured only in Avignon.
Fruit quality measurements were carried out on a minimum of 20 mature fruits per accession per watering regime harvested daily on the third to the sixth truss. All the fruits were weighed (FW, g) and their firmness was measured with a Durofel device (FIR). Only in Avignon, fruits were pooled in three groups in each watering regime. Half of the fruits of each pool were used to assess dry matter weight (DMW, %), pH, and soluble solid content (SSC, °Brix). From the second half of the fruit replicates, pericarps were crushed in liquid nitrogen and assayed for total vitamin C content (VitCFM) according to the microplate method described in Stevens et al. (2006), for sugar content (glucose and fructose) according to the enzymatic method described in Gomez et al. (2007), and for organic acid content (malic and citric) according to the HPLC method reported in Wu et al. (2002). The different metabolite concentrations were expressed relative to fresh matter (g 100 g -1 of FM) and relative to dry matter (g 100 g -1 of DM). Yield (g per plant) was computed by multiplying average fruit FW by average fruit number per plant.

Plant genotyping and SNP filtering
The GWA population was genotyped using the Tomato Infinium Array developed within the SolCAP project (http://solcap.msu. edu/) (Hamilton et al., 2012;Sim et al., 2012). The maximum rates of missing data were fixed at 25% per accessions and 10% per SNP. A minor allele frequency threshold of 0.04 was applied to discard markers with very rare alleles according to Aulchenko et al. (2007). After filtering, the set of markers was constituted of 6100 SNPs. Prior to any genetic analysis, the remaining missing genotypes were replaced by the allele frequency of the major allele. The SNPs were renamed according to their positions on the tomato genome (SL2.50), as S01_58000085 at base pair 58 000 085 on chromosome 1 (Supplementary Table S2).

Statistical analysis of the phenotypic data
All statistical analyses were performed using R (R Development Core Team, 2012). Because fewer and different traits were measured in Agadir experiments, data from both locations were analysed separately (Pearson correlations for the common trait means available in Supplementary Table S4-all significant). Prior to the ANOVAs and when distributions were skewed, phenotypic data were normalized using Box and Cox transformations. The ANOVAs were performed according to the following model: Y ijkl was the phenotypic value of accession j from genetic group i in watering regime k, µ the overall mean, Gr i the fixed effect of genetic group i, Gr i (G j ) the fixed effect of accession j nested in genetic group i, W k the fixed effect of watering regime k, and e ijkl the residual error effect. No significant microenvironment pattern was identified and we chose not to include any spatial effect in the model. When the interaction Gr×W was significant, we computed a Tukey's post-hoc test to compare the means. Then, in both watering regimes, restricted maximum likelihood estimates of the genetic and residual variances (σ 2 G and σ 2 e ) were computed with a second linear model: Y ijk= µ+Gr i +Gr i (G j )+e ijk (Gr j fixed, G i and e ijk random). Broad-sense heritabilities (H 2 ) were calculated under both watering regimes as the ratio between the genetic variance and the total phenotypic variance: H 2 =σ 2 G /σ 2 Total , with σ 2 Total =σ 2 G +1/n×σ 2 e (with n the number of replicates per accession). Spearman coefficients estimated the correlations between H 2 and σ 2 G under drought and control conditions for the same trait. Average values per accession in each watering regime and location were used for subsequent analyses. Plasticity was computed on the accession means as: ∆ki=(D ki -C ki )/C ki , with ∆ki the plasticity value for trait k and accession i, D ki the mean of trait k under drought condition for accession i, and C ki the mean of trait k under control condition for accession i.

Construction of kinship and structure matrices
We performed a principal co-ordinate analysis (PCoA) on the genotype matrix. The co-ordinates of the accessions on the first three components are available in Supplementary Table S3 and displayed graphically in Supplementary Fig. S1. A kinship matrix (K) based on identity by state among the 6100 SNPs was estimated.

GWA mapping
Average values for each trait following the transformation giving the least skewed distribution were used in the mapping models. GWASs were performed using correction for population structure (PCoA) and modeling genetic variance with the kinship matrix (K). Two mixed models were implemented.
First, the bivariate multitrait mixed model (MTMM) developed by Korte et al. (2012) to take into account the correlation structure of multienvironment data sets and increase the detection power was implemented. The MTMM approach includes two different tests: (i) the 'global test' compared a model including only the genotype effect with a null model to identify markers with common effect between watering regimes ('constitutive QTLs'); and (ii) the 'G×W test' compared a full model with a model including only the genotype effect to identify markers with an interactive effect between the watering conditions ('interactive QTLs'). SNPs with a P-value <10 -4 were considered as significant. From each test, the percentage of variation explained by the marker (individual PVE for each significant marker) was computed.
Secondly, the univariate multilocus mixed model (MLMM) developed by Segura et al. (2012) to increase the detection power for polygenic characters was used to identify associations for each trait under each watering regime ('specific QTLs') and for the ∆ values ('interactive QTLs'). We implemented a new model selection criterion in the MLMM framework to allow for a more permissive detection threshold to compromise between type I (false-positive) and type II (false-negative) errors, while limiting the number of cofactors selected to avoid overestimation of the P-values due to the relatively small size of the population. Models with a maximum of five cofactors all having a raw P-value <10 -4 were retained. From the optimal model selected, the percentage variation explained by the selected markers (global PVE for all the significant markers) was computed for each trait.
For all the QTLs identified, we computed phenotypic effects under both watering conditions as: (Minor allele mean−Major allele mean)/2. Among the interactive QTLs, we distinguished between 'antagonist QTLs' (effect changing direction according to the watering regime) and 'differential QTLs' (effect changing intensity according to the watering regime).

Linkage disequilibrium estimation and confidence interval definition
To define intervals around QTLs, we used a strategy based on LD between pairs of markers inspired from Cormier et al. (2014). We used the r 2 estimator implemented in the package 'genetics' (Warnes and Leisch, 2012) to assess LD between marker pairs. First, we performed LD calculation between 100 000 randomly chosen pairs of unlinked loci (on different chromosomes). The 95th percentile of the unlinked-r 2 distribution equal to 0.28 was considered as the critical LD threshold. Then, for each significant marker, we computed LD with all the markers upstream and downstream on the same chromosome. We defined the lower (upper) boundary of the interval as the last marker downstream (upstream) on the chromosome that presented an LD with the significant marker above the 'critical LD' threshold. For the QTLs detected with the MTMM procedure, when two markers presented a LD higher than the LD threshold, we considered them as a unique QTL. The number of genes within each interval was identified from the tomato genome (ITAG2.4).

Comparison between linkage and association QTLs and identification of candidate genes
For the comparison with the QTLs detected in the recombinant inbred lines (RILs) grown under the same conditions and phenotyped for the same traits , we projected the QTLs detected in both populations onto the tomato genome (SL2.50). In the comparison, we considered related traits as a single trait: pH, malic acid, and citric acid contents were grouped as 'acids', and SSC, glucose, and fructose contents as 'sugars'. Besides, whatever the QTL type ('interactive', 'constitutive', or 'specific') and the location of the trial, we considered that a single QTL was present when the intervals overlapped between RIL and GWA QTLs.
We then focused on the QTLs for vitamin C, sugar, and acid content including <100 genes to identify putative candidate genes with a reasonable confidence. Under those QTLs, we refined the set of candidates by selecting the genes expressed in tomato fruits according to gene expression data published by the Tomato Genome Consortium (2012). Then, we examined their functional annotations and focused on genes with annotations corresponding to related functions. Finally, we screened the polymorphism data obtained through the whole-genome resequencing of four accessions of our GWA population chosen to represent a large range of the molecular variability present in small fruit tomato : Cervil (13.3× sequence depth), Criollo (8.1×), LA1420 (12.5×), and Plovdiv (12.2×). First, we considered the nucleotide variants with moderate (non-synonymous polymorphisms in coding regions) to high (modification of splice sites or start/stop codons) effect on the protein sequence (detected using SnpEff; Cingolani et al., 2012). Then, the predicted impacts of the variants on the protein function were assessed using the web interfaces of PROVEAN (http://provean.jcvi. org/seq_submit.php) (Choi and Chan, 2015).

Dissection of the phenotypic variations in the GWA population
In the variance analysis, the part of the total variation attributed to the genotype effect was predominant (35-80%, all P-values <0.001) compared with the one attributed to the genetic group (0-15%, all P-values <0.05) and the watering regime (0-28%, significant for 17 traits), except for leaf length in Agadir and stem diameter in Avignon and Agadir ( Fig. 1; Supplementary Table S5). For those vigour traits, the watering regime represented 48-61% of the total variation.
The genetic group by watering regime interactions represented <2% of the total sum of squares for all traits and was non-significant for 12 traits. The eight significant traits were Diam.Aga, Leaf.Avi, Leaf.Aga, Ht.Avi, FW.Avi, FW.Aga, FIR.Aga, and VitCFM.Avi. Tukey's post-hoc test indicated that these interactions were mainly driven by a singular behavior of the SP group in response to water deficit ( Supplementary Fig. S2). In contrast, the genotype by watering regime interaction represented 1-19% of the total variation and was significant for all traits, except Flw.Avi, DMW. Avi, pH.Avi, and MalicFM.Avi. Interaction partitioning according to method 1 from Muir et al. (1992) indicated that the genotype by watering regime interactions were mainly due to accessions re-ranking across watering regimes (80-100%) and in a minor way to scale changes (0-20%, data not shown). The broad-sense heritabilities ranged from 30% for FructoseFM.Avi.D to 92% for FW.Avi.C. These values were correlated across watering regimes (r H 2 =0.80), as well as the genetic variances (r σ 2 G =0.99), confirming genotype re-ranking across watering regimes ( Fig. 1; Supplementary Table S5).

Impact of the water deficit on fruit quality and yield components
The RIL and GWA populations were grown in Avignon and Agadir in separate greenhouse trials over the years 2013 and 2014, while ensuring similar watering conditions (control and drought) (see Albert et al., 2016 for details concerning the RILs). On average, in both locations, water deficit impacted plant and fruit traits in the same direction in the GWA and RIL populations, with a decline in plant vigor, a decrease in yield, and a higher concentration of the metabolites in fruits (as a percentage of FM) (Table 1). However, when applying the drought treatment, FW.Avi was decreased 2-fold and Nbfruits.Avi 9-fold in the RILs (FW.Avi, -37.7%; Nbfruits, -21.7%) compared with the GWA accessions (FW.Avi, -19.0%; Nbfruits, -2.5%). It resulted in a yield decrease reaching the level of -50% in the RILs against -20% in the GWA accessions. On the other hand, SSC, DMW, and VitCFM were more strongly enhanced in the RILs (SSC, +26.3%; DMW, +30.7%; and VitCFM, +26.3%) than in the GWA accessions (SSC, +12.6%; DMW, +11.4%; and VitCFM, +12.7%).
The correlation between fruit FW in control conditions (indicator of fruit size) and ∆FW was strongly negative in the GWA accessions (Avi, r= -0.55, P=2.70 × 10 -12 ; Aga, r= -0.52, P=2.65 × 10 -10 ), as was previously noted in the RILs. This indicated greater FW loss in larger fruited accessions under drought and increased metabolite contents resulting mainly from the reduced amount of water in the fruits. Thus, the differences observed between the populations may mostly reflect differences in fruit size, with larger fruits among the RILs (8-61 g, mean=20 g, SD=9 g) compared with the GWA accessions (2-46 g, mean=13 g, SD=10 g). Nevertheless, a larger range of variation was observed among the GWA accessions for ∆Yield.Avi and ∆Nbfruits.Avi compared with the RILs (Fig. 2; Supplementary Figs S3, S4). In particular, 55 accessions exhibited an increased yield under drought in the GWA population against only two among the RILs. No noticeable geographic origin or genetic group was obvious among these 55 accessions of the GWA population (10 mixture, 43 SLC, and 2 SP).
When plotting ∆Nbfruits against ∆SSC in regard to fruit size and ∆FW.Avi, the RIL and GWA plants presented different patterns (Fig. 2). Among the RILs, only 18 accessions were present in the top right quarter of the plot corresponding to accessions with increased SSC and Nbfruits under water deficit. Besides, all the top right quarter RILs had a negative ∆FW.Avi (blue and purple color) meaning a decreased FW under drought compared with the control condition for these accessions. On the other hand, 40% of the GWA accessions were present in the top right quarter of the plot and six of them had a positive ∆FW.Avi (magenta and red color) and small to medium fruit size (FW in control from 2 g to 28 g). Similar figures were obtained when considering fruit ascorbate ( Supplementary Fig. S5), malic acid, and citric acid contents ( Supplementary Fig. S6).
When gathering the associations obtained with MLMM and MTMM, 20 associations were detected in common (same trait and same QTL type), resulting in a total of 157 associations for the 27 traits (Supplementary Tables S6-S8). Sixteen associations were detected between twice and three times with related traits ('acid' and 'sugar' traits) and/or for the same trait in the two locations. Thus, a total of 141 different associations were identified, spread unevenly over the genome (Table 2). Chromosomes carried out six (chromosomes 7 and 8) to 23 associations (chromosome 2; Supplementary Fig. S7). Thirty percent of the associations were 'constitutive' (44/141), 30% were 'control specific' (41/141), 22% were 'drought specific' (31/141), and 17% were 'interactive' (25/141). Among the interactive associations, 16 showed 'differential' effects (effect intensity changing according to watering regime) whereas nine presented 'antagonist' effects (effect direction changing according to watering regime). Up to 14, 24, and 28 different associations were mapped for vitamin C, 'acid', and 'sugar' content in fruit, respectively.

Confidence intervals and candidate gene selection under QTLs for fruit quality traits
We observed large differences in size and number of underlying genes when drawing confidence intervals around the association peaks. Eighteen QTLs mapped around the weakly recombinant centromeres covered >10 Mbp and included between 410 and 2573 genes, whereas 84 QTLs covered <5.5 Mbp and encompassed between one and 97 genes ( Supplementary Fig. S9). In the RILs grown in the same conditions , only four QTLs covered <100 genes on a total of 56 QTLs. The comparison of the QTL positions between the RIL and GWA populations resulted in a total of 11 QTLs common to both populations (Table 2), whereas 45 were specific to the RILs and 130 to the GWA population ( Supplementary Fig. S10).
To propose putative candidate genes, we focused on QTLs for vitamin C, sugar, and acid contents in fruit including <100 genes (42 among 66 QTLs) and selected in their intervals genes showing expression in the fruits according to the data from the Tomato Genome Consortium (2012). This reduced the gene list to screen for between one and 87 genes depending on the QTL intervals. Annotations were analyzed to identify genes with functions related to vitamin C, sugar, or acid metabolism under 'constitutive' QTLs and functions related to primary metabolism and/or defense against abiotic stress under 'specific' and 'interactive' QTLs. A total of 41 putative candidates were proposed for three 'constitutive' QTLs (Table 3) and 15 'interactive' or 'specific' QTLs (Table 4). Of those genes, 22 were reported to have DNA polymorphisms in the four accessions of our GWA population which were re-sequenced by Causse et al. (2013). The polymorphisms in four of those genes were predicted to change the amino acids, affecting biological function of a protein.
From the 18 dissected QTLs, 'SSC.Avi_9.1' (control specific) probably corresponded to the cloned QTL 'Brix9.2.5' controlling SSC in fruit and associated with a polymorphism in a cell   (Table 4). A second QTL ('Malic.Avi_6.3') co-localized with a previously mapped QTL for acid content in fruit in different tomato populations and for which two 'aluminumactivated malate transporter-like' genes (Solyc06g072910 and Solyc06g072920) were pointed out as putative candidate genes by Sauvage et al. (2014) (Table 3). Although these two genes presented promising polymorphisms between our four re-sequenced accessions, they displayed a very low expression in fruit (Tomato Genome Consortium, 2012; personal data) and will need further validation to be clearly associated with the phenotypes. Ten QTLs co-localized with loci identified in the RILs (Albert et al., 2016, control and drought conditions) and/or in the three tomato population analyzed by Pascual et al. (2016) (RIL, GWA, and MAGIC, control conditions) but for which no candidate gene was proposed until now, while six were present in genomic regions where, to the best of our knowledge, no QTLs for related traits were mapped thus far. In the intervals of four of them, controlling vitamin C and fructose content in a droughtspecific manner ('VitCDM.Avi_1.1', 'FructoseDM.Avi_4.1', and 'FructoseDM.Avi_10.1'), three genes coding for 'chaperone proteins dnaJ' were identified (Solyc01g105340, Solyc04g009770, and Solyc10g078560; Table 4). Five more genes coding for 'heat/cold shock proteins' (Solyc01g111280, Solyc01g111300, Solyc01g111750, Solyc04g011440, and Solyc04g011450) were identified under antagonist and drought-specific QTLs for fructose and malic acid content ('FructoseDM.Avi_1.1' and 'MalicDM.Avi_4.1'; Table 4).
Three constitutive QTLs, the first two on chromosome 7 controlling glucose and malic acid content and the third on chromosome 10 controlling fructose content, seemed particularly promising. The first two ('GlucoseDM.Avi_7.2' and 'Malic. Avi.7_2' in Table 3) shared a common interval including a gene coding for a 'phosphoenolpyruvate carboxylase' (Solyc07g062530: PEPC) and a gene coding for a 'malate dehydrogenase' (Solyc07g062650). The PEPC gene presented a non-synonymous polymorphism with a predicted impact on the protein function when comparing the four re-sequenced accessions. The third one ('FructoseDM.Avi_10.2' in Table 3) contained two genes coding for 'cell wall invertases', Lin6 (Solyc10g083290) and Lin8 (Solyc10g083300), presenting three non-synonymous polymorphisms between the re-sequenced accessions.

Table 2. Description of QTLs detected for plant and fruit traits in the GWA population through association mapping and comparison with those detected in the RIL population through linkage analysis
QTLs detected in the GWA population were classified according to their type. QTLs significant under both watering regimes are referred to as 'constitutive'. QTLs significant under one watering regime only ('control' or 'drought') are designated as 'specific'. QTLs detected with the plasticity data and/or with the interaction test are designated as 'interactive'. For each phenotypic trait and each QTL type, the number of QTLs, minimum and maximum confidence interval (CI in Mbp on genome assembly v2.5) and minimum and maximum number of genes in the interval are displayed. We considered related traits as a single trait: pH, acid malic (DM and FM), and acid citric (DM and FM) were grouped in 'acids', as well as SSC, glucose (DM and FM), and fructose (DM and FM) in 'sugars'. We gathered QTLs detected in both trial locations (Agadir and Avignon) for the same trait. For the comparison with the RIL population (results described in Albert et al., 2016), whatever the QTL type, we considered that a single QTL was present when the CI overlapped between RIL and GWA QTLs. LG Indication of QTLs confirmed in both locations, Agadir and Avignon (with the same type: 'constitutive', 'specific', or 'interactive').
c Indication of QTLs for acids and sugars confirmed with several measurement methods (pH and acid content, SSC and sugar content). Table 3. Putative candidate genes in the confidence interval around constitutive GWA QTLs for vitamin C, sugar. and acid content in fruit.

QTL(s)
Genes poorly expressed in the fruit.

Table 4. Putative candidate genes in the confidence interval around specific and interactive GWA QTLs for vitamin C, sugar and acid content in fruit
We focused on QTLs encompassing <100 genes.  . Variants which have a deleterious impact on the protein structure according to PROVEAN are indicated by '#'.

Improving fruit quality while maintaining yield in tomato under water limitation
Deficit irrigation strategies aiming to reduce non-beneficial water consumption while maximizing fruit quality and minimizing yield losses are studied in horticultural production to address environmental issues and market expectations simultaneously. It seems particularly relevant for tomato since consumers complain about lack of taste in the new varieties (Bruhn et al., 1991;Causse et al., 2010). In our trials, after a decrease in 60% of the water supply throughout plant growth, we observed on average reduced plant vigor and yield, while fruit quality was improved or stable depending on whether metabolite concentrations were expressed relative to FM or DM. This antagonistic relationship between quality and yield performances confirmed the results obtained in RILs  and the tendencies reported by other authors in tomato (Guichard et al., 2001;Kirda et al., 2004;Zheng et al., 2013), peach (Mirás-Avalos et al., 2013), or grapevine (Santesteban and Royo, 2006). Nevertheless, 50 accessions (with small to medium fruit size) had both improved fruit quality and maintained yield (or even improved) under water deficit compared with the control watering regime, although their vigor (measured through leaf length and stem diameter) was decreased. These accessions emphasized the opportunity to increase metabolite content in tomato fruits using deficit irrigation without achieving parallel limitation of the yield. In contrast, no RIL presented such a response to the water deficit treatment, and the increased sugar and acid contents observed reflected mainly concentration effects due to a decreased amount of water in fruit .
The large phenotypic variations observed mainly resulted from genotype effects (35-80%) and less from genotype by watering regime interactions (1-19%). The watering regime effect represented a significant part of the total phenotypic variability (up to 40%) only for stem diameter and leaf length. This suggests that tomato plants buffer the negative effect of water limitation by limiting their vegetative growth and reallocating the photo-assimilates to the fruits (Lemoine et al., 2013;Osorio et al., 2014).

Benefits and limits of GWA to dissect the genetic architecture of response to water deficit in tomato
Association studies aiming to identify alleles whose effects are modulated by environmental conditions are still few in plants. To date, such studies were only reported in Arabidopsis thaliana Morrison and Linder, 2014;El-Soda et al., 2015;Sasaki et al., 2015), and maize (Saïdou et al., 2014). Explicitly accounting for 'QTL by environment interactions' in QTL studies can help to discover novel genes that act synergistically with the environment, potentially leading to the identification of superior genotypes according to the environments (Des Marais et al., 2013).
We identified a total of 141 QTLs with low to medium effects. The phenotyped traits were strongly polygenic and justified the use of a multilocus GWA mapping model (MLMM: Segura et al., 2012). In particular, up to 14, 24, and 28 different QTLs were identified for vitamin C, acid, and sugar content, respectively. Among the loci identified, 51% were specific to one watering condition, 31% were constitutive and detected whatever the condition, and 18% were interactive between the watering conditions. These proportions of QTL types are relatively similar to those reported in the RILs grown in the same conditions  and in the study of Gur et al. (2011) on tomato introgression lines. However, while most of the interactive QTLs identified in the RILs presented antagonist effects, a majority of differential effects was observed in the GWA study. These discrepancies between both populations may reflect their different genetic basis: the RILs segregate between a small-and a large-fruited accession, whereas the GWA collection focuses on the polymorphisms between several diverse small-fruited accessions.
Because of the large number of markers to be used in GWA analysis, it is not straightforward to choose an appropriate significance threshold controlling for false positives while maintaining the statistical power. We thus opted for a lowered threshold of 10 -4 . If we used Bonferroni correction usually applied to exclude false positives, we should have used a significance threshold of 10 -5 . This would reduce the number of associations detected to 69 (nine 'interactive', 44 'specific', and 16 'constitutive'). With this stringent threshold, we would not have recovered some well-described tomato QTLs, such as, for example, FW11.2 and FW11.3 on chromosome 11 (fruit FW QTLs: Huang and van der Knaap, 2011;Illa-Berenguer et al., 2015). The need for more permissive thresholds in GWASs is often claimed. Strategies based on enrichment tests using known candidate genes from the literature to evaluate the false-positive rate and choose the appropriate threshold values are proposed (Atwell et al., 2010;Sasaki et al., 2015). However, these approaches are limited to well-annotated model genomes and simple traits with already well-described genetic architecture. Another solution to solve the multiple testing issues could be to use haplotypes instead of individual markers to minimize the number of tests, especially in species where the LD spans large genomic regions (Bader, 2001;McClurg et al., 2006). This has already been successfully applied in crops (Gawenda et al., 2015) and would be worth testing in tomato, but may need more markers to identify haplotypes correctly.
The projection of the QTL intervals onto the physical map of tomato allowed the comparison of QTL positions between the RIL and GWA population even though they were genotyped with different markers. This projection resulted in a total of 11 QTLs conserved between both populations. On the other hand, 45 were specific to the RIL population and 130 to the GWA population. This may seem like a relatively small number of common QTLs between the populations, but the RIL parental accessions reflected only a limited fraction of the genetic variation present in the GWA population.

Searching for candidate genes under QTLs for fruit quality traits
Our approach, combining linkage and association mapping, was powerful in recovering previously identified loci associated with fruit quality. As an example, we mapped a QTL associated with fruit fructose content on chromosome 9 which included in its interval the gene Lin5 (Solyc09g010080) known to encode a cell wall invertase affecting tomato fruit sugar content (Fridman et al., 2000). Apart from recovering previously described genes, we identified QTLs in genomic regions where QTLs associated with related traits were previously identified in other populations but for which no candidate gene was proposed until now (probably because of too large confidence intervals) or in genomic regions where, to the best of our knowledge, no QTL was reported for related traits thus far. The confidence intervals around the association peaks obtained using an LD-based approach were mostly shorter (1-97 genes for 84 intervals) compared with the intervals obtained using the RILs or introgression lines (Semel et al., 2007;Gur et al., 2011;Arms et al., 2015).
Combining publicly available expression data (Tomato Genome Consortium, 2012), exonic variants gained from resequencing of four accessions of the GWA collection  and functional analysis of the gene annotations in the confidence intervals, we proposed 41 putative candidate genes under three constitutive QTLs and 15 interactive or specific QTLs. Under the interactive and specific QTLs, genes related to protein protection (chaperone and heat/cold shock proteins), water and solute transport (aquaporins and others transporters), sugar metabolism (sucrose phosphate synthase and invertases), and hormonal signaling (auxin, gibberellin, and ethylene) were identified and may play a crucial role in responses to water deficit (Wang et al., 2003;Shinozaki and Yamaguchi-Shinozaki, 2007). Some of them presented polymorphisms with predicted impacts on the protein function when comparing the re-sequenced accessions and constitute promising targets for future functional validations.
On the bottom of chromosome 7, two QTLs, controlling glucose and malic acid content, shared a common interval including a gene coding for a 'phosphoenolpyruvate carboxylase' (PEPC) and a gene coding for a 'malate dehydrogenase'. The PEPC gene presented a non-synonymous polymorphism with a predicted impact on the protein function in the four resequenced accessions. As the PEPC is catalyzing the carboxylation of the phosphoenolpyruvate arising from glycolysis into oxaloacetate which is then converted into malate by the malate dehydrogenase or enters the Krebs cycle (Guillet et al., 2002), this gene constitutes a likely candidate. Nevertheless, although if the 'malate dehydrogenase' gene did not present any exonic SNPs in our data, it remains an interesting candidate as our four re-sequenced accessions probably did not represent the full genetic diversity present in the GWA population, and the phenotypic variations observed may result from regulation change more than modifications of the protein. On the bottom of chromosome 10, a QTL interval controlling fructose content contained two genes coding for 'cell wall invertases' (Lin6 and Lin8). Both genes presented non-synonymous polymorphisms between the re-sequenced accessions. In contrast to Lin5 on chromosome 9, Lin6 and Lin8 have not yet been associated with variation in sugar content in fruit. Cell wall invertases are extracellular hydrolases which cleave sucrose to glucose and fructose, which are then transported into the cell. They play a central role in regulating, amplifying, and integrating different signals that lead to the source-sink transition in plants.
Subsequent analyses based on either fine mapping around the candidate genes using target re-sequencing approach or functional validation, for example by genome editing, could clarify the involvement of these genes in the phenotypic variations observed.

Supplementary data
Supplementary data are available at JXB online Fig. S1. Structuration observed in the GWA population based on principal co-ordinate analysis (PCoA) on data of 6100 SNPs. Fig. S2. Box-plot of the mean distribution for the nine traits that showed a significant genetic group by watering regime interaction in the ANOVAs. Fig. S3. Distribution of the accession means for plant traits in the GWA population grown under two watering regimes.. Fig. S4. Distribution of the accession means for fruit traits in the GWA population grown under two watering regimes. Fig. S5. Relationship between plasticity of fruit number and plasticity of vitamin C content in fruit, in view of the fruit FW plasticity, in the GWA and RIL populations, respectively. Fig. S6. Relationship between plasticity of fruit number and plasticity of citric and malic acid content in fruit (relative to FW), in view of the fruit FW plasticity, in the GWA population. Fig. S7. Physical map of the QTLs detected in the GWA and RIL populations. Fig. S8. Example of co-localizations between GWA and RIL QTLs for soluble solid content and fruit FW on the bottom of chromosome 11. Fig. S9. Confidence interval (CI) sizes and numbers of genes underlying the QTLs in the GWA and RIL populations. Fig. S10. Venn diagram representing common QTLs between the RIL population (linkage mapping) and the GWA population (association mapping). Table S1. Genetic and phenotypic description of the accessions in the GWA population. Table S2. Genotypic data in the GWA population. Table S3. Principal co-ordinates analysis in the GWA population. Table S4. Correlations between Avignon and Agadir trials. Table S5. Effect of watering regime (W), genetic group (Gr), genotype nested in genetic group [Gr(G)] and the interactions [Gr×W and Gr(G)×W] on the plant and fruit traits measured in the GWA population. Table S6. QTLs identified under both watering regimes ('Control' and 'Drought') using the bivariate multitrait mixed model (MTMM) genome-wide association mapping approach. Table S7. QTLs identified under each watering regime ('Control' and 'Drought') using the univariate multilocus mixed model (MLMM) genome-wide association mapping approach. Table S8. QTLs identified for plasticity data each [(Drought-Control)/Control] using the univariate multilocus mixed model (MLMM) genome-wide association mapping approach.