Genome-Wide Transcriptome Analysis Reveals that Cadmium Stress Signaling Controls the Expression of Genes in Drought Stress Signal Pathways in Rice

Plant growth is severely affected by toxic concentrations of the non-essential heavy metal cadmium (Cd). Comprehensive transcriptome analysis by RNA-Seq following cadmium exposure is required to further understand plant responses to Cd and facilitate future systems-based analyses of the underlying regulatory networks. In this study, rice plants were hydroponically treated with 50 µM Cd for 24 hours and ∼60,000 expressed transcripts, including transcripts that could not be characterized by microarray-based approaches, were evaluated. Upregulation of various ROS-scavenging enzymes, chelators and metal transporters demonstrated the appropriate expression profiles to Cd exposure. Gene Ontology enrichment analysis of the responsive transcripts indicated the upregulation of many drought stress-related genes under Cd exposure. Further investigation into the expression of drought stress marker genes such as DREB suggested that expression of genes in several drought stress signal pathways was activated under Cd exposure. Furthermore, qRT-PCR analyses of randomly selected Cd-responsive metal transporter transcripts under various metal ion stresses suggested that the expression of Cd-responsive transcripts might be easily affected by other ions. Our transcriptome analysis demonstrated a new transcriptional network linking Cd and drought stresses in rice. Considering our data and that Cd is a non-essential metal, the network underlying Cd stress responses and tolerance, which plants have developed to adapt to other stresses, could help to acclimate to Cd exposure. Our examination of this transcriptional network provides useful information for further studies of the molecular mechanisms of plant adaptation to Cd exposure and the improvement of tolerance in crop species.


Introduction
Heavy metal ions are highly reactive and toxic to living cells, and their accumulation in plants is a major agricultural problem. Heavy metals can be classified into two categories based on their influence on plant growth. The first category includes essential mineral elements, trace amounts of which may be required for adequate growth and development, and an excess of which is toxic as soon as the concentration exceeds the threshold that plants can endure. The second category includes non-essential metals with recognized toxicity including cadmium (Cd) and lead (Pb). In particular, Cd is absorbed by the roots from the soil and transported to the shoot, negatively affecting nutrient uptake and homeostasis in plants, even at low concentrations. It is also known to adversely impact various biochemical and physiological processes including changes in the transcriptome and proteome of plants, resulting in inhibited root and shoot growth and, ultimately, reduced yield [1,2,3]. Cd pollution in arable soil has dramatically increased worldwide over the last several decades through the use of phosphate fertilizers, sludge, and irrigation water containing Cd. Furthermore, accumulation of Cd in the edible parts of plants such as seed grains places humans at a risk when ingesting them because of its highly toxic effects on human health. Thus, it is important to study the mechanisms of plant responses and defenses to Cd exposure to overcome this problem.
Cd causes oxidative stress and generates reactive oxygen species (ROS) such as superoxide radicals (O 2 2 ) and hydrogen peroxide (H 2 O 2 ), which can cause damage in various ways such as reacting with DNA causing mutation, modifying protein side chains and destroying phospholipids [4]. Many genes that respond to oxidative stress in plants, as well as genes associated with defense systems, have been studied, including induction detoxification enzymes such as glutathione S-transferase (GST), peroxidase (Prx), thioredoxins (Trx), peroxiredoxin (PrxR) and catalase (Cat) under Cd exposure, which confer Cd tolerance in plants [5,6,7,8]. For defense against Cd toxicity, chelation of metal ions by ligands such as cysteine-rich metallothioneins (MT) and phytochelatins (PC) is also induced under Cd exposure. Such molecules bind Cd 2+ ions through S-containing amino acid ligands and sequestrate the

