Phytochromes B1/B2 Are Major Regulators of Ripening-associated Epigenome Reprogramming in Tomato Fruits

Phytochrome-mediated light and temperature perception has been shown to be a major regulator of fruit development. Furthermore, chromatin remodelling via DNA demethylation has been described as a crucial mechanism behind the fruit ripening process; however, the molecular basis underlying the triggering of this epigenetic modication remains largely unknown. Here, through integrative analyses of the methylome, siRNAome and transcriptome of tomato fruits from phyA and phyB1B2 null mutants, we report that PHYB1 and PHYB2 control genome-wide DNA methylation during fruit development. The experimental evidence indicates that PHYB1B2 signal transduction is mediated by the coordinated expression of DNA methylases/demethylases, histone-modifying enzymes and chromatin remodelling factors, resulting in the transcriptional regulation of photosynthetic and ripening-associated genes. This new level of understanding provides insights into the orchestration of epigenetic mechanisms in response to environmental cues affecting agronomical traits.


Introduction
As sessile organisms, plants must constantly monitor their environment and continuously tune their gene expression to enable adaptation and survival 1 . Light is one of the main environmental cues that controls plant growth and development from seed germination to senescence 2 . Plants employ different photoreceptors to detect and respond to changes in the incident spectral composition (from UV-B to farred wavelengths), light direction and photoperiod. These photoreceptor families include (i) phytochromes (PHYs), which perceive red/far-red (R/FR) light; (ii) cryptochromes (CRYs), phototropins, and 'Zeitlupes', which sense blue/UV-A light; and (iii) the UV-B receptor UVR8 3 .
After photoreceptor activation, complex signal transduction pathways control the expression of lightregulated genes via transcriptional, posttranscriptional, and posttranslational mechanisms 2 . Several hub proteins in the light signal transduction pathway triggered by PHYs, CRYs and UVR8 have been identi ed, including transcription factors (TFs) such as PHY-INTERACTING FACTORS (PIFs) and ELONGATED HYPOCOTYL5 (HY5) as well as the ubiquitin E3 ligase complex comprising CONSTITUTIVE PHOTOMORPHOGENIC1 (COP1) 2 . Additionally, the effect of light on alternative splicing (AS) has recently been described 4 . In Arabidopsis thaliana, 1,505 genes undergo changes in the AS levels of their transcripts within 1 h of exposure to red light in a PHY-dependent manner 5 . Similarly, in the moss Physcomitrella patens, 8.4% and 8.9% of AS events rapidly respond to PHY-mediated red and blue light perception 6 . Further studies showed that the interaction between PpPHY4 and a splicing factor coregulates 70% of intron retention (IR) events in response to red light in P. patens 7 . Moreover, light controls protein localization through PHY-mediated alternative promoter selection, allowing plants to metabolically respond to changing light conditions 8 . Finally, it has been extensively reported that activated PHYs induce post-translational changes in PIF proteins, including sequestration, phosphorylation, polyubiquitylation, and subsequent degradation through the 26S proteasome-mediated pathway 3 . Although the effect of light on plant phenotypes and the plant transcriptome has been studied for decades [9][10][11] , the involvement of epigenetic regulatory mechanisms in light-dependent changes in the transcriptional landscape remains poorly addressed.
Posttranslational histone modi cations, such as acetylation and methylation, have been associated with the induction and repression of light-responsive genes 12,13 . Light-dependent enrichment of the acetylation pattern of H3 and H4 in the enhancer and promoter regions of the pea plastocyanin locus PetE has been reported 14 , and the hyperacetylation of the PetE promoter is linked to the transcriptional activity of this gene 15 . Moreover, a reduction in H3 acetylation is associated with a decrease in the expression of the A. thaliana light-responsive genes CHLOROPHYLL a/b-BINDING PROTEIN GENE (CAB2) and the RIBULOSE BISPHOSPHATE CARBOXYLASE/OXYGENASE small subunit (RBCS) 16 . Histone methylation regulates PHY-mediated seed germination in A. thaliana. Upon R light illumination, photoactivated PHYB (Pfr) targets PIF1 for proteasome-mediated degradation, releasing the expression of the JUMONJI HISTONE DEMETHYLASES JMJ20 and JMJ22. As a result, JMJ20 and JMJ22 reduce the levels of H4R3me2, which leads to the activation of the gibberellic acid biosynthesis pathway to promote seed germination 17 .
Together with histone modi cation, DNA methylation is a common epigenetic mark with a direct impact on gene expression. Nevertheless, only a few reports have speci cally addressed the effect of light stimuli on DNA methylation. Light-dependent nuclear organization dynamics during deetiolation are associated with a reduction in methylated DNA 18 . In Populus nigra, 137 genes were shown to be regulated by methylation during the day/night cycle 19 . Moreover, photoperiod-sensitive male sterility is regulated by RNA-directed DNA methylation (RdDM) in rice 20 . Finally, using a methylation-sensitive ampli ed polymorphism assay, DNA methylation remodelling was shown to be an active epigenetic response to different light qualities in tomato seedlings 21 .
Here, we investigate the impact of DNA methylation on gene expression regulation in response to PHYdependent light perception in tomato fruits. The Solanum lycopersicum genome harbours ve PHYencoding genes: PHYA, PHYB1, PHYB2, PHYE and PHYF 22 . Previous studies have shown that PHYA, PHYB1 and PHYB2 are major regulators of fruit chloroplast maturation and nutraceutical compound accumulation [23][24][25][26] . In this work, the genome-wide transcriptome and methylome were comprehensively analysed in fruits from phyA and phyB1B2 null mutants, revealing that PHY-dependent DNA demethylation is a crucial stimulus for triggering ripening-associated gene expression. Moreover, our results showed that gene-body (GB) RdDM is positively correlated with gene expression during fruit development.

