Transcriptomic changes in Cucurbita pepo fruit after cold storage: differential response between two cultivars contrasting in chilling sensitivity

Background Zucchini fruit is susceptible to chilling injury (CI), but the response to low storage temperature is cultivar dependent. Previous reports about the response of zucchini fruit to chilling storage have been focused on the physiology and biochemistry of this process, with little information about the molecular mechanisms underlying it. In this work, we present a comprehensive analysis of transcriptomic changes that take place after cold storage in zucchini fruit of two commercial cultivars with contrasting response to chilling stress. Results RNA-Seq analysis was conducted in exocarp of fruit at harvest and after 14 days of storage at 4 and 20 °C. Differential expressed genes (DEGs) were obtained comparing fruit stored at 4 °C with their control at 20 °C, and then specific and common up and down-regulated DEGs of each cultivar were identified. Functional analysis of these DEGs identified similarities between the response of zucchini fruit to low temperature and other stresses, with an important number of GO terms related to biotic and abiotic stresses overrepresented in both cultivars. This study also revealed several molecular mechanisms that could be related to chilling tolerance, since they were up-regulated in cv. Natura (CI tolerant) or down-regulated in cv. Sinatra (CI sensitive). These mechanisms were mainly those related to carbohydrate and energy metabolism, transcription, signal transduction, and protein transport and degradation. Among DEGs belonging to these pathways, we selected candidate genes that could regulate or promote chilling tolerance in zucchini fruit including the transcription factors MYB76-like, ZAT10-like, DELLA protein GAIP, and AP2/ERF domain-containing protein. Conclusions This study provides a broader understanding of the important mechanisms and processes related to coping with low temperature stress in zucchini fruit and allowed the identification of some candidate genes that may be involved in the acquisition of chilling tolerance in this crop. These genes will be the basis of future studies aimed to identify markers involved in cold tolerance and aid in zucchini breeding programs. Electronic supplementary material The online version of this article (10.1186/s12864-018-4500-9) contains supplementary material, which is available to authorized users.

modifications common to other stresses such as an accumulation of hydrogen peroxide (H 2 O 2 ) and malondialdehyde (MDA), as well as changes in the levels of endogenous abscisic acid [3,4], ethylene [5,6], polyamines [7] and in soluble sugars [8]. Moreover, previous works have demonstrated the existence of genetic variability for CI tolerance among commercial and local cultivars of zucchini squash [6,9]. The fruit of the most CI-tolerant cultivars produces less chilling-induced ethylene and accumulates lower content of H 2 O 2 and MDA, two metabolites associated with oxidative stress. Among commercial hybrids, fruit from Natura showed very low CI after 14 days of cold exposure, while fruit from cv. Sinatra presented high CI index and an important loss of fruit quality [9]. Subsequent physiological research revealed many differences between these two cultivars in some metabolic pathways involved in chilling stress and tolerance, confirming that Natura and Sinatra should be considered tolerant and sensitive cultivars to cold stress, respectively. In this sense, Natura fruit has higher soluble sugar content, higher levels of proline, lower content in chilling stress metabolites such as H 2 O 2 or MDA, and higher gene expression and activity of antioxidant defense enzymes, than Sinatra fruit [7,10,11]. In relation to the hormones involved in cold tolerance, it has been reported that the fruit of the most CI-tolerant cultivars shows a reduced induction of ethylene biosynthesis and signaling pathways under cold storage, and that the treatment with 1-MCP prevents chilling damage [6,11]. Recently we have also detected an increase in the synthesis of abscisic acid (ABA) during the first days of cold storage in the more cold-tolerant cultivar [3]. Hence, there is an increasing amount of information concerning the physiology and biochemistry of chilling in zucchini fruit, however, the molecular mechanisms underlying the response of zucchini fruit to chilling storage is limited to specific molecular pathways, including those of ethylene biosynthesis and signaling [5,6], enzymatic antioxidant system [12,13], abscisic acid synthesis and signaling, and polyamine metabolism [3,7].
In the present work, fruit of the CI-tolerant and sensitive cultivars Natura and Sinatra were stored at chilling (4°C) and non-chilling (20°C) temperature for 14 days, and their transcriptomic profiles examined by RNA-Seq. The analysis of the transcriptomic changes in response to cold stress between these cultivars contrasting in their sensitivity to chilling will provide new insight into important mechanisms and processes related to resistance against low temperature stress. Moreover, the outcomes of this study will be the basis for future studies aimed to identify markers involved in cold tolerance, which will surely improve the breeding programs of this crop.

Plant material and postharvest treatments
The commercial zucchini hybrids Natura (Enza Zaden) and Sinatra (Clause-Tezier) were grown under the same greenhouse condition in Almeria, Spain (FEMAGO S.L.). After harvest, fruits of each cultivar were stored in chambers at 4°C and 20°C during 14 days. Fruits were divided into three replicates per cultivar and storage period (0 and 14 days), each consisting in 6 fruits of similar size. After storage, weight loss, electrolyte leakage, and chilling injury-index were measured, and the exocarp tissue of each replicate was frozen in liquid nitrogen and stored at − 80°C.

