Transcriptome-based gene profiling provides novel insights into the characteristics of radish root response to Cr stress with next-generation sequencing

Radish (Raphanus sativus L.) is an important worldwide root vegetable crop with high nutrient values and is adversely affected by non-essential heavy metals including chromium (Cr). Little is known about the molecular mechanism underlying Cr stress response in radish. In this study, RNA-Seq technique was employed to identify differentially expressed genes (DEGs) under Cr stress. Based on de novo transcriptome assembly, there were 30,676 unigenes representing 60,881 transcripts isolated from radish root under Cr stress. Differential gene analysis revealed that 2985 uingenes were significantly differentially expressed between Cr-free (CK) and Cr-treated (Cr600) libraries, among which 1424 were up-regulated and 1561 down-regulated. Gene ontology (GO) analysis revealed that these DEGs were mainly involved in primary metabolic process, response to abiotic stimulus, cellular metabolic process and small molecule metabolic process. Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis showed that the DEGs were mainly involved in protein processing in endoplasmic reticulum, starch and sucrose metabolism, amino acid metabolism, glutathione metabolism, drug and xenobiotics by cytochrome P450 metabolism. RT-qPCR analysis showed that the expression patterns of 12 randomly selected DEGs were highly accordant with the results from RNA-seq. Furthermore, many candidate genes including signaling protein kinases, transcription factors and metal transporters, chelate compound biosynthesis and antioxidant system, were involved in defense and detoxification mechanisms of Cr stress response regulatory networks. These results would provide novel insight into molecular mechanism underlying plant responsiveness to Cr stress and facilitate further genetic manipulation on Cr uptake and accumulation in radish.

Transcriptome-based gene profiling provides novel insights into the characteristics of radish root response to Cr stress with next-generation sequencing Introduction Heavy metal (HM) contamination has become one of the major worldwide environmental problems affecting soil, water and human health, leading to a considerable reduction in yield and food safety (Chandra and Kulshreshtha, 2004). An elevated content of HMs in human diet constitutes a potential health risk through food chain, vegetable crops, representing an important pathway for the movement of potentially toxic metals from soil to human beings (Gupta et al., 2012). Chromium (Cr), the seventh most abundant metal on earth, is one of the most toxic HMs to living organisms and ecosystems (Katz and Salem, 1994). Cr resource is mainly from industrial process products such as leather tanning, electroplating, steel production, metal finishing, catalyst application, pigment manufacturing and metal corrosion inhibitors (Shanker et al., 2009). Trivalent Cr (III) and hexavalent Cr (VI) species are the stable forms of Cr exist in the earth's crust. Accumulation of Cr in plants can interfere with several metabolic processes such as photosynthesis, water relation and uptake of nutrients. These stresses could reduce root growth and phytomass, induce chlorosis in plants, wilting and plasmolysis in root cells and finally plant death (Hayat et al., 2011). Being a strong oxidizer, Cr (VI) at even low concentration is highly toxic and more mobile in soil/water systems, and is also known to cause allergies, irritations, and even respiratory track disorders to human beings (Shanker et al., 2005).
In recent years, the next-generation sequencing (NGS) platform including Illumina/Solexa sequencing technology has provided a powerful tool for RNA sequencing, transcriptome assembly and gene expression profiling (Mochida and Shinozaki, 2011). NGS-based RNA-Seq is a very powerful technology for transcriptomics studies, and enables us to investigate the gene activities of organisms in diverse tissues and different stages under various conditions (Geng et al., 2011). Based on genomewide expression profiles by high throughput sequencing, digital gene expression (DGE) could be employed to characterize, identify and quantify rare gene transcripts on the global transcriptome level without prior sequence knowledge (Strickler et al., 2012). In recent years, RNA-seq has provided a useful tool for identification of related genes and their expression patterns in crop pests and plant species responding to various biotic and abiotic stresses (Nachappa et al., 2012;Yu et al., 2012). Highthroughput Solexa/Illumina sequencing was successfully applied to characterize gene expression of soybean under salt, salinealkali and drought stress (Fan et al., 2013). In order to investigate head smut disease-resistance mechanisms, DGE analysis was used to identify the transcriptional changes in the roots of maize (Zhang et al., 2013). Many reports showed that the DGE approach has provided valuable tools for qualitative and quantitative gene expression analysis (Yamaguchi et al., 2010;Li et al., 2014).
Cr is absorbed by plant mainly through the roots. With transcriptomic and metabolomic techniques, some vital genes encoding PDR-like ABC transporter, multidrug resistance protein 4 and glutathione S-transferase GSTU6 were identified to be involved in response to Cr stress in rice (Dubey et al., 2010). Even though the molecular mechanisms of Cr stress effect on plants have been under scrutiny for decades, critical genes involved in response to Cr stress were not well identified and characterized in root vegetable crops. Radish (Raphanus sativus L., 2n = 2x = 18), belonging to the Brassicaceae family, is a major annual or biennial worldwide root vegetable crop especially in East Asia. However, most reports on Cr in plants have been concerned with its effects on plant growth, uptake, toxicity, translocation, and soil-plant relationships (Panda and Choudhury, 2005;Shanker et al., 2005). A comprehensive investigation on the molecular mechanism underlying Cr absorption and transport is urgently required in radish.
So far, there is no report on systematic identification and characterization of critical genes involved in Cr-responsiveness, uptake, transportation and accumulation in radish. Moreover, it is important to identify the key genes and clarify the molecular mechanism of Cr stress-response for altering the accumulation of Cr in radish plant through genetic improvements. In this study, RNA-Seq approach was employed with genome-wide transcriptional profiling of the without (control) and 600 mg/L K 2 Cr 2 O 7 treatment to identify differentially expressed genes (DEGs) under Cr stress in radish roots. The aim of this study was to identify key DEGs under Cr stress in radish roots. With the highthroughput sequencing technique applied, a large number of unigenes, Cr responsive DEGs and their expression patterns were successfully obtained, and the expression profiling of a proportion of differentially regulated genes were validated by RT-qPCR. Furthermore, based on DEGs enrichment in the corresponding pathway, a hypothetical model of Cr stress-response regulatory network in radish was proposed. This study represents a first comprehensive transcriptome-based characterization of Cr responsive DEGs in radish roots, these results would provide fundamental insights into the complex Cr-responsive gene regulatory networks, and facilitate further studies on molecular genetic mechanisms underlying plant responses to Cr stress in root vegetable crops.

Plant Material
A high-Cr-accumulation radish advanced inbred line, "NAU-RG06, " was used in the study. Germinated seeds were sown in plastic pots with matrix and soil at a ratio of 1:1 and cultivated in a growth chamber at 25 • C day/18 • C night with a 14 h light/10 h dark photoperiod. Seedlings of 30 days old were transferred into liquid medium in a plastic box and grown for 3 days. Cr treatment was carried out under the same growth conditions by soaking the root in a K 2 Cr 2 O 7 (Cr 6+ ) solution of 0 (CK) and 600 mg L −1 (Cr600) for 72 h. After treatments, at least three replicates of roots for each treatment were harvested and immediately frozen in liquid nitrogen for further use.

RNA Isolation and Illumina Sequencing
Total RNA was extracted from root samples using TRIzol reagents (Tiangen Biotech Co., Ltd., China) according to the manufacturer's instructions. The RNAs were treated with RNasefree DNase I to eliminate contaminated genomic DNA. Two radish cDNA libraries, CK and Cr600, were constructed from control and 600 mg L −1 K 2 Cr 2 O 7 treated root samples using the Illumina Paired End Sample Prep Kit. Briefly, poly (A) mRNA was enriched from total RNA using Sera-mag Magnetic Oligo (dT) Beads (Thermo Fisher Scientific, USA) and then mRNAenriched RNAs were chemically fragmented to short pieces using the fragmentation solution (Ambion, USA). Double-stranded cDNA was generated using the Superscript Double-Stranded cDNA Synthesis Kit (Invitrogen, USA). After that, the Illumina Paired End Sample Prep kit was used for RNA-seq library construction and was then sequenced using Illumina HiSeq ™ 2000.

