Metabolomics of dates (Phoenix dactylifera) reveals a highly dynamic ripening process accounting for major variation in fruit composition

Mature fruits

In the present study, 109 unique date varieties (Phoenix dactylifera) from 14 countries were collected in two separate occasions: A first collection took part in 2012 and a second one in 2013. The term variety is used here to describe a distinct phenotypic class of dates and if the same variety was collected from different countries, a different sample ID was assigned to each collected sample per country. Photos of fruits from 14 date samples collected each from a different country can be found in Fig. 1a. With each date sample, a handful of fruits were selected for pre-processing. Each fruit was weighed and the average weight was recorded for each date sample. Two fruits were halved to get a longitudinal and cross sectional view of the pericarp and seed. An international ColorChecker Color-Rendition Chart (ColorChecker Classic, X-Rite, USA) and a 20 cm ruler were positioned along the fruits on a white background under artificial light and a photograph was taken using a Canon Power Shot S100 USA camera loaded on a pre-set tripod. An example photo can be found in Additional file 1: Figure S1. RGB color values were extracted from all fruits showing on a given photo using Matlab libraries and the results were averaged for each color range separately. Readings from color charts from all processed photos were used to calibrate color measurement across the photos. Further phenotype characterization of the date samples consisted of classification into soft, semi-dry and dry types by reference to the literature as well as moisture content measurement of one representative fruit per date sample. Moisture measurement was performed for a random third of the date samples and was based on calculating the percentage of fruit weight-loss following a 116-h incubation in a 105 °C oven. A full listing of all varieties included in this study together with information on their country of production, collection point and type can be found in Additional file 2. Summary statistics for each sample collection including the count of varieties, samples and the frequency of samples per country of production are shown in Table 1-A. Overall, dates from the first sample collection were mostly from the Gulf region obtained in a fairly dried condition from shops and festivals whilst the second sample collection was dominated by North African dates obtained mostly fresh from the palm trees. For the second collection of dates, field work permissions were obtained verbally from owners of visited oases. The marketed versus fresh nature of dates between the two sample collections implies varying post-harvest conditions. All collected dates with homogenous brown color were further dried by exposing them to open air for two weeks before further processing. In general, dates were considered mature if the low moisture prevented any further change in their appearance. Notably, maturity is attained naturally with the dry class of dates but often artificially with the soft class of dates owing to intrinsically higher moisture levels (refer to background for further details).

A) Overview of the date cohorts from the first and second sample collection. Countries are denoted by their international ISO Alpha-2 code as follows: MA (Morocco), DZ (Algeria), TN (Tunisia), LY (Libya), EG (Egypt), SD (Sudan), JO (Jordan), SA (Saudi Arabia), IQ (Iraq), IR (Iran), AE (United Arab Emirates), QA (Qatar), PK (Pakistan), US (United States). Dates from the first sample collection contained largely mature dates. In contrast, 37 immature date fruits corresponding to 10 varieties were included in the second sample collection. B) Summary statistics of the metabolomics data measured from the first and second sample collection. To account for batch effect, 10 samples from the first sample collection were measured again along dates from the second sample collection. Dates from the second sample collection were only measured by Metabolon unlike dates from the first sample collection which were measured by MetaSysX and Metabolon. The median RSD (RSD = sdandard deviation/mean) from biological replicates is a combination of technical and biological variation whilst that from technical replicates only expresses technical variation

Immature fruits

Dates preprocessing and measurement protocols

The metabolic content of the date fruits from the second sample collection was measured separately a year after samples from the first collection were measured. The first collection of dates was preprocessed by MetaSysX GmbH. and measured by both MetaSysX GmbH. and Metabolon Inc., USA. Dates from the second collection were preprocessed and measured by Metabolon Inc., USA alone. The protocols for sample processing and metabolomics measurement by both MetaSysX and Metabolon are described in details in Additional file 3.

