Genome-Wide Transcriptome and Expression Profile Analysis of Phalaenopsis during Explant Browning

Background Explant browning presents a major problem for in vitro culture, and can lead to the death of the explant and failure of regeneration. Considerable work has examined the physiological mechanisms underlying Phalaenopsis leaf explant browning, but the molecular mechanisms of browning remain elusive. In this study, we used whole genome RNA sequencing to examine Phalaenopsis leaf explant browning at genome-wide level. Methodology/Principal Findings We first used Illumina high-throughput technology to sequence the transcriptome of Phalaenopsis and then performed de novo transcriptome assembly. We assembled 79,434,350 clean reads into 31,708 isogenes and generated 26,565 annotated unigenes. We assigned Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations, and potential Pfam domains to each transcript. Using the transcriptome data as a reference, we next analyzed the differential gene expression of explants cultured for 0, 3, and 6 d, respectively. We then identified differentially expressed genes (DEGs) before and after Phalaenopsis explant browning. We also performed GO, KEGG functional enrichment and Pfam analysis of all DEGs. Finally, we selected 11 genes for quantitative real-time PCR (qPCR) analysis to confirm the expression profile analysis. Conclusions/Significance Here, we report the first comprehensive analysis of transcriptome and expression profiles during Phalaenopsis explant browning. Our results suggest that Phalaenopsis explant browning may be due in part to gene expression changes that affect the secondary metabolism, such as: phenylpropanoid pathway and flavonoid biosynthesis. Genes involved in photosynthesis and ATPase activity have been found to be changed at transcription level; these changes may perturb energy metabolism and thus lead to the decay of plant cells and tissues. This study provides comprehensive gene expression data for Phalaenopsis browning. Our data constitute an important resource for further functional studies to prevent explant browning.

Our data constitute an important resource for further functional studies to prevent explant browning.

