1. Introduction

Orchidaceae is one of the largest angiosperm families, consisting of more than 800 genera and about 25, 000 species (Dressler, 1993; Ramirez et al., 2007). Orchids are well known for their tremendous diversity in plant forms and growth habits, as well as extraordinary diversity in flower morphology and high value of hybrid cultivars (Ehlers et al., 2002; Givnish et al., 2015). Orchid production has therefore become one of the most important world-wide businesses in the floriculture industry. In China, Japan, Korea, and Southeast Asia, orchids in the genus Cymbidium are economically valued for their beautiful, fragrant flowers and elegant leaves. Cymbidium ensifolium, which belongs to the subgenus Jensoa, blossoms many times a year, and is thus highly desirable in China's flower market (Su et al., 2018a; Zhang et al., 2015).

Orchids are characterized by their highly evolved flower organs, which possess three petal-like sepals in the first whorl, with one in the top or dorsal sepal, and two lower lateral sepals. In the second whorl, there are two lateral petals and a specialized bottom petal, called a lip or labellum, which is oval to triangle and often has a kinked to undulated margin. In the middle of the flower there is a pistil/stigma fused together with pollinia to form the so-called column or gynostemium. Mutations that alter flower morphology are common in the orchid family and have greatly diversified orchid flower forms (Xiang et al., 2018; Yang et al., 2017). For example, in the multi-tepal mutants, the gynostemium is replaced by a newly emerged flower, which continues to produce sepals, petals, and gynostemium centripetally (Yang et al., 2017). In the stamenoidtepal mutant, outer tepals become shorter compared with those on the standard flower (Su et al., 2018b), and in the nongynostemium mutant, the gynostemium is absent and replaced by new sepals in the center part (Yang et al., 2017). Mutants with highly specialized floral morphology are not only valuable, but also represent valuable research materials for studying the molecular mechanisms of orchid flower pattern variation. However, the large genomes and long juvenile phases of orchids have greatly limited the utility of functional genomics and gene discovery approaches that might identify genes involved in the regulating floral patterning (Cai et al., 2015; Chao et al., 2017; Yan et al., 2015; Zhang et al., 2017).

In this work, we examined a domesticated variety of C. ensifolium 'CuiYuMuDan', which we have named the "leaf-like mutant". This variety lacks typical characteristics of orchid floral organs such as zygomorphic flowers, the evolved lip, and the column; instead, these mutants produce sepal-to leaf-like structures along the inflorescence. By comparative analysis of the floral transcriptome between the wild type and the mutant, we identified 6234 differentially expressed genes (DEG) involved. GO and KEGG enrichment indicated that the unigenes responding to membranebuilding and anabolism regulation account for the majority. When referred to individual gene functions, we found 20 DEGs involved in flower development, including 13 MADS-box gene homologs, 3 CONSTANS (CO)-like genes, and other 4 unigenes that annotated as positive or negative regulator of floral meristem activity. Among them, floral homeotic MADS-box genes were investigated by phylogenetic analysis and expression validation. This study shed light on the molecular regulation of orchid flower development, particularly for floral organ defects, and provides new insights and practical guidance for improving their floral patterning regulation and molecular breeding.

2. Materials and methods 2.1. Plant materials and growth conditions

C. ensifolium wild-type plants (also known as 'XiaoTaoHong', the most widely known commercial cultivar in China) and the domesticated variety 'CuiYuMuDan', were artificially cultivated and collected from the cultivation base of the Environmental Horticulture Research Institute, Guangdong Academy of Agricultural Sciences, China. All plants were grown and maintained in pots in a greenhouse at day/night temperatures of 26/23 ℃ under a 16-h light/8-h dark photoperiod.

2.2. Library construction and illumina sequencing

