Phenotyping

The tomato collection was phenotyped for five nutritional and five fruit quality traits. The former group included metabolites with antioxidant activity, such as ascorbic acid (AsA), β-carotene (β-C) cis-lycopene (c-LYC), trans-lycopene (t-LYC) and phenolics (PHE), whereas the latter consisted of dry matter (DMW) and fresh fruit weight (FW), pH, soluble solids content (SSC) and titratable acidity (TA). Detailed information on phenotyping performed for each trait and genotype is reported in Additional file 1.

Heritability values calculated on the two years of phenotypic characterization were higher than 0.5 for all traits except than for cis-lycopene (Table 1). Therefore, phenotypes data were averaged over the two years, and the minimum, maximum and mean values are reported in Table 1, together with the coefficient of variation (% CV). A large range of variation was found for all traits, as also shown in Figure 1. In particular, in the figure is clearly evident that for β-carotene the genotype E71 represents an outlier, since it exhibited a value of 25 μg g⁻¹ FW compared to 1.99 μg g⁻¹ FW mean value of the whole population. Indeed, the genotype E71 corresponds to the variety Caro Red, which was specifically selected for this trait [22]. Consequently, in order to prevent bias, the genotype E71 was excluded from subsequent analyses. As for the other traits, variability estimated by the coefficient of variation ranged from values of approximately 10% to 50%, with only one trait (pH) showing a very low CV value (2.87) and one trait (FW) exhibiting a very high CV value (90.9%).

Trait	H 2	Min	Max	Mean	CV%
Ascorbic Acid (mg 100 g−1 FW)	0.56	22.40	51.23	33.59	17.33
β-carotene (μg g−1 FW)a	0.75	0.11	7.79	1.99	60.02
Trans-lycopene (μg g−1 FW)	0.61	0.19	193.17	85.50	54.97
Cis-lycopene (μg g−1 FW)	0.48	0.00	8.60	3.40	56.13
Phenolics (mg GAE 100 g−1 FW)	0.52	25.06	86.23	48.73	20.47
Dry matter weight (g 100 g−1 FW)	0.53	6.50	10.93	8.56	12.47
Fresh Weight (g)	0.87	5.35	313.50	65.92	90.92
pH (pH units)	0.64	4.08	4.79	4.35	2.87
Soluble solids content (Brix)	0.87	5.13	8.90	6.60	11.54
Titratable acidity (g c.a. 100 mL−1 juice)	0.66	0.27	0.75	0.47	21.54

Heritability (H ²), minimum (min), maximum (max) and mean values, and coefficient of variation (CV%) are shown for each trait.

^aAll data reported for β-carotene are referred to the whole collection except than genotype E71 (see Figure 1).

The Pearson correlation coefficients (r) among traits (Additional file 2) showed a positive value between t-LYC and c-LYC (r = 0.89) and a negative value between pH and TA (r = −0.70). In addition, AsA, PHE and TA were positively correlated with SSC and negatively correlated with FW. PHE content was also negatively correlated with t-LYC and c-LYC (r = −0.38 and −0.49, respectively), whereas it was positively correlated with AsA (r = 0.55).

Genotyping and population structure

Association mapping

AsA content was associated with markers 2383 and 7588, which map on chromosome 3 spanning a region of 150 kbp, and with marker 1241 on chromosome 5. For markers on chromosome 3, genotypes with major alleles showed an increasing AsA level, compared to genotypes with minor alleles (Additional file 9). By contrast, for marker 1241 the minor allele incremented the phenotype. For β-carotene, the analysis revealed significant associations for markers 2022, 2025 and 2028 mapping on chromosome 1. Each markers explained approximately 20% of the phenotypic variation and the minor alleles in all cases contributed to enhance values. Markers 3525 and 3526, co-localized on chromosome 3, and marker 3104 mapping on chromosome 10, were associated with t-LYC with R² values of 0.175 and 0.150, respectively. In all cases, the major alleles showed a very high effect with respect to the corresponding minor alleles. PHE was associated with markers 354 on chromosomes 8 and 4365 on chromosome 11. In both cases, the minor alleles increased the metabolite content.

