West Virginia 106 cherry tomato plants were grown in a multi-span Venlo-type greenhouse in 5l pots (potting compost P3 Tref, Tref EGO substrates BV). Water and nutrients were supplied to the plants using a drip irrigation system to maintain 20–30% drainage. Light intensities of 300–700 photosynthetically active radiation were obtained over the culture period. Flowers were mechanically pollinated three times a week and side shoots removed as they appeared. The 35S-RNAi lines used have been previously described [AO, Solyc04g054690 (Garchery et al., 2013), GLD, Solyc10g079470 (Alhagdow et al., 2007), MDHAR (or MD), Solyc09g009390 (Gest et al., 2013a)]. In each case the specific transcript levels and enzyme activity were checked and shown to decrease. The individual lines chosen were the following: GLD5.3, AO15.1, and MD5 with a wild-type (WT) untransformed control. Plants homozygous for the transgene from the T2 generation were used. A total of 15 plants per line were used, 10 for fruit sampling and 5 for destructive phenotypic measurements. All 15 plants were used for non-destructive phenotypic measurements.

Plant height from the pot to the top of the apical meristem was measured weekly for 7 weeks on 10 plants per line. The growth rate was calculated as cm gained per day. The average leaf area of a fully expanded leaf (one leaf per plant, 10 plants) was calculated by scanning the leaf and calculating the area of the image using ImageJ software. Fruit characteristics (number of days between anthesis and orange maturity stage, and fruit weight) were calculated on at least six fruits per plant from 10 plants. Fruit weight, diameter, firmness measurements (compression and maximal force for skin rupture) and color were measured on a minimum of 20 fruits per line immediately after harvesting (ripe fruits). The weight loss of seven fruits per line, related to conservation properties, was measured after 5 days post-harvest at 20°C. Texture measurements were performed using a multi-purpose texturometer (TaPlus, Lloyd). The global firmness of each fruit was measured by a uniaxial compression test with a 50 mm flat disk (test-speed 10 cm.min⁻¹) and expressed as the pressure (kPa) required for a 3% deformation of fruit height. The maximal force (N) corresponding to skin rupture was recorded using a flat cylinder tip (2 mm diameter) for a 7 mm penetration into the equatorial part of the fruit (test-speed 20 mm.min⁻¹). The color was determined using a CR-400 chromameter (Minolta Co. Ltd., Osaka, Japan) and expressed in the CIE 1976 L^*a^*b^* Color Space coordinates (illuminant D65, 0_ view angle, illumination area diameter 8 mm). The chroma a^* axis is a good indicator for tomato ripening, varying from green (−60) to red (+60).

Statistical analysis on phenotypic characteristics was carried out using XLstat (Addinsoft, Paris, France) using a non-parametric comparison of means (Kruskal Wallis) with correction (Dunn). Different letters in the tables indicate significant differences (5% significance level).

Fruit were tagged at anthesis and harvested 20 days later (for ascorbate measurements and proteome analysis only) or when orange. Fruit were harvested during the growing season (from September to November) from trusses 2, 3, and 4 from 10 plants using three fruit/truss/plant to make 3 pools of 30 fruits per line. Fruit pericarp including the epidermis was retained, frozen in liquid nitrogen and stored at −80°C before being ground to a fine powder in liquid nitrogen for the following analyses.