The cDNA libraries of wild-type C. ensifolium 'XiaoTaoHong' and the mutant plant 'CuiYuMuDan' were prepared from an equal mixture of RNAs isolated from mature flowers of five individual plants, respectively. Two replications were included in each sample to create four multiplexed cDNA libraries. Briefly, mRNAs were purified from total RNA using the Oligotex mRNA Midi Kit (QIAGEN, Germany) and quantified using a Nano-Drop 2000 spectrophotometer (Thermo Scientific, USA) to generate the cDNA library according to the Illumina manufacturer's instructions as in our previous work (Yang and Zhu, 2015). MRNAs were purified from total RNA and fragmented to approximately 200 bp. Subsequently, the collected mRNAs were subjected to first strand and second strand cDNA synthesis followed by adaptor ligation and enrichment with a low-cycle according to the instructions of the TruSeq® RNA HT Sample Prep Kit (Illumina, USA). The purified library products were evaluated using the Agilent 2200 TapeStation and Qubit®2.0 (Life Technologies, USA) and then diluted to 10 pM for cluster generation in situ on the HiSeq3000 pair-end flow cell followed by sequencing (2 × 100 bp). An average of more than 75 million reads was generated for each sample.

2.3. Analysis of differentially expressed genes

Gene expression levels were calculated by FPKM values: FPKM =[total transcript fragments/mapped fragments (millions)] × transcript length (kb). Significant differences in gene expression between wild type and mutant were determined using edgeR. The false discovery rate (FDR) was applied to identify the threshold of the P-value in multiple tests; with FDR × 0.05 and|log₂ ratio|> 1 (two-fold change) as the threshold for a significant difference in gene expression.

