R2R3-MYB genes coordinate conical cell development and cuticular wax biosynthesis in Phalaenopsis aphrodite.

Petals of the monocot Phalaenopsis aphrodite (Orchidaceae) possess conical epidermal cells on their adaxial surfaces, and a large amount of cuticular wax is deposited on them to serve as a primary barrier against biotic and abiotic stresses. It has been widely reported that subgroup 9A members of the R2R3-MYB gene family, MIXTA and MIXTA-like in eudicots, act to regulate the differentiation of conical epidermal cells. However, the molecular pathways underlying conical epidermal cell development and cuticular wax biosynthesis in monocot petals remain unclear. Here, we characterized two subgroup 9A R2R3-MYB genes, PaMYB9A1 and PaMYB9A2 (PaMYB9A1/2), from P. aphrodite through the transient overexpression of their coding sequences and corresponding chimeric repressors in developing petals. We showed that PaMYB9A1/2 function to coordinate conical epidermal cell development and cuticular wax biosynthesis. In addition, we identified putative targets of PaMYB9A1/2 through comparative transcriptome analyses, revealing that PaMYB9A1/2 acts to regulate the expression of cell wall-associated and wax biosynthetic genes. Furthermore, a chemical composition analysis of cuticular wax showed that even-chain n-alkanes and odd-chain primary alcohols are the main chemical constituents of cuticular wax deposited on petals, which is inconsistent with the well-known biosynthetic pathways of cuticular wax, implying a distinct biosynthetic pathway occurring in P. aphrodite flowers. These results reveal that the function of subgroup 9A R2R3-MYB family genes in regulating the differentiation of epidermal cells is largely conserved in monocots and dicots. Furthermore, both PaMYB9A1/2 have evolved additional functions controlling the biosynthesis of cuticular wax.

Introduction

Plants have developed a range of specialized cells on the epidermis of different organs, such as pavement cells on the leaves or conical cells on the petals, for environmental adaptation (Noda et al., 1994; Javelle et al., 2011). Petals are often the most visible part of a flower and are likely to play essential roles in successful pollination. Conical epidermal cells serve as a feature of petals with various functions in enhancing color brightness, reflecting sunlight, reducing wettability, and increasing pollinator attraction (Whitney et al., 2011). The MIXTA and MIXTA-like genes, members of the subgroup 9A R2R3-MYB gene family in eudicots (Brockington et al., 2013), function to regulate the differentiation of various epidermal cell classes, including cotton fibers (Machado et al., 2009; Walford et al., 2011), trichomes (Perez-Rodriguez et al., 2005; Plett et al., 2010), and conical epidermal cells (Noda et al., 1994; Baumann et al., 2007; Di Stilio et al., 2009; Lashbrooke et al., 2015). In contrast, in monocots, little is known about the biological functions of subgroup 9A members of the R2R3-MYB gene family in controlling epidermal cell fate. DcMYBML1, a MIXTA-like gene from Dendrobium crumenatum, can rescue the large extra-branched trichome phenotype when ectopically expressed in Arabidopsis (Arabidopsis thaliana) noeck mutant plants. However, D. crumenatum flowers, in which the DcMYBML1 transcripts are expressed, do not possess trichomes (Gilding and Marks, 2010).

In addition to the evolution of diverse structures and functions of epidermal cells, innovations have occurred in the cuticle coating on the surface of epidermal cells of aerial plant parts to cope with terrestrial environments. The plant cuticle is a continuous hydrophobic layer in which wax is embedded in (intracuticular layer) or overlaid on (epicuticular layer) a polymerized cutin matrix to form hydrophobic layers (Kunst and Samuels, 2009). Furthermore, hydroxyl and dicarboxylic fatty acids, terpenoids, phenylpropanoids, flavonoids, and glycerol can also be cuticular components (Kolattukudy, 2001; Beisson et al., 2012; Yeats and Rose, 2013). The outer layer of the cuticle is termed the cuticle proper and is mainly composed of cuticular waxes. It serves as a water-permeable barrier to prevent the evaporation of water from the epidermis and to provide the primary protective barrier against biotic (pathogens and pests) and abiotic (UV radiation and desiccation) agents (Yeats and Rose, 2013; Martin and Rose, 2014). Cuticular wax in A. thaliana consists of very long-chain fatty acids (VLCFAs: ≥20 carbon atoms) and their derivatives, including aldehydes, alkanes, secondary alcohols, ketones, and wax esters (Kunst and Samuels, 2009). The biosynthetic pathways of cuticular wax have been reported. The C16 and C18 fatty acids synthesized in the plastid are further elongated into VLCFAs by fatty acid elongase complex on the endoplasmic reticulum (Kunst and Samuels, 2009). VLCFAs are further modified via the alcohol-forming pathway (also referred to as the acyl reduction pathway) or alkane-forming pathway (also referred to as the decarbonylation pathway; Millar et al., 1999; Lee et al., 2016). The alcohol-forming pathway is responsible for primary alcohol and wax ester production, in which fatty acyl-CoA reductase (FAR3/CER4) functions to reduce VLCFAs into primary alcohols that are subsequently catalyzed by a bifunctional wax synthase/acyl-CoA:diacylglycerol acyltransferase to form wax esters (Rowland et al., 2006; Li et al., 2008). On the other hand, the alkane-forming pathway is responsible for alkane production via the ECERIFERUM 1 (CER1), CER3, and cytochrome b5 protein complex using aldehydes as precursors; the alkanes are further oxidized into secondary alcohols and ketones by midchain alkane hydroxylase (Greer et al., 2007; Bernard et al., 2012). Finally, these VLCFAs and their derivatives are secreted and deposited on the surface of epidermal cells to form cuticular wax layers (Kunst and Samuels, 2009).

