Identification of the Key Regulatory Genes Involved in Elaborate Petal Development and Specialized Character Formation in Nigella damascena (Ranunculaceae)

Extensive transcriptomic and functional analyses of elaborate petal development and specialized character formation in Nigella damascena point to the mechanisms underlying lateral organ diversification in plants.


INTRODUCTION
Petals are highly specialized lateral organs that lie between the outer sepals and inner stamens and/or carpels in a flower with a dimorphic perianth and play key roles in plant-animal interactions. Petals display spectacular arrays of diversity in size, shape, color, and function and are increasingly being used as a model system for the study of plant organogenesis and evolution (Irish, 2008;Huang and Irish, 2016). In Arabidopsis (Arabidopsis thaliana; Brassicaceae) and many other petalous taxa (such as snapdragon [Antirrhinum majus], petunia [Petunia hybrida], columbine [Aquilegia], peonies [Paeonia], roses [Rosa], lilies [Lilium], and orchids [Orchidaceae]), petals, like any other types of lateral organs of plants, are formed through at least three highly conserved developmental processes: initiation, growth, and maturation. Genes, pathways, and networks involved in the proliferation, expansion, and differentiation of cells, as well as those specifying the adaxial-abaxial, proximal-distal, and lateral-medial polarities of the entire organ, are therefore indispensable (Irish, 2008;Huang and Irish, 2016;Walcher-Chevillet and Kramer, 2016;Shan et al., 2019). Yet, unlike many other types of lateral organs, petals are usually brightly colored and/or unusually shaped, suggesting that, in addition to the common themes of lateral organs, there are developmental processes that are specific to petals (such as determination of petal identity; generation of the sometimes highly specialized, three-dimensional structures; and/or formation of intriguing color patterns; Irish, 2008Irish, , 2017Huang and Irish, 2016;Moyroud and Glover, 2017). Hence, identifying the regulatory genes, pathways, and networks involved in petal morphogenesis is not only the prerequisite for elucidating how petals were made through development and evolution but also is helpful for understanding the mechanisms underlying lateral organ diversification.
Notably, almost all our knowledge about the mechanisms underlying petal organogenesis was gained from studying species with simple petals (such as those of Arabidopsis). In flowering plants, however, there are many species whose petals are highly specialized, either having lobes, teeth, fringes, or appendages along their margins, or possessing scales, spurs, pockets, or other types of modifications on their adaxial/abaxial sides, or both (Figures 1A to 1H;Endress and Matthews, 2006). These petals, also known as elaborate petals, have been recorded in hundreds of families of at least 23 orders of angiosperms and are generally believed to have played key roles in the adaptive evolution of the corresponding lineages (Endress and Matthews, 2006;Yao et al., 2019). Interestingly, in addition to the modification of their basic structures through ventral or marginal elaboration ( Figure 1I), many elaborate petals also bear some newly gained characters (such as the labellum of the orchid Ophrys speculum, which mimics a female wasp [Colpa aurea]; Endress and Matthews, 2006;Yao et al., 2019). Uncovering the developmental mechanisms underlying the organogenesis of elaborate petals therefore will not only pave the way to elucidating the mechanisms underlying petal elaboration but also help in understanding the general patterns of new character origination and flowering plant evolution.
The genus Nigella (Ranunculaceae) appears to be an excellent system for the study of petal elaboration and its underlying mechanisms, for four reasons. First, petals of several species in this genus have long been regarded as being elaborate, possessing not only a rod-like stalk, a laminar upper lip, and a bifurcate lower lip but also several highly specialized morphological characters, such as long hairs, short trichomes, a nectary, and a pair of pseudonectaries (Figures 1J to 1L;Strid, 1970;Tamura, 1995;Zhao et al., 2011). This attribute makes it possible to study both marginal/ventral elaboration and new character origination. Second, within the genus, the complexity of the elaborate petals shows considerable variation among species so that the evolutionary histories of many features can be inferred once a reliable phylogenetic tree is available (Zohary, 1983;Bittkau and Comes, 2005;Jaros et al., 2018;Yao et al., 2019). Third, in addition to elaborate petals, species of Nigella also possess petaloid sepals that are almost indistinguishable from the simple petals of many other species (including Arabidopsis) in both morphology and function ( Figure 1J). This juxtaposition allows us to compare the elaborate petals with simple petals within the same flower. Fourth, almost all species of the genus are easy to collect, cultivate, and manipulate for which virus-induced gene silencing (VIGS) and many other technologies are applicable Liao et al., 2020).
Recently, using species of Nigella as materials, some interesting progress has been made to understand the process, pattern, and mechanisms of petal elaboration. For example, it has been shown that, like all other types of lateral organs, the highly specialized elaborate petals of N. damascena have a very simple starting point-a dome-like primordium. During development, as the upper lip, long hairs, pseudonectaries, ridges, and short trichomes were initiated successively, the complexity of the petal increased gradually . In addition, when the ancestral states of important characters were reconstructed, it becomes evident that the evolutionary elaboration of Nigella petals was a gradual process, involving not only modifications of pre-existing structures but also originations of new characters  (NidaAGL6), and NidaSEPALLATA1 (NidaSEP1), NidaSEP2, and NidaSEP3 (Gonçalves et al., 2013;Zhang et al., 2013;Wang et al., 2015), as well as those involved in pseudonectary development, such as NidaYABBY5 (NidaYAB5; Liao et al., 2020), have been identified and functionally characterized. Genes controlling the formation of many other characters, however, remain largely unknown.
Here, using N. damascena as a model, we conducted a detailed time-course transcriptomic analysis of petals to identify the genes and programs that are potentially involved in elaborate petal development. By correlating the expression pattern of gene with the morphogenetic events during petal development, we successfully identified 30 candidate genes associated with the formation of key morphological characters. Functional analyses further confirmed that a class I homeodomain-leucine zipper (HD-Zip) family transcription factor gene, Nigella damascena LATE MERISTEM IDENTITY1 (NidaLMI1), plays important roles in the development of short trichomes and bifurcation of the lower lip. Our results not only elucidate the differences between simple and elaborate petals in molecular bases but also pave the way to a comprehensive understanding of the mechanisms underlying lateral organ diversification and the new character origination.