Briefly, with MetaSysX, 50 mg of the peel and flesh of the date fruits were flash frozen in liquid nitrogen and extracted according to standardized procedures [23]. The dried metabolite extracts were measured with a Waters ACQUITY Reversed Phase Ultra Performance Liquid Chromatography (RP-UPLC) coupled to a Thermo-Fisher Exactive mass spectrometer which consists of an ElectroSpray Ionization source (ESI) and an Orbitrap mass analyzer. C8 and C18 columns were used for the lipophilic and the hydrophilic measurements, respectively. Chromatograms were recorded in Full Scan MS mode (Mass Range [100–1500]) [23]. Chromatograms from the UPLC-FT-MS runs were analyzed and processed using the software REFINER MS® 7.5 (Genedata, Switzerland). The data were further filtered and analyzed using in-house software tools (refer to Additional file 3). The samples were also measured using the Agilent Technologies GC coupled to a Leco Pegasus HT mass spectrometer which consists of an EI ionization source and a TOF mass analyzer. Column: 30 meters DB35; Starting temp: 85 °C for 2 min; Gradient: 15 °C per min up to 360 °C. NetCDF files exported from the Leco Pegasus software were imported into “R”. The Bioconductor package TargetSearch was used to transform retention time to retention index (RI), to align the chromatograms, to extract the peaks and to annotate them by comparing the spectra and the RI to the GMD [24, 25]. Obtained data from both platforms were normalized according to sample weight and to the measurement day to minimize process error over the course of many days of measurement.

With Metabolon, date samples were prepared and extracted according to the standard solvent extraction method by Metabolon Inc. [26]. The UPLC/MS/MS analysis was based on the Waters ACUITY ultra performance liquid chromatography (Waters Corporation, USA) and the ThermoFischer Scientific Orbitrap Elite high-resolution accurate mass spectrometer (Thermo Fischer Scientific Inc., USA) equipped with a heated electrospray ionization (HESI) source and an Orbitrap mass analyzer. The dried sample extracts for the LC positive and LC negative mode were reconstituted in acidic and basic LC- compatible solvents. Two independent injections were performed on each sample using separate dedicated columns. The mass spectra analysis alternated between MS and data dependent MS² scans using dynamic exclusion. With GC/MS, the samples were further dried under vacuum desiccation for an entire day and derivatized under dried nitrogen using bistrimethyl-silyl-trifluoroacetamide (BSTFA). The GS/MS analysis was based on a Thermo Finnigan™ TRACE™ DSQ™ (ThermoFinnigan, USA) fast-scanning single –quadrupole mass spectrophotometer using electron impact ionization source. The GC column was 5 % phenyl and the temperature ramp range was from 40 to 300 °C in a time span of 16 min. The raw data files from both platforms were extracted using the in-house informatics system (refer to Additional file 3). A reference library maintained by Metabolon Inc. [27], consisting of chemical standards with retention time, retention index, mass to charge ratio (m/z) and chromatographic data including MS/MS spectral data was used to identify metabolites in experimental samples as detailed in [28]. In this study, the samples were analyzed over a span of two or three days, and therefore data normalization step was performed to correct variation from instrument inter-day tuning differences.

Measurement experimental design

With the first collection of dates containing 62 date samples, the MetaSysX measurement was done in triplicates yielding a total of 186 measured metabolic profiles (Table 1-B). With Metabolon, 34 samples were measured in duplicates whilst the 28 remaining as singletons, amounting to 96 measured metabolic profiles (Table 1-B, Fig. 1c). For the rest of this article, we will refer to the latter as ‘DS1-bolon’ whilst the former metabolomics dataset will be referred to as ‘DS1-sysX’. Dates from the second sample collection were measured by Metabolon only and therefore the derived metabolomics data will be referred to in short as ‘DS2’. DS1-bolon and DS2 metabolomics data can be found in Additional file 4 & Additional file 5 respectively. The experimental design consisted of a singleton measurement of each of the 51 mature date samples (Table 1-B, Fig. 1c) and similarly the 37 immature fruits were each measured once. To account for batch measurement effect, 10 fruits from the first sample collection were measured again along the 88 fruits from the second collection, resulting in 98 measured metabolic profiles (Table 1-B). We distinguish between metabolomics data from the 37 immature and 61 mature date samples (inclusive of the 10 samples from the first collection) using the terms ‘DS2-immature’ and ‘DS2-mature’ respectively (Fig. 1c). The sample characteristics of DS1-sysX, DS1-bolon and DS2 as discussed here are summarized in Table 1-B. Since Metabolon measured datasets were extensively used in this paper, they are further illustrated in Fig. 1c.