Weight loss and chilling-injury index
The percentage of weight loss of each fruit was calculated as: % weight loss = (W i − W f )/W i × 100, being W i the initial fruit weight and W f the final fruit weight. Chilling injury index of the fruit surface was evaluated in fruit stored at 4°C using a subjective scale of visual symptoms previously described [1]: 0 = no pitting, 1 = slight (10% or less), 2 = medium (10-20%), and 3 = severe pitting (> 20%). CI index was determined using the following formula: Ʃ (pitting scale (0-3) × number of corresponding fruit within each class)/total number of fruit estimated.

Electrolyte leakage
Electrolyte leakage was measured as described [14]. Briefly, exocarp of zucchini fruit was separated with a vegetable peeler and 10 discs were taken from each replicate with an 11 mm diameter stainless steel cork borer. Each replicate was rinsed with 50 mL of deionized water three times for 3 min. After being incubated for 30 min and shaken at 100 rpm in 50 mL of deionized water, this solution was measured for conductivity at room temperature using a conductimeter (Consort C860 provided with a conductivity electrode Consort SK10T, Consort nv, Belgium). Total conductivity was determined after boiling the flasks for 10 min and cooling at room temperature. The electrolyte leakage was expressed as percentage of total conductivity.

Measurement of lipid peroxidation
Lipid peroxidation was determined as malondialdehyde (MDA) content using the procedure previously described [15], with some modifications. Exocarp ground in liquid nitrogen was homogenized (1:4, w/v) in 20% (w/v) trichloroacetic acid (TCA) and butylated hydroxytoluene was added to a final concentration of 0.67%. The homogenate was centrifuged at 4°C and 10,000×g for 15 min. The supernatant was mixed with 0.5% (w/v) thiobarbituric acid (TBA) in 20% TCA in proportion 1:4 (v/v). The mixture was heated at 95°C in a water bath for 30 min, cooled immediately in ice to stop the reaction, and centrifuged at 4°C and 4000×g for 10 min. Absorbance of supernatant was measured at 532 and 600 nm. MDA content was calculated by subtracting the non-specific absorption at 600 nm from the absorption at 532 nm and using a standard curve. Results were expressed as nmol MDA g − 1 of fresh weight.
Determination of H 2 O 2 content H 2 O 2 content was assayed as described [16]. Zucchini exocarp was ground in liquid nitrogen and homogenized with 0.1% (w/v) TCA (1:4, w/v). After centrifugation at 4°C and 12,000×g for 15 min, the supernatant was collected. The reaction mixture consisted of 0.25 mL supernatant, 0.25 mL 100 mM potassium phosphate buffer (pH 7) and 1 mL 1 M KI. The reaction was developed for 1 h in darkness and the absorbance measured at 390 nm. The amount of H 2 O 2 was calculated using a standard curve and expressed as μmol H 2 O 2 g − 1 of fresh weight.