A few SBG9A R2R3-MYB genes have evolved additional functions beyond regulating epidermal cell morphogenesis. For instance, in A. thaliana, the MIXTA-like gene MYB106 has been shown to regulate the morphogenesis of conical epidermal cells, trichome branching, and cuticle production (Oshima et al., 2013). MYB106 acts as a positive regulator of WIN1/SHN1; the overexpression of WIN1/SHN1 in A. thaliana results in increased cutin and cuticular wax loads on rosette leaves and flowers and cutin nanoridges decorating cauline leaves, suggesting that WIN1/SHN1 is likely to regulate the biosynthesis of cutin and cuticular wax (Aharoni et al., 2004; Broun et al., 2004). In contrast, A. thaliana 35S:SmiR-SHN1/2/3 plants, in which the three SHINE clade transcription factors (TFs; SHN1/2/3) are co-silenced by using an artificial microRNA approach, show no alteration of the cutin load on leaves but a significantly reduced cutin load on flowers compared with the wild-type (WT). However, the cuticular wax load on leaves and flowers is not significantly altered (Shi et al., 2011). In addition, WIN1/SHN1 has been demonstrated to directly positively regulate the expression of the cutin biosynthetic pathway gene Long-chain acyl-CoA synthetase 2; the increased cuticular wax load on the flowers of WIN1/SHN1-overexpressing lines is possibly an effect of increased cutin load (Kannangara et al., 2007; Oshima et al., 2013). Thus, to our knowledge, the transcriptional regulators that act directly to regulate the biosynthesis of cuticular wax during flower development have yet to be confirmed.

Phalaenopsis aphrodite flowers, and in particular their sepals and petals, strongly reflect light (i.e. glitter) and have a strong capacity for waterproofing and limiting water loss. The features are determined by the presence of cuticular waxes in the flowers (Kunst and Samuels, 2009; Supplemental Figure S1). However, the cuticular wax composition and the molecular pathway underlying cuticular wax biosynthesis remain poorly understood in Phalaenopsis spp. In addition, in monocots, the TFs that act to regulate the development of conical epidermal cells have also not been experimentally confirmed. In this study, we characterized two SBG9A R2R3-MYB genes, PaMYB9A1 and PaMYB9A2, and our data show that both PaMYB9A1 and PaMYB9A2 play pivotal roles in regulating the morphogenesis of conical epidermal cells and the biosynthesis of waxy substances in P. aphrodite petals.

Results

Morphogenesis of conical epidermal cells in P. aphrodite petals

To obtain more detailed information about the morphogenesis of conical epidermal cells in P. aphrodite, we divided the developing P. aphrodite flowers into five developmental stages, including 0.7 cm (Stage 1), 1.2 cm (Stage 2), 1.7 cm (Stage 3), and 2.2 cm (Stage 4) buds in lengths and fully bloomed flowers (Stage 5; Figure 1A). We then observed the adaxial surfaces of petals at each developmental stage by scanning electron microscopy (SEM; Figure 1B). SEM images show that at Stages 1 and 2, flat epidermal cells with a pentagon or hexagon shape are prominent on the adaxial epidermis of the petal. In addition, the dividing cells show a shallow line between independent descendant cells, suggesting that these two stages are undergoing proliferation (Figure 1B). By Stage 3, most of the epidermal cells exhibit a slight bulge at the center of the cells. After that, epidermal cells show a dome shape at Stage 4. Finally, differentiated conical epidermal cells can be observed at Stage 5 (Figure 1B). These results suggest that the differentiation of P. aphrodite petal conical cells occurs between Stages 3 and 5.

Figure 1

Morphogenesis of the petal conical epidermal cells. A, Definition of the five developmental stages of P. aphrodite flowers. B, Scanning electron micrograph of epidermal cells present on the adaxial epidermis of the petals at different developmental stages. Bars = 50 µm (C) relative expression of PaMYB9A1 and PaMYB9A2 by RT-qPCR at five distinct developmental stages during petal development. The y-axis shows relative RNA levels normalized to the expression levels of Ubi10. Data show the means ± standard deviation (sd) from three biological replicates, and each replicate comprises three technical repeats (n = 9).

Characterization of subgroup 9A R2R3-MYB TFs in P. aphrodite

MIXTA and MIXTA-like TFs belong to the subgroup 9A R2R3-MYB gene family, which has previously been reported to be involved in conical epidermal cell differentiation in plants including Antirrhinum majus, Petunia hybrid, A. thaliana, and Thalictrum thalictroides (Brockington et al., 2013). To identify homologs of the subgroup 9A R2R3-MYB gene family in P. aphrodite, we first identified all R2R3-MYB genes from the P. aphrodite genome by a BLAST-based search strategy using a conserved R2R3-MYB DNA binding domain sequence as a query (Chao et al., 2018). Amino acid sequences of the top genes identified in the BLAST search were aligned together with previously identified subgroup 9A R2R3-MYB genes and subjected to phylogenetic analysis. The phylogenetic tree shows that two P. aphrodite genes (PAXXG029600 and PAXXG123420, named PaMYB9A1 and PaMYB9A2, respectively) are classified into the subgroup 9A clade together with other MIXTA and MIIXTA-like orthologues in dicots (Supplemental Figure S2). According to the gene expression data retrieved from the Orchidstra version 2.0 database, these two genes are most highly expressed in the developing flower buds (Supplemental Figure S3; Chao et al., 2017). To investigate the association between the transcript abundance of PaMYB9A1 and PaMYB9A2 and the morphogenesis of conical epidermal cells on P. aphrodite petals, we used reverse-transcription quantitative PCR (RT-qPCR) to examine the patterns of PaMYB9A1 and PaMYB9A2 expression in the petals at different developmental stages. The RT-qPCR results showed that PaMYB9A1 transcripts accumulate at Stage 1, peak at Stage 3, and are barely detectable at Stage 5, and that PaMYB9A2 transcripts are predominantly expressed at Stages 3 and 4 (Figure 1C). The temporal patterns of both PaMYB9A1 and PaMYB9A2 expression are positively associated with the developmental period (between Stages 3 and 5) of conical epidermal cells during petal development (Figure 1, B and C), suggesting that these two genes might be involved in conical epidermal cell development.

Weighted gene coexpression network analysis of developing P. aphrodite petals

PaMYB9A1 and PaMYB9A2 exhibit distinct expression patterns during petal development, suggesting that these genes are involved in different regulatory networks (Figure 1C). Systems biology approaches were adopted to illuminate the transcriptional regulatory networks of these two TFs. The transcriptomic profiles of petals of each developmental stage were examined by RNA sequencing with biological replicates. Differentially expressed genes (DEGs) between Stages 1 and 2 (1,137 DEGs), 2 and 3 (4,108 DEGs), 3 and 4 (4,639 DEGs), and 4 and 5 (5,332 DEGs) were revealed (Supplemental Figure S4) by comparative transcriptome analysis. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs of different stages were also performed (Supplemental Figures S5 and S6). Furthermore, Weighted gene coexpression network analysis (WGCNA) was used to analyze the transcriptional regulatory networks of PaMYB9A1 and PaMYB9A2 (Langfelder and Horvath, 2008). A module is defined as a cluster of genes that have high correlation coefficients and interconnections. Thirteen distinct modules were identified and are shown in a gene dendrogram (Supplemental Figure S7). PaMYB9A1 and PaMYB9A2 are classified in the yellow and brown modules, respectively (Figure 2A). Furthermore, module-developmental stage eigengene analysis showed that the yellow module is highly associated with Stages 2–4, and the brown module is highly associated with Stages 3 and 4 (Figure 2B). In addition, an analysis of adjacent module-eigengenes indicated that the two modules show a strong relationship (Figure 2B). Hence, we speculate that PaMYB9A1 and PaMYB9A2 are involved in two distinct but associated regulatory pathways that control the differentiation of the conical epidermal cells.

