Fei He
Qingquan Liu
Li Zheng
Yaqiong Cui
Zhenguo Shen
Luqing Zheng- College of Life Sciences, Nanjing Agricultural University, Nanjing, China
Widely-spread cadmium (Cd) pollution in the soil threatens both crop production and human health. How plants deal with the excess Cd are largely unknown. To evaluate the molecular mechanism by which plants respond to Cd stress, rice seedlings were treated with two concentrations of Cd and subjected to deep RNA sequencing. Comprehensive RNA-Seq analysis of rice roots under two gradients of Cd treatment revealed 1169 Cd toxicity-responsive genes. These genes were involved in the reactive oxygen species scavenging system, stress response, cell wall formation, ion transport, and signal transduction. Nine out of 93 predicted long non-coding RNAs (lncRNAs) were detected as Cd-responsive lncRNAs due to their high correlation with the Cd stress response. In addition, we analyzed alternative splicing (AS) events under different Cd concentrations. Four hundred and seventy-six differential alternatively spliced genes with 542 aberrant splicing events were identified. GO analysis indicated that these genes were highly enriched in oxidation reduction and cellular response to chemical stimulus. Real-time qRT-PCR validation analysis strengthened the reliability of our RNA-Seq results. The results suggest that post-transcriptional AS regulation may also be involved in plant responses to high Cd stress.
Introduction
Cadmium (Cd) is one of the most toxic heavy metals to organisms. Unlike some essential heavy metals, such as copper (Cu), zinc (Zn), manganese (Mn), and nickel (Ni) that are necessary for plant growth and development, Cd is considered a non-essential metal element for plants. Widely spread Cd pollution has significantly affected human health both in terms of its direct effect on crop production and its high accumulation in the edible part of crops such as rice (Wang et al., 2015b). Previous studies have shown that excessive Cd in plants causes an accumulation of reactive oxygen species (ROS), lipid peroxidation, enzyme inactivation, and DNA and membrane damage (Hall, 2002; Boominathan and Doran, 2003), resulting in various toxicity phenotypes, such as chlorosis, wilting, growth reduction, and cell death (Sandalio et al., 2001; Rodríguez-Serrano et al., 2009). However, as signaling molecules, ROS can be produced at controlled levels and can lead to tolerance responses (Maksymiec et al., 2007; Lin and Aarts, 2012).
Plants have evolved several common mechanisms to prevent excess Cd effects, such as extracellular exudates, transport, chelation, and sequestration and repair of damaged proteins (Hall, 2002; Mendoza-Cózatl et al., 2005; Weber et al., 2006; Ahsan et al., 2009; Lin and Aarts, 2012). By inhibiting toxic metal transport, plants can reduce toxic metal influx and enhance metal removal from the cytosol (Wysocki and Tamás, 2010; Clemens et al., 2013). Several transporters involved in the acquisition, distribution and homeostasis of Cd in plants have been identified, including heavy metal ATPases (HMA), ATP-binding cassette transporters (ABC), natural resistance-associated macrophage protein (Nramp), metal transporter proteins (MTPs), and low-affinity cation transporter (LCT) (Thomine et al., 2000; Verrier et al., 2008; Yuan et al., 2012; Uraguchi et al., 2014; Slamet-Loedin et al., 2015). One recent study showed that Cd enter into rice root through OsNramp5, originally uptake transporter for essential element Mn, because of the unspecificity properties of this transporter (Sasaki et al., 2012; Yoneyama et al., 2015). The HMA family transporter OsHMA2 and OsHMA3 plays roles in translocation of Cd into shoot and sequestration of Cd into vacuole, respectively (Ueno et al., 2010; Nocito et al., 2011; Satoh-Nagasawa et al., 2012; Takahashi et al., 2012). In mature stage, OsLCT1 regulates cadmium transport into rice grains (Uraguchi et al., 2011, 2014). In Arabidopsis, the ABC transporter AtPDR8 is a cadmium extrusion pump conferring heavy metal resistance (Kim et al., 2007), whereas the phytochelatin transporters AtABCC1 and AtABCC2 mediate tolerance to cadmium and mercury (Park et al., 2012). Cadmium-inducible expression of the ABC-type transporter AtABCC3 increases phytochelatin-mediated cadmium tolerance in Arabidopsis (Brunetti et al., 2015). Glutathione (GSH) and phytochelatins (PCs) chelating represent another detoxification strategy in plant cells (Cobbett and Goldsbrough, 2002).
Microarray data provides highly useful information on the plant responses to Cd stress at the transcriptomic level. Genes related to Cd response have been identified by various microarray analyses in plants, including transcriptomic analysis of the response to Cd stress in Arabidopsis roots, comparative transcriptomic analysis of Cd-treated roots of A. thaliana and the Cd-hypertolerant metallophyte Arabidopsis halleri, a time course analysis of gene regulation under Cd stress in rice shoots and roots, a microarray-based analysis of cadmium-responsive microRNAs in rice, and an analysis of early transcriptomic responses to Cd in rice roots (Herbette et al., 2006; Weber et al., 2006; Ogawa et al., 2009; Zhao et al., 2009; Ding et al., 2013; Lin et al., 2013). These transcriptomic analyses using microarray have given rise to a view of Cd-responsive gene expression in plants. However, compared with the high-throughput sequencing strategies, the microarray based transcriptomic strategy may lead to the loss of some important candidate genes during screening (Weber et al., 2006). At present, still a few transcriptomic studies focusing on genes involved in the plant response to Cd have been conducted using high-throughput sequencing approaches (Tang et al., 2014; Peng et al., 2015; Xu et al., 2015). While some candidate transcription factors (TFs) involved in Cd tolerance have been identified using chip strategy (DalCorso et al., 2010), the regulatory mechanisms involved are still largely unknown. Expression of several genes change significantly under Cd stress, and the identification of TFs including MYB, bZIP, ethylene-responsive factor (ERF), and WRKY by transcriptomic studies in several plants, including Arabidopsis, rice and barley, suggest that transcriptional regulation plays a pivotal role in the plant response to heavy metal stress (DalCorso et al., 2010). However, the alternative splicing regulatory role as a crucial post-transcriptional mechanism in response to Cd stress is less understood.
The expression of genes is highly regulated at both transcriptional and post-transcriptional level. Post-transcriptional regulations of gene expression occur at the levels of pre-messenger RNA (mRNA) processing (capping, splicing, and polyadenylation), mRNA stability, and mRNA translation (Floris et al., 2009). Alternative splicing (AS) is a typical post-transcriptional regulation. Recent studies indicate that AS play an important role in plant in response to abiotic stresses. The regulatory role of pre-mRNA splicing in temperature (Chang et al., 2014; Capovilla et al., 2015; Filichkin et al., 2015a), nutrient deficiency (Li et al., 2013) and salt (Ding et al., 2014; Feng et al., 2015) stresses has been reported.
High-throughput sequencing strategies have provided a new perspective on unknown transcript dynamics under specific stresses. In this study, we used RNA-Seq to understand the mechanisms of Cd toxicity, cellular detoxification and protection pathways in response to Cd in rice roots. Here, we found candidate genes associated with the reactive oxygen species (ROS) scavenging system, stress related proteins, transporters, TFs, signaling transduction, cell wall metabolism, and long non-coding RNAs (lncRNAs).
The results of our study provide a new overview of transcriptomic responses to Cd stress. These responsive genes may contribute to a reduction in Cd toxicity and rice tolerance to Cd stress. Furthermore, we examined the expression patterns of lncRNAs and revealed the genome-wide dynamic changes in alternative splicing in response to Cd stress. The transcriptional and post-transcriptional alterations in rice roots under Cd stress provide a more comprehensive understanding of the plant responses to Cd stress.
Materials and Methods
Plant Materials, Growth Conditions, and Treatments
Rice (Oryza sativa spp. japonica cv. Nipponbare) seeds were sown on mesh floating in a 0.5 mM calcium chloride (CaCl2) solution and maintained for 2 days at 25–30°C in the dark, thereby inducing germination. Seedlings were transferred into Kimura B nutrient solution containing the macronutrients (mM): (NH4)2SO4 (0.18), MgSO4·7H2O (0.27), KNO3 (0.09), Ca(NO3)2·4H2O (0.18), and KH2PO4 (0.09); and the micronutrients (μM): MnCl2·4H2O (0.5), H3BO3 (3), (NH4)6Mo7O24·4H2O (1), ZnSO4·7H2O (0.4), Fe-EDTA (20), and 0.2 μM CuSO4·5H2O (Zheng et al., 2012). The pH of the nutrient solution was adjusted to 5.5, and held under normal greenhouse conditions with illumination provided by cool-white fluorescent lamps. Growth conditions were as follows: 27/24°C day/night temperatures, 60–80% relative humidity, and a 14/10-h day/night photoperiod. Fifteen-day-old rice seedlings were treated with or without 10 and 100 μM solutions of Cd (II) hydrochloride (CdCl2) for 24 h. Following Cd treatment, roots were harvested for RNA extraction and subsequent analysis. Samples were stored at −80°C if not immediately used for RNA isolation. All experiments were performed at least twice with three biological replicates each, and representative results of one experiment are shown.
RNA Isolation, RNA-Seq Library Preparation and Sequencing
Total RNA for RNA-Seq was extracted from roots using a plant RNA kit (Omega, USA). Purified RNA was analyzed either using a ND-8000 spectrophotometer (Nanodrop Technologies, Inc., Wilmington, DE, USA), by agarose gel electrophoresis, or using a 2100-Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) to determine the RNA quantity. Those RNA samples with no smear seen on agarose gels, a 260/280 ratio above 2.0, and a RNA integrity number greater than 8.0 were used.
For RNA-Seq analysis, we mixed three replication samples for each treatment into one, and total RNA samples were then sent to Genergy Biotechnology Corporation (http://www.genenergy.cn/) for sequencing. The TruSeq RNA sample preparation kit was used for mRNA purification, fragmentation, and first- and second-strand cDNA synthesis. Double-stranded cDNAs were then purified for end repair, dA tailing, adaptor ligation, and DNA fragment enrichment. The libraries were sequenced as 101-bp paired-end reads using Illumina Hiseq2500 according to the manufacturer's instructions. Illumina reads of all samples had been submitted to the Sequence Read Archive at the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/sra) under accession number SRP058434.
Validation of Gene Expression
To validate the RNA-Seq results, the expression of selected up- or down-regulated genes, long non-coding RNA genes and alternative splicing events were confirmed by quantitative reverse transcription–polymerase chain reaction (qRT-PCR) or semi-quantitative RT-PCR analysis. After 15-day-old rice seedlings were treated with or without 0.1, 1, 10, and 100 μM CdCl2 solutions of for 24 h, the roots samples were harvested for RNA extraction. After RNA extraction, genomic DNA was removed with the Rnase Free Dnase Set (Omega, USA) following the manufacturer's instructions. RNA was reverse-transcribed using a Super Reverse transcription kit (BioTeke, China). The quantitative real-time PCR was performed on a 7500 PCR system (Applied Biosystems, USA) to confirm the RNA-seq results with the primer sets shown in Supplemental Tables S1–S3. The PCR protocol was as follows: initial denaturation at 95°C for 30 s, followed by 95°C for 5 s, and 60°C for 34 s in a 40-cycle reaction, then followed by dissociation stage. ACTIN was used as an internal standard (Liu et al., 2015). 2 −ΔΔCt method calculation was used for data analysis. Moreover, we also performed qRT-PCR validations on two genes showing no expression changes upon Cd stress (Supplemental Table S31).
Mapping Reads to the Oryza sativa Genome with Combined Gene Annotation
Complete known gene/transcript annotation information for Oryza sativa was integrated from RAPDB (http://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_representative_2015-03-31.tar.gz,IRGSP-1.0_predicted_2015-03-31.tar.gz) and MSU 7.0 (ftp://ftp.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.gff3) databases using Cuffcompare method in Cufflinks (Trapnell et al., 2010; Kawahara et al., 2013; Sakai et al., 2013). For each sample, raw paired-end reads were quality- and adaptor-trimmed using Trim Galore (-q 25, -length 30). Next, clean reads were aligned to the Os-Nipponbare-Reference-IRGSP-1.0 reference genome (http://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_genome.fasta.gz) using TopHat2 with the “-G” option of the above merged reference gene GTF file (Trapnell et al., 2009).
Long Non-Coding RNA Prediction
“Cufflinks-Cuffmerge-Cuffcompare” pipeline was used to identify novel transcripts. For each sample, mapped reads were assembled using the Reference Annotation Based Transcript assembly strategy from Cufflinks utilizing the “-g” option of the above merged reference gene GTF file. Three sets of separately assembled transcripts were merged together into the final set of assembled transcripts using Cuffmerge. All transcripts were then compared to all known gene annotations of RAPDB and MSU using Cuffcompare. Unknown intergenic transcripts, intron transcripts, antisense exon transcripts overlapping within reference exons and antisense intron transcripts overlapping within reference introns were characterized as preliminary candidate lncRNA transcripts. Novel transcripts, which were single-exon or short, could be false positives due to sequencing or alignment errors. In order to effectively distinguish lncRNAs from protein-coding transcripts, two software tools named CPC and CNCI were implemented to predict plant lncRNAs. CPC (Coding Potential Calculator: http://cpc.cbi.pku.edu.cn/) was a support vector machine-based classifier which assessed the protein-coding potential of a transcript based on quality, completeness, and sequence similarity of its open reading frame to known proteins in NCBI NR database (Kong et al., 2007). CNCI (Coding-Non-Coding Index: http://www.bioinfo.org/software/cnci) was also a support vector machine-based classifier which effectively distinguished protein-coding and non-coding transcripts based on the usage frequency of adjoining nucleotide triplets between coding and non-coding region sequences (Sun et al., 2013). For either coding/non-coding classifier, the smaller than zero its predicted score was, the more reliable prediction of lncRNA transcript was. In order to reduce false positives, candidate lncRNAs were predicted by satisfying the following criteria: Exon number ≥ 2, RNA length ≥ 200 nt, CPC score < 0, and CNCI score < 0.
Differential Expression Analysis
Raw read count of each known or predicted gene was generated using HTSeq with union-count mode (Anders et al., 2014). Within the merged reference gene annotation, a small number of genes overlap, either on the same strand or on different strands. To avoid erroneous counting of sense and antisense reads, ambiguous aligned reads from overlapped regions were excluded. After normalization by estimating the dispersion factors in DESeq R package, normalized read count table was used for determining statistical significance (Anders and Huber, 2010). Differentially expressed genes exhibiting 2-fold changes and Benjamini and Hochberg-adjusted P-values (FDR) ≤ 0.05 were selected. If the gene normalized read count value was close to 0, log2 (normalized read count+1) was used to compute the fold change. In order to show expression increased or decreased trend of these differential expressed genes with the rise of Cd concentration, the log2-transformed normalized count value of each gene was standardized to the same scale with a mean of 0 and standard deviation of 1, and then a heatmap was generated using pheatmap R package. For Cd up- or down-regulated gene set enrichment analysis, agriGO was used to detect the significantly enriched GO terms compared with the genome-wide background. KEGG pathway enrichment analysis was performed using the Hypergeometric test in R. Significantly enriched GO terms and KEGG pathways were also selected by a threshold FDR (adjusted P-value) ≤ 0.05.
Differential Alternative Splicing Analysis
Alternative splicing events (AS) are classified into five broad categories including skipped exon, retained intron, alternative 5′ splice site, alternative 3′ splice site and mutually exclusive exon (Barbazuk et al., 2008). For each sample, known and novel AS event was extracted from aligned BAM file and the above merged reference gtf file using rMATS (Shen et al., 2014). rMATS also quantified AS event of which exon or intron was included by inclusion junction counts (IJC) and of which exon or intron was skipped by skipping junction counts (SJC). Expressed alternative splicing events were declared if IJC ≥ 1 and SJC ≥ 1. To identify differential alternative splicing (AS) events between control and Cd-treated conditions, rMATS evaluated whether the difference in the exon or intron inclusion level of AS event between two conditions exceeds a stringent threshold (|IncLevelDifference|≥ 0.1 and P ≤ 0.05).
Results
Alignment and Assembly of RNA-Seq Datasets
To gain additional insight into the rice transcriptomic response to environmental Cd stress, 15-day-old rice seedlings were treated with 10 or 100 μM solutions of Cd2+, or without Cd (control), for 24 h, at which point root samples were harvested and labeled as Cd+, Cd++, and control, respectively. These samples were used for 101 bp paired-end (PE) deep sequencing on an Illumina HiSeq 2500 platform. After the adaptor and low-quality sequences of pair-end reads were trimmed, in total, ~218 million clean reads (~21.5 Gb) were obtained, with an average of ~72.8 million reads (~6.9 Gb) per sample. Our workflow for analysis of RNA-Seq data is illustrated in Figure 1. An average of ~68.9 million clean reads per sample, corresponding to ~94.6% of the total clean reads per sample, were aligned to the Os-Nipponbare-Reference-IRGSP-1.0 reference genome using TopHat2. Among the mapped reads per sample, ~96.7% were uniquely aligned, and the remaining 3.3% were multiple mapped reads (Table 1). For each sample, the resulting aligned reads were then analyzed with Cufflinks, which assembles the aligned reads into transcripts. The assembled transcripts were further filtered based on greater than 0 FPKM (Fragments Per Kilobase Million) expression level. The remaining transcripts were processed separately by comparing our merged reference transcript annotation using Cuffcompare and further classified into four main categories: known reference transcript, novel isoform transcripts, antisense transcripts, and intergenic transcripts. An average 32.7% of 113,001 merged known reference transcripts per sample was expressed. Averagely, 8908 novel isoform transcripts, 157 antisense transcripts, and 754 intergenic transcripts for each sample were also detected (Table 2).
Figure 1. Flowchart of our RNA-seq data analysis workflow. High-level representation of the pipeline workflow for processing RNA-Seq data.
Table 1. Statistics of the Read Alignments in the RNA-Seq study.
Table 2. Number of transcripts detected by annotation classification reported by Cufflinks and cuffcompare.
Identification of Differentially Expressed Genes
Differential expression genes (DEGs) were identified between Cd+ or Cd++ and control samples using DESeq. A total of 1169 DEGs of Oryza sativa in response to Cd stress from the RNA-Seq data were identified. Under Cd+ treatment, 214 were up-regulated while 22 were down-regulated (Figure 2A, Supplemental Tables S4, S5). Under Cd++ treatment, 914 genes were up-regulated, and 248 were down-regulated (Figure 2A, Supplemental Tables S6, S7). These results suggest that most up- and down-regulated genes are common under the two gradients of Cd treatment. 208 out of 214 regulated genes and 21 out of 22 repressed genes regulated by Cd+ treatment were also among the up- and down-regulated genes regulated by Cd++ treatment (Figure 2B). In addition to the significantly increased number of up- and down-regulated genes after Cd++ compared with Cd+ treatment, the expression ratio also changed significantly between the two concentrations of Cd relative to the control, indicated by the Heatmap analysis (Figure 2C). We compared these DEGs with results from a previous study (Lin et al., 2013) that used an Agilent two-color Rice Oligo DNA Microarray 44 K designed for RAP-DB genes. 89 (45.4%) out of 196 Cd+ up-regulated genes and 294 (35.2%) out of Cd++ 836 up-regulated genes identified in our RNA-Seq assay were also found in the microarray datasets, while in contrast, most of the down-regulated genes varied between the two assays (Figure 2D). In addition to the differences in culture conditions and the rice cultivars we used, the differences between two-color microarray and RNA-Seq technologies may also lead to the discrepancy on some DEGs results. Real-time qRT-PCR analysis was also performed to validate some up- and down-regulated genes expression data according to the RNA-Seq results (Figure 3A).
Figure 2. Differentially expressed genes of Oryza sativa in response to Cd stress from RNA-seq data. (A) Summary of significant up- and down-regulated genes between two Cd concentration gradient treated roots. (B) Venn diagram analysis showed a total of 1169 differentially expressed genes in different conditions. (C) Heatmap showed expression increased or decreased trend of total differential expressed genes with the rise of Cd concentration after scaling each gene to a mean of 0 and standard deviation of 1. (D) Overlap of differential expressed RAPDB genes from RNA-Seq and Lin et al microarray study. Total number of RAPDB genes was 45,990. The P-value of overlap significant statistical test was 2.2E-16 and 2.2E-16 by Fisher exact test, respectively.
Figure 3. Validation of RNA-Seq data by real-time qRT-PCR or RT-PCR. (A) Validation of several differentially expressed Cd up- or down-regulated genes by qRT-PCR. (B) qRT-PCR detection of long non-coding RNAs expression under Cd stress. (C) Semi-quantitative PCR detection of alternative splicing event analysis under Cd stress. Rice seedling roots were collected 24 h after 10 μM Cd (Cd+) exposure in (A,C), or treated with 0.1, 1, 10, and 100 μM Cd for (B), see details in the Materials and Methods.
To identify possible biological processes or pathways that were altered in Cd stress treatment and between the two levels of Cd stress, Gene Ontology (GO) and KEGG pathway functional enrichment was performed using an FDR adjusted p ≤ 0.05 as the cutoff. Our analysis revealed that 208 commonly up-regulated DEGs were highly enriched in oxidation reduction, response to stimulus, response to stress, response to oxidative stress, phenylalanine metabolism, phenylpropanoid biosynthesis, protein processing in endoplasmic reticulum, cutin, suberine and wax biosynthesis and glutathione metabolism. In contrast, there were no significantly enriched biological processes and pathways in 21 commonly down-regulated DEGs. Moreover, 706 Cd++ specific up-regulated DEGs were highly enriched in oxidation reduction, response to biotic stimulus, response to stimulus, response to stress, cell wall organization or biogenesis, oxidation reduction, secondary metabolic process, lipid metabolism, regulation of transcription, phenylalanine, tyrosine and tryptophan biosynthesis, phenylpropanoid biosynthesis, alpha-Linolenic acid metabolism, and diterpenoid biosynthesis. Two hundred and twenty-seven Cd++ specific down-regulated DEGs were highly enriched in oxidation reduction, response to stress, lipid transport, lipid localization, response to stimulus, response to oxidative stress, response to chemical stimulus, phenylalanine metabolism and phenylpropanoid biosynthesis. It was noteworthy that the common significantly enriched GO terms and KEGG pathways identified by separate enrichment analysis strategy for up- and down-regulated genes were oxidation reduction (37 up- and 11 down-regulated genes), response to stimulus (66 up- and 19 down-regulated genes), response to stress (55 up- and 19 down-regulated genes), response to oxidative stress (12 up- and 8 down-regulated genes), phenylalanine metabolism (17 up- and 5 down-regulated genes), and phenylpropanoid biosynthesis (21 up- and 7 down-regulated genes). Within these functional GO terms and KEGG pathways, more genes were found to be induced than repressed, particularly, 8 up- and 5 down-regulated peroxidase genes were in common. Globally, the above results indicated that the gradual increase of Cd concentration had significant effects on cell wall biogenesis, lipid metabolism and regulation of transcription (Figure 4, Supplemental Tables S8–S11).
Figure 4. Function enrichment analysis of differentially expressed genes. (A) Gene Ontology (GO) enrichment analysis for the up-regulated genes. Only the top false discovery rate (FDR) ranked 10 enrichment of GO terms from “biological process” category were listed. (B) GO enrichment analysis for the down-regulated genes. Only “biological process” category were listed. (C) KEGG pathway analysis for the up-regulated genes. (D) KEGG pathway analysis for the down-regulated genes.
Cd Stress Regulated Transcription Factors
Given that TFs appear to have a major effect on the network of Cd-responsive genes, one objective of our work was to identify Cd regulated TFs. In total, 110 differentially expressed TF genes were identified following Cd treatment (Figure 5A, Supplemental Table S12). 20 up-regulated genes and 1 down-regulated gene were characterized as zinc finger TFs. 21 up-regulated genes were characterized as AP2-EREBP TFs. 15 up-regulated genes and 1 down-regulated gene were characterized as WRKY TFs. 14 up-regulated genes were characterized as NAC TFs. 12 up-regulated genes and 1 down-regulated gene were characterized as MYB TFs (Figure 5B). Moreover, a relatively small group of TIFY, HSF, bHLH, HB, GRAS, GNAT, LOB, MBF1, and Orphans TF family members were also present in the Cd-regulated TF family (Figure 5).
Figure 5. Differentially expressed transcription factors. (A) Summary of significant up- and down-regulated transcription factors under two different levels of Cd treatment compared with control in rice roots. (B) Gene expression pattern of five transcription factor classes, which contained most DGEs.
Identification of Long Non-Coding RNAs
A computational prediction pipeline of lncRNAs based on RNA-Seq data was performed (see Materials and Methods). In total, 93 novel genes corresponding to 122 transcripts were predicted as lncRNAs (Supplemental Table S13). According to their genomic location, 29 antisense lncRNAs were transcribed from the antisense strand of protein-coding genes, of which 13 were close to the 5′-untranslated region (UTR), 14 overlapped with middle coding exon region, only 2 were close to the 3′ UTR. 8 bidirectional lncRNAs were transcribed head-to-head with a protein-coding gene within 1 kb. The remaining 56 intergenic lncRNAs were located at least 1 kb away from a protein-coding gene. A total of 9 differentially expressed genes were detected between Cd-treated and control samples (Table 3). Among them, 5 were commonly induced and 4 were commonly repressed by Cd. Typically, read coverage signal maps of 2 up-regulated and 2 down-regulated lncRNAs were shown (Figures 6A,B). Real-time qRT-PCR validation performed in four of them confirmed their Cd-regulated expression pattern (Figure 3B).
Table 3. Differentially expressed lncRNAs between Cd-treated and control samples.
Figure 6. Read coverage signal maps and exon-intron structure of long non-coding RNAs. (A) Two up-regulated lncRNAs: XLOC_011965 and XLOC_054416. (B) Two down-regulated lncRNAs: XLOC_001126 and XLOC_048220.
Cd Regulates Alternative Splicing
To investigate the influence of Cd on alternative splicing regulation, expressed alternative splice events including skipped exon (SE), retained intron (RI), alternative 5′ splice site (A5SS), alternative 3′ splice site (A3SS), and mutual exclusive exons (MXE) in Control, Cd+ and Cd++ were identified using rMATS. A total of 8045 AS events were identified, of which 3488 AS events were commonly detected among three samples. RI, SE, and A3SS appear to be the most predominant AS events, relative to other AS event types, the most detected novel expressed SE events indicated that exon skipping greatly contributed to expressed transcript isoform diversity (Figures 7A–C). The 3488 AS events were distributed in 2415 genes. These AS genes were highly enriched in oxidation reduction, developmental process, response to hormone stimulus, signal transduction, and cell cycle (Supplemental Table S14). Using cutoff ≥ 0.1 for splicing differences and P ≤ 0.05, differential AS events were identified between Cd-treated and control samples (Figures 8, 9, Supplemental Tables S15–S24). We found 476 differential AS genes with 542 aberrant splicing events (Supplemental Table S25). Interestingly, only 8 genes were overlapped with DGEs (Supplemental Table S26). One example (Os04g0118900) corresponding to skipped exon and another (Os09g0417800) corresponding to retained intron with strong evidence of differential splicing based on their read coverage signal maps and exon-exon junction reads are shown in Figures 10, 11. Semi-quantitative RT-PCR validation performed in two selected transcripts confirmed their Cd-regulated AS pattern (Figure 3C). GO enrichment analysis indicated that 476 differential AS genes were highly enriched in oxidation reduction and cellular response to chemical stimulus (Supplemental Table S27). Two hundred and twenty-seven genes with differential RI were highly enriched in oxidation and signaling process (Supplemental Table S28). One hundred and twenty-one genes with differential SE had no significantly enriched biological processes (Supplemental Table S29).
Figure 7. Expressed alternative splicing events from RNA-seq data and comparsion analysis between Cd-treated and control. (A) Statistics of detected known and novel expressed alternative splicing events in control, Cd+, Cd++. (B) Overlap of expressed AS events between Cd-treated and control. (C) Overlap of AS genes between Cd-treated and control.
Figure 8. Differential alternative splicing events. (A) Summary of AS events with significant increased or decreased exon inclusive ratio between Cd-treated and control. (B) Venn diagram analysis showed shared 73 differentially expressed AS events in different conditions.
Figure 9. Heatmaps of exon inclusive ratio between Cd-treated and control for differential alternative splicing events. (A) Skipped exon. (B) Retained intron. (C) Alternative 3′ splice site. (D) Alternative 5′ slice site. (E) Mutually exclusive exon.
Figure 10. One differential alternative skipped exon event of Os04g0118900. The arrow showed that the middle exon was differentially spliced between Cd-treated and control. Three numbers showed how many exon inclusive reads supported the upstream splice junction, the alternative exon itself, and the downstream splice junction and how many exon skipped reads supported the skipping splice junction that directly connects the upstream exon to the downstream exon.
Figure 11. One differential alternative retained intron event of Os09g0417800. The arrow showed that the intron was differential spliced between Cd-treated and control. Three number showed how many reads supported the middle intron between the upstream exon to the downstream exon was not retained.
Discussion
Cd is one of the most toxic heavy metals; its toxicity in plants, animals and humans causes a major problem for crop production and food safety. Therefore, it is necessary to minimize the Cd concentration in edible portions of some crops, such as rice, one of the most important crops. Understanding Cd-responsive gene regulatory mechanisms can help us identify new ways to manipulate Cd accumulation and distribution in rice. High-throughput RNA sequencing and comprehensive transcriptomic analyses help us understand gene regulation and mechanisms for plant tolerance to heavy metals.
Although some studies have shown several genes involved in the response to Cd stress in rice (Ogawa et al., 2009; Zhang et al., 2012; Lin et al., 2013), a few RNA-Seq approaches have been used to investigate heavy metal stress at the whole transcriptome level. As an open system, RNA-Seq may provide an opportunity to identify novel lncRNAs and AS genes involved in the Cd stress response. Indeed, in this study, a large number of DEGs and differential AS genes was obtained, which will provide valuable information on transcriptional and post-transcriptional regulation of gene responses to Cd stress in rice.
Common Cadmium Stress Responses
Among the Cd DEGs, the overrepresented biological functional pathway genes are involved in cell wall formation, signaling transduction, ROS scavenging system etc. (Supplemental Tables S8–S11), that is consistent with previous studies (Mendoza-Cózatl et al., 2005; Ahsan et al., 2009; DalCorso et al., 2010; Lin and Aarts, 2012; Lin et al., 2013).
Cd can induce plant cell wall modifications, and the cell wall participated in preventing Cd transportation (Meyer et al., 2015; Wang et al., 2015a). Our data showed that four cell wall related-genes encoding sperm-coating glycoprotein (SCP)-like extracellular protein were up-regulated in rice roots during the Cd+ or Cd++ stress response. As a member of the pathogenesis-related (PR) protein family, SCP-like extracellular proteins are involved in cellular signaling processes with a putative Ca2+-chelating function (Fernández et al., 1997; Milne et al., 2003) and are induced during incompatible host-pathogen interactions or other stress conditions (Van Loon and Van Strien, 1999). They may play important functions in adaptation to biotic and abiotic stresses, including Cd stress. Four calmodulin encoding genes and two calmodulin-binding genes are up-regulated under Cd exposure suggesting they may directly involve in a complex network of signaling transduction that assists rice to effectively cope with Cd stress. In addition, one MAPKK gene (Os06g0191500) was detected in rice roots and was induced by Cd stress, suggest it may also be involve in signaling transduction pathway.
The ROS machine are usually activated by Cd stress, it is really the most overrepresented biological pathway present in the DEGs of our study (Supplemental Table S30). We identified 17 genes encoding GSTs, 2 genes encoding phenylalanine ammonialyase, 2 MT genes and 3 oxalate oxidase-like protein encoding genes were up-regulated in rice under Cd stress, whereas 4 Cysteine-rich polypeptide genes were down-regulated to avoid excessive consumption of cysteine, an ingredient of GSH (Supplemental Table S30). Separate enrichment analysis of GO term and KEGG pathway for up- and down-regulated genes indicated that peroxidase family genes had different expression patterns in response to Cd stress. Interestingly, 5 genes encoding class III peroxidase genes were significantly differential expressed, including 2 up-regulated (Os12g0111800, Os10g0107000) and 3 down-regulated (Os01g0294500, Os06g0490400, Os09g0323900). Similar phenomenon was also observed in one recent study, which showed that differential expression patterns of five candidate maize class III peroxidase genes under abiotic stress (Wang et al., 2015c). High Cd concentrations lead to an increase in the levels of oxygen-free radicals and glutathione peroxidase, which catalyzes the production of GSSG, thereby reducing the content of oxygen-free radicals. ROS scavenging enzymes or antioxidants, such as peroxidases and GST, can keep ROS levels under control and reduce their toxicity to plant growth.
Under stress conditions, some new proteins, such as heat shock proteins (Hsps), which are responsible for protein folding, assembly, translocation, and degradation in many cellular processes, can be induced for adaptation purposes. Hsps play a crucial role in plant tolerance to stress (Wang et al., 2004). Here, we detected 20 Hsp genes (Supplemental Table S30), most of which were up-regulated by Cd stress in rice roots. Similar results were found in rice roots of 10 and 25 μM Cd-treated and 5 μM Cu-treated samples (Ogawa et al., 2009; Lin et al., 2013), indicating that Hsps are strongly induced by heavy metal stress. Two Chitinase (Os06g0726200, Os10g0542900) was also induced, in addition to other damage-induced proteins reported previously.
Transcriptional Regulation of Cd
TFs play many roles in plant abiotic and biotic stress responses. TFs, such as WRKY, basic leucine zipper (bZIP), ethylene-responsive (ERF), and myeloblastosis protein (MYB), play a key role in controlling the expression of specific stress-related genes in response to Cd stress (DalCorso et al., 2010). WRKY, NAC, MYB, and AP2 genes were up-regulated in rice roots treated with 10 μM Cd for 3 h (Ogawa et al., 2009). In Arabidopsis thaliana roots, some genes encoding AP2 domain-containing protein, C2H2 zinc finger protein, HSF, and MYB TFs are specifically induced by Cd2+ but not Cu2+ (Weber et al., 2006). The expression analysis of six poplar C2H2 ZFPs showed their involvement in response to cold, salt, osmotic, and mechanical stresses (Gourcilleau et al., 2011). Wheat ZAT7 and ZAT12 belonging to C2H2 zinc finger protein (ZFP) family were shown to improve ROS tolerance in Arabidopsis and may improve Al-induced oxidative stress tolerance in wheat (Ali-Benali et al., 2013). An AP2/ERF transcription factor, RAP2.6, participates in ABA, salt and osmotic stress responses in Arabidopsis (Zhu et al., 2010). A stress-responsive NAC gene isolated from rice, SNAC2, can improve stress tolerance in rice (Hu et al., 2008). Wheat overexpressing a TaNAC6 gene exhibit higher tolerance to stress (Xue et al., 2011). MYB genes can also respond to one or more stress treatments, such as high salt, PEG and exogenous ABA (Zhang et al., 2012). A number of WRKY genes are also activated by heavy metal lead stress in Arabidopsis (Liu et al., 2009). WRKY38 from barley is involved in cold and drought stress responses (Marè et al., 2004). TcWRKY53 from Cd-treated Thlaspi caerulescens is strongly induced by NaCl, drought, cold, and salicylic acid (SA) stresses (Wei et al., 2009). Five up-regulated WRKY genes including OsWRKY71 (Os02g0181300), OsWRKY28 (Os06g0649000), OsWRKY22 (Os01g0820400), and OsWRKY42 (Os02g0462800) were involved in plant defense response (Liu et al., 2007; Abbruscato et al., 2012; Chujo et al., 2013; Han et al., 2014). One up-regulated OsWRKY72 (Os11g0490900) was involved in the abscisic acid signal and auxin transport pathway (Yu et al., 2010). These results indicate that TFs may be directly involved in abiotic stress responses. Currently, we are performing functional analyses of several TFs for their putative roles in heavy metal tolerance.
Long Non-Coding RNAs
Long non-coding RNAs (lncRNAs) are a class of RNAs of more than 200 nt that do not encode proteins. The roles of lncRNAs in gene activation are poorly understood; some studies found that lncRNAs were key regulators in the cell and are involved in plant stress responses (Wilusz et al., 2009; Xin et al., 2011). In this study, the 9 lncRNA genes expressed showed good correlation with Cd concentration (Table 3). It has been determined that one type of lncRNA transcribed by RNA polymerase V in Arabidopsis can provide binding sites associated with other factors and regulate chromatin modification (Rowley et al., 2013). Using computational analysis and an experimental approach, 125 putative lncRNAs involved in the stress response in wheat and non-conserved in plants were identified (Xin et al., 2011). From these results, we can speculate that lncRNAs play an important role in plant resistance to various environment stresses. However, the biological functions of these lncRNAs have yet to be determined. Moreover, some of the lncRNAs are antisense to protein-coding genes in the 5′- or 3′-untranslated region (UTR); therefore, it is more likely that they function by regulating the expression of other genes, as well as processing into microRNAs to suppress downstream genes (Ariel et al., 2015). Therefore, the identification of Cd-regulated lncRNAs suggests that post-transcriptional regulation may also be involved in heavy metal stress responses.
Alternative Splicing
Alternative splicing (AS) is a vital regulatory mechanism used to increase transcriptome and proteome diversity and represents adaptation during evolution and stress (Staiger and Brown, 2013; Filichkin et al., 2015b). There are more examples of the positive effect and biological function of AS during stress. To our knowledge, there is limited information on the effects of heavy metal stress on AS regulation, especially in the case of global transcriptomic profiles induced by treatment with certain heavy metals. Our data presented here emphasize the importance of AS in the response to stress adaptation, including that to heavy metals. One recent study showed that Zn regulation of AS of one major facilitator superfamily (MFS) transporter is essential for Zn toxicity (Remy et al., 2014). Under excess Zn concentrations, intron retention of ZIF2 in the 5′-UTR resulted in a more stable transcript and led to increased production of proteins that function in Zn detoxification. This example is highly indicative of the involvement of functional AS in heavy metal stress responses. In the present study, three AS-regulated splicing factors under Cd stress (Os09g0491756, Os02g0146400, Os04g0118900; Supplemental Table S26), of which the typical arginine/serine-rich (SR) splicing factor gene (Os04g0118900) with significantly increased intron inclusive level indicated that several AS self-regulatory circuits of splicing regulatory factors play a unexpected role in in plant Cd stress response (Figure 10). It would be interesting to investigate the consequence on genome-wide gene expression using genetic material of the SR gene in future study. We also found that one WRKY transcription factor gene (Os09g0417800) was differential alternative spliced with significantly decreased exon inclusive level under different levels of Cd stress (Figure 11). We can speculate that alternative splicing of the WRKY transcription factor may lead to target regulation of downstream genes that is close associated with plant responses to Cd stress.
Conclusions
Recent studies on the transcriptomes of plants under Cd stress using microarray analysis have identified Cd stress-responsible metabolic pathways, including ROS scavenging enzymes, heavy metal chelators and transporters. Our present study utilizing RNA-Seq will provide a more comprehensive analysis of transcriptome under Cd stress, and will show that, in addition to transcriptional regulation, post-transcriptional regulation may also play a vital role in plant adaptation to environmental Cd stress. Moreover, our findings will not only provide us a more complete picture of the regulatory network of the Cd stress response in plants, but will also provide us with the potential to produce Cd stress-tolerant crops.
Author Contributions
Manuscript draft: FH, QL, and LZ; Analyzing data: FH and LZ; Experiment: QL, LZ, and YC; Conception and supervision of the research: ZS and LZ.
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.
The reviewer Yingchun Xu declares that, despite being affiliated with the same institute as the authors of the manuscript, the review process was handled objectively.
Acknowledgments
This work was supported by research grants from the Project of the National Natural Science Foundation of China (No.41401351 and 31172021), “the Fundamental Research Funds for the Central Universities (KYTZ201402, KYRC201302)” and the Natural Science Foundation of Jiangsu Province (BK20140682). We thank Bin Bai (Genergy Biotechnology Corporation) for excellent technical help in RNA-seq experiment.
Supplementary Material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2015.01136
Supplemental Tables S1–S3. Primers used for qRT-PCR analysis of Cd up- or down-regulated genes. Primers used for Semi-quantitative PCR detection of alternative splicing event analysis. Primers used for qRT-PCR detection of lncRNAs expression.
Supplemental Tables S4–S5. Up- and down-regulated DEGs by Cd+.
Supplemental Tables S6–S7. Up- and down-regulated DEGs by Cd++.
Supplemental Tables S8–S11. Gene Ontology and KEGG Pathway enrichment analysis results for commonly up-regulated or down-regulated DEGs between Cd+ and Cd++. Gene Ontology and KEGG Pathway enrichment analysis results for Cd++ specific up-regulated or down-regulated DEGs.
Supplemental Table S12. Transcription factors regulated by Cd.
Supplemental Table S13. Predicted lncRNAs using RNAseq data.
Supplemental Table S14. Gene Ontology and KEGG Pathway enrichment analysis results for detected AS genes.
Supplemental Tables S15–S19. Differential AS events between Cd+ and control, successively including skipped exon, retained intron, alternative 5′ splice site, alternative 3′ splice site and mutually exclusive exon.
Supplemental Tables S20–S24. Differential AS events between Cd++ and control, successively including skipped exon, retained intron, alternative 5′ splice site, alternative 3′ splice site and mutually exclusive exon.
Supplemental Table S25. Functional description for differential AS genes.
Supplemental Table S26. Gene Ontology enrichment analysis results for overlapped 8 genes between differential AS genes and DEGs.
Supplemental Table S27. Gene Ontology and KEGG Pathway enrichment analysis results for differential AS genes.
Supplemental Table S28. Gene Ontology enrichment analysis results for overlapped AS genes with differential RI.
Supplemental Table S29. Gene Ontology enrichment analysis results for overlapped AS genes with differential SE.
Supplemental Table S30. Cd stress responsive genes identified by previously reports, including Glutathione S-transferases, Peroxidase, PAL, Metalloendoproteinas, Oxalate oxidase, Cysteine-rich peptide, Terpene synthase, Terpenoid synthase, Heat shock protein, Chitinase, transporter.
Supplemental Table S31. qRT-PCR validation of two genes with no changes in expression under Cd+ stress.
References
Abbruscato, P., Nepusz, T., Mizzi, L., Del Corvo, M., Morandini, P., Fumasoni, I., et al. (2012). OsWRKY22, a monocot WRKY gene, plays a role in the resistance response to blast. Mol. Plant Pathol. 13, 828–841. doi: 10.1111/J.1364-3703.2012.00795.X
Ahsan, N., Renaut, J., and Komatsu, S. (2009). Recent developments in the application of proteomics to the analysis of plant responses to heavy metals. Proteomics 9, 2602–2621. doi: 10.1002/pmic.200800935
Ali-Benali, M. A., Badawi, M., Houde, Y., and Houde, M. (2013). Identification of oxidative stress-responsive C2H2 zinc fingers associated with Al tolerance in near-isogenic wheat lines. Plant Soil. 366, 199–212. doi: 10.1007/s11104-012-1417-y
Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 11:R106. doi: 10.1186/gb-2010-11-10-r106
Anders, S., Pyl, P. T., and Huber, W. (2014). HTSeq–A Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169. doi: 10.1093/bioinformatics/btu638
Ariel, F., Romero-Barrios, N., Jéqu, T., Benhamed, M., and Crespi, M. (2015). Battles and hijacks: noncoding transcription in plants. Trends Plant Sci. 20, 362–371. doi: 10.1016/j.tplants.2015.03.003
Barbazuk, W. B., Fu, Y., and McGinnis, K. M. (2008). Genome-wide analyses of alternative splicing in plants: opportunities and challenges. Genome Res. 314, 587–1378. doi: 10.1101/gr.053678.106
Boominathan, R., and Doran, P. M. (2003). Cadmium tolerance and antioxidative defenses in hairy roots of the Cadmium hyperaccumulator, Thlaspi caerulescens. Biotechnol. Bioeng. 83, 158–167. doi: 10.1002/bit.10656
Brunetti, P., Zanella, L., De Paolis, A., Di Litta, D., Cecchetti, V., Falasca, G., et al. (2015). Cadmium-inducible expression of the ABC-type transporter AtABCC3 increases phytochelatin-mediated cadmium tolerance in Arabidopsis. J. Exp. Bot. 66, 3815–3829. doi: 10.1093/jxb/erv185
Capovilla, G., Pajoro, A., Immink, R. G., and Schmid, M. (2015). Role of alternative pre-mRNA splicing in temperature signaling. Curr. Opin. Plant Biol. 27, 97–103. doi: 10.1016/j.pbi.2015.06.016
Chang, C. Y., Lin, W. D., and Tu, S. L. (2014). Genome-wide analysis of heat-sensitive alternative splicing in physcomitrella patens. Plant Physiol. 165, 826–840. doi: 10.1104/pp.113.230540
Chujo, T., Miyamoto, K., Shimogawa, T., Shimizu, T., Otake, Y., Yokotani, N., et al. (2013). OsWRKY28, a PAMP-responsive transrepressor, negatively regulates innate immune responses in rice against rice blast fungus. Plant Mol. Biol. 82, 23–37. doi: 10.1007/s11103-013-0032-5
Clemens, S., Aarts, M. G., Thomine, S., and Verbruggen, N. (2013). Plant science: the key to preventing slow cadmium poisoning. Trends Plant Sci. 18, 92–99. doi: 10.1016/j.tplants.2012.08.003
Cobbett, C., and Goldsbrough, P. (2002). Phytochelatins and metallothioneins: roles in heavy metal detoxification and homeostasis. Annu. Rev. Plant Biol. 53, 159–182. doi: 10.1146/annurev.arplant.53.100301.135154
DalCorso, G., Farinati, S., and Furini, A. (2010). Regulatory networks of cadmium stress in plants. Plant Signal. Behav. 5, 663. doi: 10.4161/psb.5.6.11425
Ding, F., Cui, P., Wang, Z., Zhang, S., Ali, S., and Xiong, L. (2014). Genome-wide analysis of alternative splicing of pre-mRNA under salt stress in Arabidopsis. BMC Genomics. 15:431. doi: 10.1186/1471-2164-15-431
Ding, Y., Chen, Z., and Zhu, C. (2013). Microarray-based analysis of cadmium-responsive microRNAs in rice (Oryza sativa). J. Exp. Bot. 64, 4271–4287. doi: 10.1093/jxb/ert240
Feng, J. L., Li, J. J., Gao, Z. X., Lu, Y. R., Yu, J. Y., Zheng, Q., et al. (2015). SKIP Confers Osmotic Tolerance during Salt Stress by Controlling Alternative Gene Splicing in Arabidopsis. Mol. Plant. 8, 1038–1052. doi: 10.1016/j.molp.2015.01.011
Fernández, C., Szyperski, T., Bruyère, T., Ramage, P., Mösinger, E., and Wüthrich, K. (1997). NMR solution structure of the pathogenesis-related protein P14a. J. Mol. Biol. 266, 576–593. doi: 10.1006/jmbi.1996.0772
Filichkin, S. A., Cumbie, J. S., Dharmawardhana, P., Jaiswal, P., Chang, J. H., Palusa, S. G., et al. (2015a). Environmental stresses modulate abundance and timing of alternatively spliced circadian transcripts in Arabidopsis. Mol. Plant. 8, 207–227. doi: 10.1016/j.molp.2014.10.011
Filichkin, S., Priest, H. D., Megraw, M., and Mockler, T. C. (2015b). Alternative splicing in plants: directing traffic at the crossroads of adaptation and environmental stress. Curr. Opin. Plant Biol. 24, 125–135. doi: 10.1016/j.pbi.2015.02.008
Floris, M., Mahgoub, H., Lanet, E., Robaglia, C., and Menand, B. (2009). Post-transcriptional regulation of gene expression in plants. Int. J. Mol. Sci. 10, 3168–3185. doi: 10.3390/ijms10073168
Gourcilleau, D., Lenne, C., Armenise, C., Moulia, B., Julien, J. L., Bronner, G., et al. (2011). Phylogenetic study of plant Q-type C2H2 zinc finger proteins and expression analysis of poplar genes in response to osmotic, cold and mechanical stresses. DNA Res. 18, 77–92. doi: 10.1093/dnares/dsr001
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
Han, M., Kim, C. Y., Lee, J., Lee, S. K., and Jeon, J. S. (2014). OsWRKY42 represses OsMT1d and induces reactive oxygen species and leaf senescence in rice. Mol. Cells. 37, 532–539. doi: 10.14348/molcells.2014.0128
Herbette, S., Taconnat, L., Hugouvieux, V., Piette, L., Magniette, M. L., Cuine, S., et al. (2006). Genome-wide transcriptome profiling of the early cadmium response of Arabidopsis roots and shoots. Biochimie 88, 1751–1765. doi: 10.1016/j.biochi.2006.04.018
Hu, H., You, J., Fang, Y., Zhu, X., Qi, Z., and Xiong, L. (2008). Characterization of transcription factor gene SNAC2 conferring cold and salt tolerance in rice. Plant Mol. Biol. 67, 169–181. doi: 10.1007/s11103-008-9309-5
Kawahara, Y., de la Bastide, M., Hamilton, J. P., Kanamori, H., McCombie, W. R., Ouyang, S., et al. (2013). Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice 6:4. doi: 10.1186/1939-8433-6-4
Kim, D. Y., Bovet, L., Maeshima, M., Martinoia, E., and Lee, Y. (2007). The ABC transporter AtPDR8 is a cadmium extrusion pump conferring heavy metal resistance. Plant J. 50, 207–218. doi: 10.1111/j.1365-313X.2007.03044.x
Kong, L., Zhang, Y., Ye, Z. Q., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi: 10.1093/nar/gkm391
Li, W. F., Lin, W. D., Ray, P., Lan, P., and Schmidt, W. (2013). Genome-wide detection of condition-sensitive alternative splicing in Arabidopsis roots. Plant Physiol. 162, 1750–1763. doi: 10.1104/pp.113.217778
Lin, C. Y., Trinh, N. N., Fu, S. F., Hsiung, Y. C., Chia, L. C., Lin, C. W., et al. (2013). Comparison of early transcriptome responses to copper and cadmium in rice roots. Plant Mol. Biol. 81, 507–522. doi: 10.1007/s11103-013-0020-9
Lin, Y. F., and Aarts, M. G. (2012). The molecular mechanism of zinc and cadmium stress response in plants. Cell. Mol. Life Sci. 69, 3187–3206. doi: 10.1007/s00018-012-1089-z
Liu, Q., Zheng, L., He, F., Zhao, F. J., Shen, Z., and Zheng, L. (2015). Transcriptional and physiological analyses identify a regulatory role for hydrogen peroxide in the lignin biosynthesis of copper-stressed rice roots. Plant Soil. 387, 323–336. doi: 10.1007/s11104-014-2290-7
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
Liu, X., Bai, X., Wang, X., and Chu, C. (2007). OsWRKY71, a rice transcription factor, is involved in rice defense response. J. Plant Physiol. 164, 969–979. doi: 10.1016/j.jplph.2006.07.006
Maksymiec, W., Wójcik, M., and Krupa, Z. (2007). Variation in oxidative stress and photochemical activity in Arabidopsis thaliana leaves subjected to cadmium and excess copper in the presence or absence of jasmonate and ascorbate. Chemosphere 66, 421–427. doi: 10.1016/j.chemosphere.2006.06.025
Marè, C., Mazzucotelli, E., Crosatti, C., Francia, E., and Cattivelli, L. (2004). Hv-WRKY38: a new transcription factor involved in cold-and drought-response in barley. Plant Mol. Biol. 55, 399–416. doi: 10.1007/s11103-004-0906-7
Mendoza-Cózatl, D., Loza-Tavera, H., Hernández-Navarro, A., and Moreno-Sánchez, R. (2005). Sulfur assimilation and glutathione metabolism under cadmium stress in yeast, protists and plants. FEMS Microbiol. Rev. 29, 653–671. doi: 10.1016/j.femsre.2004.09.004
Meyer, C. L., Juraniec, M., Huguet, S., Chaves-Rodriguez, E., Salis, P., Isaure, M. P., et al. (2015). Intraspecific variability of cadmium tolerance and accumulation, and cadmium-induced cell wall modifications in the metal hyperaccumulator Arabidopsis halleri. J. Exp. Bot. 66, 3215–3227. doi: 10.1093/jxb/erv144
Milne, T. J., Abbenante, G., Tyndall, J. D., Halliday, J., and Lewis, R. J. (2003). Isolation and characterization of a cone snail protease with homology to CRISP proteins of the pathogenesis-related protein superfamily. J. Biol. Chem. 278, 31105–31110. doi: 10.1074/jbc.M304843200
Nocito, F. F., Lancilli, C., Dendena, B., Lucchini, G., and Sacchi, G. A. (2011). Cadmium retention in rice roots is influenced by cadmium availability, chelation and translocation. Plant Cell Environ. 34, 994–1008. doi: 10.1111/j.1365-3040.2011.02299.x
Ogawa, I., Nakanishi, H., Mori, S., and Nishizawa, N. K. (2009). Time course analysis of gene regulation under cadmium stress in rice. Plant Soil. 325, 97–108. doi: 10.1007/s11104-009-0116-9
Park, J., Song, W. Y., Ko, D., Eom, Y., Hansen, T. H., Schiller, M., et al. (2012). The phytochelatin transporters AtABCC1 and AtABCC2 mediate tolerance to cadmium and mercury. Plant J. 69, 278–288. doi: 10.1111/j.1365-313X.2011.04789.x
Peng, H., He, X., Gao, J., Ma, H., Zhang, Z., Shen, Y., et al. (2015). Transcriptomic changes during maize roots development responsive to Cadmium (Cd) pollution using comparative RNAseq-based approach. Biochem. Bioph. Res. Co. 464, 1040–1047. doi: 10.1016/j.bbrc.2015.07.064
Remy, E., Cabrito, T. R., Batista, R. A., Hussein, M. A., Teixeira, M. C., Athanasiadis, A., et al. (2014). Intron retention in the 5' UTR of the novel ZIF2 transporter enhances translation to promote zinc tolerance in Arabidopsis. PLoS Genet. 10:e1004375. doi: 10.1371/journal.pgen.1004375
Rodríguez-Serrano, M., Romero-Puertas, M. C., Sparkes, I., Hawes, C., Luis, A., and Sandalio, L. M. (2009). Peroxisome dynamics in Arabidopsis plants under oxidative stress induced by cadmium. Free Radical Biol. Med. 47, 1632–1639. doi: 10.1016/j.freeradbiomed.2009.09.012
Rowley, M. J., Böhmdorfer, G., and Wierzbicki, A. T. (2013). Analysis of long non-coding RNAs produced by a specialized RNA polymerase in Arabidopsis thaliana. Methods 63, 160–169. doi: 10.1016/j.ymeth.2013.05.006
Sakai, H., Lee, S. S., Tanaka, T., Numa, H., Kim, J., Kawahara, Y., et al. (2013). Rice Annotation Project Database (RAP-DB): an integrative and interactive database for rice genomics. Plant Cell Physiol. 54, e6. doi: 10.1093/pcp/pcs183
Sandalio, L. M., Dalurzo, H. C., Gómez, M., Romero−Puertas, M. C., and del Río, L. A. (2001). Cadmium−induced changes in the growth and oxidative metabolism of pea plants. J. Exp. Bot. 52, 2115–2126. doi: 10.1093/jexbot/52.364.2115
Sasaki, A., Yamaji, N., Yokosho, K., and Ma, J. F. (2012). Nramp5 is a major transporter responsible for manganese and cadmium uptake in rice. Plant Cell. 24, 2155–2167. doi: 10.1105/tpc.112.096925
Satoh-Nagasawa, N., Mori, M., Nakazawa, N., Kawamoto, T., Nagato, Y., Sakurai, K., et al. (2012). Mutations in rice (Oryza sativa) heavy metal ATPase 2 (OsHMA2) restrict the translocation of zinc and cadmium. Plant Cell Physiol. 53, 213–224. doi: 10.1093/pcp/pcr166
Shen, S., Park, J. W., Lu, Z. X., Lin, L., Henry, M. D., Wu, Y. N., et al. (2014). rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc. Natl. Acad. Sci. U.S.A. 111, E5593–E5601. doi: 10.1073/pnas.1419161111
Slamet-Loedin, I. H., Johnson-Beebout, S. E., Impa, S., and Tsakirpaloglou, N. (2015). Enriching rice with Zn and Fe while minimizing Cd risk. Front. Plant Sci. 6:121. doi: 10.3389/fpls.2015.00121
Staiger, D., and Brown, J. W. (2013). Alternative splicing at the intersection of biological timing, development, and stress responses. Plant Cell. 25, 3640–3656. doi: 10.1105/tpc.113.113803
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 17, e166. doi: 10.1093/nar/gkt646
Takahashi, R., Ishimaru, Y., Shimo, H., Ogo, Y., Senoura, T., Nishizawa, N. K., et al. (2012). The OsHMA2 transporter is involved in root-to-shoot translocation of Zn and Cd in rice. Plant, Cell Environ. 35, 1948–1957. doi: 10.1111/j.1365-3040.2012.02527.x
Tang, M., Mao, D., Xu, L., Li, D., Song, S., and Chen, C. (2014). Integrated analysis of miRNA and mRNA expression profiles in response to Cd exposure in rice seedlings. BMC Genomics. 15:835. doi: 10.1186/1471-2164-15-835
Thomine, S., Wang, R., Ward, J. M., Crawford, N. M., and Schroeder, J. I. (2000). Cadmium and iron transport by members of a plant metal transporter family in Arabidopsis with homology to Nramp genes. Proc. Natl. Acad. Sci. U.S.A. 97, 4991–4996. doi: 10.1073/pnas.97.9.4991
Trapnell, C., Pachter, L., and Salzberg, S. L. (2009). TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25, 1105–1111. doi: 10.1093/bioinformatics/btp120
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Ueno, D., Yamaji, N., Kono, I., Huang, C. F., Ando, T., Yano, M., et al. (2010). Gene limiting cadmium accumulation in rice. Proc. Natl. Acad. Sci. U.S.A. 107, 16500–16505. doi: 10.1073/pnas.1005396107
Uraguchi, S., Kamiya, T., Clemens, S., and Fujiwara, T. (2014). Characterization of OsLCT1, a cadmium transporter from indica rice (Oryza sativa). Physiol. Plant. 151, 339–347. doi: 10.1111/ppl.12189
Uraguchi, S., Kamiya, T., Sakamoto, T., Kasai, K., Sato, Y., Nagamura, Y., et al. (2011). Low-affinity cation transporter (OsLCT1) regulates cadmium transport into rice grains. Proc. Natl. Acad. Sci. U.S.A. 108, 20959–20964. doi: 10.1073/pnas.1116531109
Van Loon, L. C., and Van Strien, E. A. (1999). The families of pathogenesis-related proteins, their activities, and comparative analysis of PR-1 type proteins. Physiol. Mol. Plant P. 55, 85–97. doi: 10.1006/pmpp.1999.0213
Verrier, P. J., Bird, D., Burla, B., Dassa, E., Forestier, C., Geisler, M., et al. (2008). Plant ABC proteins-a unified nomenclature and updated inventory. Trends Plant Sci. 13, 151–159. doi: 10.1016/j.tplants.2008.02.001
Wang, F., Wang, M., Liu, Z., Shi, Y., Han, T., Ye, Y., et al. (2015a). Different responses of low grain-Cd-accumulating and high grain-Cd-accumulating rice cultivars to Cd stress. Plant Physiol. Biochem. 96, 261–269. doi: 10.1016/j.plaphy.2015.08.001
Wang, M., Chen, W., and Peng, C. (2015b). Risk assessment of Cd polluted paddy soils in the industrial and township areas in Hunan, Southern China. Chemosphere 144, 346–351. doi: 10.1016/j.chemosphere.2015.09.001
Wang, W., Vinocur, B., Shoseyov, O., and Altman, A. (2004). Role of plant heat-shock proteins and molecular chaperones in the abiotic stress response. Trends Plant Sci. 9, 244–252. doi: 10.1016/j.tplants.2004.03.006
Wang, Y., Wang, Q., Zhao, Y., Han, G., and Zhu, S. (2015c). Systematic analysis of maize class III peroxidase gene family reveals a conserved subfamily involved in abiotic stress response. Gene 566, 95–108. doi: 10.1016/j.gene.2015.04.041
Weber, M., Trampczynska, A., and Clemens, S. (2006). Comparative transcriptome analysis of toxic metal responses in Arabidopsis thaliana and the Cd2+ −hypertolerant facultative metallophyte Arabidopsis halleri. Plant Cell Environ. 29, 950–963. doi: 10.1111/j.1365-3040.2005.01479.x
Wei, W., Chai, T., Zhang, Y., Han, L., Xu, J., and Guan, Z. (2009). The Thlaspi caerulescens NRAMP homologue TcNRAMP3 is capable of divalent cation transport. Mol. Biotechnol. 41, 15–21. doi: 10.1007/s12033-008-9088-x
Wilusz, J. E., Sunwoo, H., and Spector, D. L. (2009). Long noncoding RNAs: functional surprises from the RNA world. Gene. Dev. 23, 1494–1504. doi: 10.1101/gad.1800909
Wysocki, R., and Tamás, M. J. (2010). How Saccharomyces cerevisiae copes with toxic metals and metalloids. FEMS Microbiol. Rev. 34, 925–951. doi: 10.1111/j.1574-6976.2010.00217.x
Xin, M., Wang, Y., Yao, Y., Song, N., Hu, Z., Qin, D., et al. (2011). Identification and characterization of wheat long non-protein coding RNAs responsive to powdery mildew infection and heat stress by using microarray analysis and SBS sequencing. BMC Plant Biol. 11:61. doi: 10.1186/1471-2229-11-61
Xu, L., Wang, Y., Liu, W., Wang, J., Zhu, X. W., Zhang, K. Y., et al. (2015). De novo sequencing of root transcriptome reveals complexcadmium-responsive regulatory networks in radish(Raphanus sativus L.). Plant Sci. 236, 313–323. doi: 10.1016/j.plantsci.2015.04.015
Xue, G. P., Way, H. M., Richardson, T., Drenth, J., Joyce, P. A., and McIntyre, C. L. (2011). Overexpression of TaNAC69 leads to enhanced transcript levels of stress up-regulated genes and dehydration tolerance in bread wheat. Mol. Plant. 4, 697–712. doi: 10.1093/mp/ssr013
Yoneyama, T., Ishikawa, S., and Fujimaki, S. (2015). Route and regulation of zinc, cadmium, and iron transport in rice plants (Oryza sativa L.) during vegetative growth and grain filling: metal transporters, metal speciation, grain Cd reduction and Zn and Fe biofortification. Int. J. Mol. Sci. 16, 19111–19129. doi: 10.3390/ijms160819111
Yu, S., Ligang, C., Liping, Z., and Diqiu, Y. (2010). Overexpression of OsWRKY72 gene interferes in the abscisic acid signal and auxin transport pathway of Arabidopsis. J. Biosci. 35, 459–471. doi: 10.1007/s12038-010-0051-1
Yuan, L., Yang, S., Liu, B., Zhang, M., and Wu, K. (2012). Molecular characterization of a rice metal tolerance protein, OsMTP1. Plant Cell Rep. 31, 67–79. doi: 10.1007/s00299-011-1140-9
Zhang, L., Zhao, G., Jia, J., Liu, X., and Kong, X. (2012). Molecular characterization of 60 isolated wheat MYB genes and analysis of their expression during abiotic stress. J. Exp. Bot. 63, 203–214. doi: 10.1093/jxb/err264
Zhao, F. Y., Liu, W., and Zhang, S. Y. (2009). Different responses of plant growth and antioxidant system to the combination of cadmium and heat stress in transgenic and non-transgenic rice. J. Integr. Plant Biol. 51, 942–950. doi: 10.1111/j.1744-7909.2009.00865.x
Zheng, L., Yamaji, N., Yokosho, K., and Ma, J. F. (2012). YSL16 is a phloem-localized transporter of the copper-nicotianamine complex that is responsible for copper distribution in rice. Plant Cell. 24, 3767–3782. doi: 10.1105/tpc.112.103820
Keywords: rice root, cadmium, RNA-seq, alternative splicing, lncRNAs
Citation: He F, Liu Q, Zheng L, Cui Y, Shen Z and Zheng L (2015) RNA-Seq Analysis of Rice Roots Reveals the Involvement of Post-Transcriptional Regulation in Response to Cadmium Stress. Front. Plant Sci. 6:1136. doi: 10.3389/fpls.2015.01136
Received: 05 September 2015; Accepted: 30 November 2015;
Published: 21 December 2015.
Edited by:
Mohammad Anwar Hossain, Bangladesh Agricultural University, BangladeshReviewed by:
Françoise Gosti, Centre National de la Recherche Scientifique, FranceCaiyan Chen, Institute of Subtropical Agriculture Chinese Academy of Sciences, China
Yingchun Xu, Nanjing Agricultural University, China, China
Copyright © 2015 He, Liu, Zheng, Cui, Shen and Zheng. 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.
*Correspondence: Luqing Zheng, zhenglq@njau.edu.cn
†These authors have contributed equally to this work.