FW was associated with six markers when the MLM was applied using MAF > 5%: the first was 2992 on chromosome 2, which explained about 17% of the phenotypic variation. Three markers co-segregated on chromosome 8, and mapped in the same gene (Solyc08g006170.1.1). The fifth marker 2275 explained the largest phenotypic variation (R² = 0.239), and mapped 300 bp downstream to Solyc08g006170.1.1. The last marker was 1081 on chromosome 11, and it was the only confirmed using MAF >10%. Moreover, since the multiple regression analysis evidenced a great impact of the genetic structure only on FW, for this trait the association analysis was carried out also on the three separate Q sub-populations. Results confirmed that association of marker 1081 is maintained within the sub-populations (Q1, p-value = 1.90 E-05; Q2, p-value = 4.88E-05; Q3, p-value = 4.7E-02), suggesting that the association of this marker could be considered adequately robust. TA was associated with marker 955 on chromosome 1, which explained about 25% of the phenotypic variation (R²). The other five markers explaining the remaining part of phenotypic variation were marker 2032 mapping on chromosome 2, markers 443 and 3999 co-localized on chromosome 3, and markers 1010 and 1210 on chromosome 4. Finally, only one significant SNP was associated with pH, and explained 16.8% of phenotypic variation. Genotypes exhibiting minor allele for all markers associated with TA and pH had significantly higher value than genotypes with the major allele.

Finally, we evaluated the effect of different allele combinations at loci that were significantly associated with each trait (Figure 2). For each trait, mean and statistical significance among the groups of genotypes were calculated for all the allelic combinations. For AsA, four allele combinations were found. Group 1 showed the highest value (35.03 mg 100 g⁻¹ FW, average of 57 genotypes) and group 4 the lowest (27.85 mg 100 g⁻¹ FW, average of seven genotypes). For β-C, 46 genotypes in group 1 and six in group 2 had allele combinations associated with a low content and 30 genotypes with a high content. t-LYC showed three allele combinations. Four genotypes with yellow fruits belong to group 1 associated with the lowest lycopene content (4.25 μg g⁻¹ FW mean value), whereas 67 genotypes showed an allele combination associated with high lycopene content (95.63 μg g⁻¹ FW mean value). For PHE, four groups were observed. The largest was group 1, including 49 genotypes and showing the minimum amount of phenolics (44.27 mg GAE 100 g⁻¹ FW mean value), while the group associated to the maximum amount (66.38 mg GAE 100 g⁻¹ FW mean value) included eight genotypes. Eleven allele combinations were identified for titratable acidity and the one associated with the highest value (0.745 g citric ac. 100 mL⁻¹ of juice) was detected in two genotypes, whereas that associated with the lowest value (0.418 g citric ac. 100 mL⁻¹ of juice) was detected for a group of 48 genotypes. Intermediate values were detected for the other nine groups.

Nutritional traits

Concerning antioxidants traits, markers associated with AsA, β-C, t-LYC and PHE were searched for, since these are bioactive compounds exhibiting beneficial effects on human health [29]. In particular, three markers (2383, 7588 and 1241) associated with AsA were identified, which differed from those detected by Sauvage et al. [20] using a similar GWAS approach, but exploiting accessions belonging to different tomato species. Two markers we identified corresponded to genes Solyc03g112630.2.1 and Solyc03g112670.2.1 mapping on chromosome 3 and were annotated as Fas-associated factor 1-like and Genomic DNA chromosome 5 P1, respectively. The other gene (Solyc05g052410.1.1) was located on chromosome 5 and annotated as Ethylene-responsive transcription factor 1 (ERF1). The Fas-associated factor 1-like protein is involved in an apoplastic mechanism and no direct evidence was reported to correlate its function with AsA accumulation. Since no specific functions were also assigned to Solyc03g112670.2.1, it was thought that the polymorphisms identified in this region of chromosome 3 could be in LD with other candidate genes. In order to verify this hypothesis, a scan was performed of the surrounding genomic area in LD with markers 2383 and 7588. A cluster of pectinesterases (120 kbp from marker 7588), one pectate lyase (240 kbp from marker 7588) and one polygalacturonase (350 kbp from marker 7588) were detected in LD block 23 on chromosome 3. These findings suggest that the alternative D-galacturonic biosynthetic pathway could contributes to regulate AsA variation in the tomato population under study, as previously reported in tomato [30] and other species [31],[32]. In addition, concerning the ERF1 gene, Di Matteo and colleagues [30] showed that in one S. pennellii introgression line a different expression of genes associated with ethylene biosynthesis might trigger pectin degradation resulting in AsA accumulation. Taken together, these results suggest a possible regulation of genes associated with markers 2383 and 7588 (related to pectin degradation) via Ethylene Responsive Factor 1 associated with marker 1241.

Cis and trans isomers of lycopene derive from a cascade of enzymatic reactions taking place in plastids [33]. Intermediates in the first part of the pathway are cis-configured. A pro-lycopene isomerase (CrtISO) then produces all-trans-lycopenes from tetra-cis-lycopene. Subsequent reactions convert trans-lycopene into β-carotene by the action of a lycopene β-cyclase (β-Lcy). Although no associations were detected for c-LYC, three significant associations with t-LYC were identified. Markers 3525 and 3526, co-localized on chromosome 3, matched a putative metallocarboxypeptidase inhibitor (Solyc03g031480.2.1) and tyrosyl-DNA phosphodiesterase (Solyc03g031820.2.1) whereas marker 3104 on chromosome 10 did not match annotated genes. Interestingly, even if they are not directly linked to the trans-lycopene content, analysis of the genomic area highlighted the presence of a phytoene synthase 1 (Solyc03g031860.2.1) close to markers 3525 and 3526 (at 315 and 24 kbp, respectively) and a lycopene β-cyclase 2 (Solyc10g079480.1.1) at 9 kbp from marker 3104. This showed that the association mapping approach used was able to validate two candidate genes already known to be involved in the carotenoid pathway. In fact, the identified phytoene synthase 1, which catalyzes a rate-limiting step in the carotenoids pathway, corresponds to the locus “r” [34] that carries a recessive mutation conferring a characteristic yellow flesh phenotype. Four accessions in the population showed a genotype associated to the locus “r” and all have yellow flesh fruit as a consequence of a low trans-lycopene content. In addition, on chromosome 10, besides the lycopene β-cyclase, also a carotene isomerase (CrtISO, Solyc10g081650.1.1), which converts pro-lycopene to trans-lycopene [35], was localized 1.2 Mbp downstream marker 3104 (Figure 3). As concerns β-C, significant associations were found on chromosome 1 with Solyc01g087600.2.1 annotated as Protein E03H4.4, Solyc01g087670.2.1 annotated as a guanine nucleotide-binding protein, involved in blue light perception signal pathways [36], and Solyc01g087880.2.1, which has no homology with any gene of known function. These results prompted to investigate alternative genes in this region. Scanning the genomic area associated to these three markers (LD block 13 on chromosome 1), a 24-C-sterol-methyltransferase was found that is involved in steroid biosynthesis. Moreover, a cluster of three carotenoid cleavage dioxygenase 1 (CCD1) genes was also identified at 300 kbp from marker 2022. CCD1 genes cleave the carotenoid substrate at different double bonds to produce terpenoid flavour volatiles (apocarotenoids) that contribute to the overall aroma and taste of tomato fruit [37]. It is hypothesized that the variation in the carotenoid pool may depend on the metabolic flux towards the cleavage reactions to produce apocarotenoids, but further functional experiments that will validate this hypothesis are required.

Finally, two polymorphic markers significantly associated with PHE were identified: marker 354 (Solyc08g082350.2.1), which encodes for a protein of unknown function, and marker 4365 (Solyc11g010170.1.1), encoding for a LanC-like protein2, which is involved in the modification and transport of peptides in bacteria [38]. However, no well-defined functions were reported for the latter gene in plants [39], even if a probable involvement as one receptor for abscisic acid (ABA) was hypothesized [40]. No gene of the phenolics pathways was detected in the putative region in association with marker 354. By contrast, significant co-localizations (LD block 6) found with marker 4365 included transport genes encoding two 14-3-3 proteins (Solyc11g010200 and Solyc11g010470), an ABC-2 transporter (Solyc11g009100) and a MATE efflux family protein (Solyc11g010380). The involvement of these transporters in enhancing the vacuolar compartmentalization of phenolic compounds was previously reported by Gomez et al. [41] and Di Matteo et al. [42], suggesting their probable role in the metabolism of this trait. A previous work [43] also identified QTLs for phenolic content in regions of chromosome 8 and 11 close to the markers detected, confirming the involvement of these regions in phenolics control.