Data preprocessing and platform comparison

Metabolomics data, were log-transformed and scaled so that the median measurement value from each measured metabolic profile was equal to the overall median from the whole dataset. This normalization was done separately for DS1-sysX, DS1-bolon and DS2. By default, biological replicates (when available) were not combined and measurement from each replicate was treated as a separate metabolic profile. However, with few analyses, a single measurement from each date sample was required and the replicates were averaged. This will be clearly indicated where applicable. Comparison of platforms was based on average metabolite missingness level across samples and the median relative standard deviation (RSD) across biological replicates. RSD was expressed as metabolite-wise standard deviation from replicates divided by the mean. With Metabolon measurement of samples from the first collection (or DS1-bolon), data from technical replicates were available from repeated measurement of a homogenous mixture of pooled samples (refer to Additional file 3). The median RSD from these technical replicates was calculated for assessment of data quality by Metabolon.

Non-supervised PCA analysis of mature dates and quality control

The multivariate statistical analysis package SIMCA v13.0.3 was used to perform PCA on DS1-bolon, DS1-sysX and DS2-mature separately to characterize collective metabolic variation underlying significant proportions of the variance from the respective datasets. Simca default metabolite missingness threshold of 50 % was used [29]. The significance of the extracted principal components was derived from SIMCA via built-in cross validation where for each component consecutively, parts of the data are alternatingly kept out of the model then predicted [29]. Based on the PC1/PC2 two dimensional space, date samples 78-BZGZ-MA and 105-ZGHL-EG from DS2-mature located outside the Hotelling’s 95 % confidence ellipse interval were considered outliers and excluded from further analysis of the dataset [29].

SIMCA OPLS-DA and O2PLS-DA models of the dates ripening process

Metabolic signature of date ripening was modeled from analysis of the development stage dataset, or DS2-immature, a subset of the second date sample collection as follows: Initially, PCA analysis was run on measured metabolomics data to confirm the within-sample ranking of individual fruits previously set by visual assessment of the fruits’ extent of ripening (refer to the previous section). The PCA analysis revealed clusters of fruits with comparable ripening profiles across samples (more details in the results section). These clusters were used to define development stage classes that served as a training set for an OPLS-DA classifier [29, 30]. Applying the classifier on the rest of the samples in DS2 led to the calculation of class prediction scores indicative of the samples’ ripening metabolic states. For DS1-bolon, the OPLS-DA model trained on DS2-immature data was not suitable owing to likely differences between batch measurements. Also, unlike the second collection of dates, no development stage dataset was included in the first collection. Instead, we developed a strategy based on the 10 fruits from the first sample collection which were measured again along the samples from the second collection. Because the samples in question were included in both batch measurements, they will be referred to as batch 1&2 samples for the remaining parts of this article. Our strategy for predicting the ripening states of dates from the first sample collection is here described: First, we used the OPLS-DA model previously trained on the DS2-immature samples to predict the development classes of batch 1&2 samples based on their DS2 data from the same batch measurement as the training set. This class information was used to train an O2PLS-DA classifier on the same samples (batch 1&2 samples) based on their batch 1 and 2 metabolomics measured data. The O2PLS-DA procedure [29, 30] is able to identify metabolites consistently differentiating between the different classes in the training set based on multiple measurements of the training set (here from different batch reading). The integrative nature of the O2PLS-DA model meant that it could be used to calculate class prediction scores for dates from the first and second sample collection. The scores from the first sample collection served to indicate the ripening states of these date samples whilst the scores from the second collection served to optimize and validate the O2PLS-DA model by drawing a comparison to the class prediction scores for the same samples by the original OPLS-DA model (more details in Additional file 1: Figure S2). The O2PLS-DA model was only defined on Metabolon measured data.

Association analysis of PCs from mature dates with date (soft/dry) type, country of production, ripening state and color

The lm function from the statistical analysis R software version 3.1.1 was used to run the regression model ‘PC ~ date_variable’ where date_variable consisted of one of four variables: date_type, a categorical variable with two levels: Soft and dry, with semi-dry varieties assigned to the dry class (Additional file 2); date_country, an ordinal variable from ranking the sampled countries West to East; date_ripening_state corresponding to the class prediction scores calculated by the OPLS-DA and O2PLS-DA models for samples from the first and the second collection respectively and date_color, a continuous variable based on the average of the red/green/blue (RGB) color measurements. The R package maps was used to generate the geographical map in Fig. 2 depicting the dates countries’ of production.

