Upregulation of the MYB2 Transcription Factor is Associated with Increased Accumulation of Anthocyanin in the Leaves of Dendrobium bigibbum

Orchids with colorful leaves and flowers have significant ornamental value. Here, we used γ-irradiation-based mutagenesis to produce a Dendrobium bigibbum mutant that developed purple instead of the normal green leaves. RNA sequencing of the mutant plant identified 2513 differentially expressed genes, including 1870 up- and 706 downregulated genes. The purple leaf color of mutant leaves was associated with increased expression of genes that encoded key biosynthetic enzymes in the anthocyanin biosynthetic pathway. In addition, the mutant leaves also showed increased expression of several families of transcription factors including the MYB2 gene. Transient overexpression of D. biggibum MYB2 in Nicotiana benthamiana was associated with increased expression of endogenous anthocyanin biosynthesis genes. Interestingly, transient overexpression of orthologous MYB2 genes from other orchids did not upregulate expression of endogenous anthocyanin biosynthesis genes. Together, these results suggest that the purple coloration of D. biggibum leaves is at least associated with increased expression of the MYB2 gene, and the MYB2 orthologs from orchids likely function differently, regardless of their high level of similarity.


Introduction
In Orchidaceae, Dendrobium species are one of the most popular orchids known for their medicinal and commercial value in potted and cut flower industries [1]. The Dendrobium genus contains approximately 1800 species and are mainly distributed throughout Asia and the South Pacific [2]. Dendrobium catenatum (also named Dendrobium officinale), Dendrobium nobile, and Dendrobium candidum are used in herbal medicines in many Asian countries [3]. Moreover, the Dendrobium genus is known for its valuable floral traits including colors, morphologies, and scent and Dendrobium species are regarded as some of the important commercial cut flowers. A variety of Dendrobium hybrids have been created that have improved flower colors. However, limitation of genetic resources in Dendrobium limits the extent to which flower color can be modified.
D. bigibbum is an epiphytic or lithophytic orchid that contains cylindrical pseudobulbs, each having between three and five green or purplish leaves and arching flowering stems with up to 20 usually lilac-purple flowers. The D. bigibbum plants containing purple spots on their leaves are very popular in the commercial market. Although colorful leaves and flowers add significant ornamental values to orchids, our understanding of the differential pigmentation in D. bigibbum remains limited.

The Purple Mutant of D. Bigibbum Accumulates Higher Levels of Anthocyanin
We used γ-irradiated D. bigibbum rhizomes to produce a mutant that developed purple leaves in comparison to the green leaves seen on wild type (WT) plants ( Figure 1A,B) ( Figure S1). This mutant, designated as RB016-S7, was propagated through four generations of tissue culturing. To determine whether the purple coloration of the mutant's leaves was associated with anthocyanin pigmentation, we used a pH-differential-based method to quantify the anthocyanin content. The anthocyanin content in the purple leaves (11.68 mg/g dry weight) was~7.0-fold higher than in the green leaves (1.66 mg/g dry weight) ( Figure 1C). Thus, the purple coloration of RB016-S7 leaves was likely associated with the increased biosynthesis of anthocyanins.
To understand the biochemical basis of the increased anthocyanin production in the RB16-S7 mutant, we analyzed genome-wide changes in gene expression. Total RNAs from WT and RB016-S7 leaves were used to construct six cDNA libraries that were sequenced using the Illumina HiSeq 2500 platform. After filtering and quality trimming the raw reads, we obtained 47-66 million high quality reads. Using Trinity, the clean reads from the six libraries were assembled into 110,104 transcripts, To understand the biochemical basis of the increased anthocyanin production in the RB16-S7 mutant, we analyzed genome-wide changes in gene expression. Total RNAs from WT and RB016-S7 leaves were used to construct six cDNA libraries that were sequenced using the Illumina HiSeq 2500 platform. After filtering and quality trimming the raw reads, we obtained 47-66 million high quality reads. Using Trinity, the clean reads from the six libraries were assembled into 110,104 transcripts, with an average length of 1116 bp, and these were then assembled into 32,575 unigenes, with an average length of 1048 bp ( Table 1). The sequence length distribution of unigenes showed that 8373 unigenes (25.7%) ranged from 100 to 500 bp, 11,350 unigenes (34.8%) ranged from 501 to 1000 bp, and 6488 unigenes (19.91%) had lengths of more than 1500 bp ( Figure 2). The 30,714 unigenes were matched with the non-redundant (nr) database, and among these, 26,851 unigenes matched sequences from Dendrobium catenatum, followed by Phalaenopsis equestris (2263) and Apostasia shenzhenica (209) ( Figure 3A). Furthermore, this was consistent with the phylogenetic analysis carried out among native Dendrobium spp, Cymbidium spp, P. equestris and A. shenzhenica orchids, which, as expected, showed relatedness among Dendrobium spp ( Figure 3B).   Q30, base call accuracy of 99.9%; N50, the sequence length of the shortest unigene at 50% of the total genome length.

