*2.3. Analysis of Di*ff*erentially 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.


**Table 2.** Top 20 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of upregulated differentially expressed genes (DEGs).


**Table 3.** Top 20 enriched KEGG pathways of downregulated DEGs.

#### *2.4. Analysis of Anthocyanin Biosynthetic Genes in Identified DEGs*

The mutant showed an increased accumulation of anthocyanin; therefore, we used a KEGG functional enrichment to search for genes associated with anthocyanin biosynthesis among the 2513 DEGs. A total of 17 DEGs, encoding eight key enzymes, were identified, and they were three *PAL* genes (*PAL1*: denphalae05809, *PAL2*: denphalae05806, and *PAL4*: denphalae05808), two cinnamic acid 4-hydroxylase genes (*C4H*: denphalae10925 and denphalae10926), four 4-coumarate CoA-ligase genes (*4CL*: denphalae18583, denphalae22607, denphalae27156, and denphalae27157), four *CHS* genes (denphalae02657, denphalae02658, denphalae05188, and denphalae11910), and one gene each of *F3H* (denphalae02991), flavonoid 3<sup>0</sup> -monooxygenase (*F3*0*H*: denphalae11915), *DFR* (denphalae11241), and *ANS* (denphalae18276). All these DEGs were significantly upregulated in the RB016-S7 mutant compared with WT (Table 4). Among other notable genes that were upregulated in RB016-S7 were TFs that belonged to WRKY (33 genes), MYB (20 genes), bHLH (23 genes), and WD40 (1 gene) groups. Among these, *DbMYB2*, *-4*, *-30*, and *-44*, as well as *DbbHLH1*, *-62*, *-96*, *-114*, and *-148*, were highly expressed in the RB016-S7 mutant (Table S4). The expression patterns of the anthocyanin biosynthetic genes were consistent with the increased anthocyanin levels in the RB16-S7 mutant (Figure 4A,B).

*Biosynthesis* 

genes involved in anthocyanin biosynthesis.


**Table 4.** Expression profiles of anthocyanin biosynthetic genes.

FPKM, fragments per kilobase of transcript per million mapped reads. *ANS* denphalae18276 1083 70.86 180.97 2.22 1.15 FPKM, fragments per kilobase of transcript per million mapped reads.

**Figure 4.** Flavonoid**–**anthocyanin biosynthetic genes in *D. bigibbum*. (**A**) The differentially expressed genes (DEGs) between the WT and RB016-S7 mutant found in leaves are highlighted in blue. Phenylalanine ammonia-lyase, PAL; cinnamic acid 4-hydroxylase, C4H; 4-coumarate CoA-ligase, 4CL; chalcone synthase, CHS; chalcone isomerase, CHI; flavanone 3-hydroxylase, F3H; flavonoid 3′ monooxygenase, F3′H; dihydroflavonol 4-reductase, DFR; and anthocyanidin synthase, ANS. (**B**) Expression profiles determined using fragments per kilobase of transcript per million mapped reads (FPKM) values obtained from RNA-Seq data. Expression values (as FRKM) were not scaled per row to allow the visualization of original FPKM values among samples. The heatmap was generated using the R package pheatmap. **Figure 4.** Flavonoid–anthocyanin biosynthetic genes in *D. bigibbum*. (**A**) The differentially expressed genes (DEGs) between the WT and RB016-S7 mutant found in leaves are highlighted in blue. Phenylalanine ammonia-lyase, PAL; cinnamic acid 4-hydroxylase, C4H; 4-coumarate CoA-ligase, 4CL; chalcone synthase, CHS; chalcone isomerase, CHI; flavanone 3-hydroxylase, F3H; flavonoid 3 0 -monooxygenase, F30H; dihydroflavonol 4-reductase, DFR; and anthocyanidin synthase, ANS. (**B**) Expression profiles determined using fragments per kilobase of transcript per million mapped reads (FPKM) values obtained from RNA-Seq data. Expression values (as FRKM) were not scaled per row to allow the visualization of original FPKM values among samples. The heatmap was generated using the R package pheatmap.

~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

*2.5. Quantitative Real-Time PCR (qRT-PCR) Analysis of the Genes Involved in Anthocyanin* 