Figure 2

WGCNA. A, Module-stage association. The number of genes in each module is shown. Each row corresponds to a module. Each column corresponds to a developmental stage of the petal. Three biological repeats of each developmental stage are shown. Pearson correlation coefficient of each module with different developmental stages are given and colored according to the score. B, Module–module association. A high degree of correlation is indicated by red, and a lower degree of correlation is indicated by blue.

Ectopic overexpression of either PaMYB9A1 or PaMYB9A2 results in shiny leaves in transgenic tobacco W38

Several subgroup 9A R2R3-MYB genes that act to regulate the morphogenesis of varied epidermal cells have been functionally characterized by ectopic overexpression in tobacco plants (Brockington et al., 2013). To verify whether PaMYB9A1 and PaMYB9A2 function in conical epidermal cell development, we generated two transgenic tobacco lines, 35S:PaMYB9A1 and 35S:PaMYB9A2, harboring the PaMYB9A1 and PaMYB9A2 coding sequences driven by the Cauliflower mosaic virus 35S promoter, respectively (Supplemental Figure S8). The expression of PaMYB9A1 and PaMYB9A2 in transgenic tobacco plants was confirmed by reverse transcription-PCR (RT-PCR) (Supplemental Figure S9). SEM images show that ectopic expression of both PaMYB9A1 and PaMYB9A2 does not alter the morphology of leaf epidermal cells compared to that of WT tobacco (Supplemental Figure S10). Unexpectedly, 13 of 18 35S:PaMYB9A1 lines and 15 of 19 35S:PaMYB9A2 lines showed a shiny phenotype only on the abaxial surface of the mature leaves (Figure 3). Previous studies indicated that A. thaliana 35S:WIN1/SHN1 leaves have a shiny phenotype caused by the large amount of wax deposited on the epidermis (Aharoni et al., 2004); however, 35S:PaMYB9A1/2 shiny leaves showed no noticeable wax deposition on their epidermis when observed by SEM (Supplemental Figure S10). The effects of PaMYB9A1 and PaMYB9A2 activities on the cuticular wax load and composition were further evaluated.

Figure 3

Phenotypes of transgenic tobacco plants overexpressing PaMYB9A1 and PaMYB9A2. A, Adaxial surface and (B) abaxial surface of WT, 35S:PaMYB9A1, and 35S:PaMYB9A2 mature leaves. Transgenic tobacco plants overexpressing PaMYB9A1 or PaMYB9A2 show shiny phenotypes on the abaxial surface of their mature leaves. Scale bar = 1 cm.

We extracted cuticular wax from the abaxial epidermis of transgenic tobacco and WT leaves and analyzed these waxes using gas chromatography-mass spectrometry (GC–MS). GC–MS analyses revealed that, compared with the WT, the cuticular wax loads on the two types of shiny leaves were not significantly altered (Figure 4A), but the chemical composition of cuticular wax was distinctly altered (Figure 4B). The cuticular wax compositions of both types of shiny leaves exhibited a distinct decrease in the proportion of n-alkanes of ≤25C and increased amounts of C54 n-alkanes compared with the WT (Figure 4B). Altogether, these results revealed that PaMYB9A1 and PaMYB9A2 act to regulate the biosynthesis of cuticular wax in transgenic tobacco leaves.

Figure 4

Analyses of cuticular wax of WT, 35S:PaMYB9A1, and 35S:PaMYB9A2 mature leaves. A, Total wax loads and (B) cuticular wax compositions of WT, 35S:PaMYB9A1, and 35S:PaMYB9A2 mature leaves. Chloroform extracted lipids were analyzed by GC–MS after TMS derivatization. Total wax loads were normalized to the internal control (Heptatriacontane). Each wax constituent amount is normalized against the total area and shown as amounts relative to the total wax load. Chain lengths are labeled on the horizontal axis (C21–C54). Error bars represent the mean ± sd from three biological replicates. Asterisks indicate significant differences in each wax constituent between WT and 35S:PaMYB9A1 plants as well as 35S:PaMYB9A2 plants using Student’s t test (*P < 0.05 and **P < 0.01).

Overexpression of either PaMYB9A1 or PaMYB9A2 promotes cuticular wax deposition on P. aphrodite petals

To examine whether PaMYB9A1 and PaMYB9A2 also function to regulate the biosynthesis of cuticular wax in Phalaenopsis petals, as observed in transgenic tobacco leaves, we transiently overexpressed these genes in petals and observed the petal adaxial surfaces. Seven days postagroinfiltration both PaMYB9A1- and PaMYB9A2-overexpressing petals exhibited a small number of wax crystals and many circular-shaped deposits on the surface of their conical epidermal cells. These deposits were found for three independent lines of each construct and were not observed in control plants (Figure 5). On the other hand, the cuticular wax loads and compositions of the vector control, PaMYB9A1- and PaMYB9A2-overexpressing petals were determined through GC–MS analysis with three biological replicates. Compared with the vector control, the GC–MS results showed that both PaMYB9A1- and PaMYB9A2-overexpressing petals had higher total wax loads, displaying increases of 314 ± 48% and 212 ± 40%, respectively (Figure 6A). In addition, the proportions of n-alkanes of ≥C32 and C30 primary alcohols are significantly increased, with the corresponding decrease in C25, C26, and C27 primary alcohols (Figure 6B). These results revealed that both PaMYB9A1 and PaMYB9A2 act to facilitate the biosynthesis of cuticular wax in petal conical cells.

Figure 5

Transient overexpression of PaMYB9A1 and PaMYB9A2 facilitates the deposition of cuticular wax on the surface of petal conical epidermal cells. A–F, Scanning electron micrograph of conical epidermal cells present on the adaxial face of (A and D) vector control, and (B and E) PaMYB9A1-, and (C and F) PaMYB9A2-overexpressing petals. Bar = 50 µm (A–C) and 5 µm (D–F).

