The VIL gene CRAWLING ELEPHANT controls maturation and differentiation in tomato via polycomb silencing

VERNALIZATION INSENSITIVE 3-LIKE (VIL) proteins are PHD-finger proteins that recruit the repressor complex Polycomb Repressive Complex 2 (PRC2) to the promoters of target genes. Most known VIL targets are flowering repressor genes. Here, we show that the tomato VIL gene CRAWLING ELEPHANT (CREL) promotes differentiation throughout plant development by facilitating the trimethylation of Histone H3 on lysine 27 (H3K27me3). We identified the crel mutant in a screen for suppressors of the simple-leaf phenotype of entire (e), a mutant in the AUX/IAA gene ENTIRE/SlIAA9, involved in compound-leaf development in tomato. crel mutants have increased leaf complexity, and suppress the ectopic blade growth of e mutants. In addition, crel mutants are late flowering, and have delayed and aberrant stem, root and flower development. Consistent with a role for CREL in recruiting PRC2, crel mutants show drastically reduced H3K27me3 enrichment at approximately half of the 14,789 sites enriched in wild-type plants, along with upregulation of many underlying genes. Interestingly, this reduction in H3K27me3 across the genome in crel is also associated with gains in H3K27me3 at a smaller number of sites that normally have modest levels of the mark in wild-type plants, suggesting that PRC2 activity is no longer limiting in the absence of CREL. Our results uncover a wide role for CREL in plant and organ differentiation in tomato and suggest that CREL is required for targeting PRC2 activity to, and thus silencing, a specific subset of polycomb targets.