#### *2.5. 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. *Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 8 of 17

**Figure 5.** qRT-PCR analysis of 10 genes showing altered expression levels in the RNA sequencing (RNA-Seq) analysis. The genes were associated with anthocyanin biosynthesis. More specifically, A**–** J indicate the relative expression levels of *PAL1*, *PAL2*, *PAL3*, *C4H*, *4CL*, *CHS*, *F3H*, *F3′H*, *DFR*, and *ANS*, respectively. The elongation factor 1-alpha (*EF1a*) gene served as an internal control. Error bars indicate standard deviations (*n* = 3). The experiment was repeated three times with similar results. Asterisks denote a significant difference between respective WT and RB016-S7 mutant leaves samples (*t* -text, *P* < 0.0001). **Figure 5.** qRT-PCR analysis of 10 genes showing altered expression levels in the RNA sequencing (RNA-Seq) analysis. The genes were associated with anthocyanin biosynthesis. More specifically, (**A**–**J**) indicate the relative expression levels of *PAL1*, *PAL2*, *PAL3*, *C4H*, *4CL*, *CHS*, *F3H*, *F3*0*H*, *DFR*, and *ANS*, respectively. The elongation factor 1-alpha (*EF1a*) gene served as an internal control. Error bars indicate standard deviations (*n* = 3). The experiment was repeated three times with similar results. 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. 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. 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.

*Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 9 of 17

**Figure 6.** qRT-PCR analysis of 14 genes showing altered expression levels in the RNA-Seq analysis. The relative expression levels of transcription factor genes in the leaves. More specifically, **A–N** indicate the relative expression levels of *MYB1*, *MYB2*, *MYB3*, *MYB30*, *MYB44*, *MYB75*, *bHLH1*, *bHLH96*, *bHLH114*, *bHLH153*, *WRKY24*, *WRKY31*, *WRKY40*, and *WD40*. The *EF1a* gene served as an internal control. Error bars indicate standard deviations (*n* = 3). The error bars indicate SD (*n* = 3). 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). **Figure 6.** qRT-PCR analysis of 14 genes showing altered expression levels in the RNA-Seq analysis. The relative expression levels of transcription factor genes in the leaves. More specifically, (**A**–**N**) indicate the relative expression levels of *MYB1*, *MYB2*, *MYB3*, *MYB30*, *MYB44*, *MYB75*, *bHLH1*, *bHLH96*, *bHLH114*, *bHLH153*, *WRKY24*, *WRKY31*, *WRKY40*, and *WD40*. The *EF1a* gene served as an internal control. Error bars indicate standard deviations (*n* = 3). The error bars indicate SD (*n* = 3). 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).

#### *2.6. DbMYB2 Positively Regulates Anthocyanin Biosynthesis 2.6. 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 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).

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. 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.

*Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 10 of 17

**Figure 7.** Sequence alignments and phylogenetic tree of DbMYB2 from orchids with known R2R3- MYBs domain. (**A**) Alignments of the full-length deduced amino acid sequences of DbMYB2 with **Figure 7.** Sequence alignments and phylogenetic tree of DbMYB2 from orchids with known R2R3-MYBs domain. (**A**) Alignments of the full-length deduced amino acid sequences of DbMYB2 with other R2R3-MYBs present in other orchids. (**B**) Phylogenetic relationship of DbMYB2 with known anthocyanin MYB regulators from other orchid species.

other R2R3-MYBs present in other orchids. (**B**) Phylogenetic relationship of DbMYB2 with known

**Figure 8.** Transient expression of DbMYB2s were associated with increased expression of ANS, DFR, and CHS in *N. benthamiana*. (**A–C**) Real-time quantitative RT-PCR showing relative expression levels of *NbDFR*, *NbANS*, and *NbCHS* genes in *N. benthamiana* transiently overexpressing *MYB2* genes from different orchids. (**D–H**) Real-time quantitative RT-PCR showing relative expression levels of *MYB2* genes in *N. benthamiana* plants transiently overexpressing *MYB2* orthologs from different orchids. The error bars indicate SD (*n* = 4). Results are representative of three independent

anthocyanin MYB regulators from other orchid species.