De Novo Transcriptome Assembly and Annotation
Illumina pipeline was used for filtering the raw sequence reads. The 3 ′ adaptor sequence was removed from raw sequences. All low-quality tags, such as short tags (<25 nt), empty tags, and tags with only one copy number were removed. De novo transcriptome assembly was accomplished from all the clean reads with the Trinity program and unigenes were generated (Grabherr et al., 2011;Wang et al., 2013). Only sequences with perfect homology or not more than two nucleotide mismatches were considered for conservative and accurate annotation. The assembled transcripts were further used as query sequences to search against NCBI non-redundant (nr) protein and COG (Clusters of orthologous groups of proteins) databases with BLASTX. For nr annotation, the Blast2GO program was used to get GO (Gene Ontology) annotation of assembled transcripts (Conesa et al., 2005). After getting GO annotation for assembled transcripts, WEGO software was employed to fulfill GO functional classification (Ye et al., 2006).

Identification of Differentially Expressed Genes
The clean reads were mapped to the reference database which included whole radish ESTs and unigenes from NCBI and transcriptome sequences using SOAPaligner/soap2 (Li et al., 2009). RPKM (reads per kilobase of exon model per million mapped reads) were used to evaluate the gene expression value and quantify transcript levels (Mortazavi et al., 2008). Using the DEGseq program, significantly differential gene expression was identified between the CK and Cr600 libraries (Wang et al., 2010). The false discovery rate (FDR) was used to determine the threshold p value in multiple tests procedure described by Benjamini and Yekutieli (2001). In this study, a stringent of FDR ≤ 0.001 and log 2 | FC (ratio of stress/control) | ≥2 was used as the threshold to judge the significant difference of gene expression.

GO and KEGG Pathway Enrichment Analysis of DEGs
Functional classes were assigned according to GO mapping provided by the ensemble database. KEGG (Kyoto encyclopedia of genes and genomes) pathway analysis was based on the comparative results between our mapping genes and the current KEGG database (Kanehisa et al., 2008). The differentially expressed genes (DEGs) were used for GO and pathway enrichment analysis, and a corrected p ≤ 0.05 was selected as a threshold level of significance to determine enrichment in the gene sets .

RT-qPCR Analysis
Quantitative RT-PCR analysis was used to validate the expression of the candidate genes. RT-qPCR was conducted on MyiQReal-Time PCR Detection System (Bio-Rad) using a SYBR Primix Ex Taq kit (TaKaRa, Dalian, China) according to the manufacturer's instructions. Gene-specific primers and a β-actin-specific primer pair were designed according to the gene sequences using Beacon Designer 7.0 software. Each reaction was prepared in a total volume of 20 µl containing 10 µl SYBR Green mix, 2 µl diluted cDNA and 0.2 µM of each primer. Amplification was carried out with the following cycling parameters: heating for 10 s at 95 • C, 40 cycles of denaturation at 95 • C for 5 s, annealing at 58 • C for 30 s and extension at 72 • C for 10 s (Xu et al., 2012). Each sample was analyzed in triplicates and the expression values were normalized against β− actin. The molecular weight of the products was confirmed via diagnostic agarose gel and the melting curves were analyzed. Analysis of the relative gene expression data was conducted using the 2 − C T method (Livak and Schmittgen, 2001).

Radish Root Transcriptome Assembly and Transcript Function Classification
A total of 37.95M and 27.45M raw reads were generated from CK and Cr600 library with Illumina Solexa sequencing technology, respectively. Before mapping these sequencing reads to the reference genome, the reads with connectors, filtered low-length reads and low-quality reads were removed. A total of 26.38M (69.52%) and 15.03M (54.74%) 101-nt clean reads were generated in cDNA libraries of CK and Cr600, respectively ( Table 1). The sequencing reads were aligned against the combined radish reference database, and a total of 60,881 isogenes or 30,676 unigenes with 71,526,014 nt were filtered out for the differential expression analysis. Moreover, the average length of the isogene was 1,174.85 bp, and the largest and shortest isogene was 15,585 and 306 bp, respectively. The sequence length of these assembled transcripts was mainly distributed between 306 and 3000 bp with 401 and 600 bp length in the highest abundance (Figure 1).

