<figure><img src="image_1.png" /></figure>

HAL
open science

# The quest for epigenetic regulation underlying unisexual flower development in Cucumis melo

David Latrasse, Natalia Rodriguez Granados, Alaguraj Veluchamy, Kiruthiga

Gayathri Mariappan, Claudia Bevilacqua, Nicolas Crapart, Céline Camps,

Vivien Sommard, Cécile Raynaud, Catherine Dogimont, et al.

<figure><img src="image_2.png" /></figure>

## To cite this version:

David Latrasse, Natalia Rodriguez Granados, Alaguraj Veluchamy, Kiruthiga Gayathri Mariappan, Claudia Bevilacqua, et al.. The quest for epigenetic regulation underlying unisexual flower development in Cucumis melo. Epigenetics and Chromatin, 2017, 10 (22), pp.1-17. 10.1186/s13072-017-0132-6. [hal-01594682](https://hal.science/hal-01594682)

HAL Id: hal-01594682

https://hal.science/hal-01594682

Submitted on 26 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
<figure><img src="image_4.png" /></figure>

# The quest for epigenetic regulation underlying unisexual flower development in *Cucumis melo*

David Latrasse<fnref n="1" />, Natalia Y. Rodriguez-Granados<fnref n="1" />†, Alaguraj Veluchamy²†, Kiruthiga Gayathri Mariappan², Claudia Bevilacqua³, Nicolas Crapart³, Celine Camps<fnref n="1" />, Vivien Somnard<fnref n="1" />, Cécile Raynaud<fnref n="1" />, Catherine Dogimont⁴, Adnane Boualem<fnref n="1" />, Moussa Benhamed<fnref n="1" />,² and Abdelhafid Bendahmane<fnref n="1" />*

**Abstract**

**Background:** Melon (*Cucumis melo*) is an important vegetable crop from the *Cucurbitaceae* family and a reference model specie for sex determination, fruit ripening and vascular fluxes studies. Nevertheless, the nature and role of its epigenome in gene expression regulation and more specifically in sex determination remains largely unknown.

**Results:** We have investigated genome wide H3K27me3 and H3K9ac histone modifications and gene expression dynamics, in five melon organs. H3K9ac and H3K27me3 were mainly distributed along gene-rich regions and constrained to gene bodies. H3K9ac was preferentially located at the TSS, whereas H3K27me3 distributed uniformly from TSS to TES. As observed in other species, H3K9ac and H3K27me3 correlated with high and low gene expression levels, respectively. Comparative analyses of unisexual flowers pointed out sex-specific epigenetic states of TFs involved in ethylene response and flower development. Chip-qPCR analysis of laser dissected carpel and stamina primordia, revealed sex-specific histone modification of MADS-box genes. Using sex transition mutants, we demonstrated that the female promoting gene, *CmACS11*, represses the expression of the male promoting gene *CmWIP1* via deposition of H3K27me3.

**Conclusions:** Our findings reveal the organ-specific landscapes of H3K9ac and H3K27me3 in melon. Our results also provide evidence that the sex determination genes recruit histone modifiers to orchestrate unisexual flower development in monoecious species.

**Background**

Melon (*Cucumis melo*) belongs to the *Cucurbitaceae* family that includes about 800 species. Besides melon, the *Cucurbitaceae* family comprises several important vegetable crops, such as cucumber (*C. sativus*), water-melon (*Citrullus lanatus*), squash and pumpkin (*Cucurbita* spp.), and many neglected cultivated species that are major food crops in many developing countries. For

decades, melon has been serving as a reference model organism for sex determination studies, and vascular fluxes [1, 2]. Melon also exhibits extreme genetic diversity for fruit traits, including fruit ripening, fruit shape, size, flesh color, texture, sweetness and aroma. Due to its relatively late entry into the genomic era, melon is a species with a barely explored epigenetic landscape. This limits the investigation of the cellular reprogramming of gene regulatory networks that drive development, growth and evolution of species of the *Cucurbitaceae* family.

Genome function, connectivity and structure rely on several epigenetic mechanisms that dynamically regulate the conformation of DNA, and thereby, its exposure to different regulators of gene transcription—e.g., transcription factors, non-coding RNAs, transposons [3–7].

<footnote>*Correspondence: abdelhafid.bendahmane@inra.fr</footnote>

<footnote>†Natalia Y. Rodriguez-Granados and Alaguraj Veluchamy contributed equally to this work</footnote>

<footnote n="1">¹Institute of Plant Sciences Paris-Saclay (IPSL), CNRS, INRA, University Paris-Sud, University of Evry, University Paris-Diderot, Sorbonne Paris-Cite, University of Paris-Saclay, Batiment 630, 91405 Orsay, France
Full list of author information is available at the end of the article</footnote>
At the chromatin level, DNA conformation can be modi-
fied by the addition or removal of chemical groups at
specific N-terminal residues of histones [8, 9]. Some of
these modifications are in charge of genome connectiv-
ity and heterochromatin stability, while others regulate
gene expression or repression along the euchromatic
regions. The distribution and abundance of these marks,
at the gene and genome-wide scales, define an epigenetic
landscape that directs cell differentiation and identity
[10–13]. In the plant kingdom, genome-wide analyses
of epigenetic landscapes have been performed in several
species including *Arabidopsis* [14–17], maize [18–20],
rice [21–23], common bean [24] and barley [25].

Histone methylation and acetylation are epigenetics
marks that are usually associated with various biolog-
cal processes ranging from transcriptional regulation
to epigenetic silencing via heterochromatin assembly.
While histone acetylation is commonly associated with
transcription activation, histone methylation can either
promote or inhibit gene transcription [5, 11, 18, 26–29].
Histone H3 lysine 9 acetylation (H3K9ac) is an epigenetic
mark commonly found in active promoters [5, 20, 30, 31].
Genome-wide studies in *Arabidopsis* have found around
5200 non-Transposable Element (TE)-related genes that
are targeted and regulated by H3K9ac [31]. In contrast,
H3 lysine 27 trimethylation (H3K27me3) is commonly
associated with gene repression, thereby contributing
to the mitotic heritability of PRC2 (Polycomb repressive
complex)-mediated silencing of numerous developmen-
tal and homeotic genes [29, 32, 33]. In the plant kingdom,
the study of these epigenetic marks has elucidated the
epigenetic component regulating plant growth and adap-
tation to their harsh environment [34–44].

The majority of angiosperms are hermaphrodite pro-
ducing exclusively bisexual flowers. Sex determination
is a developmental evolutionary process that leads to
unisexual flowers. Monoecious species, such as melon,
exhibit male and female flowers on the same plant. Dioe-
cious species have separate male and female individuals
[45, 46].

In melon, floral primordia are initially bisexual with sex
determination occurring by the selective developmen-
tal arrest of either the stamen or the carpel primordia,
resulting in unisexual flowers [47, 48]. This sexual organ
arrest is genetically governed by the interplay of alleles
of the *andromonoecious (M)*, *androecious gene (A)* and
*gynoecious (G)* genes [47, 49–51]. The cloning and char-
acterization of *M*, *G* and *A* genes have shown that the
*gynoecious (G)* gene encodes a zinc finger transcription
factor, *CmWIP1* [47], the *andromonoecious (M)* gene
and the *androecious gene (A)* encode ethylene biosynthe-
sis enzymes: *CmACS-7* [49] and *CmACS11* [51], respec-
tively. Genetic analysis revealed a mechanistic model in

which expression of the carpel inhibitor, *CmWIP1*, is
dependent on non-expression of *CmACS11* and expres-
sion of the stamina inhibitor, *CmACS-7*, is dependent on
non-expression of *CmWIP1*. In monoecious plants, male
flowers result from non-expression of *CmACS11* that
permits *CmWIP1* expression. Female flowers develop
on branches expressing *CmACS11*, which represses the
expression of *CmWIP1*, and thus, releasing the expres-
sion of *CmACS7* inhibiting stamina development
[51–53].

Now that the identity of the androecy, monoecy and
gynoecy sex genes are revealed and genetic pathway
controlling sex determination in *Cucurbitaceae* discov-
ered, among the next challenges is to decipher how the
sex determination signals are perceived and how the
information is translated to cause organ-specific abor-
tion at the flower and the plant level. Several studies have
reported a wide-range of epigenetic processes and envi-
ronmental cues that determine sex in animals [54–57].
In plants, it is well known that sex ratio and determina-
tion can be influenced by different environmental fac-
tors such as day length and light intensity [58], water
restriction [59] and some plant hormones such as ethyl-
ene [49–51, 60–63], auxin [64] and gibberellins [65–67].
This multifactorial regulation of sex determination is still
poorly characterized from the epigenetic point of view,
where limited literature have reported the association for
instance of small RNA (sRNAs) transcription and DNA
hypomethylation with sex determination [54, 68–71].

In melon, the epigenetic control of sex determination
remains poorly characterized and is mainly attributed
to a transposon-induced DNA methylation that leads
to gynoecy [47]. To gain new information on the epige-
netic control of sex determination, we have examined the
genome-wide landscape of H3K27me3 and H3K9ac his-
tone modifications and gene expression dynamics, in five
melon organs, focusing our analysis on unisexual flow-
ers. The combination of epigenomic and gene ontology
analysis pointed out sex-specific epigenetic states of TFs
involved in ethylene response and sexual organs devel-
opment. Furthermore, using sex transition mutants, we
demonstrated that the female promoting gene represses
the expression of the male promoting genes via deposi-
tion of H3K27me3.

Results

Genome-wide landscape of H3K9ac and H3K27me3
in melon

Genome-wide studies in *Arabidopsis* [14–17], and other
plant species such as maize [19, 20], rice [21–23], and
barley [25], have evidenced the important role of H3K9ac
and H3K27me3 in gene activation and repression,
respectively. The roles of these histone modifications in
melon development remain unknown and represent a
limit to fully understand how thousands of bioprocesses
are regulated. To determine the genomic landscape and
the organ specificity of those marks, we performed ChIP-
seq analyses using H3K27me3 and H3K9ac antibodies
on different organs, namely fruits, leaves, roots, female
and male flowers. Two replicates per tissue were pro-
cessed and present a correlation coefficient of at least
0.99 (Additional file 1: Figures S1 and S2). A minimum
of 14 millions of mapped reads were obtained in the
replicates implemented (Additional file 2: Table I). The
MACS2 algorithm, which is designed to detect typical
histone peaks [72], was used to determine loci that are
significantly enriched with H3K9ac or with H3K27me3
(Additional file 1: Figures S3 and S4). We identified
18,424 H3K9ac marked genes in leaves, 18,024 in roots,
13,655 in fruits, 17,660 in male flowers and 16,348 in
female flowers. On the other hand, we identified 6486
H3K27me3 marked genes in leaves, 7741 in fruits, 6823
in roots, 6055 in male and 6054 in female flowers (Fig. 1a;
Additional file 3: Table II).

Next we analyzed the distribution of those two marks
at the chromosome and gene levels. For the sake of sim-
plicity, we will describe mainly results obtained in leaves,
and results obtained in the other organs are shown in the
additional file. The chromosomal distributions of anno-
tated genes, H3K9ac and H3K27me3 and in each organ
on the 12 chromosomes, are illustrated in Additional
file 1: Figure S3. At the chromosomal scale we observed
an enrichment of both marks in gene-rich regions, which
are localized at chromosome arms, compare to centro-
meric and pericentromeric regions, where the majority of
transposable and retro transposable elements are located
[73], consistent with the role of these histone marks in
the control of gene expression (Fig. 1b; Additional file 1:
Figures S5 and S6). In order to detail the H3K27me3
and H3K9ac distributions at the gene level, the peaks
obtained in each tissue for both modifications were ana-
lyzed and averaged within a 4000-bp region, flanking the
TSS and TES of all the annotated genes. We found that

the peak length of H3K9ac ranged from 200 bp to 700 bp
(Additional file 1: Figure S7), located preferentially at the
TSS regions, which correspond to the first nucleosomes
of the genes (Fig. 1c; Additional file 1: Figure S8). In con-
trast, H3K27me3 peaks presented an averaged length
that ranged from 2500 to 2700 bp, which covered the
entire gene body (Fig. 1c; Additional file 1: Figure S7).
Those patterns were observed in all the organs (Addi-
tional file 1: Figures S7 and S8) and were consistent with
previous studies on different plant species, highlighting
conserved aspects of the epigenetic system in the plant
kingdom. Integration of H3K9ac and H3K27me3 datasets
showed a clear anti-correlation between those two marks
(Fig. 1d). Altogether, these results showed that in melon,
as in other plant species, H3K27me3 and H3K9ac are
not enriched in intergenic regions but distributed along
the gene body, supporting the role of these two marks in
gene regulation in an exclusive manner.

Identification of organ-specific H3K27me3 and H3K9ac
landscapes

We next aimed to determine the organ-specific profile
of both H3K27me3 and H3K9ac marks, and to deter-
mine whether these chromatin modifications contribute
to the establishment of organ-specific gene expression
programs. To this end, first, for each organ, we deter-
mined the genes that are marked either by H3K9ac or by
H3K27me3 or by both. Secondly, we performed a com-
parison between all the possible paired combinations of
the five different tissues using DiffReps, a highly sensi-
tive program for the detection of differential sites from
ChIP-seq data [74]. Third, a Venn diagram was generated
to identify genes that were over-represented in the cho-
sen organ regardless of which organ it was compared to
and thus to fish out genes that are specifically marked in
this organ (Fig. 1e, g; Additional file 1: Figures S9–S12;
Additional file 4: Table III). To determine the biological
function of these genes, we performed a gene ontology
analysis (GO) using Plant MetGenMAP [75]. Organ-
specific hyperacetylated and hypomethylated genes,

(See figure on next page.)

**Fig. 1** Genome-wide distribution and tissue-specificity of H3K9ac and H3K27me3 in melon. **a** H3K27me3 and H3K9ac target genes. Number of genes presenting H3K9ac (green bar) and H3K27me3 (red bar) peaks were quantified and plotted for the five different tissues. **b** H3K27me3 and H3K9ac distribution at the chromosome level. Distribution of H3K9ac (green, upper panel), H3K27me3 (red, medium panel) in leaves and annotated genes (gray, bottom panel) are plotted along the chromosome 1. **c** Average tag density profile of H3K27me3 and H3K9ac along the gene body. ChIP-Seq densities of equal bins were plotted along the gene body and 2-kb region flanking the TSS or the TES. **d** Correlation between the genome-wide distribution of H3K27me3 and H3K9ac. Heat map representing the tag density distribution of H3K27me3 and H3K9ac across all genes and a 2 kb flank. **e** Leaf-specific H3K9 hyper-acetylated genes. Paired comparisons of leaves vs. all the organs are shown as a Venn diagram, leaf-specific hyperacetylated genes correspond to the central overlap (highlighted in green). **f** GO terms of biological processes enriched for leaf-specific hyperacetylated genes. **g** Leaf-specific H3K27 hypo trimethylated genes. Paired comparisons of leaves versus all the tissues are shown as a Venn diagram, leave-specific hypo methylated genes correspond to the central overlap (highlighted in red). **h** GO terms of biological processes enriched for leave-specific hypo trimethylated genes. *p* values for each enriched class are presented in Additional file 9: Table VIII
<figure><img src="image_5.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>e</figcaption></figure>

H3K9ac hyper-acetylated genes in leaves

<figure><img src="image_7.png" /><figcaption>g</figcaption></figure>

H3K27me3 hypo-methylated genes in leaves

<figure><img src="image_8.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_12.png" /></figure>
displayed biological processes terms consistent with
metabolic needs and physiological aspects of their corre-
sponding organ, providing evidence for the robustness of
our data (Fig. 1f, h; Additional file 1: Figures S9–S12). For
example in leaves, we found 1335 specifically hyperacety-
lated genes involved in photosynthesis, leaf development,
stomatal morphogenesis (Fig. 1f). Consistent results were
also obtained for fruits, roots, female and male flowers
(Additional file 1: Figures S9–S12).

These results suggest that both H3K9ac and H3K27me3
display an organ-specific profile crucial for cell identity
and physiology.

Relationship between gene expression and histone
modifications in melon

To connect H3K9ac and H3K27me3 histone marks with
gene expression, we generated and integrated RNA-
seq data from the five different organs. Firstly, we con-
firm that H3K9ac and H3K27me3 in *Cucumis melo* are,
respectively, associated with gene expression and gene
repression (Fig. 2a). Secondly, in order to determine
whether the level of acetylation or methylation could
be correlated with gene expression, we divided all the
protein coding genes into four quantiles based on their
expression levels and plotted their H3K9ac or H3K27me3
profile (Fig. 2b). We observed that H3K9ac occupancy
increases with expression level, while H3K27me3 showed
the opposite pattern, where it displayed a high enrich-
ment in the lowest-expressed genes but it remained unre-
sponsive to intermediary and high expression levels. This
results suggest that in *Cucumis melo* the more a gene is
marked by H3K9ac and H3K27me3, the more it will be
expressed and repressed, respectively (Fig. 2b).

To confirm the relationship between gene expression
and organ-specific landscape of H3K9ac and H3K27me3,
we identified genes specifically expressed in each organ
(Fig. 2c; Additional file 5: Table IV) and we compared
them with organ-specific H3K9ac marked genes. Inter-
estingly, in the case of leaves, we found that a signifi-
cant part of genes specifically hyperacetylated (43%) and
hypomethylated (32%) are specifically expressed in this
organ (Fig. 2d, e). Gene ontology analyses of these genes
showed a significant enrichment in photosynthesis and
response to light. In the case of roots, flowers and fruits
we observed the same tendency, suggesting that both
H3K9ac and H3K27me3 controlled organ-specific gene
expression (Additional file 1: Figure S13).

Female and male flower-specific epigenome landscapes

Beyond elucidating general aspects of the melon epige-
nome, the analyses we performed so far evidenced how
H3K27me3 and H3K9ac organ-specific distributions
can be crucial to study gene expression regulation and

organ identity. Thus, we directed epigenomic and tran-
scriptomic analyses toward the elucidation of female and
male flower-specific epigenome landscapes and their
role in sex determination. To this end, we integrated
RNA-seq and, H3K9ac and H3K27me3 ChIP-seq data of
male and female flowers in order to identify genes that
are epigenetically regulated in a sex-specific manner
(Fig. 3; Additional file 6: Table V). We observed that 22%
of genes that are hypoacetylated in female flowers com-
pared to male flower are downregulated and only 3% of
the set of hypoacetylated genes are upregulated (Fig. 3a).
Regarding hyperacetylated genes, we also found that a
significant proportion of them (28%) are upregulated in
female flowers compared to male flower and only 5% are
dow-regulated. Furthermore, gene ontology analysis of
female-specific epigenetically regulated genes showed
enrichment in biological processes such as carpel devel-
opment, response to ethylene and response to jasmonic
acid. For male-specific epigenetically regulated genes, we
observed enrichment in biological processes such as, sta-
men development and lipid metabolic process (Fig. 3a).
The same analyses were performed for hypo- and hyper-
methylated genes, obtaining similar results and con-
firming that a significant portion of sex-specific genes is
epigenetically controlled (Fig. 3b).

After analyzing the biological processes involving the
sex-specific epigenetically regulated genes, we classified
these genes based on their molecular function. We found
enrichment for transporter activity, catalytic activity,
oxidoreductase activity and interestingly, for transcrip-
tion factor (TF) activity (Fig. 3c; Additional file 1: Figures
S14–S18; Additional file 7: Table VI). Since transcription
factors are the best described and commonly associated
components of gene expression regulation, we further
focused on them and grouped them according to their TF
families (Fig. 3d). We found that the ethylene-responsive
factor (ERF) family displayed the highest number of gene
members, consistent with previous studies that reported
the role of ethylene and ethylene transduction pathways
in sex determination in the cucurbitaceae family (22). In
addition, we found other TF families such as Zinc Fin-
ger (17 genes), basic Helix-Loop-Helix (15 genes), MYB
(13 genes) and NAC families (13) (Fig. 3d). Interestingly,
we also detected 11 genes belonging to the MADS-box
TF family. Several members of this family are conserved
master regulators of flower ontogenesis [76–78]. When
we assessed the identity of these 11 genes, we observed
that all of them are potential homologs of different *Arabi-
dopsis* genes that are part of the ABC model of flower
development (Additional file 7: Table VI).

Altogether, these results highlight the organ-specificity
of H3K27me3 and H3K9ac as well as its role in the regu-
lation of key genes encoding transcription factors, which
a

<figure><img src="image_13.png" /></figure>

b

<figure><img src="image_14.png" /></figure>

c

<figure><img src="image_15.png" /></figure>

d

<figure><img src="image_16.png" /></figure>

e

<figure><img src="image_17.png" /></figure>
(See figure on previous page.)

**Fig. 2** Correlation of H3K27me3 and H3K9ac with gene expression level and organ-specificity in leaves. **a** Correlation between H3K27me3, H3K9ac and gene expression. Genes specifically marked with H3K9ac or H3K27me3 in leaves were clustered based on the log2 value of their normalized tag densities of H3K9ac, H3K27me3 and RNA-seq data. Genes displaying the highest enrichment for H3K27me3 and H3K9ac were used with their corresponding gene expression levels in a heat map. **b** Correlation of between H3K27me3, H3K9ac and gene expression level. All the melon protein-coding genes were divided in 4 quantiles according to their gene expression levels (lowest and highest expression level corresponding to *green* and *red*, respectively). For each quantile the number of H3K27me3 and H3K9ac mapped reads was averaged and plotted along the gene body and 2-kb region flanking the TSS or the TES. **c** Gene clustering analyses based on organ-specific expression. Genes were sorted depending on their expression in each organ. 5 clusters, each one corresponding to an organ, are illustrated and numbered from 1 to 5. Number of genes of each cluster is also indicated. **d** Correlation of hyper-H3K9ac with tissue-specific gene activation. Venn diagram illustrating the genes whose high and tissue-specific expression is correlated with their tissue-specific hyper-acetylation (Overlapping set of genes). A Chi-squared test was done and confirmed the observed correlation. For this set of genes, GO terms of biological processes are shown below. **e** Correlation of hypo-H3K27me3 with tissue-specific gene activation. Venn diagram illustrating the genes whose high and tissue-specific expression is correlated with their tissue-specific hypomethylation (Overlapping set of genes). A Chi-squared test was done and confirmed the observed correlation. For this set of genes, GO terms of biological processes are shown below. **p** values for each enriched class are presented in Additional file 9: Table VIII

regulate processes such as ethylene response and flower
development.

Carpel and stamen epigenetic landscapes highlight
sex-specific chromatin states of transcription factor genes

Since sex determination in melon and other cucur-
bits relies on the selective regulation of carpel and sta-
men development, we analyzed in detail some melon
TFs, which are potential homologs of *Arabidopsis* genes
involved in the control of sexual organs development.
For such genes, we first assessed the genic distribution
of H3K27me3, H3K9ac and RNA-seq reads in male and
female flowers (Fig. 4a). According to previous studies in
*Arabidopsis*, stamen development relies on the expression
of *APETELA3* (*AP3*), *PISTILATA* (*PI*) and *AGAMOUS*
(*AG*), whereas carpel development depends on the
expression of ovule-specific genes such as *SEEDSTICK*
(*STK*) and *SHATTERPROOF 1* (*SHP1*) [78]. As expected,
when we analyzed both H3K9ac and H3K27me3 epige-
netic landscapes of stamen identity genes together with
their expression, we found that they are hyperacetylated,
hypomethylated and highly expressed in the male flower
when compared to the female flower. Carpel-promoting
genes, however, presented an enrichment of H3K9ac and
high gene expression in comparison to male flowers, but

no significant differences were observed for H3K27me3,
which may result from the cell heterogeneity of the sam-
ple or from the cooperative action of additional histone
repressive marks.

Tissue heterogeneity is an important variable that
limits the interpretation of ChIP-seq data and the iden-
tification of epigenetic events occurring in a specific
cell population. This is particularly true in this con-
text of flower development, where it was reported that
organ identity relies on the whorl-specific expression
of MADS-box genes [76–78]. We therefore decided to
overcome such limitation, by performing an H3K27me3
and H3K9ac ChIP-qPCR in isolated carpels and sta-
mens. We isolated carpels from female and stamens from
male flowers through the Laser Capture Microdissec-
tion (LCM) technique [79] and performed ChIP using
H3K27me3 and H3K9ac antibodies. The immunopre-
cipitated chromatin was then amplified using primers
annealing a neighboring sequence to the TSS of *AP3*, *PI*,
*AG*, *STK-like* genes; which corresponded to the highest
peak of H3K9ac observed along the gene body. A nega-
tive control, *MELO3C115188*, that does not display any
H3K9ac or H3K27me3 signals was also used (Additional
file 1: Figure S19). In all cases, the obtained ChIP-qPCR
results were consistent with the ChIP-seq data previously

(See figure on next page.)

**Fig. 3** H3K9ac and H3K27me3 differentially marked genes between female and male flowers. **a** Venn diagram showing the relationship between H3K9 hyper-acetylated or H3K9 hypo-acetylated genes and up- or down-regulated genes in females flowers compared to male flowers (*middle panel*). Hyper-acetylated genes are predominantly up-regulated and hypo-acetylated genes are down-regulated, as highlighted in the *pink* and *blue boxes* respectively. Gene ontology analysis (GO) of the genes in the *pink* and *blue boxes* are shown in the *left panel* and the *right panel* respectively. Histograms of the values highlight the enrichment of genes compared to the reference. **b** Venn diagram showing the relationship of H3K27 hyper-methylated or H3K27 hypo-methylated genes and up- or down-regulated genes in females flowers compared to male flowers (*middle panel*). Hyper-methylated genes are predominantly down-regulated and hypo-acetylated genes are up-regulated, as highlighted in the *purple* and *blue boxes* respectively. Gene ontology analysis (GO) of the genes in the *purple* and *blue boxes* are shown in the *left panel* and the *right panel* respectively. Histograms of the values highlight the enrichment of genes compared to the reference genome (*black bars*). **p** values for each enriched class are presented in Additional file 9: Table VIII. **c** Percentage of transcription factors among the H3K9ac and H3K27me3 differentially regulated genes between male and female flowers. **d** Percentage of different known transcription factor families among the H3K9ac and H3K27me3 differentially regulated genes between male and female flowers
a

H3K9ac hyper-acetylated and up-regulated genes
in female flowers vs. male flowers

<figure><img src="image_18.png" /></figure>

Reference genome

H3K9ac hypo-acetylated and down-regulated genes
in female flowers vs. male flowers

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

b

H3K27 hyper-methylated and down-regulated genes
in female flowers vs. male flowers

<figure><img src="image_21.png" /></figure>

Reference genome

H3K27me3 hypo-methylated and up-regulated genes
in female flowers vs. male flowers

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

c

H3K9ac or H3K27me3 regulated genes (1093 genes)

Reference genome (27427 genes)

<figure><img src="image_24.png" /></figure>

d

H3K9ac or H3K27me3 regulated genes (1093 genes)

Reference genome (27427 genes)

<figure><img src="image_25.png" /></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 4** Carpel and Stamen epigenetic analysis of melon transcription factors. **a** H3K9ac and H3K27me3 profiles and RNA levels of four transcription factors in female and male flowers. Two TFs specific of the male flowers (MELO3C003778 and MELO3C0100515) and two specific of the female flowers (MELO3C002691 and MELO3C022209) were selected for epigenetic profiling in isolated carpels and stamens of melon flowers. The distributions of peaks and RNA reads obtained in our genome-wide analyses are plotted above the genome model for each TF. **b** Stamens and carpels of melon flowers were dissected by Laser Microdissection (LMD) and ChIP-qPCR assays were performed with anti-H3K9ac and anti-H3K27me3 antibodies. The *upper panels* show an example of dissections performed with the LCM technique. The *bar charts* represent the qPCR assays on immune-precipitated DNA using primers annealing a neighboring sequence to the TSS as indicated in *panel A* (purple bar)</figcaption></figure>
generated from entire flowers; however, the differences observed between stamen and carpel are more evident, presenting at least a fold change of 2.3. This allowed us to observe clear differences in H3K27me3 occupancy between male and female organs (Fig. 4b).

Altogether, these results indicate that H3K27me3 and H3K9ac target MADS-box TFs regulating carpel and stamen identity and are presented in a sex-specific manner, thereby suggesting the role of these two epigenetic marks in the regulation of their expression. Furthermore, these results showed how tissue homogeneity could represent an advantage for the elucidation of cell-specific epigenetic landscapes in melon and other plant species.

## Role of *CmACS11* and *CmWIP1* genes in sex-specific epigenome acquisition

We previously showed that sex determination in melon relies on the interplay between alleles of three sex determination genes, *M*, *G* and *A* [47, 49–51]. In monoecious melon, most of the flowers are male and female flowers develop at the youngest nodes of the growing vines (Fig. 5a). Male flowers result from the expression of the male promoting gene, *CmWIP1*. Female flowers develop on branches expressing *CmACS11*, which represses the expression of *CmWIP1*. Androecious plants result from a loss-of-function of *CmACS11* leading to expression of *CmWIP1* in all flowers on a plant. Gynoecious plants are obtained by inactivation of *CmWIP1* function [47, 51] (Fig. 5a).

To bring new insight on the role of these sex-determining genes in establishing epigenetic signatures important for sex determination in melon, we investigated the H3K27me3 genome-wide landscape in male and female flowers of wild-type monoecious melon and two androecious and gynoecious isogenic lines, carrying a loss-of-function mutation in *CmACS11* and *CmWIP1*, respectively.

To determine the proportion of genes whose differential epigenetic state between male and female flowers can depend on *CmACS11* and *CmWIP1*, we compared male and female flowers from the wild type with flowers occupying the same position on the inflorescence from androecious and gynoecious sex transition mutants, respectively (Fig. 5b; Additional file 1: Figures S20–S22). The gene lists generated by these comparisons are illustrated in a Venn diagram, where the overlapping genes correspond to *CmWIP1*—(Additional file 1: Figure S20B) and *CmACS11*-related (Additional file 1: Figure S20C) H3K27me3 changes linked to sex identity. We observed that 54% (1200/2197) of sex-specific regulated gene was under the control of the *G* locus and that 60% (1324/2197) was under the control of the *A* locus (Fig. 5b). Altogether, these data suggested that *CmACS11*

and *CmWIP1* are two master regulator controlling sex-specific epigenome acquisitions.

The androecious gene *CmACS11* plays central role in sex determination as it controls the expression the gynoecious gene *CmWIP1* [51]. Knowing from this study that *CmWIP1* is epigenetically regulated through histone methylation (Additional file 1: Figure S23), we predicted that *CmACS11* represses the expression of the *CmWIP1* by modifying its methylation status. To test this hypothesis, we analyzed the H3K27me3 status of *CmWIP1* in flowers that do or do not express *CmACS11* and in flowers impaired in *CmACS11* function. We performed an H3K27me3 ChIP followed by qPCR using 9 different primers covering the entire *CmWIP1* gene. As predicted, we observed that *CmWIP1* is hypermethylated in female buds that express *CmACS11* but not in male buds that do not express *CmACS11*. The hypermethylation of *CmWIP1* is also lost in *CmACS11* loss of function mutant (Fig. 5c). The hypomethylation correlates with the expression *CmWIP1* and the hypermethylation with the non-expression of *CmWIP1*. Based on this, we concluded that the *A* gene represses the expression of *CmWIP1* via deposition of H3K27me3.

## Discussion

Deposition of chromatin modifications is instrumental for the concerted expression or repression of thousands of genes, thereby allowing the establishment of organ or tissue specific transcriptomes [10–13]. H3K9ac and H3K27me3 are among the main chromatin marks involved in the control of gene expression [5, 20, 30–33]. In this study, we generated genome-wide and organ-specific maps of H3K27me3 and H3K9ac in melon, thereby generating a repertoire for the elucidation of mechanisms controlling the expression of the melon genome. This work will provide new information for the analysis of melon development or several agricultural traits, that have until now been mainly studied via classical genetic approaches. The distributions of H3K9ac and H3K27me3 in the melon genome displayed common chromosomal and genic distributions when compared to other plant species such as *Arabidopsis* [14–17], maize [19, 20] and rice [21–23]: H3K27me3 and H3K9ac were mainly enriched at the chromosome arms, resembling the distribution of genes, but were confined to the gene body and not spread along intergenic regions. Within the gene body, H3K27me3 distributes uniformly from the TSS to the TES, whereas H3K9ac is mainly enriched at the TSS.

As described in other plant species, we observed that H3K9ac and H3K27me3 marks display opposite correlation with gene expression. Interestingly, we found a direct proportionality between H3K9ac and gene expression levels, thereby suggesting the crucial role of H3K9ac in
<figure><img src="image_27.png" /></figure>

gene expression regulation and fine-tuning. In contrast,
H3K27me3 did not display an evident proportionality
with gene expression, being unresponsive to intermediate
gene expression levels as it has been observed for barley

[25] and maize [19]; however, it was clearly enriched at
the lowest-expressed genes, suggesting that this histone
mark is involved in the full repression of genes rather in
the fine-tuning of their expression.
Beyond their role in gene expression regulation, we observed that both H3K9ac and H3K27me3 are important in the control of organ identity and development, thus mediating the formation of different cell types and functions from a single genome. A considerable portion of H3K9ac and H3K27me3 target genes was expressed in an organ-specific manner, the differential and unique deposition of these two histone marks in each organ, creates an epigenetic landscape on genes that is crucial for its metabolic needs and physiology.

In the case of flowers, female-specific H3K9 hyperacetylated and H3K27 hypo-trimethylated genes are involved in carpel and flower development as well as in several other bioprocesses previously reported to play a role in flower feminization in cucurbits. Such is the case for the hormone ethylene, whose synthesis and perception have important roles in development of female unisexual flowers. The expression of *CmACS7* and *CmACS11*, two enzymes that are part of the ethylene biosynthesis pathway, and *CsETR1*, one of the reported ethylene receptors, have been already placed in the regulatory network promoting female flower development in cucurbits [47, 49–51, 80–83]. On the other hand, other genes related to jasmonic acid response were specifically modified with active marks in the female flower. Consistent with this, jasmonic acid (JA) was also shown to be associated with sex determination in maize [84]. Specifically, downregulated genes (H3K9 hypo-acetylated and H3K27 hyper-trimethylated) in the female flower pointed out bioprocesses such as lipid metabolism, which has been previously reported to play a role in pollen-tube growth and pollen-stigma interactions [85].

When we assessed the molecular function of sex-specific epigenetically regulated genes, we found enrichment for transcription factor, DNA and nucleotide binding activities. We, therefore, categorized these TFs according to their corresponding protein families and found that the majority of them belong to the ethylene-responsive factors [86–88]. These results suggest that histone modifications regulate ethylene transduction pathway in melon as it has previously reported in Arabidopsis [89].

Interestingly, we also found genes belonging to the MADS-box TF family, which are well recognized for their important role in flower whorl identity and is in concordance with previous reports regarding the regulation of *AP1*, *AG* and *SEP3* by the histone demethylase REF6 that removes H3K27me3 [90, 91], and polycomb-dependent regulation of *AG* [92, 93].

The elucidation of epigenetic control of the class C genes for the development of unisexual flowers demands the assessment of the whorl-specific epigenetic state of carpel- and stamen-related identity of MADS-box genes. Tissue heterogeneity has been considered a persistent

constraint for the identification of epigenetic events occurring in a specific cell population thus, masking the natural complexity of gene expression regulation. So far in the plant kingdom, few technical improvements to this limitation have been reported, especially due to the complexity of performing ChIP on small samples. In this study, we provide the first protocol to overcome such limitation that can be transferred to any plant species and demonstrate that epigenetic changes that are undetectable with standard methods using whole organs can be identified by using purified cell populations. By this method, we also demonstrated that stamen—(i.e., *PI* and *AP3*) and carpel-related (i.e., *STK*) identity genes present a contrasting H3K27me3 and H3K9ac landscapes. This suggests that the sex determination genes recruit histone modifiers to target organ identity genes in the sexual whorls of the flower, leading to unisexual flower development.

In the cucurbit sex determination pathway, the gynoecious gene, *CmWIP1*, plays a central role in unisexual flower development. From this analysis, we showed that *CmACS11* represses *CmWIP1* expression by inducing H3K27me3 at this locus (Fig. 5d). This data suggest that the local production of ethylene via *CmACS11* in companion cells is able to activate histone methyltransferase to inhibit *CmWIP1* expression in a defined spatio-temporal manner.

Detailed knowledge about spatio-temporal gene expression and epigenetic regulation dynamics is pivotal for a comprehensive understanding of development of unisexual flowers. In species recalcitrant to plant transformation, such as melon, this has been hampered mainly by difficulties in isolating sufficient amounts of tissues from distinct organ primordia for chromatin analysis. In this context, using laser capture micro-dissection, we succeeded in generating high-quality chromatin, from dissected carpel and stamina primordia, from unisexual flowers, and discovered that the sex determination genes recruit histone modifiers to orchestrate unisexual flower development in monoecious species. Genome-wide analyses of chromatin modifications in the dissected reproductive organs from unisexual and hermaphrodite flowers and from key histone modifier mutants will likely to bring new insight on the epigenetic control of unisexual flower development.

## Methods

### Plant material and growth conditions

Monoecious melon (*Cucumis melo*) cultivar charentais (*Cucumis melo* L. subsp. *melo* var *cantalupensis*) were used for H3K9ac and H3K27me3 genome-wide descriptive analyses. For each replicate, 7-day-old seedlings were transferred to individual pots and incubated in a growth
chamber for 30 days (long day conditions, temperature: 27 °C (day) and 21 °C (night), relative humidity: 60%). Young male and female flowers (3 mm length; Developmental stage 8 [52]) were collected, organized according to their location in the plant (i.e., main stem or branches) and conserved at -80 °C for epigenomic and transcriptional analyses.

For the determination of *a*- and *g*-related epigenetic changes, we used two families of EMS-treated *C. melo* lines carrying missense mutations in *CmWIP1* and in *CmACS11*, both previously reported to lead to gynoecy [47] and androecy [51], respectively. Germination, growth and flower collecting conditions were the same to the previously described for the monoecious WT plants.

## Gene expression analysis

Total RNA were extracted from leaves, root, fruit, male and female flowers using the Nucleospin RNA kit (Macherey-Nagel), according to the manufacturer's instructions. For RT-QPCR analyses, first-strand cDNA was synthesized from 1 µg of total RNA using Improm-II reverse transcriptase (A3802, Promega) according to the manufacturer's instructions. 1/25th of the synthesized cDNA was mixed with 100 nM of each primer and LightCycler® 480 Sybr Green I master mix (Roche Applied Science) for quantitative PCR analysis. Products were amplified and fluorescent signals acquired with a LightCycler® 480 detection system. The specificity of amplification products was determined by melting curves. ACT2 was used as internal control for signals normalization. Exor4 relative quantification software (Roche Applied Science) automatically calculated relative expression level of the selected genes with algorithms based on ΔΔCt method. Data were from duplicates of at least two biological replicates. The sequences of primers can be found in Additional file 8: Table VII.

For RNA-seq analysis, libraries were synthetized using NEBNext Ultra Directional RNA library Preparation Kit (NEB) according to the manufacturer's instructions. Two biological replicates were run for each tissue. Raw reads from RNA-Seq were first adaptor trimmed and quality filtered using Trimmomatic [94]. Filtered reads were mapped to the *Cucumis melo* genome v3.5.1 using TopHat v2.0.9 [95]. Transcript quantification was derived using Cufflinks v2.2.0 [96]. CuffDiff is used for differential expression analysis (*p* value 0.05; statistical correction: Benjamini Hochberg; FDR: 0.05). A cutoff of 0.5 fold up- or down-regulation has been chosen to define significant differential expression. CummeRbund v2.0.0 was used for visualization of differential analysis [97]. Functional annotation using GO terms was performed using the Plant MetGenMAP tool [75]. All *p* values are shown in Additional file 9: Table VIII.

## Chromatin immunoprecipitation experiments

ChIP assays were performed using anti-H3K9ac (Millipore, ref. 07-352) or anti-H3K27me3 (Millipore, ref. 07-449) antibodies, using a procedure adapted from Veluchamy et al. [17]. Briefly, after plant material fixation in 1% (v/v) formaldehyde, tissues were homogenized, nuclei isolated and lysed. Cross-linked chromatin was sonicated using a water bath Bioruptor UCD-200 (Diagenode, Liège, Belgium) (30 s on/30 s off pulses, at high intensity for 60 min). Protein/DNA complexes were immunoprecipitated with antibodies, overnight at 4 °C with gentle shaking, and incubated for 1 h at 4 °C with 50 µL of Dynabeads Protein A (Invitrogen, Ref. 100-02D). The beads were washed for 2 × 5 min in ChIP Wash Buffer 1 (0.1% SDS, 1% Triton X-100, 20mM Tris-HCl pH 8, 2 mM EDTA pH 8, 150 mM NaCl), 2 × 5 min in ChIP Wash Buffer 2 (0.1% SDS, 1% Triton X-100, 20 mM Tris-HCl pH 8, 2 mM EDTA pH 8, 500 mM NaCl), 2 × 5 min in ChIP Wash Buffer 3 (0.25 M LiCl, 1% NP-40, 1% sodium deoxycholate, 10 mM Tris-HCl pH 8, 1 mM EDTA pH 8) and twice in TE (10 mM Tris-HCl pH 8, 1 mM EDTA pH 8). ChIPed DNA was eluted by two 15-min incubations at 65 °C with 250 µL Elution Buffer (1% SDS, 0.1 M NaHCO₃). Chromatin was reverse-crosslinked by adding 20 µL of NaCl 5 M and incubated over-night at 65 °C. Reverse-cross-linked DNA was submitted to RNase and proteinase K digestion and extracted with phenol-chloroform. DNA was ethanol precipitated in the presence of 20 µg of glycogen and resuspended in 50 µL of nuclease-free water (Ambion) in a DNA low-bind tube.

For ChIP-qPCR experiments, fold enrichment of targets in ChIPed DNA relative to input was calculated from an average of three replicate qPCRs. The sequences of primers can be found in Additional file 8: Table VII. Positions of the amplified regions on the different loci are indicated in Fig. 5.

For ChIP-seq assays, 10 ng of IP or input DNA was used for ChIP-Seq library construction using NEB-Next Ultra II DNA Library Prep Kit for Illumina (New England Biolabs) according to manufacturer's recommendations. For all libraries, twelve cycles of PCR were used. The quality of the libraries was assessed with Agilent 2100 Bioanalyzer (Agilent), and the libraries were subjected to high-throughput sequencing by Illumina Sequencing technology.

## Computational analysis of ChIP-seq

Preprocessing of sequencing reads for quality was performed using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Filtering and trimming of reads was done using Trimmomatic with the following parameters: Minimum length of 36 bp; mean Phred quality score greater than 30; leading and trailing
bases removal with base quality below 3; sliding window of 4:15. Using Bowtie, the remaining high-quality reads were mapped onto the *Cucumis melo* genome v3.5.1 with maximum mismatch of 1 bp. Unique mapping of reads was adopted. To determine the target regions of H3K9Ac ChIP-Seq, the model-based analysis of ChIP-Seq (MACS2) was adopted (Number of duplicate reads at a location: 1; Bandwidth: 300; mfold of 5:30; q-value cutoff: 0.05) [72]. For peak detection in H3K27me3 modified regions, we used SICER with the following criteria: Window size: 200, Gap size: 600 [98]. Alignment and tag-density were inspected with IGB. HOMER was used to associate peaks to nearby genes [99]. To cluster the H3K9Ac and H3K27me3 peaks, linear normalization and clustering of tag density with Density Array method (window size 50 bp; 2 kb flanking region of genes) was performed using SeqMINER [100]. Annotation of corresponding genes was done using melonomics resource (https://melonomics.net). Average profile of coverage along the genic region (between transcriptional start sites (TSS) and TES) along with the 2 kb flanking region was plotted using NGSplot in binning mode [101]. To identify differentially enriched histone modified sites in leaves, root, fruit and flowers, we used DiffReps with the following settings: Window size 1 kb; step size: 100 bp; *p* value: 0.0001; Statistical testing method: Chi-square method) [74].

## Chromatin immunoprecipitation of Laser capture microdissected primordia

Male and female flowers of 3 mm were collected on monoecious melon plants, fixed in 1% (v/v) formaldehyde for 15 min under vacuum, and the crosslink reaction was stopped by addition of glycine (130 mM final) for 5 min under vacuum. Fixed material was then rinsed with water, dried and snap-frozen in liquid nitrogen. In a cryostat, flowers were put and oriented at the bottom of an empty cryomold and covered by cryoprotector (Tissue-Tek® O.C.T. Compound, Sakura® Finetek, VWR, France) until frozen. After 10 min at -20 °C in the cryostat (Shandon Cryotome® FSE, Excilone, France), samples were sectioned at -20 °C. To locate the targeted structures, sections were collected on glass slide and observed through the optical microscope of an XT® Arcturus Technologies microdissection system (Excilone, France). The following 4 sections (30 µm thickness) containing 10 flowers were collected on Arcturus FRAME membrane (Excilone, France). Then, the frame was put to room temperature for the sections to melt on the membrane and a glass slide was added on top of the membrane. The whole system is then placed upside down in the microdissection system. The LCM process was carried out using the XT® Arcturus Technologies

microdissection system and software. Capture was performed using Arcturus CapSure® LCM Macro Caps (Excilone, France). Carpels and stamens were identified using the 20X. FRAME Membrane Slides allowed the use of standard IR and UV lasers. UV laser was used to cut around the selected regions as IR laser was used to capture the selected regions on the macro caps without contaminating the sample with non-target-tissue. Efficiency of microdissection was evaluated by examining the CapSure® Macro caps after capture and the tissue section remaining on the slide before and after lifting off the CapSure® macro caps. To obtain enough material for downstream analysis, around 5 mm² of captured area was collected on 3 or 4 caps per sample and stored at -80 °C. Tissues were then disrupted by grinding with mortar and pestle in liquid nitrogen. Cellular debris was then resuspended in 2 ml of nuclei isolation buffer (20 mM Hepes pH8; 250 mM sucrose; 1 mM MgCl2; 5 mM KCl; 40% Glycerol; 0,25% Triton X-100), filtered through a strainer with a 63 µm pore size. Nuclei integrity was assessed by Dapi coloration and visualization under an epifluorescence microscope. Next, the samples were centrifuged 10 min at 3500g and the pellet was resuspended in 130 µl of nuclei lysis buffer (50 mM Tris-HCl pH8; 10 mM EDTA). Chromatin was sonicated using Covaris (M220 Focused-ultrasonicator; DC 5%; PIP 75Watts; CPB 200; 20 min) in 130 µl microtubes (Covaris AFA tubes) and centrifuged 10 min at 16,000g. Supernatant was then diluted ten times with ChIP dilution buffer (1,1% Triton; 1,2 mM EDTA; 16,7 mM Tris-HCl pH8; 167 mM NaCl) before immunoprecipitation with H3K9ac or H3K27me3 antibodies overnight at 4 °C with gentle shaking, and incubated for 1 h at 4 °C with 50 µL of Dynabeads Protein A (Invitrogen, Ref. 100-02D). After several washes, immunoprecipitated DNA was then recovered using proteinase K digestion, reverse crosslink and phenol-chloroform extraction as previously described. Fold enrichment of targets in ChIPed DNA relative to input was calculated from an average of three replicate qPCRs. The sequences of primers can be found in Additional file 8: Table VII. Positions of the amplified regions on the different loci are indicated in Fig. 4.

## Additional files

**Additional file 1: Figure S1.** Pearson correlation between leaves and roots ChIP-seq replicates. The color scale represents the degree of correlation between replicates. **Figure S2.** Example of melon genomic regions showing ChIP-seq signals of biological replicates in leaves. Genes are represented in blue. Comparisons of peak positions and peak intensities illustrate the high correlation between the two biological replicates. **Figure S3.** Distribution of mapped reads for H3K27me3 (red shades) and H3K9ac (green shades) along the 12 melon chromosomes of leaves, roots, fruit, male and female flowers. Local peak densities of each epigenetic mark were plotted against the genetic distance (gray) and annotation of sense
(dark blue) and antisense transcripts (light blue). **Figure S4.** Number of H3K9ac and H3K27me3 peaks identified in the different melon tissues. The computational methods MACS2 and SICER were used to determine H3K9ac and H3K27me3 target regions, respectively. **Figure S5.** H3K27me3 and H3K9ac distribution at the chromosome level in different melon organs. Distribution of H3K9ac (Green), H3K27me3 (red) and annotated genes (gray) are plotted along the chromosome 1. **Figure S6.** H3K9ac and H3K27me3 peaks distribution on chromosome 3 illustrating the enrichment of the marks in gene-rich regions located at the distal part of the chromosome.

**Figure S7.** Boxplot showing differential peak length between H3K9ac and H3K27me3 in different melon tissues. **Figure S8.** Average tag density profile of H3K27me3 and H3K9ac along the gene body in different melon tissues. ChIP-Seq densities of equal bins were plotted along the gene body and 2-kb region flanking the TSS or the TES. **Figure S9-S12.** Tissue-specific H3K9 or H3K27me3 differentially marked genes. Venn diagrams of paired comparisons of each tissue vs. all the tissues are shown on the left. The central overlap corresponds to H3K9 or H3K27me3 specifically modified genes in each tissue. Gene ontology analysis of these subsets of genes is shown on the right. **Figure S13.** Correlation of hyper- H3K9ac or hypo-H3K27me3 with tissue-specific gene activation. Venn diagram illustrating the genes whose high and tissue-specific expression is correlated with their tissue-specific H3K9 hyper-acetylation or H3K27me3 hypo-methylation (Overlapping set of genes). A Chi-squared test was done for each comparison and confirmed the observed correlation. **Figure S14.** Molecular function analysis of H3K9ac or H3K27me3 differentially regulated genes in male and female flowers. The most significant enrichment is observed for the transcription factor activity category. **Figure S15-S18.** Molecular function and biological process classification of H3K9ac and H3K27me3 differentially regulated genes in male and female flowers. Bar charts represent the number of genes in each category for the molecular function (left) and the biological process (right). **Figure S19.** ChIP-qPCR assays on microdissected stamens and carpels of a genome region that do not display any H3K9ac or H3K27me3 signals in the ChIP-seq data. **Figure S20.** Analysis of H3K27me3 deposition in *a* and *g* mutant flowers. **A.** Schematic representation of the different type of flowers used for H3K27me3 profiling in wild-type plants (AA-GG), *wip1* mutant (AA-gg) and *acs11* mutant (aa-GG) by ChIP-seq assays. **B.** Comparisons of differential H3K27me3 deposition in wild-type male (light blue) and female flowers (red) with differential H3K27me3 deposition in wild-type male (light blue) and *g* mutant female flowers in the same position (orange). **C.** Comparisons of differential H3K27me3 deposition in wild-type male (light blue) and female flowers (red) with differential H3K27me3 deposition in wild-type female flowers (red) and *a* mutant male flowers in the same position (dark blue). **Figure S21.** Proportion of H3K27me3 changes between male and female wild-type flowers controlled by the G or A locus. Hypo- and hyper-methylated genes were analyzed separately. **Figure S22.** GO analysis of G- and A-H3K27me3 dependent genes. **Figure S23.** Browser view of the H3K27me3 ChIP-seq signals on the G locus in Wt flowers and *a* mutant flowers. Only exon 2 of the *G locus* is present in the pseudomolecule of the 12 chromosomes in the 3.5.1 melon genome version. Promoter and exon 1 of the G locus are present on another pseudomolecule that contains non-anchored BAC clones (identified as chromosome 0 in the 3.5.1 melon genome version).

**Additional file 2: Table I.** Number of aligned reads for each library.

**Additional file 3: Table II.** List of H3K9ac and H3K27me3 marked genes in leaves, root, fruit, male and female flowers.

**Additional file 4: Table III.** List of H3K9ac and H3K27me3 specifically enriched or depleted genes in leaves, root, fruit, male and female flowers.

**Additional file 5: Table IV.** List of specifically over-expressed genes in leaves, root, fruit, male and female flowers.

**Additional file 6: Table V.** List of H3K9ac and H3K27me3 regulated genes in flowers.

**Additional file 7: Table VI.** List of transcription factor epigenetically regulated genes in flowers.

**Additional file 8: Table VII.** List of primers used in this study.

**Additional file 9: Table VIII.** *p*-values list of GO analyses.

## Authors' contributions

DL and NYRG performed the wet laboratory experiments, and AV, KMG and VS performed the bioinformatics analyses. CC, CD, DL, NYRG, and ABo prepared material. CB, NC, DL and NYRG were in charge of laser micro-dissection. DL, CR, NYRG, ABo, AB and MB analyzed the data. AB and MB designed the research. NYRG, DL, MB and AB wrote the manuscript. All authors read and approved the final manuscript.

## Author details

¹ Institute of Plant Sciences Paris-Saclay (IPS2), CNRS, INRA, University Paris-Sud, University of Evry, University Paris-Diderot, Sorbonne Paris-Cite, University of Paris-Saclay, Batiment 630, 91405 Orsay, France. ² Division of Biological and Environmental Sciences and Engineering, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia. ³ UMR 1313 Génétique Animale et Biologie Intégrative, Institut National de la Recherche Agronomique, 78350 Jouy-en-Josas, France. ⁴ UR 1052, Unité de Génétique et d'Amélioration des Fruits et Légumes, INRA, BP94, 84143 Montfavet, France.

## Competing interests

The authors declare that they have no competing interests.

## Availability of data and materials

The datasets generated and/or analyzed during the current study are available online GSE98054.

## Funding

This work was supported by the Plant Biology and Breeding department in INRA, the grants Program Saclay Plant Sciences (SPS, ANR-10-LABX-40), Initiative d'Excellence Paris-Saclay (Lidex-3P, ANR-11-IDEX-0003-02), L'Agence Nationale de la Recherche MELODY (ANR-11-BSV7-0024), and the European Research Council (ERC-SEXYPARTH).

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 2 March 2017 Accepted: 27 April 2017

Published online: 06 June 2017

## References

1. Nerson H. Seed production and germinability of cucurbit crops. Seed Sci Biotechnol. 2007;1:1–10.

2. Dahmani-Mardas F, Troadec C, Boualem A, Lévêque S, Alsadon AA, Aldoss AA, et al. Engineering melon plants with improved fruit shelf life using the TILLING approach. PLoS ONE. 2010. doi:10.1371/journal.pone.0015776.

3. Berger SL. Histone modifications in transcriptional regulation. Curr Opin Genet Dev. 2002;12:142–8.

4. Margueron R, Trojer P, Reinberg D. The key to development: interpreting the histone code? Curr Opin Genet Dev. 2005;15:163–76.

5. Zhou VW, Goren A, Bernstein BE. Charting histone modifications and the functional organization of mammalian genomes. Nat Rev Genet. 2011;12:7–18.

6. Rodriguez-Granados NY, Ramirez-Prado JS, Veluchamy A, Latrasse D, Raynaud C, Crespi M, et al. Put your 3D glasses on: plant chromatin is on show. J Exp Bot. 2016;67:3205–21.

7. Ramirez-Prado JS, Rodriguez-Granados NY, Ariel F, Raynaud C, Benhamed M. Chromatin architecture: a new dimension in the dynamic control of gene expression. Plant Signal Behav. 2016. doi:10.1080/15592 324.2016.1232224.

8. Lauria M, Rossi V. Epigenetic control of gene regulation in plants. Biochim Biophys Acta. 2011;1809:369–78.

9. Slotkin RK. Plant epigenetics: from genotype to phenotype and back again. Genome Biol. 2016. doi:10.1186/s13059-016-0920-5.

10. Goldberg AD, Allis CD, Bernstein E. Epigenetics: a landscape takes shape. Cell. 2007;128:635–8.
11. Jaenisch R, Bird A. Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. Nat Genet. 2003;33:245–54.

12. Schones DE, Cui K, Cuddapah S. Genome-wide approaches to studying yeast chromatin modifications. Methods Mol Biol. 2011;759:61–71.

13. Wang Z, Zang C, Rosenfeld JA, Schones DE, Cuddapah S, Cui K, et al. Combinatorial patterns of histone acetylations and methylations in the human genome. 2009;40:897–903.

14. Roudier F, Ahmed I, Bérard C, Sarazin A, Mary-Huard T, Cortijo S, et al. Integrative epigenomic mapping defines four main chromatin states in Arabidopsis. EMBO J. 2011;30:1928–38.

15. Malapeira J, Khaitova LC, Mas P. Ordered changes in histone modifications at the core of the Arabidopsis circadian clock. Proc Natl Acad Sci USA. 2012;109:21540–5.

16. Sequeira-Mendes J, Aragüez I, Peiró R, Mendez-Giraldez R, Zhang X, Jacobsen SE, et al. The functional topography of the Arabidopsis genome is organized in a reduced number of linear motifs of chromatin states. Plant Cell. 2014;26:2351–66.

17. Veluchamy A, Jégu T, Ariel F, Latrasse D, Mariappan KG, Kim SK, et al. LHP1 regulates H3K27me3 spreading and shapes the three-dimensional conformation of the Arabidopsis genome. PLoS ONE. 2016;11:1–25.

18. Shi J, Dawe RK. Partitioning of the maize epigenome by the number of methyl groups on histone H3 lysines 9 and 27. Genetics. 2006;173:1571–83.

19. Wang X, Elling AA, Li X, Li N, Peng Z, He G, et al. Genome-wide and organ-specific landscapes of epigenetic modifications and their relationships to mRNA and small RNA transcriptomes in maize. Plant Cell. 2009;21:1053–69.

20. Zhang W, Garcia N, Feng Y, Zhao H, Messing J. Genome-wide histone acetylation correlates with active transcription in maize. Genomics. 2015;106:214–20.

21. Li X, Wang X, He K, Ma Y, Su N, He H, et al. High-resolution mapping of epigenetic modifications of the rice genome uncovers interplay between DNA methylation, histone methylation, and gene expression. Plant Cell. 2008;20:259–76.

22. He G, Zhu X, Elling AA, Chen L, Wang X, Guo L, et al. Global epigenetic and transcriptional trends among two rice subspecies and their reciprocal hybrids. Plant Cell. 2010;22:17–33.

23. Du Z, Li H, Wei Q, Zhao X, Wang C, Zhu Q, et al. Genome-wide analysis of histone modifications: H3K4me2, H3K4me3, H3K9ac, and H3K27ac in *Oryza sativa* L. J Mol Plant. 2013;6:1463–72.

24. Ayyappan V, Kalavacharla V, Thimmapuram J, Bhide KP, Sripathi VR, Smolinski TG, et al. Genome-wide profiling of histone modifications (H3K9me2 and H4K12ac) and gene expression in rust (Uromyces appendiculatus) inoculated common bean (*Phaseolus vulgaris* L.). PLoS ONE. 2015. doi:10.1371/journal.pone.0132176.

25. Baker K, Dhillon T, Colas I, Cook N, Milne I, Milne L, et al. Chromatin state analysis of the barley epigenome reveals a higher-order structure defined by H3K27me1 and H3K27me3 abundance. Plant J. 2015;84:111–24.

26. Barski A, Cuddapah S, Cui K, Roh TY, Schones DE, Wang Z, et al. High-resolution profiling of histone methylations in the human genome. Cell. 2007;129:823–37.

27. Mikkelsen TS, Ku M, Jaffe DB, Issac B, Lieberman E, Giannoukos G, et al. Genome-wide maps of chromatin state in pluripotent and lineage-committed cells. Nature. 2007;448:553–60.

28. Bannister AJ, Kouzarides T. Regulation of chromatin by histone modifications. Cell Res. 2011;21:381–95.

29. Zhang X, Clarenz O, Cokus S, Bernatavichute YV, Pellegrini M, Goodrich J, et al. Whole-genome analysis of histone H3 lysine 27 trimethylation in Arabidopsis. PLoS Biol. 2007. doi:10.1371/journal.pbio.0050129.

30. Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, Hawkins RD, et al. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet. 2007;39:311–8.

31. Zhou J, Wang X, He K, Charron J-B, Elling A, Deng X. Genome-wide profiling of histone H3 lysine 9 acetylation and dimethylation in Arabidopsis reveals correlation between multiple histone marks and gene expression. Plant Mol Biol. 2010;72:585–95.

32. Schubert D, Primavesi L, Bishop A, Roberts G, Doonan J, Jenuwein T, et al. Silencing by plant Polycomb-group genes requires dispersed trimethylation of histone H3 at lysine 27. EMBO J. 2006;25:4638–49.

33. Turck F, Roudier F, Farrona S. Martin-Magniette ML, Guillaume E, Buisine N, et al. Arabidopsis TFL2/LHP1 specifically associates with genes marked by trimethylation of histone H3 lysine 27. PLoS Genet. 2007. doi:10.1371/journal.pgen.0030086.

34. Chouard P. Vernalization and its relations to dormancy. Annu Rev Plant Physiol. 1960;11:191–238.

35. Sheldon CC, Rouse DT, Finnegan EJ, Peacock WJ, Dennis ES. The molecular basis of vernalization: the central role of FLOWERING LOCUS C (FLC). Proc Natl Acad Sci. 2000;97:3753–8.

36. Amasino R. Vernalization, competence, and the epigenetic memory of winter. Plant Cell. 2004;16:2553–9.

37. Dennis ES, Peacock WJ. Epigenetic regulation of flowering. Curr Opin Plant Biol. 2007;10:520–7.

38. De Lucia F, Crevillen P, Jones AME, Greb T, Dean C. A PHD-polycomb repressive complex 2 triggers the epigenetic silencing of FLC during vernalization. Proc Natl Acad Sci. 2008;105:16831–6.

39. Crevillén P, Sonmez C, Wu Z, Dean C. A gene loop containing the floral repressor FLC is disrupted in the early phase of vernalization. EMBO J. 2012;32:140–8.

40. Lee JT. Lessons from X-chromosome inactivation: long ncRNA as guides and tethers to the epigenome. Genes Dev. 2009;23:1831–42.

41. letswaart R, Wu Z, Dean C. Flowering time control: another window to the connection between antisense RNA and chromatin. Trends Genet. 2012;28:445–53.

42. Kinoshita T, Seki M. Epigenetic memory for stress response and adaptation in plants. Plant Cell Physiol. 2014;55:1859–63.

43. Yamuroc U, Zhu J-K, Yang Z. Epigenetic modifications and plant hormone action. Mol Plant. 2016;9:57–70.

44. Álvarez-Venezas R, De-la-Peña C. Editorial: recent advances of epigenetics in crop biotechnology. Front Plant Sci. 2016. doi:10.3389/fpls.2016.000413.

45. Renner SS. The relative and absolute frequencies of angiosperm sexual systems: dioecy, monoecy, gynodioecy, and an updated online database. Am J Bot. 2014;101:1588–96.

46. Ainsworth C, Parker J, Buchananwollaston V. Sex determination in plants. Curr Top Dev Biol. 1997;38:167–223.

47. Martin A, Troadec C, Boualem A, Rajab M, Fernandez R, Morin H, et al. A transposon-induced epigenetic change leads to sex determination in melon. Nature. 2009;461:1135–8.

48. Sebastian P, Schaefer H, Telford IRH, Renner SS. Cucumber (*Cucumis sativus*) and melon (C. melo) have numerous wild relatives in Asia and Australia, and the sister species of melon is from Australia. Proc Natl Acad Sci USA. 2010;107:14269–73.

49. Boualem A, Fergany M, Fernandez R, Troadec C, Martin A, Morin H, et al. A conserved mutation in an ethylene biosynthesis enzyme leads to andromonoecy in melons. Science. 2008;321:836–8.

50. Boualem A, Troadec C, Kovalski I, Sari M-A, Perl-Treves R, Bendahmane A. A conserved ethylene biosynthesis enzyme leads to andromonoecy in two Cucumis species. PLoS ONE. 2009. doi:10.1371/journal.pone.0006144.

51. Boualem A, Troadec C, Camps C, Lemhemdi A, Morin H, Sari M-A, et al. A cucurbit androecy gene reveals how unisexual flowers develop and dioecy emerges. Science. 2015;350:688–91.

52. Bai SL, Ben Peng Y, Cui JX, Gu HT, Xu LY, Li YQ, et al. Developmental analyses reveal early arrests of the spore-bearing parts of reproductive organs in unisexual flowers of cucumber (*Cucumis sativus* L.). Planta. 2004;220:230–40.

53. Ming R, Bendahmane A, Renner SS. Sex chromosomes in land plants. Annu Rev Plant Biol. 2011;62:485–514.

54. Werren JH, Beukeboom LW. Sex determination, sex ratios, and genetic conflict. Annu Rev Ecol Syst. 1998;29:233–61.

55. Bisoni L, Battle-Morera L, Bird AP, Suzuki M, McQueen HA. Female-specific hyperacetylation of histone H4 in the chicken Z chromosome. Chromosom Res. 2005;13:205–14.

56. Angelopoulou R, Lavranos G, Manolakou P. Regulatory RNAs and chromatin modification in dosage compensation: a continuous path from flies to humans? Reprod Biol Endocrinol. 2008. doi:10.1186/1477-7827-6-12.

57. Piferrer F. Epigenetics of sex determination and gonadogenesis. Dev Dyn. 2013;242:360–70.

58. Talamali A, Bajji M, Le Thomas A, Kinet J-M, Dutuit P. Flower architecture and sex determination: how does *Atriplex halimus* play with floral morphogenesis and sex genes? New Phytol. 2003;157:105–13.
59. Freeman DC, Vitale JJ. The influence of environment on the sex ratio and fitness of spinach. Bot Gaz. 1985;146:137-42.

60. Eveland AL, Satoh-Nagasawa N, Goldshmidt A, Meyer S, Beatty M, Sakai H, et al. Digital gene expression signatures for maize development. Plant Physiol. 2010;154:1024-39.

61. Rudich J, Halevy AH, Kedar N. Ethylene evolution from cucumber plants as related to sex expression. Plant Physiol. 1972;49:998-9.

62. Sun J-J, Li F, Li X, Liu X-C, Rao G-Y, Luo J-C, et al. Why is ethylene involved in selective promotion of female flower development in cucumber? Plant Signal Behav. 2010;5:1052-6.

63. Liu H, Yang X, Liao X, Zuo T, Qin C, Cao S, et al. Genome-wide comparative analysis of digital gene expression tag profiles during maize ear development. Genomics. 2015;106:52-60.

64. Rocheta M, Sobral R, Magalhães J, Amorim MI, Ribeiro T, Pinheiro M, et al. Comparative transcriptomic analysis of male and female flowers of monoecious *Quercus suber*. Front Plant Sci. 2014. doi:10.3389/fpls.2014.00599.

65. Cheng H, Qin L, Lee S, Fu X, Richards DE, Cao D, et al. Gibberellin regulates Arabidopsis floral development via suppression of DELLA protein function. Development. 2004;131:1055-64.

66. Plackett ARG, Thomas SG, Wilson ZA, Hedden P. Gibberellin control of stamen development: a fertile field. Trends Plant Sci. 2016;16:568-78.

67. Zhang Y, Liu B, Yang S, An J, Chen C, Zhang X, et al. A cucumber DELLA homolog CsGAIP May inhibit staminate development through transcriptional repression of B class floral homeotic genes. PLoS ONE. 2014. doi:10.1371/journal.pone.0091804.

68. Janoušek B, Šíroký J, Vyskot B. Epigenetic control of sexual phenotype in a dioecious plant, *Melandrium album*. Mol Gen Genet. 1996;250:483-90.

69. Parkinson SE, Gross SM, Hollick JB. Maize sex determination and abaxial leaf fates are canalized by a factor that maintains repressed epigenetic states. Dev Biol. 2007;308:462-73.

70. Akagi T, Henry IM, Tao R, Comai L. A Y-chromosome-encoded small RNA acts as a sex determinant in persimmons. Science. 2014;346:646-50.

71. Akagi T, Henry IM, Kawai T, Comai L, Tao R. Epigenetic Regulation of the Sex Determination Gene MeG1 in Polyploid Persimmon. Plant Cell. 2016. doi:10.1105/tpc.16.00532.

72. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008. doi:10.1186/gb-2008-9-9-r137.

73. Garcia-Mas J, Benjak A, Sanseverino W, Bourgeois M, Mir G, González VM, et al. The genome of melon (*Cucumis melo* L.). Proc Natl Acad Sci. 2012;109:11872-7.

74. Shen L, Shao NY, Liu X, Maze I, Feng J, Nestler EJ. diffReps: detecting differential chromatin modification sites from ChIP-seq data with biological replicates. PLoS ONE. 2013. doi:10.1371/journal.pone.0065598.

75. Joung J-G, Corbett AM, Fellman SM, Tieman DM, Klee HJ, Giovannoni JJ, et al. Plant MetGenMAP: an integrative analysis system for plant systems biology. Plant Physiol. 2009;151:1758-68.

76. Weigel D, Meyerowitz EM. The ABCs of floral homeotic genes. Cell. 1994;78:203-9.

77. Krizek BA, Fletcher JC. Molecular mechanisms of flower development: an armchair guide. Nat Rev Genet. 2005;6:688-98.

78. Guo S, Sun B, Looi L-S, Xu Y, Gan E-S, Huang J, et al. Co-ordination of flower development through epigenetic regulation in two model species: rice and Arabidopsis. Plant Cell Physiol. 2015;56:830-42.

79. Espina V, Wulfkuhle JD, Calvert VS, Van Meter A, Zhou W, Coukos G, et al. Laser-capture microdissection. Nat Protoc. 2006;1:586-603.

80. Wang DH, Li F, Duan QH, Han T, Xu ZH, Bai SN. Ethylene perception is involved in female cucumber flower development. Plant J. 2010;61:862-72.

81. Manzano S, Martínez C, García JM, Megías Z, Jamilena M. Involvement of ethylene in sex expression and female flower development in water-melon (*Citrullus lanatus*). Plant Physiol Biochem. 2014;85:96-104.

82. Byers RE, Baker LR, Sell HM, Herner RC, Dilley DR. Ethylene: a natural regulator of sex expression of *Cucumis melo* L. Proc Natl Acad Sci USA. 1972;69:717-20.

83. Ando S, Sato Y, Kamachi S, Sakai S. Isolation of a MADS-box gene (ERAF17) and correlation of its expression with the induction of formation of female flowers by ethylene in cucumber plants (*Cucumis sativus* L.). Planta. 2001;213:943-52.

84. Acosta IF, Laparra H, Romero SP, Schmelz E, Hamberg M, Mottinger JP, et al. Tasselseed1 is a lipoxygenase affecting jasmonic acid signaling in sex determination of maize. Science. 2009;323:262-5.

85. Wolters-Arts M, Lush WM, Mariani C. Lipids are required for directional pollen-tube growth. Nature. 1998;392:818-21.

86. Riechmann JL, Heard J, Martin G, Reuber L, Jiang C, Keddie J, et al. Arabidopsis transcription factors: genome-wide comparative analysis among eukaryotes. Science. 2000;290:2105-10.

87. Heim MA, Jakoby M, Werber M, Martin C, Weisshaar B, Bailey PC. The basic helix-loop-helix transcription factor family in plants: a genome-wide study of protein structure and functional diversity. Mol Biol Evol. 2003;20:735-47.

88. Wang D, Guo Y, Wu C, Yang G, Li Y, Zheng C. Genome-wide analysis of CCCH zinc finger family in Arabidopsis and rice. BMC Genom. 2008. doi:10.1186/1471-2164-9-44.

89. Zhang F, Qi B, Wang L, Zhao B, Rode S, Riggan ND, et al. EIN2-dependent regulation of acetylation of histone H3K14 and non-canonical histone H3K23 in ethylene signalling. Nat Commun. 2016. doi:10.1038/ncomms13018.

90. Lu F, Cui X, Zhang S, Jenuwein T, Cao X. Arabidopsis REF6 is a histone H3 lysine 27 demethylase. Nat Genet. 2011;43:715-9.

91. Smaczniak C, Immink RGH, Muiño JM, Blanvillain R, Busscher M, Busscher-Lange J, et al. Characterization of MADS-domain transcription factor complexes in Arabidopsis flower development. Proc Natl Acad Sci. 2012;109:1560-5.

92. Goodrich J, Puangsomlee P, Martin M, Long D, Meyerowitz EM, Coup-land G. A Polycomb-group gene regulates homeotic gene expression in Arabidopsis. Nature. 1997;386:44-51.

93. Krizek BA, Lewis MW, Fletcher JC. RABBIT EARS is a second-whorl repressor of AGAMOUS that maintains spatial boundaries in Arabidopsis flowers. Plant J. 2006;45:369-83.

94. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114-20.

95. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25:1105-11.

96. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511-5.

97. Goff LA, Trapnell C, Kelley D. CummeRbund: visualization and exploration of Cufflinks high-throughput sequencing data. R Package. Version 2.2. 2012.

98. Zang C, Schones DE, Zeng C, Cui K, Zhao K, Peng W. A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics. 2009;25:1952-8.

99. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576-89.

100. Ye T, Krebs AR, Choukrallah MA, Keime C, Plewniak F, Davidson I, et al. seqMINER: an integrated ChIP-seq data interpretation platform. Nucleic Acids Res. 2011;39:1-10.

101. Shen L, Shao N, Liu X, Nestler E. ngs plot: quick mining and visualization of next-generation sequencing data by integrating genomic databases. BMC Genom. 2014. doi:10.1186/1471-2164-15-284.