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.

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 genome-wide 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). High-throughput Solexa/Illumina sequencing was successfully applied to characterize gene expression of soybean under salt, saline-alkali 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₂Cr₂O₇ 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 high-throughput 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.

Materials and Methods

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₂Cr₂O₇ (Cr⁶⁺) solution of 0 (CK) and 600 mg L⁻¹ (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 RNase-free DNase I to eliminate contaminated genomic DNA. Two radish cDNA libraries, CK and Cr600, were constructed from control and 600 mg L⁻¹ K₂Cr₂O₇ 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 mRNA-enriched 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₂ | 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 (Wang et al., 2013).

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

Results

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

FIGURE 1

Figure 1. The length distribution of the assembled transcripts.

TABLE 1

Table 1. Reads from RNA-Seq library sequencing.

Gene ontology (GO) analysis indicated all these genes were distributed in three categories, namely biological process, cellular component and molecular function. Totally 42,178 transcript sequences (69.28%) were assigned to 57 GO terms at the second level. Cellular process (27,602 sequences, 65.44%), metabolic process (24,851 sequences, 58.92%) and response to stimulus (13,134 sequences, 31.14%) accounted for the main part in biological process category. In relation to cellular component, there were sequences associated with cell part (31,102 sequences, 73.74%) and organelle part (12,247 sequences, 29.04%) which represented the most abundant categories. The sequences associated with binding (26,066 sequences, 61.80%) and catalytic activity (19,614 sequences, 46.50%) were dominant in the molecular function category (Figure 2).

FIGURE 2

Figure 2. Gene ontology classification of assembled transcripts.

COG analysis showed that a total of 15,590 transcript sequences with high homology were grouped into 23 functional categories (Figure 3). The first three largest categories were “general function prediction only” (3200, 20.5%), “transcription” (1842, 11.8%), “replication, recombination and repair” (1551, 9.95%), followed by “signal transduction mechanisms” (1521, 9.8%) and “post translational modification, protein turnover, chaperones” (1066, 6.8%) (Figure 3; Table S1).

FIGURE 3

Figure 3. COG function classification of unigenes.

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

FIGURE 4

Figure 4. Volcano plot of gene expression difference between Cr600 and CK libraries.

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 up-regulated transcripts were enriched, and 269 terms were down-regulated. 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 up-regulated 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).

FIGURE 5

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.

TABLE 2

Table 2. The significantly enriched pathways of up-regulated DEGs.

TABLE 3

Table 3. The significantly enriched pathways of down-regulated DEGs.

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 calcium-dependent protein kinase and calcium-binding protein genes, which were up-regulated under Cr stress in radish roots (Table 4, Table S4).

TABLE 4

Table 4. Some critical DEGs responsive to Cr stress in radish roots.

There were 94 DEGs that were of high homology with different kinds of transcription factors (TFs), such as the WRKY family (i.e., WRKY1, 6, 13, 26, 28, 33, 46, 48, 50, 70, 71, and 75), ethylene-responsive factor (ERF) family (i.e., ERF008, 012, 018, 070, 071, 088, and 098), myeloblastosis protein (MYB) family (i.e., MYB1, 4, 28, 29, 44, 46, 51, 58, 73, 95, and 108), basic leucine zipper (bZIP) family and dehydration-responsive element-binding protein-type transcription factors (DREB) in this study (Table S4). Among these DEGs, WRKY6, 26, 28, 33, 50, 71, and 75, bZIP60, ERF070, 071, 088, and 098, MYB1, 51, 108 were found to be 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 up-regulated 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 HSP20-like 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).

FIGURE 6

Figure 6. Validation of the RNA-Seq expression profiles of 12 genes randomly selected from DEGs by RT-qPCR.

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 (Wang et al., 2013).

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 stress-responsive gene regulatory network in radish roots was proposed (Figure 7).

FIGURE 7

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.

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⁻₂ in plant cells (Zeng et al., 2012). POD, CAT, and the ascorbate–glutathione cycle (APX and GR) play crucial roles in scavenging H₂O₂. 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 Dhawan and Sharma, 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 calcium-binding related proteins were involved in signal transduction in plants under HM stress (Wang et al., 2013).

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, MATE-efflux 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-molecular-weight 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 (Wang et al., 2013). 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 (Zhang et al., 2012) 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 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 bio-membranes 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.

Accession Code

The RNA SEQ raw data have been deposited in NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra) with accession numbers: SRX256970 (CK) and SRX862647 (Cr600).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was in part supported by grants from the National Natural Science Foundation (NSF) of China (31372064, 31171956), NSF of Jiangsu Province (BK20140706), Key Technology R and D Program of Jiangsu Province (BE2013429), JASTIF [CX (12)2006, CX(13)2007] and the PAPD.

Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fpls.2015.00202/abstract