Analysis of the distribution of classes of metabolites on the loading space underlying PCs from mature dates

In order to further characterize PC1, the distribution of metabolites classified into broad metabolic categories including amino acid metabolism, sugar metabolism, energy metabolism, lipid metabolism, purine and pyrimidine metabolism, secondary metabolism and vitamin metabolism was manually examined on the underlying loading value space. The latter refers to the set of loading values assigned to the metabolites by PCA analysis where each loading value expresses the correlation between the corresponding metabolite abundance profile and the PC scores. Within a broad metabolic class, sets of metabolites sharing a functional or structural feature and having comparable loading values were identified. The common feature consisted mostly of pathway co-membership, a common catalytic activity or a unifying structural theme. These sets of metabolites were mapped to subclasses within the original broad categories as follows:

Finally, a general category degradation activity and amino acid volatiles (VOC) was formulated to capture metabolites from degradation of purines, vitamins and amino acids leading to synthesis of short chain volatiles (VOCs) [8]. For the rest of the article, all afore mentioned subclasses of metabolites as well as unrefined categories vitamin metabolism, hormone metabolism, energy metabolism and degradation activity and amino acid VOC will be collectively referred to as ‘metabolite classes’. It is important to note that the analysis was restricted to Metabolon measured data.

Amino acids and related metabolites

Enrichment in free amino acids was observed in this study in dates with an early ripening profile similar to other fruit [8, 39] (Fig. 5). Amongst all detected amino acids, the levels of alanine, glutamate and aspartate declined least in dates with late ripening profile (Additional file 6), consistent with previous work looking at ripening in tomato [39]. In general, amino acids serve as building blocks for synthesis of key intermediates and end-products of the ripening process in fruits [8]. In particular the aromatic amino acids, also measured in this study, give rise to a myriad of secondary metabolites, notably color and flavor-conferring phenylpropanoids. The observed enrichment in dipeptides in date fruits with an early ripening profile in this study (Fig. 5) may be linked to protein degradation activity recruited by hormones to eliminate pre-ripening enzymes at the onset of ripening [8]. Another potential source for the dipeptides is the targeting peptide sequence, attached to pre-folded nuclear proteins, which is digested upon protein entry into organelle structures including chromoplast [40]. Chromoplasts, which are differentiated forms of chloroplasts lacking chlorophyll, act as metabolic hubbs at early ripening, necessicating constant in-flow of effector proteins from the nucleus [41]. Targeting peptide sequences contain mostly hydrophobic amino acid residues [40], consistent with the high proportion (70 %) of hydrophobic valine, leucine, isoleucine and phenylalanine amongst amino acid constituents of the dipeptides observed in this study. Interestingly, N-acetylation of chromoplast-targeted pre-folded proteins was suggested as a mechanism of organelle specificity [42]. This may account for the enrichment of N-acetylated amino acids in date samples with an early ripening profile in this study (Fig. 5); although acetylation can sometimes be a necessary intermediate reaction during metabolism of amino acids. Finally, enrichment in glutathione activity in dates with an early ripening profile (Fig. 5) is consistent with increased antioxidant activity in fruits at early ripening [43].

Primary amines and polyamines

In this study, ethanolamine, GABA, serotonin, tyramine, tryptamine and phenethylamine from the decarboxylation of serine, glutamate, 5-hydroxy tryptamine, tyrosine, tryptophan and phenylalanine were all found enriched in dates with early ripening profiles (Additional file 6). This is consistent with early ripening expression of amino acid decarboxylases leading to amine synthesis in a number of fruits [44, 45]. Tyramine and tryptamine serve as precursors for the synthesis of defense-mediating alkaloids previously detected in dates [46] and a turnover of phenethylamine into antipathogen volatiles phenylacetaldehyde and phenylethanol serves the same purpose [8, 45]. The role of serotonin in fruit ripening has not been fully investigated; however, melatonin that derives from N-acetylserotonin (a derivative of serotonin also observed in our data), has recently been found to promote various physiological aspects of ripening when given exogenously to green tomato [47]. Recently, the decrease in GABA with ripening was linked to maintaining high levels of essential glutamate and aspartate during tomato ripening [48]. Enrichment in the polyamine putrescine in dates with an early ripening profile (at the positive end of PC1) is consistent with previously reported expression of a mouse ornithine decarboxylase conjugated to a ripening specific promotor at the onset of ripening in transgenic tomato [49]. Previous work suggested a synergy between the ripening hormone ethylene and the polyamines spermine and spermidine, derivatives of putrescine [39, 50]. The potential regulatory role of putrescine may justify its co-occurrence with products from its degradation pathway in dates with an early ripening profile (Additional file 6).

