Edited by: David W. M. Leung, University of Canterbury, New Zealand
Reviewed by: Gong-yin Ye, Zhejiang University, China; Naghabushana K. Nayidu, University of Saskatchewan, Canada
*Correspondence: Liwang Liu, National Key Laboratory of Crop Genetics and Germplasm Enhancement, College of Horticulture, Nanjing Agricultural University, Weigang 1, Nanjing 210095, China
This article was submitted to Plant Biotechnology, a section of the journal Frontiers in Plant Science
†These authors have contributed equally to this work.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Radish (
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,
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,
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.,
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 K2Cr2O7 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.
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 K2Cr2O7 (Cr6+) 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.
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−1 K2Cr2O7 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.
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.
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.,
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.,
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
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
Total reads | 37,948,330 | 100 | 27,451,966 | 100 |
Connector reads | 7444 | 0.02 | 2,799,028 | 10.20 |
Low-length and low-quality reads | 11,559,006 | 30.46 | 12,227,330 | 44.54 |
Clean reads | 26,381,880 | 69.52 | 15,026,358 | 54.74 |
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
COG analysis showed that a total of 15,590 transcript sequences with high homology were grouped into 23 functional categories (Figure
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 log2ratio ≥ 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
GO analysis was performed to predict functions by mapping each DEGs into the records of the GO database. According to Bonferroni-corrected
With Pathway-based analysis specific biological functions could be assigned to genes (Figure
Drug metabolism—cytochrome P450 | ko00982 | 21 | 70 | 2.29E-09 | 7.37E-07 |
Metabolism of xenobiotics by cytochrome P450 | ko00980 | 19 | 67 | 3.65E-08 | 5.88E–06 |
Protein processing in endoplasmic reticulum | ko04141 | 61 | 472 | 3.58E-07 | 2.88E-05 |
Glutathione metabolism | ko00480 | 25 | 134 | 1.99E-06 | 0.000106593 |
Cyanoamino acid metabolism | ko00460 | 15 | 74 | 8.03E-05 | 0.002872626 |
Alpha-Linolenic acid metabolism | ko00592 | 12 | 55 | 0.00019641 | 0.006324488 |
Glucosinolate biosynthesis | ko00966 | 8 | 28 | 0.00032714 | 0.008778165 |
Glycine, serine and threonine metabolism | ko00260 | 19 | 127 | 0.0006472 | 0.015116907 |
Tropane, piperidine, and pyridine alkaloid biosynthesis | ko00960 | 9 | 39 | 0.00079123 | 0.016985096 |
One carbon pool by folate | ko00670 | 11 | 59 | 0.00144753 | 0.027417844 |
Tryptophan metabolism | ko00380 | 13 | 81 | 0.00233318 | 0.037564152 |
Pentose and glucuronate interconversions | ko00040 | 30 | 100 | 2.19E-14 | 7.48E-12 |
Starch and sucrose metabolism | ko00500 | 51 | 331 | 9.66E-11 | 1.65E-08 |
Amino sugar and nucleotide sugar metabolism | ko00520 | 48 | 327 | 1.79E-09 | 2.04E-07 |
Phenylalanine metabolism | ko00360 | 25 | 120 | 1.44E-08 | 1.23E-06 |
Phenylpropanoid biosynthesis | ko00940 | 28 | 153 | 3.97E-08 | 2.72E-06 |
Cysteine and methionine metabolism | ko00270 | 29 | 193 | 1.73E-06 | 9.86E-05 |
Glucosinolate biosynthesis | ko00966 | 9 | 28 | 1.60E-05 | 0.00068345 |
Inositol phosphate metabolism | ko00562 | 23 | 157 | 2.98E-05 | 0.00113249 |
Plant hormone signal transduction | ko04075 | 57 | 601 | 0.000123 | 0.00350661 |
Methane metabolism | ko00680 | 25 | 203 | 0.000244 | 0.0059517 |
Nitrogen metabolism | ko00910 | 18 | 127 | 0.000319 | 0.00728052 |
Glutamatergic synapse | ko04724 | 13 | 80 | 0.000565 | 0.01206837 |
Selenocompound metabolism | ko00450 | 9 | 45 | 0.000854 | 0.01717926 |
Neuroactive ligand-receptor interaction | ko04080 | 4 | 11 | 0.002514 | 0.04776096 |
Mineral absorption | ko04978 | 7 | 34 | 0.002671 | 0.04808329 |
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 (
Comp22306_c0_seq1 | 4.55530704 | 2.07E-06 | 2.46E-05 | MAPKKK19 |
Comp47424_c0_seq1 | 5.31780773 | 5.79E-05 | 0.000406339 | Putative calcium-binding protein CML44 |
Comp26037_c0_seq1 | 4.78683836 | 4.81E-10 | 1.92E-08 | bZIP transcription factor 60 |
Comp28465_c1_seq10 | 34.4390032 | 4.66E-07 | 6.86E-06 | Ethylene-responsive transcription factor ABR1 |
Comp19785_c2_seq2 | 33.6145418 | 3.57E-05 | 0.000271789 | MYB transcription factor |
Comp13434_c0_seq1 | 3.38377099 | 4.68E-06 | 4.93E-05 | WRKY DNA-binding protein 15 |
Comp23733_c0_seq1 | −34.180089 | 9.66E-06 | 9.02E-05 | Peroxidase 17 |
Comp21825_c0_seq1 | −35.687884 | 2.01E-09 | 6.63E-08 | Peroxidase 21 |
Comp24452_c0_seq1 | 4.77201196 | 6.17E-10 | 2.40E-08 | Peroxidase 34 |
Comp22278_c0_seq2 | −4.7587006 | 9.50E-05 | 0.000622163 | Peroxidase 49 |
Comp25683_c0_seq1 | −7.5525304 | 1.36E-16 | 3.97E-14 | Peroxidase 64 |
Comp26981_c0_seq20 | 4.34774293 | 3.70E-07 | 5.63E-06 | Multidrug resistance protein ABC transporter family |
Comp27391_c0_seq2 | 34.0239657 | 4.28E-06 | 4.59E-05 | ABC transporter-like protein |
Comp21544_c0_seq1 | 7.14134258 | 3.71E-13 | 3.81E-11 | MATE efflux family protein |
Comp50587_c0_seq1 | 6.15734306 | 8.76E-07 | 1.18E-05 | Ethylene-responsive transcription factor 9 |
Comp25404_c0_seq3 | 38.2011666 | 1.57E-17 | 6.33E-15 | DNAJ heat shock protein-like protein |
Comp25620_c0_seq8 | 34.0869034 | 2.94E-06 | 3.31E-05 | HSP20-like chaperone |
Comp27769_c0_seq1 | 3.86061292 | 1.96E-07 | 3.28E-06 | Glutathione synthetase |
Comp25459_c0_seq6 | 33.4053948 | 9.64E-05 | 0.000628082 | Cysteine synthase |
Comp18928_c2_seq1 | 8.72713911 | 4.38E-20 | 4.06E-17 | Glutathione S-transferase |
Comp25826_c0_seq2 | 5.81512529 | 1.37E-09 | 4.79E-08 | Phytochelatin synthase 1 |
Comp27638_c0_seq6 | −34.141249 | 1.14E-05 | 0.000103515 | Phytochelatin synthetase-like protein |
Comp26972_c1_seq1 | −3.2389198 | 8.27E-06 | 7.91E-05 | Metallothionein type 3 |
Comp14108_c0_seq2 | 34.8424264 | 4.62E-08 | 9.53E-07 | Metallothionein-like protein |
Comp20396_c0_seq1 | 33.8349318 | 1.17E-05 | 0.000105136 | Cytochrome P450 79b3 |
Comp20687_c1_seq1 | −38.63336 | 8.55E-18 | 3.82E-15 | Cytochrome P450 79f1 |
Comp28927_c0_seq1 | −6.435114 | 6.56E-14 | 8.60E-12 | Cytochrome P450 83a1 |
Comp25273_c0_seq1 | 7.08334247 | 2.61E-13 | 2.84E-11 | Cytochrome P450, family 710, subfamily A |
Comp29019_c0_seq1 | 5.02595565 | 1.09E-10 | 5.21E-09 | CYP83B1 |
Comp24070_c1_seq2 | 39.6331444 | 8.54E-22 | 1.32E-18 | Steroid sulfotransferase 4 |
Comp14258_c0_seq1 | 33.8349318 | 1.17E-05 | 1.05E-04 | Steroid sulfotransferase 3 |
Comp26378_c0_seq3 | 3.53992894 | 1.08E-05 | 9.96E-05 | Fatty acid hydroxylase 1 |
There were 94 DEGs that were of high homology with different kinds of transcription factors (TFs), such as the WRKY family (i.e.,
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
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
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.,
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 (
Chromium (Cr) is considered as a serious environmental pollutant posing a critical concern to human health through food chain (Zayed et al.,
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.,
Plants exposed to toxic heavy metal (HM) can generate reactive oxygen species (ROS), resulting oxidative stress (Panda and Choudhury,
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.,
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,
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.,
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.,
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,
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.,
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.,
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,
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.
The RNA SEQ raw data have been deposited in NCBI Sequence Read Archive (SRA,
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.
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.
The Supplementary Material for this article can be found online at: