Deciphering genetic diversity and inheritance of tomato fruit weight and composition through a systems biology approach.

Integrative systems biology proposes new approaches to decipher the variation of phenotypic traits. In an effort to link the genetic variation and the physiological and molecular bases of fruit composition, the proteome (424 protein spots), metabolome (26 compounds), enzymatic profile (26 enzymes), and phenotypes of eight tomato accessions, covering the genetic diversity of the species, and four of their F1 hybrids, were characterized at two fruit developmental stages (cell expansion and orange-red). The contents of metabolites varied among the genetic backgrounds, while enzyme profiles were less variable, particularly at the cell expansion stage. Frequent genotype by stage interactions suggested that the trends observed for one accession at a physiological level may change in another accession. In agreement with this, the inheritance modes varied between crosses and stages. Although additivity was predominant, 40% of the traits were non-additively inherited. Relationships among traits revealed associations between different levels of expression and provided information on several key proteins. Notably, the role of frucktokinase, invertase, and cysteine synthase in the variation of metabolites was highlighted. Several stress-related proteins also appeared related to fruit weight differences. These key proteins might be targets for improving metabolite contents of the fruit. This systems biology approach provides better understanding of networks controlling the genetic variation of tomato fruit composition. In addition, the wide data sets generated provide an ideal framework to develop innovative integrated hypothesis and will be highly valuable for the research community.

Introduction

Identifying the genes controlling the variation of complex traits is a key goal of evolutionary genetics and plant biology. Attempts to identify genetic variants underlying quantitative traits have been achieved by traditional linkage mapping and genome-wide association studies using molecular markers. However, resolution of quantitative trait loci (QTLs) is limited and the identification of the polymorphisms responsible for the variation is not straightforward. Furthermore, as several intermediate levels interact between the genotypes and the phenotypes, DNA sequence variation [single nucleotide polymorphisms (SNPs) or Indels] may not directly affect the traits. Intermediate molecular phenotypes such as gene expression, protein abundance, and metabolite concentration also vary in populations and are themselves quantitatively inherited (Rockman and Kruglyak, 2006). Nowadays, rapid technological advances in high-density experiments such as next-generation sequencing (NGS), RNA expression analysis through microarray or RNAseq, mass spectrometry (MS) coupled to gas chromatography (GC-MS) or to liquid chromatography (LC-MS), and nuclear magnetic resonance (NMR) metabolic profiling enable scientists to obtain large exhaustive data sets and analyse biological systems as a whole. Integration of the genome expression products at different levels should help in dissecting the genetic variation of a given quantitative trait.

Systems biology relates the variation analysed at different expression levels, from phenotype to metabolome and proteome, studying the behaviour of all the elements in a biological system (Gutierez ; Saito and Matsuda, 2010). A bottom-up systems biology approach consists of integrating ‘omic’ resources (genomic, transcriptomic, proteomic, and metabolomic) and large physiological data sets, together with statistical network analysis in order to identify candidate genes underlying phenotypes and to construct complex regulation networks (Kliebenstein, 2010). This approach was first applied to yeast by combining DNA microarrays and quantitative proteomics to describe the galactose pathway (Ideker ). It was then applied to gene expression analysis in Escherichia coli (Rosenfeld ), and to Arabidopsis by Hirai who elucidated gene to gene and metabolite to gene networks by integrating metabolomic and transcriptomic data. Systems biology has also been used to study the natural genetic variation at different levels, such as metabolomics (Keurentjes, 2009; Kliebenstein, 2009), proteomics (Stylianou ), and transcriptomics (Keurentjes ; Kliebenstein, 2009).

Tomato (Solanum lycopersicum) is the model species for the study of fleshy fruit development and composition (Giovannoni ). It is a self-pollinated species and is derived from its closest wild ancestor Solanum pimpinellifolium (Nesbitt and Tanksley, 2002). Cherry tomato accessions (Solanum lycopersicum var. cerasiforme) have an intermediate position between these two species, as their genome is a mosaic of those from S. lycopersicum and S. pimpinellifolium (Ranc ). During tomato domestication, the diversification of fruit aspect, as well as the adaptation to a wide range of environmental conditions was simultaneous with a strong reduction of molecular diversity (Miller and Tanksley, 1990; Blanca ). This lack of genetic variation in cultivated species led geneticists to study trait variation mostly in distant crosses involving wild species, and thus limited the exploitation of intraspecific variation. Today the availability of the tomato genome sequence (Tomato Genome Consortium, 2012) and of a large number of SNP markers (Sim ) allows a re-examination of the variation and inheritance of agronomical and fruit traits at the intraspecific level.

Systems biology approaches have been used in tomato to study fruit development. Carrari and Mounet analysed transcriptome and metabolome variation during fruit development. Garcia combined phenotype, metabolome, transcriptome, and proteome profiles to study genes related to ascorbate metabolism in three transgenic lines. Wang compared the transcriptome and metabolome to uncover the molecular events underlying fruit set, while Osorio compared enzyme activity, metabolite and transcript profiles to analyse the connectivity between these groups of traits in fruit ripening mutants. However, these studies were only focused on a few mutants or on the effect of introgression in S. lycopersicum of wild species alleles. Little is known about the genetic variation in metabolic, enzymatic, and proteomic profiles contributing to phenotypic trait variation in the species.

In the present study, the aim was to decipher the complex relationships between several successive levels of omic profiles to characterize the genetic variation and physiological bases of quantitative traits in tomato fruit. For this purpose, the variation of eight genotypes representing a large range of phenotypic and genotypic diversity (four S. lycopersicum and four S. lycopersicum var. cerasiforme) and four of their corresponding hybrids was compared at two stages of fruit development [cell expansion (CE) and orange-red (OR)]. Their metabolic, enzymatic, and proteome profiles were characterized. Genetic variability was analysed for all traits, and inheritance patterns of traits that were significantly different among genotypes were assessed. Relationships among traits were analysed within and between each group of traits at each stage, and networks were constructed using sparse partial least square (sPLS) regression. This systems biology approach combining proteome, metabolome, and phenotypic analysis gave insights into the diversity and relationships of quantitative traits at different levels.

Materials and methods

Plant materials

Eight tomato lines comprising four S. lycopersicum accessions (Levovil, Stupicke Polni Rane, LA0147, and Ferum) and four S. lycopersicum lycopersicum var. cerasiforme accessions (Cervil, Criollo, Plovdiv24A, and LA1420), and four of their corresponding F1 hybrids (Levovil×Cevil, Stupicke Polni Rane×Criollo, LA0147×Plovdiv24A, and Ferum×LA1420) were used in this study (details of the accessions are shown in Supplementary Table S1 and Supplementary Fig. S1 available at JXB online). Lines were selected, based on a previous molecular characterization of 360 tomato accessions (Ranc ), to include the maximum genetic diversity of the species. Xu ) genotyped these lines and the line sequenced to obtain the reference genome (Heinz1706, Tomato Genome Consortium, 2012). From the 139 SNP markers characterized, 133 were polymorphic (96%), showing the large range of molecular diversity represented by the eight lines. The range of polymorphism between the lines and the reference genome (Heinz1706) ranged from 27% to 82% (Supplementary Table S2) The genetic distances among the parents of F1 hybrids were variable. According to data of Xu ), Levovil and Cervil were the two most distant accessions (82% SNP polymorphic), followed by LA0147×Plovdiv 24A (40%), Stupicke Polni Rane×Criollo (34%), and Ferum×LA1420 (27%).

Plants were grown during spring 2010 under greenhouse conditions (16/20 °C) in Avignon (South of France). Plants were separated into two blocks; five plants per genotype were included in each block.

For proteome, metabolome, and enzymatic measurement, two stages of development, CE and OR, were selected. The CE stage was chosen as a representative stage of the growing tomato fruit; the OR stage was chosen because it is unequivocally determined and is the key step where enzyme and protein concentrations are changing and will determine the final characteristics of the fruit. The number of days after anthesis to reach cell expansion varied among genotypes depending on their fruit size. Thus CE sampling was done at 14, 20, or 25 d after anthesis for small (Cervil), medium (Criollo, Plovdiv 24A, Stupicke Polni Rane, and the four F1 hybrids), or large (LA0147, Levovil, and Ferum) fruited accessions, respectively. OR sampling was done based on fruit colour change. Three biological replicates by stage were analysed. Each replicate included 7–20 fruits from both greenhouse blocks to buffer environmental variations. Fruit pericarps were collected, immediately frozen, ground in liquid nitrogen, and stored at –80 °C until analysis. For fruit phenotypic trait measurements, five fruits were harvested from the 10 plants of each genotype at the following six stages: (i) CE stage; (ii) CE+7 d; (iii) CE+14 d; (iv) CE+21 d; (v) OR stage; and (vi) red ripe. Fruits were evaluated for fresh weight (FW), fruit diameter (FD; measured using a caliper), and dry matter content (DMC). DMC (expressed in g 100g–1 FW) was assessed after 5 d in a ventilated oven at 80 °C.

Metabolome and enzyme activity analysis

Metabolome analyses were performed at the Metabolome Facility of Bordeaux, using quantitative proton NMR (1H-NMR) profiling of polar extracts and liquid chromatography quadrupole time-of-flight tandem mass spectrometry (LC-QTOF-MS) profiling of semi-polar extracts. 1H-NMR profiling was performed as described in Deborde with minor modifications. Briefly, polar metabolites were extracted on lyophilized powder (50mg DW per biological replicate) with an ethanol–water series at 80 °C. The lyophilized extracts were titrated to pH 6 and lyophilized again. Each dried titrated extract was solubilized in 0.5ml of D2O with (trimethylsilyl) propionic-2,2,3,3-d4 acid (TSP) sodium salt (0.01% final concentration) for chemical shift calibration and EDTA (5mM final concentration for the CE and 2mM for the OR stage). 1H-NMR spectra were recorded at 500.162 MHz on a Bruker Avance III spectrometer (Bruker, Karlsruhe, Germany) using an ATMA inverse 5mm probe flushed with nitrogen gas and an electronic reference for quantification (ERETIC2). Sixty-four scans of 32 000 data points each were acquired with a 90 ° pulse angle, a 6000 Hz spectral width, a 2.73 s acquisition time, and a 25 s recycle delay. Two technological replicates were used per biological replicate. Preliminary data processing was conducted with TOPSPIN 3.0 software (Bruker Biospin, Wissembourg, France). The assignments of metabolites in the 1H-NMR spectra were made by comparing the proton chemical shifts with values of the MeRy-B metabolomic database (Ferry-Dumazet ), by comparison with spectra of authentic compounds recorded under the same solvent conditions, and/or by spiking the samples. The metabolite concentrations were calculated using AMIX (version 3.9.7, Bruker) software. The 144 1H-NMR spectra of the data set were converted into JCAMP-DX format and deposited with associated metadata into the Metabolomics Repository of Bordeaux MeRy-B (Ferry-Dumazet , http://www.cbib.u-bordeaux2.fr/MERYB/view/project/34, accessed September 2013).

LC-QTOF-MS profiling of aqueous methanol–0.1% formic acid extracts was performed from lyophilized powder (20mg in 1ml). For each biological replicate, two extractions were performed and two injections per extract were used. An Ultimate 3000 HPLC (Dionex, Sunnyvale, CA, USA) was used to separate metabolites on a reversed phase C18 column (150×2.0mm, 3 µm; Phenomenex, Torrance, CA, USA) using a 30min linear gradient from 3% to 95% acetonitrile in water acidified with 0.1% formic acid. Metabolites were detected using a QTOF mass spectrometer (Bruker). Electrospray ionization in positive mode was used to ionize the compounds. Scan rate for ions at m/z range 100–1500 was fixed at two spectra per second. Methyl vanillate was spiked in the extraction solvent and used as an internal standard. One sample was used as a quality control sample and injected after every 10 injections. Raw data were processed in a targeted manner using QuantAnalysis 2.0 software (Bruker). This resulted in eight compounds identified based on accurate mass measurement and comparison with data from Gomez-Romero .

Ascorbic acid was measured using a spectrofluorometric method and values were expressed as total ascorbate (ascorbic acid+dehydroascorbate) as previously described by Stevens . The maximum activity (V max) of 26 enzymes of primary metabolism was assayed using a robotized platform as described in Gibon and in Steinhauser . Supplementary Tables S3 and S4 at JXB online present the lists of the primary and secondary metabolites and the enzyme activities analysed.

Proteome analysis

Methods for protein extraction, two-dimensional gel electrophoresis (2-DE), and protein identification and classification were as detailed in Xu ). Briefly, proteins were extracted using the phenol extraction method developed by Faurobert . Later, proteins were separated by 2-DE. After Coomassie colloidal staining, image analysis was performed with Samespot software (version 4.1), and the normalized spot volumes were obtained. Protein identification of 424 variable spots was performed at the proteome platform of Le Moulon (Gif-sur-Yvette) using a nano-LC-MS/MS method following the procedure described in Xu ). The database search was run against the International Tomato Annotation Group (ITAG) Release 2.3 of predicted proteins (SL2.40) database (http://solgenomics.net/, accessed September 2013) with X!Tandem software (http://www.thegpm.org/TANDEM/, version 2010.12.01.1). The Fasta sequence of the identified proteins was employed to re-annotate the proteins using the Blast2GO package (Conesa ). Sequences were compared against the NCBI-nr (version April 9, 2012) database of non-redundant protein sequence using BLASTX with the default settings.

Statistical analysis and inheritance analysis

Metabolite contents and enzyme activities were expressed on a dry weight basis to be comparable. All the analyses were performed using R (R Development Core Team, 2012, http://www.R-project.org/, accessed September 2013).

Data were submitted to a two-way analysis of variance (ANOVA; P < 0.05) with genotype, stage, and interaction effect, and then to one-way ANOVA with genotype effect at each stage. In addition, to assess the mode of inheritance of the traits, one-way ANOVA was also performed with genotype effect for each cross (two parental lines and their hybrid) and stage.

Means and standard deviations were calculated for each trait (phenotypic traits, metabolite contents, enzyme activities, and protein spot volumes) in each genotype and stage. Significantly different (P < 0.05) traits in each cross were selected at each stage to estimate additive (A) and dominance (D) components of genetic variation. A is equivalent to half of the difference between two parental lines. The S. lycopersicum line was systematically the first parent in a cross. D is the difference between the hybrid value and the parental mean. The inheritance pattern of each trait was then assessed by the dominance/additivity (D/A) ratio and classified as over-recessive (OR; D/A < –1.2), recessive (R; –1.2 ≤ D/A ≤ –0.8), additive (A; –0.8 < D/A < 0.8), dominant (D; 0.8 ≤ D/A ≤ 1.2), or overdominant (OD; D/A >1.2).

Means of metabolite contents, enzyme activities, and protein spot volumes were centred and scaled to variance unit and used for the rest of the analysis. Principal component analyses (PCAs) were performed for metabolites, enzymes, and phenotypic traits, as well as for protein spot volumes for both development stages and at each stage, with the ‘pcaMethods’ package (Stacklies ). Pearson correlations and P-values were calculated between significantly variable traits at each stage. Correlations were considered to be significant when |r| >0.7 (P-value <0.01). Significant correlations were plotted using the R ‘corrplot’ package (Wei, 2012, http://CRAN.R-project.org/package=corrplot, accessed September 2013). To analyse the relationships among protein spot volumes and other traits, networks were reconstructed and visualized using sPLS correlation regression analysis with the ‘mixOmics’ package (Lé Cao ). An arbitrary threshold of 0.7 was employed for network reconstruction. Nodes represent the different traits, and edges represent the relationships between variables belonging to different levels.

Results

To represent a large range of the genetic diversity, eight tomato accessions were chosen according to previous studies (Ranc ; Xu ). The eight accessions and their four hybrids were characterized at phenotypic, metabolic, and proteomic levels. The final fruit weight of the eight parental lines and the four hybrids ranged from 5.3g to 134.4g. Fruit weights of the four hybrids were intermediate between the values of their parental lines throughout fruit development (Supplementary Fig. S2 at JXB online). Fruit diameter (Fig. 1) was highly correlated to fruit weight, as fruits were round. DMC also showed a wide range of variation (Fig. 1; Supplementary Table S5).

Fig. 1.

Fruit size (FD) and dry matter content (DMC) of the eight tomato accessions studied and their four hybrids (measured at the orange-red stage). (This figure is available in colour at JXB online.)

Metabolome, enzyme, and proteome profiles strongly differ among accessions

Metabolome profiling by 1H-NMR and LC-QTOF-MS allowed the quantification of 18 metabolites from central carbon metabolism and eight secondary metabolites (Supplementary Table S3 at JXB online). In addition, 26 enzyme activities were assessed by robotized assays (Supplementary Table S4). These analyses provided a detailed characterization of sugar, organic acid, and amino acid metabolism pathways, as well as those of glycoalkaloids and phenolic compounds (Fig. 2). Proteins were isolated following 2-DE. A total of 1230 protein spots were detected. A subset of 424 spots whose abundance was significantly different between genotypes or stages were sequenced by LC-MS/MS. A total of 422 spots were identified (Xu ). Supplementary Table S5 lists the mean and standard deviation of every trait for each genotype and stage.

Fig. 2.

Assignment of the metabolites and enzymes studied to pathways. A total of 26 metabolites are indicated in boxes with solid lines, and 26 enzymes are highlighted in boxes with dotted lines. (This figure is available in colour at JXB online.)

The 12 accessions differed for most of the metabolites and phenotypic traits according to the ANOVAs (Table 1). The means of most of the traits (27/29) were significantly different across stages. The content of glucose, fructose, citrate, asparagine, aspartate,and phenylalanine increased from CE to OR, while the other amino acids decreased. The interactions between stage and genotype were significant for 93% of 29 metabolite and phenotypic traits, and 50% of these traits showed a different trend according to the genotype. The data were thus analysed stage by stage (Table 1). A large range of variability was observed among the 12 genotypes at each stage as all the trait means were significantly different, except for that of crypto-chlorogenic acid. The fold change difference between genotypes reached values as high as 5.6 for threonine content at CE or 7.9 for malate at OR.

Table 1.

Analysis of variation for phenotypic and metabolite contents in eight tomato accessions and four F1 hybrids

Traits	Global analysis			CE/OR		CE stage				OR stage
	Fs	Fg	Fgxs	Min	Max	Fg	Min	Max	Max/min	Fg	Min	Max	Max/min
FW	***	***	***	0.10	0.41	***	1.26	30.84	24.41	***	5.33	134.36	25.22
FD	***	***	***	0.46	0.74	***	14.19	42.50	3.00	***	22.73	69.95	3.08
DMC	***	***	NS	1.12	1.42	***	6.16	10.59	1.72	***	4.95	8.81	1.78
Glc	***	***	***	0.42	0.90	***	68 996.30	193 536.82	2.81	***	158574.12	233 208.35	1.47
Suc	***	***	***	0.31	1.08	***	6417.25	11 856.59	1.85	***	7596.82	28 049.81	3.69
Fru	***	***	***	0.38	0.85	***	65 853.14	198 775.33	3.02	***	173280.21	238 532.61	1.38
Ala	***	***	***	0.79	5.51	***	288.06	1356.30	4.71	***	171.01	492.55	2.88
Asn	***	***	***	0.31	0.95	***	752.84	2461.92	3.27	***	1200.65	4593.94	3.83
Asp	***	***	***	0.20	0.41	***	579.58	1278.39	2.21	***	1564.39	3598.28	2.30
Abu	***	***	***	0.80	2.25	***	4099.47	8485.63	2.07	***	2130.82	5854.16	2.75
Gln	***	***	***	0.48	1.45	***	5976.22	26 264.35	4.39	***	6605.55	24 155.80	3.66
Ile	**	***	***	0.57	3.30	***	168.70	865.48	5.13	***	232.59	687.61	2.96
Leu	***	***	***	0.40	1.26	***	306.62	927.64	3.03	***	401.47	1042.48	2.60
Phe	***	***	***	0.35	0.94	***	1209.44	4905.71	4.06	***	1988.74	7429.73	3.74
Tyr	***	***	***	0.50	2.20	***	153.13	760.58	4.97	***	185.86	656.84	3.53
Val	***	***	***	1.15	10.95	***	206.82	1004.13	4.86	***	91.66	396.69	4.33
Thr	***	***	***	0.38	2.12	***	125.21	745.30	5.95	***	221.38	768.59	3.47
Asc	***	***	***	0.48	1.13	***	1304.64	2236.94	1.71	***	1509.86	3021.01	2.00
Cit	***	***	***	0.37	0.57	***	26 295.26	62 463.29	2.38	***	49 525.87	149 942.47	3.03
Mlt	NS	***	***	0.59	3.78	***	13 562.91	25 447.13	1.88	***	3894.13	30 910.20	7.94
Fum	***	***	***	0.88	NA	***	5.54	18.71	3.38	***	0.00	12.75	NA
Tom*	***	***	***	22.16	154.19	***	983 224.09	4 556 008.75	4.63	***	16 992.93	98 910.75	5.82
DHTom*	***	***	***	9.15	63.02	***	196 588.35	2 644 852.27	13.45	***	7239.96	80 840.49	11.17
OH Lyc*	***	***	***	0.19	1.03	***	4384.91	9953.68	2.27	***	5766.39	46972.34	8.15
5CQA*	***	***	***	0.91	1.62	***	51 444.47	111 458.40	2.17	**	43 965.78	68 607.80	1.56
3CQA*	NS	NS	NS	0.45	2.33	NS	4404.24	19 467.52	4.42	NS	5977.49	13 005.07	2.18
Nch*	***	***	***	0.00	0.30	*	408.94	2137.62	5.23	***	2788.25	684 629.20	245.54
Rut*	***	***	***	0.63	4.26	***	63 505.59	243 368.28	3.83	***	28 927.20	265 669.06	9.18
RutP*	***	***	***	0.40	1.09	***	21497.91	59 076.96	2.75	***	25 248.73	79 202.85	3.14

See Supplementary Table S2 at JXB online for metabolite abbreviations.

Fs, significance level of the ANOVA for the stage factor; Fg, significance level of the ANOVA for genotype factor; Fgxs, significance level of the ANOVA for the interaction between genotype and stage; Min, minimum average values for each variable among the 12 genotypes; Max, maximum average values for each variable among the 12 genotypes; CE/OR, ratio of cell expansion and orange-red stages (min and max) for each genotype; NA, not available; NS, non-significant (P > 0.05).

*0.01

*0.01