Secondary metabolism

The earliest sign of secondary metabolites from the phenylpropanoid pathway in our data consisted of tannins procyanidin B1 and procyanidin B2 and catechin monomers all showing maximal level in dates with an early ripening profile (Additional file 6), in accordance with the literature [8, 51]. Astringent tannin oligomers are abundant in green fruit and only lose their astringency when undergoing structural changes as ripening progresses [52]. In addition to tannins, a wide range of color and flavor flavonoids and hydroxycinnamates were observed to peak at different ranges of PC1 and some showed no correlation with PC1. In general, variance from these metabolites was poorly explained by PC1 (Additional file 6). This could be due to a much stronger influence by genetic background [8], which may contribute to unique taste and color characteristics of individual varieties. Interestingly, PC4 from DS2-mature revealed an opposing trend between classes phenylpropanoids and TCA on one hand and the accumulation of phosphorylated sugars, captured under the class glycolysis, on the other hand. This effect can be explained by energy requirement for the synthesis of phenylpropanoids through initial degradation of phosphorylated sugars during glycolysis and downstream TCA activity. A third class of secondary metabolites consisted of volatiles, major contributors to aroma in fruits. In this study, an increase in branched chain amino acid derived volatiles and hydroxycinnamate derived volatiles was observed in dates with a late ripening profile (Fig. 5 and Additional file 6), consistent with the literature [8]. Volatiles are strong attractants of seed dispersers and their sharp increase in overripe fruit could constitute a mechanism to maximize the chance of consumption before onset of senescence.

Changes to cell wall and cell membrane

Alterations in cell membrane composition have long been known to occur during ripening [53–55], but have been given little attention by the more recent literature. Key changes to membrane phospholipids during ripening include an increased desaturation level of fatty acyl chains facilitating their peroxidation. Induction of expression of a handful of desaturase isomers was found to be strongly associated with a continuous flux of linoleate and linolenate substrates of the lipoxygenase (LOX) pathway during peach ripening [56]. This pathway is known for being the mechanism of synthesis of a myriad of C6 volatile aldehydes and alcohols that contribute significantly to fruit aroma [8]. In this study, a range of mono and poly-unsaturated fatty acids including the LOX substrates and their oxidized derivative oxylipins were observed to generally plummet in dates with a late ripening profile at the negative range of PC1.

Initial excision of fatty acyl chains resulted in accumulation of lysophospholipids in dates with an early ripening profile and late accumulation of lysophospholipid degradation products in dates with a late ripening profile (Fig. 5). These included monoacylglycerols and phosphorylated head groups, free head groups and remaining lysophosphatidic acids [57] as well N-acylethanolamines, derivatives of lysophospholipids via N-acylphosphatidylethanolamine intermediates [31] (Additional file 6). These degradation products, in addition to sphingoid bases, are likely to have a signaling role in fruits and some are already known to be downstream mediators of abscisate signaling in other plant organs [58–60]. The opposing trend in lysophospholipids versus sphingoids by PC2 in both datasets is interesting and may suggest a change in signaling patterns during ripening of certain date varieties or as a response to certain external stimuli. For instance, sphingolipid signaling is known to be induced under drought conditions in plants [59].

Fruit softening during ripening is known to be associated with an increased activity of fiber degrading enzymes, many of which targeting cell wall structures [8]. Pectin is a major constituent of the plant cell wall and its hydrolysis seems to make way for synergistic disassembly of cellulose and hemicellulose polysaccharide matrices [8]. In this study, the abundance of keto-deoxyoctulosonic acid, an acidic monosaccharide located in the side chain of the rhamnogalacturonan II class of pectin peaked in dates with middle range PC1 (Fig. 6b). This served to establish a link to the Rutab phase which is characterized by maximal fruit softness (refer to results).