Figure 6

Cuticular wax analyses of vector control-containing and PaMYB9A1- and PaMYB9A2-overexpressing petals. A, Total wax loads and (B) cuticular wax compositions of vector control-containing and PaMYB9A1- and PaMYB9A2-overexpressing petals. Chloroform-extracted lipids were analyzed by GC–MS after TMS derivatization. Total wax loads are normalized to the vector control. Each wax constituent amount is normalized against the total area, and amounts relative to the total wax load are shown. Chain lengths are labeled on the horizontal axis (C22–C54). Error bars represent the sd from three biological replicates. Asterisks indicate significant differences in each wax constituent between vector controls and PaMYB9A1-overexpressing petals and PaMYB9A2-overexpressing petals using Student’s t test (*P < 0.05 and **P < 0.01).

To identify putative targets of PaMYB9A1 and PaMYB9A2, we performed a comparative transcriptome analysis between vector control and PaMYB9A1-overexpressing petals as well as PaMYB9A2-overexpressing petals to screen genes whose expression was significantly altered in response to the activity of PaMYB9A1 and PaMYB9A2. We extracted total RNA from these petals, and RNA extracted from three independent biological replicates was used for RNA sequencing. Both the results of transcriptome analysis and RT-qPCR assays confirmed that PaMYB9A1 and PaMYB9A2 were highly expressed in their corresponding transiently overexpressed petals compared with the vector control (Supplemental Figure S11, A and C). DESeq2 analyses showed that 188 and 85 up- and downregulated DEGs (log2-fold change ≥1, P < 0.01, and q-value < 0.01) were identified in PaMYB9A1-overexpressing petals, respectively, while 81 and 135 up- and downregulated DEGs were identified in PaMYB9A2-overexpressing petals, respectively (Supplemental Data Sets S1 and S2; Supplemental Figures S12 and S13). Among the upregulated DEGs, a gene encoding 3-ketoacyl-CoA synthase (KCS; PAXXG005280), which is the key enzyme involved in VLCFA biosynthesis, was identified in PaMYB9A1-overexpressing petals, and a gene encoding FAR (PAXXG145540), which is a key enzyme involved in fatty alcohol synthesis, was identified in PaMYB9A2-overexpressing petals (Supplemental Data Sets S1 and S2).

Overexpression of either PaMYB9A1-SRDX or PaMYB9A2-SRDX represses the differentiation of conical epidermal cells in petals

The overexpression of TFs containing the SRDX domain, by Chimeric REpressor gene Silencing Technology dominantly suppresses the transcription of their regulatory targets and results in a phenotype similar to that of loss-of-function mutants (Mitsuda et al., 2011). To evaluate the functional roles of PaMYB9A1 and PaMYB9A2 in conical epidermal cell differentiation in P. aphrodite petals, we constructed the chimeric repressors of PaMYB9A1 and PaMYB9A2 by fusing the SRDX repression domain to the 3′-end of their coding sequence (Supplemental Figure S8). Agrobacterium tumefaciens cells containing the empty vector or PaMYB9A1-SRDX- and PaMYB9A2-SRDX-overexpressing constructs were infiltrated into the petals at Stage 3, when conical epidermal cell differentiation begins. We observed the adaxial surfaces of the infiltrated petals at 7 d post inoculation (dpi) by SEM. In the control petals, typical conical epidermal cells could be observed on the adaxial epidermis (Figure 7A). In contrast, both the PaMYB9A1-SRDX- and PaMYB9A2-SRDX-overexpressing petals exhibited flattened epidermal cells on their adaxial face (Figure 7, B and C), suggesting that both PaMYB9A1 and PaMYB9A2 act to regulate the differentiation of petal conical epidermal cells.

Figure 7

Epidermal cell phenotypes of vector control-containing and PaMYB9A1-SRDX- and PaMYB9A2-SRDX-overexpressing petals. A, Scanning electron micrograph of conical epidermal cells present on the adaxial epidermis of the vector control petal. B and C, Scanning electron micrograph of flattened epidermal cells present on the adaxial epidermis of (B) PaMYB9A1-SRDX- and (C) PaMYB9A2-SRDX-overexpressing petals. Bar = 100 µm.

We further identified the putative targets of PaMYB9A1-SRDX and PaMYB9A2-SRDX by the same strategy described for PaMYB9A1/2 above. We first confirmed that PaMYB9A1-SRDX and PaMYB9A2-SRDX were highly expressed in the petals with transient overexpression at Stage 4 (Supplemental Figure S11, B and D). DESeq2 analyses revealed 170 and 151 up- and downregulated DEGs (log2-fold change ≥1, P < 0.01, and q-value < 0.01), respectively, in PaMYB9A1-SRDX-overexpressing petals (Supplemental Data Set S3; Supplemental Figures S12 and S13) and 913 and 664 up- and downregulated DEGs in PaMYB9A2-SRDX-overexpressed petals (Supplemental Data Set S4; Supplemental Figures S12 and S13). Among these DEGs, 14 and 23 cell wall-associated DEGs regulated by PaMYB9A1-SRDX and PaMYB9A2-SRDX, respectively, were identified (Table 1). Six DEGs, including genes encoding two 36.4-kDa proline-rich proteins, a microtubule-associated protein, a cellulose synthase-like protein, the actin cytoskeleton-regulatory complex protein PAN-1-like, and wall-associated receptor kinase 2-like, were significantly downregulated in both PaMYB9A1-SRDX- and PaMYB9A2-SRDX-overexpressing petals. In addition, three pectinesterase genes are down-regulated in PaMYB9A2-SRDX-overexpressed petals. These results suggest that PaMYB9A1 and PaMYB9A2 cooperatively control conical epidermal cell development by regulating cell wall-related genes.

Table 1

Cell wall-related down-regulated DEGs in PaMYB9A1-SRDX- and PaMYB9A2-SRDX-overexpressing petals compared with the vector control

