Histone posttranslational modifications rather than DNA methylation underlie gene reprogramming in pollination-dependent and pollination-independent fruit set in tomato.

Fruit formation comprises a series of developmental transitions among which the fruit set process is essential in determining crop yield. Yet, our understanding of the epigenetic landscape remodelling associated with the flower-to-fruit transition remains poor. We investigated the epigenetic and transcriptomic reprogramming underlying pollination-dependent and auxin-induced flower-to-fruit transitions in the tomato (Solanum lycopersicum) using combined genomewide transcriptomic profiling, global ChIP-sequencing and whole genomic DNA bisulfite sequencing (WGBS). Variation in the expression of the overwhelming majority of genes was associated with change in histone mark distribution, whereas changes in DNA methylation concerned a minor fraction of differentially expressed genes. Reprogramming of genes involved in processes instrumental to fruit set correlated with their H3K9ac or H3K4me3 marking status but not with changes in cytosine methylation, indicating that histone posttranslational modifications rather than DNA methylation are associated with the remodelling of the epigenetic landscape underpinning the flower-to-fruit transition. Given the prominent role previously assigned to DNA methylation in reprogramming key genes of the transition to ripening, the outcome of the present study supports the idea that the two main developmental transitions in fleshy fruit and the underlying transcriptomic reprogramming are associated with different modes of epigenetic regulations. ©2020 The Authors New Phytologist ©2020 New Phytologist Trust.

Introduction

The transition from flower to fruit, the so‐called fruit set, is an important determinant of crop yield and, in the face of global warming and fast‐growing world population, maintaining yield stability is becoming a major challenge. Fleshy fruit development is a genetically programmed process comprising a series of developmental transitions coordinated by a complex network of signalling pathways that trigger massive physiological, metabolic and structural changes (Pandolfini et al., 2007). Fruit set proceeds normally upon successful completion of pollination and fertilization of the mature flower. So far, the molecular nature of the signalling networks that trigger the series of subordinated programmes for fruit set remains poorly understood. Auxin and gibberellin (GA) are two central hormones involved in the flower‐to‐fruit transition, and application of both hormones to unpollinated ovaries can stimulate parthenocarpic fruit growth in several species (Gustafson, 1936; Bünger‐Kibler & Bangerth, 1982; Wang et al., 2005; Martí et al., 2007; Serrani et al., 2008; de Jong et al., 2009; Mounet et al., 2012; Garcia‐Hurtado et al., 2012). Extensive transcriptomic profiling has provided an insight into the genetic regulators underlying fruit set in tomato (Vriezen et al., 2007; Molesini et al., 2009; Wang et al., 2009b; Ruiu et al., 2015; Tang et al., 2015), however the main drivers of the transcriptomic reprogramming underlying this developmental transition remain largely undefined.

Epigenetic remodelling is a major mechanism driving transcriptomic reprogramming associated with plant developmental processes (Bouyer et al., 2011; Baulcombe & Dean, 2014; Xiao & Wagner, 2015; Lang et al., 2017; Narsai et al., 2017; Lee & Seo, 2018; Ji et al., 2019; Borg et al., 2020). In particular, histone modifications involving Polycomb group (PcG) and H3K27me3 have been shown to play key roles in developmental transitions such as the shift from seed to vegetative growth and from vegetative to reproductive phases (Pu & Sung, 2015). Histone posttranslational modifications (HPTMs) and DNA methylation at cytosine residues on the fifth carbon are the main operating modes for epigenetic regulations in addition to small RNAs and histone variants (Henderson & Jacobsen, 2007; Lauria & Rossi, 2011; Chua & Gray, 2018; Reyes et al., 2018; Zhang et al., 2018). DNA methylation has been found in all plants analysed so far (Niederhuth et al., 2016). Methylation in the CG, CHG and CHH (where H = A, T or C) sequence contexts is maintained by Methyltransferase1 (MET1) and chromomethylases (CMT3 and CMT2), respectively (Takuno et al., 2016; Niederhuth et al., 2016), whereas RNA‐directed DNA methylation pathway (RdDM) is involved in CHH methylation maintenance (Stroud et al., 2013; Matzke & Mosher, 2014; Yaari et al., 2019). Several studies have revealed the critical role of DNA demethylation in fruit ripening through de‐repression of key regulator genes such as CNR, RIN and NOR in tomato (Zhong et al., 2013; Liu et al., 2015a; Lang et al., 2017; Cheng et al., 2018; Huang et al., 2019). Also, DNA demethylase DEMETER‐like 2 (SlDML2) has been reported to actively remove methyl groups from methylated DNA during fruit ripening (Liu et al., 2015a), and genomewide DNA methylome studies indicated that hundreds of ripening‐associated genes in sldml2 mutant showed a negative correlation between changes in DNA methylation and changes in gene expression (Lang et al., 2017).

HPTMs are a major guide to the coordinated transcriptomic reprogramming associated with developmental transitions, circadian clock and plant responses to stress (Berr et al., 2011; Malapeira et al., 2012; Liu et al., 2014; Engelhorn et al., 2017; Gu et al., 2017; You et al., 2017; Lee et al., 2019; Song et al., 2019). Specific combinations of HPTMs have been correlated with the active or repressive state of chromatin in terms of gene transcription activity (Li et al., 2008; Wang et al., 2009a; Roudier et al., 2011). HPTMs have been also used to identify central regulators for leaf senescence in Arabidopsis (Ay et al., 2009; Brusslan et al., 2015) and for lipid metabolism in microalgae (Ngan et al., 2015). A widely accepted paradigm postulates that histone acetylation is associated with gene activation, whereas histone methylation can be associated either with activation or repression depending on the lysine residue and the number of added methyl groups. Despite the increasing number of studies addressing the importance of histone modifications in plants, none has been dedicated to its potential involvement in fruit set and subsequent early development. Overall, 12 histone marks have been identified in plants by chromatin immunoprecipitation combined with microarray analysis (ChIP‐chip) or deep sequencing (ChIP‐seq) (Zhang et al., 2009; Roudier et al., 2011; Sequeira‐Mendes et al., 2014). H3K9ac and H3K4me3 are two well studied histone marks that correlate with gene activation in diverse plant development processes including de‐etiolation, leaf senescence, circadian clock, shoot meristem, UV‐B treatment and abiotic stress in Arabidopsis, rice or maize (Casati et al., 2008; Charron et al., 2009; van Dijk et al., 2010; Hu et al., 2012; Malapeira et al., 2012; Zong et al., 2013; Schenke et al., 2014; Brusslan et al., 2015). By contrast, H3K27me3 histone marks are often associated with transcriptionally silenced genes (Li et al., 2007). Methylation of H3K27 is mediated by enhancer of zeste (E(Z)), the catalytic unit of Polycomb Repressive Complex 2 (PRC2), an evolutionally conserved component in living organisms (Schwartz & Pirrotta, 2007). Knockdown of tomato SlEZ2 can significantly affect reproductive development such as flower morphogenesis and fruit set efficiency (Boureau et al., 2016), and repression of FIE, another component involved in H3K27 trimethylation in tomato, results in parthenocarpic fruit formation (Liu et al., 2012). Genomewide studies have unveiled several targets of H3K9ac, H3K4me3 and H3K27me3 in various biological processes in Arabidopsis, rice and maize, yet their putative involvement in fleshy fruit development remains restricted to the role of H3K27me3 in fruit ripening (Lü et al., 2018; Li et al., 2020).