Sugar metabolism, energy and gene expression activity

In this study, sucrose and related non-reducing sugars kestose and melezitose (from DS1-bolon) were most abundant in fruits with an early ripening profile at the positive range of PC1 (Fig. 5, Additional file 6). In fruits, an increase in sucrose levels is observed at the late mature green stage and sucrose is broken down by the enzyme invertase following the onset of ripening [8]. Metabolism of released sugar monomers proceeds via the early pentose phosphate pathway, which serves to provide carbon precursor molecules for the synthesis of aromatic amino acids, secondary metabolites, vitamins and purine/pyrimidines. In parallel, flux through glycolysis and the downstream TCA cycle serves to sustain energy levels required for gene expression and anabolic reactions during ripening. The relationship between precursor non-reducing sugars and product TCA intermediates is captured by PC3 from DS1-bolon and may reflect varying glycolysis/TCA kinetics across the date cohort.

In this study, the abundance levels of energy molecules NAD+ and AMP were lowest in dates with a late ripening profile (Fig. 5, Additional file 6). This together with accumulation of the TCA cycle intermediate fumarate, members of the glycolysis pathway and ribosomal nucleosides possibly originating from ribosome degradation (Additional file 6) may indicate a diminished metabolic activity at this late stage of ripening. Interestingly, metabolites from the TCA/rRNA nucleosides classes appeared to contrast the unsaturated fatty acids and oxylipins in PC4/PC3 from DS1-bolon/DS2-mature respectively. Late synthesis of oxylipins in some species of dates may require residual TCA activity and late synthesis of key enzymes.

The accumulation of reducing sugars in dates with a late ripening profile in this study (Fig. 5) confirmed similar reports in dates at the late Tamr stage [10]. Amongst the sugars observed (also shown in Additional file 6) xylose, fucose, arabinose and glucose may derive from cell wall hydrolysis activity during ripening. Ribulose and xylulose could derive from their phosphorylated forms from the pentose phosphate pathway whilst sugar alcohols, sugar lactones and derivative acids may result from upregulation of aldose/ketose oxidoreductases in ripe dates, previously shown using a proteomics approach [61]. It is important to note that besides contributing to fruit flavor, sugar alcohols are a type of polyol osmoprotectant that may act to alleviate the impact of fruit dryness at this stage.

Vitamins and hormones

In this study, a range of vitamins has been detected including riboflavin, niacin, pyridoxine and nicotinate (Additional file 6). The enrichment in pyridoxine in dates with an early ripening profile may be attributed to its essential role in amino acid synthesis and metabolism by amino acid decarboxylases at the early phase of ripening. In contrast, the accumulation of threonate (Additional file 6), a degradation product of vitamin C, in dates with a late ripening profile is consistent with the previously described decrease in vitamin C at the Tamr stage [1].

Dates behave like climacteric fruits implying a leading regulatory role by ethylene although interplay with abscisate may be operating on certain ripening events [8]. In this study, the ripening hormone ethylene was not measured (below the mass cut-off imposed on the instrumentation) but two related metabolites were observed (Additional file 6). One is cyanoalanine, which is a conjugate of cysteine and cyanide (cyanide being a toxic byproduct of ethylene synthesis [62]) and 5-methylthioadenosine, an intermediate in the Yang cycle that replenishes the ethylene precursor SAM. Both molecules showed maximum levels in dates with an early ripening profile at the positive range of PC1 (Additional file 6). This range was previously mapped to the Khalal stage (refer to results), which follows the climacteric ethylene peak in dates [63]. A similar pattern by isopentenyl adenosine, the precursor of the cytokinin hormone zeatin, is concordant with the reported abundance profile of zeatin during tomato ripening [64]. Abscisate is an early regulator that precedes ethylene in the climacteric fruit model tomato [65, 66]. Following the peak in ethylene, abscisate was shown to increase in abundance slightly [67], which could explain the small peak in abscisate observed in this study in dates with middle range PC1 values. In summary, key hormones and related metabolites were measured as part of this work motivating future analysis of partial correlations with other measured metabolites in order to uncover regulatory mechanisms operating on distinct aspects of the ripening process in date fruit.