Functional Annotation and Classification
In the Gene Ontology (GO) analysis, 17,498 unigenes (53.71%) were assigned to three GO terms and were categorized into 41 functional groups (FDR < 0.05) ( Table S1). The GO assignments were Mean length of unigenes (bp) 1,048 N50 of unigenes (bp) 1350 Q30, base call accuracy of 99.9%; N50, the sequence length of the shortest unigene at 50% of the total genome length.

Functional Annotation and Classification
In the Gene Ontology (GO) analysis, 17,498 unigenes (53.71%) were assigned to three GO terms and were categorized into 41 functional groups (FDR < 0.05) ( Table S1). The GO assignments were
To predict and classify the gene functions, we queried all the unigenes against the evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) (v4.5) database. This database contains the functional descriptions and classifications of the orthologous proteins, including Clusters of Orthologous Groups and euKaryotic Orthologous Groups. This analysis allowed us to allocate 27,963 unigenes to 25 eggNOG classifications. Among them, the eggNOG category of functional unknown (S, 27.88%) represented the largest group, followed by signal transduction mechanisms (T, 8.58%), posttranslational modification, protein turnover, chaperones (O, 8.13%), transcription (K, 8.09%), and carbohydrate transport and metabolism (G, 5.70%) ( Figure S3).

Analysis of Differentially Expressed Genes (DEGs) Associated with Anthocyanin Biosynthesis
A total of 2513 DEGs (FDR < 0.05) were identified between the WT and RB016-S7 mutant. Compared with WT, 1870 and 706 genes were up-and downregulated in the RB016-S7 mutant, respectively ( Figure 3C; Table S3). The top 20 significant pathways for the up-and downregulated genes were selected for further analysis. The upregulated genes were mainly enriched in ribosome biogenesis, MAPK signaling pathway, plant-pathogen interaction, plant hormone signal transduction, phenylpropanoid biosynthesis, starch and sucrose metabolism, and flavonoid biosynthesis ( Table 2). The 20 significant pathways for the downregulated genes are listed in Table 3. The downregulated genes were mainly enriched in folate biosynthesis, starch and sucrose metabolism, plant hormone signal transduction, and phenylpropanoid biosynthesis.

Quantitative Real-Time PCR (qRT-PCR) Analysis of the Genes Involved in Anthocyanin Biosynthesis
To confirm the RNA-Seq data, we first selected 10 candidate genes associated with anthocyanin biosynthesis and analyzed their expression levels using qRT-PCR. The qRT-PCR analysis confirmed ~1.5-, 2.3-, and 2.4-fold higher levels for PAL1, PAL2, and PAL4, respectively, in the RB016-S7 mutant