Identification and Functional Analysis of DEGs Under Cr Stress
By comparing two Solexa libraries (CK and Cr600), a great number of differentially expressed reads were identified. To study a subset of genes that were associated with Cr stress, we analyzed the most differentially regulated tags with a log 2 ratio ≥ 2 or ≤ −2, FDR ≥ 0.001. Using these standards, 6295 isogenes (2985 unigenes) were identified to be differentially expressed. Among all the differentially expressed genes (DEGs), 2842 isogenes (1424 unigenes) were identified to be up-regulated and the rest were down-regulated (Figure 4).
GO analysis was performed to predict functions by mapping each DEGs into the records of the GO database. According to Bonferroni-corrected P ≤ 0.05, a total of 214 terms for the upregulated transcripts were enriched, and 269 terms were downregulated. In the GO cellular component class, there were more up-regulated genes annotated for cell part (GO: 0044464), intracellular part (GO: 0044424) and cytoplasmic part (GO: 0044444). While the down-regulated genes were mainly annotated with membrane (GO: 0016020), plasma membrane (GO: 0005886) and intrinsic to membrane (GO: 0031224) (Tables S2, S3). In the GO molecular function class, many up-regulated genes were annotated by catalytic activity (GO: 0003824), oxidoreductase activity (GO: 0016491) and transporter activity (GO: 0005215). The down-regulated genes were mainly annotated by binding (GO: 0005488) and ATP binding (GO: 0005524) (Tables S2, S3). In the GO biological process class, there was a higher number of DEGs annotated by response to abiotic stimulus (GO: 0009628), primary metabolic process (GO: 0044238), cellular metabolic process (GO: 0044237) and small molecule metabolic process (GO: 0044281) (Tables S2, S3).
With Pathway-based analysis specific biological functions could be assigned to genes (Figure 5). A total of 3259 DEGs were annotated by KEGG pathway analysis, and a total of 11 and 15 significantly enriched pathways were found for up-and down-regulated DEGs, respectively (Tables 2, 3). For the upregulated isogenes, the primarily enriched pathways included protein processing in endoplasmic reticulum [ko04141], glutathione metabolism [ko00480], drug metabolism -cytochrome P450 [ko00982] and xenobiotics metabolism by cytochrome P450 [ko00980] (Table 2; Figure 5). Down-regulated genes were mainly involved in starch and sucrose metabolism [ko00500], amino sugar and nucleotide sugar metabolism [ko00520] and pentose and glucuronate interconversions [ko00040] ( Table 3). From the result of pathway enrichment analysis, it could be suggested that some genes might interact with each other to play specific roles in certain biological processes (Figure 5).

