1. Plant Biology

The petal identity gene PhDEF has a major binding and regulatory action in the epidermis

  1. Emma Désert
  2. Quentin Cavallini-Speisser
  3. Evelyne Duvernois-Berthet
  4. Pierre Chambrier
  5. Patrice Morel
  6. Brice Letcher
  7. Carine Rey
  8. Jeremy Just
  9. Suzanne Rodrigues Bento
  10. Daniel Bouyer  Is a corresponding author
  11. Marie Monniaux  Is a corresponding author
  1. Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, France
  2. Laboratoire de Biologie et Modélisation de la Cellule, École Normale Supérieure de Lyon, CNRS, UMR 5239, Inserm, U1293, Université Claude Bernard Lyon 1, France
  3. CIRI, Centre International de Recherche en Infectiologie, Univ Lyon, Inserm, U1111, Université Claude Bernard Lyon 1, CNRS, UMR5308, ENS de Lyon, France
https://doi.org/10.7554/eLife.113154.2
Share this article
Cite this article
  1. Emma Désert
  2. Quentin Cavallini-Speisser
  3. Evelyne Duvernois-Berthet
  4. Pierre Chambrier
  5. Patrice Morel
  6. Brice Letcher
  7. Carine Rey
  8. Jeremy Just
  9. Suzanne Rodrigues Bento
  10. Daniel Bouyer
  11. Marie Monniaux
(2026)
The petal identity gene PhDEF has a major binding and regulatory action in the epidermis
eLife 15:RP113154.
https://doi.org/10.7554/eLife.113154.2
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

This valuable study investigates the action of a floral homeotic transcription factor in different cell types. The results provide solid evidence that PhDEF shows more extensive binding and regulatory effects in the epidermis than in the mesophyll, supporting the conclusion that cell-layer context influences homeotic gene function. While there are some limitations, they do not undermine the principal observation of layer-dependent PhDEF activity, but the interpretation of the data is complicated by uncertainties in cell-type assignment and composition in the scRNA-seq data, the cell-type-enriched rather than strictly cell-type-specific nature of the ChIP-seq data, and differences in developmental stage between the single-cell and ChIP-seq experiments. The work will be of interest to developmental biologists studying transcriptional regulation and cell identity.