Impact of light perception impairment on the fruit transcriptome
To investigate the role played by either PHYA or PHYB1/PHYB2 (hereafter PHYB1B2) in overall gene expression during fruit development, the transcriptome of fruits at the immature green (IG) and breaker (BK) stages from phyA and phyB1B2 null mutants as well as their wild-type (WT) counterpart, was determined by RNAsEq. Among the approximately 20,000 transcriptionally active loci in each biological replicate (Supplementary Table 1), 1.2% and 2.4% at the IG stage and 9.1% and 11.2% at the BK stage were identi ed as differentially expressed genes (DEGs) in phyA or phyB1B2 mutants, respectively, compared to the WT ( Fig. 1a; Supplementary Table 2). For both genotypes, the number of exclusive DEGs was signi cantly lower in the IG stage than in the BK stage; similarly, the number of genes that were commonly regulated by PHYA and PHYB1B2 was 172 at the IG stage and 785 at the BK stage (Fig. 1b). Subsequently, the altered expression of approximately 76% (23/30) of the tested genes was validated by RT-qPCR (Supplementary Table 3). Comparison with previously reported expression data for genes involved in ripening regulation, ethylene biosynthesis and signalling, and carotenogenesis further validated our RNAseq data, as 90% of the analysed genes on average showed the expected mRNA pro le at IG and BK stages. It is worth mentioning that most of the genes displayed the same transcript uctuation in the WT, phyA and phyB1B2 genotypes, though this was somewhat attenuated in the mutants (Supplementary Table 4). These results showed that PHY-meditated light perception regulates more genes in BK than in the early stages of fruit development and that PHYB1B2 has a more substantial impact than PHYA in the fruit transcriptome in both analysed stages.
A closer look at DEGs function revealed a similar distribution of loci across MapMan categories in response to phyB1B2 and phyA mutations in both developmental stages, although with distinct abundance levels ( Fig. 1c). At the IG stage, eight categories were mainly represented, including at least 2% of the DEGs identi ed in phyA and phyB1B2: photosynthesis, lipid metabolism, phytohormone action, RNA biosynthesis, protein modi cation, protein homeostasis, cell wall organization, and solute transport ( Fig. 1c; Supplementary Tables 5 and 6). It is worth highlighting the abundance of the DEGs within the photosynthesis category in the phyB1B2 mutant, among which 34 out of the 37 genes were downregulated (Supplementary Table 6). In the BK stage, at least 2% of the DEGs were related to the lipid metabolism, phytohormone action, RNA biosynthesis, protein modi cation and homeostasis, cell wall organization and solute transport categories in both genotypes ( Fig. 1c; Supplementary Tables 7 and 8). However, while phyA de ciency also affected carbohydrate metabolism and external stimuli (Supplementary Table 7), the phyB1B2 mutant showed a large number of DEGs related to the cell cycle and chromatin organization (Supplementary Table 8). Interestingly, the chromatin organization category displayed 52 DEGs, 45 of which were upregulated. These genes encode nucleosome constituent histones (H3, H4, H2A and H2B); DNA methylases/demethylases; histone post-translational modi ers such as deacetylases, methylases/demethylases, histone ubiquitination factors and histone chaperones; chromatin remodelling factors; and genes involved in RNA-independent and RNA-directed DNA methylation (Supplementary Table 8). These results led us to further investigate the impact of DNA methylation on PHY-mediated gene expression reprogramming.