References

Abdullah, S. N. A., Cheah, S. C., and Murphy, D. J. (2002). Isolation and characterisation of two divergent type 3 metallothioneins from oil palm, Elaeis guineensis. Plant Physiol. Biochem. 40, 255–263. doi: 10.1016/S0981-9428(02)01366-9

CrossRef Full Text | Google Scholar

Balusamy, S. R., Kim, Y. J., Rahimi, S., Senthil, K. S., Lee, O. R., Lee, S., et al. (2013). Transcript pattern of cytochrome P450, antioxidant and ginsenoside biosynthetic pathway genes under heavy metal stress in panax ginseng meyer. Bull. Environ. Contam. Toxicol. 90, 194–202. doi: 10.1007/s00128-012-0891-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple tesing under dependency. Ann. Stat. 29, 1165–1188. doi: 10.1214/aos/1013699998

CrossRef Full Text | Google Scholar

Chakrabarty, D., Trivedi, P. K., Misra, P., Tiwari, M., Shri, M., Shukla, D., et al. (2009). Comparative transcriptome analysis of arsenate and arsenite stresses in rice seedlings. Chemosphere 74, 688–702. doi: 10.1016/j.chemosphere.2008.09.082

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Chandra, P., and Kulshreshtha, K. (2004). Chromium accumulation and toxicity in vascular plants. Bot. Rev. 70, 313–327. doi: 10.1663/0006-8101(2004)070[0313:CAATIA]2.0.CO;2

CrossRef Full Text | Google Scholar

Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

DalCorso, G., Farinati, S., and Furini, A. (2010). Regulatory networks of cadmium stress in plants. Plant Signal. Behav. 5, 663–667. doi: 10.4161/psb.5.6.11425

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Dang, Z. H., Zheng, L. L., Wang, J., Gao, Z., Wu, S. B., Qi, Z., et al. (2013). Transcriptomic profiling of the salt-stress response in the wild recretohalophyte Reaumuria trigyna. BMC Genomics 14:29. doi: 10.1186/1471-2164-14-29

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Daurelio, L. D., Petrocelli, S., Blanco, F., Holuigueb, L., Ottadoa, J., and Orellano, E. G. (2011). Transcriptome analysis reveals novel genes involved in nonhost response to bacterial infection in tobacco. Plant Physiol. 168, 382–391. doi: 10.1016/j.jplph.2010.07.014

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Dubey, S., Misra, P., Dwivedi, S., Chatterjee, S., Bag, S. K., Mantri, S., et al. (2010). Transcriptomic and metabolomic shifts in rice roots in response to Cr (VI) stress. BMC Genomics 11:648. doi: 10.1186/1471-2164-11-648

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Fan, X. D., Wang, J. Q., Yang, N., Dong, Y. Y., Liu, L., Wang, F. W., et al. (2013). Gene expression profiling of soybean leaves and roots under salt, saline-alkali and drought stress by high-throughput Illumina sequencing. Gene 512, 392–402. doi: 10.1016/j.gene.2012.09.100

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Fariduddin, Q., Yusuf, M., Ahmad, I., and Ahmad, A. (2014). Brassinosteroids and their role in response of plants to abiotic stresses. Biol. Plantarum 58, 9–17. doi: 10.1007/s10535-013-0374-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Geng, C., Charles, W., and Liu, S. T. (2011). Overview of available methods for diverse RNA-Seq data analyses. Sci. China Life Sci. 54, 1121–1128. doi: 10.1007/s11427-011-4255-x

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Goupila, P., Souguira, D., Ferjani, E., Faure, O., Hitmi, A., and Ledoigt, G. (2009). Expression of stress-related genes in tomato plants exposed to arsenic and chromium in nutrient solution. J. Plant Physiol. 166, 1446–1452. doi: 10.1016/j.jplph.2009.01.015

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Gupta, N., Khan, D. K., and Santra, S. C. (2012). Heavy metal accumulation in vegetables grown in a long-term wastewater-irrigated agricultural land of tropical India. Environ. Monit. Assess. 184, 6673–6682. doi: 10.1007/s10661-011-2450-7

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hall, J. L. (2002). Cellular mechanisms for heavy metal detoxification and tolerance. J. Exp. Bot. 53, 1–11. doi: 10.1093/jexbot/53.366.1

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Hayat, S., Khalique, G., Irfan, M., Wani, A. F., Tripathi, B. N., and Ahmad, A. (2011). Physiological changes induced by chromium stress in plants: an overview. Protoplasma 249, 599–611. doi: 10.1007/s00709-011-0331-0

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, 480–484. doi: 10.1093/nar/gkm882

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Katz, S. A., and Salem, H. (1994). The biological and environmental chemistry of chromium. J. Appl. Toxicol. 15, 337.