**Figure 8.** Transient expression of DbMYB2s were associated with increased expression of ANS, DFR, and CHS in *N. benthamiana*. (**A–C**) Real-time quantitative RT-PCR showing relative expression levels of *NbDFR*, *NbANS*, and *NbCHS* genes in *N. benthamiana* transiently overexpressing *MYB2* genes **Figure 8.** Transient expression of DbMYB2s were associated with increased expression of ANS, DFR, and CHS in *N. benthamiana*. (**A**–**C**) Real-time quantitative RT-PCR showing relative expression levels of *NbDFR*, *NbANS*, and *NbCHS* genes in *N. benthamiana* transiently overexpressing *MYB2* genes from different orchids. (**D**–**H**) Real-time quantitative RT-PCR showing relative expression levels of *MYB2* genes in *N. benthamiana* plants transiently overexpressing *MYB2* orthologs from different orchids. The error bars indicate SD (*n* = 4). Results are representative of three independent experiments. Asterisks denote a significant difference with mock and MYB2 overexpressed plants (*t*-test, *p* < 0.0001).

from different orchids. (**D–H**) Real-time quantitative RT-PCR showing relative expression levels of *MYB2* genes in *N. benthamiana* plants transiently overexpressing *MYB2* orthologs from different orchids. The error bars indicate SD (*n* = 4). Results are representative of three independent

**Figure 7.** Sequence alignments and phylogenetic tree of DbMYB2 from orchids with known R2R3- MYBs domain. (**A**) Alignments of the full-length deduced amino acid sequences of DbMYB2 with

*Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 10 of 17

#### **3. 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*0*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–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.

#### **4. Materials and Methods**

#### *4.1. Plant Materials*

Protocorm-like bodies (PLBs) of *Dendrobium bigibbum* var. compactum were cultured on pH5.3 Hyponex medium (6.5:6.0:19:0 N:P:K; Hyponex Japan Corp., Ltd, Osaka, Japan) supplemented with sucrose (3% *w*/*v*) and agar (0.4% *w*/*v*) (Duchefa Biochemie B.V., Haarlem, The Netherlands). PLBs were cultured in 220 ml glass jars containing 30ml medium, which were closed with semipermeable plastic caps. All the cultures were maintained at 22–25 ◦C and >60% humidity. Plants were grown under white fluorescent light (PPFD = 50 µmol/m<sup>2</sup> /s) with a 16-h illumination and 8-h dark photoperiod.

## *4.2.* γ*-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.

#### *4.3. 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.

#### *4.4. 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

$$\begin{aligned} \text{Anthocyanin pigment (cyanidin } -3-\text{glucose equivalents, mg/L)}\\ = \frac{\text{A} \times \text{MW} \times \text{DF} \times 10^3}{\varepsilon \times 1} \end{aligned} \tag{1}$$

A = (*A*510 nm − *A*700 nm) pH1.0 − (*A*510 nm − *A*700 nm) pH4.5; MW (molecular weight) = 449.2 g/mol for cyanidin-3-glucoside (cyd-3-glu); DF = dilution factor; l = pathlength in cm; ε = 26,900 molar extinction coefficient, in L·mol−<sup>1</sup> ·cm−<sup>1</sup> , for cyd-3-glu; and 10<sup>3</sup> = factor for conversion from g to mg [43].

#### *4.5. 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].

#### *4.6. 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.

#### *4.7. 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.

#### *4.8. 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).

#### *4.9. 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).

#### *4.10. 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.

#### *4.11. 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,000*g* for 5 min) and resuspended to an OD<sup>600</sup> = 0.8 in a buffer containing 10 mM MES pH 5.6, 10 mM MgCl2, 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.

**Author Contributions:** Conceptualization, G.-H.L. Formal analysis, G.-H.L.; Funding acquisition, S.-Y.K. and J.-B.K.; Investigation, S.W.K.; Methodology, J.R.; Project administration, S.-Y.K., J.-B.K. and S.H.K.; Resources, S.H.K.; Writing–original draft, G.-H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a grant from the Nuclear R&D program of the Ministry of Science and ICT and the research program of KAERI, Republic of Korea.

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