Multi-omics data integration provides insights into the post-harvest biology of a long shelf-life tomato landrace

Abstract In this study we investigated the transcriptome and epigenome dynamics of the tomato fruit during post-harvest in a landrace belonging to a group of tomatoes (Solanum lycopersicum L.) collectively known as “Piennolo del Vesuvio”, all characterized by a long shelf-life. Expression of protein-coding genes and microRNAs as well as DNA methylation patterns and histone modifications were analysed in distinct post-harvest phases. Multi-omics data integration contributed to the elucidation of the molecular mechanisms underlying processes leading to long shelf-life. We unveiled global changes in transcriptome and epigenome. DNA methylation increased and the repressive histone mark H3K27me3 was lost as the fruit progressed from red ripe to 150 days post-harvest. Thousands of genes were differentially expressed, about half of which were potentially epi-regulated as they were engaged in at least one epi-mark change in addition to being microRNA targets in ~5% of cases. Down-regulation of the ripening regulator MADS-RIN and of genes involved in ethylene response and cell wall degradation was consistent with the delayed fruit softening. Large-scale epigenome reprogramming that occurred in the fruit during post-harvest likely contributed to delayed fruit senescence.


Introduction
Tomato (Solanum lycopersicum L.) has been widely used as a model system to study climacteric f leshy fruit development and elucidate the regulatory network underlying the fruit ripening process that is orchestrated by transcription factors, hormonal regulation and epigenome modifications [1]. However, it is still unclear whether fruit ripening should be regarded as a senescence process or not [2]. A better understanding of this topic can be highly beneficial for many reasons. Primarily, improving fruit shelf-life could alleviate postharvest losses that occur worldwide. Additionally, postharvest handling practices, such as refrigeration of unripe and ripe fruits, can adversely affect their quality, including nutritional and sensory attributes [3].
A reservoir of alleles for shelf-life-related traits is available in the tomato germplasm, especially in the Mediterranean long shelf-life landraces [4]. For instance, landraces from Spain, namely "de penjar" and "de Ramellet", show a delayed deterioration of the fruit associated with the alcobaça (alc) mutation in the coding sequence of NONRIPENING (NAC-NOR), a transcription factor that controls the expression of ripening-related genes [5]. In southern Italy, several long shelf-life landraces, collectively named "Pomodorino del Piennolo del Vesuvio DOP" (hereinafter referred to as "Piennolo del Vesuvio"), are traditionally grown in the area of the Vesuvius volcano according to the traditional agronomic practices established by the Protected Designation of Origin regulation. The typical preservation technique is based on the production of "piennoli": tomato bunches are harvested during the summer and then, thanks to the jointless pedicel trait, they are put on a hemp string to make a large bunch weighing 2-3 kilograms (Fig. 1a). Each "piennolo" is hung in a ventilated and dry warehouse and fruits remain firm for a period that can extend up to ten months.
The genome of a "Piennolo del Vesuvio" landrace named "Lucariello" (hereinafter referred to as LUC) has been resequenced to explore nucleotide variability associated with the long shelf-life trait [6]. Although sequence variations in ripening-related genes were identified, no clear relationship was found between genetic variants and long shelf-life.
A tomato fruit transcriptome atlas and epigenetics studies have advanced understanding of different aspects of tomato fruit biology [1,7,8]. Methylome analysis and treatment with DNA methylation inhibitor 5-Azacytidine demonstrated that DNA methylation is crucial for controlling tomato fruit ripening [9]. An additional level of ripening regulation is introduced by post-translational modifications of histones [8,10]. Indeed, it has been established that H3K27me3 is a conserved epigenetic mark associated with fruit ripening [11]. Supplementary regulation of ripening is performed by non-coding RNAs [12], including microRNAs (miRNAs) that have been extensively profiled in tomato fruits [13]. It is well known that transcript degradation and translation repression mediated by miRNAs affect key ripening genes [14]. Notably, less attention has been paid to the molecular mechanisms underlying events of tomato post-harvest being difficult even to separate fruit ripening and senescence processes [2].
As LUC is a long shelf-life landrace with available genome sequence, we adopted a multi-omics data integration approach in the attempt to fill this research gap. Our results allowed us to unveil global changes in the expression of protein-coding genes and miRNAs as well as genome methylation and histone modification dynamics, thus elucidating the integrative transcriptome-epigenome landscape in distinct phases of post-harvest.