Reference Transcriptome and DGE Profiles
To understand how elaborate petals of N. damascena were made through development, we conducted a detailed time-course transcriptomic analysis. Before that, we first obtained a reference transcriptome using mixed RNAs of floral buds, bracts, and leaves at different developmental stages, in which a total of 38,310 protein-coding genes were predicted (Supplemental Figure 1). Second, we generated digital gene expression (DGE) profiles for petals at eight developmental stages: S4, S5, S6, S7, S8/9, S10, S11, and S12 ( Figure 2A). For comparison, sepals, stamens, and carpels at four developmental stages (i.e., those in the flowers bearing S4, S6, S8/9, and S12 petals, respectively), as well as developing bracts and leaves (i.e., those that are just underneath the floral buds with S6 petals), were also included in the DGE profiles ( Figure 2A). Because the reproducibility among the biological triplicates was very high (coefficient of determination [R 2 ] > 0.85; Supplemental Figure 2), we used the average reads per kilobase of transcript per million mapped reads (RPKM) values of the triplicates as the expression level of the gene in the corresponding sample. As a result, a total of 30,790 genes were defined as expressed (RPKM $1.0) in at least one sample, of which 26,169 were expressed in petals (Supplemental Data Set 1).
To understand the general pattern of our data, we performed pairwise Pearson and Spearman correlation coefficient (PCC and SCC, respectively) analyses of the aforementioned 22 samples. We found that, for each type of organ, adjacent stages usually had higher correlation values than nonadjacent stages, suggestive of progressive changes at the transcriptomic level ( Figure 2B). The PCC/SCC correlation value between petals at S4 and those at all other stages (i.e., S5 to S12), for example, was 0.97/0. 95, 0.73/ 0.90, 0.61/0.87, 0.19/0.77, 0.04/0.52, 0.04/0.43, and 0.03/0.31, respectively. In addition, we found that the correlation values between different types of floral organs were sometimes larger than those between the same types of organs at different stages ( Figure 2B). The PCC/SCC correlation values between S4 petals and S4 stamens, S5 petals and S4 stamens, and S4 stamens and S4 carpels, for instance, were all higher than 0.90, suggestive of very high similarity in the genes expressed. Similar phenomena were observed in a principal component analysis (PCA), although the relationships among the different stages of the same types of organs became more evident ( Figure 2C). Taken together, these results suggest that at the early stages of flower development, different types of floral organs are actually not very different from each other in terms of the genes expressed; yet, during development, when different sets of genes were switched on and off, the differences become greater and greater.

