Transcriptome analysis to identify candidate genes associated with the yellow-leaf phenotype of a Cymbidium mutant generated by γ-irradiation

Leaf color is an important agronomic trait in flowering plants, including orchids. However, factors underlying leaf phenotypes in plants remain largely unclear. A mutant displaying yellow leaves was obtained by the γ-ray-based mutagenesis of a Cymbidium orchid and characterized using RNA sequencing. A total of 144,918 unigenes obtained from over 25 million reads were assigned to 22 metabolic pathways in the Kyoto Encyclopedia of Genes and Genomes database. In addition, gene ontology was used to classify the predicted functions of transcripts into 73 functional groups. The RNA sequencing analysis identified 2,267 differentially expressed genes between wild-type and mutant Cymbidium sp. Genes involved in the chlorophyll biosynthesis and degradation, as well as ion transport, were identified and assayed for their expression levels in wild-type and mutant plants using quantitative real-time profiling. No critical expression changes were detected in genes involved in chlorophyll biosynthesis. In contrast, seven genes involved in ion transport, including two metal ion transporters, were down-regulated, and chlorophyllase 2, associated with chlorophyll degradation, was up-regulated. Together, these results suggest that alterations in chlorophyll metabolism and/or ion transport might contribute to leaf color in Cymbidium orchids.


Introduction
Orchids such as Cymbidium, Dendrobium, Oncidium, and Phalaenopsis are important cash crops worldwide, and the orchid industry has contributed substantially to the economy of many Southeast Asian countries [1,2]. Because of its fragrant flowers and straight leaves, Cymbidium is a popular orchid in China, Korea, and Japan [3][4][5]. In addition to its floral diversity, PLOS

Quantitative real-time PCR (qRT-PCR) analysis
Reverse transcription (RT) and first-strand cDNA synthesis were carried out using a ReverTra Ace-α-kit (Toyobo Co. Ltd, Osaka, Japan). Quantitative RT-PCR was carried out with iQ SYBR Green Supermix (Bio-Rad, Hercules, CA, USA). Quantitative PCR was performed using CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, USA). The PCR conditions included and initial denaturation step of 95˚C for 10 min followed by 35 cycles of 95˚C for 15 s, 56-62˚C for 15 s and 72˚C for 30 s. Each sample was run in triplicates and ACTIN was used as an internal control. Cycle threshold values were calculated using Bio-Rad CFX Manager 3.1 software (Bio-Rad, Hercules, USA). Gene-specific primers used for quantitative RT-PCR are described in S1 Table. Chl and carotenoid content assay Six-month-old leaves of the WT and the S12 mutants were sampled for Chl and carotenoid determination. Amounts of Chls a and b and carotenoids were estimated following the method of Lichtenthaler [34]. The fresh leaf samples were ground using liquid nitrogen and suspended in 96% ethanol (Sigma, St. Louis, MO, USA). This extract was vortexed and placed at room temperature in dark for 24 h. The absorption spectra of the extract was measured at 470, 648.6, and 664.2 nm using a UV-1800 spectrophotometer (Shimadzu, Kyoto, Japan).

RNA sequencing and de novo assembly
The cDNA libraries were prepared from both the WT and the S12 leaves. Raw reads were trimmed by filtering out adaptor-only nucleotides that were smaller than 75 bp using Trimmomatic (v0.32) [35]. De novo assembler Trinity (v2.2.0) was used to construct large contigs from the filtered reads [36]. Trinity is a representative RNA assembler based on the de Bruijn graph algorithm. The assembly pipeline of Trinity consists of three consecutive modules: Inchworm, Chrysalis, and Butterfly. Protein coding sequences (CDSs) were extracted from the reconstructed transcripts using TransDecoder (v3.0.1), a utility included with Trinity to assist with the identification of potential coding region [37]. The prediction of coding regions is based on search of all possible CDSs, verification of the predicted CDSs by GENEID [38], and selecting the region that has the highest score among candidate sequences.