Fruit quality evaluation
LUC was characterized by oval-shaped fruits with a pointed apex, 2-3 loculi and weighing less than 20 g. Fruit ripening was apparently regular, as suggested by changes in the colour of pericarp reflecting the accumulation of carotenoids. Fully ripe fruits (RD0) harvested from the open field (OF)-grown plants showed a shelf-life that exceeded 150 days (Fig. 1b). External fruit colour index following an initial slight decrease remained unchanged at 60 (RD60) and 150 (RD150) days after harvesting (Fig. 2a).
A steady content of carotenoids was previously reported in long shelf-life tomatoes [15]. In LUC the total carotenoids even increased by 28% in RD150 vs RD60 (Table S1). Given that the colour value of LUC fruits is stable in post-harvest, being related to the lycopene content, we can hypothesize that the increase in carotenoids does not depend on the lycopene in LUC. Fruit weight showed a reduction of 18% at RD150, associated with a late fruit withering (Fig. 1b, Fig. 2b). The firmness of the fruits decreased by 31% at RD60 due to the softening process (Fig. 2c). This finding agrees with the evidence that fruit pigmentation and fruit softening are independently regulated processes [16]. Ethylene production increased immediately after breaker, peaked at 4 days after breaker and dropped to a negligible level in post-harvest (Fig. 2d). By way of comparison, we monitored fruits of "Ailsa Craig", which is a short shelf-life cultivar frequently used to study ripening and senescence behaviour [3]. Colour of AC fruits was comparable with that of LUC in the early stages of ripening, whereas significantly lower values characterized AC vs LUC at the later stages (Fig. 2a). Weight and firmness of fruits decreased in AC much faster than LUC (Fig. 2b, c). Furthermore, a greater amount of ethylene was produced by AC fruits in all ripening stages compared with LUC (Fig. 2d). Finally, the deterioration of AC fruits began ∼40 days after harvesting. In fruit samples from OF-grown LUC plants, we evaluated other quality traits such as the content of soluble solids that decreased from 8.9 to 8.1 • Brix during post-harvest (Table S1). The increase in H 2 O 2 concentration suggested that LUC fruits experienced increasing levels of oxidative stress during storage. Generally, a plethora of antioxidant molecules including carotenoids, ascorbic acid, and f lavonoids are synthesized in tomato fruits to prevent oxidative stress. Among hydrophilic antioxidants, ascorbic acid decreased at RD150, while f lavonoids increased significantly (Table  S1). The antioxidant capacity of the hydrophilic fraction proved to be constant throughout the post-harvest phases. By contrast, the antioxidant capacity of the lipophilic components showed a significant decrease during post-harvest. Overall, the total antioxidant capacity did not change significantly after harvesting (Table S1). In general, the characterization of LUC fruits provided findings consistent with its distinctive long shelf-life, including a reduced softening rate, similarly to Mediterranean long shelf-life landraces [4], and a lower amount of ethylene, as previously reported for some "de penjar" accessions [17].