Gene ID	log2 FC	P-value	Q-value	Annotation
PaMYB9A1-SRDX
PAXXG019470a	−1.09	9.57E-153	2.74E-151	36.4 kDa proline-rich protein-like
PAXXG237140	−1.21	1.15E-74	1.87E-73	36.4 kDa proline-rich protein
PAXXG190220a	−1.58	1.16E-06	3.02E-06	36.4 kDa proline-rich protein-like
PAXXG201700	−1.72	0	0	Microtubule-associated protein TORTIFOLIA1-like
PAXXG367790a	−1.12	0	0	65-kDa microtubule-associated protein 1-like
PAXXG003800a	−1.04	8.90E-34	8.30E-35	Cellulose synthase-like protein D4
PAXXG315280	−2.09	4.18E-45	4.68E-44	Vegetative cell wall protein gp1-like
PAXXG315300	−1.58	7.29E-12	2.91E-11	Vegetative cell wall protein gp1-like
PAXXG387100	−1.15	1.09E-41	1.15E-40	Xyloglucan endotransglucosylase/hydrolase
PAXXG200250	−1.2	2.77E-44	3.06E-43	Probable xyloglucan endotransglucosylase/hydrolase
PAXXG325960a	−1.23	1.20E-17	6.38E-17	Actin cytoskeleton-regulatory complex protein PAN1-like
PAXXG150460a	−1.25	3.50E-19	1.99E-18	Wall-associated receptor kinase 2-like
PAXXG282620	−.1	6.41E-17	3.30E-16	Polygalacturonase-like
PAXXG033430	−1.24	1.85E-08	5.77E-08	Altered xyloglucan 4-like
PaMYB9A2-SRDX
PAXXG319070	−1.58	0	0	Probable pectinesterase/pectinesterase inhibitor 34
PAXXG325960	−3.04	7.48E-54	2.47E-53	Actin cytoskeleton-regulatory complex protein PAN1-like
PAXXG367790	−1.05	0	0	65-kDa microtubule-associated protein 1-like
PAXXG317750	−1.39	3.47E-55	1.17E-54	65-kDa microtubule-associated protein 3-like
PAXXG019470	−2.95	0	0	36.4 kDa proline-rich protein-like
PAXXG190220	−1.07	0	0	36.4 kDa proline-rich protein-like
PAXXG003800	−1.46	4.93E-60	1.77E-59	Cellulose synthase-like protein D4
PAXXG156700	−1.79	0	0	Cellulose synthase-like protein D3
PAXXG156710	−1.01	1.38E-09	1.71E-09	Cellulose synthase-like protein D3
PAXXG260600	−1.31	0	0	Probable cellulose synthase A catalytic subunit 5
PAXXG018510	−2.04	0	0	Probable xyloglucan endotransglucosylase/hydrolase protein 23
PAXXG177360	−1.54	3.60E-13	5.24E-13	Probable xyloglucan endotransglucosylase/hydrolase
PAXXG177370	−1.91	1.22E-16	2.02E-16	Probable xyloglucan endotransglucosylase/hydrolase
PAXXG220520	−1.06	0	0	Probable xyloglucan endotransglucosylase/hydrolase protein 8
PAXXG065920	−1.61	1.30E-55	4.37E-55	Probable xyloglucan endotransglucosylase/hydrolase protein 30
PAXXG069840	−3.64	0	0	Probable xyloglucan endotransglucosylase/hydrolase protein 23
PAXXG069850	−3.45	0	0	Probable xyloglucan endotransglucosylase/hydrolase protein 23
PAXXG132100	−2.52	3.27E-52	1.06E-51	Probable xyloglucan endotransglucosylase/hydrolase protein 33
PAXXG018520	−5.07	0	0	Xyloglucan endotransglycosylase
PAXXG248380	−1.51	0	0	Probable xyloglucan glycosyltransferase 9
PAXXG153360	−1.08	1.77E-18	3.10E-18	Glucan endo-1,3-beta-glucosidase 11-like isoform X3
PAXXG226440	−1.19	0	0	Glucan endo-1,3-beta-glucosidase 14-like isoform X1
PAXXG033550	−1.75	1.37E-91	6.36E-91	Glucan endo-1,3-beta-glucosidase 10-like
PAXXG187670	−2.03	3.43E-113	1.82E-112	Endoglucanase 12-like
PAXXG033420	−1.2	4.24E-35	1.08E-34	Protein altered xyloglucan 4-like isoform X1
PAXXG093590	−1.36	2.58E-33	6.36E-33	Pectinesterase
PAXXG315920	−2.84	0	0	Pectinesterase-like
PAXXG029450	−4.22	3.27E-282	3.13E-281	Pectinesterase-like
PAXXG040880	−1.4	0	0	Probable pectate lyase 5

DEGs that are also down-regulated in PaMYB9A2-SRDX-overexpressing petals.

Discussion

In dicots, the functions of subgroup 9A R2R3-MYB gene family members, the MIXTA and MIXTA-like genes, regulate the differentiation of various epidermal cells. However, their biological function has yet to be confirmed in monocots. In this study, we validated the functions of two subgroup 9A R2R3-MYB genes, PaMYB9A1 and PaMYB9A2, from P. aphrodite through the transient overexpression of their corresponding chimeric repressors, PaMYB9A1-SRDX and PaMYB9A2-SRDX, in developing petals. Our data showed that both PaMYB9A1-SRDX and PaMYB9A2-SRDX appear to repress the differentiation of conical epidermal cells in petals. These results revealed that the function of the subgroup 9A R2R3-MYB gene family members, which control epidermal cell differentiation, is conserved between monocots and dicots.

PaMYB9A1 and PaMYB9A2 have evolved additional functions in cuticular wax biosynthesis in P. aphrodite petals

Several subgroup 9A R2R3-MYB TFs possessing functions other than controlling epidermal cell fate have been reported. For example, in A. thaliana, the myb106-2 mutant shows slightly flattened petal epidermal cells without nanoridge decoration, suggesting that MYB106 plays a principal role in regulating epidermal cell differentiation and cutin nanoridge decoration in petals (Oshima et al., 2013). In addition, in tomato, SlMIXTA-like controls cutin biosynthesis and epidermal cell patterning in the pericarp (Lashbrooke et al., 2015); in Mimulus lewisii, a MIXTA-like gene is involved in epidermal cell development and carotenoid pigmentation in flowers (Yuan et al., 2013). Here, we confirmed that both PaMYB9A1/2 coordinate conical epidermal cell differentiation and cuticular wax biosynthesis in P. aphrodite petals.

Chemical composition of cuticular wax in shiny leaves