Functional annotation
Trimmed read were annotated with protein databases including gene ontology (GO), database for annotation, visualization, and integrated discovery (DAVID) [39], Kyoto Encyclopedia of Genes and Genomes (KEGG), and eukaryotic clusters of orthologous genes (KOG) databases. GO includes biological processes (BP), cellular component (CC), and molecular function (MF). KEGG is a major recognized pathway-related database that integrates genomic, biochemical, and systemic information. DAVID is a web-accessible program that integrates functional genomic annotations with intuitive graphical displays, which highlights pathway members within the biochemical pathways provided by KEGG. KOG categories were obtained via comparisons to the KOG database using RPS-BLAST (included with BLAST v2.2.26) [40].

Identification of differentially expressed genes between wild-type and S12 mutant
Gene expression profiles were analyzed using the method of RNA-Seq following Expectation Maximization (RSEM) [41]. The unique feature of RSEM is that it does not rely on a reference genome. RSEM uses the Bowtie alignment program to align transcripts. Differentially expressed genes (DEGs) were identified using edgeR [42], a Bioconductor package based one the generalized linear model that analyzes RNA-Seq data by considering gene expression as a negative binomial. We used a false discovery rate (FDR) < 0.05 significance cut-off for multiple testing adjustments [43].

Reduced accumulation of Chls and carotenoids in the S12 mutant
The S12 mutant developed yellow-colored leaves from the seedling stage to maturity (Fig 1).
To determine if this was associated with reduced Chl or carotenoid content, we quantified Chl a and Chl b as well as total carotenoids. A significant reduction in Chl and carotenoid contents was observed in the S12 mutant leaves; Chl a, Chl b, and carotenoid levels were reduced by 85%, 78%, and 65%, respectively (Fig 2A and 2C). The ratio between Chl a and Chl b also differed nominally but significant difference between the WT and the S12 mutants. The Chl a/b ratio serves as a useful indicator of plant response to shading [44] (Fig 2B).
To determine molecular changes underlying the S12 mutation, we assayed changes in genome-wide transcript levels using RNA-Seq. To this end, equal amount of RNAs from the  WT and the S12 mutants were used to construct cDNA libraries, and then sequenced by Trimmomatic (v0.32) platform. The average length of clean reads was 308-315 bp and a total of 30.83 million (99.92%) and 27.45 million (99.94%) clean reads were generated from the WT and the S12 mutant, respectively. A total of 225,694 transcripts were assembled into 144,918 unigenes with an average length of 1,257 bp (Table 1). We used gene ontology (GO) classification to classify the predicted functions of unigenes, which were categorized into 73 functional groups (FDR < 0.05). GO assignments were divided into three categories including biological process (BP), cellular component (CC), and molecular function (MF). Integral component of membrane (15.83%), plasma membrane (13.07%), and cytoplasm (12.36%) were dominant groups in the CC category followed by chloroplast (11.93%), cytosol (8.21%), and membrane (6.07%). Predicted proteins assigned to BP category were mainly associated with transcription (6.52%), protein phosphorylation (3.51%), protein ubiquitination (2.21%), and embryo development (1.87%). In the MF category, the most heavily represented groups were linked to ATP binding (9.99%), protein binding (9.09%), metal ion binding (5.74%), and protein serine/threonine kinase activity (3.84%) (Fig 3).  Next, we mapped the assembled unigenes to the reference canonical pathway in the KEGG, including metabolism, genetic information processing, environmental information processing, and cellular processes (http://www.kegg.jp/kegg/pathway.html). The 144,918 unigenes were assigned to 22 KEGG sub-pathways (Fig 4). Theses pathways included KEGG orthology (KO) entries for metabolism (734 KOs), genetic information processing (69 KOs), environmental information processing (14 KOs), and cellular processes (108 KOs) ( Table 2).

The implicated role of ion transport and Chl catabolism in the S12 phenotype according to DEG analysis
Gene expression analysis identified a total of 2,267 DEGs (FDR < 0.05) between the WT and the S12 mutant. Among these genes, 724 genes were up-regulated and 529 were down-regulated ( Fig 5). The DEGs were categorized into 27 functional groups in GO classification (FDR < 0.05). Predicted proteins assigned to biological process (BP) were mainly associated with single-organism process, which corresponded to the largest group. Cell and membrane terms were dominant among cellular components (CC). Those assigned to molecular function (MF) were mainly linked to ATP binding and transport activity (Fig 6). Especially, when DEG GO and total GO were compared in CC and MF, percentage of membrane group and transporter activity group were higher in DEG GO.
The DEGs were further classified into 22 functional categories using the KOG system. The largest KOG cluster was carbohydrate transport and metabolism (G), followed by secondary metabolites biosynthesis, transport and catabolism (Q), post-translational modification (O), energy production and conversion (C), and general function (R) (Fig 7). The KOG categories were further divided into multiple classes including cellular processes and signaling (273), information storage and processing (90), metabolism (705), and poorly characterized (135) classes (S2 Table). Next, we analyzed gene pathway assignment of DEGs to KEGG pathway. As a result, DEGs were assigned to seven sub-pathways (Fig 8), which in turn served as a template for assaying specific metabolic processes associated with the S12 mutant. To understand the biological roles of the screened major targets, we used DAVID to perform the GO biological process enrichment analysis. The DAVID, the most widely used on-line tool for functional classification, relies on a partitioning approach that groups genes together on the basis of their functional similarities. As predicted, a major category of DEGs was involved in membrane and chloroplast functions (Table 3). These DEGs included ABC transporters that were up-regulated in the S12 mutant and genes associated with ion transport were down-regulated  (Tables 4 and 5). Interestingly, the metal ion transporters have been associated with leaf color and photosynthesis [45,46]; their encoding genes were thus considered as candidates for the genes responsible for the phenotype seen in the S12 mutant. We next identified genes associated with Chl metabolism to assay association between transcript and Chl levels. A previous study in the model plant Arabidopsis revealed that 16 genes are involved in the conversion of glutamyl-tRNA to Chl [47]. In addition, genes associated with several major Chl catabolites have been identified, including Chlorophyll a-b binding protein 13 (CAB13), non-yellow coloring 1 (NYC1), NYC1-like (NOL), chlorophyllases (CLHs), pheophorbide oxygenase (PAO), and red Chl catabolite reductase (RCCR) [48]. No significant differences in expression levels were observed between the WT and the S12 mutant for genes involved in Chl biosynthesis. However, elevated expression levels of CLH2, a gene involved in Transcriptome analysis to identify genes associated with the yellow-leaf phenotype in Cymbidium mutant  Chl degradation, was observed in the S12 mutant (Fig 9). To validate the RNA-Seq data, we used qRT-PCR to measure the expression levels of 16 genes, including genes associated with ion transport and Chl biosynthesis and degradation. The qRT-PCR analysis confirmed~18 fold higher levels of CLH2 in the S12 mutant. In comparison, CAB13 and seven ion transporters, including two metal ion transporters, were down-regulated in the S12 mutants, by 1.65-and 1.3-4.05-fold, respectively (Fig 10O and 10A-10I). The expression levels of five other ion Transcriptome analysis to identify genes associated with the yellow-leaf phenotype in Cymbidium mutant Transcriptome analysis to identify genes associated with the yellow-leaf phenotype in Cymbidium mutant transporter genes were~4.5-52.3-fold higher in the S12 mutants (Fig 10J-10N). Together, these results suggest that the reduced levels of Chl in the S12 mutant could be a result of increased Chl degradation and/or impaired ion transport associated with Chl biosynthesis.

Discussion
In this study, we used the γ-ray-based mutagenesis procedure to isolate a leaf-color mutant in Cymbidium, which showed a notable decreases in Chl and carotenoid levels. The RNA-Seq analysis identified 2,267 DEGs, including 724 up-regulated and 529 down-regulated in the S12 mutant. A functional classification of these genes allowed us to identify the chlorophyllase gene, CLH2, which was induced~18-fold in the S12 mutant. The CLH2-encoded enzyme facilitates Chl catabolism, and increased levels of CLH2 have been associated with reduced Chl contents in a number of plants [49][50][51][52][53][54][55][56]. These results are consistent with a previous study [57], which suggested that the yellow-striped leaves of C. sinense variants were associated with increased Chl degradation. Notably, C. sinense variants characterized by Zhu et al. [57] showed~1-3-fold higher expression levels of CLH2 compared with the~18-fold increase seen in the S12 mutant. Zhu et al. also observed a~1-3-fold increase in the expression levels of RCCR [57], which acts downstream of CLH2 in the Chl degradation pathway [57]. However, unlike the C. sinense variants characterized by Zhu et al. [57], the S12 mutant showed normal expression levels of RCCR. Consistent with previous studies [18,[23][24][25][26]57], the S12 mutant also showed reduced levels of carotenoids, a phenotype commonly observed in Chl-deficient mutants. The S12 mutant also contained an altered Chl a to Chl b ratio, in which Chl b is required to stabilize the light-harvesting protein complex [58]. This, in turn, correlated with a reduced expression level of NYC1 in the S12 mutant, which encodes a chlorophyll b reductase that catalyzes the conversion of Chl b to Chl a. Thus, a reduced Chl content and decreased Chl a/b ratio in the S12 mutant may be associated with fewer light-harvesting antenna complexes. In contrast, the S12 mutant showed normal Transcriptome analysis to identify genes associated with the yellow-leaf phenotype in Cymbidium mutant expression levels of all 16 genes associated with Chl biosynthesis. This further suggests that the reduced Chl content in S12 leaves was likely associated with Chl catabolism and not biosynthesis. Metals such as iron (Fe), copper (Cu), manganese (Mn), zinc (Zn), and magnesium (Mg) play key roles as cofactors in photosynthesis: Fe, as a cofactor for three photosynthetic electron transfer chain complexes; Cu, as a cofactor for thylakoid lumen electron transport protein plastocyanin; Mn, for photosystem II functions; Zn, plays a key role in the catalysis of chloroplastic ß-carbonic anhydrase enzyme; and Mg, in the center of the Chl ring [29,45,59]. The transitions of these metal ions are regulated to maintain cellular homeostasis, and excess metal ions cause oxidative stress because of their deleterious interactions with oxygen [29,59]. In this study, we identified diverse metal ion transporters among down-regulated DEGs. In particular, the down-regulated expression levels of CCH (Cu-ion transporter) and YSL9 (Fe-ion transporter) were confirmed by qRT-PCR (Table 5). The Arabidopsis CCH gene is a functional homolog of the yeast (Saccharomyces cerevisiae) gene Anti-oxidant 1, which, when mutated, results in a reduced Fe-uptake capability [60]. Senoura et al. [61] reported that OsYSL9 mainly localizes in the plasma membrane and transports Fe(II)-nicotianamine and Fe(III)-deoxymugineic acid. The expression of OsYSL9 is repressed in the leaves under Fe-starvation conditions. The down-regulated expressions of CCH and YSL9 could affect the Fe ion transition in the cell, which could be correlated with the reduced Chl content in the S12 mutant. NRAMP proteins are involved in Fe homeostasis [33], and the expression levels of AtNRAMP1, 3, and 4 are up-regulated in response to Fe deficiency [62,63]. Additionally, AtNRAMP3 and 4 remobilize vacuolar Mn in leaves and have important roles in photosystem II functions [64]. In cyanobacteria, the ABC-type Mn transport complex is induced under Mn-starvation conditions [65]. In the S12 mutant, the reduced Fe content is expected because of the down-regulated expression levels of CCH and YSL9 (Table 5). Additionally, the up-regulated expression levels of NRAM2 and the ABC transporter family were identified (Table 4), which was consistent with previous reports [62][63][64][65].
In the present study, we found that seven genes involved in ion transport, including two metal ion transporters (CCH and YSL9), were down-regulated, and CLH2, associated with Chl degradation, was up-regulated in the yellow leaf-color mutant, S12. This provides useful information for understanding Chl biosynthesis and degradation in Cymbidium. In addition, these results show that γ-ray-based mutagenesis can be employed as a useful tool to generate genetic diversity among orchid species.
Supporting information S1