Polar metabolites were extracted from lyophilized powder with a series of hot ethanol/water solutions and quantified by ¹H-NMR as previously described (Mounet et al., 2007). For extract preparation and NMR acquisition parameters, special care was taken to allow absolute quantification of individual metabolites through addition of 5 mM ethylene diamine tetra-acetic acid disodium salt solution to improve the NMR signal resolution and quantification of organic acids, adequate choice of the NMR acquisition parameters, and use through an electronic reference for quantification (Moing et al., 2004). Unknown compounds were quantified in arbitrary units. The lyophilized extracts were titrated with KOD to apparent pH 6 in 400 mM potassium phosphate buffer solution in D₂O and lyophilized again. The ¹H-NMR spectra were recorded on each dried titrated extract solubilized in 0.5 mL D₂O with the sodium salt of (trimethylsilyl)propionic-2,2,3,3-d4 acid (TSP) at a final concentration of 0.01% for chemical shift calibration, at 500.162 MHz on an Avance spectrometer (Bruker BioSpin, Wissembourg, France) using a 5 mm inverse probe. Sixty-four scans of 32K data points were acquired with a spectral width of 6000 Hz, acquisition time of 2.73 s, 90° pulse angle and recycle delays of 25 s. The raw ¹H NMR spectral profiles have been deposited, with associated metadata, into the Metabolomics Repository of Bordeaux MeRy-B (http://services.cbib.u-bordeaux2.fr/MERYB/public/PublicREF.php?REF=T07001).

Ascorbate and dehydroascorbate concentrations in fruit powder were analyzed according to previous protocols based on the reduction of an iron-dipyridyl complex leading to a color change detectable spectrophotometrically with a filter at 525 nm (Stevens et al., 2006) with modifications (Garchery et al., 2013).

Individual polyphenols were analyzed using 200 mg of lyophilised powder from the orange fruit (3 pools), to which of 25 μL of taxifolin (1 mg.mL⁻¹, Extrasynthèse, Genay France) was added as an internal standard. Three successive extractions with 70% ethanol were carried out before evaporation to dryness in a vacuum concentrator. The pellet was resuspended in 300 μL water, 700 μL methanol and filtered with a PTFE syringe filter. Ten micro liters were analyzed by reversed-phase HPLC (Shimadzu Prominence system coupled with a Diode Array Detector (HPLC-DAD), Merck LichroCART® 250-4 column; phase Superspher® 100 RP-18 end capped 4 μm) using a gradient elution with ortho-phosphoric acidified water (pH 2.6) (A) and methanol (B) as solvents. The program was as follows: 0 min 15% B, 10 min 15% B, 20 min 21% B, 32 min 21% B, 42 min 23% B, 52 min 23% B, 107 min 60% B; flow rate 0.5 mL.min⁻¹ from 0 to 44 min, 0.5 mL.min⁻¹ to 0.7 mL.min⁻¹ from 44 to 46 min, 0.7 mL.min⁻¹ from 46 to 60 min. Phenolic compounds were monitored at 280 nm. Chlorogenic acid and rutin were quantified using external calibration curves of standard compounds (Extrasynthèse, Genay, France); other compounds were expressed in arbitrary units.

Protein extraction, two-dimensional gel electrophoresis, and protein identification and classification were carried out as previously described (Xu et al., 2013). Briefly, proteins were extracted using the phenol extraction method developed by Faurobert et al. (2007). Proteins were then separated by two-dimensional gel electrophoresis. After Coomassie colloidal staining, image analysis was performed with Samespot software (version 4.1), and the normalized spot volumes were obtained. Of the 108 spots revealed, 49 were eliminated as not being variable. Protein identification of the 59 variable spots was performed at the proteome platform of Le Moulon (Gif-sur-Yvette) using a nano-LC-MS/MS method following a previously described procedure (Xu et al., 2013). The database search was run against the tomato unigene database version 4 held by the Solanaceae Genome Network (http://solgenomics.net/; Tomato_200607_build_2) with X!Tandem software (http://www.thegpm.org/TANDEM/, version 2010.12.01.1). Sequences were subsequently converted to Solycs (ITAG version 2.4) by Blast analysis. This conversion was confirmed by re-analyzing the mass spectra or by searching for similarity between the protein peptide sequences and the ITAG2.4 version of the genome by TBLASTN.

Microarray experiments were performed on the three biological replicates (3 pools) and technical duplicates with the dyes reversed (dye swap). The TOM2 cDNA microarray [CGE, Boyce Thompson Institute, Cornell University, http://bti.cornell.edu/CGEP/CGEP.html, (Fei et al., 2011)] containing 12,160 70-mer oligonucleotides corresponding to 11,862 unique and randomly selected transcripts of the tomato genome were used. Sample preparation and microarray hybridization and analysis were carried out as described (Mounet et al., 2009). Briefly, one μg of DNA-free total RNA was amplified with the MessageAmpTM aRNA Ambion Kit (Ambion/Applied Biosystems), labeled using the CyScribe Post-Labeling Kit (GE Healthcare) and the cy3- and cy5-labeled cDNAs were hybridized to TOM2 microarrays using an automatic hybridization station HS 4800 (Tecan) as described. Microarray slides were scanned with a Genepix 4000 B fluorescence reader (Axon Instruments) using Genepix 4.0 Pro image acquisition software. The photomultiplier tube voltage was adjusted to 640 V for Cy3 and 700 V for Cy5. Spot flagging was carried out by Genepix (missing spots) and by visual inspection of the images in order to exclude saturated and heterogeneous spots. The raw data were then submitted to LIMMA R Package v2.12.0 (Ritchie et al., 2015) for data visualization, normalization as described (Alhagdow et al., 2007). Background was first subtracted from the raw median intensity values using the normexp function. Two normalization algorithms were successively applied to the raw data (Yang et al., 2002): print-tip loess (default parameters) as the within-array normalization, followed by a quantile normalization (default parameters) as between-array normalization. The normalized data were then analyzed as described below. To convert to current Solyc identifiers, each of the TOM2 probes was blasted against the tomato genome (ITAG2.3, minimum alignment length 40 bp, minimum identity 84%), which had been sequenced in the meantime, and following this 7633 unique Solycs were retained.

Differential expression and transcriptome comparisons were carried out using the R/maanova package [https://www.bioconductor.org/packages/release/bioc/html/maanova.html (R Development-Core-Team, 2015)]. Means and standard errors were calculated for each gene in each line. A simple model taking into account the genotype (line) was used. A matrix of comparisons was used for the pairwise comparison of the transcriptome datasets first of each transgenic line with WT and then with each other. For each comparison, a threshold for the FDR corrected p-value of 0.001 was applied. Following initial tests, we added a sub-hypothesis to the comparison matrix to test whether the AO-GLD transcriptome and the AO-MDHAR transcriptomes were inverted i.e., the wild-type expression was the average of the gene expression in the two transgenic lines being compared. Expression differences were expressed as log2 of the Fold Change between the two lines. The linear model was tested using 1,000 iterations and then an FDR adjustment of the p-values obtained. A threshold for the corrected p-value of 0.001 was set to identify the core gene set for the network visualization.

Weighted gene correlation network analysis (WGCNA) was carried out using the normalized expression data for the 7633 genes with a “trait” file including the 22 metabolites and 40 proteins using the WGCNA package in “R” (Langfelder and Horvath, 2008). No outliers were detected after sample clustering so the whole dataset was used. A soft threshold of 14 and minimum module size of 25 were set. Networks were visualized using Cytoscape version 3.2.0 (Cline et al., 2007) (www.cytoscape.org) using the node and edge output file from the WGCNA package with a threshold set to 0.52 (to fix a number of nodes and edges suitable for visualization). The dataset used contained the 182 differentially expressed genes (differentially expressed in the pairwise comparisons of the three types of transgenic lines), the 22 metabolites and 40 proteins. The number of genes was limited to facilitate network creation and to give weight to the proteins and metabolites. The layout was first designated as degree sorted circle layout and then a degree layout was chosen to aid node identification. Node size is proportional to the degree of connectivity.

Gene ontology (GO) enrichment on the whole transcriptome signatures generated from pairwise comparisons of the lines was carried out using the R package “goseq” (Young et al., 2010) with the corrected p-value false detection rate set to 5%. The package requires the log of fold-change and p-value from the R/maanova analysis. No length correction was required as the data was from microarrays. Tomato gene ontology for the 182 differentially expressed genes was carried out using the online tools available from the Boyce Thompson Institute, Cornell University (http://bioinfo.bti.cornell.edu/cgi-bin/MetGenMAP/home.cgi, Boyle et al., 2004).

The ROSMETER bioinformatic tool allows evaluation of transcriptome signatures linked to reactive oxygen species (ROS) (Rosenwasser et al., 2013). The tool consists of a web-based interface whereby vector-based correlations between stored transcriptome data, chosen for its ROS signature, and the provided data of interest are calculated. The stored data are microarray-based data from wild type and Arabidopsis mutants for the following genes: ascorbate peroxidase, the flu mutant, mitochondrial alternative oxidase, catalase 2 and superoxide dismutase subjected or not to the following stresses: ozone, high light, hydrogen peroxide, drought, rotenone, methyl viologen and 3-aminotriazole as described in the original research paper. The data from the RNAi lines normalized to wild-type (fold change and p-value) were compared with this series of transcriptomes arising from different conditions or mutations in Arabidopsis thaliana. In order to compare the transcriptomes, the tomato Solycs were first converted to their closest Arabidopsis homolog and submitted to the ROSMETER interface. In the heatmaps, red represents a positive correlation and green a negative one.

The three RNAi transgenic lines (AO: ascorbate oxidase; GLD: galactonolactone dehydrogenase and MDHAR: monodehydroascorbate reductase) and wild-type plants grown in the greenhouse were subjected to simple phenotypic analyses to measure plant height, growth rate and final leaf area, the duration of fruit ripening (anthesis to orange fruit) and fruit weight were also measured. The results for the plant phenotypes are shown in Table 1: AO and GLD RNAi lines had a slower growth rate and a smaller final plant height compared to wild-type plants. Fruit from the same two transgenic lines also took longer to reach maturity (Table 2). Fruit from the MDHAR silenced lines were significantly smaller than fruit from the other lines, including wild-type (Tables 2, 3), as per previous results (Truffault et al., 2016). Fruits from the GLD RNAi line were also smaller in diameter (Table 3). Harvested fruits were analyzed for fruit firmness, color, and post-harvest weight loss (Table 3). No significant differences were found between the lines in terms of compression or color but the AO and MDHAR lines required less force to penetrate the skin compared to wild type. MDHAR RNAi fruit also showed greater post-harvest weight-loss compared to the WT, AO, and GLD lines.

Metabolome analysis was carried out on the fruit pericarp samples using ¹H-NMR profiling of polar extracts. Ascorbate (Table 4) and major fruit polyphenols (Table 5) were analyzed by specific methods. Ascorbate was significantly increased in orange fruit pericarp of MDHAR lines as previously observed (Gest et al., 2013a) and dehydroascorbate was significantly different in AO and MDHAR lines, although this difference was the opposite of what would be expected given the activities of the enzymes; GLD and WT had intermediate values. Ascorbate data for fruit 20 days after anthesis is presented in Supplementary Table S1, unlike orange fruit, no significant differences are seen. The polyphenol that varied most in the transgenic lines was chlorogenic acid, which is also the most prevalent polyphenol, and its stereoisomer cis-chlorogenic acid: as for dehydroascorbate, levels were significantly different in AO and MDHAR lines, AO lines containing less chlorogenic acid than wild-type whereas the MDHAR lines contained more (Table 5 and Figure 2). Levels of these polyphenols in the GLD RNAi fruit were not significantly different to wild-type. The metabolites detected and quantified by ¹H-NMR are shown in Figure 2 (and Supplementary Table S2). Metabolites that differentiate the lines included the amino acids alanine and tyrosine which were highest in MDHAR and lowest in AO lines (significant for alanine) and the hexose sugars glucose and fructose which showed a tendency in contrast to the aforementioned amino acids to be highest in AO lines.

The differentially expressed protein spots led to the identification of 34 proteins by mass spectrometry (Table 6). Differences between the lines are presented as a heatmap of the log2 of the ratio with the WT (Figure 3; Table 6, raw data in Supplementary Table S3). The line most affected in terms of its fruit proteome is MDHAR, for which up-regulated proteins included a 60S ribosomal protein, a class I heat-shock protein and a cysteine protease inhibitor. Down-regulated proteins included a heat-shock protein, (hsp70 homolog), lactoylglutathione lyase and an amidase hydantoinase/carbamoylase family protein. The proteome of the GLD line did not show significant differences to the wild-type, despite visually resembling the ascorbate oxidase profile. A number of proteins showed opposite regulation in AO and MDHAR lines, these included chaperone DnaK (MDHAR down-regulated, AO up-regulated), a probable carbonic anhydrase, a protein disulphide isomerase and a metacaspase 7, a protein with cysteine endopeptidase activity and a positive mediator of cell death which is AO down-regulated and MDHAR up-regulated. One protein was up-regulated in the AO line which is an F1F0-ATPase inhibitor protein. The data for fruit 20 days after anthesis is shown in Supplementary Figure S1: only two proteins showed statistically significant differences in spot volumes between wild-type and the MDHAR line.

The normalized transcriptome data for the four lines was subjected to a weighted-gene correlation network analysis to identify clusters of correlated genes and compare them with different phenotypes (metabolite and protein levels). We identified 38 modules of between 26 and 851 genes (Supplementary Figure S2). The number of genes in each module was too small to give significant differences following GO enrichment analysis but revealed co-regulation of metabolites or proteins with genes in modules or groups of modules. For example, the amino acids were positively correlated with the gene modules “dark red,” “blue,” “green-yellow” and “yellow.” The full data set with module membership is found in Supplementary Table S4.

Each of the three transgenic lines was compared with wild-type and then the lines were compared two by two. For the pairwise comparisons of the transgenic lines, the hypothesis that the gene expression in the wild-type was the average of the gene expression in each transgenic line was tested using the maanova model (see Materials and methods). Volcano plots for each pairwise comparison showing the log2 of the fold-change in gene expression between the lines and the associated p-value (−log10 transformed) are shown in Figure 4. The AO-MDHAR and AO-GLD comparisons generated the highest number of differentially expressed genes (comparisons 4; 290 genes and comparison 6; 330 genes). The comparisons WT = (transgenic X + transgenic Y)/2 are comparisons 5, 8, and 9. If the hypothesis is true then the p-values become non-significant and the −log10 of the p-values gives low values which is the case for the comparisons 5 [WT = (AO + MDHAR)/2] and 8 [WT = (AO + GLD)/2] but not the case for the comparison 9 [WT = (GLD + MDHAR)/2]. The whole dataset including means, standard error, p and F-values and the fold-change with p-values for each comparison is available in Supplementary Table S5.

The AO-MDHAR and AO-GLD comparisons generated the highest number of differentially expressed genes and the expression in the WT is always intermediate. Comparison of these two gene sets showed 182 genes common to each comparison, defining a core set of genes for further analysis. Independently of the comparison and the direction of change, differentially expressed genes were part of one “molecular function” category: structural constituent of ribosome (GO:00037) and two “biological process” categories: ribosome biogenesis and translation (GO:00422 and GO:00064). These genes were up-regulated in the AO line and down-regulated in the MDHAR and GLD RNAi lines. The core set of 182 genes were analyzed using the Solanaceae Genome Network GO-enrichment tool and this analysis gave further categories associated with translation, ribosomes, unfolded protein binding, GTPase regulator activity and nucleotide binding. The categories “structural constituent of ribosome,” “small GTPase regulator activity” and “oxidoreductase activity acting on NADH or NADPH” were more than ten-fold enriched in the cluster compared to the genome (Table 7).

The 182 core gene set was used in network analysis with the metabolites and the protein dataset. A topological overlap matrix using the WGCNA script was produced and an edge file for network visualization using Cytoscape. The network visualized with a degree-sorted circle layout, is shown in Figure 5 with the top eleven hub genes. The most highly connected genes included two or three ribosomal proteins, an observation which is compatible with those above, and two heat-shock proteins (including DnaJ) and a heat-shock transcription factor.

Our tomato transcriptome data was compared using ROSMETER, a bioinformatic tool which unites transcriptome data from different stress experiments or mutants in Arabidopsis linked to specific ROS production in specific organelles (see Materials and methods for more information). Its use is based on the hypothesis that the chemical identity and sub-cellular localization of ROS leave a specific imprint on the transcriptomic response (Rosenwasser et al., 2013). We used the RNAi AO and RNAi MDHAR datasets to compare with the Arabidopsis transcriptomes. Our first observation is that independently of the stress and/or organelle, the transcriptome responses are inversed when these two lines are compared (Figure 6). The particular conditions to be highlighted include mitochondrial stress: in this case, the application of rotenone gives a strong negative correlation with the RNAi MDHAR and a positive correlation with RNAi AO 3 h after application. Negative correlations are also found for RNAi MDHAR with Arabidopsis catalase mutants after 3 and 8 h of high light, whereas the correlations are positive for RNAi AO. Positive correlations between RNAi MDHAR exist with 30 min of methyl viologen and generally the knockout ascorbate peroxidase mutant subjected to high light. Again for RNAi AO, in contrast to RNAi MDHAR, a negative correlation is found for 30 min of methyl viologen treatment.

The transcriptomes revealed a switch linked to the activity of ascorbate oxidase on one side of the switch and the enzymes MDHAR/GLD acting in the other direction, always with wild-type gene expression in the middle. There is a tendency for genes to be up-regulated when ascorbate oxidase is down-regulated. The transcriptome inversion is also seen in datasets from Arabidopsis provided by the ROSMETER tool. The switch seen between AO and MDHAR RNAi lines is capable of distinguishing the conditions found in some of the ROSMETER transcriptome data, although there are limits to this analysis because correlations can be based on a small number of genes (a minimum of 45 according to the original paper). The ROSMETER analysis does not reveal links with hydrogen peroxide metabolism because correlations between our lines and the public data are inversed when comparing ascorbate peroxidase and catalase mutants, both of which in theory should contain more hydrogen peroxide as a result of decreased detoxification. The strong response to rotenone is interesting because rotenone is an inhibitor of mitochondrial complex I, which is also where the galactonolactone dehydrogenase protein is located (Schimmeyer et al., 2016). Also in the paper by Rosenwasser et al. (2013), rotenone-induced stress is also positively correlated with developmental senescence, a process which has parallels with fruit ripening.

The data seem to indicate the existence of a signal affecting transcription closely linked to GLD activity and the redox state of the ascorbate pool. The nature of this signal is unknown but a candidate is monodehydroascorbate, as it is the molecule that links ascorbate oxidase and monodehydroascorbate reductase activity, which could act as a primary signal given (i) the short half-life of the molecule, a prerequisite for a role in signaling and (ii) the proton transfer that occurs when the radical is reduced back to ascorbate (Bielski et al., 1981) which gives potential for protein modification. Most oxidation reactions with ascorbate are indeed single electron reactions (Heber et al., 1996) generating the radical monodehydroascorbate, this radical is considered to be common and a marker for stress. Another possibility is that redox changes in mitochondria act as a retrograde signal (Rhoads and Subbaiah, 2007); this could be a possibility given the involvement of GLD. Interestingly a dehydroascorbate reductase (DHAR) gene (Solyc05g054760) is one of the 182 differentially regulated genes and is down-regulated in GLD and MDHAR lines and up-regulated in the AO line. Compensatory effects, either by enzymes such as DHAR, or other isoforms of the enzymes under study, cannot be ruled out.

The transcriptional signature and network analysis suggest a role for ascorbate in control of translation and protein folding, at least in the non-photosynthetic tissue studied. Some evidence from the literature has established links between ascorbate and these cellular mechanisms: a good example is the posttranscriptional negative feedback control that ascorbate exerts over its synthesis via the GDP-l-galactose phosphorylase enzyme (Laing et al., 2015) which appears to affect ribosome density on the 5′UTR of the mRNA. Ribosomes and translation are also regulated by stress and redox processes (Halbeisen and Gerber, 2009; Khandal et al., 2009; Mustroph et al., 2009; Gerashchenko et al., 2012; Gismondi et al., 2014; Moore et al., 2016). Protein folding is also affected by environmental challenges: the unfolded protein response results from a period of stress and leads to either up-regulation of proteins required for folding or degradation or in extreme cases cell death if prolonged (Ruberti et al., 2015; Wan and Jiang, 2016). Functional hypotheses such as these would need to be tested. For example, purification of polysomes from the available transgenic lines could lead to information on the RNAs undergoing translation and/or associated proteins and ribosome subunits. Considering the fruit context, this experiment would be especially informative in ripening fruits subjected to high temperature stress or in ripe fruits subjected to stressful postharvest storage conditions that may result in protein mis-folding such as high CO₂ (Rothan et al., 1997) or chilling temperatures (Stevens et al., 2008).

Heat-shock proteins also figure amongst the hub genes identified: three heat-shock proteins are identified: SolycHsfB2a (Solyc03g026020), a Class II small heat-shock protein Le-HSP17.6 (Solyc08g062340) and a chaperone protein DnaJ (Solyc03g123560). Heat-shock proteins are known to be involved in numerous biological phenomena including stress responses, circadian clock regulation, flowering time, and hypocotyl elongation (Charng et al., 2007; Mittler et al., 2012; Kolmos et al., 2014). The hub gene protein DnaJ (hsp40) stimulates the ATPase activity of DnaK (an hsp70 homolog which is interestingly differentially accumulated in the AO and MDHAR proteomes). These proteins bind to unfolded polypeptides to prevent their aggregation, for example under conditions of heat-shock. The DnaK protein is of interest because it has been detected in the proteome of the Arabidopsis mutant for GDP-L-galactose phosphorylase (vtc2), a major control point for ascorbate synthesis, as being up-regulated after 3 and 5 days of high light compared to the wild-type (Giacomelli et al., 2006). This parallels the situation in the AO line (up-regulation) and is opposite to the MDHAR line phenotype (down-regulation) meaning that ascorbate deficiency mirrors the RNAi AO phenotype. Another hub protein, Arabinogalactan, is involved in many biological processes including germination, cell extension, fertilization, binding, and release of calcium and is linked to the auxin response (Ellis et al., 2010). Only one of the hub genes has already been characterized from tomato (Table 7) which is the tft8 protein, a member of the 14-3-3 protein family involved in defense responses (Roberts and Bowles, 1999).

The major phenotypic differences seen in the transgenic lines are related to growth: in previous studies RNAi AO lines tend to show improved yield stability (Garchery et al., 2013) whereas the opposite is true for GLD and MDHAR RNAi lines (Alhagdow et al., 2007; Truffault et al., 2016), a phenomenon which is correlated to the changes in the fruit transcriptomes and could be related to fruit growth. Explanations for the phenotypes might be found in the identities of the hub genes which are broadly involved in signaling translation and protein (re)folding. We also noted the presence of many ATP-dependent proteins. The yield phenotypes could therefore reflect differences in cellular energy levels.

As many studies have observed, RNA and protein abundance are often weakly correlated (Halbeisen and Gerber, 2009; Petricka et al., 2012) and our data corroborate this: the quantity of the differentially accumulated proteins is generally not linked to their transcriptional levels (Table 6). For the observable proteins, the GLD proteome is not significantly affected: the GLD proteome actually is remarkably visually similar to the AO proteome which is surprising as the transcriptome data showed contrasted differential expression of genes for these two lines. The metabolic profile is similar to the proteome profile in that several metabolites show opposite changes compared to wild-type between AO and MDHAR RNAi lines, this includes dehydroascorbate, chlorogenic acid, and alanine. The changes seen in dehydroascorbate are counter-intuitive when compared to what is expected based on enzyme activity. We, with other independent groups, have previously reported increases in ascorbate levels on reduction of MDHAR activity (a peroxisomal and cytosolic isoform; Haroldsen et al., 2011; Gest et al., 2013a) which implies feedback regulation. Links have been shown between polyphenols such as chlorogenic acid or anthocyanins and ascorbate (Page et al., 2012); also oxidation of chlorogenic acid is correlated with oxidation of ascorbate to dehydroascorbate (Takahama, 1998) and ascorbate can regenerate chlorogenic acid from its corresponding radicals (Yamasaki and Grace, 1998). Our results show less chlorogenic acid in lines with lower ascorbate oxidase. We should point out that the genotype used in this study (West Virginia 106) is particularly rich in chlorogenic acid, the wild-type contains 3.25 mg/g dwt of chlorogenic acid whereas reported concentrations range from 0.24 to 3.69 mg/100 g fresh weight (between ~0.024 and 0.37 mg/g dwt; Slimestad and Verheul, 2005; Peng et al., 2008). The differences in alanine are also of interest as this metabolite is produced directly from pyruvate and therefore has close links with glycolysis, the oxidative pentose phosphate pathway and the TCA cycle and cellular energy production. Alanine is a product of anaerobic metabolism resulting from high rates of amino acid interconversion (Rocha et al., 2010b). Unlike lactate and ethanol, alanine accumulation does not have detrimental side effects but maintains the glycolytic flux (Rocha et al., 2010a).

In conclusion, this study in a non-photosynthetic tissue reveals that the activity of enzymes controlling the oxidation state of the ascorbate pool are correlated to a few metabolic changes and some changes in protein quantity, particularly those involved in protein folding. The relatively low number of metabolic changes compared to gene expression probably reflects buffering of the system. In addition to this post-transcriptional regulation, a transcriptional response is seen which is also linked to the activity of the final biosynthesis enzyme, galactonolactone dehydrogenase. Plants silenced for this enzyme show a transcriptional response similar to that of MDHAR silenced plants and opposite to ascorbate oxidase silenced plants implying signaling from the ascorbate pool. Network analysis reveals hub genes involved in fundamental cellular mechanisms including cellular signaling, translation and protein synthesis as well as the heat-shock response.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.