Transcriptome reprogramming during post-harvest
To investigate variation in gene expression profiles of the fruit in post-harvest, RNA-seq was performed at three time points (i.e. RD0, RD60, RD150). Each time point was analysed in triplicate and principal component analysis (PCA) revealed homogeneity between replicates (Fig.   S1). In total, 397 439 858 reads were sequenced, 99.7% of which were successfully aligned back onto the LUC genome (Table S2). A total of 15 844, 14 755 and 15 808 genes were expressed (FPKM ≥1) in fruit pericarp at RD0, RD60 and RD150, respectively. FPKM-normalized values for each transcript were hierarchical-clustered and used to generate the heatmap that clearly shows a distinct expression pattern at each time point (Fig. 3a). To validate RNA-Seq data, RT-qPCR analysis was carried out on eight genes involved in different biological processes, thus showing that expression profiles obtained with the two methods were comparable (Fig. S2). Differentially expressed genes (DEGs) were identified by combining three statistical approaches (false discovery rate, FDR ≤ 0.05; Methods S1). Only DEGs called by all the three methods were considered for downstream analysis ( Fig.  S3 and Table S3 at doi:10.17632/jnhhh4v9zm.1). Overall, 5460 and 7619 DEGs were identified in RD60 vs RD0 and RD150 vs RD60 comparisons, respectively, with 2642 (48%) up-and 2818 down-regulated genes in RD60 vs RD0 and 4060 (53%) up-and 3559 down-regulated genes in RD150 vs RD60 (Fig. 3b). In the latter comparison the total number of DEGs was higher than in RD60 vs RD0, with the amount of up-regulated genes greater than that of downregulated ones. In addition, 403 and 546 were respectively the common up-and down-regulated DEGs between the two comparisons.
Gene ontology enrichment analysis was performed to identify enriched (i.e. over-represented) GO terms associated with DEGs using the whole transcriptome set as a background reference. For each GO category a Fisher's exact test (Bonferroni-corrected, FDR ≤ 0.05) was performed. Looking at up-regulated genes in RD60 vs RD0 the most enriched GO terms (>10%) were "signal transduction", "response to heat", "mRNA splicing, via spliceosome", "photosynthesis, light-harvesting" and "exocytosis" (Fig. S4a). As for up-regulated genes in RD150 vs RD60, the over-represented terms were "carbohydrate transport" (21%) and "retrograde transport, endosome to Golgi" (17%) (Fig. S4b). None of the enriched GO terms among up-regulated genes were shared between the two comparisons. Overall, these results indicated a largescale transcriptional reprogramming, which occurred in the fruit during post-harvest.
We compared expression data generated for LUC with those publicly available for AC to provide additional clues on the genes expressed in post-harvest. The expression data retrieved from the FruitENCODE database refers to different developmental and ripening stages [11]. PCA revealed that LUC samples at the red ripe stage (RD0) and post-harvest stages (i.e. RD60 and RD150) unexpectedly grouped with AC samples in the early ripening stages (Fig. 4).
We calculated the correlation between gene expression profiles and PC1 of PCA, which explained 51.9% of the variance observed in LUC and AC samples. A total of 1084 genes showed a high Pearson's correlation coefficient (≥ 0.75) (Table S4). Among these genes, we recognized ripening-related genes [1] including (i) known regulators such as MADS-RIN, SBP-CNR, NAC-NOR, TAGL1, AP2a; (ii) genes involved in hormone biosynthesis and signaling pathways; and (iii) genes encoding cell wall degrading enzymes. Among ripening regulators, MADS-RIN (Solyc05g012020) showed the highest Pearson's correlation value (0.96). MADS-RIN is a key transcription factor that regulates the ripening process in tomato [18] by targeting several classes of genes including those involved in ethylene biosynthesis and cell wall metabolism [19,20]. Out of the 356 known targets of RIN, 289 showed a positive correlation (average value = 0.72) with PC1 of PCA. We examined the expression of genes involved in ethylene biosynthesis, perception and signaling that had a positive correlation with PC1 of PCA (Table S4). They mostly appeared less expressed in LUC vs AC (Fig. 5a). A similar behaviour was shown by genes encoding components of the core abscisic acid (ABA) signaling pathway [21] (Fig. S5). Given the prolonged firmness of LUC fruits, we examined genes related to cell wall degradation (Table  S4). Genes that were up-regulated at breaker+7 (B + 7) and B + 10 in AC appeared mostly down-regulated at RD0 and RD60 in LUC (Fig. 5b). An increase of expression at RD150 was likely associated with the fruit softening.
We wondered if genes related to the above-reported classes and correlated with PC1 of PCA (Table S4) were regularly expressed during the ripening in LUC fruit. Therefore, we examined by RT-qPCR the expression of representative genes by using pericarps of fruits collected at different ripening stages (MG, B, B + 4, B + 8) from LUC and AC plants grown under the same conditions ( Fig. 6). In AC, MADS-RIN showed a steady increase in expression from MG to B + 8, confirming previous findings [1,22]. Conversely, MADS-RIN expression was significantly lower in LUC and decreased after the peak at  B + 4. ACS2 (Solyc01g095080) that is a target of MADS-RIN and is involved in ethylene biosynthesis [23] showed a significant increase of expression in LUC at B + 4 and dropped-off afterwards.
The expression of ethylene receptors ETR3 (Solyc09g07 5440), corresponding to the Never Ripe locus, and ETR4 (Solyc06g053710), were not significantly different. We examined the ABA receptor SlPYL1 (Solyc08g082180) and the SNF1-related protein kinase SlSnRK2.8 (Solyc04g012160). Although SlPYL1 was expressed to a reduced extent during ripening in both LUC and AC, it showed a significant reduction in LUC vs AC at B + 4. SlSnRK2.8 expression was apparently lower in LUC vs AC but the difference was not statistically significant. Finally, we examined the expression of Lipoxigenase B (SlLoxB, Solyc01g099190) and Pectate Lyase (SlPL, Solyc03g111690) involved, respectively, in membrane breakdown and cell wall degradation and both leading to a shelf-life extension when silenced [24,25]. Both are targets of MADS-RIN and mirrored its expression pattern in LUC vs AC during ripening. Their expression dropped-off downstream the B + 4 peak in LUC, while steadily increasing till B + 8 in AC.
Overall, the expression profiles of ripening-related genes in LUC suggested that ripening and senescence processes were impaired in this landrace. Given the importance of epigenetics modifications in tomato fruit ripening, we investigated whether different miRNAs as well as epi-marks could play a role in determining the long shelf-life phenotype.