[Editors' note: this paper was reviewed by Review Commons.]

https://doi.org/10.7554/eLife.113154.2.sa0

Abstract

In flowering plants, floral organ identity is specified by the combinatorial action of homeotic genes. While the role of these genes in the early specification of organ identity is well established, their late function throughout floral organ development and in specific cell types is much less characterized. In particular, since plant organs are structured in clonally independent cell layers, whether and how homeotic identity interacts with cell layer identity is unknown. We have previously identified cell layer-specific mutants for the petal identity gene PhDEF in petunia flowers, resulting in drastically different petal phenotypes whether PhDEF is expressed in the petal epidermis or in the mesophyll. In this study, using a combination of single-cell RNA-Seq and chromatin immunoprecipitation on phdef cell layer-specific mutants, we find that PhDEF regulates a different set of target genes in the petal epidermis and mesophyll, with a major regulatory action in the epidermis. We uncover a high diversity of binding profiles in PhDEF target genes, with a complex combination of layer-specific or non-specific binding sites, and a much more prominent binding of PhDEF in the epidermis than in the mesophyll. Our study highlights that floral homeotic genes like PhDEF can have different regulatory actions in different cell contexts, here different cell layers, demonstrating that cell layer identity indeed influences the regulatory processes underlying homeotic identity.

Introduction

In flowering plants, homeotic genes define the identity of floral organs in a combinatorial manner, as outlined in the ABC model (Schwarz-Sommer et al., 1990; Coen and Meyerowitz, 1991; Bowman and Moyroud, 2024). These genes encode master transcription factors, most of them from the MADS-box family, that regulate the expression of a wide set of genes throughout floral organ development until organs are mature. Additionally, floral organs are structured in different cell layers that are inherited from the three clonally distinct layers (L1, L2, and L3) of the shoot apical meristem (Satina et al., 1940; Stewart and Burk, 1970; Meyerowitz, 1997; Scheres, 2001). The action of homeotic genes is therefore superimposed on a pre-established cell layer identity, raising the question whether homeotic function is dependent or influenced by cell layer identity.

Petals are often key players in pollinator attraction, and in the case of our focus genus Petunia, this is typified by the bright pigmentation or strong scent emission from the conical cells of the limb epidermis, while the long tube observed in some petunia species limits nectar access to adapted pollinators (Galliot et al., 2006; Rodrigues et al., 2018). In the model organism Petunia x hybrida (petunia), petal identity is defined by B-class MADS-box regulators only: PhDEF (Petunia x hybrida DEFICIENS), PhGLO1 (GLOBOSA1), and PhGLO2 (Morel et al., 2017; Vandenbussche et al., 2004; Rijpkema et al., 2006). The encoded proteins form PhDEF/PhGLO1 and PhDEF/PhGLO2 obligate heterodimers that translocate to the nucleus, bind DNA, and regulate gene transcription, as was reported for their homologs in Antirrhinum and Arabidopsis (Tröbner et al., 1992; Schwarz-Sommer et al., 1992; McGonigle et al., 1996). These complexes activate the expression of the PhDEF, PhGLO1, and PhGLO2 genes themselves, resulting in an auto-activation loop (Tröbner et al., 1992; Schwarz-Sommer et al., 1992; Vandenbussche et al., 2004). This loop, that is conserved in all eudicots, buffers stochastic noise, leading to a sharp transition to petal identity, and ensures the maintenance of high expression levels of these genes until petals are mature (Winter et al., 2002; Lenser et al., 2009). Other proteins, in particular the E-class SEPALLATA MADS-box proteins (Pelaz et al., 2000; Pelaz et al., 2001), are part of higher-order complexes with the B-class heterodimers to define petal identity (Immink et al., 2003; Ferrario et al., 2003; Morel et al., 2019). B-class genes, after defining petal identity, also participate in their development and in the differentiation of their cell types. Indeed, experiments of temporal inactivation of B-class genes in Arabidopsis thaliana have revealed that, at intermediate stages of flower development, petals easily revert to a sepal-like identity (Wuest et al., 2012), suggesting that expression of B-class genes needs to be maintained over the course of petal development for normal petal identity specification. Similarly, inhibition of B-class gene expression at later stages of Antirrhinum majus flower development causes a reduction in petal size, conical cell size, and scent emission (Manchado-Rojo et al., 2012), and a reduction in scent production and emission was also recently reported in petunia petals when PhDEF expression is repressed at anthesis (Bednarczyk et al., 2025). This suggests that petal identity remains partly uncommitted until intermediate-to-late stages of flower development, and that B-class genes act to maintain this late specification of identity through the control of specific cell fates and traits.

In this study, we investigate the late role of the petal identity regulator PhDEF in the two cell layers of the petunia petal, stemming from our previous identification of cell layer-specific phdef mutants (Chopy et al., 2024). In petunia, petals derive from the L1 and L2 meristematic layers that will form the petal epidermis and mesophyll, respectively (Satina and Blakeslee, 1941; Chopy et al., 2024). While the phdef mutant flowers display a homeotic conversion of petals into sepals, the layer-specific phdef mutation results in two contrasting phenotypes (Chopy et al., 2024; Figure 1A): flowers with a phdef mutant epidermis (that we named star) form a normal tube but a small, star-shaped, and unpigmented limb. Meanwhile, flowers with a phdef mutant mesophyll (which we named wico for wide corolla) have a much shorter tube than normal, but the limb is normally shaped and pigmented, although the petals are often pink rather than red. These flowers derive from the excision of the transposable element dTph1 (inserted in the first exon of the PhDEF gene and knocking-out its function) in a unique cell layer, as we previously described (Chopy et al., 2024). The star and wico phenotypes show that PhDEF plays distinct roles in the two layers of the petal and drives partly independently limb and tube development. This material allows for the identification of cell layer-specific regulatory processes driven by PhDEF without the need to resort to tissue dissection. Here, to further understand how PhDEF achieves these cell layer-specific roles, we identified the cell layer-enriched target genes of PhDEF by characterizing layer-specific transcriptional profiles obtained by single-cell RNA sequencing (scRNA-Seq) and layer-specific binding profiles obtained by chromatin immunoprecipitation sequencing (ChIP-Seq), in wild-type (WT), star, and wico petals. We found that, although PhDEF is similarly expressed in the two layers of the petal, PhDEF binds to and regulates the expression of many more target genes in the epidermis than in the mesophyll. The target genes of PhDEF display a combination of layer-specific and non-specific binding sites in their regulatory regions, suggesting a complex pattern of transcriptional control among cell layers. We also found differential enrichment of motifs under epidermal and mesophyll-specific peaks, suggestive of different protein partners for PhDEF in the two petal cell layers. This study offers a novel cell layer-specific perspective on transcriptional regulation and shows that homeotic regulators such as PhDEF are influenced by the pre-established cell layer identity.

Figure 1 with 7 supplements see all
Download asset Open asset
Single-cell RNA-sequencing separates epidermal and mesophyll cell types in the petunia petal.

(A) Representative petunia wild-type (WT), star, wico, and phdef flowers from an upper (top picture) and side view (bottom picture). The limb and tube are indicated. Scale bar = 1 cm. Below the flower pictures, a schematic petal primordium is depicted with L1 (future epidermis) and L2 (future mesophyll) cell layers, and PhDEF expression in the WT, star, wico, and phdef primordia is indicated in orange. (B) Overview of the experimental protocol to produce petal protoplasts ready for isolation with the 10 X Genomics Chromium device. Upper picture: a WT flower cut into small fragments in the cell wall-digesting solution (scale bar = 1 cm). Bottom picture: isolated protoplasts from a WT petal after a 5 hr-long digestion process, viewed by bright-light microscopy (scale bar = 50 µm). (C) Uniform Manifold Approximation and Projection (UMAP) plot of 11,632 WT petal cells sequenced for their transcriptome, after integration from two biological replicates. Clusters of cells are color-coded and are ordered from the biggest to the smallest, based on the number of cells. (D) Dotplot of key marker genes for petal identity, stamen and carpel identity, epidermal identity, and pigmentation. The size of the dot represents the percentage of cells that express a given gene, and the color scale indicates the average expression level of the gene across all cells in a given cluster. st., stamen; ca., carpel. (E) Barplot of the number of genes enriched in the epidermis (magenta) or in the mesophyll (green), as defined by varying the log2FoldChange (log2FC) threshold and the layer specificity threshold (% of cells expressing the gene in the epidermis − % of cells expressing the gene in the mesophyll). The Kolmogorov-Smirnov two-sided test (KS test) was applied to compare the gene number distributions between epidermis- and mesophyll-enriched genes at a given layer-specificity threshold; n.s., non-significant, *p<0.05, **p<0.005, ***p<0.001.

Results

Single-cell RNA sequencing reveals separate clusters of epidermal and mesophyll cell types in the petunia petal

The star and wico phenotypes (Figure 1A) suggest that PhDEF regulates a distinct set of target genes in the two cell layers, contributing to either tube or limb development. We first aimed to identify the genes that exhibited enriched or non-enriched expression in the two WT petal layers, and for this we sequenced the transcriptome of single cells from WT petunia petals. To this end, the petals were isolated by carefully removing the stamen filaments that are fused to the base of the tube, and we performed a 5 hr-long protoplast isolation of the whole corolla (tube and limb) by cell wall-digesting enzymes (Figure 1B). While we were originally aiming to isolate protoplasts from petals at an early developmental stage (stage 8 as defined in Chopy et al., 2024) to capture the dynamics of gene regulation during petal development, our protocol only allowed us to isolate a sufficient amount of protoplasts for scRNA-Seq in fully differentiated petals. A visual inspection of protoplasts under the microscope showed that the proportion of epidermal (pigmented with pink anthocyanins) and mesophyll (non-pigmented and/or chloroplastic) cells recovered was not significantly different from the proportion of epidermal and mesophyll cells observed in petal cross-sections (Figure 1—figure supplement 1), indicating that protoplasting allowed to efficiently recover cells from the two petal layers. Protoplasts were subjected to transcriptome sequencing with the 10 X Genomics Chromium technology (Figure 1—figure supplement 2), resulting in the transcriptome from 11,632 cells from WT mature petals in two biological replicates (Figure 1C, Pearson’s correlation coefficient r=0.96 between pseudo-bulk datasets of the two replicates).

The B-class petal identity genes PhDEF, PhGLO1, and PhGLO2 were found to be broadly expressed in all eight clusters (Figure 1D), while the stamen and carpel identity genes PhTM6, FBP6, and pMADS3 (Rijpkema et al., 2006; Heijmans et al., 2012) were hardly detected, indicating the absence of contamination by stamen tissue. In order to assign identity to each cluster, we looked at known marker genes from the literature, identified their putative homologs in petunia (Supplementary file 1), and observed their expression in each cluster (Figure 1D, Figure 1—figure supplement 3). We also explored Gene Ontology (GO) term enrichment for the best marker genes of each cluster (Figure 1—figure supplement 4, Supplementary file 1), and we determined the limb or tube origin of respective clusters by RT-qPCR on a chosen set of cluster markers (Figure 1—figure supplement 5). The specific expression of the sucrose transporter SWEET11, and several UmamiT amino acid transporters and sulfate transporters, clearly marked cluster 5 as a group of vascular cells (Kim et al., 2021; Figure 1—figure supplement 3C). Similarly, the strong expression of several anthocyanin biosynthesis genes allowed the identification of cluster 6 as pigmented cells, hence cells from the upper limb epidermis (Bombarely et al., 2016; Figure 1D). Other known epidermal markers, such as ketoacyl-CoA synthases (KCS, wax synthesis enzymes) (Hegebarth and Jetter, 2017) or members of the Homeodomain Glabrous (HDG) family (Nakamura et al., 2006), identified clusters 0 and 4 as two other epidermal clusters (Figure 1D). Cluster 0 likely groups cells from the lower and upper limb epidermis; indeed, these cells faintly express pigmentation genes (Figure 1D); marker genes from this cluster are strongly expressed in limb tissue (Figure 1—figure supplement 5A); and in situ hybridization for the cluster marker IAA14 indicated a clear expression in both upper and lower limb epidermis (Figure 1—figure supplement 5C). Consequently, this suggests that cells in cluster 6 represent only a fraction of the upper limb epidermal cells, that is the ones most strongly expressing pigmentation genes, perhaps the most distal cells from the limb that are the latest to differentiate. The highly specific expression of TERPENE SYNTHASE1 (PhTPS1), a tube-specific gene involved in volatile emission inside the flower (Boachon et al., 2019), suggested that cluster 4 contains cells from the tube epidermis (Figure 1—figure supplement 5A,B). The strong and specific expression of several histone genes and S-phase genes (Figure 1—figure supplement 3B, D) identified the smallest cluster (7) as a cell-cycle state cluster (i.e. dividing cells, independent of a specific cell identity), as often found in single-cell transcriptomics (Zhang et al., 2021; Chen et al., 2025). Finally, we assigned a mesophyll identity to the remaining clusters (1, 2, and 3), since they expressed genes strongly associated with photosynthetic activity for cluster 3 (Figure 1-figure supplement 3) and with water transport for cluster 2 (Figure 1—figure supplement 4, Supplementary file 2); they did not express any of the epidermal or pigmentation markers (Figure 1D), and these clusters contained cells from either tube or limb origin (Figure 1—figure supplement 5A). We verified that the petunia homologs to genes defined as epidermal or mesophyll-enriched in tobacco petals also displayed a similar enrichment in our clusters (Zhang et al., 2021; Figure 1—figure supplement 3E). Taken together, although there remains some uncertainty about the precise spatial location of each cluster, we are confident about their belonging to either the epidermis or the mesophyll layer.

Protoplasting is known to induce strong transcriptomic changes that could impact the clustering and the assignment of cell identities (Denyer and Timmermans, 2022). To estimate the effect of protoplasting on gene expression in our samples, we performed bulk RNA-Seq on protoplasted vs. intact WT petal tissue, and we observed a high correlation between the gene expression levels in the two conditions (r=0.84). Still, a high number of differentially expressed genes (DEGs) were identified (5571 and 5891 genes upregulated and downregulated by the protoplasting process, respectively, Supplementary file 3). We explored the effect of removing these DEGs on the clustering from the scRNA-Seq data (Figure 1—figure supplement 6): the Uniform Manifold Approximation and Projection (UMAP) plots exhibited a high degree of similarity in their shape irrespective of the removal of the DEGs, and the same clusters were recovered with only slight variation in the number of cells in each cluster. Therefore, we decided to keep all genes and all cells in our analysis.

In summary, our assignment of cluster identities revealed that the three mesophyll clusters are not clearly distinct in the UMAP space, in contrast to the three epidermal clusters that have distinct transcriptomic profiles. Importantly, there is a clear separation of epidermal and mesophyll cell types in the UMAP space, as previously observed in scRNA-Seq data from the petals of mature tobacco flowers (Kang et al., 2022), showing that these layers have distinct transcriptomic signatures overall. Since our main interest was in the layer-specific functions of PhDEF, we next merged the three epidermal and the three mesophyll clusters together, setting apart the vasculature.

The petal epidermis expresses more specific genes than the petal mesophyll

We first explored the difference in gene expression between WT petal layers. We computed gene differential expression between the epidermis and the mesophyll, and we examined the layer-enriched genes, that is with a higher expression level in one layer and with a high specificity of expression in one layer (higher percentage of cells expressing the gene in this layer). We noticed that there were consistently more genes whose expression is specifically enriched in the epidermis than in the mesophyll (Figure 1E). For instance, with log2FoldChange >0.75 and layer-specificity of expression >10% (see Materials and methods), there were 695 genes enriched in the epidermis compared to 151 genes enriched in the mesophyll (Supplementary file 4), and this asymmetry was observed at varying log2FoldChange and layer-specificity thresholds (Figure 1E). Some of the GO terms enriched in layer-enriched genes were consistent with the function of these layers (Figure 1—figure supplement 7, Supplementary file 2): the epidermis-enriched genes were associated with GO terms related to the metabolism of surface lipids or the phenylpropanoid pathway (an extensive network that leads to the production of lignin, anthocyanins, and flavonols), while the mesophyll-enriched genes were associated with several GO terms related to photosynthesis. Overall, this initial analysis of layer-enriched genes showed that there are more epidermis-enriched genes than mesophyll-enriched genes in the WT petal.

PhDEF regulates a different set of genes in the two petal layers

To explore how petal identity is defined by PhDEF in a layer-specific fashion, we generated scRNA-Seq data from 3875 and 3737 cells from star and wico mature petals, respectively. Due to the difficulty in collecting high numbers of flowers from these genetic chimeras, we were only able to produce a single biological replicate for each. We integrated the WT, star, and wico scRNA-Seq data using Harmony (Korsunsky et al., 2019), allowing for the identification of epidermal and mesophyll clusters in all samples (Figure 2A and B). The pseudo-bulk transcriptomes displayed a higher correlation between WT replicates (Pearson’s correlation coefficient, r=0.96) than with the other genotypes (r=0.9 between WT and wico, r=0.93 between WT and star), as expected. The levels of expression of the main petal identity genes PhDEF and PhGLO1 were consistent with our assumptions (Figure 2C–E): in WT petals, both genes were expressed evenly in all cell types, while in star petals we observed high expression in the mesophyll clusters and low expression in the epidermal clusters, and the converse was observed in wico petals. This is particularly apparent when looking at PhGLO1 expression, that is higher than PhDEF expression at this stage (Figure 2D).

Figure 2 with 2 supplements see all
Download asset Open asset
PhDEF regulates a different set of genes in the petal epidermis and mesophyll.

(A) Uniform Manifold Approximation and Projection (UMAP) plots of 11,632 WT, 3875 star, and 3737 wico petal cells sequenced for their transcriptome, after integration and clustering. Here, epidermal clusters (purple) and mesophyll clusters (green, without the vasculature) were merged. (B) Percentage of cells from the different epidermal and mesophyll identities identified from WT, star, and wico petal scRNA-Seq. (C) Per-cell expression levels of PhDEF and PhGLO1 in WT, star, and wico displayed on UMAP plots. (D) Dotplot of PhDEF and PhGLO1 expression level (color-coded) and percentage of cells expressing the gene (coded in the size of the dot) in the epidermis and mesophyll cells of WT, star, and wico petals. (E) Pearson’s correlation coefficient (r) plot between the pseudo-bulk transcriptomes from the mesophyll, epidermis, and vasculature cells extracted from the WT, star, and wico scRNA-Seq data. A Fisher r-to-z transformation test indicates that all correlation coefficients are significantly different (p<0.05). (F) Venn diagram of the number of epidermis (purple) and mesophyll (green) DEGs, revealing the number of epidermis- and mesophyll-specific DEGs, as well as common DEGs (i.e. differentially expressed in both layers), based on WT, star, and wico scRNA-Seq data. The percentage of activated (light gray) and repressed (dark gray) genes is also displayed. (G) Differential expression (log2FC = log2(FoldChange)) of the common DEGs in the epidermis or in the mesophyll, showing that almost all DEGs are either activated or repressed in both layers. Three genes of interest (PhDEF, PhGLO1, and PhGLO2) are displayed as color points. (H) Ten most enriched Gene Ontology (GO) terms for biological processes in epidermis-specific (left), mesophyll-specific (middle), and common (right) DEGs, after GO term redundancy reduction and sorting by p-value.

When comparing the star and wico petals with WT petals, we did not uncover new clusters, which one might have expected from altered cell identity, but we observed changes in the proportions of cells in each cluster (Figure 2B). First, in star, a marginal number of upper limb epidermal cells were recovered, which corresponds to the small secondary L1-revertant sectors that are frequently observed in star petals (Chopy et al., 2024; Figure 2—figure supplement 1) and that were thus later removed from the dataset. Apart from that, two major epidermal clusters (tube epidermis and limb epidermis) were strikingly depleted in cell numbers in star as compared to their WT counterparts. We interpret this as the result of two processes: (1) the lower proportion of limb tissue, that is mostly constituted of epidermal cells (Reale et al., 2002; Chopy et al., 2024), in star petals than in WT petals; (2) the alteration of epidermal cell identity by the phdef mutation, which might cause epidermal cells to cluster with mesophyll cells. In line with the latter, the pseudo-bulk transcriptome of the star epidermis was slightly more similar to the one from the WT mesophyll (r=0.92) than to the one from the WT epidermis (r=0.89; Figure 2—figure supplement 2).

In wico petals, the main epidermal clusters were recovered in proportions higher than in star petals, but still much lower than in WT petals (Figure 2B). Surprisingly, the small cluster of strongly pigmented limb epidermal cells, which has a very clear transcriptional signature in WT petals (cluster 6 in Figure 1C), was not found in wico flowers. We speculate that this is because the wico flowers have a reduced pigmentation as compared to WT petals (Figure 1A, Figure 2—figure supplement 1), which we previously found to be caused by the 6 bp insertion left by the transposon excision from the PhDEF locus, slightly altering the function of PhDEF in activating pigmentation (Chopy et al., 2024). Consistently, we found 13 genes, out of the 42 associated with anthocyanin biosynthesis and its regulation, to be differentially expressed in the wico epidermis as compared to the WT epidermis. The depletion in epidermal cells in the wico UMAP space could therefore result from this slightly altered epidermal identity. In line with that, the pseudo-bulk transcriptome of the wico epidermis was as similar to the one from the WT epidermis (r=0.9) than from the WT mesophyll (r=0.9; Figure 2—figure supplement 2). In contrast, in both star and wico, the mesophyll resembled more the WT mesophyll (r=0.94 and 0.92, respectively) than the WT epidermis (r=0.88 and 0.84, respectively; Figure 2—figure supplement 2). This global inspection of layer-specific star and wico transcriptomes suggested that alteration of epidermal identity tends to shift the transcriptome towards a mesophyll-like transcriptome, while alteration of mesophyll identity did not have the converse effect; therefore, the mesophyll identity appears as the default cell identity in the petal.

Next, we explored in detail the effect of phdef mutation on gene expression in a given petal layer. Since the star flowers are only mutant for PhDEF in the epidermis, we defined epidermal-specific DEGs as differentially expressed in the star epidermis as compared to the WT epidermis (with |log2FC|>1, adjusted p-value <0.01 and specificity of expression >10 %, see Materials and methods). A major limitation of this definition is that, as previously shown, star epidermal cells might cluster with mesophyll cells because their transcriptional identity has been altered. This would lead to an artificial depletion of star epidermal cells in the analysis, which would result in the underestimation of epidermal-specific DEGs. Conversely, mesophyll-specific DEGs were defined as the genes differentially expressed in the wico mesophyll as compared to the WT mesophyll. We found that 286 genes were DEGs in both layers, whereas 628 were exclusive DEGs in the epidermis and 364 were exclusive DEGs in the mesophyll (Figure 2F, Supplementary file 4). Hence, even accounting for a possible underestimation of epidermal-specific DEGs, we found more epidermis-specific DEGs than mesophyll-specific DEGs. These genes were either activated (75%) or repressed (25%) by PhDEF (Figure 2F), and there was a slightly higher proportion of genes whose expression was activated by PhDEF in the epidermis (78.5%) than in the mesophyll (69.3%, Figure 2F). The vast majority (99.3%) of the genes commonly regulated by PhDEF in both layers, hereafter named common DEGs, were regulated in the same direction in the two layers, that is either activated (75.2% of genes) or repressed (24.8% of genes, Figure 2G). Hence, the direction of regulation by PhDEF is not heavily influenced by the petal layer in which it acts. As expected, PhDEF, PhGLO1, and PhGLO2 were found among the 287 common DEGs, activated by PhDEF in the two petal layers (Figure 2G, Supplementary file 4). GO terms enriched in epidermis-specific DEGs were related to redox processes (glutathione metabolism, cell redox homeostasis), while those enriched in mesophyll-specific DEGs were strongly associated with photosynthesis, and various different GO terms were enriched in common DEGs (Figure 2H, Supplementary file 2).

In summary, our scRNA-Seq dataset from layer-specific phdef mutant flowers indicates that PhDEF regulatory action is largely different between the two cell layers, as the number of common DEGs is lower than the number of layer-specific DEGs, and this is consistent with the drastically different star and wico phenotypes. We also found about two times more epidermal-specific than mesophyll-specific DEGs.

PhDEF binds to more target genes in the petal epidermis than in the petal mesophyll

The scRNA-Seq transcriptomic data confounds the effect of direct gene regulation by PhDEF, and the indirect gene regulation that results from it, which might be extremely different in terms of numbers of genes, direction, and intensity of the transcriptional regulation among the pathways affected. To circumvent this limitation, we then attempted to identify PhDEF direct targets in the two petal layers. For this, we performed chromatin immunoprecipitation followed by sequencing (ChIP-Seq) in WT, star, and wico petals, using a custom antibody directed against the endogenous PhDEF protein (Chopy et al., 2024). After a set of preliminary tests to evaluate ChIP efficiency (Figure 3—figure supplement 1 and Materials and methods), we performed the experiment on two biological replicates, on flowers at stage 8 when the tube is about half its final size, in an attempt to capture most gene regulation events underlying both limb and tube development (Chopy et al., 2024). After sequencing and mapping the reads on the Petunia axillaris genome (Bombarely et al., 2016), we called the peaks for each replicate separately, then computed the Irreproducible Discovery Rate (IDR) (Li et al., 2011) to identify reproducible peaks and removed peaks detected in the input chromatin (see Materials and methods), which resulted in 3134, 4276, and 2216 peaks for WT, wico, and star, respectively (Supplementary file 5, Supplementary file 6). We observed that PhDEF bound to similar regions of the genome in WT, wico, and star, with a generally higher read coverage in WT and wico samples than in star samples (Figure 3A and B, Supplementary file 6). Importantly, this is consistent across the two ChIP replicates for each genotype (Supplementary file 5), and is not due to a more variable PhDEF expression in the mesophyll as compared to the epidermis (Figure 3—figure supplement 1), which could have caused a general lower enrichment of PhDEF binding in the mesophyll. This indicates that binding of PhDEF to the genome is more frequent in the epidermis than in the mesophyll, but overall takes place in the same regions of the genome.

Figure 3 with 2 supplements see all
Download asset Open asset
PhDEF binds to more target genes in the epidermis than in the mesophyll.

(A) Metaplot of PhDEF binding sites in WT, wico, and star petals, with the distance from the Transcription Start Site (TSS) and the Transcription Termination Site (TTS). (B) Heatmap of the read coverage of all peaks detected in WT samples, and of corresponding regions of the genome in wico and star samples. The peaks are sorted according to the highest to lowest read coverage in WT #2, so that the same regions of the genome are in the same line in all libraries. The read coverage is color-coded, and position 0 represents the start of the peak. (C) Read coverage for PhDEF binding to the genomic regions of PhDEF and PhGLO1 in WT, star, and wico petals. Peaks that pass our pipeline for reproducibility (see Materials and methods) are displayed as thick red lines. Predicted MADS-binding sites (MADS-bs) are indicated by blue lines, while no TCP-bs were found in these sequences. Although the peak in the promoter region of PhDEF does not pass the IDR threshold, its read coverage profile strongly suggests that PhDEF binds to this position. (D) Intersection of the number of gene-associated peaks from WT, star, and wico ChIP-Seq data, and definition of the possible PhDEF binding profiles. The intersection occasionally resulted in the artificial duplication of peaks, and once corrected, this marginally changed the total number of peaks identified for each genotype in the intersection. (E) Distribution of PhDEF binding profiles across its target genes. (F–G) PhDEF binding profile over the genomic locus of AN1 (F) and EOBII (G). Predicted MADS-binding sites (MADS-bs) are indicated by blue lines, and predicted TCP-binding sites (TCP-bs) are indicated by orange lines.

In total 2017, 2771, and 1136 peaks were associated with genes for WT, wico, and star, respectively (Supplementary file 6, Supplementary file 7). Most of these peaks were located in promoter regions (defined as 6 kb upstream of the start codon, 34.5% for WT, as compared to a random 14.9% occupancy, p=2.06e–153, hypergeometric test) or terminator regions (defined as 6 kb downstream of the stop codon, 20.1% for WT, p=6.44e–13) of genes, with a small and non-enriched fraction found in gene bodies (9.8% on average, p=0.51; Figure 3A, Supplementary file 6). These proportions were similar in WT and wico, but star displayed a lower proportion of peaks located in the promoter regions (24.5% for star vs. 34.5% and 36.2% for WT and wico, respectively, Supplementary file 6), suggesting that binding of PhDEF to promoter regions is more frequent in the epidermis than in the mesophyll. Out of the genes bound by PhDEF in WT, we found a substantial proportion (37.0%, p=1.20e–118, hypergeometric test) also differentially expressed in the previously published phdef-151 mutant transcriptome (Chopy et al., 2024), which is consistent with the proportions that have been reported from previous ChIP-Seq assays on MADS-box proteins (Wuest et al., 2012; Chen et al., 2018; Kaufmann et al., 2010; Kaufmann et al., 2009). In particular, we could detect clear binding of PhDEF on its own promoter, and on the promoter of the other B-class gene PhGLO1 (Figure 3C), that are two established regulatory targets of PhDEF, detected as differentially expressed both in bulk RNA-Seq and in scRNA-Seq.

We next examined the intersection of gene-associated peaks (including promoter, gene body and terminator regions) between WT, wico and star flowers (Figure 3D, Supplementary file 7). Strikingly, and consistently with the high number of peaks detected in wico, there were many more epidermal-enriched peaks (i.e. found in WT +wico or wico alone, 2041 peaks in 1896 genes) than mesophyll-enriched peaks (i.e. found in WT +star or star alone, 407 peaks in 398 genes). A large number of peaks (1,176 in 1,110 genes) were also shared between layers (i.e. found in WT, star +wico or WT +star + wico). When shared, most peaks were higher in wico than in WT than in star (see Figure 3C for PhGLO1 for instance), supporting again that the binding of PhDEF to its target sites is more frequent in the epidermis than in the mesophyll. Hence for most peaks, PhDEF binding in the WT sample can be interpreted as the average of frequent PhDEF binding in the epidermis and rare or absent PhDEF binding in the mesophyll, explaining why some peaks are detected in wico but not in WT. The 317 peaks uniquely found in star were also weakly bound in the epidermis, and the 1,147 peaks uniquely found in wico were also weakly bound in the mesophyll (Figure 3—figure supplement 2), showing that differential binding between layers is quantitative, with a relative enrichment in the epidermis or in the mesophyll that can differ across genomic regions.

While most of PhDEF target genes (91.8% on average) were associated with a single ChIP peak, others displayed between 2 and 4 binding sites for PhDEF (Supplementary file 6). To our surprise, these different binding sites over a single target gene were not necessarily occupied by PhDEF similarly in the two petal layers: for instance, the gene encoding the key regulator of anthocyanin biosynthesis ANTHOCYANIN1 (AN1) (Spelt et al., 2000; Koes et al., 2005) displayed three binding sites for PhDEF in its regulatory regions, two of them being epidermal-specific and one of them being shared between petal layers (Figure 3F). A similar binding profile, with one epidermal-specific and one shared peak, was observed for the gene encoding the MYB-R2R3 transcription factor EMISSION OF BENZENOIDSII (EOBII), a major regulator of scent emission and petal maturation (Spitzer-Rimon et al., 2010; Spitzer-Rimon et al., 2012; Chopy et al., 2023; Figure 3G). These complex binding profiles accounted for a rather small number of genes (153 genes with epidermis-specific +shared peaks, 23 genes with mesophyll-specific +shared peaks), compared to the vast number of genes displaying a single type of PhDEF binding profile (1694 genes with only epidermis-specific peaks, 977 genes with only shared peaks and 325 genes with only mesophyll-specific peaks). Hence, the PhDEF binding profile can be complex, with layer-specific or -aspecific binding sites, which might reflect differential transcriptional regulation between cell layers.

Epidermis DEGs display more epidermal and shared binding sites for PhDEF

We next wanted to explore further how different PhDEF binding profiles are associated with the expression profiles of its target genes by combining our transcriptomic and ChIP-Seq assays on WT, star, and wico flowers. We compared our ChIP-Seq results with either the scRNA-Seq dataset or with the bulk RNA-Seq dataset on WT, star, and wico petals that we previously obtained (Chopy et al., 2024). These two datasets have complementary benefits: the scRNA-Seq dataset provides a layer-specific view of gene differential regulation, while the bulk RNA-Seq dataset has been generated at the exact same stage as ChIP-Seq and better captures lowly expressed genes. Both bulk RNA-Seq and scRNA-Seq datasets overlap with the ChIP-Seq dataset more than by chance (Figure 4—figure supplement 1, p=4.73e–13 and p=1.37e–04 respectively, hypergeometric test).

Using bulk RNA-Seq data, we explored the distribution of PhDEF binding profiles for DEGs in wico and/or in star petals at stage 8 (Figure 4A), and compared it with the distribution we had previously computed for all PhDEF target genes (Figure 3E). Overall, distributions were very similar, with a majority of epidermis-specific peaks and shared peaks, as previously observed. A significantly different distribution of PhDEF binding profiles was observed for star-specific DEGs (Figure 4A, Chi2 goodness-of-fit, p=0.039), with more genes displaying epidermis +shared binding sites (test of equal proportions, one-sided, p=0.037). No significant deviation from the expected distribution is found for wico-specific DEGs (Chi2 goodness-of-fit, p=0.080) nor for common DEGs (p=0.064), but wico-specific DEGs were enriched for mesophyll-specific peaks (test of equal proportions, one-sided, p=0.024). To explore whether layer-specific regulation profiles might be blurred in the bulk RNA-Seq data, since it mixes the transcriptome of both cell layers, we analyzed PhDEF binding profile in the layer-specific DEGs identified by scRNA-Seq (Figure 4B). Again, distributions of PhDEF binding profiles were highly similar for epidermal-specific, mesophyll-specific, or shared DEGs (Figure 4B, Chi2 goodness-of-fit, p=0.079, 0.28, and 0.94, respectively), but we could find an enrichment of epidermis +shared binding sites in epidermis-specific DEGs (test of equal proportions, one-sided, p=0.033). Overall, and keeping in mind the limitations of these analyses, these results show that epidermis DEGs tend to display more epidermis-specific and shared binding sites for PhDEF, suggesting that layer-specific binding and transcriptional regulation are partially linked, at least in the epidermis.

Figure 4 with 3 supplements see all
Download asset Open asset
Intersection between PhDEF binding and regulatory profiles, and motif enrichment under PhDEF binding sites.

(A) Distribution of PhDEF binding profiles across genes differentially expressed (DEGs) in star only, wico only, or both in star and wico (common DEGs), as compared to WT, from bulk RNA-Seq of petals at stage 8. Stars indicate a significant deviation from expected PhDEF binding profiles, either for the whole distribution (red star, Chi2 goodness-of-fit test) or an enrichment of individual categories (black star, test of equal proportions, one-sided). Binding categories with less than 5 genes were not tested for significant enrichment. (B) Distribution of PhDEF binding profiles across epidermis-specific, mesophyll-specific, or common (in both layers) DEGs, obtained from WT, star, and wico scRNA-Seq by comparing layer-specific transcriptomes. (C) Selected motifs enriched under PhDEF ChIP-Seq peaks for WT, wico, and star petals, bs = binding site. The full list of motifs detected is in Figure 4-figure supplement 2.

Different motifs are enriched under layer-specific PhDEF binding sites

A possible mechanism by which PhDEF regulates different target genes in the two petal cell layers could be its interaction with different protein partners, which would modify its DNA-binding specificity. This should be readily detectable in the enriched motifs for transcription factor binding under ChIP-Seq peaks. When examining these motifs, we found that gene-associated ChIP-Seq peaks for WT, wico, and star were all enriched in CArG-boxes (Figure 4C, Figure 4—figure supplement 2), the typical MADS-box binding site. This likely reflects PhDEF direct binding to these regions but might also be the result of other MADS-box proteins binding, together with PhDEF as dimers or tetramers. Indeed, PhDEF is known to interact with the B-class proteins PhGLO1 and PhGLO2 (Vandenbussche et al., 2004), and with the E-class proteins FBP2 and FBP5 (orthologous to Arabidopsis SEPALLATA3 (SEP3) and SEP1/2, respectively) (Ferrario et al., 2003; Morel et al., 2019). We examined the cell layer-specific expression of these genes in our WT scRNA-Seq dataset and found that the B- and E-class genes were overall homogeneously expressed between the two petal layers or hardly detected (Figure 4—figure supplement 2A). It is therefore unlikely that different floral quartets would form in the epidermis and in the mesophyll. Apart from CArG boxes, we also found enrichment for TCP (teosinte branched 1/cycloidea/proliferating cell nuclear antigen factor 1) binding sites in WT and wico gene-associated peaks, but not in star (Figure 4C, Figure 4—figure supplement 2). TCP transcription factors are known to be partners of MADS-box proteins (Bemer et al., 2017; Wanamaker et al., 2017), and also to regulate the same genes as the ones regulated by MADS-box proteins (Zhao et al., 2020), explaining why their motifs are occasionally found in co-occurrence under ChIP peaks (Kaufmann et al., 2009; Goslin et al., 2023). The fact that this motif is not found enriched under star peaks could be due to two reasons: either the number of peaks in star was too low to actually detect such enrichment, or TCP proteins could be partners of PhDEF or coregulators of PhDEF target genes only in the petal epidermis. As a support to the second hypothesis, we detected a strong enrichment of predicted TCP-binding sites under epidermal-specific peaks and shared peaks, but only a weak one under mesophyll-specific peaks, as compared to random sequences (Supplementary file 8). For gene-associated peaks, 25.2% of epidermal-specific peaks, and 25.0% of shared peaks, had at least one predicted TCP binding site, in contrast to only 11.8% of mesophyll-specific peaks (Supplementary file 8). We also examined the expression of petunia TCP genes (Zhang et al., 2020) in our WT petal scRNA-Seq dataset and found a clear enrichment for PaTCP4a in the epidermis, while other TCP genes displayed low expression and/or similar expression in the two cell layers (Figure 4—figure supplement 3B). The expression of this gene increases during petal development, is lower in star than in wico flowers at stage 8 (Figure 4—figure supplement 3C), and is also higher in the star and wico epidermis than in their corresponding mesophyll tissue (Figure 4—figure supplement 3D). The expression of PaTCP4a was also reported to be higher in P. axillaris flowers with large petals than in flowers with small petals (Zhang et al., 2020). PaTCP4a therefore appears as a good candidate to be an epidermal-specific partner of PhDEF in the petal, which remains to be verified experimentally.

Discussion

In this work, we investigated the cell layer-specific functions of the petal identity regulator PhDEF in petunia, based on the initial observation that epidermal- or mesophyll-specific expression of PhDEF yielded drastically different petal phenotypes. We found that PhDEF regulates a different set of genes in the two petal layers, with a major binding and regulatory action in the epidermis. We also identified target genes with complex binding profiles, and overall the association between binding and regulatory profile of PhDEF target genes was limited, but significant for genes bound and regulated in the epidermis, displaying more epidermal +shared peaks than expected. This suggests that PhDEF can achieve layer-specific regulation of its target genes through layer-specific binding, in particular in the epidermis. The reasons for the otherwise low association between binding and regulatory profiles are likely manifold: PhDEF binding does not always entail target gene regulation; PhDEF can act either as an activator or as a repressor of gene transcription; long-distance looping between PhDEF binding sites (as has been shown for MADS-box proteins ; Mendes et al., 2013; Melzer and Theissen, 2009) creates complex binding profiles that exceed the environment of a single gene; and most importantly, many other transcription factors than PhDEF regulate the expression of its target genes.

Globally, gene transcriptional regulation was consistently more prominent in the epidermis than in the mesophyll. First, we observed that the epidermis expresses more specific genes than the mesophyll, which, to our knowledge, had not been reported before. Hence, the mesophyll might be regarded as a default identity for the petal cells, upon which the regulation of a large set of genes is appended for the specification of epidermal identity. Second, we found that phdef mutation leads to the differential expression of more genes in the epidermis than in the mesophyll. Third, we observed that PhDEF strikingly binds to more binding sites and target genes in the epidermis than in the mesophyll. This prominent role of PhDEF in the petal epidermis aligns with the multiple specific features found in this tissue, and in epidermal tissues in general: trichomes, stomata, glandular cells, or conical cells; and the very specific functions that the epidermis can fulfill for pigmentation, volatile emission, cuticle biosynthesis, or defense against biotic and abiotic stress (Robinson and Roeder, 2015; Cavallini-Speisser et al., 2021; Riglet et al., 2021; Zuch et al., 2022). In contrast, the mesophyll appears less differentiated, with cells visually quite similar and fulfilling less diverse functions than the epidermis (Cavallini-Speisser et al., 2021; Borsuk et al., 2022; Zhang and Ambrose, 2024), although this might be biased by the low accessibility of this tissue, and the reduced interest it has attracted. Several epidermal-specific genes have been described over the years, in particular the early specifiers of embryonic epidermal identity Meristem Layer 1 (ML1) or Protodermal Factor 2 (PDF2) from the previously cited family of HDG transcription factors (Lu et al., 1996; Abe et al., 2003; Bhatia et al., 2021; Schrick et al., 2023), or the genes acting in later epidermal-specific functions, such as stomata (Zuch et al., 2022) or trichome (Schellmann and Hülskamp, 2005) formation, cuticle synthesis (Fich et al., 2016), or pigmentation (Yeo and Moyroud, 2025). In contrast, a much smaller number of genes specific to the L2 layer have been described (Bhatia et al., 2021; Yadav et al., 2014), and these genes are all expressed post-embryonically (Bhatia et al., 2021), suggesting the possibility that the sub-epidermal identity does not require specific regulators for its determination. Therefore, establishing epidermal identity might truly require the regulation of several specific genes, starting from the canvas of a default mesophyll identity. Linking this prominent action of PhDEF in the epidermis with the star and wico phenotypes suggests that building the limb, which is mostly constituted of epidermal tissue, requires the regulation of a larger set of genes than building the tube.

We have two main hypotheses explaining why PhDEF adopts a different binding and regulatory behavior in the epidermis and in the mesophyll: First, PhDEF might have access to different protein partners in the two petal layers, and since binding to different protein partners entails changes in DNA-binding specificity, which has been demonstrated for MADS-box proteins (Smaczniak et al., 2017; Käppel et al., 2023), this would lead to the regulation of different target genes. Our identification of a TCP-binding motif enriched under epidermal-specific PhDEF binding sites suggests that MADS-TCP complexes could contribute to the specific regulation of PhDEF epidermal targets. However, it would account for 25% of PhDEF target sites at most. Second, PhDEF might have access to different target genes in the two petal layers due to different chromatin accessibilities. Recently, single nuclei ATAC-Seq (Assay for Transposase-Accessible Chromatin) assays in maize and sorghum leaves revealed that epidermal and mesophyll chromatin accessibilities were indeed very different, with more accessible regions in the epidermis than in the mesophyll in both species (Swift et al., 2024). MADS-box proteins have also been shown to interact with chromatin remodelers (Liu et al., 2011; Smaczniak et al., 2012), which might participate in the reinforcement of differential chromatin accessibility between layers. The fact that the difference between epidermal and mesophyll binding of PhDEF is mostly quantitative, with more frequent binding in the epidermis, seems in favor of an elevated chromatin accessibility in the epidermis. Apart from these two main non-mutually exclusive scenarios, a major aspect that remains unexplored in our study is the presence and activity of the PhDEF protein in the petal cell layers. The mere enrichment of PhDEF in the epidermis could explain its prominent regulatory role in this layer, which would then open investigations on the regulation of PhDEF translation, PhDEF protein movement between layers (Urbanus et al., 2009; Urbanus et al., 2010), or PhDEF protein degradation.

Overall, our study provides a cell-layer-specific view on transcriptional regulation, directly accessible here thanks to the star and wico periclinal chimeras. Other studies than ours have reported that regulators could have different functions in different cell layers (Carpenter and Coen, 1990; Carpenter and Coen, 1995; Bouhidel and Irish, 1996; Jenik and Irish, 2001; Sessions et al., 2000; Filippis et al., 2013), for other floral homeotic genes but not exclusively, which suggests that cell layer identity likely impacts transcription factor activity as a general rule.

Materials and methods

Plant material and culture conditions

Request a detailed protocol

Plants were grown in a culture room in long day conditions (16 hr day at 22 °C, 8 hr night at 18 °C, 75-Valoya NS12 LED bars, light intensity: 130 μE, 60% humidity). All plant material derives from the Petunia x hybrida R27 line containing the active dTph1 transposable element. The wico and star flowers spontaneously arose from phdef-151 homozygous mother plants, by the cell layer-specific excision of dTph1, resulting in branches carrying either wico or star flowers that were repeatedly obtained from several different phdef-151 individuals, as described in Chopy et al., 2024. These branches were subsequently maintained by cuttings.

Histology

Request a detailed protocol

Petal cross-sections and toluidine blue staining were performed as described in Chopy et al., 2024. Cell types were counted as indicated in Figure 1-figure supplement 1, from four different cross-section images per limb and per tube (n=1038 limb cells and 2736 tube cells counted in total).

Bulk RNA-Seq

Request a detailed protocol

From WT petals at stage 12 (Reale et al., 2002; Chopy et al., 2023), a fraction of the corollas was directly flash-frozen; the rest was protoplasted as described below and the protoplast pellet was flash-frozen. RNA extraction was performed with Sigma’s Spectrum Plant Total RNA Kit with On-Column DNase I Digestion, following manufacturer’s recommendations. RNA integrity and quantity were determined using a Bioanalyzer RNA 6000 Nano assay (Agilent), libraries were prepared with the CORALL RNA-Seq Library Prep (Lexogen) and sequenced with an Illumina NextSeq500 (single-end reads, 84 bp). Reads were mapped on the P. axillaris transcriptome (Bombarely et al., 2016) as described in Chopy et al., 2024. Differential gene expression was computed with DESeq2 (Love et al., 2014; 1.34.0) in RStudio (R 4.4.1, RStudio 2024.04.2).

RT-qPCR

Request a detailed protocol

For RNA extraction, limb or tube tissue was collected from WT flowers at anthesis in three biological replicates. For limb samples, the limb from one petal from three different flowers constitutes one replicate; for tube tissue, the whole tube from three different flowers constitutes one replicate. RNA extraction, Reverse Transcription (RT), quantitative PCR (qPCR), and result analysis were performed as described previously (Morel et al., 2019). Primer sequences are available in Supplementary file 9.

In situ hybridization

Request a detailed protocol

In situ hybridization was performed as previously described (Chopy et al., 2024), using f3h mutant flowers with white petals (Yu et al., 2021), to avoid detecting an aspecific signal in pigmented epidermal cells. The primers used to synthesize the antisense probe for IAA14 can be found in Supplementary file 9.

Single-cell RNA-Seq

Petal protoplast isolation

Request a detailed protocol

Protoplast isolation was performed in sterile conditions. WT, star, and wico corollas (ca. 15, 30, and 30 flowers per replicate, respectively) were briefly soaked in 70% ethanol and in 0.5% bleach for 30 s then rinsed three times with sterile water. Sterile corollas were transferred to a Petri dish containing 2 mL of Digestion Mix (0.4% macerozyme R-10, 0.8% Cellulase Onozuka R-10, w/v in TEX Buffer [3.1 g/L Gamborg B5 salts, 500 mg/L MES, 750 mg/L CaCl2*2H2O, 250 mg/L NH4NO3, 136.9 g/L Sucrose, pH 5.7]) and cut into ca. 0.5 cm² pieces using a new scalpel blade for each corolla to reduce tissue wounding. 10 mL of Digestion Mix was added, and digestion was performed for 5 hr at 26 °C in the dark, with gentle orbital agitation (20 rpm) for the last 15 min. After filtration through a 40-µm mesh, the volume was adjusted to 25 mL using 0.4 M sucrose (492 mOsm/kg-1 H2O) and the mixture was centrifuged for 10 min at 100 × g with acceleration 2/9 and deceleration 0/9 with a swing-out rotor. The underlying buffer was removed using a peristaltic pump (Gilson MINIPULS Evolution) connected to a sterile Pasteur pipette without perturbing the protoplast layer, at a rate of ca. 100 µL/s. After adjusting again the volume to 25 mL using 0.4 M sucrose, the whole process was repeated twice. Protoplast concentration and viability was assessed using a Kova slide and staining with 1% Evans Blue dye solution (w/v in 0.4 M Sucrose) as described in Fernandez-Da Silva and Menendez-Yuffa, 2006. Cell types were estimated as indicated in Figure S1, from four different pictures taken from a single protoplast isolation event (n=3338 cells counted in total).

From single-cell transcriptome sequencing to clustering

Request a detailed protocol

Protoplast suspensions were adjusted to 345, 480, 590, and 560 cells/µL with Mannitol-BSA (0.44 M Mannitol, 0.1% BSA (w/v), 498 mOsm.kg-1 H2O) for WT (first replicate), star, wico, and WT (second replicate), respectively. Protoplasts were loaded in the 10 x Genomics Chromium chip and libraries were prepared using the Chromium Next GEM Single Cell 3ʹ Reagent Kits v3.1 (dual index) kit following manufacturer’s instructions. Libraries were sequenced on an Illumina NovaSeq 6000 SP or S1. Read quality was checked with FastQC (Andrews, 2010; v0.12.1) and reads were aligned on the P. axillaris v1.6.2 transcriptome (annotation from Bombarely et al., 2016 transferred to a genome improved by HiC Chopy et al., 2024; Dudchenko et al., 2017, which was the best available version of a Petunia genome at that time). CellRanger count (v7.0.1, 10 X Genomics) was used to filter reads, count barcodes and UMIs, and generate the HDF5 matrix. The second WT replicate was split and sequenced in two technical replicates on the same Chromium chip, and read counts were aggregated with the Cell Ranger function aggr. Normalized gene expression, UMAP plots, and cluster markers were computed with Seurat 5.1.0 (Hao et al., 2024; Satija et al., 2015) in R (R Development Core Team, 2024) (v4.3.3) as follows, using default parameters unless stated otherwise: cells with less than 200 genes and genes detected in less than three cells were removed; data was normalized with NormalizeData, and highly variable features were identified with FindVariableFeatures and scaled with ScaleData; data dimension was reduced by Principal Component Analysis with RunPCA; Doublets were counted and removed with DoubletFinder (McGinnis et al., 2019); the number of significant Principal Components (PC) was determined with the function JackStraw (Chung and Storey, 2015); UMAP dimension reduction (McInnes et al., 2020) was performed with the function RunUMAP using the appropriate number of PCs. For the WT sample, MultiK (Liu et al., 2021) was used to explore the optimal number of clusters; a nearest-neighbor graph was constructed with the function FindNeighbors; clustering was computed with the function FindClusters, using the parameters determined by MultiK; cluster markers were identified with the function FindMarkers. For the identification of clusters, orthology to A. thaliana genes was defined as best or second best reciprocal Blast hit. A similar pipeline, up to clustering, was applied to the star and wico scRNA-Seq samples.

Data integration

Request a detailed protocol

To estimate the effect of protoplasting on WT scRNA-Seq clusters, genes differentially expressed after protoplasting with different log2FoldChange values (absolute log2FoldChange >1, 1.5, 2, or 2.5) were removed from the scRNA-Seq dataset, to generate four datasets. These were integrated with the original WT dataset with Harmony (Korsunsky et al., 2019; implemented in Seurat), and cell proportions per cluster were compared. To compare the WT, star, and wico samples, the two biological replicates of WT scRNA-Seq were first integrated with Harmony; then the WT, star, and wico samples (after computing the UMAP) were integrated together with Harmony, to ensure an unbiased representation of each genotype for cluster identification after integration.

Layer-specific gene differential expression

Request a detailed protocol

On the WT, star, and wico datasets after integration, the three epidermal clusters and the three mesophyll clusters were merged together by redefining their identities, while the vasculature and replicating/dying cells were kept untouched. For the star sample only, the small cluster of limb epidermal cells (secondary L1-revertants) was removed from the epidermal cluster. Only genes with more than 50 total counts in the dataset were kept for further analysis. Layer-specific DEGs were computed with the function FindMarkers, and subsetted for an absolute log2FoldChange >0.75 and layer-specificity>10% (computed as pct.1 – pct.2, pct.1 being the percentage of cells from group 1 expressing the gene, and pct.2 from group 2). DEGs in the star epidermis and wico mesophyll, as compared with the corresponding tissue in WT, were computed similarly and subsetted with an absolute log2FoldChange >1, adjusted p-value <0.01, and pct.1–pct.2 > 10%. These cut-offs were chosen in an attempt to capture a meaningful number of genes, while being stringent for their layer-specificity of expression. Specifically, Seurat proposes a default cut-off of 0.25 for the absolute log2FoldChange, that we considered to be too permissive for our own data. The cut-off for layer-specificity of expression was chosen at 10%, as the lowest cut-off for which the numbers of epidermis- and mesophyll-enriched genes were statistically different in the WT petal. Venn diagrams were generated with the help of Interactivenn (Heberle et al., 2015).

Chromatin immunoprecipitation and sequencing (ChIP-Seq)

ChIP assay and sequencing

Request a detailed protocol

Chromatin immunoprecipitation on WT, star, and wico petals at stage 8 was performed as previously described (Bouyer et al., 2018; Chopy et al., 2024), using polyclonal antibodies directed against PhDEF devoid of its highly conserved DNA-binding domain. Initial enrichment tests by qPCR on two positive control genes (PhDEF and Pos2) and one negative control gene (Neg1) allowed selecting the appropriate amount of sonicated chromatin, set at 25 µL, 100 µL, and 100 µL for WT, star, and wico, respectively. Four replicate ChIP experiments were performed for each genotype, and enrichment was measured by ChIP-qPCR (Figure S10) to select the two best replicates to send for sequencing. The quality of IP and INPUT (WT chromatin) samples was assessed with Tapestation 4150 HS D5000 (Agilent). Libraries were prepared with MicroPlex Library Preparation Kit v3 (Diagenode) and around 1 ng of each library was PCR-amplified (12 cycles), following Diagenode’s recommendations. Libraries were analyzed TapeStation 4150 HS D5000 and quantified with Qubit 4.0 with the Qubit dsDNA HS Assay Kit (Thermofisher), then mixed at equimolar ratio for sequencing with an Illumina NextSeq 500 (paired-end, 2×76 bp, dual indexing). Read quality was assessed with FastQC (Andrews, 2010) and MultiQC (Ewels et al., 2016), and reads were trimmed with fastx-trimmer (Hannon, 2010) to remove the first 6 nucleotides on the 5’ end and the last nucleotide on the 3’ end of each read, resulting in 2×68 bp reads. Amplification bias was removed as previously described (Leemans et al., 2023). Reads were mapped on the P. axillaris genome (v1.6.2 superscaffolded by HiC Bombarely et al., 2016; Dudchenko et al., 2017; Chopy et al., 2024), which was the most recent version of the genome available at the time data analysis was initiated. Mapping was performed with Bowtie2 (Langmead and Salzberg, 2012; v2.4.2) with default parameters, and only reads with mapping quality above 20 were kept. The effective genome size was estimated at 9.99e8 with khmer (Crusoe et al., 2015), with a kmer size of 68.

Peak detection and annotation

Request a detailed protocol

We explored different strategies for peak calling, but the standard MACS2 pipeline was not capturing PhDEF binding on its own promoter, which we used as our positive control, likely because the input sample was too noisy. Therefore, we explored different thresholds for peak calling with MACS2 (Zhang et al., 2008; False Discovery Rate = FDR) and for the IDR (Li et al., 2011) between IP replicates. Applying a FDR of 0.1 and an IDR of 0.1 to IP samples, and a FDR of 0.05 to INPUT samples, allowed to capture visually striking binding events in the IP and to remove spurious peaks in the INPUT in a satisfactory way. Peaks detected in the INPUT sample were removed from the reproducible IP peaks if they overlapped by at least 25% reciprocally, or if INPUT peaks were included in IP peaks, which was computed with a custom script using a combination of bedtools (Quinlan and Hall, 2010; v2.26.0) commands (intersect, window, subtract). Such customization of peak detection has been already applied to other transcription factors (Durand et al., 2023). Annotation of genes close to peaks was performed with Python. In cases when a peak was assigned to several possible genes, the gene with the shortest distance to the peak was retained. Peak intersections were performed independently using bedtools intersect on (i) the complete peak sets from WT, star, and wico and (ii) the subset of peaks located within 6 kb upstream of the transcription start site (TSS) or 6 kb downstream of the transcription termination site (TTS) of annotated genes, hereafter referred to as gene-associated peaks. Shared peaks were defined by a reciprocal overlap of at least 25%, in order to identify marked regions in one, two, or three genotypes. The distance of 6 kb, which was retained for the rest of the analysis, was chosen based on the metaplots, and because it encompasses 80% of the ChIP peaks. Motifs under the peaks were explored with MEME-ChIP (v5.5.7; Machanick and Bailey, 2011) with a minimal size of 6 bp and a maximal size of 15 bp against the JASPAR 2024 (Rauluseviciute et al., 2024). CORE non-redundant v2 database, and selected motifs (WT motifs from Figure 3B) were used to predict binding sites on individual genomic sequences with FIMO (Grant et al., 2011). Random genomic sequences, with the same size and size distribution as test files, were generated from the P. axillaris genome, for gene regions only or for the whole genome, to test for enrichment in TCP-binding sites. Peaks and gene models were visualized with the Integrative Genome Viewer (Robinson et al., 2023). For heatmaps of read coverage, coverage files (bigwig files) were calculated from BAM files using bamCoverage from the package deepTools2 (Ramírez et al., 2016; v3.5.4). Heatmaps were built using computeMatrix and plotHeatmap functions using a window of 3 kb beginning at the start of the peak.

GO term database improvement and enrichment analysis

Request a detailed protocol

Protein sequences were predicted from the P. axillaris genome annotation transferred from the published version v1.6.2 (Bombarely et al., 2016) to the Hi-C superscaffolded one using GFFread (Pertea and Pertea, 2020) and TransDecoder.LongOrfs (TransDecoder, 2025); with parameter ‘-m50’. A total of 56,351 predicted protein sequences were obtained. We then applied custom Perl scripts (refer to Data availability section) to select the most plausible CDS for each predicted mRNA, resulting in 31,626 predicted protein sequences. These sequences were then compared with UniProtKB protein database version 2024_06 (UniProt Consortium, 2023). The complete Swiss-Prot database of curated proteins (containing 44,534 plant sequences and 571,609 non-plant sequences) and the plant subset of the noncurated database TrEMBL (containing 20,024,921 sequences) were used. The comparison was performed using BlastP (Gish, 2026) program (AB-Blast v3.0) with parameters ‘W=3 Q=7 R=2 matrix = BLOSUM80 B=200 V=200 E=1e-6 hitdist = 60 hspsepqmax = 30 hspsepsmax = 30 sump postsw’. The BLAST output was filtered using custom Perl scripts to retain only those matches with a log10(e-value) no lower than 75% of the best log10(e-value). The GO annotations for proteins contained in the Uniprot KnowledgeBase were downloaded as a GPA file (v222, 2024-08-01) from the GOA database (Huntley et al., 2015) hosted at EBI. The GO database (Ashburner et al., 2000) was retrieved via the QuickGO web interface (Binns et al., 2009) in October 2024; GO term relationships and UniprotKB GO annotations were subsequently loaded into a custom PostgreSQL database (v13.13) for easy querying. Custom Perl scripts were written to associate each predicted petunia protein sequence to the GO terms associated with all its matched proteins or with at least five matched proteins. The comprehensive list of GO annotations for the entire Petunia proteome is accessible in the project dataset as a text file (see Data availability section), with the list of underlying Blast evidences. Enrichment analysis for test genes was performed using all genes expressed in the WT petal bulk RNA-Seq as background genes. A hypergeometric test (R version 4.2.2 R Development Core Team, 2024) was applied to assess the significance of enrichment of each subset (Pavlidis et al., 2004; Falcon and Gentleman, 2007). GO terms enriched (p-value <0.01 with a hypergeometric test) and with at least two genes in the pathway were analyzed with REVIGO (Supek et al., 2011; Tiny mode) to reduce redundancy. Enriched GO terms for biological processes were sorted by their log10(p-value), and the 10 best terms were displayed in barplots generated with ggplot2 (Wickham, 2016) in R.

Data availability

All data and code are publicly available as of the date of publication. Bulk RNA-Seq raw fastq files and normalized counts after DESeq2 have been deposited at GEO (GSE336494). Single-cell RNA-seq raw fastq files, h5 matrices, and rds images after integration have been deposited at GEO (GSE290697). ChIP-Seq fastq and BED files have been deposited at GEO (GSE308255). Scripts for bulk RNA-Seq analysis, GO database update, and GO term enrichment, and scRNA-Seq analysis are available at https://gitbio.ens-lyon.fr/rdp/phdef_flower_layer (Cavallini-Speisser et al., 2026). ChIP-Seq BAM and BAI files, ChIP-Seq and scRNA-Seq FastQC htmlHTML reports, and functional annotation of the Petunia axillaris proteome are available at https://entrepot.recherche.data.gouv.fr/dataverse/PhDEF_Flower_layer.

The following data sets were generated
    1. Cavallini-Speisser Q
    2. Monniaux M
    (2026) NCBI Gene Expression Omnibus
    ID GSE336494. Cell layer-specific functions of the petal identity regulator PhDEF in petunia II.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE336494
    1. Cavallini-Speisser Q
    2. Monniaux M
    (2025) NCBI Gene Expression Omnibus
    ID GSE290697. Cell layer-specific functions of the petal identity regulator PhDEF in petunia.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE290697
    1. Desert E
    2. Monniaux M
    (2025) NCBI Gene Expression Omnibus
    ID GSE308255. Cell layer-specific functions of the petal identity regulator PhDEF in petunia [ChIP-seq].
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE308255
    1. Just J
    2. Cavallini-Speisser Q
    (2025) Recherche Data Gouv
    Petunia axillaris proteome functional annotation.
    https://doi.org/10.57745/GUUBSE
The following previously published data sets were used
    1. Chopy M
    (2023) NCBI BioProject
    ID PRJNA951505. Transcriptome sequencing of petunia petals and cell layer-specific petal identity mutants.
    https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA951505

References

  2. Software
    1. Andrews S
    (2010)
    FastQC: a quality control tool for high throughput sequence data
    FastQC.
  35. Website
    1. Gish WA
    (2026) AB-BLAST
    Accessed March 17, 2020.
    https://blast.advbiocomp.com
  38. Website
    1. Hannon GJ
    (2010) FASTX-Toolkit
    Accessed January 15, 2023.
    http://hannonlab.cshl.edu/fastx_toolkit/
    1. Quemener E
    2. Corvellec M
    (2013)
    SIDUS—the solution for extreme deduplication of an operating system
    Linux J 3:3.
  80. Software
    1. R Development Core Team
    (2024) R: a language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
    http://www.r-project.org

Article and author information

Author details

  1. Emma Désert

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Data curation, Formal analysis, Investigation, Visualization, Methodology, Writing – review and editing
    Contributed equally with
    Quentin Cavallini-Speisser
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0001-3087-5559

  2. Quentin Cavallini-Speisser

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Data curation, Software, Formal analysis, Validation, Investigation, Methodology, Writing – review and editing
    Contributed equally with
    Emma Désert
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0008-1989-0310

  3. Evelyne Duvernois-Berthet

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Software, Formal analysis, Visualization, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1442-6603

  4. Pierre Chambrier

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5390-1345

  5. Patrice Morel

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared

  6. Brice Letcher

    Laboratoire de Biologie et Modélisation de la Cellule, École Normale Supérieure de Lyon, CNRS, UMR 5239, Inserm, U1293, Université Claude Bernard Lyon 1, Lyon, France
    Contribution
    Software, Formal analysis, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8921-6005

  7. Carine Rey

    CIRI, Centre International de Recherche en Infectiologie, Univ Lyon, Inserm, U1111, Université Claude Bernard Lyon 1, CNRS, UMR5308, ENS de Lyon, Lyon, France
    Contribution
    Software, Formal analysis, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4871-8110

  8. Jeremy Just

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Software, Formal analysis, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0842-9808

  9. Suzanne Rodrigues Bento

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Methodology
    Competing interests
    No competing interests declared

  10. Daniel Bouyer

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Project administration, Writing – review and editing
    For correspondence
    daniel.bouyer@ens-lyon.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7438-5551

  11. Marie Monniaux

    Laboratoire de Reproduction et Développement des Plantes, Université de Lyon, ENS de Lyon, UCB Lyon 1, CNRS, INRAE, Lyon, France
    Present address
    Univ. Lille, CNRS, UMR 8198, Evo-Eco-Paleo, Lille, France
    Contribution
    Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    marie.monniaux@univ-lille.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6847-3902

Funding

Agence Nationale de la Recherche

https://doi.org/10.67599/anr-19-ce13-0019
  • Marie Monniaux

Agence Nationale de la Recherche

https://doi.org/10.67599/anr-21-ce12-0036
  • Daniel Bouyer

Agence Nationale de la Recherche (ANR-21-ESRE-00016)

  • No recipients declared.

Agence Nationale de la Recherche (ANR-21-CE20-0007-02)

  • Daniel Bouyer

Agence Nationale de la Recherche (ANR-21-CE12-0036-01)

  • Daniel Bouyer

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Francesco Quattrocchio and Shuangjiang Li for sharing their protocol for protoplast isolation; Cyril Dégletagne and the Cancer Genomics Platform from the Centre de Recherche en Cancérologie de Lyon for support in the scRNA-Seq experiment; Nicolas Dalle and Annick Dubois for assistance in chromatin immunoprecipitation; Benjamin Gillet and Sandrine Hugues from the sequencing platform of the Institut de Génomique Fonctionnelle de Lyon for library preparation and sequencing of the transcriptomes of this study; the Master of Biology of the Ecole Normale Superieure de Lyon for their participation in funding and analyzing the ChIP-seq experiment; and Sergio Sarnataro from Spatial-Cell-ID for advice on scRNA-Seq analysis. We also thank three anonymous reviewers for their constructive criticism that greatly helped improve this manuscript. This work was supported by a grant to QC-S and MM from the Agence Nationale de la Recherche (grant ANR-19-CE13-0019, FLOWER LAYER), by a grant to ED from the French Ministry of Higher Education and Research, by grants to DB from the Agence Nationale de la Recherche (grants ANR-21-CE12-0036-01 and ANR-21-CE20-0007-02), and by the EquipEx +Spatial-Cell-ID under the "Investissements d'avenir" program (ANR-21-ESRE-00016). We gratefully acknowledge support from the CBPsmn (PSMN, Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. The platform operates the SIDUS solution (Quemener and Corvellec, 2013) developed by Emmanuel Quemener.

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.113154. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2026, Désert, Cavallini-Speisser et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 267
    views
  • 22
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Emma Désert
  2. Quentin Cavallini-Speisser
  3. Evelyne Duvernois-Berthet
  4. Pierre Chambrier
  5. Patrice Morel
  6. Brice Letcher
  7. Carine Rey
  8. Jeremy Just
  9. Suzanne Rodrigues Bento
  10. Daniel Bouyer
  11. Marie Monniaux
(2026)
The petal identity gene PhDEF has a major binding and regulatory action in the epidermis
eLife 15:RP113154.
https://doi.org/10.7554/eLife.113154.2

Share this article

https://doi.org/10.7554/eLife.113154

Categories and tags