Quality traits

Major fruit quality traits of interest for both the fresh market and processing tomatoes include fruit size, shape, total solids, colour, firmness, ripening, pH, titratable acidity, soluble solids content and dry matter. In this study, a large number of associations with FW and TA were found, only one association with pH and no association with SSC and DWM.

FW is a quantitatively inherited trait controlled by up to 28 QTLs, even though QTL analyses in previous studies revealed that most (67%) phenotypic variation in fruit size could be attributed to six major loci (fw1.1, fw1.2, fw2.1, fw2.2, fw3.2 and fw11.3) localized on chromosomes 1, 2, 3 and 11 [44]-[47]. The present study confirmed only one of the above loci (fw11.3).

Indeed, on chromosome 11 marker 1081 matched Solyc11g071840.1.1, annotated as a calmodulin binding protein, and was located in the LD block 40 that spans an interval of 167 kbp. This region contains both a portion of the fw11.3 locus (starting 20 kbp downstream of marker 1081) and the fas-YABBY locus (24 kbp upstream of the marker), previously hypothesized to determine fruit size [48]. Therefore, the findings here reported not only confirmed the involvement of locus fw11.3 in FW variation but also restricted the region of putative candidate genes with respect to the previously identified region of 149 kbp, which included 22 predicted genes [48]. Indeed, considering the LD block strategy used, marker 1081 was strongly associated only with a portion of fw11.3 of around 70 kbp (52,301,894 – 52,365,467), including eight predicted genes. This region must be enriched in SNPs to locate precisely one or more responsible polymorphisms associated to the trait and further investigation should be carried out to fine map the putative gene responsible for the trait variation.

Among important quality traits in tomato, TA influences shelf-life of processed tomato and low pH values reduce the risk of pathogen growth in tomato products [49]. One locus associated with pH and five loci associated with TA were evidenced. Only marker 2246 that matches Solyc11g017070.1.1 on chromosome 11 (encoding for an eukaryotic translation initiation factor 3 subunit 2) proved associated with pH trait. An ATPase and a Sucrose_H+ symporter co-localized with this marker; they are proton pumps responsible for acidifying cellular compartment and might correlate with this trait. Moreover, previous studies identified a pH QTL in this region [44],[50] supporting our hypothesis and confirming the possible involvement of this region in regulating the pH level in tomato fruit.

The most significant association with TA was with marker 955, which explains around 25% of the variation. The marker matched the predicted gene Solyc01g107550.2.1, which encodes for a methylthioribose kinase, an enzyme involved in recycling of methionine through the methylthioadenosine (MTA) cycle. The marker did not fall in any significant LD block on chromosome 1 and for this reason a narrow area around the marker was investigated to look for candidate genes. Four annotated genes not related to the trait under study and various unknown genes were identified, leaving the mechanism responsible for TA variation of this locus still obscure. Other significant associations were found on chromosome 3 for markers 443 and 3999 matching Solyc03g083440.2.1 (Glutamate synthase) and Solyc03g093310.2.1 (F-box family protein), respectively, in a region where previous studies detected QTLs for TA [44],[50],[51]. Markers 443 and 3999 are 1 Mbp away from one another and do not fall in the same LD block. Marker 443 falls in the 280 kbp LD block 15 that includes at least 40 genes, none of which might be so far involved in TA determination. On the other hand, for marker 3999 no significant LD block was inferred and few putative co-localizations could be highlighted. In particular, an aquaporin (Solyc03g093230.2.1) was identified as a potential candidate. The involvement of the aquaporin gene family in the modulation of fruit acidity was hypothesized in tomato antisense experiments, indicating a strong effect of this protein on the sugar/organic acid ratio in fruit [52]. An additional locus for TA was found on chromosome 4 between markers 1210 and 1010. This region includes about 30 genes and we pointed out a purple acid phosphatase (Solyc04g005450), a metalloenzyme that hydrolyses phosphate esters and anhydrides under acidic conditions [53]. No relevant co-localizations were instead found for marker 2032 on chromosome 2.

Plant material