Identification of miRNAs and their targets
Single-end sequencing of small non-coding RNAs was performed in fruit samples at three-time points (i.e. RD0, RD60, RD150). About 99% of the reads were successfully aligned back onto the reference LUC genome. Overall, 4582 unique miRNAs were identified and ∼1% was identical to known miRNAs recorded in miRBase Release 22.1 (Table S5).
MicroRNAs were annotated and their anticorrelated targets were identified (Tables S6-S8). The function of target genes was described with GO terms belonging to the "biological processes" domain ( Fig. S6). Enrichment analysis highlighted over-representation of the GO categories "lipid biosynthesis" and "isoprenoid metabolism" in RD150 vs RD60. Networks were built to depict the relationships among DE miRNAs and their anticorrelated targets within DEGs in RD60 vs RD0 and RD150 vs RD60 (Figs. S7, S8). As expected, we found several ripening-related genes among the predicted targets of miRNAs. For instance, miRNA_3127 was up-regulated at RD60 vs RD0. It targets the transcription factor SlAP2a (Solyc03g044300), a member of the APETALA2/Ethylene Response Factor (AP2/ERF) family, which was downregulated at RD60 vs RD0. Expression of SlAP2a was found to be regulated also by miRNA_3648, a member of miR172 family, which is known to regulate genes belonging to the AP2 family [14]. Finally, miRNA_2340 was down-regulated at RD150 vs RD60, while its target the ethylene-insensitive 3-like gene (EIL5, Solyc01g014480), likely involved in fruit senescence, was consistently up-regulated at the same time point.

Genome-wide DNA methylation patterns
A single-base resolution map of DNA methylation across RD0, RD60, and RD150 was generated for LUC fruits to investigate DNA methylation dynamics during postharvest. For each sample at least 88 M paired-end reads were produced. About 83% of the reads were aligned back onto the reference LUC genome, covering 87.4% of the genome. On average, sequenced methylomes had ∼19fold coverage per DNA strand. We compared the DNA methylomes in RD60 vs RD0 and RD150 vs RD60 to detect the differences in methylated Cs and identify the differentially methylated regions (DMRs). During post-harvest, the amount of methylated Cs increased from 1 805 000 in RD60 vs RD0 to 2 574 000 in RD150 vs RD60, while the demethylated Cs fell from 1 342 000 to 945 000 (Fig. 8a).
In RD60 vs RD0, we identified a total of 599 640 DMRs, 55% of which were hypermethylated (hyper-DMRs) ( Table S9). In RD150 vs RD60, the amount of DMRs raised (673153) as well as the frequency of hyper-DMRs (69%) ( Table S9). Hyper-DMRs occurred predominantly in the CHH context at the frequency of 66% and 72% in RD60 vs RD0 and RD150 vs RD60, respectively (Fig. S9). Taken together, our findings pointed out that LUC fruits underwent a global increase in DNA methylation during post-harvest, especially in CHH context. Furthermore, we analysed DMR composition and distribution along with genes and their f lanking sequences. We found that DMRs were preferentially located in intergenic regions (Fig. 8b). In addition, DMRs associated with protein-coding genes showed different patterns across sequence contexts. For instance, hyper-and hypo-DMRs at CGs gradually increased in 5 and 3 flanking regions. By contrast, DMRs at CHGs and CHHs primarily occurred within the gene bodies, whereas they dramatically decreased in the regions flanking the genes (Fig. 8c). To explore whether the increased DNA methylation observed during post-harvest was attributable to enhanced DNA methyltransferase activity or reduced action of DNA demethylation, we examined the expression levels of DNA methyltransferase and demethylase genes. Among the seven methyltransferases functionally characterized in tomato (SlMET, SlDRM1L, SlDRM1L1, SlDRM5, SlMET1L, SlMET3L, SlMETL), we detected changes in transcriptional activity only for SlDRM1L1 (Solyc10g078190), which increased its expression, and SlMET3L (Solyc12g100330) that decreased (Fig. 8d). Regarding the DNA demethylase genes, SlDML2 (Solyc10g083630), which encodes a fruit-specific demethylase [1], progressively reduced its expression during post-harvest (Fig. 8d). These results suggested that the increase in DNA methylation from RD0 to RD150 could be dependent on the dual activity of demethylation decrease and methylation increase.