Background
Plant tissue culture is an important tool in both basic and applied studies as well as in commercial application, Such as for propagation of orchids, especially in production of Phalaenopsis. The Phalaenopsis genus of popular ornamental plants belongs to the family Orchidaceae and is mainly used in tissue culture because of its rapid propagation. However, explant browning in Phalaenopsis tissue culture presents a major problem in producing regenerated Phalaenopsis in culture. To date, little is known about the browning mechanisms of Phalaenopsis. Therefore, establishing effective and economical methods to reduce explant browning in tissue culture of Phalaenopsis will provide a valuable tool for plant propagation and transgenic manipulation.
Recent advances in sequencing technologies have enabled genomic-scale sequencing projects for many model organisms. These projects include de novo transcriptome analysis and reference mapping of expressed transcripts [14][15][16]. Using RNA-sequencing (RNA-seq) techniques, Mellidou et al [17] identified different expressed gene (DEGs) mainly involved in lipid metabolism, secondary metabolism, and cell wall modifications in apple fruit browning disorder. The energy-related and stress-related genes were also altered during apple fruit browning development. To investigate pear fruit surface brown (SBS), transcriptome analysis showed that up-regulated the expression of genes related to oxidative phosphorylation, phenolic compound synthesis and PPO [18]. These studies provide a genomics basis for botanists to understand the molecular mechanisms of enzymatic browning.
Recently, OrchidBase (http://lab.fhes.tn.edu.tw/est) has been established from 37,979,342 sequence reads collected from 11 in-house Phalaenopsis orchid cDNA libraries using multiple sequencing techniques [19]. A total of 1,233,823 unique sequences were obtained from Phalaenopsis aphrodite using Roche 454 and Illumina/Solexa high-throughput sequencing technologies [20]. Based on these accomplishments in Phalaenopsis genome research, in the current study, we carried out transcriptome analysis and expression profiling of Phalaenopsis leaf explants during browning.
Here, we used Illumina short-read sequencing for de novo transcriptome assembly and analysis of Phanalenopsis hybrid: Konggangjinli (A red flower cultivar, S1 Fig) explants at the early stage of browning. We constructed a mixed library from 0-d-cultured (control), 3-d-cultured (prior to browning, browning rate 0%), and 6-d-cultured explant (brown rate 100%) (S2 Fig). After sequencing, we used BLAST to compare these reads to the NCBI database and Orchid-Base (http://lab.fhes.tn.edu.tw/est) to determine their encoded proteins. The genes and isoforms were annotated and functionally mapped to Gene Ontology (GO) terms. Furthermore, we constructed 3 libraries from the above-mentioned three kinds of explant for expression profile analysis during Phalaenopsis explant browning. The reads were first mapped to the transcriptome and OrchidBase and then analyzed to identify differentially expressed genes (DEGs) among 0 d, 3 d and 6 d. DEGs were categorized into three groups: explants cultured for 3 d vs. 0 d (before explant browning), explants cultured for 6 d vs. 0 d (after explant browning), and explants cultured for 6 d vs. 3 d. Next, we confirmed the differential expressions of interesting genes by quantitative real-time polymerase chain reaction (qRT-PCR). This study sheds fresh light on explant browning at a genome-wide scale and will facilitate the functional analysis of browning-related genes in Phalaenopsis.

Materials and Methods
Plant material [21]. BLAST2GO was used to assign putative functionalities, GO terms, and KEGG (Kyoto Encyclopedia of Genes and Genomes) based metabolic pathways [22]. Furthermore we conducted sequence similarity analysis of protein sequences between the assembled transcripts and Arabidopsis thaliana, Nicotiana tabacum, Glycine max, Lycopersicon esculentum, Triticum aestivum, Oryza sativa, and equestris of the Phalaenopsis database. In this analysis, the Phalaenopsis equestris database was downloaded from the NCBI database (http://www.ncbi.nlm.nih.gov/ bioproject/53345), by BLASTn search with a cut-off (E-value of 10 -5 ). Arabidopsis and Oryza were chosen as they represent dicot and monocot model species, respectively.

Illumina PE library preparation for expression profile analysis and transcriptome mapping for bioinformatics analysis
To examine expression profiles, we then constructed and sequenced three independent libraries from 0 d, 3 d, and 6 d explants. Each cDNA library was constructed as described above and then sequenced by Illumina HiSeq 2000 (1 × 50 bp read length). After removing adaptor sequences, empty reads, and low-quality sequences with unknown sequences "N", the high quality clean reads were mapped to our Phalaenopsis transcriptome reference database using the Bowtie software with default parameters. Then the transcript abundances were quantified by SEM (http://dewelabbiostat.wise.edu/rsem). The differentially expressed genes (DEGs) between the three samples were identified using the number of mapped reads as EdgeR inputs (http://bioconductor.org/packages/release/bioc/html/edgR.html). In this algorithm, we used FDR (False Discovery Rate) 0.05 and the absolute value of |log2FC|>1 as the threshold to judge the significance of differential gene expression. Here, FDR was used to determine the threshold p-value in multiple tests and analysis through manipulating the FDR value, where p-value corresponded to the differential gene expression test.
Three comparisons were conducted: explants cultured for 3 d vs. explants cultured for 0 d (3 d vs. 0 d), explants cultured for 6 d vs. explants cultured for 0 d (6 d vs. 0 d), and explants cultured for 6 d vs. explants cultured for 3 d (6 d vs. 3d). A Venn diagram was constructed using the Venny tool software [23].
The three groups of DEGs were mapped to GO terms. The DEGs were aligned to the GO database for GO functional enrichment analysis using GO tools software with FDR corrections (http://githb.com/tanghaibao/goatools). Statistical significance was considered at P<0.05. We also mapped all DEGs to the KEGG database and searched for significantly enriched KEEG pathways at P<0.05 level using KOBAS software (http://kobas.cbi.pku.edu.cn/home.do). The sequence similarity analysis of protein sequences between the DEGs with Arabidopsis thaliana and Oryza sativa were also conducted as mentioned above. Homologous protein domains for DEGs were identified by searching against the Pfam database using Hmmscan Validation of gene expression by quantitative real-time PCR (qRT-PCR) qRT-PCR was performed to confirm the RNA-seq data. Total RNA was extracted for cDNA preparation as described above. Total RNA (1 μg) was reverse transcribed into single-stranded cDNA using the Primescript RT reagent kit (TaKaRa, Dalian, China). qRT-PCR was implemented using the SYBR premix Ex Taq kit (TaKaRa, Dalian, China), with first-strand cDNA as the template. The actin gene from Phalaenopsis was used as an internal control. The relative quantitative method (2 -ΔΔCt ) was used to calculate the fold change of the target genes. The primers employed in the qRT-PCR are listed in S1 Table. The results of these reference genes were compared to fragments per kilobase of exon model per millions mapped reads (FPKM) estimates.

Phalaenopsis transcriptome assembly and annotation
To examine the Phalaenopsis transcriptome, we collected RNA from 0-d-cultured, 3-d-cultured, and 6-d-cultured explants and constructed one mixed cDNA library. We used the tagged cDNA library for 2×100 bp paired-end sequencing on a single lane of the Illumina Hiseq2000 and generated 118,996,000 raw paired-end reads resulting in 12 gigabases of sequence. After cleaning and quality control, this resulted in 79,434,350 reads (7,544,398,171 bp) ( Table 1). Using the RNA seq de novo assembler Trinity, we obtained 31,708 isogenes with a mean size of 1108.07 bp and lengths from 351 bp to 9,265 bp (S2 Table). Using BLASTn sequence similarity search with cut-off E-value 10 -5 , about 90% of isogenes (28,267/31,708) have significant sequence similarity to the Phalaenopsis EST database (data not shown). Only 10% of isogenes (3441) were not hit, and were thought to be specific sequence in our study. We used this reference transcriptome for our subsequent analysis of expression profiles during explant browning.
About 84% assembled transcripts (26,565/31,708, S3 Table) were annotated by using BLASTp and BLASTn algorithms to search against the NR, gene, string and Nt databases (E-value 10-5). A BLASTx sequence similarity search with cut-off E-value 10 -5 , showed that Phalaenopsis proteins have significant sequence similarity to the proteomes of other plants ( Table 2). For example the Phalaenopsis proteome showed substantial amino acid sequence identity to Arabidopsis thaliana (61.71%) and Oryza sativa (63.78%) proteomes. The function of the Phalaenopsis sequences hit in Arabidopsis thaliana and Oryza sativa were shown as supporting information in S4 and S5 Tables. Most of the transcripts had most-significant hits to the Phalaenopsis equestris database, and only 14.34% of the sequences did not show significant identity in the Phalaenopsis equestris database. These do not found homologous sequences suggesting that they may novel proteins in our study.

Phalaenopsis expression profile changes during explant browning
To examine the expression of Phalaenopsis genes in browning, we next constructed three cDNA libraries and used 2×100 bp paired-end sequencing to obtain sets of expression reads for explants cultured for 0, 3, and 6 d. After removing the adaptor sequences and low-quality sequences, we mapped the clean reads for each sample (0, 3, and 6 d) to the assembled reference transcriptome of Phalaenopsis (Table 3). We used Bowtie (http://bowtie-bio.sourceforge.net/index.shtml) to map reads, used RSEM to assign reads to genes and isoforms, and used edgeR to standardize the reads. We refer to the normalized reads in fragments per kilobase of exon model per millions mapped reads (FPKM), and we used the number of genes with FPKM values !0.5 to estimate the total number of genes expressed during Phalaenopsis explant browning(S6 Table). We found 19,156 expressed genes in the explants cultured for 0 d, representing almost 60% of the annotated transcriptome. We also found 22,185 genes (70%) were expressed in 3 d, and 22,902 (72%) in 6 d ( Table 3).

Identification of genes differentially expressed during Phalaenopsis explant browning
We compared our three time points (0, 3, and 6 d) and identified 698 differentially expressed genes (DEGs) (Fig 1, S7 Table) with FDRs of <0.05. We identified 267 DEGs before Phalaenopsis explant browning in the comparison of 3 d vs. 0 d, 169 (644 isoforms) up-regulated and 98 (523 isoforms) down-regulated. We identified 534 DEGs after explant Phalaenopsis browning in the comparison of 6 d vs. 0 d, 314 (921 isoforms) up-regulated and 220 (641 isoforms) down-regulated. Finally, we only identified 159 DEGs in the comparison of 6 d vs. 3 d, 76 (748 isoforms) up-regulated and 83 (553 isoforms) down-regulated. These results showed that many more genes increased in expression before and after explant browning. Also, the proportion of DEGs was highest in the 0 d vs. 6 d comparison and lowest in the 3 d vs. 6 d comparison.
To identify shared and unique DEGs during explant browning, we generated a Venn diagram based on the three comparisons (Fig 2). Only 11(44 isoforms) DEGs were shared in the three groups. Furthermore, DEGs that have the best homolog in Arabidopsis thaliana and Oryza sativa were also studied (S8, S9 and S10 Tables). In the comparison of 3 d vs. 0 d, 6 d vs. 0 d 6 d vs. 3 d, there were 31, 49and 38 isoforms shown sequences identity in Arabidopsis thaliana, and  173, 239,173 isoforms hit in Oryza sativa respectively. The number of homolog in three groups showed low identity in Arabidopsis thaliana comparison with Oryza sativa. Based on the homologs in Arabidopsis thaliana, we found multiple DEGs involved in biosynthesis of secondary metabolites, nucleotide metabolism, carbohydrate metabolism, and amino acids metabolism, such as: NAD (P)-binding rossmann-fold superfamily protein, nucleoside diphosphate kinase family protein, protein kinase superfamily protein, glucose-6-phosphate dehydrogenase 2, caffeoyl-CoA 3-O-methyltransferase.
The highly expression of DEGs before and after Phalaenopsis explant browning was shown in Table 4. The results showed that the expression of pathogenesis-related protein 10c was in higher level and was significantly elevated after explant browning. Notably, glutathione S-transferase (comp16940_0, comp23612_0, comp18470_0, comp23433_0), which involved in Glutathione metabolism, increased by 10-700 times after explant browning, It was also worth mentioning that the expressions of nonsymbiotic hemoglobin (comp23970_0), bibenzyl synthase (chalcone synthase, comp22408_0), 2-hydroxyacyl-CoA lyase (comp22023_0), Sadenosyl-L-methionine synthase (comp19689_0), cytochrome P450 71D11-like (comp20309_0), which appear to be involved in peroxisome, cysteine and methionine metabolism and secondary metabolism, respectively, were very low in the control samples, but increased significantly after Phalaenopsis explant browning occurred, the expression increase by 30-600 folds, of which nonsymbiotic hemoglobin expression increase more than 300 folds. These results demonstrated that there is a close relationship between metabolism and browning. The expression of RNA-dependent RNA polymerase (comp23180_c0), was in high level in control explant and was significantly decreased before and after explant browning. Three DEGs involved in energy metabolism, such as V-type proton ATPase 16 kDa proteolipid subunit-like (comp 16560_c0), which take part in oxidative phosphorylation, putative chlorophyll a/b-binding protein (comp 23175_c0), type III chlorophyll a/b-binding protein (comp20543_c0) which relate to photosynthesis expression reduced before and after explant browning, respectively.
To identified DEGs that were induced before and after Phalaenopsis explant browning, the result showed that there were 36 and 69 DEGs induced before and after Phalaenopsis explant browning, respectively (S11 Table). The top 10 DEGs were showed in Tables 5 and 6. We found that cytokinin oxidase(comp21576), probable mannitol dehydrogenase (comp21741),4-coumarate:CoA ligase(comp20466), which involved in zeatin biosynthesis, phenylpropanoid biosynthesis, fatty acid elongation, respectively, were highly induced before and after Phalaenopsis explant browning.

Gene ontology (GO) analysis for Phalaenopsis transcriptome and DEGs
To characterize the Phalaenopsis transcriptome, we used GO assignments to classify the functions of the predicted Phalaenopsis genes. Using Blast2GO, we grouped the 21,384 annotated sequences into 45 sub-categories within the three main GO categories of biological process, cellular component and molecular function (Fig 3). Based on the biological processes, we established 15 categories. We found 6,318 genes (52.99%) involved in cellular processes, 6,167 genes (51.72%) in metabolic processes and 2, 420 genes (20.29%) in responding to stimulus. In addition, we found 14 clusters of genes with roles in cell formation, including cell parts (8,275 genes, 69.41%), organelles (6,885 genes, 57.75%), and membrane (3,431 genes, 28.79%). Based on their putative molecular functions, we also found 23 categories, including genes involved in catalytic activity (5,642 genes, 47.32%) and binding activity (5, 534 genes, 46.42%).
Once we had characterized the transcriptome, we next assigned GO categories to the DEGs and compared them with those in the transcriptome. We found a similar distribution of GOannotation function of DEGs in the three comparisons (S3 Fig). For biological processes, before and after explant browning, most DEGs were assigned to metabolic process, followed by cellular process and responding to stimulus. For cellular component, the top three categories for DEGs were cell part, cell, and organelles. For molecular functions, the top categories for DEGs were catalytic activity and binding. Furthermore, GO enrichment analyses showed that before Phalaenopsis explant browning, the DEGs for 3 d vs. 0 d mainly function in the protein-DNA complex, and nucleosome organization. However, after Phalaenopsis explant browning, the DEGs for 6 d vs. 0 d mainly function in metabolic processes, such as small molecule metabolism, oxoacid metabolism, and carboxylic acid metabolism (S12 Table). We also analyzed the GO enriched functional categories in each of the three pairwise comparisons, which condensed and compressed the genes by removing categories that did not show as significantly different, and uses a false color heatmap to indicate up-or down-regulated classes (Fig 4). Genes encoding enzymes for cytokinin dehydrogenase, oxidoreductase, acting on CH or CH2 groups, with an iron-sulfur protein as acceptor, 4-hydroxy-3-methylbut-2-en-1-yl diphosphate synthase activity were strongly   up-regulated in 3 d and 6 d explants compared to control (0 d) (Fig 4). Also, we found no significant change for 6 d vs. 3 d in GO function enrichment analyses.

Pathway analysis for Phalaenopsis transcriptome and DEGs
Pathway-based analysis can help elucidate the biological functions of genes. To identify the biological pathways that are active in Phalaenopsis, we mapped all the unigenes to the reference canonical pathways in KEGG. A total of 14,099 unigenes were assigned to 123 KEGG pathways (S13 Table), where 'Metabolic pathways' (1,975 unigenes, 14%) is dominant, followed by 'Biosynthesis of secondary metabolites '(791 unigenes, 5%) and 'Spliceosome '(309 unigenes, 2%) ( Fig 5). KEGG classification analysis identified transcripts involved in signal transduction, carbohydrate metabolism and lipid metabolism have a high degree in Phalaenopsis (Fig 6). The distribution of KEGG pathway of Phalaenopsis genes was shown in Fig 7. The top 30 pathways of Phalaenopsis transcripts mapped to the Arabidopsis thaliana and Oryza sativa KEGG were showed that the purine metabolism, spliceosome, pyrimidine metabolism and RNA transport were high hit in Oryza sativa, while metabolic pathways and biosynthesis of secondary metabolites were high mapped in Arabidopsis thaliana. About 293,230 and 214 transcripts mapped in the spliceosome, RNA transport and ribosome in Phalaenopsis were tooken over the top 3 (Fig 7).
To characterize the functional consequences of gene expression changes associated with explant browning, we performed pathway enrichment analysis on DEGs based on the KEGG database with p <0.05 as the threshold. We found that 18 of 90 pathways demonstrated significant changes in 3 d vs. 0 d, 25/102 pathways in 6 d vs. 0 d, and 11/98 pathways in 6 d vs. 3 d (S14 Table). Comparison of the three groups of DEGs revealed that phenylpropanoid biosynthesis, flavonoid biosynthesis, phenylalanine metabolism and terpenoid backbone biosynthesis demonstrated significant changes in all three comparisons. There were seven pathways of carbohydrate metabolism, seven pathways of amino acid metabolism, five of lipid metabolism, four pathways of secondary metabolism significantly changed before and after Phalaenopsis explant browning (Table 7).

Protein domain analysis for Phalaenopsis transcriptome and DEGs
Conserved protein domains were identified in the Phalaenopsis unigenes against the Pfam database by using Hmmscan software with cutoff E-value 0.1, 12,341 of 26,565 assembled unigenes matched entries that are corresponding to 3,021 different domains/families (S15 Table). The 20 most abundant domains identified are provided in Table 7. The top 3 abundant domains were identified included Pkinase (protein kinase domain) with 568 unigenes, followed by Pkinase_Tyr (protein tyrosine kinase) with 284 unigenes and RRM_1(RNA recognition motif) with 202 unigenes. The PPR domain proteins and cytochrome P450 families were also highly represented.
To specifically identify DEGs protein domains in the Pfam database, comparison of Pfam searches between three groups of DEGs indicated an overall similarity of identified domains. The top 20 domains were shown in Fig 8. Protein kinase domain, cytochrome P450, methyltransferase domain, Protein tyrosine kinase and ABC transporter were most highly represented in three groups of DEGs.   Validation of the expression profile data using real-time PCR Based on KEGG pathway enrichment analysis and previous studies [4][5][6][11][12], we selected eleven DEGs and subjected them to quantitative real-time PCR (qRT-PCR) to validate the results of the expression profile analysis (Fig 9). We found a high correlation (R 2 = 0.8841) between RNA-seq and qRT-PCR results (S1 Table). For example, the transcript levels of CHS (chalcone synthase), 4CL(4-coumaroyiÁCoA ligasc), and F3'H (flavonoid3'-hydroxylase)were high in 6 d explants, and DFR (dihy-droflavonol 4-reductase)decreased in 3 and 6 d explants (Fig 9h and 9i). ATPase γ subunit, ATPase α subunit, psbB, and psbD were down-regulated in the 3 and 6 d explants, consistent with the results of expression profile analysis. Only the psbB transcript level increased in the 6 d explants, in contrast to the sequence data (Fig 9e and 9k).
As expected [4,11], the transcripts of PPO, POD, and PAL were up-regulated during browning (Fig 9I). Our results indicate that our RNA-seq data reliably identified browning-related genes in Phalaenopsis.

Transcriptome and expression profile analysis of DEGs during Phalaenopsis explant browning
Using the Illumina HiSeq 2000 next generation sequencing technology, we obtained transcriptome sequences from explants that were cultured for 0, 3, and 6 d. We assembled more than 79,434,350 reads into approximately 31,708 isogenes and annotated 26,565 proteins. Homology searches showed that 63.78 and 61.71% isogenes have significant similarities to Oryza sativa and Arabidopsis thaliana protein sequences and about 85.66% of isogenes show identity with the publicly available Phalaenopsis equestris. Our data show more sequences with significant identity to known sequences in the plant proteome than the sequences identified by Hsu et al. [24]. Also, we annotated and functionally mapped 21, 348 sequences to GO terms. Gene annotation can help predict potential genes and their functions at the transcriptome level. In this study, we found predominant gene clusters involved in the cellular and metabolic processes of biological processes, the binding and catalytic activities of molecular functions, as well as the cell, cell part, and organelles categories of cellular components. Our results were different from those obtained in Phalaenopsis aphrodite, in which the identified genes were predominantly involved in the biosynthesis of macromolecules and in nitrogen, protein, and nucleic acid metabolism [20]. This difference may be due to the fact that we constructed our cDNA library only from Phalaenopsis leaf tissue. The DEGs among the explants cultured for 0, 3, and 6 d probably relate to their metabolism and functions. We found that 267 and 534 genes exhibited significant differences in expression before and after (3d and 6d) Phalaenopsis explant browning, respectively. Among them 36 and 69 DEGs were induced before and after explant browning, respectively. The three comparisons shared only 11 genes in common, indicating that the patterns of gene expression differed over time during explant browning. The identified DEGs were involved in biological processes, cellular components, and molecular functions GO terms. Pfam search showed that protein kinase domain, cytochrome P450, methyltransferase domain, protein tyrosine kinase and ABC transporter were most highly represented in three groups of DEGs. Protein kinase and its subcategories tyrosine-protein kinase were known to regulate the majority of cellular pathways. Cytochrome P450 is a large plant protein family that might contribute to extensive modifications of various secondary compounds. These results revealed DEGs with multiple functions, suggesting that 3 d and 6 d were critical stages for explant browning, when many substances related to browning were produced. DEGs involved in polyphenol synthesis and phenylpropanoid pathways before and after Phalaenopsis explant browning Previous studies indicated that phenol oxidase activity results in explant browning and the degree of explant browning positively correlates with the contents of phenolic compounds [2][3][4][5]25]. The accumulation and oxidation of phenolic compounds causes enzymatic browning. Also, reduction of polyphenol synthesis and inhibition of the phenylpropanoid pathway repressed browning [26][27][28]. PAL is the first dedicated enzyme in the phenylpropanoid pathway and converts phenylalanine into trans-cinnamic acid, providing a substrate for further synthesis of phenolic compounds. PPO and POD also function in enzymatic browning. Both enzymes have common diphenolic substrates leading to brown products. In our previous studies, we found that tannin filled in the vascular bundles of browning tissues [29]. Based on previous work and this study, we selected seven genes for qPCR: PhCHS, Ph4CL, PhF3'H, and PhDFR, which function in flavonoid biosynthesis, and PhPPO, PhPOD, and PhPAL, which function in phenylalanine metabolism, and which also showed significant changes in the KEGG pathway enrichment analysis in this study. Consistent with our expectations, PhPPO, PhPOD, and PhPAL were up-regulated in Phalaenopsis explant browning [2,4,8,30] and qRT-PCR confirmed the expression profiling obtained in this study.
Anthocyanins, a class of flavonoids that is responsible for the colors in fruits and most flowers of higher plants, are major water soluble pigments [31,32] and exhibit important physiological functions, such as antioxidative [33,34]. The degradation of anthocyanins and/or the oxidation of phenolics caused by polyphenol oxidase (PPO) results in an enzymatic browning reaction of fruits and vegetables [35,36]. Genes in anthocyanin biosynthesis such as PhPAL, PhCHS, Ph4CL, PhF3'H were up-regulated in 6 d explants of Phalaenopsis indicates that anthocyanin may play important roles in the early stages of Phalaenopsis explant browning. Previous studies have shown that disruptions of oxidation and reduction may induce tissue browning, leading to reactive oxygen species (ROS) bursts, which would affect the membrane integrity [37][38][39]. ROS which were induced in tissue culture is eliminated by biosynthesis of anthocyanin. DFR is a key enzyme in the catalysis of the stereospecific reduction of dihydroflavonols to leucoanthocyanidins in anthocyanin biosynthesis. Over-expression DFR in the sweet potato enhanced scavenging of reactive oxygen radicals in plants under stressful conditions [40]. PhDFR showed opposite expression profile compared to that of PhPAL, Ph4CL, PhCHS and PhF3'H in the early stages of Phalaenopsis explant, so anthocyanin biosynthesis may be effected. This may lead to imbalance of ROS metabolism and membrane lipid peroxidation, and destroy of cell membrane integrity.
Meanwhile, Cytochrome P450 plays an important role in flavone biosynthesis and has been well characterized. About 115 different cytochrome P450 isogenes were detected in our RAN sequencing data. The cytochrome P450 71D11-like was highly expression after Phalaenopsis explant browning (Tables 4 and 8). Cytochrome P450 expressed only after lotus tissue browning [41], indicating that browning was beneficial for secondary metabolite production.
DEGs involved in photosynthesis, energy synthesis and other metabolism before and after Phalaenopsis explant browning Photosystem II (PSII) is a multi-subunit, pigment-protein complex localised in the chloroplast thylakoid membranes. The PSII complex is composed of a number of proteins including D1 (PsbA), D2 (PsbD), CP43 (PsbC), CP47 (PsbB), apoproteins of cytochrome b559, and the extrinsic manganese-stabilizing protein, along with a number of small proteins [42]. In our KEGG pathway enrichment analysis, we found that DEGs involved in photosynthesis also changes significantly in browning explants, so we selected PhpsbB, and PhpsbD for qRT-PCR. F-ATPases (also known as F1F0-ATPase, or H(+)-transporting two-sector ATPase) (EC:3.6.3.14) function as the prime producers of ATP in mitochondria, chloroplasts and bacterial plasma membranes, using the proton gradient generated by oxidative phosphorylation (mitochondria) or photosynthesis (chloroplasts). Thus, we also selected the genes encoding the PhATPase γ subunit, and APhTPase α subunit for qRT-PCR. Our qRT-PCR results were in accordance with the RNA-seq expression profiles in which the transcript levels decreased in the 3 and 6 d explants. By contrast, our qRT-PCR analysis showed that the transcript levels of PhpsbB increased in the 6 d explants. We have found that during explant browning, chloroplasts shrink to dark-colored structures and chlorophyll b content declines [25,29]. In this study, we found that the abundance of PhpsbB and PhpsbD declined before and after Phalaenopsis explant browning. Consistent with this results, two kinds of chlorophyll a/b binding protein were found to be decreased significantly during explant browning (Table 4). To our knowledge, this is the first report to show that photosynthesis genes are associated with explant browning. Liu et al [18] also confirmed that photosynthesis changes during pear fruit browning. These metabolic processes will disrupt the energy metabolism providing a molecular context for the formation of ROS in browning tissue. We conceive that the wound-induced natural defense response leads to accumulation of toxic compounds that impair photosynthesis in browning plant tissues and ultimately damage or kill plant cells.
In plant tissues, ATP content is important for regulating anti-oxidase expression, such as superoxide dismutase, catalase, and ascorbate peroxidase [43]. Jiang et al. [40] also found that ATPase was differentially expressed in Lotus tuber browning tissue using 2-DE and MALDI-TOF-TOF for protein profile analysis. Proteins related to nucleotide metabolism: adenosine kinase, and nucleoside diphosphate kinase were down-expressed after Lotus tissue browning. Their down-regulations indicated that substance metabolism was out of control after browning and the abilities of energy synthesis and utilization declined, accompanied by disorder of ATP energy metabolism. Senescence and stress can also lead to lower ATP production and higher accumulation of reactive oxygen species, resulting in membrane system disorder and triggering browning [42]. Moreover, using LCQ MS/MS for protein expression analysis in explants cultured for 3 d, we linked ATPase with explant browning [12]. In line with this, we herein verified that the transcript levels of the genes encoding the PhATPase γ and PhATPase α subunits decreased in the 3 d and 6 d explants, implying that PhATPase has important functions during explant browning. The expression of V-type proton ATPase 16 kDa proteolipid subunit-like also decreased significantly before and after explant browning (Table 4). Down-regulation of genes encoding ATPase subunits indicated that energy synthesis and utilization declined. Less energy means incapacity to repair membranes, in turn, loosening of cellular compartmentalization allows the release of PPOs and phenols to the cytosol, triggering the actual discolouration reactions. Thus, we consider that ATP metabolism might affect explant browning in Phalaenopsis.
Before and after Phalaenopsis explant browning, DEGs involved in amino acid metabolism, carbohydrate metabolism, lipid metabolism were significantly changed (Table 7). A total of 22 and 44 DEGs involved in amino acid metabolism before and after explant browning, indicating that expressions of some proteins resulted in promotion of the browning reaction. After Lotus tissue browning, Jiang et al found the expression of proteins relating to amino acid metabolism: cysteine synthase and glutamate-ammonia ligase reduced, respectively [40]. Related to lipid metabolism the expression type 1 non-specific lipid transfer protein reduced during explant browning, implying that lipid might be related to plant stress tolerance. Many DEGs involved in lipid degradation pathways were significantly up-regulated in apple fruit browning disorder [17]. It is understandable that changes in these metabolic processes will disrupt the energy metabolism and cell compartments, resulting oxidation of phenolic compounds and influence the explant browning.

Conclusions
Our RNA-seq analysis of changes in gene expression substantially increases the number of available equestris in Phalaenopsis and provides comprehensive information on the mechanism of Phalaenopsis explant browning at a genome-wide scale. Most of the DEGs we identified were up-regulated during Phalaenopsis explant browning. More importantly, transcripts of the genes that control secondary metabolism, energy metabolism, carbohydrate metabolism, and lipid metabolism were differentially expressed in explant browning. In addition, in this study, we identified many new browning-related genes. More detailed and extensive analyses of the identified genes will be essential to further decipher their possible roles in explant browning. Thus, our data provide genome-wide information that enables further exploration of explant browning.