Differentially expressed genes (DEGs) were annotated using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses following Yang and Zhu (2015) and Yang et al. (2017). First, all DEGs were mapped to GO terms (or KEGG pathways) in the databases (http://www.geneontology.org/ or http://www.genome.ad.jp/) and gene numbers were calculated for every term (or pathway) (Conesa et al., 2005; Kanehisa and Goto, 2000). This was followed by a hypergeometric test to find the significantly enriched terms in DEGs compared to the genome background. For this purpose, we used the corrected P-value ≤ 0.05 as a threshold for the significantly enriched GO terms and Q-value ≤ 0.05 as a threshold for KEGG pathways.

2.4. Real-time quantitative RT-PCR (qRT-PCR) analysis

One μg of total RNA was quantified with a NanoDrop 2000 spectrophotometer (Thermo Scientific) and used for reverse-transcription following an oligo(dT) primed cDNA synthesis protocol (Fermentas). The resulting cDNA was subjected to relative quantitative PCR using Bio-Rad CFX-96 RealTime PCR System (Bio-Rad, USA) in a final volume of 20 μl containing 2 μl of cDNA and 10 μl of SYBR premix Ex-taq^TM (Takara, Japan). To compare gene expression levels between accessions, ubiquitin was used as an internal control for normalization. A relative quantitative method (DDCt) was used to evaluate the quantitative variation. For each reported result at least three independent biological samples were subjected to a minimum of three technical replicates. The primers were designed with Primer 5.0 software and are listed in Table S1.

2.5. Multiple sequence alignment and phylogenetic analysis

A total of 108 complete protein sequences were aligned using MUSCULE, including 7 CeMADS genes from C. ensifolium, 26 AsMADS from Apostasia, 28 DcMADS from Dendrobium catenatum, 25 EpMADS from Erycina pusila and 22 PEQUMADS from Phalaenopsis equestris. Based on the alignment, a Maximum likelihood (ML) phylogenetic tree was constructed using MEGAX with 1000 bootstrap value and Jones-Taylor-Thornton method (Kumar et al., 2016). Genbank accession number of MADS-box genes used in phylogenetic analysis are listed in Table S2.

2.6. Scanning electron microscopy

Dissected apices of mature flowers were fixed in a solution of 3% glutaraldehyde and 2% formaldehyde for 24 h. Samples were then dehydrated in acetone, critical-point dried in liquid CO₂, and mounted on stubs and sputter coated with 25 nm gold. Samples were examined using a JSM-6360LV (JEOL) scanning electron microscope.

3. Results 3.1. Morphology of Cymbidium ensifolium leaf-like flower mutant

The wild-type flowers of C. ensifolium 'XiaoTaoHong' have three sepals in the first whorl and three petals in the second whorl. Together these are called the tepals. Two of the lateral petals have the same shape, and the bottom petal which is called the labellum (or lip) is highly evolved, with an ovate to triangular shape (Fig. 1A and B). The male and female reproductive organs are highly fused to form a gynostemium, which evolved through the complete fusion of the style, stigma and staminal filament, and has four pollinia on a semi-circular viscidium (Fig. 1C and D). In comparison, these typical characteristics of orchid flowers are absent from the mutant plant 'CuiYuMuDan' and differentiation of floral organs is severely inhibited. Instead of differentiated floral organs, the leaf-like flower mutant has sepal-to leaf-like structures along the inflorescence (Fig. 1E-G). To identify differences in the shapes of wild type and mutant floral organs, we observed perianth epidermal cells with a scanning electron microscope. In wild-type flowers, the epidermal cell shapes and ultrastructural features of wild-type sepals, petals, lips, and columns vary (Fig. 2A-H). On both the adaxial and abaxial epidermis, sepal cells are polygonal; petal cells are hemispherical; lip cells are spherical; and column cells are papillae/cupola/conica. In contrast, mutant sepal and petal-like structures (Fig. 1, G) have no polygonal to ellipticalshaped cells; instead, consistent with the abnormal smooth surface of the leaf-like flower structures, mutant flowers have flat epidermal cells (Fig. 2, I & L).

3.2. Transcriptome sequencing and differentially expressed genes

Transcriptome sequencing generated 77, 520, 950 total reads for the wild-type 'XiaoTaoHong' (WT_Xth) and 85, 119, 788 for the leaf-like flower mutant 'CuiYuMuDan' (Mut_Cymd), with corresponding clean reads of 76, 749, 454 (99.00%) and 84, 486, 426 (99.26%) (Table 1). Currently, no C. ensifolium reference genome sequence exists; thus, we assembled all clean reads (25, 748, 328, 939 bp) de novo and obtained a dataset of 44, 480 unigenes, which had a mean length of 891 bp (NCBI Genbank accession number SRR9117944, SRR9117943 and Table S1). This data set is smaller than the available ESTs derived from C. ensifolium flowers in our previous report (Yang and Zhu, 2015). We therefore mapped all sequencing reads to 111, 892 unigenes reported in our previous work. The total mapped reads for WT_Xth and Mut_Cymd were 68, 700, 479 (89.51%) and 60, 198, 127 (80.08%), respectively. This suggests that most transcriptionally active genes were captured in our initial transcriptome data set.

To identify the transcriptional responses in the mutant flower, we determined FPKM values for each gene and calculated the number of DEGs (FPKM ≥1). Overall, 6234 transcripts showed significantly altered expression in the mutant, with 3578 up-and 2656 down-regulated (Table S3). These findings suggest that these DEGs may play a role in leaf-like flower development.

Numbers inside the parentheses indicate the percentage of total reads. Total mapped:Number of reads mapped to the reference dataset; Multiple mapped:Number of reads with multiple alignment sites; Uniquely mapped:Number of reads with single alignment site.

3.3. Annotation statistics and functional enrichment

GO enrichment analysis revealed that leaf-like flower mutant DEGs can be categorized into 66 functional groups. Of these, 18 GO terms were related to 'biological process', 28 were related to 'cellular component', and 20 were related to 'molecular functions. Within each of the three main categories of the GO classification scheme, the dominant subcategories were 'regulation of transcription & DNA-templated' (278), 'nucleus' (1286) and 'protein binding' (377), respectively. In addition, we noticed that a relatively high-percentage of genes came from the GO terms of 'cytoplasm', 'plasma membrane', 'integral component of membrane' and 'chloroplast' in the 'cellular component' categories (Fig. 3). After normalization to the total number of assigned genes for each of the 66 GO terms, we found significant enrichment occurred in 'secondary shoot formation', 'respiratory chain', 'pentose-phosphate shunt', and 'oxidoreductase activity' (Fig. 4A).

We also conducted an enrichment analysis using KEGG path-ways, which assigned 4239 DEGs to 134 pathways. Of these, the majority of hits fell into 'plant-pathogen interaction (184)' and 'plant hormone signal transduction' (182), followed by 'starch and sucrose metabolism' (149) and 'carbon metabolism' (146) (Table S4). However, the most significant enrichment was detected in the 'citrate cycle' (Fig. 4, B). These findings indicate that the 'citrate cycle' pathway, which functions as the ultimate metabolic pathway of the three major nutrients (carbohydrates, lipids and amino acids), is involved in the regulation of leaf-like flower development. More details of the pathway annotations for significant hits in unigene sets are provided in Table S4.

3.4. Differential expression of transcription factors

Transcription factors play a critical role in the control of plant reproductive development. In our study, a total of 152 putative transcription factors were detected with significant changes in expression between wild-type and the mutant, with 110 up-and 42 down-regulated, respectively (Table S5). We found that the genes that showed the greatest differences in expression were transcription factors in the bHLH family (24), followed by MYB (15), WRKY (14), and MIKC families (13) suggesting that transcription factors in these gene families play a crucial role in C. ensifolium flower development (Fig. 5). We also noted expression differences in NAC, TCP, SBP and TALE transcription factor families. Expression differences in these developmentally important transcription factor families suggest that genes in these families may play roles in leaf-like flower mutant defects in floral organ development.

3.5. Differentially expressed genes involved in floral development

We specifically searched for the orthologs of genes and gene families that are known to be involved in flower initiation and the ABCDE flower developmental model of Arabidopsis. According to the unigene annotations, we identified 20 genes involved in flower development. Eleven of these gene were up-regulated, whereas nine genes were down regulated. These DEGs associated with floral development included 13 unigenes which were identified as MADS-box gene homologs, 3 CONSTANS (CO)-like genes, and other four unigenes that were annotated as positive or negative regulators of floral meristem activity, including CLAVATA3 (CLV3), WUSCHEL (WUS) and PERIANTHIA (PAN) (Somssich et al., 2016; Wang et al., 2017; Uemura et al., 2018). Expression of two genes homologous to WUS increased by 1.49-and 1.09-fold respectively, whereas two unigenes with CLV3-like and PAN-like sequences, were expressed at lower levels in the mutant (Table 2). According to reports on Arabidopsis, WUS is required to maintain stem cells in an undifferentiated state. The size of the WUS expression domain controls the size of the stem cell population through WUS indirectly activating the expression of CLV3 in the stem cells and CLV3 repressing WUS transcription through the CLV1 receptor kinase signaling pathway (Guo et al., 2018; Uemura et al., 2018). Thus, our findings that, in the mutant, WUS-like sequences are expressed at high levels and CLV3- like genes are expressed at low levels is consistent with sustained differentiation of the inflorescence axis.

3.6. Confirmation of differential expression of floral homeotic genes

MADS-box genes are known to play key roles in flower development; we therefore analyzed MADS-box genes that were differentially expressed in the leaf-like flower mutant. We used 50RACE to sequence complete MADS-box genes (Table S6) and then constructed a maximum likelihood (ML) phylogenetic tree of 108 protein sequences from four orchid species. Seven Cymbidium ensifolium MADS-box genes were clustered with floral homeotic genes isolated from other orchid species:two AP1/FUL homologs, two SEP homologs, and one putative ortholog each of AGL6, AP3 and MIKC^* genes (Fig. 6).These gene sequences could be used as an important resource for future studies on floral development and flower organ formation.

We used real-time quantitative RT-PCR to validate the differential expression of these seven MADS-box genes. A-class AP1 homologous genes TDN27680 and TDN23829 were expressed 12.3- and 41-fold higher, respectively, in the mutant than in the wildtype. Similarly, AGL6-like TDN30954 and MIKC^* gene TDN23236 were expressed 8.5-and 2.5-fold higher, respectively, in the mutant than that in the wild-type flowers (Fig. 7). In contrast, B-class AP3 homolog TDN23421 and E-class SEP-like sequences TDN29274 and TDN28990 were expressed 0.9-to 0.3-fold less in the mutant than in the wild-type. The differential expression of these floral homeotic genes between mutant and wild-type flowers is consistent with our initial analysis of differential gene expression.

4. Discussion

Leaf-like flower variations occur frequently in angiosperms. The genes that are thought to regulate these variations (e.g., LEAFY, AGL24, SVP, phytoplasma effector SAP54) have been well documented in several plant species, including Arabidopsis, rice, and tomato (Ferrario et al., 2004; Fornara et al., 2008; Li et al., 2016; MacLean et al., 2011; Schultz and Haughn, 1991; Yu et al., 2004). Numerous economically valuable Cymbidium varieties have been produced that have diverse inflorescence architectures, including leaf or flower color and floral patterning. Nevertheless, few studies have examined the molecular mechanisms of floral organ development and variation in Cymbidium spp.

During the past decade, genomics approaches have been widely applied to horticultural crops, and several Cymbidium EST libraries have been reported (Sun et al., 2016; Yang and Zhu, 2015; Yang et al., 2017, 2019). Application of next-generation sequencing technologies has greatly advanced the study of orchid genomics. Recently we published a transcriptome data set containing 111, 892 C. ensifolium transcript clusters. In this study, we dissected the C. ensifolium leaf-like flower mutant for the first time. Comparative transcriptome analysis identified 6234 DEGs, most of which are involved in cell division and anabolism. Gene Ontology terms enriched for defective flowers were related to 'mitochondrion', 'plasma membrane', 'nucleus', and 'RNA transport'; whereas, KEGG pathways were related to 'Starch and sucrose metabolism', 'Carbon metabolism', 'Biosynthesis of amino acids', 'Protein processing in endoplasmic reticulum', 'Phenylpropanoid biosynthesis', 'RNA transport', and 'Amino sugar and nucleotide sugar metabolism'. The differential expression of genes involved in plant cell expansion and energy metabolism is consistent with the developmental processes that underlie flower development, i.e., cell division, membranebuilding, and regulation of anabolism.

We also identified transcription factors that correspond to floral organ defects in the mutant. Specifically, we identified 110 transcription factors that are up-regulated in the leaf-like flower mutant and 42 transcription factors that are down-regulated. The majority of these genes belong to the transcription factor family bHLH followed by MYB, WRKY, and MIKC gene families. The differential expression of these genes in mutant flowers suggest that they play crucial roles in C. ensifolium flower development. In addition, expression differences were detected in developmentally important transcription factor families such as NAC, TCP, SBP, and TALE. Taken together, these findings suggest that flower development is regulated by a complex regulatory network of genes. For example, the NAC genes CUC1 and CUC2 have been shown to be involved in petal initiation and petal number determination; furthermore, other NAC members may function in different stages of petal development, such as petal initiation and petal expansion (Cucinotta et al., 2018; Spinelli et al., 2011). TCP proteins have important functions in plant development, such as control of leaf shape and flower symmetry, and the suppression of lateral shoot branching (Nicolas and Cubas, 2016a, 2016b). In addition, TALE genes are required for meristem maintenance and proper patterning of organ initiation, and play a critical role in the diversity of flower and leaf form (Furumizu et al., 2015; Tao et al., 2018). The regulatory relationship between SPL proteins and the meristem identity AP1/FUL-like genes has also been well documented (Chen et al., 2018; Liu et al., 2017).

We further characterized differentially expressed genes from the MIKC type MADS-box transcription factor family, including two AP1/FUL homologs, two SEP homologs, one putative ortholog of AGL6, AP3, and MIKC^*. Quantitative RT-PCR confirmed that in the mutant A-class AP1 homologous genes TDN27680 and TDN23829, AGL6-like TDN30954, and MIKC^* gene TDN23236 were highly expressed, whereas expression of B-class AP3 homolog TDN23421 and E-class SEP1/3-clade genes TDN29274 and TDN28990 decreased.

Previous studies on the regulation of flower organ development have noted that perianth patterning in orchids is mediated by ABCE model-related MADS-box genes. For example, during flower patterning in Phalaenopsis, PeMADS6/PeMADS2-5 (PI-like/AP3- like) form important heterodimeric complexes; in Dendrobium DcOSEP1/DcOPI/DcOAP3A or DcOAP3B (SEP-like/PI-like/AP3-like) form multimeric proteins; similarly, in Oncidium, AGAMOUS-LIKE6- like (AGL6-like) and AP3-like proteins form multimeric proteins; and in Habenaria radiata, HrDEF-C3/HrAGL6-C2 form complexes (Tsai et al., 2004; Chang et al., 2010; Hsu et al., 2015; Mitoma et al., 2019). Notably, SEP genes serve as glue factors to interact with the mass of multi-MADS-box proteins and facilitate the formation of various protein complexes during sepal, petal, and lip determination. Pan et al. (2014) revealed that in Phalaenopsis dimeric, trimeric and tetrameric protein-protein interactions form between PeSEP and other MADS-box proteins, such that participation of these four SEP paralogs in Phalaenopsis varies during floral transition, morphogenesis, and ovule development. For example, PeSEP3 has a significant effect on floral morphology. Virus-induced silencing of PeSEP3 alters the epidermal identity of tepals and the contents of anthocyanin and chlorophyll, causing tepals to become leaf-like organs. In our work, the expression of CeSEP1/3-clade genes TDN29274 and TDN28990, the orthologs of PeSEP1 and PeSEP3, respectively, was reduced by 30-90%. However, the other two CeSEP genes showed equal or slightly higher expression in the mutant (data not shown). This result agrees well with research from Phalaenopsis that shows that defects of CeSEP1/3-clade genes contributes to the leaf-like flower phenotype in the mutant. Furthermore, these results lend support to the idea that SEP paralogs differ in their ability to regulate floral organ specificity.

This is the first study to identify effectors of leaf-like flower development. These transcriptome data sets provide a significant resource for the discovery of genes related to floral patterning and improve our understanding of the floral morphological diversity of the orchid family.

Declaration of Competing Interest

We would like to submit the enclosed manuscript entitled "Transcriptome Analysis Reveals Clues into Leaf-Like Flower Development in Chinese Orchid C. ensifolium", which we wish to be considered for publication in Plant Diversity. No conflict of interest exits in the submission of this manuscript, and the manuscript is approved by all authors for publication. I would like to declare on behalf of my co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part.

In this work, we dissected a domesticated variety of C. ensifolium 'CuiYuMuDan', which continuous to produce sepal-to leaf-like structures and named as leaf-like flower mutant accordingly. Using comparative transcriptome analysis, we identified 6234 differentially expressed genes (DEGs) mainly involved in membranebuilding, anabolism regulation, and plant hormone signal transduction implying their potential roles involved in floral development deficiency in the mutant. When annotated with plant transcription factor database, we found 152 differentially expressed transcription factors and the most differences occurred in the families of bHLH, MYB, MIKC, and WRKY. Moreover, when referred to individual gene functions, we found 20 DEGs involved in flower development, including MADS-box genes, CLAVATA3 (CLV3), WUSCHEL (WUS) and PERIANTHIA (PAN). Among them, 7 floral homeotic genes were further investigated by phylogenetic analysis and expression validation, which displayed distinctive spatial expression patterns and significant changes between the wild type and the mutant. Our result shed light on the molecular regulation of orchid flower development, particularly represents the first report on transcriptome analysis of floral organ deficiency in C. ensifolium, and provides new insights and practical guidance for improving floral patterning regulation and molecular breeding of Chinese orchid plants. So we believe that our findings could be of interest to the readers of Plant Diversity, and we hope that the editorial board will agree on the interest of this study.

We deeply appreciate your consideration of our manuscript, and we look forward to receiving comments from the reviewers. If you have any queries, please don't hesitate to contact me at the address below.

Acknowledgements

This research was funded by grants from National Key R & D Program (2018YFD1000404), the National Natural Science Foundation of China (31672184), the Natural Science Foundation of Guangdong Province (2017A030312004), Guangzhou Science and Technology Project (201707010307, 201904020026), Innovation Team of Modern Agricultural Industry Technology System in Guangdong Province (2019KJ121), and the Guangdong Academy of Agricultural Sciences Discipline Team Construction Project (201612TD, 2017A070702008, 201721).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2019.12.001.