Characterization of Cr Stress-Responsive Genes
Signal sensing and transduction were important participants in response to heavy metal stress. Once entered the plant cell, Cr induces reactive oxygen species (ROS) which were responsible for the activation of MAPK (Mitogen-activated protein kinases) kinase cascade. MAPKs including MAPKKK, MAPKK, and MAPK, were involved in response to various environmental, hormonal and developmental stimuli. Ca/CaM is the signal transduction cascade system, which is associated with to heavy metal stress. In this study, a total of eight DEGs had been identified to be highly homologous to  genes encoding MAPKs (MAPKKK18, MAPKKK19, MAPK3, MAPK16, MAPK17, MAPK18, MAPK19, and MAPK20). In total, 10 DEGs including CML42, CML44, and CML50, were calciumdependent protein kinase and calcium-binding protein genes, which were up-regulated under Cr stress in radish roots ( Table 4,  Table S4).
High Cr concentration would induce ROS and result in oxidative stress. Accordingly, plants would produce antioxidant enzymes and non-enzymic antioxidants to alleviate ROS damages. In the current, the transcripts encoding peroxidase34, dehydroascorbate reductase (DHAR) and glutathione reductase (GR) were identified to be up-regulated, while peroxidase17, 21, 49, and 64, superoxide dismutase (SOD) and gamma-tocopherol methyltransferase genes were down-regulated in response to Cr stress ( Table 4, Table S4).
Metal transporters and chelating compounds play significant roles in coping with HM stress by exclusion, chelation and compartmentalization. A total of 69 DEGs were identified as members of different metal transporter families, which were mainly related to sulfate transporters, heavy metal transport/detoxification domain-containing protein, ATP binding cassette (ABC), ZRT, IRT-like proteins (ZIPs), Natural resistance-associated macrophage proteins (Nramps) families, cation diffusion facilitators (CDFs) and MATE efflux family protein ( Table 4, Table S4). Almost all the genes encoding MATE efflux family proteins identified in this study were up-regulated under Cr stress. Four DEGs including three upregulated and one down-regulated, were found to be homologous with genes encoding metallothioneins (MTs) under Cr stress. In addition, one phytochelatin (PC) related gene (comp25826_c0) was up-regulated in radish roots in the present study ( Table 4,  Table S4).
Furthermore, Heat shock proteins (HSPs) and Cytochrome P450 (CYP450) genes were also identified in radish roots response to Cr stress. It was reported that HSPs were induced by various kinds of stresses, such as high temperature, heavy metal and drought stress (Dubey et al., 2010). A total of 10 HSP genes are up-regulated in response to Cr stress including two HSP20like chaperone genes, four heat shock transcription factor (HSF) genes, and two DNAJ HSP-like protein genes, HSP21 and HSP22. Of the four down-regulated HSPs three were in the HSP-like family and one HSP23. CYP450 genes were involved in a group of secondary metabolite production that helped the plants to resist against various stressful conditions (Balusamy et al., 2013). For example, CYP79F1, CYP83A1, and CYP79B3 involved in glucosinolate biosynthesis were differentially expressed in radish roots under Cr stress. A total of 24 CYP450 genes including 17 up-regulated and seven down-regulated genes were differentially expressed under Cr stress ( Table 4, Table S4).

RT-qPCR Validation
To evaluate the validity of Illumina sequencing and to further confirm the patterns of differential gene expression, a subset of 12 genes highly expressed in the Cr treatment were selected and detected by RT-qPCR analysis with gene-specific primers. These 12 genes were mainly related to defense and detoxification mechanisms including signaling protein kinases (MPK19), TFs (WRKY33, Bzip44), metal transporters (ABCA3, ZFP9, ZFP4), chelate compound biosynthesis (GSTU19, GST1, PCS1), antioxidant system (PX49, APX5) and HSP (Hsc70-2). The expression patterns from RT-qPCR indicated a general agreement with those from the Solexa sequencing results (Figure 6; Table S5).

Discussion
Chromium (Cr) is considered as a serious environmental pollutant posing a critical concern to human health through food chain (Zayed et al., 1998;Shanker et al., 2005). Exploration of the molecular mechanisms underlying plant response to heavy metal (HM) stress, particularly those facilitate effectively controlling HM uptake, translocation and accumulation, has been a focus work for plant genetic manipulation and genomic research .
Radish is an economically important vegetable crop with an edible taproot. The mechanism of radish response to Cr stress has not been clarified at the molecular level before. A global analysis of transcriptome could facilitate the identification of critical gene expression and regulatory mechanisms in plant response to biotic and abiotic stresses such as cold (Shan et al., 2013), salinity (Dang et al., 2013) and HM stress (Yamaguchi et al., 2010;Yu et al., 2012;Li et al., 2014). In the present study, with a transcriptome profiling of radish root, some critical DEGs responsive to Cr stress were discovered, and the hypothesis model of Cr stressresponsive gene regulatory network in radish roots was proposed (Figure 7). FIGURE 5 | The key candidate genes and their interaction with major pathways under Cr stress of radish roots. Red and green boxes represent up-and down-regulated candidate genes, respectively. Blue and gray boxes represent gene definition and major pathways, respectively.