The overexpression of both PaMYB9A1 and PaMYB9A2 in tobacco caused a shiny phenotype on the abaxial surface of mature leaves. However, unlike overexpression in P. aphrodite petals, there were no apparent epicuticular wax crystals deposited on the epidermal cells of the shiny tobacco leaves, and the wax load was not distinctly altered on the shiny leaves (Figures 4–6). These results revealed that the shiny phenotype was caused by an altered cuticular wax composition in 35S:PaMYB9A1/2 tobacco leaves and not by increased wax deposition as observed in A. thaliana.

The cuticular wax composition of P. aphrodite petals reveals an unknown cuticular wax biosynthetic pathway

The cuticular wax composition varies in different plants and different tissues/organs across the same plant. In this study, we identified the chemical composition of cuticular wax deposited on P. aphrodite petals for the first time; it is mainly composed of two VLCFA derivatives, n-alkanes (the major constituents) and primary alcohols. In contrast, the cuticular wax of A. thaliana flowers is composed of fatty acids, primary and secondary alcohols, n-alkanes, branch-alkanes, and ketones (Shi et al., 2011). The P. aphrodite petals have a more simplified chemical composition of cuticular wax than A. thaliana flowers.

The well-known biosynthetic pathways of cuticular wax and the alcohol-forming and alkane-forming pathways, primarily derived from previous studies in A. thaliana. In these studies, even-chain VLCFAs were used as precursors to produce even-chain primary alcohols and wax esters through the alcohol-forming pathway, and to produce even-chain aldehydes and odd-chain n-alkanes, secondary alcohols, and ketones through the alkane-forming pathway (Seo et al., 2011). The two biosynthetic pathways can fully explain how the VLCFA-derived constituents of cuticular wax are synthesized in A. thaliana stems. However, the cuticular wax of P. aphrodite petals contains large proportions of even-chain n-alkanes (70.4%) and odd-chain primary alcohols (26.1%) as the main constituents. The cuticular wax on A. thaliana flowers also contains small amounts of even-chain alkanes as the constituents (Shi et al., 2011). These observations suggest that an unknown biosynthetic pathway is responsible for producing even-chain alkanes and odd-chain primary alcohols in flowers. In A. thaliana, an odd-chain VLCFA (C29:1 fatty acid) is one of the constituents of cuticular wax deposited on flowers, suggesting that odd-chain VLCFAs could be generated in flowers (Shi et al., 2011). Therefore, one possibility is that odd-chain VLCFAs also serve as precursors to generate even-chain alkanes and odd-chain primary alcohols through well-known wax biosynthetic pathways in flowers. However, the biosynthetic pathway for producing even-chain alkanes, odd-chain alcohols, and odd-chain VLCFAs in flowers remains unclear and requires further investigation.

PaSHN is not a downstream target of PaMYB9A1 or PaMYB9A2

In A. thaliana, the MIXTA-like gene MYB106 acts as an activator of WIN1/SHN1 to regulate the differentiation of conical epidermal cells and the decoration of cutin nanoridges in petals (Aharoni et al., 2004; Shi et al., 2011; Oshima et al., 2013). We found that only the WIN1/SHN1 homologous gene PaSHN (PAXXG128820) was present in the P. aphrodite genome and was predominantly expressed in petals at developmental Stages 2–4 (Supplemental Figure S14). However, based on the comparative transcriptome analysis, the expression levels of PaSHN were not altered significantly in PaMYB9A1-, PaMYB9A2-, PaMYB9A1-SRDX-, and PaMYB9A2-SRDX-overexpressing petals compared to vector controls, suggesting that, unlike in A. thaliana, PaMYB9A1 and PaMYB9A2 are not upstream regulators of PaSHN. The results are consistent with petal conical epidermal cells without nanoridge decoration in P. aphrodite.

Wax synthesis and cell wall-related genes are putative direct targets of PaMYB9A1

The application of chromatin immunoprecipitation (ChIP) followed by PCR- or sequencing-based methods is a well-accepted approach for the high-throughput screening of direct interactions between DNA and proteins to identify the putative targets of TFs (Schmidt et al., 2009). However, it is challenging to avoid mistargeted reactions (false-positive results). For example, 4,282 putative target genes of the MADS-box TF SEPALLATA3 have been identified (false discovery rate (FDR) < 0.001) by ChIP followed by ultrahigh-throughput Solexa (Illumina, San Diego, CA, USA) sequencing (Kaufmann et al., 2009). The 4,282 genes account for approximately one-sixth of the total number of genes in A. thaliana, suggesting the potential inclusion of many mistargeted false-positive results. In this study, we tried to identify the putative direct targets of PaMYB9A1/2 by a combined strategy of comparative transcriptome analysis and detailed temporal expression analysis of genes at each petal developmental stage.

We hypothesized that (1) the expression levels of direct targets of PaMYB9A1/2 would be significantly increased in PaMYB9A1/2-overexpressing petals and would be significantly decreased in PaMYB9As-SRDX-overexpressing petals and that (2) PaMYB9A1/2 and their direct targets should have similar temporal expression patterns during petal development. Therefore, we assessed the levels of genes expressed in PaMYB9A1-, PaMYB9A1-, PaMYB9A1-SRDX-, and PaMYB9A2-SRDX-overexpressing petals and compared them to those of the vector control to identify DEGs by DESeq2 analysis. However, according to our hypothesis, 35 and 14 targets of PaMYB9A1 and PaMYB9A2, respectively, met condition 1 above, and only 6 of 35 targets of PaMYB9A1 met condition 2 as well; the number of candidate direct PaMYB9A1 targets that are up-regulated in PaMYB9A1-overexpressing petals and downregulated in PaMYB9A1-SRDX-overexpressing petals was reasonably narrowed down to 35 genes (Supplemental Figure S15). Among these genes, 6 genes (17.1%; 6/35) exhibited an expression pattern similar to PaMYB9A1 and were also classified into the yellow module by WGCNA (Figure 8). Among these putative direct targets of PaMYB9A1, one gene (PAXXG005280) encodes a KCS 11-like protein that is involved in the wax biosynthetic pathway (Lee et al., 2009). In addition, two genes, PAXXG201700 and PAXXG201700, encode the microtubule-associated protein TORTIFOLIA1-like and cellulose synthase-like protein D4, which likely function in cell wall remodeling (Rajangam et al., 2008). However, the other 29 putative target genes (82.9%; 29/35) exhibited different temporal expression patterns than PaMYB9A1 during petal growth, suggesting that these genes are likely false-positive results (Figure 8). These results revealed that this strategy is efficient in screening out the direct targets of TFs from the putative targets identified by comparative transcriptome analysis.

Figure 8