Genome-wide patterns of histone modifications
ChIP-Seq was performed to trace global patterns of activating (H3K4me3, H3K9K14ac) and repressive (H3K27me3) histone modifications (HMs) in LUC fruit pericarp at RD0, RD60, and RD150. The sequencing generated ∼24.4 Gb reads (Table S10) that were processed to identify regions of ChIP enrichment (i.e. peaks). For each sample, we counted the amount of identified peaks and the amount of bases associated with the histone marks (Table S11). On average, 6.3%, 4.3% and 6.2% of the genome was associated with H3K4me3, H3K9K14ac and H3K27me3 modifications, respectively. In all cases, HMs closely resembled the expected distribution. Indeed, H3K4me3 and H3K9K14ac showed a peak immediately after the putative transcriptional start site, while H3K27me3 signals covered wider regions (Fig. S10).
The most represented CS was E5 (Fig. S11a). As expected, E6 (H3K27me3-only regions) was strongly associated with repetitive sequences (∼7% of total repeats) (Fig. S11b). By contrast, CSs marked by H3K4me3 and/or H3K9K14ac, namely E1, E2 and E3, were found to be associated with promoters (42%) and gene bodies (51%) (Fig. S11c-d). Generally, these CSs referred to the most expressed genes, whereas E4 and E6, marked by H3K27me3, were associated with genes with low levels of expression (Fig. 9b). Based on CS changes between time points, we estimated that 12% of chromatin is dynamic during post-harvest. One-third of CS changes were related to gain/loss of H3K27me3 mark as shown by transition from E5 to E6 and vice versa (Fig. 9c, d).
To investigate transcriptional regulation in postharvest, differences in HMs between stages were highlighted through the call of differentially HM Regions (DHMRs) around the transcriptional start sites. Following the pairwise comparisons of HM enrichment profiles (i.e. RD60 vs RD0 and RD150 vs RD60) we observed an overall increase in H3K4me3 marked transcriptional start site regions at RD60 and RD150 (Fig. 10). Indeed, the balance between gain and loss sums to ∼1.3 Mb considering both comparisons. On the other hand, H3K9K14ac showed a net loss of 1.7 Mb during post-harvest, with an increase  Among DEGs we tracked the expression of genes encoding histone lysine demethylases with Jumonji Cdomain and we found that three genes, SlJMJ6, SlJMJ7, and SlJMJ12, were up-regulated in RD150 vs RD60. Since it was demonstrated that SlJMJ6 demethylates H3K27me3 [26] our findings were consistent with H3K27me3 loss at RD150.

Mechanisms of gene regulation in post-harvest
During post-harvest of LUC fruits, ∼50% of DEGs were found to be engaged in at least one epi-mark change besides being a miRNA target in ∼5% of cases (Table S12, S13). For DNA methylation in the CG, CHG and CHH context, 1944 and 2687 DEGs were associated with DMRs located at promoter regions in RD60 vs RD0 and RD150 vs RD60, respectively. These findings corroborate the possibility that these genes can be subjected to epi-regulation during post-harvest. Gene Ontology enrichment analysis performed on DEGs associated with DMRs showing hypermethylation (Tables S14-S19 at doi:10.17632/jnhhh4v9zm.1) suggested that the increase of DNA methylation, occurring during postharvest, could contribute to long shelf-life by regulating different pathways (Tables S20-S25). Based on the known components of tomato ripening regulatory circuit [1], we found that expression of transcription factors such as TAGL1 (Solyc07g055920) and SBP-CNR (Solyc02g077920), which were down-regulated in post-harvest, was consistent with the increase of DNA methylation at promoter regions.
As far as histone marks H3K4me3, H3K9K14ac, and H3K27me3 are concerned, we found that 986 and 1509 DEGs were associated with DHMRs in RD60 vs RD0 and RD150 vs RD60, respectively (Table S12, S13). Analysis of over-represented GO categories (FDR ≤ 0.05) indicated that DEGs associated with the H3K9K14ac changes (Table  S26) were mostly related to nucleic acid metabolism (RD60 vs RD0) (Table S27). In case of DEGs engaged in H3K4me3 changes (Table S28), GO terms dealing with cell membranes appeared to be enriched in the aging tissues (RD150 vs RD60) (Table S27). Among the subset of DEGs positively correlated with the loss of the activation mark H3K4me3, DNA demethylase SlDML2 (Solyc10g083630) showed a high r value (r = 1) (Table  S27). SlDML2 expression decreased during post-harvest consistently with the general trend of DNA methylation increase. Our findings suggested that SlDML2 could be regulated by H3K4me3 loss in post-harvest. As for H3K27me3 DHMRs, we focused on anticorrelated DEGs and particularly on those associated with the loss of the repressive mark H3K27me3. Indeed, post-harvest is characterized by the peculiar decrease of this histone mark (Table S29). Among DEGs, SlASR3, ABA stress ripening 3 (Solyc04g071590), is negatively correlated (r = −0.94) with H3K27me3. It is reported that ASR3 positively regulates fruit ripening and stress responses [27].
Finally, we identified interesting DEGs by examining the co-occurring epi-marks, including DNA methylation and histone post-translational modifications as well as miRNAs in RD60 vs RD0 and RD150 vs RD60 with respect to the DEGs (Fig. 11a, b). For instance, we found two candidates, i.e. abscisic acid-induced MYB1, SlAIM1 (Solyc12g099120) and Solyc01g096600 which will be investigated for their role in post-harvest.