Candidate Genes Associated with Morphological Characters
To identify the genes responsible or indispensable for the formation of key morphological characters on N. damascena petals, we conducted more detailed investigations of the DGE profiles. We started with six essential and well-characterized developmental programs in lateral organ development: organ initiation, organ size determination, organ polarity specification, nectary development, trichome differentiation, and anthocyanin biosynthesis (Hay and Tsiantis, 2010;Yang and Ye, 2013;Albert et al., 2014;Fukushima and Hasebe, 2014;Huang and Irish, 2016;Min et al., 2019). We compared the expression patterns of the component genes in petals with those in sepals and stamens (Figure 4; Supplemental Figure 5). In all these programs, there were genes whose expression patterns in petals were largely consistent with those in sepals and stamens ( Figure 4B), suggestive of functional conservativeness. A small number (23 of 66) of genes, however, showed obviously different expression patterns in petals, having either distinct peaks or much higher levels of expression in petals than in sepals or stamens (Figures 4B and 4C; Supplemental Figure 5), suggestive of possible functional diversification. Based on its expression pattern, and because its counterparts in Arabidopsis, PLETHORA1 (PLT1) and AINTEGUMENTA-LIKE6 (AIL6), function to promote cell proliferation in floral organs (Krizek, 2009;Wuest et al., 2012), we hypothesize that NidaPLT1-1 (peaked at S5) is required for the formation of the nectar chamber, where an increased rate of cell division was observed at S5. Similarly, because ectopic expression of polarity genes can sometimes lead to regional outgrowths (Fukushima and Hasebe, 2014), we propose that the polarity genes NidaKAN1 (peaked at S7) and NidaYAB5 (peaked at S8/9) are required for the formation of pseudonectaries. Consistent with this notion, a recent study has demonstrated that knockdown of NidaYAB5 can indeed led to the complete loss of pseudonectaries (Liao et al., 2020).
By using the same strategy and, more importantly, by correlating the expression patterns of genes with key morphological events in petal development, we identified many other candidate genes ( Figure 5; Supplemental Figure 6). The bifurcation of the lower lip, for example, is very likely controlled by NidaCUC1 (peaked at S7), a NAM, ATAF1/2, and CUC2 (NAC) family member whose ortholog has been reported to control leaf dissection in   S6 S7 S8/9 S10 S11 S12 S4 S6

NidaLMI1 as an Essential Regulator for Petal Elaborations
To test the reliability of our predictions, we conducted expression and functional experiments for one of the predicted candidate genes, NidaLMI1 (Supplemental Figure 7). It has been reported that, as a class I HD-Zip family transcription factor gene, LMI1 is involved in the formation of dissected leaves in Cardamine hirsuta and Gossypium hirsutum, tendrils in Pisum sativum, and trichomes in Cucumis sativus (Hofer et al., 2009;Vlad et al., 2014;Zhao et al., 2015;Andres et al., 2017). In N. damascena, NidaLMI1 was preferentially expressed in petals, with the highest expression being detected in S8/9 petals ( Figure 5), suggesting that its expression has been expanded into petals. To understand the precise spatiotemporal expression of NidaLMI1, we conducted detailed in situ hybridization assays. The expression of NidaLMI1 was undetectable in shoot apical meristem, leaf primordia, or floral primordia before S6 of petal development (Supplemental Figure 8). At S6, the expression signals of NidaLMI1 were first detectable in the newly formed junction of the forked lamina (Supplemental Figure 8) and then spread to the area just above the developing nectar chamber, thereby bridging the two pseudonectaries at S7 ( Figure 6A). Subsequently, the expression signals of NidaLMI1 were restricted to the regions where short trichomes would arise (at S8/9) and the developing short trichomes (at S11 and S12; Figures 6B and 6C). Clearly, the expression of NidaLMI1 not only coincides with the formation of short trichomes but also correlates with the bifurcation of the lower lip.
To understand the function of NidaLMI1, we knocked down its expression using the VIGS technique with two Tobacco rattle virus (TRV) constructs, NidaLMI1_V1 and NidaLMI1_V2 (Supplemental Figure 9). Compared with the TRV2-treated plants (i.e., mock), the TRV2-NidaLMI1_V1-and TRV2-NidaLMI1_V2-treated plants with strong phenotypic changes were very similar to each other. Specifically, the two lobes of the lower lip were completely fused and the short trichomes disappeared from the chewing surface of the lower and upper lips, while leaves and other floral organs remained largely unaffected (Figures 6D to 6K; Supplemental Figure 9). In plants with moderate phenotypic changes, the two petal lobes were partially fused and the short trichomes were sparsely distributed (Supplemental Figure 9). In accordance with their morphological features, very low expression of NidaLMI1 was detected in flowers and petals with strong and moderate phenotypic changes (Supplemental Figure 9). These results not only confirmed the predication that NidaLMI1 is required for the differentiation of short trichomes but also suggest that it is involved in the bifurcation of the lower lip.