Temporal expression patterns of putative direct target genes of PaMYB9A1 during petal development. The heatmap shows the expression patterns of PaMYB9A1 and 30 putative direct-target gene candidates of PaMYB9A1 with FPKM > 10 at least in one developmental stage. Data are from the three biological replicates for each developmental stage. Asterisks indicate the genes that are classified into the yellow module that includes PaMYB9A1 and its putative direct target genes. The color scale indicates the expression levels (log2-fold change) relative to Stage 1 based on FPKM values for upregulated genes (red), genes that show no expression difference (pale yellow), and downregulated genes (blue).

Materials and methods

Plant materials and growth conditions

Phalaenopsis aphrodite subsp. formosana seedlings were purchased from Ciean Po International Co., Ltd. (Ping Tung, Taiwan). All P. aphrodite plants were grown at a high temperature of 32°C/30°C with a 12-h photoperiod to maintain vegetative growth. The plants were shifted to a low ambient temperature of 24°C/18°C with a 12-h photoperiod to induce spiking. After A. tumefaciens infiltration, the P. aphrodite plants were grown at a constant temperature of 25°C with a 12-h photoperiod. On the other hand, the Nicotiana tabacum L. cv W38 plants were grown at 25°C/20°C with a 12-h photoperiod.

RNA isolation and RT reactions

For P. aphrodite, samples were collected and ground in liquid nitrogen. RNA was extracted using TEN buffer (7.8 mM Tris, 15.6 mM Na2EDTA, 156 mM NaCl, 4% w/v SDS (sodium dodecyl sulfate), 16 mM DTT, pH 9.0) following phenol–chloroform extraction and precipitated, at 13,000 rpm for 40 min at 4°C using a final concentration of 2 M LiCl after overnight incubation at 4°C. For tobacco, total RNA from leaf samples was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. A total of 1 μg of RNA was reverse-transcribed in a 30-μL RT reaction using the First-strand cDNA Synthesis kit (Promega, Madison, WI, USA) and oligo-dT primer, followed by the addition of 70 μL ddH2O to the synthesized cDNAs.

RT-qPCR

RT-qPCR was carried out with an ABI 7500HT sequence detection system (Applied Biosystems, Foster City, CA, USA) using SYBR Green Master Mix (Applied Biosystems) following the manufacturer’s instructions. Twenty microliters of the RT-qPCR reaction contained 1 μL cDNA, 0.25 mM primers, and 10 μL of KAPA SYBR FAST master mix (KAPA Biosystems, Wilmington, MA, USA). RT-qPCR was performed using triplicate technical replicates, and the experiments were repeated with RNA isolated from three independent biological samples. Relative quantification analysis was carried out using ABI SDS version 1.4 software (Applied Biosystems). Transcript levels were normalized to the endogenous control gene Ubi10 using the comparative cycle threshold method (Lu et al., 2007). The primers used for RT-qPCR are listed in Supplementary Table S1.

Construction of the PaMYB9A1, PaMYB9A2, PaMYB9A1-SRDX, and PaMYB9A2-SRDX overexpression vectors

The PaMYB9A1 and PaMYB9A2 coding sequences were amplified using the primer pair attb1-PaMYB9A1F/attb2PaMYB9A1R and attb1-PaMYB9A2F/attb2- PaMYB9A2R (Supplemental Table S1). The PCR fragments were then cloned into pDONR221 using the Gateway BP reaction kit (Thermo Fisher Scientific, Carlsbad, CA, USA) according to the cloning protocol to generate pDONR221-PaMYB9A1 and pDONR221-PaMYB9A2. pH2GW7 was incubated with pDONR221-PaMYB9A1 and pDONR221-PaMYB9A2 using the Gateway LR reaction kit (Thermo Fisher Scientific, Carlsbad, CA, USA) according to the cloning protocol to generate pH2GW7-PaMYB9A1 and pH2GW7-PaMYB9A2. pH2GW7-PaMYB9A1-SRDX and pH2GW7-PaMYB9A2-SRDX were constructed using the same process described for pH2GW7-PaMYB9A1 and pH2GW7-PaMYB9A2 except that the primers attb2PaMYB9A1R and attb2PaMYB9A2R were used instead of attb2PaMYB9A1R-SRDX and attb2PaMYB9A2R-SRDX, respectively.

Nicotiana tabacum L. cv W38 transformation

Constructs pH2GW7-PaMYB9A1 and pH2GW7-PaMYB9A2 were transformed into A. tumefaciens LBA4404 by electroporation and grown on YEB (10 g/L yeast extract, 10 g/L peptones, 5 g/L NaCl and 15 g/L agar; pH 7.0) plates with 50 µg/mL rifampicin and 100 µg/mL spectinomycin. Tobacco leaf discs were transformed according to the method developed by Horsch et al. (1985).

Agroinfiltration

Agroinfiltration was carried out as previously described (Lu et al., 2012) with some modifications. Briefly, A. tumefaciens C58C1 (pTiB6S3ΔT)H competent cells were transformed with overexpression constructs using an electroporation system (Bio-Rad Laboratories, Hercules, CA, USA). Then, A. tumefaciens was incubated in 10 mL YEB at 28°C until the optical density at 600 nm (OD600) reached 0.8–1.0. After centrifugation (4,000 rpm, 10 min, in room temperature), A. tumefaciens was suspended in 20 mL AB-MES medium (17.2 mM K2HPO4, 8.3 mM NaH2PO4, 18.7 mM NH4Cl, 2 mM KCl, 1.25 mM MgSO4, 100 μM CaCl2, 10 μM FeSO4, 50 mM MES, 2% glucose (w/v), and pH 5.5) with 200 μm acetosyringone and cultured with shaking at 28°C overnight. The overnight culture was centrifuged (4,000 rpm, 10 min, at room temperature), the supernatant was removed, and the pellet of A. tumefaciens was suspended in 2 mL of infiltration medium containing 50% MS medium (1/2 MS salt supplemented with 0.5% sucrose (w/v), pH 5.5), 50% AB-MES and 200 μm acetosyringone. Suspensions of A. tumefaciens cells showing OD600 values of 0.4 and 0.8 were infiltrated into Stages 3 and 4 petals, respectively, by using insulin syringes (60 mm × 31 G).

Phylogenetic analysis

The amino acid sequences of all R2R3-MYB family genes in P. aphrodite and well-characterized MIXTA and MIXTA-like genes were first aligned using ClustalW, and phylogenetic analysis was performed using the maximum likelihood method with the JTT model by MEGA version 10 (Kumar et al., 2018). Bootstrap values were calculated with 1,000 replicates.