Discussion
In this study, we investigated the transcriptome and epigenome context of tomato fruit in a landrace (LUC) that is characterized by a long shelf-life. Like most of the Mediterranean long shelf-life landraces, LUC has maintained genetic identity thanks to in situ conservation and it is a valuable genetic resource to improve shelf-life, fruit quality and resilience in tomatoes [4,28].
Multi-omics analysis performed in this study revealed that thousands of genes were differentially expressed at the red ripe stage and during post-harvest (60 and 150 days after harvesting) and about half of DEGs were potentially epi-regulated. Among the regulators necessary for the overall ripening process, MADS-RIN can be the key player in the gene regulatory network leading to delayed fruit ripening in LUC. Previously, it was ruled out that MADS-RIN could be related to long shelf-life in "de penjar" landraces since the rin mutation was absent in those genotypes [5,18]. Recently, genetic polymorphisms were detected in the MADS-RIN region in long shelflife germplasm but their possible relationship with long shelf-life needs to be proven [28]. In LUC, no sequence variation was detected in the coding sequence either in MADS-RIN or in the other ripening regulators except for TAGL1 showing a small deletion along with a SNP [6].
Evidence gathered in this study highlighted that MADS-RIN expression was significantly down-regulated at RD0 in LUC vs AC and progressively declined during post-harvest. About 81% of the genes reported to be direct targets of MADS-RIN [19,20] were down-regulated at RD0, at least. It is known that MADS-RIN is necessary to complete the fruit ripening program as the blocking of its expression leads to the down-regulation of many important ripening-related genes [16,22,29]. For instance, ethylene-related genes including ACC synthases, ACS2 and ACS4, along with signaling factors are targets of MADS-RIN [9,19,20,30], which is also involved in a positive feedback loop leading to ethylene production [11]. Recently, it has been reported that MADS-RIN is required for autocatalytic system II of ethylene biosynthesis [23]. Our study revealed that a large number of genes involved in ethylene biosynthesis, perception and signaling, including ACC synthases (ACSs) and oxidases (ACOs), receptors (ETRs) and response factors (CTRs and ERFs), were significantly down-regulated in LUC at RD0 consistently with the down-regulation of MADS-RIN at the same stage. Ethylene-related genes maintained their repressed status during post-harvest with the exception of a subset of genes showing up-regulation at RD150 likely associated with fruit senescence. Down-regulation of MADS-RIN and ethylene-related genes could explain the downstream effects on genes encoding cell wall degradation enzymes in LUC. Indeed, down-regulation of one or more family members of each type of cell wall modifying enzymes, including polygalacturonases and pectate lyases, occurred at RD0 until RD60. We hypothesize it is likely associated with delayed fruit softening. MADS-RIN and ethylene are required for the expression of genes involved in cell wall degradation, thereby suggesting that both induce normal softening during ripening [23,31]. MADS-RIN has contrasting effects on tomato fruit softening as evidenced by the different allelic mutants [16,23]. Delayed fruit softening was observed in rin mutants with truncated/chimeric RIN protein [16,18,22]. Taken together, our findings indicate that long shelf-life in LUC is associated with the incomplete fruit ripening program that has beneficial effects, including carotenoid accumulation and delayed fruit softening.
Long shelf-life is likely inf luenced by the balance of other hormones in addition to ethylene. ABA-related genes were reported to be involved in tomato fruit ripening [32]. A higher level of ABA observed at the onset of fruit senescence in some "de penjar" landraces was related to higher retention of water in fruits [17]. Since LUC fruits exhibit a reduced water loss, we monitored the expression of genes involved in ABA biosynthesis and signaling pathway, including PYL receptors, PP2C phosphatases, and SNF1-related protein kinases (SnRK2s). Up-regulation of genes encoding key enzymes within the ABA biosynthetic pathway, such as ZEP and NCED, combined with suppression of genes for uridine diphosphate glucosyltransferases (UGTs) involved in ABA homeostasis regulation [33] was observed, suggesting enhancement of ABA level in postharvest. This is consistent with the view that ABA acts during late ripening as a negative regulator of both ethylene biosynthesis and fruit ripening, thus playing a role in shelf-life extension [34]. The variegated expression pattern of PYLs, PP2Cs and SnRK2s during post-harvest is difficult to interpret since their function in ripening is not fully understood, yet. It was reported that the overexpression and suppression of SlPYL9 resulted in opposite ripening phenotypes, accelerated and delayed, via altered expression of ABA signaling genes and that SlSnRK2.8 was differentially expressed in PYL mutants [35].
Given that epigenome reprogramming controls tomato fruit ripening with effects on long shelf-life [1], genomewide analysis of miRNAs, DNA methylation, and histone post-translational modifications (HMs) was performed in this study. We found that about half of DEGs were associated with differentially methylated regions and/or enrichment of HMs and/or potentially targeted by differentially expressed miRNAs during the post-harvest phases. It is known that miRNAs together with other noncoding RNAs participate in the regulation of fruit ripening and senescence [12]. They are involved in different pathways, including those related to ethylene synthesis and fruit softening [36,37]. About 5% of DEGs were estimated to be potential targets of miRNAs in LUC during post-harvest. Network analysis performed in this study shows how miRNAs and in silico identified target genes within DEGs were involved in different pathways underlying long shelf-life.
DNA methylation increased in LUC as the fruit progressed from RD0 to RD150. This is probably due to the interplay between different methylases and the SlDML2 demethylase down-regulation. Since expression of SlDML2 was regularly activated during ripening in both LUC and AC (data not shown), DNA hypomethylation mainly due to increased expression of the SlDML2 could have occurred during ripening in LUC as previously reported in AC [9]. Extension of shelf-life described in SlDML2 knock-down and knock-out mutants occur via DNA hypermethylation and repression/activation of different classes of ripening-related genes [38,39]. Specifically, it has been suggested that a regulatory loop occurs between DNA methylation and ripening transcription factors, including MADS-RIN [38]. Although MADS-RIN was not associated with any DMR in our experiment, we cannot exclude that DNA hypermethylation inf luenced MADS-RIN down-regulation no later than RD0.
Histone modifications also play an important role during post-harvest. Our results provided evidence that the repressive histone mark H3K27me3 is lost at RD150. The biological significance of this loss could lie in a wide activation of genes involved in the fruit senescence program. H3K27me3 is required for regulating the core ripening genes with the aim of keeping the autocatalytic ripening loop under strict developmental control [11]. Previously, it has been found that MULTICOPY SUP-PRESSOR OF IRA1 (SlMSI1), encoding a Polycomb-Group (PcG) protein that acts in Arabidopsis through the Polycomb Repressive Complex2 [40] driving trimethylation of H3K27, inhibited fruit ripening by negatively regulating a large set of ripening-related genes when overexpressed in tomato [41]. A functional study on SlJMJ6 encoding a histone lysine demethylase with Jumonji Cdomain capable of demethylating H3K27me3, demonstrated that the removal of the repressive H3K27me3 mark from ripening-related genes including MADS-RIN and SlDML2 is a key feature of tomato fruit ripening [26]. In LUC, the expression of SlJMJ6 along with SlJMJ7 and SlJMJ12 being all up-regulated at RD150 is consistently associated with H3K27me3 loss. It is worth mentioning here that both SlJMJ6 and SlJMJ12 have mutations in the LUC genome but their relationship with long shelf-life has not yet been proven [6].
Based on our integrative approach including DEGs and co-occurring epi-marks, we identified two interesting candidates. Abscisic acid-induced MYB1, SlAIM1 (Solyc12g099120), was not expressed at RD0, but it gradually increased its expression to ∼700-fold across post-harvest, acquiring the activation histone marks H3K4me3 and H3K9K14ac in combination with DNA methylation increase in promoter and decrease in gene body. To the best of our knowledge, the expression of this gene was never reported in tomato fruit. However, the role of SlAIM1 as an ionic/oxidative stress tolerance factor has clearly been established [42]. Another interesting gene, Solyc01g096600, increased its expression in post-harvest. DNA methylation loss occurred at CG in the promoter, while H3K4me3 loss affected its gene body. It encodes a component (subunit 1) of the THO/TREX (TRanscription-EXport) complex, which is functionally conserved among eukaryotes and contributes to different mRNA biogenesis steps, including splicing and export to the cytoplasm [43]. The complex selectively carries out a regulation of gene expression through the transcription-export machinery. Arabidopsis homolog, THO1/HPR1, showed transcription elongation selectivity for RTE1 (REVERSION TO ETHYLENE SENSITIVITY1) expression, thereby affecting the consequent ethylene-induced leaf senescence [44,45].
On the basis of the main results gathered from our work, two different scenarios associated with early and late post-harvest can be hypothesized in "Piennolo del Vesuvio" landrace "Lucariello" (Fig. 12).
Firstly, the down-regulation of MADS-RIN and ethylenerelated genes has a repression effect on the genes involved in cell wall degradation, thereby delaying fruit softening. Afterwards, up-regulation of ethyleneand ABA-related genes in association with epigenetic changes lead to the beginning of fruit senescence program. In summary, we present an integrative multiomics approach, hitherto unperformed in tomato postharvest, providing new molecular insights into the transcriptome and epigenome landscape underlying long shelf-life. Currently, additional studies are ongoing on the candidates identified in this work, to eventually elucidate their function. Moreover, a segregating population is being developed to find associations between the genetic variants previously identified in LUC [6] and the expression and/or the epigenetic variations observed in this study. We believe these efforts will contribute further to the understanding of the long shelf-life processes. Finally, the exploitation of a valuable genetic resource along with our publicly available datasets will facilitate future works for tomato breeding.

