Data files were imported into Acuity 4.0 for normalization and analysis. Only spots with intensity values higher than the background plus two standard deviations of the background median, in at least one channel, were used for analyses. Before normalization, the median local feature background was subtracted. Data were normalized by Lowess with a centered print-pin tip using the Acuity default values. To generate the raw data to be used for the expression analysis, a Lowess M Log Ratio was used as the expression value, and patterns with more than 80% of non missing values were selected. In all, 3350 probes met the threshold for hybridization quality . The data sets supporting the results of this article are available in the Array express repository, [E-MEXP-3902].Differentially expressed genes were identified from the raw dataset using the Significance Analysis of Microarray software. Missing values were imputed by 10-Nearest Neighbors Imputer algorithm, with 100 blocked permutations and a random seed value set by default in the program. PCA and 2D-hierachical cluster analyses were performed on the significant data using Acuity . A principal component analysis was calculated for those factors explaining 100% of variance. For calculations spots with missing values were replaced with the average values across the arrays.

Profiles with the same shape pattern were centered around the mean value across arrays, large pot with drainage to avoid the effect that the magnitude of response might have on the average profile. For the hierarchical cluster, a Pearson correlation centered on 0 was used as a similarity metrics. A complete linkage was used to link clusters together to produce the tree. Transcripts and/or samples were ordered in the clusters according to their contribution to principal component 1 of the PCA performed with the same dataset. The ChillPeach genes were classified into 34 distinct functional categories and 702 specific processes by extensively reviewing the literature and by searching in reference databases . Functional enrichment on a ranked list of genes was performed with, a local, customized version of ‘catscore.pl’ Perl script described in Cheung et al., using a two-tailed Fisher exact t-test with adjusted p-value cut-off of 0.05.Correlations were calculated by the Pearson product moment correlation method using Matlab 2007 . P values below 0.01 were selected for statistical significance. A statistical significance level of 1% was assessed with the correlation coefficients over 0.8. Those genes whose expression profiles contained 100% of data points in the samples analyzed were used to calculate correlations. The complete list of the microarray wide gene expression correlations with the Mealiness Index are listed in Table S3. Functional enrichment is performed as indicated above. The 96.96 dynamic arrays were obtained from the Fluidigm Corporation and were used to set up four sets of qRT-PCR reactions of 64 cDNA preparations corresponding to 32 samples: 15 genotypes in the M stage and/or CS1 samples and 5 pools .

Two biological replicates were included in each array for all the 15 genotypes and pools, each one representing at least three different fruits. Two replicated 96.96 Fluidigm dynamic arrays were used. For the Fluidigm analysis, 72 genes were selected from our microarray results . Oligo pairs for selected genes were obtained using the Primer Express version 2.0 software . To design primers, the following conditions were used: Tm 58–60uC, GC content 20–80%, primer length 20–22 base pairs and an amplicon size of 140–150 bp. A virtual PCR was carried out for each oligo pair obtained with the ‘primersearch’ program from the EMBOSS open software suite, using the full set of known peach sequences as potential template sequences. The interrogated peach sequence databases included the ChillPeachDB, ESTreeDB and GDR_Prunus sequences. Only the oligo pairs yielding a single PCR product from each unique gene, based on the sequence assembly of all the known Prunus sequences, were considered. When more than one specific oligo was obtained for a gene, the oligo pair with the lowest penalty value , and which mapped most of the 39 end of the gene, was selected using custom Perl scripts. Three genes were selected to normalize qRT-PCR results on the basis of low variability in the chillpeach microarray under all conditions analyzed in this paper: a gene with unknown function , an ABC1 family protein and, an esterase/lipase/thioesterase gene They were validated by qRT-PCR as described in [45]. The comparative DDCt method, as described by in Livak and Schmittgen, was used to confirm a flat pattern throughout the samples. For the Fluidigm analysis, the cDNA synthesized from total RNA following standard methods was diluted to 1:10 using the DA Assay Loading Buffer .