Plant material consisted of 96 tomato genotypes, including Italian and Latin American landraces, and vintage and modern varieties collected from seed banks in Italy and worldwide. In detail, accessions were derived from Italian breeders’ collections, Regional and National Italian Institutions (Regione Campania, ARCA2010 Cooperative for Agriculture, MiPAF-CRA Centre for Research in Agriculture, University of Naples Federico II), and International Institutions (Plant Genetic Resources Unit, USDA, Tomato Genetics Resource Center, Davis, USA, Hebrew University of Jerusalem, Israel). All genotypes were grown during the seasons 2011 and 2012, according to a randomized complete block design with three replicates (10 plants per replicate), in field plots at the Agricultural Experiment Station of CRA-ORT in Battipaglia (Salerno, Italy). Field trials were conducted in accordance with the Italian legislation. All genotypes were subject to phenotypic and genotypic analyses. Out of 96 genotypes, six were excluded from phenotyping, since they produced few fruits or segregated for some morphological traits. Five additional genotypes were filtered out since were considered unreliable due to the large number of low quality SNP scores. Therefore, the number of samples finally used for association mapping was reduced to 85.

Phenotyping

Chemical and physical traits were evaluated on ten fruits harvested at red ripe stage for each biological replicate as recommended in the SCAR Agro-Food Tomato Working Group [54]. Traits included fresh weight (FW), total dry matter weight (DMW), soluble solids content (SSC), titratable acidity (TA) and pH.

Metabolic analyses were carried out on collected fruits stored at −80°C (three replicates of 5–8 fruits per accession): ascorbic acid (AsA), β-carotene, trans- and cis-lycopene, and phenolics content were measured as described below. AsA content was determined as reported by Di Matteo et al. [30] with minor modifications. Briefly, 500 mg of frozen powder were added to 300 μL of ice-cold 6% trichloroacetic acid (TCA) in 2 mL Eppendorf tubes. Samples were vortexed and left on ice for 15 min, and then centrifuged for 15 min at 25,000 g at 4°C. The supernatant was transferred to a clean tube and 20 μL were used for the assay, as described in the manuscript cited above. The AsA content was expressed as mg per 100 g⁻¹ of FW.

Extraction and analysis of carotenoids were carried out on 5 g of fruit pericarp, according to Ishida et al. [55]. Reversed Phase-High Performance Liquid Chromatography (RP-HPLC) analysis was performed through a Waters E-Alliance HPLC system constituted by a 2695 separations module with quaternary pump, autosampler, and a 2996 photodiode array detector; data were acquired and analyzed with Waters Empower software. The chromatographic separations were performed at a flow rate of 0.8 mL min⁻¹ and at 0.005 AUFS (Absorbance Units Full Scale) by using a reversed phase, analytical polymeric C30 column (250 × 4.6 mm i.d.; 3 μm particle diameter; YMC, Wilmington, NC, USA). Results were expressed as μg g⁻¹ of FW. Solvents used for sample preparation and extractions were of analytical grade, while those for HPLC analysis (methyl-t-butyl ether, methanol, ethyl acetate and tetrahydrofuran) were of HPLC grade; all were obtained from Merck (Darmstadt, Germany). Trans-carotenoid standards (lycopene and β-carotene) used in HPLC analyses were purchased from Sigma Chemical Co (Sigma-Aldrich Company, St. Louis, MO, USA).

Total phenolics content was assayed using a modified procedure of the Folin–Ciocalteu test [56]. In brief, 250 mg of frozen ground tissue were homogenized in a mortar with pestle and extracted using 1 mL of 60% methanol. Samples were left on ice for 3 min in the dark. Crude extracts were transferred into a 15 mL tube and volume was increased to 5 mL by adding 60% methanol. The samples were centrifuged at 3000 g for 5 min; afterwards, 62.5 μL of the supernatant, 62.5 μL of Folin–Ciocalteu’s reagent (Sigma) and 250 μl of deionized water were mixed and incubated for 6 min; 625 μL of 7.5% sodium carbonate and 500 μL of deionized water were added to the samples and incubated for 90 min at room temperature in the dark. Absorbance was measured at 760 nm. The concentration of total phenolics was expressed in terms of mg of gallic acid equivalents (GAE) per 100 g⁻¹ of FW.

