Cd accumulation during Cd-stress
Under the control condition, node I (as “N”) in the two cultivars (21.05 mg/kg DW in “Yuzhenxiang” as “y” and 10.25 mg/kg DW in “Xiangwanxian No. 12” as “X”) had higher Cd accumulation compared with panicle node (as “P”) (2.17 mg/kg DW in “y” and 1.40 mg/kg DW in “X”; p<0.01, Fig. 1). The Cd treatment increased the Cd accumulation in all tissues, especially that in the node I (56.43 mg/kg DW in “y” and 44.25 mg/kg DW in “X”). Grains of “X” cultivar had lower Cd accumulation (both in control and Cd treatment) than that in “y” cultivar. These data confirmed that “y” was a high Cd accumulation cultivar, and node I had higher capacity in Cd sequestration.
Summary of the mRNA-seq and miRNA-seq
Twenty-four cDNA libraries were constructed and a total of 1111.34 M clean reads, with 89.80% average mapping rate (76.45%–92.37%), to the rice reference genome (Additional file 2: Table S1). Principle component analysis (PCA) (Fig. 2A) and sample-to-sample clustering analysis (Fig. 2B) showed that the samples of the same tissue (“P” or “N”) from the same cultivar of control (“C”) and Cd treatment (“T”) were clustered together. Twenty-four small cRNA libraries generated total 350.83 M clean reads (112.57 M unique reads) with an averaged mapping rate of 82.01% to the O. sativa reference miRNA in miRbase (http://www.mirbase.org/cgi-bin/mirna_summary.pl?org = osa). The alignment rate in miRBase of each sample sequenced ranged from 0.62% to 2.37% (average 1.49%,Additional file 2: Table S1). PCA (Fig. 2C) and sample-to-sample clustering analysis (Fig. 2D) showed that samples of the same tissue (“P” or “N”) were grouped together.
DEGs between high and low Cd accumulation cultivars
A total of 4,535 DEGs were identified by pairwise comparison of Cd-treated (T) vs. untreated control (C) of node I (N) and panicle node (P) in “X” and “y” cultivar, respectively (Fig. 3A-B and Table 1).The results showed that the number of down-regulated genes in “X” cultivar was greater than that in “y” cultivar. GO and KEGG enrichment based on these down regulated genes were performed to identify the main biological processes. The results showed that GO terms of “transporter activity” and “response to stimulus” were significantly enriched in “X”. By contrast, only a few genes were differentially expressed in these two GO terms in “y” cultivar (Fig. 4A). Results of hierarchical clustering analysis of the 84 common DEGs, including 69 genes down regulated after Cd treatment in “X”, showed an enrichment in “transporter activity” (Fig. 4B and Additional file 3: Table S2), whereas another 74 common DEGs, including 62 genes down regulated after Cd treatment in “X”, showed an enrichment related to “response to stimulus” (Fig. 4C and Additional file 4: Table S3). Notably, most of the DEGs in “X” were down-regulated by Cd treatment, and most down-regulated genes in “X” by Cd treatment were unchanged in “y”, especially for the DEGs associated with “response to stimulus” (Fig. 4C).
Among all the DEGs related to “transporter activity”, the iron-regulated transporter 1 (IRT1, OS03G0667500) and metal transporter Nramp5 (natural resistance-associated macrophage protein 5, Mg2+ and Cd transporter, Mn and Cd uptake protein, OS07G0257200) were noteworthy for they had higher expression in the panicle node compared with node I (Additional file 3: Table S2), which indicated that Cd treatment increased OsIRT1 and OsNramp5 in the panicle node of the two cultivar. For node I, OsIRT1 and OsNramp5 increased in “X”, while decreased in “y” cultivar (Additional file 1: Fig. S1 A). Although the expression of OsIRT1 and OsNramp5 increased after Cd treatment, overall it remained at a low level in “X” cultivar.
OsNRT1.5A (OS02G0689900, nitrate transporter 1.5A) and OsVIT2 (OS09G0396900, Vacuolar Iron Transporter 2) played a key role in the upward transport of Cd. OsNRT1.5A showed a higher expression level in node I than panicle node in both cultivars. In “X” cultivar with Cd treatment, OsNRT1.5A was up- and down-regulated in node I and panicle node, respectively. For the expression of OsNRT1.5A in “y” cultivar, it was up regulated in panicle node. OsVIT2 was up regulated following Cd treatment in node I in both cultivars, but it was down regulated in panicle node in the “X” cultivar. Cd stress had little effects on the expression level of OsVIT2 in panicle node in “y” cultivar (Additional file 1: Fig. S1A and Additional file 3:Table S2). With Cd treatment, 6 aquaporin genes (PIPs) (OS02G0629200, OS04G0233400, OS02G0666200, OS04G0559700, OS07G0448100 and OS07G0448800) showed similar expression profiles (Additional file 1: Fig. S1B). The expression of 6 PIPs members was higher in “X” than in “y” cultivar under control condition, and down regulated after Cd treatment in “X”, but with little changes in “y”. The expression of PIPs did not show significant difference between node I and panicle node in two rice cultivars.
Among the DEGs related to “response to stimulus”, heat shock transcription factor (HSF) A2d/B2c genes (including OS03G0161900/OS09G0526600), light-harvesting chlorophyll a-b binding protein (LHC-II) genes (e.g. OS02G0197600), and genes encoding heat shock protein (HSP71.1/70/20; including OS03G0276500, OS01G0840100 and OS06G0253100) showed similar changes in the panicle node and node I of the two rice cultivars (Additional file 4: Table S3). Higher expression levels of aforementioned genes were found in “X” than in “y” cultivar (Additional file 1: Fig. S1C). Cd treatment decreased the expression of all these genes in the “X” cultivar, but increased OsHSF-A2d/B2c, OsLHC-II and OsHSP71.1 in the panicle nodeof “y” cultivar. The distinct expression profiles and responses of aforementioned genes to Cd treatment between the two rice cultivars are likely in part account for the differential Cd accumulation.
Expression of known Cd-responsive genes
In order to decipher the Cd-stress responses in panicle node and node I, we analyzed the expression profiles of 52 Cd-responsive genes previously reported in the literature. Among these genes associated with Cd-stress in rice, including metallothionein 1 (OsMT1) [16], cadmium tolerant 1 (OsCDT1),, OsCDT2 [17, 18], OsMTP1 [7], cation diffusion facilitator (OsCDF1),, ATP-binding cassette transporter multidrug resistance protein 1 (OsMRP1/ABCC1) [19, 20] showed higher expression levels in the two rice cultivars. Furthermore, they expressed higher in the panicle node of the two cultivars compared with node I (Fig. 5). Cd enhanced the expression of OsCDT1 and OsABCC1 in the node I of “X” cultivar. Another gene serine hydroxymethyltransferase 1 (OsSHM1),, which showed relatively higher expression level, was down regulated only in the panicle node of “X” cultivar by Cd stimulus (Fig. 5). In addition, the nicotinamine synthase 3 (OsNAS3) gene and dense and erect panicle 1 (OsDEP1) genes, which had low expression levels, showed relatively high expression in node I of the two cultivars. Cd treatment enhanced the expression of OsNAS3 in node I and that of OsMT1d in both nodes (Fig. 5). Other Cd-responsive genes including OsMT1f, OsYSL15, OsIRT2, and OsGST4 had low expression levels in the two cultivars and were found Cd unresponsive in our study.
Identification of differentially expressed DEmiRNAs
A total of 70 non-overlapping DEmiRNAs were identified from panicle node and node I in the two rice cultivars (Fig. 6A and Additional file 5: Table S4). Most DEmiRNAs were up-regulated by Cd treatment. There were 12 common DEmiRNAs among all the pairwise comparisons (Fig. 6B), of which six miRNAs (osa-miR398b, osa-miR408–3p, osa-miR408–5p, osa-miR528–5p, osa-miR528–3p and novel50_mature) were up-regulated by Cd treatment (Fig. 6C). In addition, only one common DEmiRNA (osa-miR528–3p) was identified as Cd up-regulated in the panicle node of both cultivars. There were six common DEmiRNAs (osa-miR408–3p, osa-miR408–5p, osa-miR528–5p, osa-miR398b, osa-miR166–5p and novel50_mature) in node I of both cultivars, among which, only osa-miR166–5p showed a different expression pattern in node I of “X” cultivar and panicle node of “y” cultivar. In order to reveal the differential expression profiles of microRNAs more comprehensively, we screened the other miRNA family members (how you select these members) (Fig. 6D). Osa-miR398b, osa-miR408–3p, osa-miR528–3p and osa-miR528–5p in panicle node showed a higher expression level than node I with Cd treatment. We then constructed the DE miRNA-mRNA regulatory network (Fig. 7) based on the 15 miRNA (12 common and 3 family members) and the targets among the down-regulated DEGs. In the regulatory network, a HSP member (OS06G0253100) was regulated by two miRNAs including osa-miR528–5p and osa-miR5493. OsMYB5P (OS02G0624300), OsbZIP18 (Os02g0203000), and OsERF141 (Os02g0638650) were regulated by osa-miR528–5p and novel13_mature, indirectly. Another bZIP member OsbZIP23 (Os02g0766700) was the target gene of osa-miR1846a/b/c–5p; SNF1-related protein kinase 1 subfamily protein (SnRK) gene (OS02G0178000, OsSnRK1) was regulated by osa-miR528–3p. In addition, OsAAE3 (OS04G0683700) was regulated by both osa-miR408–3p and novel50_mature (Fig. 7 and Additional file 1: Fig. S1D).
qRT-PCR verification of RNA-seq data
A total of 14 genes were selected for qRT-PCR analysis. The fold-changes found in qRT-PCR and RNA-seq were highly consistent, and the correlation coefficient was 0.6 (Fig. 8).