Quantitative Real-Time PCR (qRT-PCR) Analysis of the Genes Involved in Anthocyanin Biosynthesis
To confirm the RNA-Seq data, we first selected 10 candidate genes associated with anthocyanin biosynthesis and analyzed their expression levels using qRT-PCR. The qRT-PCR analysis confirmed 1.5-, 2.3-, and 2.4-fold higher levels for PAL1, PAL2, and PAL4, respectively, in the RB016-S7 mutant compared with the WT. In addition, the qRT-PCR analysis showed that C4H, 4CL, and CHS were induced~2-, 3.3-, and 10-fold in the RB016-S7 mutant compared with the WT. Similarly, F3'H, F3H, DFR, and ANS were induced~1.2-, 7.2-, 5.0-, and 3-fold in the RB016-S7 mutant compared with the WT ( Figure 5). The qRT-PCR data were consistent with results obtained from the RNA-Seq data. Thus, the purple pigmentation in RB016-S7 may be associated with the increased expression levels of genes involved in anthocyanin biosynthesis. Asterisks denote a significant difference between respective WT and RB016-S7 mutant leaves samples (t -text, P < 0.0001).
Next, we analyzed the expression levels of regulatory genes associated with anthocyanin biosynthesis. The RNA-Seq dataset showed that DbMYB2, -30, and -44 were highly upregulated in the RB016-S7 mutant compared with the WT, while the expression of DbMYB75 was not significantly different from in WT plants. Notably, the qRT-PCR analysis was only able to confirm a ~13-fold induction in DbMYB2, while the expression of DbMYB30, -44, and -75 remained at WT levels.
Comparisons of expression levels of genes encoding bHLH TFs showed that only DbbHLH1 was expressed at higher levels in RB016-S7 than WT plants. In comparison, DbbHLH96, -114, and -153 showed WT-like expression levels. DbWD40, which showed a 67.97% identity to the Arabidopsis ortholog AtTTG1, had a WT-like expression level [32]. A recent report also suggests roles for WRKY TFs in anthocyanin biosynthesis. RNA-Seq data showed that several WRKY TFs were highly expressed in the RB016-S17 mutant compared with WT. However, the qRT-PCR analysis was only able to confirm ~1.5-3-fold inductions of DbWRKY24, WRKY31, and WRKY40 genes ( Figure 6). Thus, only a select group of TFs were upregulated in the mutant plant, and these, in turn, could play roles in the regulation of genes involved in anthocyanin biosynthesis. Asterisks denote a significant difference between respective WT and RB016-S7 mutant leaves samples (t-text, p < 0.0001).
Next, we analyzed the expression levels of regulatory genes associated with anthocyanin biosynthesis. The RNA-Seq dataset showed that DbMYB2, -30, and -44 were highly upregulated in the RB016-S7 mutant compared with the WT, while the expression of DbMYB75 was not significantly different from in WT plants. Notably, the qRT-PCR analysis was only able to confirm a~13-fold induction in DbMYB2, while the expression of DbMYB30, -44, and -75 remained at WT levels.
Comparisons of expression levels of genes encoding bHLH TFs showed that only DbbHLH1 was expressed at higher levels in RB016-S7 than WT plants. In comparison, DbbHLH96, -114, and -153 showed WT-like expression levels. DbWD40, which showed a 67.97% identity to the Arabidopsis ortholog AtTTG1, had a WT-like expression level [32]. A recent report also suggests roles for WRKY TFs in anthocyanin biosynthesis. RNA-Seq data showed that several WRKY TFs were highly expressed in the RB016-S17 mutant compared with WT. However, the qRT-PCR analysis was only able to confirm 1.5-3-fold inductions of DbWRKY24, WRKY31, and WRKY40 genes ( Figure 6). Thus, only a select group of TFs were upregulated in the mutant plant, and these, in turn, could play roles in the regulation of genes involved in anthocyanin biosynthesis. Results are representative of two independent experiments. Asterisks denote a significant difference between respective WT and RB016-S7 mutant leaves samples (t-test, p < 0.0001).