RNA extraction
Total RNA was extracted from the exocarp of 6 fruits for each replicate as reported [17], treated with RNAse-Free DNAse and purified using RNasy ® MiniElute™ Cleanup columns (Qiagen, Hilden, Germany). The quality and quantity of RNA was determined by agarose gel electrophoresis and NanoDrop Lite spectrophotometer (Thermo Fisher Scientific, MA, USA).
Sequencing data processing and gene expression analysis RNA samples were sequenced using Illumina Hiseq2000 at Boyce Thompson Institute (Ithaca, NY, USA). The quality of the single reads generated by Illumina was checked using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). In order to obtain highquality data, the raw reads were pre-processed and trimmed using the software NGS_CRUMBS (https:// bioinf.comav.upv.es/ngs_crumbs/). Through the different utilities the adapters used during the sequencing process were removed, as well as low quality sequences with a Phred quality score Q < 20 and ambiguous sequences with N. Using bowtie2 [18], the high quality reads were mapped against Cucurbita pepo transcriptome v3.0 which is included in the genome version 4.1 but was not available when this study was performed [19]. The expression levels were calculated and normalized by the FPKM method with RSEM [20]. Differential expression transcripts were identified using DESeq2 package [21] of the bioconductor package [22,23]. Transcripts with an adjusted padj (p-value adjusted for multiple comparisons using Bejamini-Hochberg method) < 0.05 and a log2 fold change (FC) ± 1.5 based in three biological replicates were considered as DEGs. Principal component and clustering analysis were performed with Mev software [24].

Gene ontology (GO) terms enrichment analysis
BlastoGO software (v2.8.0) [25] was used for GO term differential analysis using the Cucurbita pepo transcriptome v3.0 GO annotation which contains 24,402 annotated unigenes. GO terms enrichment for each data set was calculated by a binomial test model with FDR cut off of 0.05.

Gene expression analysis by qRT-PCR for RNA-Seq validation
The expression patterns of 10 random DEGs identified by RNA-Seq in this study were validated by quantitative RT-PCR. Primers pairs for each gene (Additional file 1: Table S1) were designed using Primer3 web tool (http:// bioinfo.ut.ee/primer3-0.4.0/primer3/). Total RNA was extracted as above. First-strand cDNA was synthesized from 1 μg total RNA using Maxima Reverse Transcriptase (Thermo Fisher Scientific, Rockford, IL, USA). For qRT-PCR, amplifications were run in a 96-well-plates iCycler iQ thermal cycler (Bio-Rad) using iQ SyBr Green Supermix (BioRad). Quantification was performed with the iCycler iQTM associated software (Real Time Detection System Software, version 2.0). The relative gene expression was calculated using non-stored fruit as the calibration sample. EF-1α was used as the internal reference gene for normalizing the transcript profiles following the 2 -ΔΔCt method [26].

Statistical analysis
The experimental design was completely randomized. Data were subjected to analysis of variance (ANOVA) or unpaired t-test using Statgraphics Centurion XVI (Statpoint Technologies, Inc., Warrenton, VA, USA). When appropriate, means were separated by Tukey's HSD test and differences at p < 0.05 were considered significant.

Fruit physiological parameters and incidence of chilling injury
Natura and Sinatra fruit were stored at 4 and 20°C during 14 days and postharvest quality parameters including percentage of weight loss, CI index, electrolyte leakage, lipid peroxidation (as MDA content), and H 2 O 2 were recorded (Table 1). Storage at 4°C was effective in reducing the fruit weight loss observed in zucchini fruit stored at 20°C. On the other hand, chilled fruit from Natura and Sinatra showed greater membrane permeability, lipid peroxidation, and H 2 O 2 content than non-chilled fruit, which presented values similar to fresh harvested fruit in both cultivars. In spite of the clear effect of chilling, the extent of the changes on quality parameters was more pronounced in Sinatra than in Natura. Our data confirmed that cold-stored Sinatra fruit had a greater loss of quality, showing higher weight loss, CI index, electrolyte leakage, lipid peroxidation, and H 2 O 2 content than Natura fruit.

Differential gene expression in the cold stored fruit of the two cultivars
The molecular network regulating zucchini fruit response to chilling was studied by performing a RNA-Seq analysis from exocarp of Natura and Sinatra fruit before and after 14 days of storage at 4°C and 20°C. A total of 146 million single reads with an average of 8.1 million reads per sample were generated by Illumina Hiseq2000. After pre-processing and trimming, 133.4 million high quality reads were obtained (an average of 7,412,696 per sample) (Additional file 2: Table S2) and mapped against Cucurbita pepo transcriptome v3.0 [19]. To analyze the complexity of the transcriptomic data and to cluster samples according to their gene expression profile, we first performed a principal component analysis (PCA) over the expression data of the 18 biological samples (Fig. 1). The analysis showed that in all conditions the gene expression profile of the three independent biological replicates clustered together; thus the experiment was considered reliable for further analysis. Furthermore, the PCA revealed that at harvest both cultivars presented a different gene expression pattern. Within each cultivar, fresh-harvested fruit clustered away from stored fruit in PC1, which explains 42.9% of the variation. On the other hand, the differences between 4°C and 20°C stored fruit in Sinatra are explained by both, PC1 and PC2 (38% of the variation). Interestingly, the changes in the gene expression pattern between Natura chilled and non-chilled fruit were smaller, clustering close together with little separation in either axes.
Gene expression was compared using a pairwise analysis. The comparison of the transcriptome of Natura and Sinatra fruit at harvest (Natura FH/Sinatra FH) revealed 503 DEGs, 262 of which were up-regulated and 241 down-regulated in Natura versus Sinatra ( Table 2). The changes in the transcriptomic profiles of the two cultivars were then identified by comparing gene expression values in 4°C and 20°C stored fruit (cold-stored/ control) in each cultivar independently ( Table 2). In Natura, there were 5636 DEGs in cold-stored fruit, from which 2522 were up-regulated (45%) and 3114 downregulated (55%). Similarly, in cold-stored fruit of Sinatra there were 6623 DEGs, of which 2845 genes were upregulated (43%) and 3778 genes were down-regulated (57%). To validate the results of transcriptomic profiling, the expression of 10 DEGs was analyzed by quantitative RT-PCR. Linear regression analysis was conducted using Log2 fold change obtained by each approach (Additional file 3: Figure S1). Results showed a significant positive correlation between both methods (R 2 = 0. 956), thus validating the transcriptomic results obtained by RNA-Seq. DEGs of each cultivar were compared using a Venn diagram to identify specific genes related to cold tolerance or sensitivity (Fig. 2). Results showed that 2682 DEGs were specific to Natura and 3669 were specific to Sinatra. Regarding up-regulated DEGs, 1221 were specific of Natura and 1568 were exclusively identified in Sinatra (45% and 43% respectively). Down-regulated genes differentially regulated in Natura were 1461; while in the case of Sinatra, 2101 were specific, representing 54% and 57% of the DEGs in each cultivar, respectively. The percentage of specific DEGs induced and repressed in response to chilling was similar in both cultivars. On the other hand, Natura and Sinatra shared 2954 DEGs, from which 1264 were induced (43%) and 1640 were down-regulated (55%) in both cultivars. Moreover, there was a group of common DEGs showing opposite cold-specific regulation, i.e. 37 DEGs were up-regulated in Natura and down-regulated in Sinatra, while 13 DEGs were repressed in Natura and induced in Sinatra. The complete list of specific and common cold-regulated genes as well as those found differentially expressed between cultivars at harvest is available in Additional file 4: Table S3.
Functional analysis of the differential gene expression in fruit of the two cultivars To explore the biological functions of the cold-induced and repressed DEGs of each cultivar, a gene ontology (GO) enrichment analysis was conducted. Overrepresentation of GO terms was evaluated to correlate different biological processes (BPs), molecular functions (MFs), and cellular components (CCs) with cultivar-dependent chilling response in the exocarp of zucchini fruit.
Functional analysis of the differential gene expression in fruit of the two cultivars before storage Regarding the comparison between Natura and Sinatra fruit at harvest, 5 GO categories (4 BPs and 1 CC) were up-regulated and 3 (2 BPs and 1 MF) were downregulated in Natura compared to Sinatra. The BPs represented by the largest number of DEGs up-regulated in Natura fruit compared to Sinatra fruit before storage were 'response to oxidative stress' , 'response to cold' and 'fatty acid biosynthetic process' , while the three GO categories down-regulated in these fruit were the BPs 'Glycoxylate metabolic process' and 'glycoxylate cycle' and the MF 'Isocitrate lyase activity' (Fig. 3, Additional file 5: Table S4).
Functional analysis of the differential gene expression in cold-stored fruit of the two cultivars Figure 4 shows the overrepresented GO terms in DEGs that were specific for the cold tolerant cv. Natura or for the cold-sensitive cv. Sinatra, and in DEGs that were common to both cultivars. In Natura-specific upregulated DEGs several GO terms were found that belong to abiotic-and biotic-stress responses (16% of up-regulated DEGs) such as response to 'metal ion' , 'salt stress' , 'temperature stimulus' , and 'defense response to bacterium' (Fig. 4a). Interestingly, three BPs upregulated in Natura, i.e. 'response to misfolded protein' , 'proteasome core complex assembly' , and 'proteasomal ubiquitin-dependent protein catabolic processes' , were related to proteolysis (Additional file 6: Table S5). Other BPs overrepresented in Natura were 'glycolysis' , 'toxin catabolic process' , and 'hydrogen peroxide catabolic process'. The functional analysis of up-regulated DEGs in Natura also revealed that the largest MF categories in this cultivar were 'threonine-type endopeptidase activity' , 'transferase activity' , and 'protein domain specific  Fig. 2 Venn diagram of the differentially expressed genes in fruit from Natura and Sinatra exposed to cold storage (4°C vs 20°C). Bold-italic numbers depict the DEGs specifically up-regulated due to cold storage in each cultivar. Italic numbers depict the DEGs specifically down-regulated due to cold storage in each cultivar binding' (Additional file 7: Figure S2A), while 'cytosol' , 'chloroplast stroma' , and chloroplast envelope' were the CCs represented with the largest number of DEGs (Additional file 8: Figure S3A). On the other hand, only two enriched GO terms resulted from the functional analysis of down-regulated DEGs specific of Natura. These GO terms were the BP 'meristem maintenance' and the CC 'cell periphery' (Fig. 4a and Additional file 8: Figure S3). The former is integrated by DEGs that control the expression of the meristem genes and regulates the response to hormones such as auxins, while the latter includes DEGs whose function are expressed in plasma membrane or cell wall. The functional analysis of cold response in Sinatra fruit revealed that there were no overrepresented GO terms in upregulated genes, while there were many terms overrepresented in cold-repressed DEGs (Additional file 9: Table S6). Regarding BPs, 21% of down-regulated DEGs in Sinatra were represented in GO terms related to exposure to biotic and abiotic stress conditions such as 'defense response to other organism' , 'response to salt stress' , 'response to cadmium ion' , 'response to cold' , 'response to abscisic acid stimulus' , 'response to wounding' , and 'response to water deprivation' (Fig. 4b). The MFs more overrepresented in down-regulated DEGs were related to ion binding (Additional file 7: Figure S2B) and the largest CCs were the 'plasma membrane' , 'plasmodesma' , and 'chloroplast thylakoid membrane' (Additional file 8: Figure S3B).
Up-regulated transcripts common to both cultivars were only enriched in three BPs; 'cellular response to nutrient levels' , 'heat acclimation' , and 'negative regulation of endopeptidase activity' (Fig. 4c). The most overrepresented MFs were related to oxidation-reduction reactions, 'monooxygenase activity' , and 'oxidoreductase activity' (Additional file 7: Figure S2C). On the other hand, the functional analysis of cold-repressed DEGs common to both cultivars (Fig. 4c, Additional file 10: Table S7) revealed that most of the overrepresented BPs were related to response to several stimuli ('response to hormone stimulus' , 'response to other organism' , and 'response to osmotic stress'), cell wall and cuticle biosynthesis ('pectin metabolic process' , 'cuticle development' , and 'wax biosynthetic process'), and photosynthesis and exposition to light ('photosynthesis, and light harvesting' , 'non-photochemical quenching' and 'anthocyanin accumulation in tissues in response to UV light'). The largest MF in down-regulated DEGs from both cultivars was 'transferase activity' (Additional file 7: Figure S2C) and the CCs were 'plastid envelope' , 'apoplast' , and 'plant-type cell wall' (Additional file 8: Figure S3C).

Transcription factors
Around 289 and 314 DE transcription factor (TF) genes were identified in Natura and Sinatra fruit, respectively, when 4°C and 20°C stored fruit of each cultivar was compared (Additional file 4: Table S3). From these, 47 and 60 were specifically up and down-regulated in Natura respectively; most of them belonged to the zinc finger, NAC, and MYB families. With respect to the cold-sensitive Sinatra, the number of specific downregulated TFs was much larger (96) than the upregulated ones (30) and they were part of the zinc finger, WRKY, and AP2/EREBP (including ERFs and CBF/ DREBs) families. From these DE TF genes, those belonging to overrepresented BPs identified in our functional analysis are listed in Table 3.

Putative candidate genes for cold tolerance in zucchini fruit
In order to identify genes more likely to be related to acquisition of chilling tolerance in zucchini, a careful analysis of all DEGs in our dataset was carried out. The selected genes were those showing a differential expression pattern between Natura and Sinatra and were differentially expressed in cold-stored fruit compared to fruit stored at 20°C and/or compared to freshly-harvested fruit. The selected candidate genes were grouped according to the following metabolic pathways: carbohydrate and energy metabolism, lipid metabolism, peptide transport, transcription, and signal transduction (Table 4). Their expression patterns and their possible role in cold-tolerance in zucchini fruit will be discussed below.

Discussion
Zucchini fruit is known to be sensitive to cold storage; although the degree of chilling susceptibility is very dependent on the cultivar [6,9]. Present fruit quality parameters in Natura and Sinatra fruit after cold storage confirmed previous reports indicating that chilling affects the fruit of the two cultivars, although Natura fruit were more tolerant to chilling and withstood cold better than Sinatra fruit. Furthermore, the chilling tolerance induced by diverse postharvest treatments, including preconditioning at moderate temperature before cold storage [12] or individual shrink wrapping [11], correlates with a better fruit antioxidant status. However, a holistic overview of the molecular bases of the zucchini response to chilling is still lacking.
Transcriptional bases for the differential response of fruit of the two zucchini cultivars to cold storage The PCA analysis showed a differential transcriptomic response between cultivars; Sinatra fruit stored at 20°C were clustered away from fruit stored at 4°C , whereas chilled and non-chilled fruit from Natura were grouped together. Likewise, the number of DEGs in cold-stored fruit compared to fruit stored at 20°C was higher in cv. Sinatra than in cv. Natura. This result suggests larger transcriptomic changes in the cold-sensitive Sinatra after cold storage. On the other hand, 2954 DEGs were shared among Sinatra and Natura, reflecting that some of the chilling responses to cold storage are common   in fruit of both cultivars. In this sense, the ratio of up-and down-regulated genes was also similar in both cultivars, meaning that the number of repressed genes was slightly higher than the number of induced genes. Transcriptomic studies in other species suggest that this ratio of differential gene expression in response to cold may not be related to acclimation but rather to a species dependent response. Thus, in several species such as Arabidopsis thaliana or Camellia sinensis, it has been reported that in response to cold the number of up-  regulated genes was larger [27,28], whereas in other species such as Populus simonii or Solanum lycopersicum this number was either lower or the same [29,30].

Molecular mechanisms related to chilling tolerance in zucchini fruit
The functional analysis of DEGs in zucchini fruit before storage revealed that the response to abiotic stresses ('response to cold' , 'response to oxidative stress' and 'regulation of the circadian clock by temperature') was up-regulated in Natura fruit when compared to Sinatra. Interestingly, the most overrepresented BPs in cold-induced genes from Natura fruit after cold storage were also those related to abiotic ('response to metal ion' , 'response to salt stresses, and 'response to temperature stimulus') and biotic ('defense response to bacterium') responses to stress conditions. On the other hand, in the sensitive cultivar Sinatra, BPs related to biotic ('defense response to other organism') and abiotic ('response to salt stress' , 'response to cadmium ion' , 'response to cold' , 'response to abscisic acid stimulus' , 'response to wounding' , and 'response to water deprivation') stresses were also overrepresented but the genes were down-regulated. Our data suggest that this differential regulation may be related to the degree of sensitivity of zucchini fruit cultivars to low temperature and may be fundamental in preventing the negative effect of cold storage in the more chilling-tolerant cultivar. Furthermore, our results support previous reports showing a cross-talk between different stresses, as is the case of cassava apical shoots [31], asparagus bean seedlings [32], and table grapes [33] exposed to chilling. Other BPs that were highly enriched from the analysis of specifically up-regulated DEGs in Natura were those involved in the maintenance of energy and redox status, including 'glycolysis' , 'fatty acid β-oxidation' and 'pentosephosphate shunt'. Among these DEGs were pyruvate kinase (CUUC9946, CUUC104231), phosphoglycerate kinase (CUUC91403), glyceraldehyde-3-phosphate dehydrogenase (CUUC117427), ribulose-phosphate 3-epime rase (CUUC104592), phosphoribulokinase (CUUC 62669), citrate synthase (CUUC110803), malate dehydrogenase (CUUC107944), and fructose-1,6-bisphosphatase (CUUC111898, CUUC110408). Enhanced expression of these genes may result in activation of carbohydrate catabolism and therefore in a rise of the energy supply for Natura fruit, which agrees with previous reports indicating a higher pool of ATP and high energy status on fruit of this cultivar [7]. Similar findings were revealed by a proteomic analysis in leaves of a cold-tolerant maize genotype [34]. Likewise, Cai and coworkers [35] showed that cold-stored grapes treated with salicylic acid had a better postharvest quality and presented a higher accumulation of an important number of proteins belonging to carbohydrate and energy metabolism pathways. Natura may also increase the ATP production by enhancing respiration, since a delta subunit of two mitochondrial ATP synthase genes (CUUC101699 and CUCC101701) was specifically upregulated in Natura chilled fruit.
The BPs toxin and hydrogen peroxide catabolic processes were also overrepresented in the chilling-tolerant cultivar. These results suggest that concomitant to an increase in the defense response to stress, the trigger of some mechanism of detoxification may allow Natura fruit to avoid or reduce the oxidative stress associated with cold storage. Several DEGs in this group encode detoxification enzymes such as catalase (CUUC61495), peroxidase 2 (CUUC106379), or NADPH-dependent thioredoxin reductase 3-like (CUUC65543), all of them related to ROS scavenging. This differential expression in Natura would explain the lower H 2 O 2 levels measured in this cultivar; i.e. cold-stored Natura fruit maintained H 2 O 2 levels similar to fresh or 20°C stored fruit, but much lower than those observed in the cold-stored fruit from the sensitive cultivar Sinatra. This is also in accordance with previous reports indicating that the chilling tolerance induced by different postharvest treatments correlated with a high catalase activity and gene expression in zucchini fruit [12,36].
The RNA-Seq data also point to protein degradation as a key mechanism in the acclimation of cold-tolerant fruit to low temperature. Different processes related to proteolysis ('response to misfolded protein' , 'proteasome core complex assembly' , and 'proteasomal ubiquitindependent protein catabolic processes') were specifically induced in cold-tolerant fruit after 14 days at 4°C. The ubiquitin-dependent protein degradation involves first the ubiquitination of the target protein by the ubiquitin enzymes in a multi-step process (E1, E2, and E3 enzymes), followed by the degradation of the modified protein by the 26S proteasome, a large multi-catalytic endopeptidase complex [37]. It is generally accepted that the ubiquitin-proteasome system (UPS) allows cells to respond rapidly to intracellular signals and changing environmental conditions such as drought, salinity, and cold stress [38]. The UPS function in stress responses is most communly accomplished through targeting and degradation of a negative regulator in response to a stimulus enabling the activation of signaling pathways required for tolerance of the perceived stress [37]. In Natura, many genes encoding for different subunits of the 26S proteasome (CUUC61601, CUUC111088, CUUC114894, etc.) together with two E3 ubiquitin ligases (CUUC62721 and CUUC60741) were up-regulated in cold-stored fruit, reflecting the importance of the UPS in facilitating the response to chilling conditions. ' Alternative oxidase' (AOX) is another interesting functional category overrepresented in Natura fruit (CUUC91321, CUUC91320, CUUC91319). The alternative pathway of electron transport limits ROS production in mitochondria and has been described as a mechanism to prevent chilling damage in cold-sensitive species [39]. In other freshly harvested fruit such as sweet pepper and tomato, different treatments that increase chilling tolerance such as methyl jasmonate or salicylate also induced the expression of AOX and correlated with CI resistance [40,41]. Similarly, it is likely that Natura fruit could ameliorate the adverse effect of chilling by an increase in the transcription of the AOX genes.

Molecular mechanisms related to chilling sensitivity in zucchini fruit
The chilling sensitive cultivar Sinatra showed the opposite profile; no enriched GO terms were found in the upregulated DEGs, but 60 different GO categories were identified among cold-repressed DEGs. Many of these GO terms were related to response to biotic, abiotic and endogenous stimuli including 'response to cold' , 'response to salt stress' or 'response to ABA stimulus'. This enrichment with cold-repressed genes indicates that, contrary to what happens in Natura, chilling-sensitive fruit are unable or have a lower capacity to cope with the adverse environmental conditions imposed. In particular, several of the down-regulated genes belonging to these BPs were related to Ca 2+ signaling such as calcium-dependent protein kinases (CUUC118571, CUUC99129), calmodulin like-proteins (CUUC93980, CUUC65863), or calmodulin-domain protein kinases (CUUC118569) (Additional file 6: Table S6), revealing a general repression of Ca 2+ signaling in cold-sensitive fruit. These results suggest that similarly to other species, under cold stress conditions Ca 2+ signaling must be crucial for the acquisition of cold-tolerance [42,43].
It is important to highlight that one of the first symptoms of chilling is related to damage of the plasma membrane, explaining why some measurements such us electrolyte leakage and lipid peroxidation are considered as good biochemical markers of chilling damage in zucchini fruit [9]. It is notable that nearly 18% of specifically cold-repressed DEGs in Sinatra fruit are associated specifically with the plasma membrane. These DEGs are mainly transporters and protein kinases, which may indicate a loss of functionality due to damage in coldsensitive fruit after 14 days of storage at low temperature.

Candidate genes associated with cold tolerance in zucchini fruit
In an attempt to identify genes that may be involved in the mechanisms that control cold tolerance in zucchini fruit, a series of genes that showed differential expression among cultivars were selected. As described previously, energy and carbohydrate metabolism must play an essential role in the tolerance of zucchini to low temperature. From all the DEGs found in the BPs related to maintenance of energy and redox status, the expression of a malate dehydrogenase (MDH) (CUUC107944) was specifically induced in fruit of the cold-tolerant cultivar. MDH reaction is involved in central metabolism and redox homeostasis between organelle compartments [44]. In transgenic apple plants, the overexpression of a cytosolic malate dehydrogenase improved their cold and salt tolerance [45], so it is possible that the activation of this enzyme would increase the redox response in zucchini fruit.
Another candidate gene to regulate zucchini response to chilling is the one encoding for phosphatidylinositol:ceramide inositolphosphotransferase 2-like (CUUC60795), which catalyzes the transfer of the phosphorylinositol group from phosphatidylinositol to phytoceramide, an essential step in sphingolipid biosynthesis. These molecules have been proposed to play important roles in signal transduction, membrane stability, host-pathogen interactions, and stress responses [46]. The expression of this gene was strongly down-regulated in the sensitive cultivar; however, the role of these lipids in the response of zucchini fruit to low temperature is unknown. In Arabidopsis, their accumulation in response to cold is increased in wild-type (WT) plants with respect to the cold-sensitive double mutant sld1sld2, although in the slightly more cold-tolerant mutant atbi-1 the levels of sphingolipids were similar to those of WT [47]. Additionally, the virusinduced silencing in tomato plants of a sphingolipid Δ8 desaturase (SlSLD) involved in sphingolipid synthesis induced severe chilling damage after a 4°C treatment.
Previous studies showed that peptide transport is involved in stress tolerance. It was reported in Arabidopsis that PTR3-like gene was overexpressed in response to abiotic [48] and biotic stress [49], and the AtPTR3 knockout mutant showed higher sensitivity to bacterial pathogen infection and salt stress, suggesting that this transporter protects against biotic and abiotic stress [49]. In zucchini, we also found that the peptide transporter PTR3-A-like (CUUC107903) and PTR2-like (CUUC89270 and CUUC89268) were specifically upregulated by low temperature in the chilling tolerant cv. Natura.
The signal transduction mechanisms that follow plant stress perception are usually triggered or mediated by TFs. Therefore, TF genes are good candidates for regulating chilling tolerance in zucchini fruit. A large number of TFs have been related to cold stress as well as to other unfavorable conditions in different species [50][51][52]. Among them, the MYB76-like TF gene (CUUC97743) from zucchini was only detectable in fruit of the cold-tolerant cv. Natura, and its expression increased drastically after cold storage. MYB TFs are crucial in regulating the network responses to abiotic and biotic stresses [53]. CUUC97743 encodes a polypeptide of R2R3 MYB type that shows a high homology with MYB TFs linked to processes such as epidermal cell differentiation and trichome development in cucumber [54,55].
Other TFs specifically up-regulated in the tolerant cultivar Natura were a DELLA protein GAIP (CUUC105134) and AP2/ERF domain-containing transcription factor (CUUC104934), whereas five ERFs were specifically down-regulated in the sensitive cultivar, including a C-repeat binding factor (CBF) (CUUC92773). APETALA2 (AP2)/ethylene-responsive-element-binding protein (ER EBP) is a large family of TFs unique to plants that have been implicated in plant responses to stresses such as cold and drought [56,57]. The family includes ERF (ethylene responsive factors), and DREB (dehydration responsive element binding proteins) involved in ethylene-related responses. The best known cold regulatory signaling pathway is that mediated by CBF/DREB1 [42], a small subfamily of DREB transcription factors which activates the expression of cold-responsive (COR) genes [58]. The expression of Arabidopsis CBF1, CBF2, and CBF3 is induced shortly after exposure to low temperature, and their overexpression promotes freezing tolerance in Arabidopsis [59] and other species including tomato, rice, and potato [42], Furthermore, CBF/DREB1s appear to induce the accumulation of DELLA (nuclear growth-repressing proteins), which could be responsible for the growth retardation observed in CBF/DREB1s overexpressing plants [60]. Interestingly, a DELLA protein GAIP-coding gene (CUUC105134) was specifically up-regulated in Natura, while a DELLA protein-coding gene (CUUC81379) was specifically down-regulated in Sinatra concomitant with the down-regulation of the CBF, CUUC92773, in this cold-sensitive cultivar.
Regarding C 2 H 2 -type zinc finger transcription factors, ZAT10-like (CUUC61049) showed a very interesting expression profile during zucchini postharvest. In Natura its expression was specifically up-regulated in response to low temperature, while in Sinatra this gene was down-regulated with respect to fruit stored at 20°C, reaching similar expression levels to fresh harvested fruit. Zhu and coworkers [61] also described the expression of a ZAT10-like protein gene that was induced in mandarin after 60 days of cold storage. Arabidopsis plants overexpressing ZAT10 exhibited growth retardation and enhanced tolerance to drought, salt, osmotic, heat, and oxidative stress, as well as photo-inhibitory light [62][63][64]. Although Mittler and coworkers [62] observed that knockout and RNAi mutants of Zat10 were more tolerant to osmotic and salinity stress, it has recently been reported that this gene functions as a positive regulator in osmotic stress tolerance, withZAT10 phosphorylation being required for its function in Arabidopsis [65]. We propose that the transcription factors identified in our study are promising candidate genes for controlling chilling tolerance in zucchini, especially those specifically induced by low temperature in the coldtolerant cultivar Natura (MYB76-like, CUUC97743; AP2/ ERF-like, CUUC104934; and ZAT10-like, CUUC61049).
An important group of genes related to stress sensing and signal transduction also increased their transcription in Natura compared to Sinatra. Among them, five different genes encoding 14-3-3 like proteins showed higher expression in Natura fruit after 14 days of storage at 4°C (CUUC94909, CUUC94908, CUUC94906, CUUC94907, CUUC111254). It has been reported that phosphorylated like proteins play an important role in abiotic and biotic stress response pathways by interacting and modulating the activity of target proteins [66], or they may interact with components of hormone signaling pathways, such as the ABA signaling pathway [67], that is known to be active under temperature and other stresses. In zucchini, the tolerance to cold storage during fruit postharvest has been associated with an increase in ABA synthesis [3], and higher expression of 14-3-3-like proteins could be related to this behavior.

Conclusions
In this work, transcriptomic changes that take place in zucchini after cold storage have been compared in two contrasting cultivars for cold tolerance, Natura and Sinatra. The main response of the cold-tolerant cv. Natura was an induction of the mechanisms common to different stress conditions, whereas that of the cold-sensitive cv. Sinatra was a down-regulation of the same mechanisms. This study also highlights the crucial role of some pathways including carbohydrate and energy metabolism, as well as the regulation of transcription and signal transduction in the acquisition of cold tolerance in zucchini during long-term storage. The data suggest the importance of protein trafficking and degradation in the adaptation of the cold tolerant fruit to low temperature. Among the molecular networks related to chilling tolerance that have been detected by functional analysis of RNA-Seq data, different candidates genes have been selected; these genes could be useful as markers for selection of new lines and hybrids in the current breeding programs of zucchini.

Additional files
Additional file 1: Table S1. Primers pairs used to perform quantitative RT-PCR. (DOC 35 kb)