To date, the correlation between histone modifications and reprogramming of gene transcription has been investigated in many plant species, such as Arabidopsis (Charron et al., 2009; Lafos et al., 2011; Brusslan et al., 2015; You et al., 2017; Lee et al., 2019), maize (Rossi et al., 2007; Casati et al., 2008), rice (Li et al., 2008; Liu et al., 2015b), poplar (Li et al., 2019) and moss (Widiez et al., 2014). Such information is still lacking for many economically important crop species like tomato, a model system for fleshy fruit research. The flower‐to‐fruit transition represents a major developmental transition suited to investigate the role of epigenetic variation in fruit set and to address the comparative contribution of DNA methylation and HPTM to the global changes in the level of gene transcription. Whether changes in DNA methylation have a similar critical role in fruit set as in fruit ripening, remains an open question. Our study expanded the current view of the epigenetic regulation underlying transcriptional reprogramming in tomato and provides new insight into the mechanisms controlling a biological process that has a decisive impact on yield and quality of a major crop.

Materials and Methods

Plant materials and sampling

Tomato plants Solanum lycopersicum L. cv Micro‐Tom. were grown under standard culture chamber conditions (14 h : 10 h, 25°C : 20°C, day : night cycle, 80% relative humidity, 250 mol m−2 s−1 intense light).

Ovary samples were collected at anthesis stage and were referred to as 0 d postanthesis (0DPA) when stamens were still loosely enclosed by petals. Fruits at 4 d postanthesis (4DPA) corresponded to the whole developing fruit organ (Xiao et al., 2009; Pattison et al., 2015). Fruit samples collected at 4 d after IAA treatment (4IAA), corresponded to flowers emasculated 1 d before anthesis (to avoid accidental self‐pollination), then the ovaries were treated over the next 4 d with 10 μl 500 μM indole‐3‐acetic acid. Each biological replicate corresponded to a pool of at least 50 ovaries or young fruits from 25 plants.

Chromatin immunoprecipitation and sequencing

ChIP assays were performed as described previously (Gendrel et al., 2005) with minor modifications. Tissues corresponding to 0DPA, 4DPA and 4IAA were cross‐linked by vacuum infiltration for 15 min in 1% formaldehyde. To ensure efficient crosslinking, 4DPA and 4IAA fruits were cut in half before crosslinking. Crosslinking was stopped by adding glycine (0.125 M final concentration) and incubation under vacuum infiltration for an additional 5 min. After washing twice with cold 1× PBS solution, samples were thoroughly dried between paper towels, snap‐frozen in liquid nitrogen and stored at −80°C. To perform the ChIP assay, c. 1 g of cross‐linked tissue was ground to a fine powder in liquid nitrogen and then the chromatin was sheared by Diagenode Bioruptor sonication (5 runs of 10 cycles for 30 s ‘ON’ and 30 s ‘OFF’). The size of the sonicated chromatin was within the range of 100–500 bp. Next, 10 μl of sonicated supernatant was kept aside as input. For each sample (120 μl supernatant), dilution buffer was added to bring the final volume to 1.5 ml and, depending on the histone mark, either 5 μl of H3K9ac or H3K4me3 rabbit polyclonal antibody (Millipore, Darmstadt, Germany; Cat. no. 07‐352; Lot no. 2586454; Cat. no. 07‐473; Lot no. 2430389) or 8 μl of H3K27me3 rabbit polyclonal antibody (Millipore; Cat. no. 07‐449; Lot no. 2475696) was added before incubation overnight (4°C, 10 rpm). For the control experiment without histone mark antibodies, 5 μl of nonimmunised rabbit IgG antibody (Millipore; Cat. no. 12‐370; Lot no. 2426484) were added. For the empty control (Mock), no antibody was added. Afterward, 50 μl of protein A/G agarose beads (Pierce™ Protein A/G UltraLink™ Resin; Thermo Scientific, Waltham, MA, USA; Cat. no. 53133) were added and the samples incubated for 3 h at 4°C. Elution was carried out as described previously (Gendrel et al., 2005). For each sample 10 ng immunoprecipitated DNA was used for library construction and sequencing. Two biological replicates corresponding to independent plant materials were used for ChIP assays, one replicate for genomewide sequencing and the second one for ChIP‐qPCR experiments.

ChIP‐sequencing was performed at the GeT‐PlaGe core facility (INRA Toulouse). Sequencing libraries were prepared using the TruSeq ChIP Library Preparation Kit for Illumina Sequencing. Sequencing was performed on an Illumina HiSeq 3000 system with the Illumina SBS HiSeq 3000 reagent kits (2 × 150 nt paired‐end reads). The enrichment of DNA fragments (% input) was validated by quantitative real‐time PCR using primers listed in Supporting Information Table S1.

Genome bisulfite treatment and sequencing

Genomic DNA was isolated from 0DPA ovaries, 4DPA and 4IAA fruits using the Wizard® Genomic DNA Purification Kit (Promega, Madison, WI, USA). Here, 10 μl of total genomic DNA were used for Whole‐Genome Bisulfite Sequencing (WGBS). WGBS was performed at the Beijing Genomics Institute (BGI, Hong Kong, China). WGBS libraries were prepared according to Bioo Scientific’s protocol using the Bioo Scientific NEXTflex Bisulfite Library Prep Kit for Illumina Sequencing. The bisulphite treatment was performed using the EZ DNA Methylation‐Gold Kit (Zymo Research, Irvine, CA, USA) followed by PCR. Library quality was assessed using an Advanced Analytical Fragment Analyser and libraries were quantified by qPCR using the Kapa Library Quantification Kit. Two biological replicates were performed.

Sequencing was performed on a Illumina HiSeq 3000 system using Illumina SBS HiSeq 3000 reagent kits (2 × 150 nt paired‐end reads; see Table S2 for the number of reads).

RNA sample preparation and sequencing

For each sample, total RNA was extracted from c. 200 mg tissue, using the TRIzol RNA Isolation Kit (Thermo Fisher Scientific). After DNA removal (DNA‐free™ DNA Removal Kit, Ambion, Austin, TX, USA), RNA was purified and the quality checked using an Agilent 2100 bioanalyser. Only samples with a RIN > 8.6 were used for Illumina sequencing. Three biological replicates were performed for each sampling stage. Paired‐end RNA‐sequencing (2 × 125 nt) was carried out at the GeT‐PlaGe core facility (INRA Toulouse) using a TruSeq Illumina SBS Kit V4 and a HiSeq 2500 platform.

Gene ontology analysis