DISCUSSION
In this study, by conducting detailed transcriptomic and functional analyses on N. damascena, we explored the mechanisms underlying elaborate petal development and specialized character formation. We found that, in addition to the genes and programs indispensable for simple petal development (such as initiation, identity determination, growth, polarity specification, and maturation), there are genes and programs that are not normally expressed in simple petals (such as the formation of lobes, nectaries, and trichomes), suggestive of uniqueness of the N. damascena petals. By analyzing functionally important developmental programs and correlating the expression patterns of genes with morphological changes of petals, we also identified 30 candidate genes that likely play key roles in the development of elaborate petals (especially marginal/ventral elaboration) and the occurrence of morphological events. Our results not only provide a portrait of elaborate petal development but also elucidate the mechanisms underlying lateral organ diversification in plants.

Mechanisms Underlying Marginal and Ventral Elaboration
It has been proposed that elaborate petals can arise through marginal elaboration or ventral elaboration, or both, with the former being the formation of lobes, teeth, fringes, or appendages along the margin of the petal. In leaves, the formation of lobes, teeth, and leaflets has been generally believed to result from regional expression of cell division inhibitors, such as CUC-like genes of the NAC family, LMI1-like genes of the HD-Zip family, and CINCINNATA (CIN)-like genes of the TCP family (Crawford et al., 2004;Blein et al., 2008;Hofer et al., 2009;Koyama et al., 2011;Sicard et al., 2014;Vlad et al., 2014;Andres et al., 2017). In N. damascena, the bifurcate lobes of the lower petal lip are produced through marginal elaboration at S6. Coincidently, at S6, the expression of NidaLMI1 is restricted to the junction region of the forked lower lip. The fact that knockdown of NidaLMI1 led to the generation of intact rather than bifurcate lower petal lips further confirmed that NidaLMI1 is responsible for the bifurcation of the lower lip. In addition, we found that NidaCUC1 shows a relatively higher expression level in S4 to S7 petals relative to sepals and stamens (Supplemental Figure 5), suggestive of a possible role in marginal elaboration. Taken together, these results suggest that, similar to the formation of lobes, teeth, and leaflets in leaves, marginal elaboration of the N. damascena petals also requires the regional expression of cell division inhibitors.
Unlike marginal elaboration, ventral elaboration refers to the formation of lips, scales, spurs, or other protruded structures on the ventral (i.e., abaxial) side of the petal (Endress and Matthews, 2006;Yao et al., 2019). The mechanisms underlying ventral elaboration are still unclear; yet, the data suggest that it may be related to the reactivation of the meristem program or rearrangement of the expression of adaxial/abaxial genes, or both (Walcher-Chevillet and Kramer, 2016). In Antirrhinum majus and Linaria vulgaris (Plantaginaceae), reactivation of the KNOXI genes, a key regulator of meristem development and maintenance, can cause the production of spur-like outgrowths on petals (Golz et al., 2002;Box et al., 2011). The formation of the peltate leaves of Tropaeolum majus (Tropaeolaceae) and traps of the carnivorous plant Utricularia gibba (Lentibulariaceae), however, has likely been caused by the rearrangement of the expression of the adaxial/ abaxial genes (Gleissberg et al., 2005;Whitewoods et al., 2020). In this study, we found that, at S4, when the upper lip starts to emerge through ventral elaboration, genes of the KNOXI-CUC module, such as NidaSTM1/2, NidaBP, and NidaCUC1, show elevated expression in petals compared to sepals and stamens, suggestive of possible reactivation of the program (Figure 4). Presumably, it was the reactivation of the KNOXI and CUC genes that led to the formation of the upper lip. Nevertheless, we still cannot exclude the possibility that expression rearrangement of adaxial/abaxial genes also contributed considerably to this process, because the abaxial gene NidaYAB5 is also expressed in the adaxial side of the upper lip Liao et al., 2020). More studies are needed to elucidate the exact roles of these genes. (A) to (C) Results of in situ hybridization of NidaLMI1 in petals at S7 (A), S8/9 (B), and S11/12 (C). For petals at S7 and S8/9, series transverse (T in blue) and longitudinal (L in orange or gray) sections are presented, with their relative positions being shown with dashed lines in the left cartoons (adaxial view for the lower lip and abaxial view for the upper lip). For petals at S11 and S12, only longitudinal sections are presented. In the cartoons, the black dashed lines indicate the edge of the nectar chamber and the pink regions represent the expression domains of NidaLMI1 in petals. Hybridization results of the sense probe are presented as negative controls. Bars 5 100 mm.