Cuticular wax analysis

Cuticular waxes were extracted and analyzed as previously described with slight modifications (Kurdyukov et al., 2006). Plant materials (5 × 10 cm2 of tobacco leaf or 3 g fresh weight of P. aphrodite petal) were cut and immersed in chloroform for 30 s at room temperature. The resulting solution of cuticular waxes was spiked with 10 mg of heptatriacontane (Fluka Chemie AG, Buchs, Switzerland) as an internal standard. Chloroform soluble extracts were then evaporated under a continuous gas nitrogen flow. Compounds were derivatized by the addition of 100 µL of a pyridine: N,O-bis-trimethylsilyl-trifluoroacetamide (Macherey-Nagel, Allentown, PA, USA) solution (1:1, v:v) for 20 min at 100°C before GC–MS analysis. The composition was analyzed using a capillary gas chromatograph with a DB-5 MS column (length 30 m, inner diameter 0.25 mm, film thickness 0.25 μm, Agilent, Santa Clara, CA, USA), with helium applied at a linear velocity of 0.69 mL/min as the carrier gas, and detected on a mass spectrometric detector (70 eV, m/z 50–700, GCMSQP2020, Shimazu, Kyoto, Japan). The GC–MS program was as follows: injection at 250°C, holding for 5 min, an increase of 3 µL/min to 290°C, holding for 10 min, an increase of 2 µL/min to 300°C, and holding for 10 min. Total wax loads were normalized to the internal control (Heptatriacontane). Each wax constituent amount is normalized against the total area and shown as amounts relative to the total wax load. Statistical significance of the measurements, from three biological replicates, was determined using a Student’s t test.

WGCNA

Three biological replicates were collected from developmental Stages 1–5 petals and their total RNA was extracted for RNA sequencing. RNA sequencing was entrusted to BGI company, and paired-end sequencing was carried out with a BGISeq400 system at Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China). The sample datasets were filtered to remove low abundance genes that did not reach a maximum FPKM value of 2 in all samples and a value of 0 in at least 13 samples. FPKM values (normalized) of transcripts obtained from the fifteen samples were used to generate the coexpression networks using the WGCNA (version 1.47) package. An adjacency matrix was constructed using a soft threshold power of 9. Network interconnectedness was measured by calculating the topological overlap using the TOMdist function with a signed TOM Type. Average hierarchical clustering using the hclust function was performed to group the genes based on the topological overlap dissimilarity measure (1-TOM) of their connection strengths. Network modules were identified using a dynamic tree cut algorithm with a minimum cluster size of 50 and merging threshold function at 0.25. To identify biologically significant modules, the eigengene (first principal component) for each module was plotted using ggplot2 in R. To identify hub genes of the modules, the module membership (MM) for each gene was calculated based on the Pearson correlation between the expression level and the module eigengene. Genes within the module with the highest MM were highly connected within that module.

SEM

For quick SEM, samples were frozen immediately in liquid nitrogen and placed on the observation platform. After vacuum treatment, the samples were observed by SEM (Hitachi TM3000, Tokyo, Japan). For cryo-SEM, the samples were dissected and rapidly frozen in the cryo chamber (Quorum PP2000TR) at −210°C for at least 15 s. They were then transferred to the preparation chamber of a cryo-SEM system (Quanta 200, FEI), where the frozen petals were sputter-coated with gold (2 min at −130°C) and then transferred to the observation chamber and observed at −190°C.

Comparative transcriptome analysis

Three biological samples were collected from vector control, PaMYB9A1- and PaMYB9A2-overexpressing petals, and their total RNA was extracted for RNA sequencing. RNA sequencing was entrusted to BGI, and paired-end sequencing was carried out with a BGISeq400 system at Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China). After sequencing, low-quality raw reads (containing >5% unknown nucleotides and >50% low-quality (q-value ≤ 10) bases) were filtered by SOAPnuke. The high-quality clean reads were mapped to the P. aphrodite genome (Chao et al., 2018) by Hierarchical Indexing for Spliced Alignment of Transcripts. Gene expression was calculated and normalized to FPKM (Fragments Per Kilobase per Million) values. DESeq2 was used to identify the DEGs (log2fold change >1 and P ≤ 0.001, FDR (q-value) < 0.01) through comparison with vector control. Enrichment analysis of the GO and KEGG of the DEGs was analyzed by phyper in R.

Sequence read archive numbers

Sequencing data used in this article can be found in the National Center for Biotechnology Information (NCBI) databases (https://www.ncbi.nlm.nih.gov/) under the following BioProject numbers: PRJNA637697, PRJNA728444, PRJNA728478. The genome data of P. aphrodite can be found in Chao et al. (2018).

Accession numbers

Sequence data from this article can be found in the Orchidstra databases (http://orchidstra2.abrc.sinica.edu.tw/orchidstra2/index.php) under the following accession numbers: PaMYB9A1 (PAXXG029600) and PaMYB9A2 (PAXXG123420).

Supplemental data

The following supplemental materials are available in the online version of this article.

The P. aphrodite flower.

Phylogenetic tree of P. aphrodite R2R3-MYB proteins.

Expression profiles of PaMYB9A1 and PaMYB9A2 genes in different organs.

Scatter plots depicting differentially expressed genes (DEGs).

Gene ontology (GO) analysis.

KEGG pathway enrichment analysis.

Hierarchical clustering tree showing coexpression modules classified by WGCNA.

Diagram of the constructs used in this study.

Expression of PaMYB9A1 and PaMYB9A2 in putative transgenic tobacco plants.

SEM images of the leaf epidermis of WT, 35S:PaMYB9A1, and 35S:PaMYB9A2 transgenic tobacco plants.

Expression profiles of PaMYB9A1 and PaMYB9A2 in different agroinfiltrated petals.

Gene ontology (GO) analysis.

KEGG pathway enrichment analysis.

Expression profile of PaSHN.

Comparative transcriptome analysis of empty vector-containing and PaMYB9A1-, and PaMYB9A1-SRDX-overexpressing petals.

List of primers used in this study.

List of DEGs in PaMYB9A1-overexpressing petals.

List of DEGs in PaMYB9A2-overexpressing petals.

List of DEGs in PaMYB9A1-SRDX-overexpressing petals.

List of DEGs in PaMYB9A2-SRDX-overexpressing petals.

Click here for additional data file.