DbMYB2 Positively Regulates Anthocyanin Biosynthesis
Increased expression of DbMYB2 in the RB016-S17 mutant suggested that MYB2 could positively regulate expression of anthocyanin genes and thereby anthocyanin levels. This is further supported by an earlier study that showed that Dendrobium hybrid MYB2 positively regulated anthocyanin pigmentation in Dendrobium petals. Amino acid alignment of DbMYB2 with DhMYB2 BS No.3 [33] showed ~92% identity. Likewise, amino acid alignment of MYB2 orthologs from D. hybrid, D. candidum, D. nobile, and Cymbidium sinense showed ~80%, ~80%, ~62%, and 63% identity, respectively, with DbMYB2 ( Figure 7A). The amino acid alignment showed that the R2R3 repeat region was highly conserved among various MYB proteins ( Figure 7A). Phylogenetic analysis between these MYB proteins placed DbMYB2, DhMYB2, and DcMYB2 in the same clade ( Figure 7B).
To determine whether increased expression of DbMYB2 positively regulated expression of anthocyanin genes, we expressed MYB2 genes from D. bigibbum, D. candidum, D. nobile, D. hybrid, and C. sinense in Nicotiana benthamiana and evaluated expression of N. benthaminana genes ANS, DFR, and CHS, which are associated with anthocyanin biosynthesis. All the MYB2 genes showed varying levels of increased expression at 36 h post-agroinfiltration ( Figure 8D-H). Interestingly, however, only transient expression of DbMYB2 was associated with increased expression of ANS, DFR, and CHS in N. benthamiana ( Figure 8A-C). These results strongly suggest that DbMYB2 positively regulates expression of genes associated with anthocyanin biosynthesis, and that higher anthocyanin levels in the RB016-S17 mutant are likely due to higher expression levels of DbMYB2. Inability of other MYB2 orthologs to increase expression of ANS, DFR, and CHS suggests that, regardless of their homology, the MYB2 orthologs function differently.

DbMYB2 Positively Regulates Anthocyanin Biosynthesis
Increased expression of DbMYB2 in the RB016-S17 mutant suggested that MYB2 could positively regulate expression of anthocyanin genes and thereby anthocyanin levels. This is further supported by an earlier study that showed that Dendrobium hybrid MYB2 positively regulated anthocyanin pigmentation in Dendrobium petals. Amino acid alignment of DbMYB2 with DhMYB2 BS No.3 [33] showed~92% identity. Likewise, amino acid alignment of MYB2 orthologs from D. hybrid, D. candidum, D. nobile, and Cymbidium sinense showed~80%,~80%,~62%, and 63% identity, respectively, with DbMYB2 ( Figure 7A). The amino acid alignment showed that the R2R3 repeat region was highly conserved among various MYB proteins ( Figure 7A). Phylogenetic analysis between these MYB proteins placed DbMYB2, DhMYB2, and DcMYB2 in the same clade ( Figure 7B).
To determine whether increased expression of DbMYB2 positively regulated expression of anthocyanin genes, we expressed MYB2 genes from D. bigibbum, D. candidum, D. nobile, D. hybrid, and C. sinense in Nicotiana benthamiana and evaluated expression of N. benthaminana genes ANS, DFR, and CHS, which are associated with anthocyanin biosynthesis. All the MYB2 genes showed varying levels of increased expression at 36 h post-agroinfiltration ( Figure 8D-H). Interestingly, however, only transient expression of DbMYB2 was associated with increased expression of ANS, DFR, and CHS in N. benthamiana ( Figure 8A-C). These results strongly suggest that DbMYB2 positively regulates expression of genes associated with anthocyanin biosynthesis, and that higher anthocyanin levels in the RB016-S17 mutant are likely due to higher expression levels of DbMYB2. Inability of other MYB2 orthologs to increase expression of ANS, DFR, and CHS suggests that, regardless of their homology, the MYB2 orthologs function differently.