PHYs regulate the epigenome pro le
The global pro le of methylated cytosines (mCs) in the epigenome of tomato fruits was assessed by whole-genome bisul te sequencing in the IG and BK stages for phyA, phyB1B2 and WT genotypes. In agreement with previous reports 27,28 , regardless of the genotype and fruit stage, the greatest total number of mCs was located in the CHH context, followed by the CG and CHG contexts, while the methylation level was highest in the CG (80%) context followed by the CHG (67%) and CHH (23%) contexts (Supplementary Table 9, Supplementary Fig. 1). For further comparisons, we selected only cytosines with coverage > 10X, and except for chromosome 9 in the transposable elements TEs enriched region, all samples met this cutoff. In all contexts, the highest cytosine density was associated with gene-rich euchromatic regions located at chromosome arm ends ( Supplementary Fig. 1). Conversely, in symmetrical contexts (CG and CHG), the highest methylation rates were found across pericentromeric regions enriched in TEs and in the CHH context in gene-rich regions associated with a higher density of sRNAs ( Supplementary Fig. 1). The comparison of the methylation status between the two fruit stages showed that ripening-associated demethylation 27 occurs mainly in the CG context, especially in gene-rich regions, and that it is impaired in phyB1B2 mutant BK fruits ( Supplementary Fig. 1).
The subsequent comparison between genotypes revealed global epigenome alteration in phy mutants in all contexts analysed. The most remarkable observation was the presence of considerable hypermethylation in all contexts across gene-rich regions in BK-stage fruits from phyB1B2 (Fig. 2a). In contrast, phyA exhibited hypermethylation in CHG and CHH contexts associated with TE-rich regions (Fig. 2a), suggesting that different PHYs control DNA methylation across speci c genomic regions through distinct regulatory mechanisms. Interestingly, PHY-associated hypomethylation was exclusively detected in the CG context of gene-rich regions in IG-stage fruits from phyA and in the CHH context of TErich regions for BK-stage fruits from phyB1B2. In summary, these data revealed that both PHYA and PHYB1B2 regulate the global methylome, but PHYB1B2 has a greater impact on ripening-associated methylation reprogramming across gene-rich genomic regions in tomato fruits.
To investigate the relationship between PHY-regulated cytosine methylation and gene expression, we rst identi ed genes with differentially methylated promoters (DMPs, 2 kb upstream of the transcription start site) in all three contexts. Interestingly, associated with the massive alteration previously observed, the pattern of DMPs varied with the mC context, stage and genotype (Fig. 2b, Supplementary Tables 10 and 11). Regarding the CG context, whereas the phyA mutant showed virtually the same frequency of hyperand hypomethylated promoters in the two stages, the status of hypermethylated promoters in phyB1B2 increased over 60% from the IG to BK stage, while the number of loci with hypomethylation decreased 50% (Fig. 2b, Supplementary Table 12). In contrast, phyA showed a greater number of hypermethylated promoters in the CHG context in the IG stage than in the BK stage, while the levels in the WT and phyB1B2 mutant remained similar upon ripening (Fig. 2b, Supplementary Table 13). In the CHH context, the number of hypermethylated promoters decreased in both genotypes from the IG to BK stages (Fig. 2b, Supplementary Table 14).
These results indicate that PHY de ciency results in massive promoter hypermethylation in both the IG and BK stages of tomato fruit development. Moreover, they reinforce the role of PHYB1B2 in ripening-associated demethylation and its putative effect on gene expression.
Effect of PHY-mediated differential methylation on the transcriptome To assess whether the differential methylation of gene promoters affects mRNA levels, we crossed data from DEGs and DMPs between genotypes. Supplementary Fig. 2 shows scatter plots of promoter methylation vs mRNA fold changes for comparisons of the two genotypes at the two examined developmental stages in the three mC contexts. The most evident result was that among the thousands of loci with identi ed DMPs (Fig. 2b), only hundreds of the loci were also differentially expressed (Supplementary Table 15) (0.7% for IG phyA, 1.6% for IG phyB1B2, 5.6% for BK phyA and 7.4% for BK phyB1B2), raising an intriguing question about the biological signi cance of the extensive change in the methylation pattern observed in the mutants. In contrast, the percentages of the DEGs showing DMPs were 73% for IG phyA, 76% for IG phyB1B2, 72% for BK phyA and 75% for BK phyB1B2 ( Supplementary  Fig. 2). Many more DEGs with DMPs were observed in BK than in IG fruits and in phyB1B2 than in the phyA genotype. The functional categorization of these genes revealed a similar category distribution to the DEGs (Fig. 1c, Supplementary Table [16][17][18][19]. At the IG stage, there were seven categories in which at least 2% of the loci showed DMPs and differential expression in both genotypes: photosynthesis, phytohormone action, RNA biosynthesis, protein modi cation and homeostasis, cell wall organization and solute transport, whereas phyB1B2 additionally impacted lipid metabolism (Fig. 1c). In the BK stage, the categories in which at least 2% of the DEGs showed DMPs were lipid metabolism, phytohormone action, RNA biosynthesis, protein modi cation and homeostasis, cell wall organization and solute transport-related functions in both genotypes, while only phyA impacted carbohydrate metabolism and external stimuli, and only phyB1B2 affected photosynthesis, chromatin organization and cell cycle categories.
Interestingly, in the comparison of the IG and BK stages, 42.5%, 34.2% and 18.8% of the DMPs were associated with DEGs, while 79.5%, 76.6% and 71.5% of the DEGs showed differences in promoter methylation in WT, phyA and phyB1B2, respectively (Supplementary Fig. 3). All of these data demonstrated that the PHY-dependent mRNA pro le is profoundly affected by promoter methylation; however, massive genome-wide PHY-induced methylation reprogramming has a still uncharacterized role beyond the regulation of mRNA accumulation. Moreover, promoter methylation has a greater effect on gene expression regulation during BK than in the IG stage. Additionally, the data showed that PHYB1B2 has a more extensive in uence on gene expression regulated via promoter methylation than PHYA, reinforcing the above conclusion that PHYB1B2 affects CG ripening-associated demethylation ( Supplementary Fig. 2).

The sRNAome is altered by PHY de ciency
To assess the involvement of RdDM in PHY-mediated transcriptome regulation, the sRNAome was analysed in fruits at the IG and BK stages from both mutants and the WT genotype (Supplementary   Table 20a). A total of 28,314 clusters of sRNAs were identi ed across the whole genome in at least one of the samples, including 7,984 in gene bodies, 7,863 in promoter regions, 7,966 in TEs and the remaining 4,501 across intergenic regions ( Supplementary Fig. 1, Supplementary Table 20b). The methylation level was evaluated for each sRNA cluster-targeted genomic region (sCTGR), and as previously observed for promoter regions, a higher proportion of hypermethylation was observed in BK fruits from phyB1B2 in the CG symmetrical context. Moreover, the greatest number of differentially methylated sCTGRs was observed in the asymmetrical context CHH. (Fig. 3a, Supplementary Table 20g-j).
sCTGR methylation levels and sRNA accumulation data were intersected, and among a total of 154, 318, 267 and 257 differentially accumulated sRNA clusters (Supplementary Table 20c-f), 88, 154, 99 and 82 also showed differential methylation changes (> 5%) in phyA IG, phyB1B2 IG, phyA BK and phyB1B2 BK fruits, respectively (Fig. 3b, (Supplementary Table 20g-j), showing a strong association (P < 0.005) of the two datasets. Intriguingly, this positive association was not observed in the transition from the IG to BK stages ( Supplementary Fig. 4), suggesting that the global methylation changes via RdDM could be attributed to PHY de ciency. Moreover, a clear disturbance in sRNA accumulation was observed in phyB1B2, since almost no clusters with less sRNA accumulation were observed in BK compared to the IG stage ( Supplementary Fig. 4). We further analysed whether this positive correlation observed between sRNA accumulation and sCTGR methylation impacted gene expression levels. Notably, regardless of the fruit developmental stage, changes in the accumulation of sRNA located in gene bodies (GBs), and not in the promoter (P) region, were positively correlated with the mRNA level (Fig. 3c, Supplementary Table 20kn). Among these loci, two interesting examples were identi ed: the well-known ripening-associated genes RIPENING INHIBITOR (RIN, Solyc05g012020 29 ) and FRUITFULL2 (FUL2, Solyc03g114830 30 ), which showed higher expression in phyB1B2 at the IG stage (Fig. 4a) and higher sRNA accumulation and sCTGR methylation across their GBs (Fig. 4b) compared to WT. The premature expression of these TFs was in agreement with the previously reported anticipation of ripening onset in the phyB1B2 mutant 31 . Altogether, these ndings revealed (i) impaired RdDM in BK fruits of phyB1B2, indicated by the absence of clusters with less sRNA accumulation, and (ii) that GB RdDM is an important mechanism that positively regulates gene expression in a PHY-mediated manner during fruit development.

PHYB1B2-dependent methylation regulates fruit chlorophyll accumulation
The categorization of DEGs associated with differential promoter methylation revealed prominent representation of the photosynthesis category in the fruits of the phyB1B2 mutant at the IG stage (Fig. 1c). Among the 32 genes, 22 were downregulated and hypermethylated in the promoter region (Supplementary Tables 6 and 17). Most of these genes encode chlorophyll-binding proteins, structural photosystem proteins and chlorophyll biosynthetic enzymes. This might at least partly explain the reduction of 50% in the total chlorophyll level observed in phyB1B2 IG fruits (Fig. 5a). Detailed promoter analysis of the chlorophyll-related enzymes PROTOCHLOROPHYLLIDE OXIDOREDUCTASE 3 (POR3, Solyc07g054210), and two CHLOROPHYLL A/B BINDING PROTEINs (CBP, Solyc02g070990 and CAB-3c, Solyc03g005780) showed that the reduced mRNA levels of these three genes in phyB1B2 (Fig. 5b) were correlated with the presence of hypermethylated regions in the promoters. Interestingly, several of these hypermethylated sequences overlapped with HY5 and PIF photomorphogenic TF binding sites, such as Gbox (CACGTG), CA hybrid (GACGTA), CG hybrid (GACGTG) and PBE-box (CACATG) motifs 32 (Fig. 5c). These results suggest that the transcription of genes involved in chlorophyll metabolism and the photosynthetic machinery in tomato fruits is regulated by the PHYB1B2-mediated methylation status of their promoter regions in addition to the PHY-mediated post-translational regulation of HY5 and PIF protein levels 2 .
The methylation-mediated regulation of fruit ripening is PHYB1B2-dependent In their seminal study, Zhong et al. (2013) 27 revealed that the extensive methylation in the promoter regions of ripening-associated genes gradually decreases during fruit development. Interestingly, RNA biosynthesis, which includes transcription factors, was the most abundant functional category among the DEGs that showed DMPs (Fig. 1c). Thus, we examined a set of ripening-associated master transcription factors: RIN, Solyc05g012020, NON-RIPENING (NOR, Solyc10g006880 33 ), COLORLESS NORIPENING (CNR, Solyc02g077920 34 ) and APETALA2a (AP2a, Solyc03g044300 35 ). The evaluation of the promoter regions clearly showed that while their methylation level decreased from the IG to BK stage in the WT genotype, they remained highly methylated in phyB1B2 (Fig. 6a). As previously observed, several hypermethylated regions are closely linked to HY5 and PIF binding sites. The maintenance of high methylation levels in the promoters of these key regulatory genes at the onset of fruit ripening was highly correlated with their transcriptional downregulation at the BK stage (Fig. 6b).
Carotenoid accumulation is probably the most appealing and best investigated trait of tomato fruits; in agreement with previous ndings 26 , ripe phyB1B2 fruits showed a ve-fold reduction in carotenoid content compared to WT (Fig. 7a).
With the aim of evaluating whether this effect is a consequence of the methylation-mediated regulation of carotenoid biosynthesis genes, we further analysed the promoters of PHYTOENE SYNTHASE 1 (PSY1, Solyc03g031860), PHYTOENE DESATURASE (PDS, Solyc03g123760), 15-CIS-ζ-CAROTENE (ZISO, Solyc12g098710) and ZETA-CAROTENE DESATURASE (ZDS, Solyc01g097810), which, with the exception of PDS, were hypermethylated in phyB1B2 BK fruits (Supplementary Table 11). The mC pro le con rmed the presence of hypermethylated regions in all four promoters, which were predominantly located near light-dependent TF binding motifs in the PSY1, PDS and ZISO promoters (Fig. 7b), explaining the reduced mRNA levels of these genes observed in phyB1B2 (Fig. 7c).
RIN is one of the main TFs controlling ripening-associated genes by directly binding to their promoters.
RIN binding occurs in concert with the demethylation of its targets 27 . To examine whether RIN binding site methylation could be affected by the phyB1B2 mutation in the ripening-related master transcription factors and carotenoid biosynthesis gene promoters, we mapped the available RIN ChIP-seq data 27 and performed de novo motif discovery ( Supplementary Fig. 5). Interestingly, the levels of mCs associated with RIN targets were higher in the NOR, CNR and AP2a promoters in phyB1B2 than in WT. Moreover, the RIN promoter itself was hypermethylated across the RIN binding site in phyB1B2 BK fruits, suggesting a positive feedback regulatory mechanism (Fig. 6a). Finally, in the phyB1B2 mutant, the PSY1, PDS, ZISO and ZDS promoters showed higher methylation overlapping with RIN target binding motifs, indicating that the upregulation of carotenoid biosynthesis genes during tomato ripening is dependent on the PHYB1B2mediated demethylation of RIN target sites.
Altogether, our ndings showed that PHYB1B2 is a major player in fruit ripening by controlling the promoter demethylation of master transcriptional regulators and carotenoid biosynthesis genes.

Discussion
The dynamic methylation pattern during tomato fruit development has been demonstrated to be a critical ripening regulation mechanism 27,28 . DNA demethylation, mainly in the CG context, triggers the activation of genes involved in ripening and is required for pigment accumulation and ethylene synthesis 27,36 . Simultaneously, the dynamic epigenome during fruit development is strictly regulated by environmental cues 37 . The prevailing model establishes PHYs as major components involved in the coordination of fruit physiology with the ever-changing light and temperature environment 23,26 . Thus, we explored the link between fruit epigenome reprogramming and these well-established light and temperature sensors.
Our data clearly showed that phyA and phyB1B2 de ciencies modi ed the epigenome pro le through methylome and sRNAome reprogramming. In particular, PHY-mediated DMPs and GB methylation impacted the transcriptome pattern, affecting tomato fruit development and demonstrating that a portion of the ripening-associated demethylation previously reported 27 is dependent on active PHYs, especially PHYB1B2. However, the massive alteration of methylation patterns observed in PHY mutants suggests the existence of a still unclear genome regulatory mechanism. Whether it is related to light-induced AS 3 and/or PHY-mediated alternative promoter usage to control protein localization 8 , as mentioned above, remains to be explored.
The phyA and phyB1B2 mutants showed a positive correlation between cluster sRNA accumulation, target methylation in GB and mRNA levels. In angiosperms, GB methylation has been associated with constitutively expressed genes 38,39 ; however, PHY de ciency intriguingly seems to deregulate this mechanism, affecting temporally and spatially regulated genes. The RIN and FUL2 examples analysed here clearly showed that sRNA accumulation and methylation were mainly located near transposable elements (TEs) (Fig. 4). It is known that the insertion of TEs within GB can disrupt gene expression; thus, methylation-mediated TE silencing and GB methylation are evolutionarily linked 38 . The enhancement of TE-associated DNA methylation in GB (Fig. 3c) and the absence of clusters with less sRNA accumulation in BK compared to the IG stage in phyB1B2 ( Supplementary Fig. 4) can be explained by the overexpression of canonical RdDM genes. Solyc12g008420 and Solyc06g050510 encode homologs of RNA-DEPENDENT RNA POLYMERASE (RDRP) and the associated factor SNF2 DOMAIN-CONTAINING PROTEIN CLASSY 1 (CLSY1), respectively, both of which were upregulated in BK fruits from phyB1B2 plants (Supplementary Table 2). These proteins are key players in canonical RdDM in A. thaliana required for PolIV-dependent sRNA production 40 . Similarly, Solyc09g082480 and Solyc03g083170, which were also upregulated in phyB1B2 BK fruits, are homologs of A. thaliana RNA-DIRECTED DNA METHYLATION 1 (RDM1) and DEFECTIVE IN MERISTEM SILENCING 3 (DMS3), respectively. The protein products of these genes, together with DEFECTIVE IN RNA-DIRECTED DNA METHYLATION 1 (DRD1), form the DDR complex, which enables RNA Pol V transcription 41 . To our knowledge, this is the rst report to associate PHY-mediated sRNA accumulation and DNA methylation with mRNA levels in plants.
Several pieces of evidence have shown that PHYB1B2 has a more substantial impact on genome regulation than PHYA. For example, BK fruits from phyB1B2 displayed (i) a large number of DEGs associated with chromatin organization (Fig. 1c); (ii) overall promoter hypermethylation in the CG context (Fig. 2b); (iii) the highest number of DEGs associated with DMPs ( Supplementary Fig. 2); and (iv) half the number of DMPs associated with DEGs between the IG and BK stages compared to the WT ( Supplementary Fig. 3). These massive epigenomic alterations in phyB1B2 led us to look more closely at genes related to chromatin organization that showed alterations in mRNA levels.
The chromomethylase SlMET1L (Solyc01g006100) (also referred to as SlCMT3 42 ) displays the highest transcript abundance in immature fruits, which declines towards the fully ripe stage 43 . Hence, in line with the higher level of DNA methylation, our transcriptome analysis showed that SlMET1L was upregulated in phyB1B2 BK fruits. Conversely, the upregulation of the SlROS1L demethylase (Solyc09g009080 43 ; also referred to as SlDML1 44 ) was also observed in phyB1B2 BK fruits, which might seem contradictory at rst glance. However, it has been reported that the Arabidopsis thaliana ROS1 gene promoter contains a DNA methylation monitoring sequence (MEMS) associated with a Helitron transposon, which is methylated by AtMET1, positively regulating AtROS1 gene expression 45 . Similarly, SlROS1L harbours two transposable elements within its promoter and showed a higher methylation level in phyB1B2 than in the WT genotype, suggesting a similar regulatory mechanism in tomato ( Supplementary Fig. 6, Supplementary Table 15).
The tomato homologue of A. thaliana DECREASED DNA METHYLATION 1 (DDM1, Solyc02g085390) showed higher mRNA expression in phyB1B2 mutant BK fruits than in their WT counterparts. DDM1 is a chromatin remodelling protein required for maintaining DNA methylation in the symmetric cytosine sequence 46 , which can be associated with the observed methylation differences in the CG context in phyB1B2 (Fig. 2a).
Several histone modi ers showed altered expression in BK fruits from the phyB1B2 mutant (Supplementary Table 8). The methylation of lysine residues 9 and 27 on histone H3 is associated with repressed genes. Histone lysine methyltransferases are classi ed into ve groups based on their domain architecture and/or differences in enzymatic activity 47 . The BK fruits of the phyB1B2 mutant displayed three differentially expressed lysine methyltransferases: Solyc03g082860, an upregulated H3K27 Class IV homologue, and Solyc06g008130 and Solyc06g083760, two H3K9 Class V homologs that showed lower and higher expression than in WT fruits, respectively. Histone arginine methylation is catalysed by a family of enzymes known as protein arginine methyltransferases (PRMTs). Solyc12g099560, a PRMT4a/b homologue, was upregulated in phyB1B2 BK fruits. Interestingly, in A. thaliana, PRMT4s modulate key regulatory genes associated with the light response 48 , reinforcing the link between the PHYB1B2 photoreceptor and epigenetic control. Finally, tomato histone demethylases have been recently identi ed. SlJMJ6, whose expression peaks immediately after the BK stage, has been characterized as a positive regulator of fruit ripening by removing the H3K27 methylation of ripening-related genes, and SlJMJ6-overexpressing lines show increased carotenoid levels 49 . SlJMJC1 (Solyc01g006680), which exhibits the same expression pattern 49 , is downregulated in the phyB1B2 mutant, suggesting that this gene might exhibit a similar regulatory function to its paralogue, inducing ripening in a PHYB1B2dependent manner (Figs. 6 and 7).
Histone deacetylation plays a crucial role in the regulation of eukaryotic gene activity and is associated with inactive chromatin 40 . Histone deacetylation is catalysed by histone deacetylases (HDACs). Fifteen HDACs were identi ed in the tomato genome 50 . Among these HDACs, SlHDA10 (Solyc01g009120) and SlHDT3 (Solyc11g066840) were found to be downregulated and upregulated in phyB1B2 BK fruits, respectively. SlHDA10 is localized in the chloroplast, and its transcript is highly expressed in photosynthetic tissues 50 ; whether SlHDA10 deacetylates chloroplast proteins by silencing photosynthesis-related genes remains to be determined. Although SlHDT3 is mainly expressed in immature stages of fruit development and its expression declines with ripening, its silencing results in delayed ripening and reduced RIN expression and carotenogenesis. On the other hand, the expression level of SlHDT3 is increased in ripening-de cient mutants such as Nr or rin 51 . Our results showed that phyB1B2 mutant fruits displayed higher expression of SlHDT3 and reduced RIN transcript levels at the BK stage, suggesting reciprocal regulation between these two factors (Fig. 8a). During the IG stage, SlHDT3 is highly expressed, contributing to the epigenetic inhibition of ripening. The reduction in SlHDT3 expression towards BK releases DNA methylation and both directly and indirectly upregulates RIN. Additionally, PHYB1B2 inhibits the expression of SlHDT3 and induces that of RIN, thereby inducing or repressing DNA methylation, respectively. Finally, RIN maintains SlHDT3 downregulation in BK 51 . This mechanism precisely tunes ripening-related epigenetic reprogramming and contributes to explaining the high methylation levels observed in the phyB1B2 mutant (Fig. 2).
The link between DNA methylation levels and tomato fruit ripening-associated gene expression has been previously reported 27,28 , but the stimuli and the molecular mechanisms underlying this relationship remained unknown. The integrated analysis of the experimental evidence together with previous gene functional studies in tomato and A. thaliana allowed us to propose that PHYB1B2 is an important triggering factor for chromatin remodelling and, consequently, transcriptional regulation during fruit development. PHY signal transduction is mediated by the coordinated expression of DNA methylases/demethylases, histone-modifying enzymes and chromatin remodelling factors resulting in the induction of photosynthesis and ripening-related genes in immature and breaker fruit stages, respectively (Fig. 8b). The vast reservoir of data released here brings a new level of understanding about how epigenetic mechanisms orchestrate the response to light and temperature uctuations affecting important agronomical traits in eshy fruits.

Material And Methods
Plant material, growth conditions and sampling phyA and phyB1B2 phytochrome mutants in the Solanum lycopersicum (cv. MoneyMaker) genetic background were previously characterized [52][53][54] . Tomato seeds were grown in 9L pots containing a 1:1 mixture of commercial substrate and expanded vermiculite, supplemented with 1 g L − 1 of NPK 10:10:10, Five replicates per genotype were cultivated. Fruits were sampled at the immature green (15 mm diameter), mature green (when the placenta displays a gelatinous aspect), breaker (beginning of ripening process when the fruit shows the rst yellowish colouration) and red ripe (7 days after the breaker stage) stages. All fruits were harvested at the same time of day with four biological replicates (each replicate was composed of a single fruit per plant). The columella, placenta, and seeds were immediately removed, and the remaining tissues were frozen in liquid nitrogen, ground and freeze-dried for subsequent analysis.  Table 1) and were used for statistical analysis.

Reverse transcription quantitative PCR (RT-qPCR)
Total RNA extraction was performed with the ReliaPrep™ RNA Cell and Tissue Miniprep System (Promega), and cDNA synthesis was conducted with SuperScript™ IV Reverse Transcriptase (Invitrogen). The primers used for qPCR are listed in Supplementary Table 21. RT-qPCR was performed in a QStudio6 Absolute uorescence data were analysed using LinRegPCR software to obtain Ct and primer e ciency values. Relative mRNA abundance was calculated and normalized according to the ΔΔCt method using EXPRESSED and CAC as reference genes 57 .

MethylC-Seq analysis
Methyl C sequencing was performed as described in a previous report 58 Table 9).
The Bioconductor package methylKit 60 was used for the detection of methylation levels across the analysed regions: promoters (2 kb upstream of transcription start site) and sRNA cluster-targeted genome regions (sCTGRs). Only Cs with 10X coverage were considered. Methylation differences with an FDR < 0.05 in each comparison (WT vs phyA; WT vs phyB1B2) were recorded as differentially methylated promoters (DMPs) or differentially methylated sCTGRs. Differential methylation in the CG, CHG and/or CHH context was considered to exist if the region contained at least 10 differentially methylated Cs in the corresponding context. Finally, for the comparison of global methylation levels between genotypes, only common Cs with at least 10X coverage in all samples were analysed. sRNAome pro le sRNA extraction and quality parameters were determined from three independent biological replicates of green and breaker stage fruits from each genotype, as described above in the "Transcriptional pro le" section. After RNA integrity con rmation, libraries were prepared using a TruSeq Small RNA Library Prep and sequenced using the Illumina HiSeq 4000 platform to generate a read length of 50 bp. The raw sequencing reads that were generated were quality trimmed with Trimmomatic 55 Table 20a). All libraries were aligned to genome version SL3.0 using ShortStack v3.8.1 61 with default parameters (allowing the distribution of multimapping reads according to the local genomic context). Then, the de novo identi cation of clusters of sRNAs was performed for all libraries, and individual counts for each library and cluster were obtained using the same software.