The Role of Antioxidant System in Plant
Response to Cr Stress Plants exposed to toxic heavy metal (HM) can generate reactive oxygen species (ROS), resulting oxidative stress (Panda and Choudhury, 2005). Oxidative stress can lead to lipid peroxidation in plants and disrupt the cellular functions, resulting in inhibition of plant growth and development (Yadav, 2010). To fight against oxidative stress, plants have evolved a complex antioxidant system for elimination of ROS to protect the cell from damage. Actually, the defense system including enzymatic and non-enzymatic antioxidants was considered as an important mechanism of metal detoxification and tolerance in plants. SOD is the first line of defense against ROS-mediated damage and could scavenge toxic O − 2 in plant cells (Zeng et al., 2012). POD, CAT, and the ascorbate-glutathione cycle (APX and GR) play crucial roles in scavenging H 2 O 2 . Cr has the redox behavior which exceeds that of other metals including Co, Fe, Zn, and Ni (Panda and Choudhury, 2005). In the current study, some functional genes encoding peroxidase17, 21, 34, 49, and 64, SOD, DHAR, tocopherol and GR, were identified in the radish root exposed to Cr ( Table 4, Table S4). It was interesting that almost all the DEGs related to peroxidase were down-regulated under Cr stress, which might result from the cellular injury including bio-membrane lipid peroxidation under Cr stress.
Glutathione is also an important antioxidant involved in cellular defense against toxicants. GSH is a redox buffer protecting the cytosol and other parts of cells against ROS which are induced by biotic and abiotic stresses (Zitka et al., 2013). In this study, one cysteine synthase (comp25459_c0), two glutathione synthetase (GSS) and 17 out of 19 glutathione S-transferase (GSTs), were up-regulated in radish roots following Cr stress ( Table 4, Table  S4). These findings provided further evidence for the role of GSH in antioxidant metabolism. Furthermore, GSH, serving as a precursor, could participate in the synthesis of phytochelatins (PCs).  Therefore, it could be concluded that antioxidant system plays a fundamental role in plant cellular detoxification.

Signal Transduction and TFs in Plant Response to Cr Stress
When encountered with extra-cellular stimuli or to prevent HM entrance from susceptible sites into the protoplast, the cell walls could activate a variety of specific stress-responsive signaling proteins, such as MAPKs and calcium-binding related proteins. It was reported that MAPKs were important mediators in signal transmission, connecting the perception of external stimuli to cellular responses (Sanchita et al., 2014). The calmodulin (CaM) system could regulate HM ion transport, gene expression, and stress tolerance (DalCorso et al., 2010). In this study, MAPK18, MAPK20, and other calcium-binding related protein genes were expressed differentially in the radish root under Cr exposure. These genes were also found to be expressed differentially in radish under Pb stress, indicating that MAPKs and calciumbinding related proteins were involved in signal transduction in plants under HM stress .
Once phosphorylated, MAPKs enter into the nucleus through the phosphorylation of TFs to regulate gene expression. A series of TFs including bZIP, WRKY, ERF, and MYB have been shown to play important roles in regulating the expression of specific stress-related genes (Thapa et al., 2012). In this study, different kinds of TFs including up-regulated WRKY, ERF, MYB and bZIP60 were identified in radish response to Cr stress ( Table 4,  Table S4). It could be concluded that the plants would activate a series of TFs to regulate corresponding transcriptional processes to alleviate the phytotoxicity of HMs.

Transporter and Chelate Compound in Plant Response to Cr Stress
It was reported that a large number of transporter families including ABC, Nramps, ZIPs, CDFs, and MATE efflux family proteins involved in HM uptake, transport, distribution and plant tolerance to HM stress (Dubey et al., 2010;Lan et al., 2013). A large number of transporters including sulfate transpoters, MATEefflux family proteins, ABC transporter family proteins, heavy metal-associated domain containing proteins, were differentially expressed after challenge with Cr in the rice root (Dubey et al., 2010). Wang et al. (2013) reported that in radish response to Pb stress, the main differentially expressed metal transpoters were ABC, Nramps, ZIPs, and CDFs, which were also identified in the present study ( Table 4, Table S4). These results indicate that there would be similar mechanisms underlying the uptake, transport and detoxification of HMs among different plant species. It was found that PCs and MTs could chelate metal ions to form metal-chelate compounds and sequester them in the vacuole (Sharma and Chakraverty, 2013). MTs are low-molecularweight cysteine-rich metal binding peptides. Under HM stresses, different MT types have been reported in plants including oil palm (Abdullah et al., 2002), poplar (Kohler et al., 2004), and radish . In the present study, the genes encoding MT3 and MT-like proteins were also found under Cr stress in the radish root. Moreover, two phytochelatins (PCs) related genes were also identified ( Table 4, Table S4). It was reported that PCs have a high affinity in binding HMs, and PC-HM complexes are then transported into the vacuole, leading to inactivation of HMs in the plant cells (Zeng et al., 2012).