Sample preparation
Rice (Oryza sativa ssp. japonica cv. Nipponbare) seeds were germinated and grown by hydroponic culture in nutrient media [1.425 [23] in a growth chamber at 28uC and 70-80% humidity in a 16h light/8h dark cycle. After 10 days, the seedlings of uniform size and growth were subjected to Cd stress treatment by transferring them to a similar medium with 50 mM CdSO 4 . The plants were maintained under Cd stress conditions for 120 h and then details of plant growth were recorded and sampling was performed as described previously [21]. Total RNA was extracted from all tissue samples using an RNeasy Plant Kit (Qiagen, Hilden, Germany) according to manufacturer's instructions. Construction of 13 cDNA libraries (2 tissues, 3 conditions, and 2-3 replicates) from total RNA using a TruSeq RNA sample preparation kit and sequencing with the Illumina Genome Analyzer IIx (Illumina Inc., San Diego, CA, USA) was performed according to the manufacturer's protocols.
Sequencing and mapping of short reads onto the rice genome More than two biological replicates for each set of conditions were highly correlated (coefficient . 0.92), and reads from the same treatment were merged for subsequent analysis. Trimming of Illumina adaptor sequences and low-quality bases (Q , 20) at the 59 and 39 ends of each read was performed by Cutadapt [24] (http://code.google.com/p/cutadapt/) and a custom-made program. These pre-processed reads were mapped to the IRGSP-1.0 genome assembly (http://rapdb.dna.affrc.go.jp/) to reconstruct the transcript structures by a series of programs; Bowtie for shortread mapping [25], TopHat for defining exon-intron junctions [26], and Cufflinks for gene structure predictions [27]. To estimate the expression levels of each transcript, all pre-processed reads were mapped to the Os-Nipponbare-Reference-IRGSP-1.0 genome assembly (http://rapdb.dna.affrc.go.jp/) by Bowtie with default parameters [28]. The expression level for each transcript was calculated as RPKM (Reads Per Kilobase exon Model per Million mapped reads)-derived read counts [29] based on the number of uniquely mapped reads that overlapped with exonic regions. The resulting RNA-Seq data were deposited in the DDBJ Sequence Read Archive (Accession No. DRA001092).

Analysis of the responsive transcripts and alternative splicing patterns
We performed G-test to detect differentially expressed transcripts in control and Cd treatments based on the statistical null hypothesis that the proportions of mapped reads to the transcripts are the same between the two conditions. The frequency distribution of transcripts was determined by constructing 262 contingency tables with variables corresponding to the number of mapped and unmapped reads on a given transcript in control and each Cd treatment, respectively. A false discovery rate (FDR, 0.01) was used in multiple hypothesis testing to correct for multiple comparisons. When calculating fold changes, 1 was added to avoid division by 0. Gene Ontology (GO) terms were assigned to each transcript from RAP-DB for each GO category. Enrichment of GO terms in the biological process category was evaluated by Fisher's exact test with a FDR threshold of 5% for responsive transcripts in major clusters at 1 h and 24 h after Cd treatment. The results were plotted as 2log10 of FDR values in a heatmap. Microarray data (GSE6901, http://www.ncbi.nlm.nih.gov/geo/) were used to compare the responsive transcripts between Cd exposure (24 h) and other abiotic stresses (drought, salt and cold). The responsive transcripts ($2-fold or # 0.5-fold) in each treatment were used for Venn diagram analysis using the ''venn'' function of the R base package gplot version 2.10.1. For analysis of alternative splicing patterns in Cd-responsive genes, RPM (Reads per Million mapped reads) values of the splice sites in 5,222 representative loci were calculated using reads mapped on the sites. The splicing patterns in sites with RPM values . 1 were compared between control and Cd exposure samples.

qRT-PCR
The expression of Cd-upregulated metal ion transporter genes in root and shoot samples were confirmed by quantitative RT-PCR analysis. Rice seeds were germinated and grown in water in a growth chamber. After 10 days, the seedlings were subjected to different stress treatments by transferring them to water containing different reagents ( for Zn exposure). Ten-day-old rice seedlings transferred to water were used as a control for the stress treatment. Total RNA was extracted from samples collected after 24 h of each treatment. After DNase I (Takara, Shiga, Japan) treatment, first-strand cDNA was synthesized using the Transcriptor First Strand cDNA synthesis kit (Roche, Basel, Switzerland) according to the manufacturer's protocol. The resulting cDNA was used for PCR amplification in the LightCycler 480 system (Roche, Basel, Switzerland) with each primer set. The detection threshold cycle for each reaction was normalized using Ubiquitin1 with 59-CCAGGACAAGATGATCTGCC-39 and 59-AAGAAGCT-GAAGCATCCAGC-39 as primers. Three technical replicates for each treatment were used for analysis.

Changes in plant morphology under Cd exposure
In the hydroponically cultured rice, growth was greatly affected by Cd exposure (Figure 1). In Cd-exposed rice, growth retardation of the shoot was apparent after 24 h, many dark spots appeared on all leaves after 48 h, the leaves turned yellow and the leaf tips of the seedlings began to wilt after 72 h, and after 120 h of 50 mM cadmium exposure all leaf blades were curled completely and the seedlings were shrunken and wilting, which is attributable to supply of Cd ( Figure 1). This concentration or higher concentrations of Cd have been previously shown to elicit robust physiological responses and gene expression as acute toxic responses in rice seedlings [3,30,31]. The wilting occurred gradually compared with drought treatment for 24 h, in which symptoms started to appear after 2 h treatment and plants were completely dried up after 24 h, in the same growth chamber ( Figure S1). The detoxification processes of the plant are insufficient to cope with the toxic metal beyond a 10 mM dose, and wilting and reduction of plant fresh weight are accompanied by decreased leaf conductance and increased stomatal closing [32]. It has been suggested that part of the fatal damage to plants from Cd exposure occurs through drought stress.

Cd treatment induced overall changes in gene expression in rice
To ascertain our hypothesis at the transcriptional level, we generated transcriptome profiles of the early response to Cd exposure using RNA-Seq during plant growth, particularly at 1 h and 24 h after 50 mM Cd treatment and before treatment (0 h). For each set of conditions, an average of approximately 48.3 million (90.7%) quality-evaluated reads were mapped to the rice genome sequence and used for further analysis (Table S1) Table S2), suggesting the effect of Cd exposure was enhanced gradually up to 24 h. Several MTs, Prxs and heat shock proteins (Hsps) were strongly upregulated among the 20 genes with the greatest relative expression at 24 h. Table 1 shows strongly upregulated (top 20) transcripts without probes in the Rice 44K Microarray in roots and shoots at 24 h. We also found transcription factors (e.g. AP2-EREBP, NAC) that may function as regulatory factors under Cd exposure among these novel Cd- responsive transcripts, most of which had uncharacterized functions. Furthermore, the responsive transcripts included 1.84-7.98% unannotated transcripts predicted by the Cufflinks program in each particular tissue/time. A validation experiment on the responsive unannotated transcripts was performed by qRT-PCR analysis [21]. These novel transcripts may represent interesting novel targets to understand Cd regulatory pathways in rice. Moreover, among 44,519 representative loci on the rice genome (IRGSP-1.0), alternative splicing isoforms were confirmed in 5,222 (11.7%) representative loci. Changes in the alternative splicing patterns of these loci were examined against Cd exposure and 2,873 loci (6.5%) showed different patterns. RNA-Seq is a promising tool for analyzing alternative patterns that microarrays are unable to, and approximately ,48% of rice genes showed alternative splicing patterns [33]. This suggests that Cd responsive expression is also regulated by alternative splicing mechanisms with respect to the timing of the response or tissue specificity. In conclusion, RNA-Seq data were far superior to data derived from the microarray and that Cd treatment affected many genes and caused drastic changes in gene expression in rice.

Functional characterization of Cd-responsive transcripts
To investigate the functions of cadmium stress-responsive transcripts, we performed Gene Ontology (GO) analysis and found an overrepresentation of specific GO keywords that denote involvement with processes triggered by stress, including metal ion transport (GO:0030001) (upregulated, root), response to stress (GO:0006950) (upregulated, root and shoot), trehalose biosynthetic process (GO:0005992) (upregulated, shoot), DNA replication (GO:0006260) (downregulated, root), DNA repair (GO:0006281) (downregulated, root), translation (GO:0006412) (downregulated, root and shoot), and photosynthesis (GO:00015979) (downregulated, shoot) ( Figure S2). These enriched GO terms among the upregulated transcripts indicate that RNA-Seq was successful in identifying Cd-responsive genes. The enriched GO terms among downregulated transcripts result in growth retardation ( Figure 1). These results imply that Cd controls a significant part of the defense to stress and plant growth.
Among the GO keywords identified, response to stress (GO:0006950) included many drought stress responsive genes such as the Arabidopsis cor47 homolog Dip1 [34], Rab 21 ( = Rab16A) (responsive to abscisic acid 21) genes [35], and other Rab 16 genes [36] in roots and shoots at 24 h under Cd exposure (Table S3). The signaling pathways in drought, high-salinity and low temperature stress show high levels of cross-talk among abiotic stresses in Arabidopsis [37]. Our results suggested that expression of these genes responsive to drought stress (also implying a relationship to highsalinity and low temperature stresses) was affected by Cd exposure, but a relationship between Cd and drought stresses at the transcription level has not been reported as far as we know.

Cd exposure triggered upregulation of drought stress responsive genes
We investigated the expression of well-characterized drought stress-related transcription factors (TFs) and their downstream genes under Cd exposure to ascertain the relationship between Cd and drought stresses.
i) DREB/CBF TFs. Among the DREB/CBF TFs, which contain drought-responsive elements (DRE), DREB1A/CBF3 and DREB1G were drastically upregulated in roots after 24 h of Cd exposure and DREB1C/CBF2 was drastically upregulated in roots after 1 h ( Figure 3A). Over-expression of DREB1A/CBF3 improved tolerance to drought, high-salinity and low temperature, and resulted in growth retardation under non-stress conditions, suggesting functional similarity to Arabidopsis [38,39]. Overexpression of DREB1G significantly elevated tolerance to drought [40]. As the transcriptional networks mediated by DREB1 are conserved in plants and DREB1 regulates many downstream target genes under abiotic stress [37], part of the signaling pathway related to DREB1 might be enhanced under Cd exposure in rice.
ii) bZIP TFs. Several bZIP domain TFs were characterized that regulated expression of ABA-responsive genes in rice. The abiotic stress-inducible genes bZIP23 [41] and ABI5 [42] were upregulated in roots at 1 h and 24 h during Cd exposure, respectively ( Figure 3B). bZIP23 and ABI5 are thought to function as key components in ABA-dependent transcriptional networks in rice, suggesting the existence of gene expression regulation through ABA-dependent pathways under Cd exposure. The existence of an ABA-dependent pathway was supported by strong upregulation of NCED3 and NCED9 in roots (Table S3). In our BLASTP search, the Arabidopsis proteins AtNCED3 [43] and AtNCED9 [44], key enzymes of ABA synthesis during drought   stress and seed development, respectively, were the top hits to NCED3 and NCED9 in rice. ABA is known to accumulate rapidly and confer stress tolerance by inducing many genes under drought [43,44,45,46]. Synthesis of ABA following Cd stress can lead to decreased leaf conductance as Cd 2+ toxicity perturbs the plantwater relationship [32] (Figure 1). LEA3-1 was upregulated in bZIP23 overexpressors [41], and here it was also upregulated in roots and shoots (Table S3). LEA genes often appear to be ABAdependent [47], and their proteins may function in protecting macromolecules such as enzymes, protein complexes and membranes [48]. Overexpression of LEA3-1 in rice can significantly improve relative yield under drought stress [49].
iii) NAC TFs. Several NAC TFs, including NAC5 [50] and SNAC1 [51], are induced by drought, high-salinity and lowtemperature. NAC5 was upregulated in roots at 1 h and shoots at 24 h, and SNAC1 was upregulated in roots at 24 h during Cd exposure ( Figure 3C). NAC5 is regulated through an ABAdependent pathway [52]. LEA3, which was identified as a NAC5 downstream gene [50,53], and a rice ERD1 (early responsive to drought 1) homologue (OsERD1) were also upregulated in roots and shoots (Table S3). The expression of OsERD1 is probably regulated by SNAC1 through the NAC recognized sequence (NACRS) promoter [51,54]. Over-expression of SNAC1 can significantly improve drought resistance and confer strong tolerance to high-salinity stress [51].
iv) Other TFs. Many other stress-related TFs, including AP2/ERFs (AP37, AP59), C2H2 zinc finger (ZFP252), TIFY (TIFY11) and MYB (Myb4), were also upregulated under Cd exposure (Table S3). Overexpression of AP37, AP59 [55] and ZFP252 [56] results in drought and high-salinity tolerance. Overexpression of ZFP252 in rice increased the amount of soluble sugars and free proline [56]. The accumulation of soluble raffinose and the upregulation of galactinol synthase genes, key enzymes of raffinose synthesis, have been reported in drought-, high salinityand low temperature-treated Arabidopsis plants [57]. In our study, galactinol synthase and raffinose synthase genes were strongly upregulated in roots under Cd exposure (Table S3). Cd treatment has also been shown to increase the level of raffinose in Arabidopsis [58]. The AtGolS3 galactinol synthase gene is controlled by DREB1A through DRE and DRE-like cis-elements in Arabidopsis [57,59]. These results suggest raffinose is a metabolite for adaptation to Cd and might be produced by enhancing a DREB-related pathway in rice. P5CS, a key enzyme of osmoprotectant proline synthesis, was upregulated in shoots under Cd exposure (Table S3). Accumulation of proline as an enzyme protectant has been reported in rice seedlings exposed to Cd [60,61], suggesting that proline is also a metabolite for adaptation to Cd. Overexpression of TIFY11a [62] resulted in high-salinity and drought tolerance, and overexpression of Myb4 [63] resulted in tolerance to low temperature, so upregulation of TIFY11 and Myb4 may also function in tolerance to Cd exposure (Table S3). Thus, many reports suggest that overexpression of stress-inducible TFs can increase abiotic stress tolerance to drought, high-salinity, or low temperature in plants, so upregulation of such TFs could be useful for developing transgenic crops with enhanced tolerance to Cd stress. Some (DREB, bZIPs and NCED) were more upregulated in roots than in shoots under Cd exposure up to 24 h and the tendency of their expression patterns was confirmed by qRT-PCR analysis, which was different to the tendency under drought stress for 3 h ( Figure S3, Table S4). As OsDREB1A and OsDREB1B were included in the 50 genes with the greatest relative expression in roots exposed to 10 mM Cd-treatment for 3 h [20], roots might be affected more directly and earlier by Cd exposure. In conclusion, our data clearly indicated that Cd affects the expression of genes in the drought-related signaling pathway, especially in roots, because the Cd exposure treatment was performed in hydroponic culture.
We also confirmed the upregulation of Jacalin1, LOX (lipoxygenase), BBTI 1 (Bowman Birk trypsin inhibitor 1), Receptor kinase containing LRR repeats, PSLS (Phospho sulfolactate synthase), Hsp70, PP2Ca and PP2Cb in one or more of the tissue/treatment combinations (Table S3). The upregulation of these genes was ascertained to help plants acclimate to stress conditions, such as drought in AtCBF3-or AtABF3-overexpressing transgenic rice [64]. These results suggest that Cd exposure more or less perturbs the expression of genes in drought, high salinity-and low-temperature stress signaling pathways, which results in similar morphological changes in response to both Cd and drought (Figure 1, Figure S1). Even though comparative analysis using microarray data revealed that 28.0-54.3% of the responsive transcripts ($2-fold or # 0.5-fold) responded to other abiotic stresses (drought, high salinity or low-temperature) in roots and shoots ( Figure S4), Cd-responsive transcripts responded to drought/high-salinity stresses more than twice as much as to low temperature stress ( Figure S4), suggesting there are higher levels of cross-talk between Cd and drought/high-salinity stresses than in the cross-talk between Cd and low temperature stress as expected from observation of wilted seedlings under Cd exposure ( Figure 1). We also identified that 9% of Cd-upregulated transcripts were commonly upregulated among the four stresses, including SNAC1, ZEP252 and Myb4 (Table S5), which may confer general adaptation under abiotic stress.
Upregulation of genes in the defense system to Cd exposure and RNA-Seq revealed the regulated genes of multigenic families Next, we investigated the expression of antioxidant and detoxification enzymes to confirm that our expression profiles reflected the response to Cd exposure.
The production of ROS is unavoidable under various biotic and abiotic stresses including Cd exposure. The accumulation of ROS is largely counteracted by an intricate antioxidant defense system [6]. The enzymatic scavengers GST [65] (Pfam accession: PF02798.15, PF13417.1, PF13410.1, PF00043.20, http://pfam. sanger.ac.uk/) and class III Peroxidase (Prx) (E.C. 1.11.1.7) [66] have large and complex gene families that control key metabolic steps in many eukaryotic systems. Among 80 GST genes, the responsive genes tended to be more upregulated in roots at 24 h compared with the other tissue/treatment combinations. In particular, GSTU12 (shoot), GSTU19 (root and shoot), GSTU24 (root and shoot), GSTU30 (shoot), GSTU31 (root) and GSTU39 (root) were strongly upregulated [fold change (FC) . 20] ( Figure 4A). GSTU19, GSTU24 and GSTU39 were found to be upregulated in the roots of rice seedlings under arsenate exposure using a microarray platform [67], suggesting that part of the Cd defense might be similar to the As defense in the roots of rice. Enhanced GST with glutathione peroxidase activity in transgenic tobacco increased glutathione-dependent peroxide scavenging and alterations in glutathione and ascorbate metabolism that led to reduced oxidative damage [68]. Class III Prx genes catalyze the reduction of H 2 O 2 by transferring electrons to various donor molecules such as lignin precursors and secondary metabolites [69,70]. Prx22, Prx59, Prx62, Prx86, Prx110, Prx111 and Prx131 among 138 Prxs were upregulated strongly (FC . 20) in shoots after 24 h of Cd exposure ( Figure 4A). Prx86 and Prx111 were strongly upregulated in the plant defense against a gall midge [71] and against Xanthomonas oryzae and Magnaporthe grisea [72]. As many i) Antioxidant enzymes for ROS-scavenging.
studies have described the diverse functions of Prx in Arabidopsis [8,73], these genes might be expressed in a time and tissue specific manner in rice.
Genes encoding various kinds of ROS-scavenging enzymes, such as glutaredoxin (Grx) [74], Trx [75], PrxR [76], monodehydroascorbate reductase (MDAR) [77], alternative oxidase (AOX) [78], Cat (PF00199), and alpha-dioxygenase (a-DOX1) [79] were also upregulated (FC . 3) (Figure S5), suggesting activation of several antioxidative systems in the plant cells under Cd exposure. Most of these genes tended to be upregulated more at 24 h compared with 1 h after Cd exposure. a-DOX1, which plays a role in protecting tissues from oxidative damage [80], was the most upregulated (,94-fold) in shoots at 24 h under Cd exposure. a-DOX1 is well-known to be upregulated and show increased activity under oxidative-and heavy metal-stresses [80]. Tissue-specific upregulation was observed for some genes, such as Os01g0194600 (Grx) in roots and a-DOX1 in shoots. Even though the genes were from the same family, we observed tissue-specific expression for different PrxR and AOX genes. One of the two genes in each family was upregulated more in roots compared with shoots and the other gene showed an opposite upregulation pattern ( Figure S5). Respiratory burst oxidase homolog (Rboh: NADPH oxidase) genes that produce ROS [6,81] were also upregulated ( Figure S5). These results indicate that the treatment evokes ROS in rice cells. However, while the gene expression profiles of the whole multigenic family have not been well investigated under Cd exposure, the regulated genes in the multigenic family were elucidated clearly by RNA-Seq analysis. Through the use of gene transfer technology to manipulate the activities of these antioxidant enzymes, ROS-scavenging systems in plants might be used to increase plant stress tolerance. ii) Detoxification enzyme through chelation. The chelation of heavy metals by particular ligands is a well-characterized mechanism of Cd detoxification that has evolved in plants. The two best-characterized heavy metal-binding ligands in plant cells are PC and MT, which are involved in detoxification [9]. PC is The expression of Cd-upregulated heavy metal ion transporters was investigated in various liquid media containing different kinds of metal ions by qRT-PCR analysis. The x-axis shows treatments and the y-axis shows relative expression. Transcript expression levels were normalized using an internal control (ubiquitin 1) and plotted relative to expression in water (control) at 24 h in roots (R) and shoots (S). The transcripts were classified into four groups based on their expression patterns. doi:10.1371/journal.pone.0096946.g005 synthesized from reduced glutathione (GSH) in a transpeptidation reaction [82]. PC synthase (PF05023.9) and GSH synthase (PF04107, PF03917, PF03199) genes did not show much response to Cd exposure in this study, but cad1, a PC deficiency mutant in Arabidopsis, is Cd sensitive [83]. MTs have been found in diverse organisms including mammals, plants, and fungi as well as some prokaryotes and function in maintaining homeostasis of essential metals and metal detoxification [1,9]. MTs are rich in cysteine residues that coordinate multiple metal atoms such as cadmium, zinc and copper. Eleven MT genes were found to be differentially expressed during growth and development, in various tissues and during biotic and abiotic stresses [84,85]. Here, 9 of these 11 MT genes were strongly or weakly upregulated in at least one set of conditions ( Figure 4B). Interestingly, two sets of three MT genes (Os12g0567800, Os12g0568200, Os12g0568500 and Os12g0570700, Os12g0571000, Os12g0571100) reside on chromosome 12 within 40-kb and 15-kb regions, respectively. As the clustering of MTs is well-known in mouse and human [86], these might have evolved from gene duplications in rice. These clustered genes were upregulated at different levels in at least two of the tissue/treatment combinations ( Figure 4B). Os12g0568200 showed the strongest upregulation (821-fold) among them, in roots at 24 h under Cd exposure ( Figure 4B). Simultaneously, cis-regulatory changes and minor changes in the regulatory regions of genes may occur during evolution [87,88], resulting in diversity of gene expression levels. Members of the rice MT gene clusters differed in their tissue expression patterns, suggesting that each gene may perform different functions in specific tissues. The responsive transcripts of antioxidative and detoxification enzymes differed in their expression patterns among metal stresses ( Figure S6, Table  S4). Their upregulation indicated that our expression profiles under Cd exposure obtained by RNA-Seq analysis were reliable.

Responses of various metal transporter genes under Cd and other metal stresses
To confirm whether other ions affect the expression of Cdresponsive genes, we investigated the expression of metal transporter genes from the HMA, MatE (multi antimicrobial extrusion protein) [89], Zip, CDF, NRAMP, and PDR families, as well as LCT1 [12] and LCT1 homologues (Os06g0576200, Os06g0674000), because at least one gene from each of these families and LCT1 have been reported to function in Cd transport. As LCT1 did not have any Pfam domain for heavy metal binding, we searched for homologues in RAP DB using a BLASTP search (threshold of 1e-03). These gene families contain specific metal ion binding domains [PF01554 (MatE), PF08370 (PDR_assoc), PF01545 (Cation_efflux), PF02535 (Zip), PF00403 (HMA), PF01566 (Nramp)], which may function in Cd transport. In total, 88.9% of the transporters (168 transcripts) were responsive to Cd exposure, with 35.7% (60 transcripts) in the HMA family, 29.8% (50 transcripts) in MatE, 10.7% (18 transcripts) in Zip, 7.7% (13 transcripts) in Cation_efflux and Nramp, 6.0% (10 transcripts) in PDR_assoc and 2.4% (4 transcripts) being LCT1 and LCT1 homologue. PDR9 [11], LCT1 [12] and HMA3 [13], which are important of Cd translocation in rice, were identified as the responsive transcripts. Transcripts upregulated more than 5-fold in one or more of the tissue/treatment combinations are shown in Table 2. Os02t0585200 containing an HMA domain was the most upregulated in roots (87-fold), and Os01t0609900 containing a PDR_assoc domain was the most upregulated in shoots (53-fold) at 24 h during Cd exposure. These results emphasized the potential of the RNA-Seq strategy to reveal novel Cd-responsive transporter transcripts in rice. It should be mentioned that the expression levels were diverse among the gene families, suggesting they may have differentiated functions in transporting various metal ions.
Next, the expression of 12 randomly selected upregulated transporter genes was analyzed by qRT-PCR in seedlings treated with various metals, such as Zn, Mg and Pb ( Figure 5, Table S4). The seedlings were grown in water for 10 days to avoid priming effects before transferring to various media. The expression patterns of the Cd-upregulated transporter transcripts were classified into four groups: three transcripts (MatE, HMAs) upregulated significantly under Cd exposure (group 1), three transcripts (HMAs) upregulated significantly under both Cd and the other treatments (group 2), two transcripts (MatE, Zip) upregulated under the other treatments but not Cd (group 3), and four transcripts (Cation efflux, MatE, PDR_assoc, HMA) upregulated under general stress (group 4), though their expression might have been affected by ion concentration, treatment time, the balance of other ions, tissue and growth stage. The transcripts of group 1 might respond specifically to Cd. In our conditions, the transcripts of group 2 were not only upregulated by Cd, but also by the essential metal Ni, suggesting they might function in Ni transport. The nickel transporter TjZnt in the Ni hyperaccumulator Thlaspi japonicum has been reported to have Cd 2+ transport ability [90]. The transcripts of group 3 did not respond to Cd in water substituted for nutrient-rich medium, suggesting they were easily affected by other metal ions. The transcripts of group 4 might respond to general stress. The expression of most of these was affected to different degrees by exposure to other metal ions, suggesting specific systems for transporting Cd may have not developed in rice because Cd is a non-essential metal for the plant.

Conclusions
In summary, the sequencing of mRNAs uncovered new genes involved in signal transduction, antioxidation, detoxification and metal transport, and enabled us to develop new hypotheses for the transcriptional network underlying the response to Cd exposure. A highlight was the discovery of a relationship between Cd stress and drought stress. The data suggested that the Cd stress signaling pathway is involved in controlling drought stress signaling pathways that might confer to tolerance to Cd exposure. It is quite possible that the network underlying Cd stress responses and tolerance, which plants have developed to adapt to other stresses, could help to acclimate to Cd exposure because it is a non-essential metal. We have summarized the scheme of signal transduction for acclimation to Cd exposure in rice in Figure 6. As ROS are a toxic byproduct of aerobic metabolism, but also act as signaling molecules in the complex signaling network of cells [5], the generated ROS might act as messengers to activate gene expression in both Cd and abiotic stress pathways. Understanding the relationships between the transduction of different stress signals is useful to develop transgenic rice that show enhanced stress tolerance. Establishing the exact composition and organization of the transcriptional network underlying the response to Cd exposure will provide a robust tool for manipulating the stress tolerance of crops in the future. Figure S1 Changes in rice morphology after 24 h drought stress. Rice seedlings grown by hydroponic culture in nutrient media were subjected to drought stress treatment by transferring them to a case without media. The rice seedlings began to show signs of wilting in shoots after 1 h and the changes gradually became more prominent, such that after 24 h of drought stress the shoots were completely wilted compared with control rice not removed from the nutrient media. (TIFF) Figure S2 Identification of GO terms enriched in Cdresponsive transcripts. Significant GO terms identified by GO enrichment analysis based on the most enriched biological processes associated with variations under Cd exposure are shown in a heatmap of responsive transcripts (left: upregulated transcripts, right: downregulated transcripts). The bar with red-black gradation indicates the level of significance of GO enrichment with the extremes representing statistically significant (red) and non-significant (black) GO terms. (TIFF) Figure S3 qRT-PCR analysis of Cd-responsive genes that may function in drought. The expression of Cdresponsive drought-related genes was investigated under Cd exposure up to 24 h (black) and drought at 3 h (gray) in qRT-PCR analysis. The x-axis shows treatments and the y-axis shows relative expression. Transcript expression levels were normalized using an internal control (ubiquitin 1) and plotted relative to expression in non-treated samples (control) in roots (R) and shoots (S). (TIFF) Figure S4 Venn diagram analysis of Cd and other stress responsive transcripts. The resulting four-way Venn diagrams for roots and shoots show the number of transcripts responsive ($2-fold or # 0.5-fold) to Cd (24 h), drought, high-salinity and low temperature relative to the control (0 d). (TIFF) Figure S5 Expression analysis of gene families that may function in defense against Cd stress. The graph shows the expression of ROS-scavenging enzyme genes and respiratory burst oxidase homolog (Rboh) genes under Cd exposure in RNA-Seq analysis. The x-axis shows genes and y-axis shows relative expression. The black bar shows the relative expression in roots at 1 h, the red bar roots at 24 h, the light blue bar shoots at 1 h and the blue bar shoots at 24 h. (TIFF) Figure S6 Expression patterns of Cd-upregulated antioxidative and detoxification enzymes in various medium conditions by qRT-PCR analysis. The expression of Cdupregulated antioxidative and detoxification enzymes was investigated in various liquid media containing different kinds of metal ions by qRT-PCR analysis. The x-axis shows treatments and the y-axis shows relative expression. Transcript expression levels were normalized using an internal control (ubiquitin 1) and plotted relative to expression in water (control) at hour 24 in roots (R) and shoots (S). The transcripts were classified into four groups based on their expression patterns. (TIFF)