Molecular Bases of Other Morphological Events
In addition to ventral and marginal elaboration, there are other morphological events that play key roles in the making of elaborate petals. Some of these events have significantly changed the developmental trajectory of the petals and eventually led to the formation of new characters (such as nectaries, pseudonectaries, short trichomes, and long hairs). The exact mechanisms underlying the occurrences of these morphological characters are still unclear, except that a recent study has shown that ectopic expression of one of the abaxial genes, NidaYAB5, may be responsible for the formation of pseudonectaries through regional thickening; knockdown of NidaYAB5 led to the complete loss of pseudonectaries (Liao et al., 2020). In this study, we not only confirmed the importance of NidaYAB5 in petal development but also identified the genes tightly associated with the formation of other characters. The nectary, for example, is likely controlled by NidaSTY1, a STY family gene whose ortholog in Aquilegia coerulea is indispensable for nectary development (Min et al., 2019), whereas the formation of the ridges between pseudonectaries and nectar chamber may be regulated by the NidaTCP5 gene ( Figure 5).
In addition to conical cells, petals of N. damascena also possess short trichomes and long hairs, both of which are unicellular structures. Long hairs, which start to emerge at S6, are ;1200 mm in length and sparsely distributed on the lobes and pseudonectaries of the lower lip, whereas short trichomes, which start to initiate at S10, are ;100 mm in length and densely distributed on the chewing surface of the upper and lower lips. Recently, it has been suggested that while long hairs on petals likely existed before the diversification of the genus Nigella, short trichomes evolved during the evolution of the genus . In this study, we proposed that NidaGL3 may be responsible for the differentiation of long hairs whereas NidaLMI1 is indispensable for the emergence of short trichomes. When the expression of NidaLMI1 was knocked down by using the VIGS technique, short trichomes disappeared, whereas long hairs remained unaffected. This differential effect confirms that long hairs and short trichomes are related but clearly different morphological characters that are possibly controlled by different regulators. Interestingly, it appears that, like conical cells, the initiation and development of short trichomes also require the proper functioning of NidaMIXTA, the ortholog of the well-known conical cell inducer MIXTA, for three reasons. First, during development, short trichomes and conical cells were not distinguishable until S10 . Second, like NidaLMI1, the expression of NidaMIXTA was very low at the early stages (before S7) of petal development but started to increase at S8/9, when both short trichomes and conical cells started to initiate ( Figure 5). Third, ectopic expression of MIXTA into tobacco (Nicotiana tabacum) can trigger the formation of excess trichomes on leaves and floral organs (Glover et al., 1998). These results suggest that NidaLMI1 and NidaMIXTA have something in common in determining the differentiation of short trichomes, although the exact relationship between the two genes awaits to be tested.