Discussion
Anthocyanins are pigments that confer color to various plant parts [34]. The color is determined by the composition and concentration of pigments, which vary greatly among plant species [35]. Cyanidin-3-glucoside is a major anthocyanin found in most plants [36]. Other common anthocyanin pigments present in plants include delphinidin, pelargonidin, peonidin, malvidin, and petunidin. Earlier studies on Dendrobium orchids primarily focused on anthocyanin profiles in flowers and stems, which contain pelargonidin, cyanidin, peonidin, delphinidin, and/or malvidin [24]. In contrast, we were only able to detect malvidin in the leaves of D. bigibbum (data not shown), and its levels were associated with increased purple pigmentation in the RB016-S7 mutant's leaves. Thus, anthocyanin pigments present in leaves versus flowers and stems might be associated with the specific genes expressed in these tissues. We determined that the increased anthocyanin accumulation in RB016-S7 was associated with increased expression levels of PAL, CHS, F3 H, and DFR genes that are involved in anthocyanin biosynthesis. Although the anthocyanin biosynthetic genes are well-conserved, the timing, level, and spatial distribution of anthocyanin biosynthesis are primarily determined by TFs.
A recent study offered useful insights into the functions of WRKY TFs in anthocyanin biosynthesis [37], which in turn regulates the MYB/bHLH/WD40complex [27,[38][39][40]. Our analysis also identified 33 WRKY-, 20 MYB-, 23 bHLH-and 1WD40-encoding genes that were differentially expressed in the RB016-S7 mutant leaves. It is possible that these TFs regulate the expression of one or more genes involved in the anthocyanin pathway (Table S5). An example of complex regulation underlying anthocyanin biosynthesis includes the feed-forward loop mechanism in which TFs regulate each other and jointly regulate target genes [28]. In Dendrobium hybrid petals, DhMYB2 and DhbHLH1 TFs play regulatory roles in the anthocyanin biosynthetic pathway [33]. Consistent with this finding, we determined that the expression levels of DbMYB2 and DbHLH1 were significantly higher in the RB016-S7 mutant. Notably, the A. thaliana and M. truncatula MYB2 proteins act as transcriptional repressors of anthocyanin biosynthesis, and the overexpression of either MYB2 abolishes anthocyanin biosynthesis [28,29]. Likewise, heterologous overexpression of Malus domestica MYB3 (MdMYB3) in Nicotiana tabacum is associated with increased anthocyanin biosynthesis [41]. Thus, orthologs of MYB and possibly other genes involved in anthocyanin biosynthesis may play opposite roles in different plants. This was further evident in our analysis, which showed that increased expression of MYB2 in RB016-S7 plants positively correlated with anthocyanin biosynthesis. This was further consistent with our result that heterologous overexpression of D. bigibbum MYB2 in N. benthamiana led to increased expression of ANS, DFR, and CHS genes. Interestingly, transient overexpression of MYB2 orthologs from other Dendrobium spp. or C. sinense did not alter expression of ANS, DFR, and CHS genes. These results strongly suggest that DbMYB2 positively regulates expression of genes associated with anthocyanin biosynthesis, even though DbMYB2 showed high levels of homology with other MYB2 orthologs. Thus, subtle changes in MYB2 sequence are likely sufficient to alter their function. It is possible that MYB2 in other orchids could regulate anthocyanin biosynthesis genes by serving as a part of the bigger complex that contains other factors like bHLHs or WD40. Deciphering the exact biochemical functions of various TFs involved in anthocyanin biosynthesis in D. bigibbum will require more detailed analyses of these proteins.

γ-Irradiation of in Vitro Shoot Cultures
Six-month-old in vitro regenerated shoots of approximately 3 cm in length were exposed to γ-radiation using 60Co γ-irradiator (60 Gy/24 h) at the Korea Atomic Energy Research Institute, Jeongeup, Korea [42]. The first vegetative generation in which treatment was performed was referred to as M1V1. The study continued until the fourth generation (M1V4) to confirm the stability of the induced traits. The purple-colored leaf mutant RB016-S7 was obtained and its physiological traits were analyzed.

RNA Extraction and Quantitative Real-Time PCR (qRT-PCR) Analysis
Total RNA was extracted from the WT and the RB016-S7 plants using an RNeasy plant mini kit (Qiagen, Hilden, Germany), following the manufacturer's instructions. The RNA concentration and quality of each sample were determined using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose electrophoresis, respectively. The cDNA was transcribed from 500 ng of total RNA using a ReverTra Ace-α-kit (Toyobo Co. Ltd, Osaka, Japan). The qRT-PCR was performed with a CFX96 touch real-time PCR detection system (Bio-Rad, Hercules, CA, USA) using iQ™ SYBR ® Green supermix (Bio-Rad, Hercules, CA, USA). The D. bigibbum actin gene was used as an internal control, and the 2 −∆∆Ct method was used to analyze differential expression levels.
Cycle threshold values were calculated using CFX Manager 3.1 software (Bio-Rad, Hercules, CA, USA). Gene-specific primers are listed in Table S5.

Measurement of Total Anthocyanin Content
After harvesting the leaves, the samples were freeze-dried and subjected to solvent extraction using a solution of 85% ethanol acidified with 15% 1.5 N HCl. The samples were incubated at 4 • C for 24 h. Samples were diluted in two buffer solutions: potassium chloride buffer 0.025 M (pH 1.0) and sodium acetate buffer 0.4 M (pH 4.5). Absorbance was measured via spectrophotometer at 510 and 700 nm after 15 min of incubation at room temperature, respectively. Absorbance was calculated as Anthocyanin pigment (cyanidin − 3 − glucoside equivalents, mg/L) = A×MW×DF×10 3

RNA-Seq Analysis, De Novo Assembly, and Unigene Generation
The cDNA libraries were prepared independently from both WT and RB016-S7 leaves. Low quality and duplicated reads, as well as adapter sequences, were removed from RNA-seq raw data using Trimmomatic with default parameters [44]. The de novo assembly was performed using Trinity (ver. 2.8.4) with default parameters [45]. Afterwards, redundant sequences were removed from the assembled transcript sequences using cd-hit-est (ver. 4.7) with a similarity threshold of 90% (i.e., removing similar sequences sharing more than a 90% identity), generating nr transcript sequences. Protein coding sequences (CDSs) were predicted and extracted from the nr transcript sequences using TransDecoder (ver. 5.5) with a parameter of selection of the longest CDS by comparison with Pfam database. The collection of extracted CDSs was designated as the unigene set and used for further analyses. The completeness of the unigene set was validated by analysis with Benchmarking Universal Single-Copy Orthologs (ver. 3.1.0) [46].

Functional Annotation of Unigenes
Sequences homologous to unigenes were identified using BLASTP analyses (cutoff e-value 1e-5) against the NCBI nr protein database. The GO terms, and eggnog (ver. 3.0) and KEGG pathways, were assigned to the unigenes based on BLASTP results using the Blast2GO program (ver. 5.2.5). Conserved domains in the unigene sequences were identified using InterProScan program (ver. 5.34-73.0) with default parameters. In addition, a KEGG pathway analysis was also performed with the KEGG Automatic Annotation KAAS Server using the single-directional best hit method and searching against representative gene sets from both eukaryotes and monocots.

Expression Profiling of Unigenes
Trimmed high-quality RNA-seq reads were mapped on the unigene sequences using BWA (ver. 0.7.17-r1188) [47] and then, RNA reads mapped on unigene sequences were counted using SAMtools (ver. 1.9) [48]. Fragments per kilobase of transcript per million mapped reads values were calculated using the number of RNA-seq reads mapped on unigene sequences and used for the expression profiling of unigenes.

Identification of DEGs between the WT and RB016-S7 Mutant
The bioconductor package DESeq (ver. 1.22.1) was used to identify DEGs between samples [49].
Genes showing over two-fold expression changes with p-values of less than 0.05 were considered DEGs. The GO enrichment analysis was performed for the DEGs using Fisher's exact test with an adjusted p-value of 0.05 in the Blast2GO program (ver. 5.2.5).

Analysis of Unigenes Involved in Anthocyanin Biosynthesis
For anthocyanin biosynthesis, unigenes assigned to the anthocyanin biosynthetic reference pathway from the KEGG pathway analysis were first selected. In addition, BLASTP searches (cutoff e-value: 1e-10) were performed using 40 Arabidopsis genes involved in anthocyanin biosynthesis [50] as queries, and then, unigenes with high similarity levels (≥ 60% identity and ≥ 80% alignment length) to the query sequences were selected as candidate unigenes that could be involved in anthocyanin biosynthesis. Expression values (fragments per kilobase of transcript per million mapped reads) for the genes were retrieved from expression profiles of the unigenes set and used for generating heatmaps using the R -package pheatmap (ver. 1.0.12).

Cloning of Orchid MYB Genes
The full-length MYB2 cDNA (denphalae23719) was PCR-amplified from D. bigibbum leaves and cloned into the Gateway binary vector pMDC32 vector, under the 35S CaMV promoter. The primers used for amplification of MYB2 sequences are listed in Table S5. All the amplified products were sequenced.

Agroinfiltration of N. Benthamiana
The pMDC32-MYB2 plasmids were transformed into A. tumefaciens strains LBA4404 via the freeze-thaw method [51]. Agrobacteria were grown in the LB medium supplemented with 100 mg/L kanamycin and incubated at 28 • C with shaking. Bacteria were pelleted by centrifugation (14,000g for 5 min) and resuspended to an OD 600 = 0.8 in a buffer containing 10 mM MES pH 5.6, 10 mM MgCl 2 , and 200 µM acetosyringone. Cultures were then incubated for 2-4 h at room temperature. Bacteria were infiltrated into the underside of N. benthamiana leaves using a needleless 1 ml syringe. The agroinfiltrated plants were kept in the growth chamber maintained at 23 • C with a 16-h photoperiod.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/16/ 5653/s1. Figure S1. Images of D. bigibbum. Morphological phenotypes of typical wild type (WT) and the RB016-S7 mutant. Figure S2. Distribution of annotated sequences based on a Gene Ontology (GO) analysis. The GO functional classification assigned 17,498 unigenes to 41 subcategories under the three main GO categories of biological process, cellular component, and molecular function. The x-axis indicates the subcategories and the y-axis indicates the number of unigenes in each category. Figure S3. Numeric distribution of eggNOG annotations of unigenes. Letters on the x-axis refer to the categories on the right. The y-axis indicates the number of unigenes in the corresponding eggNOG category. Figure S4. Distribution of annotated sequences as assessed using a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The x-axis indicates enriched KEGG pathways, and the y-axis represents the number of unigenes within each KEGG pathway. (A) Metabolism; (B) genetic information processing; (C) environmental information processing; (D) cellular processes; (E) organismal systems. Table S1. Significantly enriched GO terms. Table S2: Significantly enriched KEGG pathways. Table S3: Clusters of annotated GO terms in the biological process category enriched in up-and downregulated genes between WT and the RB016-S7 leaves. Table S4: Expression profiles of the regulatory genes for leaf color in D. bigibbum. Table S5: Gene-specific primers used for quantitative real-time PCR and Agrobacterium transient assay. Table S6: MYB2 gene sequence from Dendrobium orchids and Cymbidium sinense for Agrobacterium transient assay.

Conflicts of Interest:
The authors declare no conflict of interest.