Google Scholar

Kohler, A., Blaudez, D., Chalot, M., and Martin, F. (2004). Cloning and expression of multiple metallothioneins from hybrid poplar. New Phytol. 164, 83–93. doi: 10.1111/j.1469-8137.2004.01168.x

CrossRef Full Text | Google Scholar

Lan, H. X., Wang, Z. F., Wang, Q. H., Wang, M. M., Bao, Y. M., Huang, J., et al. (2013). Characterization of a vacuolar zinc transporter OZT1 in rice (Oryza sativa L.). Mol. Biol. Rep. 40, 1201–1210. doi: 10.1007/s11033-012-2162-2

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Li, J. M., Liu, B., Cheng, F., Wang, X. W., Aarts, M. G. M., and Wu, J. (2014). Expression profiling reveals functionally redundant multiple-copy genes related to zinc, iron and cadmium responses in Brassica rapa. New Phytol. 203, 182–194. doi: 10.1111/nph.12803

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Li, R., Yu, C., Li, Y., Lam, T. W., Yiu, S. M., Kristiansen, K., et al. (2009). SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25, 1966–1967. doi: 10.1093/bioinformatics/btp336

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Liu, T., Liu, S., Guan, H., Ma, L., Chen, Z., Gu, H., et al. (2009). Transcriptional profiling of Arabidopsis seedlings in response to heavy metal lead (Pb). Environ. Exp. Bot. 67, 377–386. doi: 10.1016/j.envexpbot.2009.03.016

CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2^-ΔΔC_T Method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mochida, K., and Shinozaki, K. (2011). Advances in omics and bioinformatics tools for systems analyses of plant functions. Plant Cell Physiol. 52, 2017–2038. doi: 10.1093/pcp/pcr153

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Mortazavi, A., Williams, B. A., Mccue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNASeq. Nat. Methods 5, 621–628. doi: 10.1038/nmeth.1226

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Nachappa, P., Levy, J., and Tamborindeguy, C. (2012). Transcriptome analyses of Bactericera cockerelli adults in response to “Candidatus Liberibacter solanacearum” infection. Mol. Genet. Genomics 287, 803–817. doi: 10.1007/s00438-012-0713-9

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Narusaka, Y., Narusaka, M., Seki, M., Umezawa, T., Ishida, J., Nakajima, M., et al. (2004). Crosstalk in the responses to abiotic and biotic stresses in Arabidopsis: analysis of gene expression in cytochrome P450 gene superfamily by cDNA microarray. Plant Mol. Biol. 55, 327–342. doi: 10.1007/s11103-004-0685-1

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Panda, S. K., and Choudhury, S. (2005). Chromium stress in plants. Plant Physiol. 17, 95–102. doi: 10.1590/S1677-04202005000100008

CrossRef Full Text | Google Scholar

Rai, V., and Mehrotra, S. (2008). Chromium-induced changes in ultramorphology and secondary metabolites of Phyllanthus amarus Schum & Thonn – an hepatoprotective plant. Environ. Monit. Assess 147, 307–315. doi: 10.1007/s10661-007-0122-4

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Sanchita, Dhawan, S. S., and Sharma, A. (2014). Analysis of differentially expressed genes in abiotic stress response and their role in signal transduction pathways. Protoplasma 251, 81–91. doi: 10.1007/s00709-013-0528-5

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Shan, X. H., Li, Y. D., Jiang, Y., Jiang, Z. L., Hao, W. Y., and Yuan, Y. P. (2013). Transcriptome profile analysis of maize seedlings in response to high-salinity, drought and cold stresses by deep sequencing. Plant Mol. Biol. Rep. 31, 1485–1491. doi: 10.1007/s11105-013-0622-z

CrossRef Full Text | Google Scholar

Shanker, A. K., Cervantes, C., Tavera, H., and Avudainayagam, S. (2005). Chromium toxicity in plants. Enviorn. Int. 31, 739–753. doi: 10.1016/j.envint.2005.02.003

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Shanker, A. K., Djanaguiramanb, M., and Venkateswarlu, B. (2009). Chromium interactions in plants: current status and future strategies. Metallomics 1, 375–383. doi: 10.1039/B904571F

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Sharma, J., and Chakraverty, N. (2013). “Mechanism of plant tolerance in response to heavy metals,” in Molecular Stress Physiology of Plants, eds G. R. Rout and A. B. Das (India: Springer), 289–308.