Implications for the Origination of New Characters
It has been proposed that co-option is one of the most important mechanisms underlying the origination of new characters (Shubin et al., 2009;Arthur, 2011). By definition, co-option is the deployment of a gene or program with a well-established function to a new function over the course of evolution (Arthur, 2011). The origins of zygomorphic flowers in multiple flowering plant lineages and horns in beetles, for example, can be explained by the cooption of the CYCLOIDEA (CYC)-like genes and the limb program, respectively (Shubin et al., 2009;Hileman, 2014). In the genus Nigella, elaborate petals are made through both marginal/ventral elaboration and the successive originations of several new characters, such as the aforementioned nectary, pseudonectaries, long hairs, and short trichomes . In this study, we found that the originations of at least some new characters may have also been the results of co-option. The short trichomes, for instance, have likely evolved through co-option of the LMI1 gene, because the LMI1 orthologs in other species (including Arabidopsis) are usually not expressed in petals. In fact, a recent study has shown that the nectaries likely arose from cooption of the STY-like genes into the nectariferous petals before the divergence of the Ranunculaceae and Berberidaceae (Min et al., 2019). If these hypotheses are correct, it would be interesting to know how these genes and programs were co-opted into petals.
It has been suggested that the origination of one character sometimes can induce, or stimulate, the origination of other characters, as is the case for the feathers and wings of birds (Prum, 1999;Prum and Brush, 2002). If this prediction is correct, it implies that different characters may have connections in terms of the underlying mechanisms. Interestingly, in this study, we found that the NidaLMI1 gene has dual functions, that is, in both the bifurcation of the lower lip and the production of short trichomes; when the expression of NidaLMI1 was knocked down, short trichomes disappeared and the two lobes of the lower lip failed to bifurcate ( Figure 6). In Nigella, the lobes on the lower lip and short trichomes are the new characters originated at different stages of evolution . It is therefore very likely that LMI1 has been co-opted twice during evolution, with the first being related to the bifurcation of the lower lip and the second being associated with the production of short trichomes. While this hypothesis awaits to be tested, there are already some data that support it. For example, it has been shown that the LMI1 orthologs are involved in the formation of dissected leaves, tendrils, and trichomes in Ca. hirsuta/G. hirsutum, P. sativum, and Cu. sativus, respectively (Hofer et al., 2009;Vlad et al., 2014;Zhao et al., 2015;Andres et al., 2017). Hence, it would not be surprising if it were co-opted in Nigella petals to control the formation of lip lobes and short  trichomes. The relationship between the two new characters (i.e., whether the formation of the former is the prerequisite for the latter), however, remains unclear. A phylostratigraphic study on the DGE data and functional dissection of the key candidate genes will probably help elucidate the molecular mechanisms underlying, as well as the relationship between, these two successively originated new characters.

Problems, Challenges, and Future Directions
It should be mentioned that although the candidate genes identified in this study are very likely key regulators of petal development in N. damascena, it is still too early to say that the mechanisms underlying petal elaboration have been elucidated, for three reasons. First, except for NidaLMI1 and NidaYAB5, the expression patterns and functional properties of all other candidate genes have not yet been studied. Even for the genes with functional studies, the results were based on the VIGS rather than mutant, transgenic, or other more reliable analyses (Gonçalves et al., 2013;Wang et al., 2015;Liao et al., 2020;this study). Second, the candidate genes were identified because their expression patterns in petals match the occurrence of key morphological event and because their orthologs in other species have been reported to perform relevant functions ( Figure 5; Supplemental Figure 6). Yet, because members of the same coexpression module tend to have the same or very similar functions, it is quite possible that other genes, especially those whose orthologs in other species have not yet been functionally characterized or do not perform relevant functions, actually play more important roles than the ones we identified. Third, and most importantly, the candidate genes were interesting because they have interesting expression patterns, yet the reasons why they have such interesting expression patterns are still unclear. In the future, while detailed expression and functional investigations could and should be conducted on the identified candidate genes, special attention should also be paid to the other genes of the corresponding coexpression modules and the factors that shape their expression patterns. It is hoped that, by doing this, we can not only identify more regulators of elaborated petal development but also elucidate the molecular bases of key morphological events.
In spite of these caveats, our study clearly demonstrates that detailed time-course transcriptomic analyses, followed by necessary expression (e.g., in situ hybridization) and functional (especially VIGS) investigations, may be a very promising strategy in elucidating the molecular basis of plant evolution. Functional studies on nonmodel systems can be extremely challenging; yet, many plants with interesting and/or important traits (including, but not limited to, the elaborate petals) are nonmodel species for which transgenic or even genetic manipulations are very difficult. A high-quality transcriptome, therefore, can act as an excellent starting point for studies when the reference genome of the species is lacking. Next, by conducting detailed DGE and bioinformatic analyses, genes, pathways, and networks responsible for, or at least tightly related to, the formation of the trait under consideration may be identified, and functional verification (if necessary) can be done by using RNA in situ hybridization, VIGS, and/or other technologies. More importantly, when this strategy is applied to multiple species, it will become possible to understand the commonness and peculiarity of a trait and elucidate the molecular underpinnings of its evolution.