Statistical analysis for RNAseq and sRNAome
Genes/sRNA clusters with read/count numbers smaller than two per million were removed. Read/count values were normalized according to the library size factors. Statistical analyses were performed with edgeR from Bioconductor® 62,63 using a genewise negative binomial generalized linear model with the quasi-likelihood test 64 and a cutoff of the false discovery rate (FDR) ≤ 0.05.

Gene functional categorization
The DEGs were functionally categorized with MapMan application software 65 followed by hand-curated annotation using MapMan categories.

In silico regulatory motif predictions and RIN ChIP-seq analyses
To predict the transcription factor binding motifs, we used the complete collection of 530 plant transcription factor-binding sites (TFBS) modelled as position frequency matrices (PFMs) from the JASPAR 2020 database 66 . Putative binding sites were obtained by scanning the whole S. lycopersicum 3.00 genome with Fimo 67 , p-value < 1e − 5.

Carotenoid and chlorophyll analysis
Carotenoid and chlorophyll extraction was performed from aliquots of 20 mg dry mass, sequentially with 100 µL of a saturated solution of NaCl, 200 µL of dichloromethane and 500 µL of hexane:diethyl ether 1:1 (v/v). The pellet was extracted three additional times with 500 µL of 1:1 (v/v) hexane:diethyl ether. All supernatant fractions were combined, completely dried by vacuum and suspended in 200 µL of acetonitrile. Chlorophyll, phytoene, phyto uene, lycopene, β-carotene and lutein levels were determined via HPLC with a photodiode array detector 69 .

