Genome-Wide Transcriptome Analysis of Cadmium Stress in Rice

Rice growth is severely affected by toxic concentrations of the nonessential heavy metal cadmium (Cd). To elucidate the molecular basis of the response to Cd stress, we performed mRNA sequencing of rice following our previous study on exposure to high concentrations of Cd (Oono et al., 2014). In this study, rice plants were hydroponically treated with low concentrations of Cd and approximately 211 million sequence reads were mapped onto the IRGSP-1.0 reference rice genome sequence. Many genes, including some identified under high Cd concentration exposure in our previous study, were found to be responsive to low Cd exposure, with an average of about 11,000 transcripts from each condition. However, genes expressed constitutively across the developmental course responded only slightly to low Cd concentrations, in contrast to their clear response to high Cd concentration, which causes fatal damage to rice seedlings according to phenotypic changes. The expression of metal ion transporter genes tended to correlate with Cd concentration, suggesting the potential of the RNA-Seq strategy to reveal novel Cd-responsive transporters by analyzing gene expression under different Cd concentrations. This study could help to develop novel strategies for improving tolerance to Cd exposure in rice and other cereal crops.


Introduction
Cadmium (Cd) is a widespread heavy metal pollutant that is highly toxic to living cells. Accumulation of the nonessential metal Cd in plants is a major agricultural problem. Specifically, Cd is absorbed by the roots from the soil and transported to the shoot, negatively affecting nutrient uptake and homeostasis in plants, even in very small amounts. Many agricultural soils have become contaminated with Cd through the use of phosphate fertilizers, sludge, and irrigation water containing Cd. Cd exposure inhibits root and shoot growth and ultimately reduces yield. Furthermore, Cd accumulation in the edible parts of plants such as seed grains places humans at a risk because of its highly toxic effects on human health. Reducing the Cd concentration in plants below the maximum level indicated by the Codex Alimentarius Commission of FAO/WHO [1] is necessary to avoid negative impacts 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, which can cause damage in various ways such as reacting with DNA causing mutation, modifying protein side chains, and destroying phospholipids [2]. Various biochemical and physiological processes associated with defense systems are active in plants under Cd exposure. Many genes such as glutathione S-transferase (GST) for detoxification and cysteine-rich metallothioneins (MT) for defense against Cd toxicity respond to Cd stress in plants and might confer Cd tolerance in rice. Transporters with heavy metal binding domains are key factors for root uptake of Cd from soil and efflux pumping of Cd at the plasma membrane; however, the manner in which these genes respond to low Cd concentrations has not been well investigated in rice.
In a previous study, we investigated the gene expression of rice plants (Oryza sativa L. cv. Nipponbare) under a high Cd concentration using the RNA-Seq platform. A clear and detailed view of the transcriptomic changes triggered by Cd exposure is important to understand the gene expression network of the basal response to Cd stress. This could not be obtained from past studies using the microarray platform, but RNA-Seq can accurately quantify gene expression levels over a broad dynamic range with high resolution and sensitivity [3]. We found that drought stress signaling pathways were activated under Cd exposure through the responses of many drought-related genes [4]. Thus, the recently elucidated scaffolding mechanisms for Cd signaling pathways are complex but of great importance. In this study, we performed rice transcriptome analysis under different low Cd concentrations using the RNA-Seq platform to deepen our understanding of Cd responses.  [5]. After 10 days, seedlings of uniform size and growth were subjected to Cd stress treatment by transferring them to a similar medium with 0.2, 1, or 50 M Cd. These values were chosen based on a report that the total dissolved Cd in 64 fields with Cd-contaminated soils ranged from 0.03 to 182 g/L [6] in previous experiences. The plants were maintained under Cd stress conditions for 14 d. Root and shoot samples were collected at approximately 9:00 AM, frozen in liquid nitrogen, and stored at −80 ∘ C until subsequent analyses. Total RNA was extracted from both root and shoot samples using an RNeasy Plant Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Construction of 34 cDNA libraries (2 tissues, 4 conditions, 2 treatments, 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 as a part of establishing TENOR (Transcriptome Encyclopedia of Rice, http://tenor.dna.affrc.go.jp/) [7]. The resulting RNA-Seq data were deposited in the DDBJ Sequence Read Archive (Accession number DRA000959).

Identification of Differentially Expressed Transcripts.
The biological replicates (2-3) for each set of conditions were highly correlated (coefficient > 0.95), so reads from the same treatment were merged for subsequent analysis. Trimming of Illumina adaptor sequences and low-quality bases ( < 20) by Cutadapt [8] and mapping of preprocessed reads to the IRGSP-1.0 genome assembly (http://rapdb.dna.affrc.go.jp/) were performed as described previously [9]. To estimate the expression levels of each transcript, all preprocessed reads were mapped to the IRGSP-1.0 genome assembly by Bowtie with default parameters [10]. The expression level for each transcript was calculated as the RPKM-(Reads per Kilobase Exon Model per million mapped reads-) derived read count [11] based on the number of uniquely mapped reads that overlapped with exonic regions. A -test was performed to detect differentially expressed transcripts in the control and Cd treatments based on the statistical null hypothesis that the proportions of mapped reads to the transcripts were the same between the two conditions. 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.

Hierarchical Clustering and Gene Ontology Enrichment
Analysis. The Cd-responsive transcripts in root and shoot were used for hierarchical clustering analysis. We used the heatmap.2 in the R package gplots (version 2.11.0) to perform clustering analyses of transcripts. The scores were used to compare significant changes in gene expression. A Gene Ontology (GO) term was assigned to each transcript based on the GO annotations for biological process, molecular function, and cellular component in RAP-DB. GO enrichment was evaluated by Fisher's exact test with a FDR threshold of 5% for responsive transcripts in the biological process category of each cluster. The results were plotted as − log 10 of FDR values in a heatmap.

Low Cd Concentration Exposure of Rice Plants and Growth
Retardation during the Treatment. We used rice plants grown in hydroponic culture, which enabled us to control the Cd exposure easily. High Cd concentration exposure has been previously shown to elicit robust physiological responses and gene expression as acute toxic responses in rice seedlings [12][13][14]. Growth retardation of the shoot was slightly visible after 1 d (data not shown), the leaves turned yellow and the leaf tips of the seedlings began to wilt after 4 d, and all leaf blades were curled completely and the seedlings were wilting after 10 d under high Cd concentration (50 M) exposure ( Figure 1). While no visible symptoms were observed in shoots under low Cd concentration exposure (0.2 and 1 M Cd) after 1 d, growth retardation occurred gradually compared with the control, with symptoms starting to appear after 7 d. Plants in the same growth chamber exposed to different Cd concentrations showed clear growth differences after 10 d ( Figure 1). Even after 28 d, the seedlings under low Cd concentration exposure did not show yellow leaves or wilting (data not shown). These results suggested that high Cd concentration exposure causes fatal damage to plants while low Cd concentrations lead to growth retardation ( Figure 1), which is supported by the fact that plant detoxification processes are insufficient to cope with this toxic metal beyond a 10 M dose [15].

Gene Expression Profiles under Low Cd Concentration Exposure in Rice.
We next analyzed the transcriptome profiles of the response to Cd exposure using RNA-Seq during plant growth, at 1, 4, and 10 d after Cd treatment, and before treatment (0 d). For each set of conditions, an average of approximately 15.1 million (92.2%) quality-evaluated reads (total 211 million) were mapped to the rice genome sequence and used for further analysis (   Figure 2). Fifty-one transcripts including GST, MT, Prx (peroxidase), and heat shock proteins were upregulated more than 20fold among the upregulated transcripts in roots at 1 M Cd ( Table 1). Induction of detoxification enzymes against oxidation stress such as GST and Prx under Cd exposure might be associated with the defense system that confers Cd tolerance to plants [16][17][18] even at low Cd concentrations. The cysteine-rich MT might function as a ligand for chelation of metal ions to defend against Cd toxicity [19]. The DREB/Crepeat binding factor (CBF) specifically interacts with the DRE/CRT cis-acting element and controls the expression of many stress-inducible genes in plants [20]. The activation of gene expression in several drought stress signal pathways under Cd exposure has been reported [4]. Five heat shock proteins (Hsps) were strongly upregulated in roots under 1 M Cd, with the greatest relative expression at 1 d (Table 1). These genes may contribute to cellular homeostasis by protecting macromolecules such as enzymes, protein complexes, and membranes under Cd exposure. This result suggested that the roots of hydroponically cultured rice might be affected more directly and earlier by Cd exposure. There was a difference between the low Cd concentrations in that no Hsps were strongly upregulated in roots at 0.2 M Cd (Table 1), suggesting that the effect of this condition might be small or show time lag. In shoots, 15 and 11 transcripts were upregulated more than 20-fold among the upregulated transcripts under 0.2 and 1 M Cd, respectively (Table S2). Nine transcripts including Nramp1 (natural resistance-associated macrophage protein) were upregulated under both 0.2 and 1 M Cd (Table S2). In Arabidopsis, Nramp1 localizes to the plasma membrane and functions as a high-affinity transporter for manganese (Mn) uptake [21], while OsNramp5 uptakes Mn and Cd [22]. Transporters with heavy metal binding domains are often capable of transporting several metals, such as Fe, Zn, Mn, and Cd, because of their low substrate specificity [23][24][25][26]. We found that upregulation of a HLH DNA-binding domain containing transcription factor (Os04g0301500) in both roots and shoots peaked at 4 d under 0.2 M Cd; this protein may function as a regulatory factor under Cd exposure ( Table 1, Table S2). The number of downregulated transcripts in roots peaked at 4 d after Cd exposure, while the number in shoots gradually increased under low Cd concentration exposure (Figure 2). A few dozen transcripts were downregulated less than 0.05-fold among the downregulated transcripts in roots and shoots under Cd exposure (Table S2). Therefore, a small part of transcripts were strongly up-or downregulated among several thousand responsive transcripts under low Cd concentration exposure. Largescale changes in gene expression occurred in rice under Cd exposure, even at low concentrations, possibly because Cd is a nonessential metal for the plant.
To obtain a functional annotation of responsive transcripts under Cd exposure, we used GO biological process categories. The responsive transcripts in shoot and root were clustered into several groups based on their expression patterns. GO enrichment analysis was performed using clustered transcripts assigned by GO  were also included in both clusters. These may function in plant growth. Thus, these correspond to the observed changes in phenotype (Figure 1), which clearly validated the RNA-Seq expression profiling data obtained from rice tissue under Cd stress condition. However, the pattern of gene expression is quite complex and would require more detailed analysis.

Constitutively Expressed Genes Responded Differently under Low Cd Concentration to High Cd Concentration.
As many genes responded to both low and high Cd concentrations [4], we assessed the effect of the stress degree on rice seedlings through the expression of constitutively expressed genes. We investigated the expression of 18 genes annotated by the RAP that were expressed constitutively in 39 tissues collected throughout the life cycle of the rice plant from two varieties according to 190 Affymetrix GeneChip Rice Genome Arrays, in addition to four genes annotated by the RAP that have frequently been used as internal controls in expression analyses [27]. The results showed that the expression of more than half of them fluctuated drastically (>2 or <2) in roots or shoots after 1 d of high Cd concentration exposure (Figure 3). This drastic response may be partly because RNA-Seq can accurately quantify gene expression levels over a broad dynamic range with high resolution and sensitivity [10,28,29]. However, our results suggest that their expression is greatly affected by strong stress, even though they are expressed constitutively across the developmental course. Note that a high Cd concentration can cause fatal damage to rice seedlings, such as by affecting homeostasis, which corresponds to the observed changes in phenotype (Figures 1 and 3).

Comparative Gene Expression Analysis between Low and High Cd Concentrations Reveals Novel Cd-Responsive Transporters.
We investigated the expression of metal transporter genes containing metal ion binding Pfam domains [PF01554 (MatE), PF08370 (PDR assoc), PF01545 (Cation efflux), PF02535 (Zip), PF00403 (HMA), and PF01566 (Nramp)] that may function in Cd transport under Cd exposure. The expression of 183 transport transcripts was compared between low and high Cd concentration treatments in roots and shoots at 1 d, because Cd uptake from the hydroponic culture and efflux pumping are initial responses to Cd exposure ( Figure 4, Table S3). The transcripts tended to be more responsive in roots and shoots under higher Cd concentration exposure. This result indicated the potential of the RNA-Seq strategy to reveal novel Cd-responsive transporters by analyzing gene expression under exposure to different Cd concentrations. The responsive transcripts might function in roots at the early stage of Cd exposure. No transcripts were upregulated more than 3-fold in shoots under low Cd exposure ( Figure 4, Table S3), suggesting that the effect takes more time to appear in shoots. Os03g0667500 (Zip, root) encoding iron-regulated transporter 1 (IRT1) was upregulated more than 5-fold under low Cd concentrations but responded only slightly under the high Cd concentration. IRT1s often transport Cd because of their low substrate specificity [24][25][26]30]. Os02g0585200 (HMA, root), Os03g0152000 (HMA, root), Os0g0584800 (HMA, root), Os01g0609900 (PDR assoc, shoot), and Os01g0609300 (PDR assoc, shoot) showed the highest (32-fold) upregulation under high Cd concentration exposure and responded only slightly to low Cd concentrations (Table S3). The balance between Cd and various other metal ions in the hydroponic culture might affect the expression of these genes, because specific systems for transporting Cd may have not developed in rice as it is a nonessential metal. The effects of other ions on the expression of transporters [4] and responsive genes associated with defense systems against Cd (Supplementary Figure S2) have been indicated.

Conclusions
We generated gene expression profiles for rice seedlings grown under low Cd concentrations. Phenotypic observations and constitutive gene expression indicated that low Cd concentrations cause growth retardation but are far from being fatal in rice. Several genes associated with defense systems were strongly upregulated; the expression of metal ion transporter genes tended to correlate with Cd concentration and GO enrichment analysis of the clustered genes based on their expression patterns, suggesting that our transcriptome profiles reflect responses to Cd in rice. Our data also suggest that it could be dangerous to eat plants that do not show specific Cd pollution symptoms growing in soil contaminated by small amounts of Cd. Establishing the exact composition and organization of the transcriptional network underlying the response to Cd exposure will provide a robust tool for improving crops in the future, for example, by creating low Cd uptake plants.  [27] suggested the following genes as candidates for constitutive expression: glycine-rich RNA-binding protein (Os12g0632000), expressed protein (Os06g0686700), profilin (Os06g0152100), ADP-ribosylation factor (Os05g0489600), triosephosphate isomerase (Os01g0147900), glycine-rich RNA-binding protein (Os03g0670700), peptidyl-prolyl cis-trans isomerase (Os02g0121300), endothelial differentiation factor (Os08g0366100), ubiquitin monomer (Os06g0681400), protein translation factor SUI1 (Os07g0529800), GAPDH (Os08g0126300), polyubiquitin (Os02g0161900), protein elongation factor (Os02g0519900), translation initiation factor (Os03g0758800), ubiquitin-conjugating enzyme (Os01g0819500), GTP-binding nuclear protein (Os05g0574500), peptidyl-prolyl isomerase (Os02g0760300), and 60S ribosomal protein L31 (Os02g0717800). Their paper also introduced the following genes that have frequently been used as internal controls in expression analyses: elongation factor1-alpha (Os03g0177500), ubiquitin fusion protein (Os03g0234200), GAPDH (Os02g0601300), and tubulin beta-6 chain (Os01g0805900).