Plant Materials and Growth Conditions
Seeds of Nigella damascena were sown in nutrient soil and grown in a growth chamber under a 12-h-light (4280 Lux)/12-h-dark photoperiod at 24°C and 60% relative humidity.

Generation of the Reference Transcriptome and the DGE Profiles
Total RNAs were separately extracted from floral buds, bracts, and leaves at different developmental stages using the SV Total RNA Isolation System (Promega) and then equally mixed to generate a pool for library construction. Paired-end 100-bp-long reads were sequenced by Illumina HiSeq2000 (Novogene), and the clean reads were assembled into transcripts using Trinity (Grabherr et al., 2011) after the adaptors and the lowquality reads were removed. Next, an all-by-all BLASTN was performed to group transcripts with different alternative splice forms into unigenes. Subsequently, the unigenes were annotated by sequence similarity comparisons against the nonredundant protein database in National Center for Biotechnology Information and SwissProt using BLASTX (E value < 1e-5). If no best hit was found, the unigenes were annotated through ESTScan (Iseli et al., 1999), and only protein sequences longer than 100 amino acids were retained. Finally, the reference transcriptome of N. damascena was generated (Supplemental Figure 1).
For the generation of the DGE profiles, 22 samples, each with three biological replicates, were subjected to total RNA extraction as described above. In total, 66 libraries were constructed independently for single-end 100-bp-long reads sequencing on Illumina HiSeq2000. The clean reads of resultant 66 DGE profiles were separately mapped to the reference transcriptome by TopHat2 (Kim et al., 2013) and then the mapped reads were counted by HTSeq v0.5.4p3 (Anders et al., 2015) to calculate the RPKM values (Mortazavi et al., 2008). The quality of all the 66 DGE profiles was reflected by read utilization rates (82.40 to 87.49%) and gene coverage (74.84 to 90.35%; Supplemental Data Set 8). The reproducibility among the biological triplicates was evaluated by the pairwise coefficient of determination of the RPKM values, which was calculated in R using log 10transformed RPKM values [i.e., log 10 (RPKM 1 1)] as inputs.

PCC, SCC, and PCA Analyses
The overall comparisons of the expression profiles between samples were conducted by PCC, SCC, and PCA analyses. PCC and SCC were performed under the cor.test function in R using the methods of Pearson and Spearman, respectively. PCA was performed under the prcomp function in R using the original RPKM values as inputs.

Identification of Specifically and Preferentially Expressed Genes and GO Enrichment Analysis
A gene was defined as expressed in a specific type of organ if its highest RPKM is $1.0 in the sampled stages, and its expression level was determined in terms of the highest RPKM value. A gene was considered as specifically expressed in a specific organ if its RPKM is $1.0 only in this organ but <1.0 in any other organs. A gene was regarded as preferentially expressed in a specific organ if its RPKM in this organ is at least twofold higher than those in any other organs. The same criterion was used to identify the genes specifically and preferentially expressed in the petals at different stages. For the GO enrichment analysis of specifically and preferentially expressed genes, the protein sequences were BLAST against the Arabidopsis (Arabidopsis thaliana) nonredundant protein database with an E value cutoff <1e-10 (Ashburner et al., 2000). Next, the GO terms of each gene were adopted in terms of its best hit, and the GO enrichment of a cluster of genes was performed using the agriGO program (Tian et al., 2017) with the false discovery rate threshold of 0.05.

Gene Coexpression Analysis
The coexpression modules were identified by R package WGCNA (Zhang and Horvath, 2005;Langfelder and Horvath, 2008) based on the RPKM data of 13,395 genes whose expression showed coefficient of variation >0.5 across the eight stages. We first created a matrix of pairwise PCCs between all pairs of genes and constructed an adjacency matrix by raising the coexpression measure, 0.5 1 0.5 3 PCC, to the power b 5 14, which was interpreted as an optimal soft threshold of the correlation matrix. The resultant adjacency matrix was then used to calculate the topological overlap, based on which genes were hierarchically clustered. The Dynamic Hybrid Tree Cut algorithm was used to cut the hierarchical clustering tree. Clusters with fewer than 20 genes were merged into their closest larger neighbor cluster. Each module was summarized by the first principal component of the scaled gene expression profiles (referred to as module eigengene). Module-stage associations were qualified by PCC analysis, where each module was represented by its module eigengene, and each stage was represented with a numeric vector of "1" for a specific stage and "0" for all the other stages. The same strategy as described above was used for the GO enrichment analysis of the coexpression modules.