For each year, the normal distribution of data was verified using a Shapiro-Wilk test. Six (β-C, PHE, FW, pH, SSC, TA) out of 10 traits were not normally distributed and were log₁₀ transformed before performing the association mapping analysis. In order to test the consistency of phenotypic characterization over years, heritability values were calculated as reported in [15]. Moreover, year and genotype effects were assessed by two-factor analysis of variance. Since the genetic effect over the two years was more significant than the genotype x year interaction for all traits, associations were calculated by using means over years. Pearson correlation coefficients were calculated among all pairs of traits. Analyses were carried out with the R program [57].

SNP genotyping

Samples were genotyped using a tomato array built in the framework of the Solanaceae Coordinated Agricultural Project (SolCAP) from NIFA/USDA and based on the ILLUMINA Infinium Technology. The SolCAP tomato panel includes 7,720 markers constructed on eSNPs deriving from six tomato genome sequences. Details of the SolCAP SNP discovery pipeline are described in Hamilton et al. [2] and Sim et al. [6]. Information about SolCAP SNPs are available in the SGN database (http://solgenomics.net/).

For each accession, genomic DNA was extracted from fresh, young leaf tissue using the DNeasy Plant Mini kit (QIAGEN, Valencia, CA) according to the manufacturer’s recommendations. DNA quality and concentration were evaluated on agarose gel and spectrophotometrically using the Nanodrop instrument (Thermo Fisher Scientific, Wilmington, USA). Genotyping was conducted at the Genomix4Life S.r.l (http://www.genomix4life.com) using 250 ng of DNA per accession following the manufacturer’s protocol for the Illumina Infinium assay. Intensity data and SNP calls were performed by GenomeStudio version 1.7.4 (Illumina Inc., San Diego, CA, USA). SNPs were called using the Infinium chip cluster file based on the SolCAP tomato collection and a manual classification was implemented when the default clustering was not clearly defined. In addition, quality and reproducibility were tested using duplicated DNA samples.

Data analysis

The set of SNPs was filtered in order to perform molecular analyses. Markers with more than 10% missing genotypes and with minor allele frequency (MAF) <5% were removed. The physical position of all SNPs on the 12 tomato chromosomes was obtained from the SGN database (http://solgenomics.net).

Linkage disequilibrium (LD) values were calculated on SNP set with MAF greater than 5%. Pairwise r ² between markers was calculated for each chromosome using TASSEL v.4.0 [58]. Values of r ² were plotted against physical distance and LD decay was inferred via locally weighted scatterplot smoothing (LOWESS) using R software and testing smoothing parameters fixed to 0.2. In order to determine the distance of LD decay, for each chromosome different r ² baseline values (0.2, 0.3 and 0.5) were tested [59]. A Bayesian population classification was carried out using STRUCTURE 2.3.3 [60]. Since population structure estimates assume unlinked markers, a sub-set of 600 markers with an average distance of 600 kbp was chosen to perform STRUCTURE analysis. STRUCTURE runs were carried out with a length of burn-in and MCMC (Markov chain Monte Carlo) of 100,000 each. Twelve independent runs were conducted allowing K (number of populations) varying from 2 to 13. Optimal K was inferred by using the Evanno et al. [61] transformation method. The influence of the population structure on the phenotypic variation of considered traits was assessed by multiple regression analysis. Associations between genotypes and phenotypes were determined using the General Linear Model (GLM) and the Mixed Linear Model (MLM) in TASSEL v4.0. For both models the structured association model (Q model) was performed [62] and two different thresholds of MAF (>5% and >10%) were used. In the association analysis, we considered the kinship matrix based on the SNP data in the model of MLM, and the population structure covariates detected in the tomato accessions with the STRUCTURE analysis. The significance of association between traits and markers was estimated by using an adjusted P value (Bonferroni correction) and the threshold for the association was set to <5.88×10⁻⁴ (0.05/85). In addition, in order to estimate LD blocks according to the definition of Gabriel et al. [63], the HAPLOVIEW software [64] was used with the following parameters: MAF >0.05; Hardy–Weinberg P-value cut-off, 0; percentage of genotyped lines >0.75.

Data availability

The list of genotypes including details on source and distribution, and the SNP genotype dataset are deposited on the LabArchives repository hosted http://dx.doi.org/10.6070/H4TT4NXN.