Google Scholar

Strickler, S. R., Bombarely, A., and Mueller, L. A. (2012). Designing a transcriptome next-generation sequencing project for a nonmodel plant species. Am. J. Bot. 99, 257–266. doi: 10.3732/ajb.1100292

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Thapa, G., Sadhukhan, A., Panda, S. K., and Sahoo, L. (2012). Molecular mechanistic model of plant heavy metal tolerance. Biometals 25, 489–505. doi: 10.1007/s10534-012-9541-y

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Vázquez, M. N., Guerrero, Y. R., González, L. M., and Noval, W. T. (2013). Brassinosteroids and plant responses to heavy metal stress. Open J. Metal 3, 34–41. doi: 10.4236/ojmetal.2013.32A1005

CrossRef Full Text

Wang, L. K., Feng, Z. X., Wang, X., Wang, X. W., and Zhang, X. G. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. doi: 10.1093/bioinformatics/btp612

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Wang, Y., Xu, L., Chen, Y. L., Shen, H., Gong, Y. Q., Limera, C., et al. (2013). Transcriptome profiling of radish (Raphanus sativus L.) root and identification of genes involved in response to lead (Pb) stress with next generation sequencing. PLoS ONE 8:e66539. doi: 10.1371/journal.pone.0066539

CrossRef Full Text | Google Scholar

Xu, Y., Zhu, X., Gong, Y., Xu, L., Wang, Y., and Liu, L. (2012). Evaluation of reference genes for gene expression studies in radish (Raphanus sativus L.) using quantitative real-time PCR. Biochem. Biophys. Res. Co. 424, 398–403. doi: 10.1016/j.bbrc.2012.06.119

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Yadav, S. K. (2010). Heavy metals toxicity in plants: an overview on the role of glutathione and phytochelatins in heavy metal stress tolerance of plants. South Afr. J. Bot. 76, 167–179. doi: 10.1016/j.sajb.2009.10.007

CrossRef Full Text | Google Scholar

Yamaguchi, H., Fukuoka, H., Arao, T., Ohyama, A., Nunome, T., Miyatake, K., et al. (2010). Gene expression analysis in cadmium-stressed roots of a low cadmium-accumulating solanaceous plant, Solanum torvum. J. Exp. Bot. 61, 423–437. doi: 10.1093/jxb/erp313

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z. J., et al. (2006). WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34, 293–297. doi: 10.1093/nar/gkl031

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Yu, L. J., Luo, Y. F., Liao, B., Xie, L. J., Chen, L., Xiao, S., et al. (2012). Comparative transcriptome analysis of transporters, phytohormone and lipid metabolism pathways in response to arsenic stress in rice (Oryza sativa). New Phytol. 195, 97–112. doi: 10.1111/j.1469-8137.2012.04154.x

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Zayed, A., Lytle, C. M., Qian, J. H., and Terry, N. (1998). Chromium accumulation, translocation and chemical speciation in vegetable crops. Planta 206, 293–299. doi: 10.1007/s004250050403

CrossRef Full Text | Google Scholar

Zeng, F. R., Qiu, B. Y., Wu, X. J., Niu, S. Z., Wu, F. B., and Zhang, J. P. (2012). Glutathione-Mediated alleviation of chromium toxicity in rice plants. Biol. Trace Elem. Res. 148, 255–263. doi: 10.1007/s12011-012-9362-4

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Zhang, S. P., Xiao, Y. N., Zhao, J. R., Wang, F. G., and Zheng, Y. L. (2013). Digital gene expression analysis of early root infection resistance to Sporisorium reilianumf. sp.zeae in maize. Mol. Genet. Genomics 288, 21–37. doi: 10.1007/s00438-012-0727-3

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Zhang, X. M., Zhao, L., Larson-Rabin, Z., Li, D. Z., and Guo, Z. H. (2012). De novo sequencing and characterization of the floral transcriptome of Dendrocalamus latiflorus (Poaceae: Bambusoideae). PLoS ONE 7:e42082. doi: 10.1371/journal.pone.0042082

CrossRef Full Text | Google Scholar

Zitka, O., Krystofova, O., Hynek, D., Sobrova, P., Kaiser, J., Sochor, J., et al. (2013). “Metal transporters in plants,” in Heavy Metal Stress in Plants, eds D. K. Gupta, F. J. Corpas, and J. M. Palma (Berlin; Heidelberg: Springer), 19–41.

Google Scholar