Plant material
Plants of LUC were grown in the open field (OF) on the slopes of the Vesuvius volcano (Campania, Italy) by using a randomized design with 9 blocks and 3 replications. Each block consisted of 20 plants. Fruits were harvested from 2 nd truss at the red ripe stage corresponding to 8 days post-breaker (RD0) and stored in a nonrefrigerated warehouse to be collected at 60 (RD60) and 150 (RD150) days after harvesting. Pericarp from handdissected fruits was frozen in liquid nitrogen, ground by a blender (Waring) and stored at −80 • C until processed. For Chromatin Immunoprecipitation Sequencing (ChIP-Seq), pericarp was pre-processed before liquid nitrogen treatment.
LUC and "Ailsa Craig" (AC) were also grown in a controlled environment at 25 • C with a daily light of 16 hours to collect fruits for comparative analysis. Fruit samples were harvested at the different ripening stages according to Shinozaki et al. [7].

Fruit physico-chemical analysis
Fruit colour was measured using a Minolta Chromameter CR-400 (Konica Minolta Inc., Osaka, Japan). Colour was assessed at two equatorial regions and two regions at the apex and base of each fruit (5 fruits per stage). Readings were taken based on CIE L * a * b * colour system with L * (whiteness or darkness), a * (redness/greenness) and b * (blueness/yellowness). Colour index (CI) was calculated from the equation CI = (2000·a * )/[L * · (a * 2 + b * 2 ) 1/2 ] [46]. Fruit firmness was estimated using a durometer AGROSTA® 100X (Agro-Technologie, Compainville, France) equipped with a 5-mm f lat-head cylindrical probe. The probe was pressed 2.5 mm to the fruit surface and the peak force (N) was recorded.
The content of total soluble solids ( • Brix) from fruit juice (30 fruits) was measured with a digital refractometer equipped with a temperature compensation system (Atago ® Model 3405, Tokyo, Japan). Additional chemical analyses are detailed in Methods S1.
Ethylene was measured on fruits at different ripening and post-harvest stages (5 fruits per stage) based on the procedure described by Nguyen et al. [47] with minor modifications. Each fruit was placed in a sealed gas-tight 100-mL jar for 17 hours. One ml of sample from each headspace was analysed by a gas chromatograph GC Clarus ® 580 (PerkinElmer, Waltham, MA, USA) equipped with a hydrogen flame detector (FID) and a TG-IBOND Alumina column (30 m × 0.53 mm × 10 μm, Thermo Scientific, Waltham, MA, USA) to quantify the ethylene released.
All the above-mentioned analyses were performed with fruits collected from plants of LUC and AC grown at the same conditions.

RNA extraction and sequencing library preparation
Total RNA of nine fruit samples (3 stages x 3 biological replicates) from open-field-grown LUC was extracted and purified according to Di Matteo et al. [48]. RNA amount and quality were assessed using a NanoDrop ND-1000 Spectrophotometer and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Indexed libraries were prepared by Genomix4Life Ltd (Baronissi, Salerno, Italy) from 1 μg/ea of purified RNA with the TruSeq Stranded mRNA and the TruSeq Small-RNA Sample Prep Kit (Illumina, San Diego, CA) according