Gene ontology (GO) analysis of selected genes was performed using the R package goseq (Young et al., 2010). All GO terms used in our study were obtained from the tomato genome website (https://solgenomics.net/). Gene length bias existing in RNA‐seq was taken into account when the enrichment of the GO category was computed. An over‐represented P‐value of < 0.05 was used to select significantly enriched GO categories.

Identification of putative orthologues in the tomato genome

Genes were functionally categorised based on orthologues from the well studied model plant Arabidopsis, with manual re‐assignment according to the tomato genome (Tomato Genome Consortium, 2012) and NCBI annotations (Table S4). Local Blastp analysis was performed to obtain putative orthologues by blasting with the Arabidopsis protein database (TAIR10) with an E‐value < 1E−20 and maximum selection of three targets. After removing the redundancy, the putative orthologues with the highest score was selected for further analysis.

Generation of CRISPR/Cas9 mutants and genotyping

Constructs for CRISPR/Cas9 mutagenesis was performed following standard protocols (Belhaj et al., 2013; Brooks et al., 2014). Briefly, two single‐guide (sg)RNAs were designed ahead of the coding sequence for the SET domain of the target genes, using the CRISPR‐P server (http://cbi.hzau.edu.cn/cgi‐bin/CRISPR) (Lei et al., 2014). All vectors were assembled using the GoldenGate cloning system (Werner et al., 2012). sgRNA1 and sgRNA2 were first cloned into the Level 1 vectors pICH47751 and pICH47761 driven by the Arabidopsis U6 promoter, respectively. Level 1 constructs pICH47732‐NOSpro::NPTII, pICH47742‐35S:Cas9, pICH47751‐AtU6pro:sgRNA1 and pICH47761‐AtU6::sgRNA2 were then assembled into the Level 2 destination vector pAGM4723. All sgRNA sequences are listed in Table S1. For genotyping of the transgenic lines, genomic DNA was extracted using the ReliaPrep™ gDNA Tissue Miniprep System (Promega). CRISPR/Cas9‐positive lines were further genotyped for mutations using the primers listed in Table S1.

Results

Transcriptomic profiling of the flower‐to‐fruit transition

To investigate the extent to which epigenetic modifications were associated with the transcriptomic reprogramming underlying the flower‐to‐fruit transition, the same plant material samples were used for WGBS, chromatin immuno‐precipitation coupled to deep sequencing (ChIP‐seq) and for genomewide transcriptomic profiling by RNA‐seq. To test first whether pollination‐dependent and ‐independent fruits set shared similar genetic reprogramming, ovaries of emasculated flowers at the anthesis stage (0 day postanthesis, 0DPA) and young developing fruit initiated by pollination (4DPA) or by auxin (4IAA), were collected and subjected to RNA‐seq profiling (Fig. 1a). Auxin‐treated ovaries underwent active growth, reaching the same size as pollination‐induced fruit at 4DPA and 9DPA, by contrast with the control ovaries treated with mock solution that failed to grow (Fig. 1a). Deep sequencing generated 23–33 million reads depending on the sample, with 79–85% of the reads being uniquely mapped to the S. lycopersicum genome (ITAG2.3; Table S2). Gene count analysis showed that biological replicates from each stage/condition clustered together, with 4DPA and 4IAA samples displaying a high degree of similarity while being clearly distinct from 0DPA tissues (Fig. 1b). In total, transcripts corresponding to 28 466 genes, representing 82% of the 34 727 annotated tomato genes, were detected in at least one of the three samples. Using a threshold fold change ≥ 2 and adjusted P‐value < 0.01, allowed the assignment of 7582 genes (3845 upregulated and 3737 downregulated) as differentially expressed genes (DEGs) in pollination‐induced fruit set and 5447 genes (2962 upregulated and 2485 downregulated) in auxin‐induced fruit (Table S5). Notably, 4219 DEGs were shared between auxin and pollination‐induced fruit (Fig. 1c) suggesting that auxin triggered fruit set through a genetic reprogramming largely similar to that induced by pollination/fertilisation.

Fig. 1

Genomewide transcriptomic profiling of tomato genes during fruit set process by natural pollination and auxin‐inducing. (a) Fruit initiation programme in tomato (cv Micro‐Tom). Unpollinated ovaries at 0 d postanthesis (0DPA) and young developing fruits at 4DPA and at 4IAA were sampled for RNA‐sequencing. Bars, 1 cm. (b) Cluster dendrogram of gene counts in RNA‐sequencing samples. Distance matrix of gene counts from all RNA‐seq libraries were implemented by DESeq2. Dendrograms were generated by hierarchically clustering samples based on distance values. The darker blue indicates a closer distance. (c) The number of differentially expressed genes during fruit set. Fold change > 2 and P‐value < 0.01. (d) GO enrichment analysis for DEGs regulated by pollination or auxin during fruit set. Selected significant enriched biological processes (Benjamini–Hochberg (BH) adjusted over‐represented P‐value < 0.05) were annotated in the figure. Genes with P‐value < 0.05 were selected for GO analysis. (e) Differential expression of cell cycle related genes during fruit set. Log2 of expression were indicated as gene expression in y‐axis, and x‐axis refers to the gene name annotated either based on tomato referenced studies or, when missing, according to the best corresponding orthologue in Arabidopsis. Gene names were ordered by the group including common DEGs (blue shaded), auxin‐specific DEGs (yellow shaded) and pollination‐specific DEGs (purple shaded). Solid line in the graph represents natural pollination while dash line represents auxin treatment. Genes with significant differential expression are marked by asterisks (**, fold > 2 and 0.001 < P‐value < 0.01; ***, fold > 2 and P‐value < 0.001). (f) Differential expression of histone modification related genes during fruit set.

Biological processes enriched during fruit set

GO analysis of the DEGs, performed separately for upregulated and downregulated genes in pollination‐dependent and auxin‐induced groups, indicated that 152 and 128 GO terms were significantly enriched in 4DPA and 4IAA samples, respectively (Fig. 1d; Table S6). Among these, 83 biological processes were commonly enriched in upregulated genes, including those related to cell proliferation and differentiation, photosynthesis and hormone regulation. Out of the 96 cell division‐related genes identified in the tomato genome, 50 were upregulated, reflecting the active ongoing cell division at early stages of fruit development (Fig. 1e). The data indicated that pollination and auxin‐induced fruit set trigger the same subset of cell division genes, although pollination seemed to recruit a larger set of genes during this developmental transition.

RNA‐seq profiling also revealed that a large number of DEGs related to ‘DNA methylation’, ‘histone modification’, ‘histone lysine methylation’, ‘H3‐K9 methylation’ and ‘histone phosphorylation’, were significantly enriched (Fig. 1d). To capture all epigenetic‐related genes present in the tomato genome we performed a genomewide Blast search to identify new genes related to histone modification, DNA (de)methylation, chromatin remodelling and PcG, based on the phylogenetic relationship with those described in Arabidopsis. The outcome of this search extended the total number of tomato genes putatively involved in epigenetic processes from 111 to 213, of which 137 belonged to the histone modification category including PcG, 42 to DNA methylation and 34 to chromatin remodelling (Table S4). Interestingly, among 20 DEGs corresponding to histone modifiers 12 were significantly upregulated in 4DPA, of which five were also upregulated in 4IAA fruit (Fig. 1f).

Transition from flower to fruit is associated with a decrease in global DNA methylation

Study of methylomes already reported the critical role of DNA methylation in the transition to ripen of tomato fruit (Zhong et al., 2013), however its potential contribution to fruit set remains unexplored. To gain insight into the dynamics of DNA methylation during fruit set induced by either pollination or auxin treatment, ovaries of emasculated flowers (0DPA) and 4DPA or 4IAA young fruit were collected, genomic DNA was extracted and then subjected to WGBS analysis at single‐base resolution. Over 120 M paired‐end reads were produced for each sequencing library, covering around 70% of the S. lycopersicum genome (ITAG2.3) (Table S2). Each methylome was sequenced with an average of 11‐fold coverage resulting in sequencing depth comparable with DNA methylomes previously reported for tomato (Zhong et al., 2013; Lang et al., 2017; Lü et al., 2018).

For each sample (0DPA, 4DPA and 4IAA) two biological replicates were treated separately giving rise to similar data (Fig. S1a) and revealing that on average, 83% of CG cytosine contexts (c. 530 million), 58% of CHG (c. 523 million) and 12% of CHH (c. 2.93 billion) were methylated in 0DPA samples. The three types of methylated cytosine (5mC) displayed a slight decrease at 4DPA, with an average methylation level in CG, CHG and CHH cytosine contexts reaching 80%, 55% and 10%, respectively (Fig. 2a). A smaller, but not significant, decrease was also observed in 4IAA fruits. In all samples, the cytosine methylation levels for all cytosine contexts were at similar magnitudes as those reported in ripening fruit (Zhong et al., 2013). To identify the differentially methylated regions (DMRs), the bins method was used to compare the 0DPA ovary either to 4DPA (4DPAvs0DPA group) or 4IAA fruit (4IAAvs0DPA group). When considering 30%, 20% and 10% differences in DNA methylation level for CG, CHG and CHH contexts, respectively, 140 148 DMRs representing 0.29% of the tomato genome were identified in 4DPAvs0DPA group. Of these, 2% (2831) corresponded to CG‐DMRs, 14.2% (19 835) to CHG‐DMRs and 83.8% (117 482) to CHH‐DMRs (Fig. 2b; Table S7). In total, 88 919 DMRs (0.21% of the tomato genome) were identified in the 4IAAvs0DPA group, among which 2.1% (1889) corresponded to CG‐DMRs, 16.7% (14 853) to CHG‐DMRs and 81.1% (72 177) to CHH‐DMRs (Fig. 2b). Most DMRs in both 4DPA or 4IAA corresponded to hypomethylation (Fig. 2c). Notably, only a small proportion of DMRs overlapped between the two groups (Fig. 2d) suggesting that the observed changes in DNA methylomes might reflect distinct features between pollination‐induced and auxin‐induced fruit set. The lower number of DMRs associated with auxin‐induced fruit might reflect, at least partly, the absence of genetic reprogramming related to seed development in these parthenocarpic fruit. The distribution of DMRs in the tomato genome indicated they were mostly located in transposable element rich regions for all cytosine contexts (Fig. S1b).

Fig. 2

Differentially methylated regions (DMRs) during pollination‐dependent and auxin‐induced tomato fruit set. (a) Global DNA methylation levels in each cytosine context. (b) Total number of DMRs identified during pollination‐induced (#4DPAvs0DPA) and auxin‐induced (#4IAAvs0DPA) fruit set using ‘bins’ method (window size, 100 bp and gap size, 50 bp). (c) Percentage of hyper‐ and hypo‐methylated DMRs. (d) Overlapping of DMRs between pollination‐dependent and auxin‐induced fruit set, in CG, CHG and CHH sequence context. (e) Correlation between DEGs and DMRs during pollination‐dependent (left panel) and auxin‐induced (right panel) fruit set. X‐axis represents log2‐transformed differentially methylated (DM) fold and y‐axis represents log2‐transformed expression change. Spearman correlation coefficients (R) and corresponding P‐value in each case were noted in each panel. Promoter genes containing ≥ 2 DMRs in their promoter regions were removed for correlation analysis. (f) Distribution of DMRs in genic region including 2 kb upstream promoter, exon, intron and 2 kb downstream terminator. (g) Change in expression of promoter DMR genes during pollination‐dependent (upper panel) and auxin‐induced (lower panel) fruit set. Promoter DMR genes correspond to all methylated cytosines in 2 kb promoter regions in all types of sequence context. The number of DE genes (fold change > 2 or < −2 and P‐value < 0.01) is mentioned the red area of the circle for upregulated and green area for downregulated.

DMRs weakly correlated with changes in gene expression during fruit set

Cross‐referencing the RNA‐seq and Methylome data revealed no correlation between DMR and DEGs for both pollination‐induced and auxin‐induced fruit set under all three types of cytosine contexts (Fig. 2e). Although 44% of genic DMRs corresponded to promoter regions (Fig. 2f), the overwhelming majority (71%) of these did not correspond to DEGs. Moreover, about half of upregulated DEGs were not associated with promoter hypomethylation, and half of downregulated DEGs did not display promoter hypermethylation (Fig. 2g). Notably, genes related to hormones GA and auxin, known to be instrumental to the fruit set process, further illustrated that changes in cytosine methylation were not necessarily associated with gene expression reprogramming underlying the flower‐to‐fruit transition (Fig. S2). Indeed, SlGA20ox1, a GA synthesis gene, was upregulated in both pollination and auxin‐induced fruit set while undergoing hypomethylation in one case and hypermethylation in the other case (Fig. S2). The lack of correlation between DEGs and DMRs is also illustrated by SlARF9A, another critical regulator of fruit set, which showed differential regulation but no significant change in DNA methylation at the promoter level (Fig. S2).

RNA‐seq profiling indicated that, with the exception of the low expressing SlDML2, none of the tomato DML‐like genes involved in DNA demethylation showed significant increase in expression during fruit set either triggered by natural pollination or by auxin treatment (Fig. S3). Conversely, a number of genes related to the RdDM pathway (13 out of 23 in total) were downregulated during both pollination‐induced and auxin‐induced fruit set, suggesting that reduced DNA methylation activities via the RdDM pathway may contribute to the global DNA hypomethylation observed during the flower‐to‐fruit transition. By contrast, SlMET1 and SlCMT3‐like1/2, encoding enzymes necessary for maintenance of CG and CHG methylation, displayed increased expression at 4DPA and 4IAA (Fig. S3), likely to be related to their putative role in maintaining DNA methylation following replication during the extensive cell division characterising early stages of fruit development.

Genomewide mapping of H3K9ac, H3K4me3 and H3K27me3 histone marks

Given the high representation within DEGs of histone modifier genes (Fig. 1d) known to mediate posttranslational modifications of histone H3, including acetylation of lysine 9 residues (H3K9ac) and trimethylation of lysine 4 (H3K4me3) and lysine 27 (H3K27me3), we investigated the genomewide distribution of these three histone marks by ChIP‐seq approach. The number of sequencing reads corresponding to immuno‐precipitated chromatin and control input samples (Table S2) ranged from 48 to 80 million, and on average 95% of the reads mapped to the S. lycopersicum reference genome. The enrichment of the three marks was validated for nine selected regions (genes) by targeted real‐time PCR using independent 0DPA and 4DPA samples, giving results that were fully consistent with the ChIP‐seq data (Fig. S4).

The total number of peaks corresponding to H3K9ac or H3K4me3 was higher than that corresponding to H3K27me3 (13 072 to 14 308) in both 0DPA and 4DPA samples (Fig. 3a) and, depending on the sample, these histone mark‐associated regions covered 4.4% to 6.8% (36–56 Mb) of the tomato genome sequence (Table S8). The breadth of the covered regions (Fig. 3b) was higher for H3K27me3 (median length 2500–2900 bp) than for H3K4me3 (1653–1722 bp) and H3K9ac (1748–1763 bp). The three marks were mainly associated with gene‐rich euchromatin regions, even though H3K27me3 displayed higher association with heterochromatin (Fig. S5a). More than 80% of the regions associated with H3K9ac and H3K4me3 mapped to genes, while < 40% of H3K27me3 were located in genic regions (Fig. 3c). Further inspection of histone mark‐associated regions revealed that most peaks (>74% of H3K9ac, >89% of H3K4me3, and >88% of H3K27me3) were limited to single genes, and only a minor fraction encompassed two or more genes (Fig. S5b). In 0DPA tissues, 53% of the annotated tomato genes were associated with H3K9ac, 54% with H3K4me3 and 19% with H3K27me3, and these proportions only slightly increased at 4DPA stage (60%, 56% and 19%, respectively).

Fig. 3

Genomewide identification of histone modified regions during tomato fruit set. (a) Number of identified regions for H3K9ac, H3K4me3 and H3K27me3. (b) Length distribution of the peak regions. The outliers are dotted with black colour. The histone marks are displayed in yellow (H3K9ac), red (H3K4me3) and blue (H3K27me3). (c) Frequencies of peaks associated with genic and nongenic regions. A region spanning 1.5 kb upstream of the annotated transcription start site (TSS) to 0.5 kb downstream of transcription end site (TES) was designated as a genic region. Genic + repeats, the genic regions overlapped with TE. (d) Average association profile of input (grey), H3K9ac (yellow), H3K4me3 (red) and H3K27me3 (blue) in genic regions at 0DPA and 4DPA. The gene set is adapted from publicly available RNA end‐sequencing data (Zhong et al., 2013) which defines the TSS and TES. Mean counts within 100‐bp window covering 2.5 kb upstream to 2.5 kb downstream the TSS and, 2.5 kb upstream to 25 kb downstream the TES were extracted and plotted. (e) Global view of gene expression and histone mark association at 0DPA. 23 852 genes were filtered by association either with H3K9ac, H3K4me3 or H3K27me3. The expression levels were used as anchors to sort genes. The occupancy of histone marks in the gene region spanning from 2 kb upstream to 2 kb downstream of CDS was represented and visualised by DeepTools.

To investigate the association pattern of the three histone marks at the gene level, the average read count corresponding to each mark was examined using a large set of genes (16 982), the annotation of which has been refined by RNA end‐sequencing (Zhong et al., 2013). H3K9ac and H3K4me3 showed a similar high enrichment in a narrow region downstream of the Transcription Start Sites (TSS) by contrast with H3K27me3 marks that showed lower intensity and were located throughout the gene bodies, as well as upstream of TSS and downstream of the transcription end site (TES) regions (Fig. 3d,e).

To determine the putative correlation between histone marks and the expression level of individual genes, we profiled the ChIP signal intensity in genic regions using all genes that were associated with at least one of the three histone marks (23 852 genes in total). The data revealed that H3K9ac and H3K4me3 marks co‐occurred in 95% of the cases in the same set of genes. Moreover, a genomewide view of gene expression levels and histone mark distribution revealed that these two marks strongly correlated with high levels of gene expression and confirmed their high enrichment at the 5′‐end of genes (Fig. 3e). Profiling the three epigenetic marks at 4IAA fruits by ChIP‐seq using the same parameters as indicated for 4DPA fruits, identified 24 938 regions associated with H3K9ac, 21 827 with H3K4me3 and 13 982 with H3K27me3 (Fig. S6a; Table S8). These numbers were similar to those observed for pollination‐induced fruits (4DPA) and so were the occupancy of all histone marks associated regions (Fig. S6b), the number of histone mark‐associated genes (Fig. S6c,d) and the histone mark distributions in genic regions (Fig. S6e). Overall, the association patterns of the three histone marks are consistent with known patterns in plants and other eukaryotic organisms (Zhang et al., 2007, 2009; Wang et al., 2009a; He et al., 2010; Zhou et al., 2010; Lafos et al., 2011; Veluchamy et al., 2015).

Pollination‐dependent and auxin‐induced fruit set share similar histone posttranslational modifications

Plotting the histone mark association and the expression changes revealed a strong correlation between H3K9ac and H3K4me3 marks and DE genes in both pollination‐dependent and auxin‐induced fruit sets, whereas no such correlation was observed with H3K27me3 histone marks (Fig. 4a). Genomewide profiling of the dynamic changes of the three histone marks revealed that 13 451 genes were differentially associated (DA) with H3K9ac, 11 333 with H3K4me3 and 3973 with H3K27me3 (Fig. 4b; Table S9). A high proportion of H3K9ac and H3K4me3 DA genes displayed differential expression in both pollination and auxin‐induced fruit set, while this proportion was lower for H3K27me3 DA genes (Fig. 4b). Assessing the expression level for genes showing significant change in histone marks indicated that enrichment in H3K9ac or H3K4me3 marks corresponded to enhanced gene transcript levels, while a decrease in these marks was associated with reduced gene expression (Fig. 4c,d). The close similarity between the DEGs induced by auxin and by pollination, prompted the investigation into whether the two types of fruit set also shared similar patterns of HPTMs. Markedly, a high proportion of DA genes associated with pollination‐dependent fruit set overlapped with DA in auxin‐induced fruit initiation (Fig. 4e). While these data further emphasised the strong correlation between the change in histone marks and the transcriptomic reprogramming underlying the flower‐to‐fruit transition, they also clearly indicated that auxin‐induced and pollination‐dependent fruit set relied mainly on a common set of genes.

Fig. 4

Correlation between differentially associated marks (DA) and differentially expressed genes (DE) during the flower‐to‐fruit transition in tomato. (a) Correlation between DE and DA during fruit set. X‐axis represents log2‐transformed DA and y‐axis represents log2‐transformed expression change (DE). The spearman correlation coefficients (R) and the corresponding P‐value were mentioned in each case. (b) DA and DE genes in pollination and auxin‐induced fruit set. The number of DA (light green colour bar) and DE (dark green colour bar) genes (P < 0.01) were indicated above the bars. ‘P’ and ‘A’ represent 4DPAvs0DPA and 4IAAvs0dDPA group, respectively. (c) Average expression level of DA genes based on fold change ≥ 2 and P‐value < 0.01. The number of DA genes is indicated at the top of each panel. Wilcox test was used to compare means in 0DPA and 4DPA samples. *, P < 0.05; ***, P < 0.001. (d) Density profile of histone marks (DA fold change > 2 and P‐value < 0.01) and gene expression change represented at the gene level in 0DPA and 4DPA tissues. (e) Genes differentially associated with histone marks (P < 0.01) overlapping between pollination‐dependent and auxin‐induced fruit. Only DA regions corresponding to single genes were considered. Per cent refer to the fraction of genes undergoing histone post‐translational modifications (HPTM) in pollination (P) and auxin (A) induced young fruit.

Gain or loss of histone marks correlates with changes in the expression of important fruit set‐related genes

Interestingly, plotting DAs and DEGs allowed the identification of a set of highly enriched GO terms common to pollination‐dependent and auxin‐induced fruit set (Fig. S7), indicating that they might be essential for the initiation of the fruit development programme. Genes involved in important biological processes underlying the fruit set transition underwent concomitant changes in both histone marks and gene expression (Fig. S7; Tables S9, S10). Of particular interest, a large number of genes related to hormones, like auxin and gibberellin, two hormones known to be critical for fruit initiation, was differentially expressed and DA with histone marks (Table S11). Out of 112 auxin‐related genes identified in the tomato genome, 44 were differentially expressed during either pollination‐induced or auxin‐induced fruit set, of which 32 (73%) underwent consistent differential association with at least one transcriptionally permissive histone mark (Table S11). Genes involved in all aspects of auxin metabolism and responses (Fig. 5a) were affected by these changes including auxin synthesis (tryptophan aminotransferases and flavin monooxygenases), transport (SlPINs, SlLAXs and SlPILSs) and signalling (13 Aux/IAAs and 5 Auxin Response Factors). Gibberellin‐related genes also underwent dramatic changes in their status, with 17 out of 34 displaying differential expression between 0 and 4DPA or 4IAA. Of these 17 DEGs, 14 were DA with at least one histone mark (Tables S8, S10) and, among these, 11 genes displayed consistent enrichment in H3K9ac or H3K4me3 transcriptionally permissive histone marks, including KAO2‐Like2, GA20OX1 and GID1B‐like (Fig. 5b). Strikingly, the overwhelming majority of ethylene‐related DEGs (46 out of 61) were downregulated (Table S11), including those involved in ethylene biosynthesis (four ACC synthases and five ACC oxidases) and response (five ETRs, one GRL, one CTR, five EIN‐like, three EBFs and 23 ERFs). The vast majority of these ethylene‐related DEGs (78.6%) showed decreased association with at least one transcriptionally permissive histone mark as exemplified by ACO4, ETR4 and EIL3 genes (Fig. 5c). Genes related to other hormones also underwent major changes in transcript levels associated with differential histone marks (Table S11), supporting the idea that the fruit set process required intricate multiple hormone signalling, strongly associated with changes in histone marks.

Fig. 5

Differential expression and differential histone mark associations of auxin‐related, gibberellin‐related and ethylene‐related genes during tomato fruit set. Differential mark association and gene expression of (a) auxin synthesis, transport and signalling genes, (b) gibberellin‐related genes and (c) ethylene‐related genes were visualised in Integrative Genomics Viewer (IGV). Histone mark associations are marked blue (upper panel) and gene expression red (lower panel).

A large proportion of DEGs involved in processes essential for initiation of fruit development, such as cell division, hormone signalling and seed development, exhibited concomitant changes in H3K9ac and H3K4me3 transcriptionally permissive histone marks (Table 1; Table S9). Most cell division‐related genes that displayed significant change in their transcript levels also showed consistent enrichment in H3K9ac and H3K4me3 marks in both pollination‐dependent (84%) and auxin‐induced fruit set (86%). Likewise, hormone‐related genes that displayed differential expression also underwent consistent changes in H3K9ac and H3K4me3 histone marks (Table 1). The consistency between differential expression and HPTMs was also observed for genes related to seed and embryo development although, not surprising, the number of DEGs in this category was much higher in pollination‐dependent (456) than in auxin‐induced (282) fruit set, this latter process leading to seedless fruit.

Table 1

Transcriptomic changes and histone marks reprogramming of genes involved in processes known to be essential for fruit set.

Group	DEG a	DA b  + DEG
		Consistent change c	Inconsistent change d
		Nb (%)	Nb (%)
4DPAvs0DPA
Cell division	50	42 (84)	1 (2)
Hormone‐related	220	159 (72)	21 (10)
Embryo/seed e	456	220 (48)	33 (7)
4IAAvs0DPA
Cell division	35	30 (86)	2 (6)
Hormone‐related	155	133 (86)	20 (13)
Embryo/seed	282	108 (38)	31 (11)

DEG differentially expressed genes with P‐value < 0.01, fold change > 2.

DA genes associated with gain/loss of H3K9ac or H3K4me3 histone marks with P‐value < 0.01.

Consistent change corresponds to enhanced expression associated with gain of H3K9ac or H3K4me3 marks, or to decreased expression associated with loss of H3K9ac or H3K4me3 marks.

Inconsistent change corresponds to enhanced expression associated with loss of H3K9ac or H3K4me3 marks, or to decreased expression associated with gain of H3K9ac or H3K4me3 marks.

Cluster of 12 and 19 genes from previous tissue‐specific transcriptomic data (Pattison et al., 2015).

Overall, the outcome of the study revealed that a high fraction of DEGs (79–90%) underwent changes in transcriptionally permissive histone marks, while only a minor fraction (23–28%) displayed changes in cytosine methylation in their promoter region (Table 2). These data clearly support the notion that the transcriptomic reprogramming underlying fruit set is more strongly associated with HPTM than with rearrangement in cytosine methylation.

Table 2

Correlation between the two modes of epigenetic modifications and the transcriptomic changes during fruit set.

Group	Promoter DMR a /DEG b	DA c /DEG
4DPAvs0DPA	28%	79%
4IAAvs0DPA	23%	90%

Percentages for promoter DMR/DEG refer to percentage of genes displaying increased DNA methylation with decreased gene expression level, or decreased DNA methylation level with increased gene expression level. DA/DEG percentages refer to changes in active histone marks with concomitant changes in expression level.

DMR, differentially methylated regions (P < 0.01) in 2 kb promoter regions at all cytosine sequence contexts. Promoters associated with different changes of DMRs were removed before analysis.

DEG, differentially expressed genes with P‐value < 0.01, fold change> 2;

DA, genes differentially associated with gain/loss of H3K9ac or H3K4me3 histone marks with P‐value < 0.01.

Histone modifier genes putatively involved in tomato fruit set

A large proportion (31 out of 39) of the epigenetic‐related genes assigned as DEGs exhibited changes in histone marks (Table S9; Fig. 6a); mining the expression dataset from different tomato cultivars in the TomExpress (http://gbf.toulouse.inra.fr/tomexpress) database identified six differentially expressed histone modifiers that displayed an expression pattern associated with the fruit set transition (Fig. 6b). Further analysis of their expression levels in various tomato tissues and development stages revealed that SlSDG5, SlSDG16, SlSDG27, SlSDG30 and SlPRMT8 genes encoding putative histone methyltransferases were upregulated at the transition phase, whereas SlJMJ17 histone demethylase gene displayed a decreased expression (Fig. S8). Because SlSDG27, SlSDG5 and SlSDG16 showed expression patterns that fully matched the flower‐to‐fruit transition phase, we addressed the functional significance of these genes in tomato fruit set using genome editing via CRISPR‐Cas9 technology. No mutant lines could be generated for SlSDG5, suggesting that the mutation in this gene might be lethal or detrimental to the regeneration process during plant transformation. Several SlSDG16 KO lines were obtained (Fig. S9) but none of these mutants gave rise to detectable phenotypes, suggesting a potential functional redundancy within this large gene family. By contrast, two different mutations in SlSDG27 were obtained, both resulting in a frameshift predicted to produce a truncated protein with a deletion of the complete SET domain and a partial deletion of the PHD domain (Fig. 6c). Consistent phenotypes related to fruit set were observed in these mutants exhibiting pollination‐independent fruit formation (Fig. 6c). However, only heterozygous lines for these SlSDG27 KO‐mutants could be selected, presumably due to the lethal effect of this mutation at the homozygous state. The parthenocarpic fruit produced by these lines suggested the potential ability of SlSDG27 to trigger the flower‐to‐fruit transition independently from successful fertilisation of the flower. However, the seedless fruit did not allow obtaining progenies from these lines to undertake further characterisation.

Fig. 6

Histone modifier genes during tomato fruit set. (a) Differential expression and histone mark association of epigenetic‐related genes during fruit set. The upper panel of the graph corresponds to the expression profile, where x‐axis refers to the gene name annotated either based on tomato referenced studies or, when missing, according to the best corresponding orthologue in Arabidopsis. Genes with significant differential expression were marked by asterisks (*fold > 2 and 0.01 < P‐value < 0.05; **fold > 2 and 0.001 < P‐value < 0.01; ***fold > 2 and P‐value < 0.001). The lower panel corresponds to the heatmap of DEGs and DAs. The pink blocks indicate an increase in gene expression (fold > 2 and P‐value < 0.01) or histone mark association (P < 0.01); the green blocks indicate a decrease in gene expression (fold > 2 and P‐value < 0.01) or histone mark association (P < 0.01). The blue blocks indicate that at least two DA regions are found in the same gene and these DAs show both gain and loss of histone marks (P < 0.01). (b) Heatmaps showing absolute mean normalised expression values. The expression patterns of 19 DE histone modifier genes (left panel) were assessed using expression datasets available in TomExpress platform which includes several developmental stages and tissue types from wild species S. pimpinellifolium (middle panel) and Solanum lycopersicum (right panel). Of these, six genes showing consistent expression patterns among various experimental contexts are highlighted red to signify the potential importance of these histone modifiers in the control of the fruit set. (c) Knockout of SlSDG27 in tomato leads to seedless fruit formation in CR#SDG27 plants generated by CRISPR/Cas9 gene editing. Two guide RNAs (sgRNA1 and sgRNA2; green bars) were designed for editing the target gene. Protospacer‐adjacent motif (PAM) are indicated in blue letters. Two independent mutations within SlSDG27 sequence are shown in red and sequence gaps represented by dashes. The predicted truncated proteins are schematically illustrated to show they are impaired in PHD and SET functional domains. The lower panel shows parthenocarpic fruits in heterozygous lines representative of the phenotypes displayed by the two SlSDG27 mutants. Bars, 1 cm.

Discussion

Epigenetic control operates either through repositioning of histone marks or changes in DNA methylation, however the respective contributions of the two mechanisms to the dynamic changes of gene expression underlying various plant growth and development processes have been only partially explored. The flower‐to‐fruit transition provides a case study that is well suited for comparatively addressing the contribution of DNA methylation and HPTMs to the transcriptomic reprogramming underlying a developmental transition that has a major impact on crop yield. The outcome of the present study supported a scenario in which histone modification associated with transcriptional reprogramming underpinned the initiation of tomato fruit development. This is different from previous reports emphasising the prominent role of DNA methylation in reprogramming key genes during the transition to ripening of tomato fruit (Zhong et al., 2013; Liu et al., 2015a; Lang et al., 2017). Notably, HPTMs are strongly associated with the global transcriptomic changes during the flower‐to‐fruit transition in both pollination‐dependent and ‐independent fruit set.

Our study indicated that, during fruit set, DEGs are mainly associated with changes in H3K9ac or H3K4me3 marking, in line with the dynamic histone modifications reported in several developmental transition processes such as Arabidopsis meristem differentiation (Lafos et al., 2011), shoot apical meristem to inflorescence meristem transition (You et al., 2017), floral morphogenesis (Engelhorn et al., 2017), diurnal rhythms regulation (Song et al., 2019) and stomatal guard cell differentiation (Lee et al., 2019). By contrast, our data showed that changes in H3K27me3 marks correlated weakly with differential expression of individual genes, similar to the situation during the transition from shoot apical meristem to inflorescence meristem in rice and Arabidopsis in which, in most instances, gain or loss of H3K27me3 did not show correlation with changes in gene expression (Liu et al., 2015b; You et al., 2017). Moreover, in rice lines downregulated in H3K27me3 methytransferase SDG711 expression, most genes that displayed loss of H3K27me3 mark did not show increased expression levels (Liu et al., 2015b). Similarly, in Arabidopsis, only a small number of genes displayed a change in expression upon loss of H3K27me3 marks in PRC mutants or gain of H3K27m3 in demethylase mutants (Lafos et al., 2011; Lu et al., 2011; Shu et al., 2019). This supports the view that depletion of H3K27me3 alone is not sufficient to promote gene expression (Lafos et al., 2011; Lu et al., 2011; Wang et al., 2016; Shu et al., 2019).

Although H3K27me3 repressive mark did not show clear correlation with changes in gene transcript levels, it may still play an active role in gene expression reprogramming during tomato fruit set as suggested by the downregulation of PcG components, SlEZ2 and FIE, that leads in both cases to altered fruit set in tomato (Liu et al., 2012; Boureau et al., 2016). Our data clearly indicated that increased gene expression largely correlated with gain of H3K9ac and H3K4me3 marks, irrespective of the change in H3K27me3 marks. This sustains the view that H3K9ac and H3K4me3 marks represented the main HPTMs events associated with gene expression reprogramming during the flower‐to‐fruit transition. It is however likely that other histone marks might also be involved in the transcriptional reprogramming underlying fruit set, given that c. 20% of the expressed genes were devoid of any of the three histone marks taken into consideration in our study. Conversely, it is also possible that some active genes were not necessarily associated with chromatin regulation.

The data indicated that H3K9ac and H3K4me3 marks might act synergistically to promote gene transcription, while most H3K27me3‐marked genes exhibited low or no expression and were devoid of the two transcriptionally permissive histone marks (Fig. 3e). Highly expressed genes showed extremely low association with H3K27me3, but high enrichment in H3K9ac and/or H3K4me3 suggested that H3K9ac and H3K4me3 marks were mutually exclusive of H3K27me3. Some genes, in which the three marks co‐exist, showed variable expression in different tissues of the ovary and which might reflect the mixed nature of the tissues used in the ChIP‐seq experiments. Although it cannot be ruled out that genes detected as associated with antagonistic histone marks may correspond to two versions of the same gene bearing distinct marks in two different cell or tissue types. However, sequential ChIP performed at the cellular level provided proof that several bivalent genes or regions can truly exist in plants and other organisms (Sequeira‐Mendes et al., 2014). It has been reported that the trithorax group (TrxG) and PcG, known to mediate H3K4 and H3K27 histone modifications, respectively, worked antagonistically to activate or repress overlapping sets of genes (Schuettengruber et al., 2007; Pu & Sung, 2015). In this regard, Arabidopsis AGMOUS and FLOWERING LOCUS C genes were shown to be a common target for H3K4me3‐TrxG and H3K27me3‐PcGs, thus conferring to these histone modifiers key roles in the shift from embryo to seedling or from vegetative to reproductive phases (Saleh et al., 2007; Pien et al., 2008). Some critical genes have been reported to carry bivalent marks in Arabidopsis and potato (Saleh et al., 2007; Berr et al., 2010; Sequeira‐Mendes et al., 2014; Zeng et al., 2019) and, accordingly, our data revealed the concomitant occurrence of H3K4me3, H3K9ac and H3K27me3 histone marks in a number of tomato homeobox transcription factor genes such as HB‐BELL, HB‐WOX, HB‐KNOX, HB‐HD‐ZIP and zf‐HD involved in vegetative (Reiser et al., 1995; Byrne et al., 2003; Mele, 2003; Du et al., 2009) or reproductive organ development such as TAGL1, TAGL11 and RIN MADS‐box (Itkin et al., 2009; Martel et al., 2011) (Table S12).

Several studies have reported the critical role of DNA demethylation during the ripening of tomato fruit through de‐repression of key ripening regulator genes such as CNR, RIN and NOR (Manning et al., 2006; Li et al., 2008; Zhong et al., 2013; Liu et al, 2015). This process is mediated by SlDML2 DNA demethylase, which was shown to exhibit high expression levels during the fruit ripening transition. By contrast, our data showed that SlDML2, as well as three other DNA demethylase genes, displayed extremely low expression levels during the flower‐to‐fruit transition, suggesting a minor contribution of DNA demethylation to changes in transcriptomic reprogramming during the fruit set process. Altogether, our study indicated that change in DNA methylation was weakly associated with the fruit set transcriptomic programme, although a small set of genes involved in the fruit set transition might be controlled by DNA methylation. In line with this idea, recent DNA methylome studies in soybean and Arabidopsis seeds showed that changes in cytosine methylation did not correlate with changes in transcript levels of genes extremely important for seed development and germination (Lin et al., 2017; Kawakatsu et al., 2017). Taking together, it seems that different developmental transitions in tomato fruit might be associated with different modes of epigenetic remodelling, even though HPTMs have not been thoroughly addressed in the case of fruit ripening,

The ChIP‐seq data revealed a net enrichment for the three histone marks in 4DPA and 4IAA compared with 0DPA samples (Fig. 4a), consistent with the decreased expression of histone deacetylase and the upregulation of histone methyltransferase genes revealed by RNA‐seq during fruit set. However, impairing the expression of histone methyltransferase genes that displayed expression patterns matching the fruit set transition in tomato failed to provide clear clues to the putative role of these histone modifiers in the fruit set process. In some cases, knockout mutations in these genes might be lethal or impair the regeneration of essential plant organs during the genetic transformation process, thus precluding obtaining corresponding mutants. In some other cases, functional redundancy among members of gene families may obstruct the observation of the phenotypes. Nevertheless, heterozygous lines bearing a knockout mutation within the SlSDG27 gene showed pollination‐independent fruit formation, suggesting a potential role for this gene in triggering the flower‐to‐fruit transition.

Overall, the present study sheds new light on the main events and regulatory mechanism underlying the fruit set transition, and provides novel targets for the design of breeding strategies, aiming to ensure yield stability in the face of climate change by targeting not only genetic variation but also epigenetic regulation.

Author contributions

GH and PF performed the experiments. MBe helped in setting ChIP analysis. GH, AD and MZ, performed the bioinformatic analysis. GH, EM and KW analysed the data. PG, BH and KW provided input and expertise. MB and MZ conceived and directed the project, and GH, MB, KW and MZ wrote the manuscript.

Data availability

Datasets supporting the conclusions of this article are available (study PRJEB19602 and PRJEB38607) from the European Nucleotide Archive (http://www.ebi.ac.uk/ena/data/) with the following accession numbers: ERS1572540 to ERS1572558 for RNA‐seq analysis; ERS1572559 to ERS1572570, for ChIP‐seq analysis; ERS4597399 to ERS4597404 for bisulfite‐seq analysis. All data processing details are provided in Methods S1.

Fig. S1 Cluster dendrogram of DNA cytosine methylation level in 0DPA, 4DPA and 4IAA samples.

Fig. S2 DNA methylation profiles in promoter of GA and auxin‐related genes differentially expressed during fruit set.

Fig. S3 Expression profile of genes involved in the regulation of DNA methylation during fruit set process.

Fig. S4 ChIP‐qPCR validation of ChIP‐seq data.

Fig. S5 Genomic distribution of H3K9ac, H3K4me3 and H3K27me3.

Fig. S6 Identification of histone modified regions in 4IAA sample.

Fig. S7 GO terms enriched in RNA‐seq and ChIP‐Seq datasets.

Fig. S8 Mean normalised expression of the six histone modifier genes selected.

Fig. S9 Generation of SlSDG16 knockout lines by CRISPR/Cas9 technology.

Methods S1 Data processing of RNA‐seq, ChIP‐seq and RNA‐seq data.

Click here for additional data file.

Table S1 Primers used in this study.

Click here for additional data file.

Table S2 Read mapping summary for RNA‐seq, ChIP‐seq, BS‐seq libraries.

Click here for additional data file.

Table S3 DA regions identified by MAnorm method.

Click here for additional data file.

Table S4 List of genes related to epigenetic regulation in tomato.

Click here for additional data file.

Table S5 List of differentially expressed genes (DEGs).

Click here for additional data file.

Table S6 List of GO enrichment biological processes from DEGs (P < 0.05).

Click here for additional data file.

Table S7 Identification of DMRs in CG, CHG and CHH sequence context.

Click here for additional data file.

Table S8 Size of associated region for each mark and percent of the tomato genome covered by the three histone marks.

Click here for additional data file.

Table S9 Differential expression of important fruit set‐related genes and their differential association with histone marks.

Click here for additional data file.

Table S10 Gene ontology analysis of DEGs and active mark‐associated‐DA genes.

Click here for additional data file.

Table S11 Hormone‐related DEGs and DA genes.

Click here for additional data file.

Table S12 Putative bivalent genes associated with both H3K9ac/H3K4me3 and H3K27me3.

Please note: Wiley Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Click here for additional data file.