RNA-Seq Analysis of Rice Roots Reveals the Involvement of Post-Transcriptional Regulation in Response to Cadmium Stress

Abstract

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.

Keywords: RNA-seq; alternative splicing; cadmium; lncRNAs; rice root.

Figure 1

Flowchart of our RNA-seq data analysis workflow. High-level representation of the pipeline workflow for processing RNA-Seq data.

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.

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.

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.

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.

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.