CYP450s and HSPs May Be Involved in Plant Response to Cr Stress
CYP450s play crucial roles in biosynthesis of a variety of endogenous lipophilic compounds such as fatty acids, sterols, phenylpropanoids, terpenoids, and phytoalexins brassinolides, which could enhance the plant tolerance under stresses (Narusaka et al., 2004;Vázquez et al., 2013;Fariduddin et al., 2014). It was reported that CYP450s were induced by abiotic stress, such as As (Chakrabarty et al., 2009) in rice, Pb (Liu et al., 2009) in Arabidopsis, Cr (Dubey et al., 2010) in rice and dehydration stress in cabbage (Yu et al., 2012), and biotic stress including insecticide stress in ladybirds  and bacterial infection in tobacco (Daurelio et al., 2011). The present results showed that a total of 24 CYP450 genes were differentially expressed in the radish root under Cr stress. It was interesting that biotic and heavy metal FIGURE 6 | Validation of the RNA-Seq expression profiles of 12 genes randomly selected from DEGs by RT-qPCR.
FIGURE 7 | A regulatory model predicted from the genes expressed differentially during Cr stress in radish roots. A green box represents a gene enrichment site and corresponding gene or protein name. A blue box represents a gene enrichment site and corresponding gene category. stress response in plants may have a similar mechanism (Rai and Mehrotra, 2008).
Furthermore, HSPs, which are molecular chaperones, have the function of repairing proteins under stress conditions. Hsp90-1 was found to be highly accumulated in tomato under Cr and As stress (Goupila et al., 2009). Three up-regulated and five down-regulated HSP genes were identified under Cr stress in rice (Dubey et al., 2010). In this study, a similar result was obtained and a total of 14 HSP genes were identified in the radish root under Cr stress. It could be inferred that HSPs play vital roles in enhancing the plant tolerance to Cr stress.

Cr Stress-Response Regulatory Networks in Radish
Plants have developed a variety of complex physiological and genetic mechanisms to cope with HM toxicity. In general, there are several strategies to defense and regulate HM stress including immobilization, exclusion, chelation and compartmentalization (Hall, 2002;Sharma and Chakraverty, 2013). Based on the DEGs enrichments with the technique, Cr stress-response regulatory network was put forward in radish roots (Figure 7). As shown in Figure 7, when a plant cell was exposed to Cr, oxidative stress would be induced. The oxidative stress could generate lipid peroxidation and cause severe damages to biomembranes in plant cells (Panda and Choudhury, 2005). In order to decrease excess levels of oxidative damage, several metabolites and enzymes including GPX, SOD, APX, and CAT were activated. After sensing the signal, the plant cells activated some Cr stress responsive signaling molecules such as MAPKs and CMLs, which then guided TFs (WRKY, bZIP, ERF, MYB families) in the nucleus to regulate the expression of some functional genes. These genes may encode transporters (ABC, ZIP, Nramps, and CDFs) and chelating compounds (MTs and PCs), which would alleviate Cr toxicity via diverse mechanisms in plant cells (Figure 7; Table 4; Table S4).

Conclusion
The comprehensive transcriptome-based characterization of Cr responsive differentially expressed genes (DEGs) was firstly conducted in radish roots. This result revealed that the expression of 1424 unigenes was up-regulated, and 1561 unigenes were down-regulated in radish roots under Cr stress. Most of the DEGs were relative to antioxidant system, signal transduction and TFs, transporters and chelate compounds biosynthesis. In addition, a model of Cr stress-response regulatory network in the radish root was proposed. These results would provide useful information for dissecting the molecular genetic mechanism underlying Cr uptake, sequestration, translocation and detoxification in root vegetable crops.