Introduction
Polycomb Repressive Complex 2 (PRC2) is a conserved complex that represses gene expression by trimethylating lysine 27 of histone H3 proteins (H3K27me3) [1][2][3]. PRC2 activity counteracts, and is counteracted by, the transcription-promoting functions of trithorax-group proteins [4]. The core PRC2 is composed of 4 subunits. In plants, some of these subunits are encoded by small gene families, allowing the formation of multiple, distinct complexes. Different plant PRC2 complexes have been shown to regulate specific developmental processes such as endosperm development, flowering time and flower development [2,3]. As PRC2 complexes do not have DNA binding domains, they are recruited to target loci by interacting proteins [2,[5][6][7][8][9][10]. One of the most characterized PRC2-regulated processes in Arabidopsis is the induction of flowering in response to prolonged cold, termed vernalization. In response to vernalization, PRC2 promotes flowering by silencing the flowering inhibitor FLC. The vernalizationspecific VRN-PRC2 complex is recruited to FLC by complexing with PHD proteins from the VERNALIZATION INSENSITIVE 3-LIKE (VIL) family [7,[11][12][13]. In Arabidopsis, the VIL family consists of 4 members, including VIN3 and VRN5. Vernalization induces VIN3 expression, while VRN5 is expressed constitutively. VIL proteins also repress additional members of the FLC family during vernalization, and VRN5 and VIL2 are also involved in other flowering pathways [7,11,14,15].
VIL proteins have been identified from several plant species [16][17][18][19][20][21][22][23][24][25]. They have been shown to promote flowering in all tested species, including species that do not have an FLC ortholog and/or do not respond to vernalization. In rice, the OsLF and OsLFL1 genes encode transcription factors that inhibit flowering and have been identified as VIL targets [18,20]. A VIN3 ortholog has also been identified in tomato [24]. While the vast majority of research on VIL proteins concerned their involvement in flowering induction, several studies reported additional developmental effects. For example, Arabidopsis vrn5 mutants had increased leaf curling, increased numbers of petals, and distorted siliques [13]. In rice, leaf inclination2 (lc2) mutants had an altered leaf angle, curled leaves and severe sterility, and other rice VIL genes were found to affect spikelet development, branching and grain yield [16,17,19,26]. Silencing the Brachypodium distachyon BdVIL4, which is similar to VIN3, led to increased branching [21]. Pepper cavil1 mutants affect leaf development, apical dominance and branching [25]. However, the knowledge about the involvement of VIL proteins in these and other developmental processes is limited, and their role in compound-leaf development has not been explored. In addition, it is not clear whether VIL proteins recruit PRC2 mainly to targets involved in the induction to flowering or whether they have broader roles in plant development.
Tomato plants have compound leaves, which are composed of multiple leaflets [27]. The elaboration of compound leaves depends on slow maturation of the developing leaf, which enables an extended organogenesis activity at the leaf margin, during which leaflets are formed [28][29][30][31][32]. Leaflets are formed by differential growth at the leaf margin, where regions of blade growth are separated by intercalary regions of growth inhibition [33]. Auxin has been shown to promote growth and its response is inhibited in the intercalary domains [34][35][36][37][38][39][40]. Mutations in the tomato gene SlIAA9/ENTIRE (E), encoding an auxin-response inhibitor from the Aux/ IAA family that specifies the intercalary domain, result in simplified leaves due to ectopic blade growth in the intercalary domain [34,41,42].
Here, a screen for suppressors of the e simple-leaf phenotype identified the crawling elephant (crel) mutant, which substantially suppresses the ectopic blade growth of e. We found CREL to encode a tomato VIL gene, related to Arabidopsis VIL1/VRN5. crel mutants affect many aspects of tomato development, including plant and organ maturation. Comparison of H3K27me3 modifications between wild type and crel plants showed that CREL is required for H3K27me3 deposition at approximately half of the H3K27me3-enriched sites found in wildtype plants. Therefore, CREL promotes maturation throughout the plant life by promoting selective deposition of H3K27me3 and gene silencing at a subset of PRC2 targets.

crawling elephant (crel) mutants suppress entire (e) and have very compound leaves
entire (e) mutants, mutated in a tomato Aux/IAA gene, have simplified leaves in comparison to the wild-type compound leaves [34,40,42,43] (Fig 1A and 1B). To identify genes that are involved in compound-leaf development, we generated an Ethyl Methane Sulfonate (EMS) mutant population in the background of e, and screened for suppressors of the e simplified leaf phenotype. This screen identified the crawling elephant-1 (crel-1) mutant as a strong e suppressor. e crel-1 double mutants had distinct, clearly separate primary leaflets, and occasionally had secondary leaflets, in contrast to the mostly entire leaf shape of single e mutants (Fig 1A-1C). To characterize the unique crel-1 phenotype, we backcrossed crel-1 to the parental line (Solanum lycopersicum M82), and identified single crel-1 F2 individuals ( Fig 1D). Leaves of single crel-1 mutants were much more compound than wild-type leaves, with a similar number of primary leaflets but many more secondary leaflets than the wild type. In contrast to wild typeleaves, crel-1 leaves also had tertiary leaflets (Fig 1A, 1D, 1J and 1K). Therefore, crel mutants suppress the e simplified-leaf phenotype, and form many more leaflets in both the wild type and the e backgrounds. Interestingly, previously identified e suppressors such as slmp and slarf19a,b had a reduced number of leaflets in the respective single mutants [39].
We identified several additional crel alleles from the Menda EMS and fast neutron mutant population [44], in the M82 background, and confirmed allelism by complementation tests. The alleles showed a range of phenotypic severities, including a diverse increase in leaflet number (Fig 1E-1H, 1J and 1K). Similar to crel-1, crel-2 also suppressed the e simplified leaf phenotype (S1A-S1D Fig).

CREL acts during relatively late stages of leaf development
To investigate the timing of the effect of crel mutants on leaf development, we compared leaf development between wild type and crel-1 plants. Early stages of leaf development were very similar between wild type and crel-1 plants when similar developmental stages were compared, although the terminal leaflet expanded earlier in crel-1 ( Fig 1I). However, the rate of leaf and leaflet initiation was much slower in crel-1 mutants than in the wild type (S1E Fig). At later stages of leaf development, when wild type leaves stopped generating new leaflets, crel-1 and crel-2 leaves continued to form leaflets more than a month later ( Fig 1K). Therefore, crel leaves develop slower than the wild type, and while the terminal leaflet appears to differentiate early, overall leaf differentiation is substantially delayed in crel mutants.

PLOS GENETICS
To further characterize this effect of crel on leaf development, and understand the timing and developmental context of CREL action, we analyzed its genetic interaction with mutants that affect the developmental window of leaflet morphogenesis. Leaflets are formed during the morphogenesis stage of leaf development, which follows leaf initiation and precedes leaf expansion and differentiation [27,28,45,46]. The elaboration of compound leaves depends on an extended morphogenesis stage. The CIN-TCP transcription factor LANCEOLATE (LA) and the MYB transcription factor CLAUSA (CLAU) promote maturation and differentiation and thereby restrict the morphogenetic window [29,[47][48][49][50]. La-2 is a semi-dominant mutant in which LA is expressed precociously due to a mutation in the miR319 binding site. This accelerates leaf differentiation and leads to a simple leaf form (Fig 2A and 2C). La-2 was epistatic to crel-1 (Fig 2B-2D), indicating that the morphogenetic window in La-2 is terminated before the timing of CREL action, in agreement with the relatively early effect of LA and late effect of CREL on leaf development. Leaves of loss-of-function clau mutants have an extended morphogenetic window, leading to a substantial increase in leaf complexity and leaflet number ( Fig 2E) [48,51]. crel-1 clau double mutants had very complex leaves (Fig 2E and 2F), suggesting that CREL acts in parallel with CLAU to restrict leaf elaboration and promote maturation. Removing the activities of both regulators leads to prolonged, extensive leaflet morphogenesis. Together, these results suggest that CREL acts in relatively late stages of leaf development to promote maturation and differentiation.  [52]; EdII-GUS is a stabilized form of E fused to the GUS reporter; ARF10m is a mutant form of ARF10 that is mutated in the miR160 binding site. A-D-epistasis of La-2/+, in which differentiation is accelerated, over crel. E, F-enhancment of crel by clau, in which differentiation is delayed. G-J-enhancment of genotypes with narrow leaves due to reduced auxin response by crel. K, L-suppression of the ectopic blade growth of

PLOS GENETICS
The suppression of the e phenotype by crel raised the question of whether CREL is also involved in the differential growth at the leaf margin that leads to the formation of separate leaflets. To address this question, we crossed crel to mutants affected in auxin-mediated blade growth. Ectopic expression of a stabilized form of E resulting from a mutation in domain II of the E (IAA9) protein (EdII) resulted in leaflet narrowing [42], (Fig 2G). This effect was strongly enhanced in the crel-1 background (Fig 2H), suggesting that E inhibits and CREL promotes blade expansion, but they act in at least partially parallel pathways. Similarly, crel-1 enhanced the narrow blade phenotype resulting from ectopic expression of a miR160-resistant ARF10, a negative regulator of blade expansion (FIL>>ARF10m, Fig 2I and 2J). In agreement, crel-1 suppressed the ectopic blade growth resulting from ectopic expression of miR160, which negatively regulates ARF10 and 4 additional ARF proteins (FIL>>miR160, Fig 2K and 2L) [40]. The suppression of FIL>>miR160 was more prominent in later leaves than in early leaves (S2 Fig). Interestingly, leaflet number was reduced in FIL>>ARF10 crel-1 relative to both single mutants, and crel-1 FIL>>EdII-GUS plants were extremely small with almost no leaflets, suggesting that extreme repression of lamina growth leads to a reduction in leaflet formation and overall growth. Together, these results suggest that CREL promotes blade growth during compound-leaf development, and acts at least partially in parallel to auxin.

CREL is a VRN5 homolog
To identify the CREL gene, we genetically mapped the crel-1 mutation using an F2 mapping population from a cross between the crel-1 mutant, in the Solanum lycopersicum M82 background, and S. pimpinellifolium. crel-1 was mapped to chromosome 5. Further mapping was hampered by an introgression of S. pimpinellifolium sequences in the M82 line in this region [53]. We therefore used RNA-seq to identify possible causative mutations in crel-1 and crel-2, which led to the identification of mutations in the gene Solyc05g018390 in both crel-1 and crel-2. In crel-1, a G to A substitution at position 4264 from the transcription start site (TSS) led to a stop codon in exon III. The fast neutron allele crel-2 contains a 12,826-bp-long deletion, which results in the elimination of exon I and II and part of exon III (Fig 3A). Sequencing the Solyc05g018390 gene in two additional crel alleles identified a 1-bp deletion in the first exon at position 322 from the TSS in crel-3, and an A to T substitution in position 3630 leading to a stop codon in the third exon in crel-5 ( Fig 3A). We therefore concluded that Solyc05g018390 is CREL. CREL is predicted to encode a plant homeodomain (PHD) finger protein (Fig 3A and  3B). It is most similar to the Arabidopsis VRN5 gene.

CREL is expressed in expanding blades
We characterized the expression of CREL in the fifth leaf produced by the plant, at different developmental stages, to examine how its expression correlates with its activity. CREL was expressed throughout leaf development, with relatively low expression in apices containing the SAM and very young P1-P3 primordia. Later, its expression was gradually upregulated, peaking at P6/P7 (Fig 3C). To spatially localize CREL in leaf primordia, we cloned a 2960-long CREL promoter and used it to generate a CREL driver line in the transactivation system [52,56]. In developing leaves, the CREL promoter drove expression in leaf margins. In agreement with the qPCR experiment, expression appeared lower in young primordia, and increased from P4 on. Expression was mainly visible in expanding regions of the leaf margin, the terminal leaflet at P4, and the expanding leaflets at P6 and on (Fig 3D-3G). The expression appeared to follow the basipetal differentiation wave of the leaf, with strong expression first appearing in the terminal leaflet, which is the first to expand and differentiate, and then progressing basipetally in expanding leaflets. This leaf expression pattern is compatible with the crel leaf phenotype, which starts to differ from the wild type around the P5 stage (Fig 1I), and with the genetic interactions showing that crel affects leaf maturation and blade formation (Fig 2).

CREL promotes multiple aspects of plant maturation and differentiation
In addition to their effect on leaf differentiation and patterning, additional differentiation processes were also delayed and/or impaired in crel mutants. crel plants failed to maintain an

PLOS GENETICS
upright position and the plants exhibited a sprawling growth habit. Similar to the effect on leaf shape, this phenotype developed at a relatively late stage of plant development (Fig 4A and  4B). To understand the basis for the "crawling" phenotype, we dissected developing wild-type and crel-2 stems at successive developmental stages. We sectioned the internode between the cotyledons and the first leaf from different plants grown together, between the ages of 3 and 10 weeks. In three-week-old plants, crel-2 stems were narrower than the wild type with nearly normal although slightly less developed vascular bundles ( Fig 4C, 4D, 4I and 4J). crel-2 vasculature continued to develop slower than the wild type, and ceased maturation and differentiation prematurely. This resulted in a thin and undeveloped xylem in crel-2 stems. Specifically, crel-2 stems failed to complete a vascular cylinder, and had only partial secondary xylem development ( Fig 4E-4H, 4K-4N). The reduction in supporting tissue likely contributes to the reduced strength of the crel stem.
Root vasculature development was also delayed and impaired in crel mutants. crel-2 roots were narrower than wild-type roots, and their vascular tissue developed slowly and failed to reach full differentiation (Fig 5A-5F). To investigate the effect of crel on the root system as a whole, wild type and crel-2 plants were grown hydroponically, and root volume and length were calculated at successive times. Root volume and length were reduced in crel-2 plants, and the difference increased with time, although the difference was statistically significant at one of the time points only (Fig 5G and 5H). Therefore, CREL plays an important role in root development and differentiation.
The crel mutation also affected flowering time and flower development. crel-1 and crel-2 mutants flowered much later than the wild type, after producing 12-13 leaves, compared to 6

PLOS GENETICS
leaves in the wild type ( Fig 6A). Mature crel flowers were not fully developed, had short and distorted organs and were sterile (Fig 6B and 6C). Early flower development was similar between wild type and crel plants, except for the sepals that were curled backwards in crel, resulting in an open bud where the inner organs were not covered by the sepals. However, crel flower organs ceased development and growth prematurely (Fig 6D-6G). Therefore, crel mutants were delayed in multiple developmental pathways. In some cases such as flower, stem, and root development, these organs failed to properly differentiate, while in others, such as leaf development and flowering time, they differentiated substantially slower than the wild type. Overall, crel plants had aberrant plant and organ structure, which resulted in weak and sterile plants.
In Arabidopsis, VIL proteins affect flowering time by silencing the expression of the flowering repressor FLC. Three tomato genes were included in the FLC/MAF clade in a phylogenetic analysis of tomato MADS-box genes, SlMBP8 (Solyc12g087830), SlMBP15 (Solyc12g087820) Therefore, SlMBP8 may be involved in the late flowering phenotype of crel, but further analysis is required to assess its role.

PLOS GENETICS
this prediction, we performed ChIP-seq for the H3K27me3 modification in shoot apices of 4-week-old wild-type and crel-2 plants (S2 Table). In wild-type tissue, H3K27me3 was found to be enriched in a characteristic pattern over gene bodies, with crel-2 mutants showing a decrease in average H3K27me3 over genes (Fig 7A). To determine whether this decrease was due to a small loss at all H3K27me3 sites or a large loss at a smaller number of sites, we identified all enriched sites found in two replicates for each genotype and compared these reproducibly enriched sites between wild-type and crel-2. Of the 14,789 reproducibly enriched regions in WT, 7,469 of these sites were essentially depleted of H3K27me3 in crel-2 mutant plants ( Fig  7B and 7D; crel-lost H3K27me3 sites). Interestingly, 1,295 sites showed increases in H3K27me3 in the crel-2 mutant (Fig 7C and 7E; crel-gained H3K27me3 sites). The vast majority of these sites are normally enriched for H3K27me3 in wild-type, suggesting that in the

PLOS GENETICS
absence of CREL, excess PRC2 activity is directed to these sites that were only modestly enriched with H3K27me3 in wild-type. The crel-lost and crel-gained H3K27me3 site coordinates as well as associated genes are described in S3 Table. To examine the effects of H3K27me3 loss on the transcriptome, we performed RNA-seq in multiple replicates on shoot apices of 4-week-old wild-type plants and plants with 3 independent crel alleles (Fig 8). RNA was prepared from a total of seven biological replicates of crel mutant

PLOS GENETICS
plants and compared to three biological replicates of wild-type S. lycopersicum. Differential gene expression analysis revealed that 431 genes were significantly upregulated in all crel mutants and 265 were significantly downregulated compared to the wild-type (Fig 8A and S4 Table). Gene ontology analysis utilizing agriGO [59] showed significant enrichment among these upregulated genes for the biological function "Transcription Factor Activity" (FDR = 8.1e-6) and other terms associated with transcriptional regulation (S5 Table). The genes downregulated in crel mutants showed no significantly enriched GO terms using the same parameters.
To understand the relationship between chromatin alterations and transcriptional changes in crel, we mapped the crel-lost and crel-gained H3K27me3 sites to genes. The 7,469 crel-lost sites were associated with 4,242 genes, 133 of which were also upregulated in crel mutants (of 431 total upregulated genes) (Fig 8B and S6 Table). This reinforces that loss of CREL results in changes in gene expression by affecting H3K27me3 deposition at a subset of genes. While not all genes associated with crel-lost H3K27me3 peaks showed significant gene expression changes under the stringent expression criteria used, when statistical significance is not considered the majority of the genes that lose H3K27me3 in crel have log2 fold changes greater than 0 (S5 Fig). As many aspects of plant and organ maturation and differentiation are delayed and/or impaired in crel mutants, we asked whether the molecular signature of maturation and differentiation is affected in crel.   [32] defined sets of "morphogenesis" and "differentiation" genes based on transcriptomic data of successive stages of leaf development as well as mutants that affect morphogenesis and differentiation. We examined whether there is an enrichement in "differentiation" and "morphogenesis" signatures among genes that showed altered expression and/or altered H3K27me3 modifications in crel. This analysis showed that "differentiation" genes were significantly enriched among the genes that were downregulated in crel, with 4 times more "differentiation" genes then expected by random overlaps between groups being present in the subset of genes that were significantly down-regulated in crel (p-val < 8.509e-15, S7 Table).

Discussion
VIL proteins have been shown to affect flowering in several plant species, by repressing the expression of flowering repressors, such as FLC in Arabidopsis. In addition to their effect on flowering, VIL genes were found to affect an array of developmental processes in different species. This work identifies the tomato VIL gene CREL as a mediator of diverse developmental processes, via the modulation of H3K27me3 modifications in many genes throughout the tomato genome, likely by recruiting PRC2 complexes to a subset of their target genes.

CREL promotes plant and organ maturation
So far, VIL genes have been mainly implicated in flowering time [7,14,16,18,[20][21][22]25,60]. Here, we uncover a much broader role for this gene family in plant development, as revealed from the phenotypes and the effect on H3K27me3 modification. crel mutants are affected in many aspects of plant maturation and differentiation in addition to the delay in flowering time. crel mutants flower late and have delayed leaf maturation, resulting in an extended leaf morphogenesis and more compound leaves. Interestingly, while flowering and leaf maturation eventually occur in crel, stem, root, and flower differentiation are impaired in crel and these organs do not reach full differentiation and function. CREL accumulates relatively late during leaf development, thus enabling prolonged morphogenesis. Recently, a growth-rate dependent mechanism of controlling VIN3 accumulation in the cold has been described [61]. It would be interesting to understand the mechanism by which CREL expression is delayed during organ maturation to enable timed maturation and differentiation.
Other genes involved in the induction of flowering were also shown to affect maturation and differentiation in additional developmental aspects. For example, the tomato flowering inducer SFT, the ortholog of FT, was shown to promote leaf maturation and affect stem differentiation [62]. Recently, SFT was shown to specifically affect secondary cell wall biosynthesis in tomato stems [63]. FT was also proposed to promote maturation and termination in additional species [64,65]. The pepper cavil1 mutants, impaired in the pepper CREL ortholog, have reduced vascular development but wider stems [25]. Therefore, both precocious and delayed differentiation impairs the final form and function of stems. In Cardamine hirsuta, plant maturation and flowering was shown to be coordinated with age-dependent changes in leaf shape in plants with variable FLC activity [66]. In contrast to CREL and SFT, which promoted all aspects of plant and organ maturation, tomato CIN-TCPs were shown to promote leaf maturation but delay plant maturation, while AP1/FUL MADs BOX genes promoted plant maturation and delayed leaf maturation [67,68]. Overall, similar to CREL, genes that have been implicated mainly in the promotion of flowering in Arabidopsis were found to promote a wide range of differentiation and maturation aspects.
The involvement of CREL in plant and organ differentiation is in agreement with a role in mediating PRC2 activity. A common function of PRC2 genes in plants is the maintenance of a differentiated state, and prc2 mutants in both mosses and seed plants have phenotypes related to dedifferentiation and overproliferation [1,2,69]. Therefore, CREL may aid in recruiting PRC2 to differentiation-related target genes.
VIL genes from other species have also been shown to affect other developmental processes in addition to flowering. [13,16,19,21,25,26]. Interestingly, beside the common effect on flowering time, the specific developmental effects only very partially overlap among these species. This suggests that the VIL family may be used as a tool for developmental innovations, recruiting an existing tool to different processes. Specifically interesting in this respect is the comparison between tomato and pepper, which are closely related species that differ in several key developmental features, such as flowering architecture and leaf shape. cavil1 mutants have reduced vasculature development, increased plant and organ size, increased branching and reduced angle of axillary branches [25]. Interestingly, only some of these additional phenotypes overlap with crel. In contrast to the simple leaves of pepper, tomato leaves are compound, with several orders of leaflets. This is correlated with faster differentiation of the pepper leaf, similar to tomato La-2/+ mutants [29]. The current work revealed an important role for CREL in the development of the compound leaf, with an effect on both the rate of differentiation and leaf patterning (Figs 1 and 2), further supporting the notion that VIL genes have been recruited to diverse, partially species-specific processes.
We propose that, in addition to its general effect on growth and differentiation, CREL also promotes blade growth in developing tomato leaves. crel mutants suppress the ectopic blade outgrowth of e mutants and miR160-overexpressing plants (Figs 2K and 2L, S2). Furthermore, crel enhances the narrow-leaf phenotype caused by overexpression of E or miR160-targeted ARFs. In addition, CREL expression is elevated in later stages of leaf development when the blade begins to expand, and the CREL promoter shows high expression in growing regions of the leaf margin (Fig 3D-3G). The genetic interactions between crel and auxin-related mutants suggest these auxin mediators and CREL act via independent pathways to regulate blade growth. Therefore, CREL likely promotes blade growth either downstream of auxin or through an at least partially parallel pathway. As most effectors of compound-leaf development have been shown to affect either the organ-level differentiation rate (for example LA and CLAU) or local differential growth (E, CUC, SlMP) [27], it is interesting that CREL appears to affect both aspects.

CREL affects H3K27me3 modifications throughout the tomato genome
Only a handful of VIL targets have been identified so far, most of which are related to their role in promoting flowering. In Arabidopsis, FLC and FLC-related proteins are targeted by different VIL proteins in specific flowering pathways [11,14]. In rice, the flowering inhibitors OsLF and OsLFL1 were identified as VIL targets [18,22,70]. In addition, a cytokinin catabolism gene from the CKX family and the bud-outgrowth inhibitor OsTB1 have been identified as a VIL target in rice [19,26]. The microRNA miR156 was proposed as a target of BdVIL4 in Brachypodium [21], and several putative targets have been proposed to mediate the effect on flowering of pepper VIL1 [25]. Here we found that of the 14,789 sites enriched with H3K27me3 in wild-type plants, approximately half (7,469) lose nearly all H3K27me3 in crel mutants. The identification of such a global effect on H3K27me3 modification in crel mutants indicates that this VIL protein has a widespread role in polycomb silencing. Together with the pleiotropic phenotypic effect, this suggests that VIL proteins are involved in a wide range of developmental processes, and play a central role in recruiting PRC2 complexes to many targets genome wide. The similarly pleiotropic effect of Cavil1 mutants in pepper, together with its effect on gene expression [25], suggests that this is also true in other species.
Interestingly, the widepsread reductions in H3K27me3 across the genome in crel are also associated with gains in H3K27me3 at a smaller number of sites (1,295) that normally have modest levels of the H3K27me3 mark in wild-type plants. This phenomenon of H3K27me3 redistribution in PRC2-targeting mutants has been previously observed in both plants and animals [8,71] and suggests that sufficient PRC2 activity is a limiting factor in H3K27me3 accumulation under normal circumstances. However, we are aware that, prior to mapping the CREL binding sites, we cannot rule out that CREL plays a role in targeting PRC2 activity to many loci while directly inhibiting it at others. This distinction awaits future experimental testing.

A conserved role for VIL genes in promoting flowering
VIL genes have been shown to affect flowering time in many species where mutations or silencing of these genes have been described, including both dicots and monocots [7,14,16,18,20,22,25]. In Arabidopsis, VIL proteins promote flowering by recruiting PRC2 to flowering repressors from the FLC family, thus facilitating the deposition of the repressive chromatin modification H3K27me3. Different Arabidopsis VILs act to induce flowering in specific combinations of flowering pathways, timing and target genes [11,14]. Interestingly, while tomato plants do not require vernalization for flowering, crel mutants are late flowering. Similarly, VIL genes promote flowering in additional species that do not require vernalization for flowering, or in species with a different vernalization mechanism. Rice VIL proteins including LC2 were shown to act by repressing OSLFL1 and OsLF, two flowering repressors unrelated to FLC [18,22,70]. Therefore, while the effect on flowering and possibly the molecular mechanism are conserved, the target genes differ among species [25]. While the FLC homolog SlMBP8 (Solyc12g087830.1) shows reduced H3K27me3 enrichment in crel-2 (S3 Fig), further research is required to understand whether it mediates the crel late-flowering phenotype.

Plant material and growth conditions
Tomato plants (Solanum lycopersicum cv M82) were germinated and grown in a controlled growth room or in a commercial nursery for four weeks. Then the seedlings were transferred to a greenhouse with natural daylight and 25˚C/ 20˚C day/night temperature, or to an open field with natural daylight and temperature. crel-1 was isolated in this work by a mutant screen in the e-3 background (Berger 2009, Ben Gera 2012), as described below. crel-2-crel-5 are from the mutant populations described by Menda et al., [44]. The transactivation system, described previously [52,56], was used to characterize the CREL promoter and for leaf-specific expression. This system consists of driver lines and responder lines. In the driver lines, specific promoters drive the expression of the synthetic transcription factor LhG4, which does not recognize endogenous plant promoters. In the responder lines, a gene of interest or a reporter is expressed downstream of the E.coli operator, recognized by LhG4 but not endogenous plant transcription factors. A cross between a driver and a responder line results in a plant (designated PROMOTER>>GENE) expressing the gene of interest/marker under the control of the specific promoter. La-2, clau, the FIL driver line and the ARF10, miR160 and EdII-GUS responder lines have been previously described [35,40,47,51,52,72].

Generation of a mutant population, screening and identification of crel-1
Around 750 entire-4 (e-4) seeds were treated with the mutagenic substance Ethyl-Methane Sulfonate (EMS, Sigma m0880) at a concentration of 0.6% for 10 hours. Around 50 seeds underwent a control treatment without exposure to EMS. The treated seeds were sown in a commercial nursery and the seedlings (M1 generation) were transferred to a greenhouse. M1 plants were self-pollinated to increase the number of seeds per plant. M2 seeds were collected separately from each of around 650 M1 plants. M2 progeny (around 40 seeds per family) were grown in an open field, and screened frequently during the season for mutants that affect the development of the leaf, flower and fruit. e crel-1 was identified in this screen as a recessive mutant segregating 1:3 in an M1 family. Single crel-1 mutants were generated by a cross between e crel-1 and wild type and identification of single crel-1 mutants in the F2 generation. crel-1 was then back-crossed three times to M82 for further characterization.

Allelism tests and genetic interactions
As crel mutants are sterile, allelism tests were performed by crossing heterozygote siblings. Progeny of a cross between two heterozygous alleles segregated ¼ mutant progeny. Similarly, genetic interactions between crel and other mutants or transgenic lines were generated by crossing heterozygous crel siblings with the respective mutant or transgenic line.

Identification of the CREL gene
An F2 mapping population was generated by crossing crel-1/+ plant, in the M82 background to Solanum pimplinellifolium, and collection of F2 progeny from individual F1 plants. Initial mapping with 30 F2 plants showing the crel-1 phenotype, and 50 mapping markers developed by Revital Bronstein, Yuval Eshed (Weizmann Institute) and Zach Lippman (CSHL) and spread along the tomato genome, identified linkage to 3 markers on chromosome 5. Fine mapping of 120 crel-1 F2 individuals and additional markers located the gene to a region between markers zach 43.2 and jose 58.1 dcap, located between bases 43,123,344 and 58,170,500 on chromosome 5. Further mapping was hampered by an introgression of S. pimpinellifolium sequences in the M82 line in this region [53] We therefore used RNAseq of wild type, crel-1 and crel-2 plants to identify polymorphism in these alleles.
For RNAseq, shoot apices containing the SAM and the 4 youngest primordia were collected from 14-day-old M82, crel-1 and crel-2 plants, in which L4 (the 4 th leaf produced by the plant) was at the P4 stage. RNA was extracted using the RNeasy micro kit (Qiagen), using the manufacturer's instructions. Two biological replicates were used for M82 and crel-2, and one biological replicate for crel-1. Sequencing libraries were prepared according to the Illumina TruSeq RNA protocol and sequenced on an Illumina HiSeq2000 platform at the Genome Center of the Max Planck Institute for Plant Breeding Research. We obtained between 21,3 and 28,3 million 96-bp single-end reads per library (average of 25,8 million). Reads were aligned to the S. lycopersicum reference sequence v2.40 using TopHat v2.0.6 [73] with the following parameters: -max-insertion-length 12 -max-deletion-length 12 -g 1 -read-gap-length 12 -read-editdist 20 -read-mismatches 12 -no-coverage-search -read-realign-edit-dist 0 -segment-mismatches 3 -splice-mismatches 1. To detect polymorphisms between the crel mutants and wildtype M82, biological replicates from each genotype were merged. Then, duplicated reads were removed using default settings in Picard (http://broadinstitute.github.io/picard/), indels were realigned using GATK v2.2-8, and variants called in all samples simultaneously using default parameters in GATK v2.2-8 [74]. Next, we estimated the effect of each variant in annotated transcripts (ITAG 2.3) using ANNOVAR [75]. Variants in the candidate region in chromosome 5 determined by QTL analysis were evaluated manually.

Phenotyping and imaging
Characterization of early leaf development and rate of leaf initiation. Plants were sown, germinated and grown in a growth chamber. Every two weeks, the number of leaves and leaf primordia were counted from six plants from each genotype. The fifth leaf (L5) was photographed by a stereoscope (Nikon SMZ1270) and its developmental stage determined. Six different plants were used for each time point.
Quantification of leaf complexity. Leaves 5, 7 and 9 were marked at the time of their emergence from the shoot apex. Then, the number of leaflets was counted every 7-14 days for each marked leaf. Primary, intercalary, secondary and tertiary leaflets were counted. At least 9 plants were counted for each genotype.
SEM characterization of flower development. Flowers from different developmental stages of each genotype were collected and their petals removed using a stereoscope, placed on a microscope stub with a carbon strip and analyzed with Hitachi TM3030 Plus SEM.
Phenotypic quantification and statistical analysis. For the quantification of the number of leaves to flowering, plants were grown in a greenhouse, and with the appearance of the first flower, the number of leaves formed before the flower were counted. At least 9 biological repeats, each consisting of a single plant, were quantified. The experiment was repeated twice, once with plants germinated in a commercial nursery and once with plants germinated in a growth chamber. For the quantification of root phenotypes, seedlings were grown hydroponically in Hoagland nutrient solution (pH 6.5), in a growth room set to a photoperiod of 12/12-h night/days, light intensity (cool-white bulbs) of *250 μmol m −2 s −1 , and 25˚C. After 28, 34, 38 and 43 DAG the roots of 3 plants of each genotype were scanned and analyzed using a flatbed scanner (Epson 12000XL, Seiko Epson, Japan) and root architecture was analysed using WinRhizo software (Regent Instruments Ltd., Ontario,Canada). Statistical analysis was performed using the JMP software (SAS Institute, http://www.sas.com). Means and p values were calculated using the Student's t-test or the Dunnett's test, as indicated in the figures.
Histological characterization of stem tissues. Ten to 50 day-old plants were free-hand dissected using a double-sided razor blade. 1-2-mm-long sections were dissected from up to 5 cm below the node. Sections were dehydrated in acetic acid: ethanol [1:10] for 1 hour and then stained directly with TBO (0.01% aqueous, sigma). Images of early developmental stages were captured using Nikon a SMZ1270 stereoscope equipped with a Nikon DS-Ri2 camera and NIS-ELEMENTS software.
Confocal imaging. For analysis of pCREL:nYFP expression, dissected whole-leaf primordia were placed into drops of water on glass microscope slides and covered with cover slips. The pattern of YFP expression was detected by a confocal laser scanning microscope (CLSMmodel SP8; Leica), with the solid-state laser set at 514 nm for excitation and 530 nm for emission. Chlorophyll expression was detected at 488nm for excitation and 700nm for emission. ImageJ software was used for analysis, quantification, and measurements of captured images. Images were adjusted uniformly using Adobe Photoshop CS6. Tomato stems and primary roots were cut to sections of 200um and 300um width respectively using Leica VT1000 vibratome and were and cleared using ClearSee [76], cell wall staining was performed using SR2200 (Renaissance Chemicals) prior to mounting and visualization using 405nm laser.

Cloning and plant transformation
The CREL promoter was generated by amplifying 3000 bp upstream of the CREL ATG from genomic DNA and cloned upstream to LhG4 generating the pCREL:LhG4 driver line.
Plant transformation and tissue culture were performed as described in Israeli et al 2019 [39]. At least five independent kanamycin-resistant transgenic lines from each transgene were genotyped and, in the case of pCREL:LhG4, crossed to an OP:YFP stable line to generate pCREL>>YFP. Three lines from each transgene or resultant cross were examined, and a representative line was selected for further analysis.

Phylogenetic analysis
Phylogenetic analysis was performed using full-length protein sequences of the tomato, Arabidopsis, rice and pepper VIL gene family or the FLC/MAF MADS-BOX clade. The sequences were obtained from the Sol Genomics Network (SOL, https://solgenomics.net/), The Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org/) and the Plaza tool (https:// bioinformatics.psb.ugent.be/plaza/). Sequences were aligned and phylogenetic trees were constructed using MEGA X [54,55], using using a Maximum Likelihood method, branch lengths represent the expected number of substitution per site.

RNA extraction and qRT-PCR analysis
RNA was extracted using the Plant/Fungi Total RNA Purification Kit (Norgen Biotek, Thorold, ON, Canada) according to the manufacturer's instructions, including DNase treatment. cDNA synthesis was performed using the Verso cDNA Kit (Thermo Scientific, Waltham, MA, USA) using 1 μg of RNA. qRT-PCR analysis was carried out using a Corbett Rotor-Gene 6000 real-time PCR machine, with SYBR Premix for all other genes. Levels of mRNA were calculated relative to EXPRESSED (EXP) [77] or TUBULIN (TUB) [78] as described [29]. Primers used for the qRT-PCR analysis are detailed in S1 Table. ChIP-seq procedures WT and crel plants were grown on soil under 16 hrs of light/8 hrs dark cycles for 28 days after germination. Shoot apices, (0.8 g for each replicate and two replicates per genotype) containing approximately three visible expanding leaves, were harvested and fixed in 1% formaldehyde + 0.2% Silwet L-77 for 17 minutes under vacuum. Glycine was then added to a final concentration of 0.125 M and tissue was placed under vacuum for an additional 5 minutes, followed by washing several times in water. ChIP was performed on the fixed tissue using the procedure of Gendrel et al. [79]. For each ChIP reaction, we used 2 ug of a rabbit polyclonal antibody against H3K27me3 (Millipore, catalog #07-449). Input and ChIP DNAs were converted to Illumina sequencing libraries using the Accel-NGS 2S Plus DNA library kit according to the manufacturer's instructions (Swift Biosciences). Libraries were sequenced on an Illumina NextSeq 500 instrument using 50-nt single end reads at the University of Georgia Genomics and Bioinformatics Facility.

ChIP-seq data processing
Raw reads were mapped to the SL3.0 build of the tomato genome using Bowtie2 [80] with default parameters. Raw mapped reads were then processed using Samtools [81] to retain only those with a mapping quality score greater than or equal to 2 and PCR duplicates were removed using Samtools "rmdup" function. Enriched regions (peaks) for H3K27me3 were then identified for each replicate using the "Findpeaks" function of the HOMER package [82]. Further analyses only considered peaks that were identified in both replicates for each genotype and reads were scaled to the lowest read depth using Samtools "view -s" function. Unique WT and crel sites were detected by using the Bedtools "subtract" function with the -wa parameter, where only peaks that were completely lost in either WT or crel-2 were retained. These unique peaks from WT and crel-2 were mapped to their nearest TSS using the "annotatePeaks. pl" function from the HOMER package [77]. See S2A Table for more information on read count numbers after each of these steps. For normalization and visualization, quality-filtered reads in.bam format were converted to bigwig format using the "bamcoverage" script in deep-Tools 2.0 [83] with a bin size of 1 bp and RPKM normalization. Heat maps and average plots displaying ChIP-seq data were also generated using SeqPlots visualization software [doi: 10. 12688/wellcomeopenres.10004.1].

RNA-seq library preparation and sequencing
Wild-type and crel mutants shoot apices were grown and harvested as for ChIP-seq, and RNA was prepared using the Qiagen RNeasy Mini Plant kit according to the manufacturer's instructions. This included 3 biological replicates of wild-type tissue and 7 biological replicates from crel mutants (from 3 indpendent alleles: crel-1, crel-2, and crel-3). RNA-seq libraries were prepared using the Zymo RiboFree Total RNA Library Kit according to the manufacturer's instructions and libraries were sequenced using paired-end reads on an Illumina NovaSeq 6000.

RNA-seq data processing
Adapter sequences from raw RNA-Seq reads were trimmed using the "ILLUMINACLIP" function from the Trimmomatic software package using standard paired-end parameters. Trimmed reads were then mapped to the SL3.0 build of the tomato genome using the Spliced Transcripts Alignment to a Reference (STAR) aligner [84] with default parameters. Mapped reads were annotated to genes using the "featureCounts" function from the Subread software package. Differential gene expression analysis was calculated by using DeSeq2 package [85] comparing all crel mutant samples (3 independent alleles, 7 total biological replicates) together against wild-type data (3 biological replicates). See S4 Table for DeSeq2 output. For our analysis, only genes with adjusted p-values that met our Bonferroni-corrected p-value threshold of 1.395673e-06 (0.05/total # of genes detected; 35825) were considered. From there, differentially expressed genes were sorted into three bins: upregulated, downregulated, if log2 fold change from WT was over 0.6, less than -0.6 or from -0.59-0.59, respectively. See S2 Table for information on read processing numbers.

Enrichment analysis
The representation factor determines whether genes from one list (list A) are enriched in another list (list B), assuming that these genes behave independently. The representation factor is defined as: (number of genes in common between both lists) � (number of genes in the genome)/(number of genes in list A) � (number of genes in list B). p-values were calculated using exact hypergeometric probability [86]. Venn Diagram shows the overlap between all transcripts detected by RNA-seq with log2 fold change equal to or less than 0 (yellow) and transcripts with log2 fold change over 0 (blue) compared to genes associated with H3K27me3 peaks unique to WT (crel-lost sites; green) and unique to crel-2 (crel-gained sites; red). (TIF) S1 Table. Primers used in this work. (DOCX) S2 Table. Sequence read numbers for ChIP-seq. Two biological replicates of ChIP-seq for H3K27me3 were performed on shoot apices of WT and crel1-2 plants. The table indicates for each biological replicate the number of total sequencing reads obtained, the number and percentage mapping to the tomato genome, and the total number of reads remaining after filtering for mapping quality. Reads with a MAPQ score of 2 or greater were used for further analysis. (XLSX) S3 Table. Genes associated with crel-gained and crel-lost H3K27me3 sites and the peak coordinates for these sites.  Table. Gene ontology analysis of all Upregulated genes in crel and all genes associated with crel-lost H3K27me3 sites. (XLSX) S6 Table. Overlap between: 1. crel upregulated genes and genes associated with crel-lost H3K27me3 sites; 2. crel downregulated genes and genes associated with crel-increased H3K27me3 sites. (XLSX) S7 Table. Anlaysis of enrichemnt of "differentiation" and "morphogenesis genes in the lists of crel-upregulated genes, crel-downregulated genes, genes associated with crelincreased H3K27me3 siets and genes associated with crel-lost H3K27me3 sites.