The Nanoflex 4-IFC Controller and the BioMark Real Time PCR system by the Fluidigm Corporation were used to run the dynamic arrays under the standard conditions employed at the General Hospital lab, Valencia, Spain. The cycling program consisted of 10 min at 95uC followed by 40 cycles of 95uC for 5 sec and 1 min at 60uC. PerlqXpress was used to calculate ‘‘fold expression values’’ from the Ct values obtained directly from the BioMark Real-Time PCR Analysis Software . Briefly, PerlqXpress filter outliers within a sample, corrected differences in background control levels, centers and scales data. The mean centered and scaled Ct values were transformed into relative quantities using the exponential function with the efficiency of PCR reaction as its base. For each gene the RQ was corrected using a normalization factor. FC is calculated by dividing normalized RQ to reference sample in each biological replicate . Mean, standard deviation, and coefficient variation were calculated for each replicate. Replicates were filtered by the coefficient variations. At least 4 good replicates were used to calculate ‘‘fold expression change’’ values. To extend the validity of the results obtained in the pools to individual lines, for which we had individual MI index values, qRT-PCRs were performed on 15 individual peach genotypes from the popDG progeny . For each gene pair in a predefined expression set, the Pearson correlation coefficients between their expression profiles in the individual Pop-DG siblings were obtained by Gitools 1.8.2. A gene was selected as consistent and was confirmed over the individual lines when it correlated with a predefined expression pattern.Harvest maturity, a factor known to influence mealiness, was tested before cold treatments to ensure all fruits were in the same maturity stage. Table S1 shows there were no significant differences in firmness, SSC and TA between genotypes. This indicates that at harvest, both populations were at the same physiological stage and differences in the subsequent cold response can be mainly attributed to the cold sensibility without significant distortions owing to lack of adequate maturity stage. To asses the effect of the cold stress on peach fruits from siblings of the Pop-DG population, a subset of the cold stored fruits were ripened for 2–3 days at 20uC and mealiness was evaluated as the proportion of measured fruits with mealiness or Mealiness index . Figure 1A shows the average MI of pools of fruits grouped according to their sensitivity to develop mealiness. The pool S had higher MI as compared with pool LS after the same cold storage times , although tend to converge after increasing cold storage periods, drainage collection pot indicating a non complete tolerance of fruits LS. The difference was more pronounced after one week of cold storage at 5uC, where the mealiness symptoms were already visible in the pool S but not the pool LS . No significant differences in the frequency of other CI symptoms were observed between pools S and LS . Thus the characteristic feature, differentiating the cold response of the pools, was their sensitivity to develop mealiness. Given that the proportion of mealy fruits increased with the time of cold storage, our hypothesis is that despite mealiness was not showing until fruit was allowed to ripe, relevant molecular changes may had already started to occur during cold storage.Bulked segregant analysis in combination with the Chillpeach expression microarray was used to compare the transcriptomes of peach fruits from the S and LS Pop-DG siblings. In all, 3350 transcripts showed a significant difference in expression levels at least for one condition using two criteria: a false discovery rate ,5% and a p-value ,0.05.

The principal component analysis of the complete dataset variance is seen in Figure 1B. PC1 clearly separated fruit samples which came directly from cold storage , from fruits M and R . The proximity between fruits M and R, if compared to CS, indicated that the effect of cold storage on the peach transcriptome was much greater than that of ripening. PC2 separated fruits M from R. Both, fruits S and LS seemed to follow parallel ripening programs, but fruits LS showed delayed or less intense ripening transcriptomic changes than fruits S. It should be noted that Pop-DG siblings in each pool were selected on the basis of their cold response, as revealed by the MI after shelf life ripening, so it is not surprising that some differences in the ripening programs may exist. In addition, PC2 roughly separated cold stored samples according to the eventual increase in the MI of the fruit should they be submitted to shelf life after cold storage . According to this component, fruits from the pool S stored for 1 or two weeks have achieved a pattern of ripening similar to fruits R . This may indicate that during cold storage at 5uC, some internal ripening may result in chilling sensitivity and in a shortened shelf life. The loading plots for PC2 revealed 37 genes among with were genes previously reported in the regulation of temperature responses, including the transcriptional factor CBF, GASA5 and SCR2. Thus the transcript levels contributing to component PC2 may be relevant for the development of a tolerance mechanism in cold, which could affect the way cold storage interrupted or slowed down the ripening program and eventually how fruits ripen afterward. The bidimensional hierarchical cluster analyses revealed a similar sample separation to that obtained with PCA . Furthermore, 2D_HCA segregated CS1-LS from the rest of the cold-stored fruits , according to the fact, that if fruits CS1-LS ripen, they do not develop mealiness. These results indicate that from the molecular point of view one week of cold storage is a critical time with maximum differences expected to be found at this stage between fruits S and LS, including any CI associated trait. This global analysis also revealed that, although the expression profiles were generally similar between the S and LS pools of fruits, there were qualitative and quantitative differences . To further describe the cold response mechanism from a global point of view and its possible relation to eventual CI, we conducted a functional enrichment analysis of the 11 resulting clusters from 2D-HCA . The most over represented functional category in cluster CS-glob8, containing genes up-regulated during cold storage in both fruit pools, was RNA transcription regulation, which comprised 47 genes . In this category, we found several transcription factors whose orthologs were up-regulated during cold acclimation in Arabidopsis and some were assigned to specific cold acclimation regulons . The other functional category enriched in CS-glob8 was with 37 genes, secondary metabolism, a functional category previously associated with cold tolerance. In addition, and in agreement to the higher tolerance of fruits LS, structure maintenance proteins and an antioxidant system were among the functional categories over represented in differential clusters CS-glob7 and CS-glob 10 , both highly induced in the pools of fruits LS as compared to fruits S . Moreover, cluster CS-glob 9 was enriched in RNA translation and protein assembly, energy production, and trafficking machinery and membrane dynamics , indicating that these activities can be enhanced in fruits LS. This suggests that some kind of cold adaptation was activated in both S and LS peach fruits during cold storage. The genes in cluster CS-glob 2 were down-regulated in both S and LS fruits , and were enriched in glycolysis/pentose phosphate pathway, the photorespiratory pathway and organelle division . Lowered expression levels of carbohydrate metabolism and glycolytic genes correlated to cold sensitivity in Arabidopsis. However, the extensive down-regulation of primary metabolism, together with the down-regulation of posttranscriptional, posttranslational and protein degradation , was probably associated with the relative higher cold tolerance of fruits LS.