Phylogenetic Analysis
Coding sequences of the class I HD-Zip family genes were retrieved through BLAST searches against publicly available databases by using Arabidopsis genes as queries (Supplemental Table 1). Sequence alignment and adjustment were conducted as previously described by Yu et al. (2016). The alignable nucleotide matrix was used for phylogenetic reconstruction in PhyML 3.0 using the maximum-likelihood method (Guindon et al., 2010). The general time reversible (GTR) 1 I 1 G model was applied, and 1000 bootstrap replicates were performed. The matrix in a PHYLIP format and the tree topology in a Newick format can be found in Supplemental Data Sets 9 and 10, respectively.

Gene Expression and Functional Studies
To amplify the cDNA sequence of NidaLMI1, the total RNA was reversetranscribed into cDNA using SuperScript III First-Strand Synthesis System (Life Technologies) and used as the template. Amplified fragments of expected length were purified and cloned into pEASY-T3 cloning vector (TransGen) for sequencing. The expression patterns of NidaLMI1 were revealed by mRNA in situ hybridization, following the procedures as described by Wang et al. (2015).
The function of NidaLMI1 was investigated through the VIGS technique. Two specific fragments, which spanned the 39 and 59 ends of NidaLMI1, respectively, were amplified, purified, and independently introduced into the TRV2-based PY156 vector to generate constructs NidaLMI1_V1 and NidaLMI1_V2 for VIGS experiments (Supplemental Figure 9). Wild-type plants were treated at least two rounds (Supplemental Table 2). Flowers showing phenotypic changes were photodocumented. The morphology and micromorphology of flowers with visible phenotypic changes were investigated as described by Yao et al. (2019). The silencing efficiency of VIGS treatments was determined by quantitative RT-PCR (qRT-PCR) as described by Wang et al. (2015). The cDNAs from individual floral buds or petals at S8/9 of the wild-type and VIGS-treated plants were used as templates. The primers used in this study are provided in Supplemental  Table 3. The significant difference of the relative expression levels between the mock and the lines showing strong phenotypic changes was evaluated using the wilcox.test with the two-sided function in R. To avoid the offtarget effect, the sequence specificities of the in situ probe, VIGS fragments, and qRT-PCR primers were evaluated by BLASTN against the reference transcriptome with an E value cutoff < 1e-10 (Supplemental Table 4).

Accession Numbers
The RNA-sequencing data have been deposited in National Center for Biotechnology Information Short Read Archive with the accession number PRJNA638226. Sequence data from this article can be found in the GenBank data library under the following accession number: NidaLMI1 (MT780182). Supplemental Table 1. Gene information used for the phylogenetic analysis.
Supplemental Table 2. Information on VIGS experiments. Table 3. Primers used in this study.

Supplemental
Supplemental Table 4. Evaluation of sequence specificities of the in situ probe, VIGS fragments and qRT-PCR primers.

ACKNOWLEDGMENTS
We thank Shin-Han Shiu and Wengen Zhang for valuable comments and members of the Kong laboratory for helpful discussions. We also thank editors and three anonymous reviewers for valuable suggestions. This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant XDB27010304) and the Chinese Academy of Sciences (grant ZDBS-LY-SM022).

AUTHOR CONTRIBUTIONS
R.Z., H.S., and H.K. designed the research. R.Z. carried out the sample collection. P.W. performed the de novo assembly of the reference transcriptome of N. damascena. R.Z. and X.F. analyzed the DGE data with the help of G.X. R.Z. and C.Z. conducted mRNA in situ hybridization, VIGS, and qRT-PCR experiments with the help of H.L. X.Y., X.D., and Y.Y. performed the scanning electron microscopy experiment. J.C. made the virtual clay models. R.Z., X.F., H.S., E.M.K., and H.K. wrote the article.