Competing interests
The authors declare no competing interests.  elements (TEs) and mC in all contexts (mCG, mCHG, mCHH) for the WT genotype. Global methylation changes for phyA and phyB1B2 in comparison with the wild type (WT) at the immature green (IG) and breaker (BK) stages are shown (bin size, 1 Mb). Gene and TE densities were estimated according to the number of nucleotides covered per million. The methylation levels in the CG, CHG and CHH contexts are 40-90%, 25-80% and 10-30%, respectively. The mC difference was relative to the corresponding WT fruit stage within a -5% (hypomethylated) range +5% (hypermethylated). (b) Number of genes with differentially methylated promoters (DMPs, 2 kb upstream transcription start site) in phyA, phyB1B2 and both mutants. Hyper-and hypomethylation are indicated by grey and darker-coloured bars, respectively. DMPs show statistically signi cant differences (FDR < 0.05) relative to WT. Boxplots show changes in the accumulation of cluster sRNAs in promoter (P, 2 Kb upstream of the 5' UTR end) and gene body (GB) regions for up-and downregulated DEGs. Asterisks indicate statistically signi cant differences by the Wilcoxon-Mann-Whitney test (** p<0.0001). All results represent the comparison of phyA and phyB1B2 to the wild type in immature green (IG) and breaker (BK) fruit stages.

Figure 4
Methylation across promoter and gene body regions differentially affects gene expression. (a) Relative expression of RIPENING INHIBITOR (RIN) and FRUITFULL 2 (FUL2) in immature green (IG) and mature green (MG) fruits from phyB1B2 determined by RT-qPCR. Red dots indicate data from RNA-seq in the same stage. Expression levels represent the mean of at least three biological replicates and are relative to the wild type (WT). Asterisks indicate statistically signi cant differences by the two-tailed Student´s t test compared to WT (* p<0.05). (b) Differential gene body methylation (green bars) and sRNA accumulation (black bars) within RIN and FUL2 in IG fruits from the phyB1B2 and WT genotypes.  genotypes. Green and orange indicate total mC in immature green (IG) and breaker (BK) fruits, respectively. HY5 and PIF transcription factor binding motifs are denoted with arrows. Thick blue lines indicate RIN binding sites according to ChIP-seq data27. (b) Relative expression from the RT-qPCR assay of genes encoding master ripening transcription factors in BK and red ripe (RR) fruits from phyB1B2. Red dots indicate data from RNA-seq in the same stage. Expression levels represent the mean of at least three biological replicates and are relative to WT. Asterisks indicate statistically signi cant differences by twotailed Student's t test compared to WT (* p<0.05).